MOCCA-SURVEY database I. Accreting white dwarf binary systems in globular clusters -- III. Cataclysmic variables -- Implications of model assumptions

TL;DR

This study analyzes the population of cataclysmic variables in globular clusters using advanced models, revealing extremely low duty cycles and the influence of initial binary populations and model assumptions on CV properties.

Contribution

It provides new insights into how different initial conditions and model parameters affect the formation and observable characteristics of CVs in globular clusters.

Findings

Present-day CV duty cycles are extremely low (<0.1%)

Initial binary populations significantly influence CV properties

Bright CVs are typically recently formed, while faint CVs are older

Abstract

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.

Equations

$$\frac{{\dot{J}}_{\mathrm{CAML}}}{J}=\nu \frac{{\dot{M}}_2}{M_2}$$

$$\nu\ =\ \begin{dcases}\frac{M_{2}^{2}}{M_{1}(M_{1}+M_{2})},&\text{classical \cite[citep]{(\@@bibref{AuthorsPhrase1Year}{King_1995}{\@@citephrase{, }}{})}}\\ &\\ \frac{0.35}{M_{1}},&\text{empirical \cite[citep]{(\@@bibref{AuthorsPhrase1Year}{Schreiber_2016}{\@@citephrase{, }}{})}}\end{dcases}$$

$$E_b=\alpha \Delta E_{\mathrm{orb}}$$

$$E_b=\frac{GM(M-M_c)}{\lambda R}.$$

$$E_b=-{\int }_{M_c}^M\frac{Gm}{r(m)}dm+{\alpha }_{\mathrm{rec}}{\int }_{M_c}^M{\epsilon }_{\mathrm{rec}}(m)dm$$

Full text

††pagerange: MOCCA-SURVEY database I. Accreting white dwarf binary systems in globular clusters – III. Cataclysmic variables – Implications of model assumptions–MOCCA-SURVEY database I. Accreting white dwarf binary systems in globular clusters – III. Cataclysmic variables – Implications of model assumptions††pubyear: 2016

1 INTRODUCTION

Cataclysmic variables (CVs) are interacting binaries composed of a white dwarf (WD) undergoing stable mass transfer from a main sequence (MS) star or a brown dwarf (BD) (e.g. Warner, 1995; Knigge et al., 2011). They are expected to exist in non-negligible numbers in globular clusters (GCs) that are natural laboratories for testing theories of stellar dynamics and evolution.

CVs in GCs have been studied by many authors, both theoretically and observationally (e.g. Knigge, 2012, and references therein). GCs are thought to play a crucial role in CV formation, since their densities are sufficiently high that dynamical encounters involving binaries should be common. Thus, in dense GCs, it is natural to expect that many CV progenitors will have been affected by dynamics in some way prior to CV formation (e.g. Shara & Hurley, 2006; Ivanova et al., 2006; Belloni et al., 2016, 2017; Hong et al., 2017).

In the first paper of this series (Belloni et al., 2016), we discussed six specific MOCCA (MOnte Carlo Cluster simulAtor) models with a focus on the properties of their present-day CV populations. In the second paper (Belloni et al., 2017), we concentrated instead on a discussion of the properties of the progenitor and formation-age populations, and how CV properties are affected by their ages.

Here we aim to complement our previous work by considering the same six initial cluster conditions that we already presented but considering one new combination of the common-envelope phase (CEP) parameters and two new approaches to the treatment of CV evolution and angular momentum loss (AML). We analyze the impact of these variables on present-day CV populations.

The paper is structured as follows. In Section 2, we describe the codes used for the cluster and CV evolution simulations. We present the suite of models analyzed in this paper in Section 3. The main results of this investigation are presented in Sections 4, 5, 6, and 7, while Section 8 contains a discussion of the main results, along with a direct comparison between our simulated CV properties and observed CV candidates in GCs. Finally, we summarize our main conclusions in Section 9.

2 METHODOLOGY

We start with a summary of the main aspects of the codes used in Belloni et al. (2016, 2017), defining as well the different dynamical groups and the different populations that we will analyze. In Section 2.2 we explain how and why we introduced a new code for simulating close WD-MS binaries and CV evolution, which is up-to-date and current. Finally, in Section 2.3 we describe the three approaches adopted in this work to treat the AML during the CV phase.

The cluster evolution is performed with the MOCCA code, which is based on the orbit-averaged Monte Carlo technique from Hénon (1971) improved by Stodółkiewicz (1986) and further improved by Giersz et al. (2013, and references therein). These last authors included in MOCCA the FEWBODY code (Fregeau et al., 2004) to handle 3 and 4-body interactions, a description of escape processes in tidally limited clusters (Fukushige & Heggie, 2000), and stellar evolution via the SSE/BSE code developed by Hurley et al. (2000, 2002). Additionally, MOCCA has been extensively tested against N-body codes (e.g. Giersz et al., 2013; Wang et al., 2016), showing good agreement.

The analysis of the CVs produced in clusters simulated by the MOCCA code was performed with the CATUABA (Code for Analysing and sTUdying cAtaclysmic variables, symBiotic stars and AM CVns) code.

A detailed description of the MOCCA and CATUABA codes, as well as all the features and main assumptions adopted in them, can be found in Belloni et al. (2016, 2017).

2.1 The influence of dynamics on CV formation

As already explained in Belloni et al. (2016, see section 3.2, for more details), CATUABA permits for the unambiguous identification of a given CV population in the cluster snapshots provided by MOCCA at various times in the cluster evolution. Based on the most relevant events in the history of each individual CV, systems are separated into three groups:

• •

Binary stellar evolution (BSE) group, which corresponds to CVs formed without any influence from dynamics.

• •

Weak dynamical interactions (WDI) group, i.e. CVs weakly influenced by dynamics.

• •

Strong dynamical interaction (SDI) group, which contains those CVs strongly influenced by dynamical interactions, including systems that underwent either exchanges or mergers.

The CATUABA code also recognizes the main formation channels for the WD-MS population that will later on become the CV population. Briefly, close WD-MS binaries can be formed through a CEP with/without weak/strong dynamical interactions, through exchanges with/without a CEP, and through mergers with/without a CEP.

Three different populations can also be distinguished:

• •

The progenitor population is defined as the population of all binaries that are CV progenitors, i.e. the population of binary systems at the time of cluster birth (i.e., $t=0$) that will become CVs.

• •

The formation-age population is defined as the population of zero-age CVs, i.e. when mass transfer starts from the donor onto the WD. Obviously, the time corresponding to the onset of the CV phase is different for every CV.

• •

The present-day population is the population of CVs at the current age (t $=$ 12 Gyrs).

2.2 Close WD-MS binary evolution and CV evolution

The most important part of MOCCA for producing CVs and related objects is the Binary Star Evolution (BSE) code from Hurley et al. (2002), which models angular momentum loss in binaries. However, its implementation is not fully up-to-date (Knigge et al., 2011; Schreiber et al., 2016, for a state-of-the-art model). Mass transfer occurs if either star fills its Roche lobe and may proceed on a nuclear, thermal, or dynamical time-scale. Prescriptions to determine the type and rate of mass transfer and the response of the primary to accretion are implemented in BSE. For the CEP, BSE assumes that the common envelope binding energy is that of the giant(s) envelope involved in the process.

The overall CV evolution can be recovered by BSE, although the CV evolution model implemented in this code is out-dated (for more details, see Section 8.3). This leads to unrealistic predictions for the location of the period minimum and the subsequent CV evolution. Additionally, BSE is unable to reproduce the orbital period gap, which is one of the most important observables related to CVs.

Since the BSE code is out-dated with respect to many aspects of CV evolution, we use instead a different code to simulate the entire population of close WD-MS binaries formed in the simulated clusters. For this purpose, we implement in CATUABA an option to extract all close WD-MS binaries (with periods less than 10 days) that are formed throughout the cluster evolution. The systems are then evolved up to 12 Gyr with the code described in Zorotovic et al. (2016).

We emphasize that in order to proceed with this method, it is reasonable to assume that once a close WD-MS binary has formed then its subsequent evolution will not be influenced by dynamics based on the results of Belloni et al. (2017), who showed that CVs are very dynamically hard, and of low mass, such that the impact parameter for interactions is small. It follows that the probability of an interaction occurring is very small (see Leigh et al., 2016, for more details regarding the interruption of binary mass transfer due to dynamical interactions). Thus, assuming that the evolution of close WD-MS binaries is not affected by dynamics is not far from what happens in simulations of realistic clusters.

2.3 Angular Momentum Loss and stability limit

It is well-known that AML in a CV determines its secular evolution. It is important to distinguish AML due to mass transfer, usually called consequential angular momentum loss (CAML, $\dot{J}_{\rm CAML}$) from AML independent of mass transfer (e.g. magnetic braking, gravitational radiation), usually called systemic AML ($\dot{J}_{\rm sys}$) (e.g. King & Kolb, 1995). Typical candidates for the CAML are: circumbinary disks (e.g. Willems et al., 2005), hydromagnetic accretion disk winds (e.g. Cannizzo & Pudritz, 1988), and frictional drag associated with nova eruptions (Schreiber et al., 2016; Nelemans et al., 2016).

Particularly interesting is the influence of nova eruptions on CV evolution. According to Schreiber et al. (2016), if the strength of CAML is inversely proportional to the WD mass, then most CVs with WD masses less than $\sim$ 0.5 ${\rm M_{\odot}}$ are dynamically unstable. The motivation for such a functional form (empirical CAML) comes from the fact that a helium (He) WD has never been observed in a CV in the Galactic field, although such WDs are observed in detached systems (Zorotovic et al., 2011). The primary mechanism thought to be responsible for such a dependence between CAML and WD mass is nova eruptions (Schreiber et al., 2016; Nelemans et al., 2016). The frictional AML produced by novae depends strongly on the expansion velocity of the ejecta (Schenker et al., 1998). For low-mass WDs, the expansion velocity is small (Yaron et al., 2005), and this may lead to strong AML by friction that makes most CVs with He WDs dynamically unstable.

The empirical CAML model is a good candidate to solve several problems related to CV evolution, like the missing low-mass WDs in CVs, the predicted versus observed space density or the period distribution (see Schreiber et al., 2016, for more details). Recently, this model has also been proposed as the explanation for the existence of single He-core WDs (Zorotovic & Schreiber, 2017). Therefore, and despite the fact that the physical mechanism behind the new empirical form of CAML is not understood at a fundamental level, we think it is appropriate to check the influence that this enhanced CAML can have on GC CVs.

Both codes used in this work for CV evolution do not follow the evolution of the mass transfer rate during nova cycles but just calculate the secular mean mass transfer rate. Mass transfer variations that may occur on shorter time scales between two nova eruptions (e.g. Shara et al., 1986, 2007; Schmidtobreick et al., 2015) are thus not considered. As shown by Schenker et al. (1998) this should not affect the stability limit which is fundamental for the predictions of binary population models. In both codes the mass loss due to nova eruption is simply taken into account by assuming that only a small fraction of the transferred mass is accreted by the WD. In BSE this fraction is 0.1 per cent (Hurley et al., 2002, Eq. 66) so the WD mass is assumed to slowly grow. Zorotovic et al. (2016) assume that all accreted matter is expelled during nova eruptions i.e. the WD mass is kept constant.

Both the BSE code and the code described in Zorotovic et al. (2016) use systemic AML given by the disrupted magnetic braking model (Rappaport et al., 1983; Schreiber et al., 2010) for the evolution of close detached binaries. This model assumes that magnetic braking is not active, or at least strongly inefficient, in systems with fully convective MS secondary stars. However, when mass transfer starts, Zorotovic et al. (2016) include CAML as an extra source of AML, which affects the orbital period evolution and the critical mass ratio for stable mass transfer. In BSE, on the other hand, CAML is assumed to be negligible and therefore is not included.

In order to check the influence of CAML on CVs in GCs, in addition to those CVs formed in MOCCA, we simulate all close WD-MS formed in our models with the code described in Zorotovic et al. (2016). We consider in this investigation two formulations for CAML, namely the classical CAML (cCAML) from (King & Kolb, 1995) and the empirical formulation (eCAML) from Schreiber et al. (2016), given below:

[TABLE]

where

[TABLE]

For both CAML approaches, we compute the stability limit imposed by the adopted formulation by equating the adiabatic mass radius exponent and the mass radius exponent of the secondary’s Roche lobe. In BSE, on the other hand, a fixed value for the critical mass ratio is assumed.

3 MODELS

In what follows we describe the main initial conditions assumed for the GCs in our different models. The most important differences concern the initial binary population (IBP), initial binary fraction and initial central density. We also explain our choice to introduce a new set of CEP parameters in order to generate six new models.

3.1 Initial cluster conditions

The IBP refers to all initial binaries in a given cluster, associated with specific distributions for their orbital parameters. In this work, we continue the analysis started in Belloni et al. (2016, 2017), who investigated models with two distinct IBPs:

• •

K models: models with 95 per cent primordial binary fractions with orbital parameter distributions that follow the Kroupa IBP, which is based on binary eigenevolution and a mass feeding algorithm (Kroupa, 2008).

• •

S models: models with low binary fractions of either 5 or 10 per cent that follow the “Standard” initial binary orbital parameter distributions. The Standard IBP assumes an uniform distribution for the mass ratio, an uniform distribution in the logarithm of the semi-major axis $\log(a)$ or a log-normal semi major axis distribution, and a thermal eccentricity distribution.

The main distributions associated with both IBPs are shown in Figures 6 and 7 in the Appendix A.

The models also differ with respect to the initial central density. We consider: sparse ($\rho_{c}\sim 10^{3}$ M⊙ pc*-3*), dense ($\rho_{c}\sim 10^{5}$ M⊙ pc*-3*), and very dense ($\rho_{c}>10^{5}$ M⊙ pc*-3*) clusters, where $\rho_{c}$ is the cluster central density.

We assume that all stars are on the zero-age main sequence when the simulations begin, and that all residual gas left over from the star formation process has already been expelled from the cluster. Additionally, all models have low metallicity ($\lesssim$ 0.001), are initially in virial equilibrium, and have neither rotation nor primordial mass segregation (e.g. Leigh et al., 2013, 2015). Finally, all models are evolved for 12 Gyr, up to the present-day.

The properties of the six models are summarized in Table 1.

3.2 Common-envelope phase

An important phase in the life of most interacting binaries (such as CVs, AM CVns, X-ray binaries, etc.) is the CEP during which time two stars share the same envelope (see Ivanova et al., 2013, for a review on common-envelope evolution). This phase leads to binaries whose orbital periods are orders of magnitude shorter than prior to this phase, due to common envelope ejection and gas dynamical friction.

Although a great deal of effort has been put in to understanding the CEP, the energy budget during this phase is not well understood, especially the role of recombination energy stored in ionized regions inside the envelope (Ivanova et al., 2015).

For convenience, the outcome of this phase is usually approximated by the balance between the change in the orbital energy and envelope binding energy, given an efficiency $\alpha$ for this process. More specifically, $\alpha$ is the efficiency with which orbital energy is used to eject the envelope, i.e. it corresponds to the fraction of the difference in orbital energy (before and after the CE) in unbinding the envelope, i.e.

[TABLE]

where $E_{b}$ is the envelope binding energy, $E_{\rm orb}$ is the orbital energy, and $\alpha$ is the fraction of $E_{\rm orb}$ used to unbind the envelope (i.e. CEP efficiency). The binding energy is usually approximated by

[TABLE]

where $\lambda$ is the binding energy parameter, which depends on the structure of the primary star.

If a fraction $\alpha_{\rm rec}$ of the recombination energy of the envelope is used to help in ejecting the envelope, then the binding energy equation is

[TABLE]

where $M_{c}$ is the giant core mass, $M$ is the giant mass, $\varepsilon_{\rm rec}$ is the specific recombination energy, and $\alpha_{\rm rec}$ is the fraction of the recombination energy used to unbind the envelope. The binding energy parameter $\lambda$ can be computed by equating equations 3 and 4. When $\alpha_{\rm rec}>0$, an additional source of energy is used in the ejection of the envelope, which implies that less orbital energy is needed to unbind the envelope.

In the six models investigated by Belloni et al. (2016, 2017), we assumed $\alpha=3$, which might not be an appropriate choice based on the study carried out by Zorotovic et al. (2010). These authors analyzed a homogeneous sample of 60 post common envelope binaries (PCEBs) identified by the Sloan Digital Sky Survey. They found that $\alpha$ should be $\sim$ 0.2-0.3, and certainly no greater than 1. In our previous works, we also used $\lambda=0.5$. As pointed out by Zorotovic et al. (2014a), in the original version of BSE, a fixed value for the binding energy parameter ($\lambda$) is assumed as input to the code. However, after improvements were made, a function to compute $\lambda$ was implemented in BSE, as described in detail by Claeys et al. (2014, see their Appendix A). In this newer implementation, a positive value of $\lambda$ represents the fraction of the recombination energy ($\alpha_{\rm rec}$) included in the calculation of the the binding energy parameter. On the other hand, for $\lambda=0$, the binding energy parameter is computed as a function of the WD progenitor (without any additional source of energy). For $\lambda<0$, the binding energy parameter is fixed and set equal to -$\lambda$. Therefore, using $\lambda=0.5$ as an input for the code, we were in fact assuming that 50% of the recombination energy contributes in the ejection process ($\alpha_{\rm rec}=0.5$).

Recent investigations have concluded that WD-MS binaries experience a strong orbital shrinkage during the CEP such that $\alpha$ should be $\sim 0.2-0.3$ without any additional source of energy, like recombination energy, or at least with a very small contribution ($\alpha_{\rm rec}\sim 0$), Zorotovic et al., 2010; Toonen & Nelemans, 2013; Camacho et al., 2014).

With the aim of exploring the influence of more realistic values of the CEP parameters, namely $\lambda$ and $\alpha$, we extend the set of models studied by Belloni et al. (2016, 2017). In this paper we simulate six new models, with exactly the same initial conditions described in Table 1, but with more reasonable choices for the CEP parameters, namely: $\alpha=1$ and $\lambda=0.0$ (i.e. $\alpha_{\rm rec}=0.0$). These assumed values imply that all of the variation in the binary orbital energy (of the giant core and the MS star) contributes to the expulsion of the CE, and that the binding energy parameter is variable and properly calculated assuming that none of the recombination energy is used to assist the expulsion of the CE. (Claeys et al., 2014, Appendix A).

To summarize, we consider in this work the six models simulated by MOCCA with the initial conditions given in Table 1, but assuming either $\alpha=3$ and $\alpha_{\rm rec}=0.5$ or $\alpha=1$ and $\alpha_{\rm rec}=0.0$. Hereafter, we will indicate the differences among the models K1, K2, K3, S1, S2, and S3, by including a subindex to indicate the assumed value of $\alpha$.

The choice of exploring different assumptions for the common envelope ejection is new to this work, and will allow us to re-assess previous results (Belloni et al., 2016, 2017) and their sensitivity to our assumptions for the CEP parameters, in addition to any differences caused by different AML approaches (Sections 2.2 and 2.3).

After having described the main codes and assumptions involved in this work as well as the 12 models analyzed, we turn our attention to our main results.

4 CLUSTER EVOLUTION AND CEP PARAMETERS

We initially calculated six models that differ mainly with respect to their IBPs, initial binary fractions and initial densities. Additionally, for this work, we simulated the same initial models, but with more realistic values for the CEP parameters. Hence, we are able to compare the global evolution of our simulated clusters when each of these two different sets of CEP parameters is adopted.

All together, we have a total of 12 models. Among these models, we find cases with black hole subsystems that form, post-core collapse clusters, intermediate-mass black holes and a cluster with a relaxation time greater than a Hubble time (initially extremely sparse).

PCEBs have shorter periods for lower values of $\alpha$, such that we find that the number of PCEBs produced in our simulated clusters is smaller for lower values of $\alpha$, since many of them merge during the CEP. A more general conclusion is: the lower is the value of $\alpha$, the longer is the initial minimum period of any PCEB progenitors (see Zorotovic et al., 2014b, for more detail). This is because any PCEBs that would have had an initial period shorter than the minimum period for PCEBs to merge during the CEP. Since the initial period is unchanged by the choice of the CEP parameters, decreasing $\alpha$ causes the number of PCEBs to decrease, thus increasing the number of mergers during the CEP.

As the majority of PCEBs in each model are WD-MS and WD-WD binaries (i.e. CV and AM CVn progenitors), which is a consequence of the adopted IMF, the number of WDs in the system is higher for clusters simulated with $\alpha=1$ and $\alpha_{\rm rec}=0.0$. During the merger process, the MS star is absorbed into the giant envelope and the outcome is a new giant. This giant will further evolve into a WD, with its type depending on the giant involved in the process. Similarly (albeit more complex), when the secondary is a WD, the outcome of the coalescence process is a giant that will likely evolve into a more massive WD (for more details on the coalescence of common-envelope cores, see Hurley et al., 2002, Section 2.7.2). Thus, instead of having the same amount of WD-MS and WD-WD binaries as in models evolved with $\alpha=3$ and $\alpha_{\rm rec}=0.5$, clusters evolved with $\alpha=1$ and $\alpha_{\rm rec}=0.0$ produce more single isolated WDs (i.e. less WD-MS and WD-WD binaries).

Clusters simulated with $\alpha=1$ and $\alpha_{\rm rec}=0.0$ produce more WD-NS binaries. This is easy to understand given how the formation of these objects might be influenced by the CEP parameters. When two WD progenitors are involved in two CEPs (one for each of the WDs), the final outcome of both CEPs (WD-WD binary) has a smaller period for $\alpha=1$ and $\alpha_{\rm rec}=0.0$ than for $\alpha=3$ and $\alpha_{\rm rec}=0.5$ (provided, of course, the binary survives the CEP). In some cases, the WDs are so close that mass transfer can proceed from a less massive He WD or C/O WD onto a more massive O/Ne/Mg WD. This can lead to the formation of a NS via accretion-induced collapse. Indeed, if the more massive WD triggers electron-capture reactions at its centre upon reaching the Chandrasekhar mass limit, then the burning which propagates outward makes central temperatures and pressures high enough to cause a collapse of the WD, inducing a contraction of the remnant (e.g. Nomoto & Kondo, 1991). These systems are interesting because they are potential progenitors of low-mass X-ray binaries and millisecond pulsars in GCs. For example, Rivera-Sandoval et al. (2015) recently discovered near-ultraviolet counterparts to five millisecond pulsars in the GC 47 Tucanae, and all of them have likely He WD companions, whose masses are $\sim$ $0.16-0.3$ M⊙ with cooling ages that lie between $\sim$ $0.1-6$ Gyr.

Despite the differences outlined above, the global cluster evolution is not affected by our choice for the CEP parameters. The time evolution of the main cluster properties remains about the same in all cases, including the rate of mass loss, the evolution of the core and half-mass radii, the rates and types of dynamical interactions, the binding energies of the resulting binaries, the binary fraction, etc. This result is expected since most PCEBs are too hard to interact on sufficiently short timescales to provide a sufficient additional source of energy to change the cluster evolution (Leigh et al., 2016).

5 WD-MS AND CV FORMATION RATES

In a previous paper (Belloni et al., 2017), we inferred that the formation rate of WD-MS binaries that will later become present-day CVs in globular clusters is strongly dependent on WD formation (due to stellar evolution). This is even the case for those WD-MS binaries formed through dynamical interactions. Additionally, given the time-delay between the WD-MS binary formation and CV formation, we inferred that the overall CV formation rate in globular clusters could be nearly constant. However, we emphasized in that work that this conclusion was based on a very small set of models.

In this section, we extend our previous analysis. Here, we bring to bear twice the number of models, each of which has more CVs relative to the old models. We again apply a two-sample Kolmogorov-Smirnov test to every pairwise combination of models. The alternative hypothesis of the test is that the paired samples are not drawn from the same theoretical CV formation rate, a null result that would support the conclusion of a model-independent CV formation rate. That is, all CV populations form in a similar way, independent of the host cluster properties.

Upon comparing models with the highest numbers of CVs (S2α=1 and K2α=1), we reject the null hypothesis at the 99.9 per cent confidence level. This suggests that the CV formation rate differs from model to model. This result is based on much better statistics than in Belloni et al. (2017). Thus, in our new sample, we are able to rule out the possibility that the CV formation rate in globular clusters is similar irrespective of the host cluster properties.

In Fig. 1, we show the cumulative distributions for the formation-age CV populations, for models S2α=1 and K2α=1. We consider all three schemes for the AML/CV evolution (MOCCA, cCAML, and eCAML). Upon comparing the different models for each AML/CV evolution scheme (same colour), we find that they are not the same.

As expected, for all models, when the eCAML formulation is adopted, a large fraction of CVs form early on in the cluster evolution ($\lesssim$ 2 Gyr). This is associated with the number of massive WDs in the population (formed at $\lesssim$ 2 Gyr)111The cluster turn-off mass at low metallicity adopted in the simulations is $\sim$ 1.5 M⊙ at $\sim$ 2 Gyr. MS stars whose masses are smaller than this are likely to become He WDs, i.e. low-mass WDs, which are not present using the eCAML scheme. and the short time-delay between WD-MS binary formation and CV formation ($<1$ Gyr). Additionally, most CVs that form early on in the cluster evolution belong to the BSE group (i.e., formed without any influence from dynamics).

We emphasize that additional correlations between the CV formation rate and the initial and present-day cluster properties might still exist. More models will be needed to address this fully. Although it is clear from our KS-tests that the CV formation rate differs from cluster to cluster, CVs are continuously being added to the population throughout the cluster evolution (Belloni et al., 2017). This is consistent with a nearly constant formation rate (in each cluster) after $\sim 1-2$ Gyr, with the “constant of proportionality” being different from cluster to cluster.

Finally, we stress that looking for theoretical correlations between CV formation rates and cluster properties is extremely important. For example, knowing this and knowing the number of bright CVs observed in a particular GC would allow us to account for the underlying observational biases and estimate the real number of CVs in such a GC, or even to use this information as a constraint for the initial GC conditions (e.g. Leigh et al., 2013, 2015).

6 NUMBER OF PRESENT-DAY CVS

In Table 2, we show the number of present-day CVs in all simulated clusters, separated according to the influence of dynamics for each formation channel. We emphasize that a few CVs formed adopting the MOCCA/cCAML/eCAML formulations for the WD-MS evolution do not survive in the cluster up to 12 Gyr. This is due to dynamical interactions or escape from the cluster or, mainly, mergers during unstable mass transfer. For the numbers shown in Table 2, only those CVs that form and survive to the present-day are included. More details on those CVs that do not survive may be found in Belloni et al. (2017, section 3.8).

In what follows, we discuss how CV production is affected by our choices for the CEP parameters (Section 6.1), AML schemes (Section 6.2), and initial cluster conditions (Section 6.3).

6.1 Influence of CEP parameters

In Section 4, we found that the number of PCEBs decreases when $\alpha$ and $\alpha_{\rm rec}$ decrease. Even though the number of PCEBs is reduced in our models with $\alpha=1$ and $\alpha_{\rm rec}=0.0$, this affects the number of CVs in the opposite way. In other words, we do not see a reduction in the number of CVs caused by the reduction in the number of PCEBs that are WD-MS binaries. Instead, the opposite occurs, i.e. the total number of CVs produced is larger for $\alpha=1$ and $\alpha_{\rm rec}=0.0$ in all models but model S3 (see Section 6.3).

This can be understood by looking at the final configuration of PCEBs. Lower values for the CEP parameters lead to shorter periods. This causes PCEBs to remain in the PCEB phase for shorter durations. In turn, this produces more CVs since the pre-CV lifetime is reduced and more systems become CVs in less than 12 Gyrs.

Regardless of the assumed AML and CV evolution prescriptions, there are no CVs produced in the Kroupa models when $\alpha=3$ and $\alpha_{\rm rec}=0.5$ are adopted. Only by decreasing $\alpha$ and $\alpha_{\rm rec}$ can the Kroupa models form CVs without any influence from dynamics. As suggested by previous works based on adopting the Standard IBP, the assumed CEP parameters should be small. This seems to be the case for the Kroupa IBP as well.

6.2 Influence of AML schemes

With respect to the CV evolution scheme, the number of CVs formed in MOCCA is smaller than the number of CVs formed adopting the cCAML formulation. This is associated with the stability criterion for dynamical mass transfer. In MOCCA, the criterion for a CV to start unstable mass transfer in a dynamical timescale is fixed ($q>q_{0}=0.695$). In the case of cCAML, this critical mass ratio depends on the mass of the secondary and it can be as great as $\sim 1.2$. Thus, many CVs that survive when cCAML is adopted merge in the pure MOCCA simulations. In the case of eCAML, the critical mass ratio also depends on the secondary mass. For low-mass secondaries it is more strict than the fixed value from MOCCA and also than the variable limit for cCAML, but it can reach extremely large values for more massive companions (Schreiber et al., 2016).

Interestingly, the relative numbers of CVs produced adopting the eCAML formulation is greater in the Kroupa models, compared to adopting the cCAML formulation. This indicates that the Kroupa IBP produces more C/O WDs (high-mass) than the Standard IBP, especially in model K$1_{\rm\alpha=1}$. This cluster is extremely sparse initially (i.e., only a few CVs form dynamically). Even so, $\sim$ 60 per cent of the CVs formed with the cCAML scheme are also present when eCAML is assumed. The same is not true for model S$1_{\rm\alpha=1}$, which is sparse as well. In this cluster, only $\sim$ 17 per cent of the CVs formed with the cCAML approach are present with the eCAML approach.

More generally, as expected, the fraction of CVs that remain with the eCAML formulation is around 15-20 per cent for the S models (Standard IBP), as pointed out by Schreiber et al. (2016). On the other hand, for the K models (Kroupa IBP), this fraction is $\sim$ 50-60 per cent. As mentioned above, this indicates that upon assuming a Kroupa IBP the production of CVs with high-mass WDs is more efficient, regardless of the formation channel.

6.3 Influence of the initial cluster conditions

In general, for both IBPs, the denser the cluster is initially, the smaller the number of CVs formed through stellar evolution alone. In fact, upon comparing the columns associated with the BSE group in Table 2, for all schemes (MOCCA, CAML and eCAML), we see a reduction in the relative number of CVs formed in the BSE group. This indicates that, for all clusters, any reduction in the relative numbers of CV progenitors is correlated with the initial cluster density, which is associated with the role of dynamical interactions in destroying CV progenitors, which is in turn related to the cluster soft-hard boundary.

The semi-major axis that defines the boundary between soft and hard binaries is proportional to $r_{h}/N$, where $r_{h}$ is the half-mass radius and $N$ is the number of objects (single + binaries) in the cluster. In addition, for clusters with similar $N$ (as is the case in the three Kroupa models), the denser the cluster (i.e. the smaller the half-mass radius), the smaller the semi-major axis that defines the soft-hard boundary. Therefore, the soft-hard boundary is inversely related to the cluster density, and as the cluster density increases, the soft-hard boundary also decreases (in semi-major axis/period). At a particular density, the soft-hard boundary will enter in to the region occupied by CV progenitors. Thus, beyond this density, more and more CV progenitors are potentially destroyed, as the density increases. This is the reason behind the above-mentioned correlation between the cluster density and the number of CV progenitors that become CVs via binary stellar evolution.

As already mentioned, the number of CVs in model S3 is smaller for lower values of the CEP parameters. This cluster is very dense and dynamics plays a huge role, due to the short orbital separation corresponding to the soft-hard boundary. For the lower values of $\alpha$ and $\alpha_{\rm rec}$, the number of WD-MS binaries in the cluster is about half the analogous numbers when the cluster is evolved with the larger CEP parameters. This is not the case for the S1 and S2 models where dynamical effects are not as strong. In these two models, for the lower CEP efficiencies, around 80% of WD-MS binaries survive the CEP in comparison to the number of WD-MS binaries obtained with the higher CEP efficiencies. In model S3, we adopt the smaller CEP efficiencies, such that a large number of CV progenitors come from initial binaries whose orbital separations are greater than the orbital separation that defines the soft-hard boundary. For these ‘wide’ binaries, the probability of interaction is larger and they can be more easily destroyed by dynamics.

Most CVs in the Standard models are formed via a CEP, which may be preceded by weak dynamical interactions. Usually, strong dynamical interactions do not play a significant role, which is directly related to the fact the Standard IBP contains more hard binaries than the Kroupa IBP (see appendix A). However, the influence of strong dynamical interactions increases slightly for the eCAML formulation. This is easy to understand based on the role of exchanges in the CV formation process, which puts more massive WDs into the CV population (e.g. Belloni et al., 2017). Thus, the probability of surviving the stability criterion imposed by the eCAML scheme is enhanced for CVs formed via strong dynamical interactions.

As already pointed out by Belloni et al. (2016, 2017), when $\alpha=3$ and $\alpha_{\rm rec}=0.5$, models that follow the Kroupa IBP produce only dynamically formed CVs, with exchanges being the main formation channel. Here we confirm that the same is obtained when CAML is included. When lower values are adopted for the CEP parameters ($\alpha=1$ and $\alpha_{\rm rec}=0.0$), the Kroupa models also have CVs formed through only stellar evolution (BSE group). Although field-like CVs can form in all Kroupa models with lower CEP efficiencies, only in the sparsest cluster (K1α=1) do they contribute significantly to the overall CV population. This is related to the fact that dynamically formed CVs are rare in sparse clusters.

7 PRESENT-DAY CV PROPERTIES

In this section, for convenience, we show only the results for models K2 and S2, for both choices of the CEP parameters. This is because these two models produce the largest numbers of CVs, which increases the statistical significance of our analysis. This is a reasonable approach since models with the same IBP tend to have similar properties in the end. However, the relative numbers of CVs formed from each channel depends on the initial cluster conditions, as we saw in Section 6.3.

In Sections 7.1 and 7.2, we discuss the CV component massses and periods. Subsequently, we address the mass transfer rate (Section 7.3), duty cycle (Section 7.4), X-ray luminosity (Section 7.5), quiescent magnitudes (Section 7.6) and spatial distribution (Section 7.7).

7.1 WD mass

The main difference in the WD mass distribution (left-hand column of Fig. 2) for the different AML formulations is the absence of low-mass WDs when eCAML is assumed. As already stated in Section 2.3, the empirical formulation for CAML is motivated by the absence of CVs with He WD primaries (Schreiber et al., 2016). Consequently, it is not a surprise that the WD mass distributions that use the eCAML scheme contain only C/O and O/Ne/Mg WDs (i.e. massive WDs). The minimum WD mass in the distribution for eCAML is $\sim$ 0.6 M⊙ (with a median value of $\sim$ 0.73 M⊙) independent of the host cluster properties, CEP parameters, formation channel or IBP. This is illustrated in Fig. 2 (left-hand column) which shows the WD mass distribution for models K2 and S2, at the moment of CV formation, i.e. at the onset of mass transfer. For the MOCCA and cCAML schemes, on the other hand, the WD mass can reach values as low as $\sim$ 0.2 M⊙.

The WD mass distribution is also affected by the IBP and CEP parameters. For the Kroupa models, the distributions have a strong peak at $\sim$ 0.7 M⊙. In the Standard models, another peak is also present at lower masses ($\sim$ 0.3 M⊙) which is absent in the eCAML approach. The peak at $\sim$ 0.7 M⊙ broadens and moves towards lower values ($\sim$ 0.6 M⊙) when smaller efficiencies for the CEP are assumed. The interquartile range of the WD mass distribution for the Standard IBP models (including all formation channels) is $[\sim 0.35;\sim 0.65]$. This indicates that the majority of the CVs formed in these models have low-mass WDs. On the other hand, upon considering the Kroupa IBP models, again for all formation channels, the interquartile range is $[\sim 0.61;\sim 0.78]$. This indicates that most CVs formed in these models have high-mass WDs. This is the main reason why the Kroupa models are less affected by the eCAML approach than are the Standard models. The differences in the interquartile range are due to the intrinsic differences between both IBPs adopted here in the region of parameter space corresponding to the CV progenitors (see Appendix A). The Kroupa IBP produces CVs with properties that more closely resemble those CVs observed in the Galactic field. But, with that said, a more thorough investigation that considers a realistic Galactic star-formation rate and observational selection effects is still to come.

CVs strongly influenced by dynamics (SDI group), mainly have massive WDs ($\sim$ 75 per cent of the CVs in the SDI group have masses greater than $\sim$ 0.6 M⊙). This is due to dynamical exchanges in which a binary composed of low-mass MS stars has one of its components replaced by either a more massive MS star or a more massive WD (Belloni et al., 2017). In both cases, the resulting CV WD mass is higher than it would otherwise have been without the exchange.

7.2 Donor mass and period

For the donor masses, we find an accumulation of systems at the very low-mass end (below the brown dwarf limit) in all the simulations, regardless of the IBP, AML prescription or CEP parameters. CVs simulated with MOCCA have in general higher secondary masses than with cCAML/eCAML. This is mainly associated with the donor mass-radius relation adopted in each code (see Section 8.3). Therefore, most CVs are period-bouncers, irrespective of the cluster properties, CEP parameters, formation channel or IBP. In general, CVs formed with a strong influence from dynamics (SDI group) have a lower fraction of period-bouncers, which implies that they have a larger probability of being observed during quiescence (due to higher mass transfer rates). However, for the Standard models, when considering either the cCAML or eCAML approaches with $\alpha=3$ and $\alpha_{\rm rec}=0.5$, the opposite occurs, i.e., for these models, the fraction of period-bouncers is higher in the set of dynamically formed CVs. This is because most CVs in the SDI group are formed at earlier times.

The period distributions are shown in the right-hand column of Fig. 2. Here we see a drastic difference in the CV evolution model implemented by MOCCA relative to the state-of-the-art model implemented in the cCAML/eCAML scheme. This is mainly due to the mass-radius relation adopted in BSE (Hurley et al., 2000, see Eq. 24), that assumes a substantial increase in the radius of the donor when its mass decreases, for secondary masses below $\sim$ 0.1 M⊙. This results in period bouncers that can reach periods longer than 3 h, which is not realistic (Knigge et al., 2011).

On the other hand, in the state-of-the-art CV evolution model, the mass-radius relation adopted (Knigge et al., 2011) assumes that below $\sim$ 0.07 M⊙ the radius decreases more slowly as the mass decreases, and at no point does it increase. Therefore, period bouncers remain close to the period minimum for longer. However, the authors caution that this mass-radius relation should not be extrapolated below $\sim$ 0.05 M⊙, although the systems should remain close to the period minimum and with very low mass transfer rates (probably undetectable).

7.3 Mass transfer rate

For the mass transfer rate distribution (left-hand column of Fig. 3), the main differences between models are associated with the CV evolution scheme and the different formation channels.

Although all schemes are concentrated towards low mass transfer rates, the peak in the distribution for systems evolved with the MOCCA scheme is at higher values than for those that incorporate CAML. As already discussed in Sect. 7.2, this is related to the mass radius relation adopted in the BSE code used in MOCCA as well as an inadequate mass transfer rate equation for the AML timescale.

For CVs formed without any or with only a weak influence from dynamics (BSE and WDI groups), most CVs have mass transfer rates smaller than $\sim 5\times 10^{-12}$ M⊙/yr, especially in the Kroupa models, for both sets of CEP parameters. For those CVs strongly influence by dynamics (SDI group), the distribution is flatter. This indicates that CVs in the SDI group should be brighter, in general, than CVs in either the BSE or WDI groups. As we will see in Section 7.6, this is in fact the case.

Applying the DNe criterion, as described in Belloni et al. (2016, see Section 3.3.1), we can separate the CVs into three groups with respect to the instability of their disks, namely CVs whose disks are hot/stable, unstable, or cold/stable.

For most CVs in our models, the disks are unstable, with very few exceptions of stable disks (less than 1 per cent). Additionally, most of the disks exhibit outbursts of Type B (i.e. inside-out) that begin in the inner parts of the disk and propagate outward (e.g. Smak, 2001; Lasota, 2001). This is due to their extremely small mass transfer rates. Thus, if such CVs do not have WDs with relatively high magnetic fields such that the inner part of their disks might be disrupted (Dobrotka et al., 2006), practically all CVs in our models are DNe and undergo outbursts. Most GC CVs seem to be DNe (e.g. Servillat et al., 2011; Knigge et al., 2011; Knigge, 2012; Webb & Servillat, 2013; Belloni et al., 2016). As we will see in what follows, the occurrence of outbursts for such CVs is rather infrequent, which makes their detection through outbursts very difficult.

7.4 Duty cycle

One important property of DN CVs with regard to their detectability is the duty cycle. This can be defined as the fraction of the DN cycle in which the DN is in outburst, i.e. the ratio between the duration of the outburst and the recurrence time. In order to compute these two quantities, we use empirical relations as described in Belloni et al. (2016, see Section 3.3.5).

7.4.1 Recurrence time

Due to the differences in the donor mass and period distributions previously discussed, the CV recurrence times with the MOCCA formulation are smaller than those for CVs formed with the cCAML and eCAML formulations. This is because the recurrence time is inversely proportional to the mass ratio, and the mass ratios for MOCCA CVs are greater. Specifically, most CVs formed in MOCCA have recurrence times shorter than $10^{4}$ days, for both sets of CEP parameters.

With the cCAML formulation, we see a clear difference with respect to the IBP. When $\alpha=3$ and $\alpha_{\rm rec}=0.5$ are adopted, the median of the distribution is $\sim 10^{4}$ days for the Kroupa models and $\sim 10^{3}$ days for the Standard models. When $\alpha=1$ and $\alpha_{\rm rec}=0.0$ are adopted, the median of the distribution for the Kroupa models is similar, while for the Standard models the median increases and becomes similar to what is seen in the Kroupa models.

With the eCAML formulation, regardless of the assumed CEP parameters and the IBP, most CVs have recurrence times longer than $\sim 10^{5}$ days. They are longer here because the mass ratios are much smaller, since these CVs are dominated by massive WDs.

In general, CVs strongly influenced by dynamics (SDI group) have longer recurrence times. This is because, usually, they are formed from more massive WDs.

7.4.2 Outburst duration

For the duration of the outburst, which is associated with the period (the longer the period, the greater is the extent of the disk and the longer is the outburst), we see a clear distinction upon comparing the MOCCA and cCAML/eCAML formulations. As noted before, most CVs formed in the cCAML and eCAML simulations have similar periods, which is close to the period minimum (in turn, the outbursts have similar durations). On the other hand, CVs in MOCCA simulations have much longer periods (up to $\sim 3-4$ h), which makes the outbursts last longer. Even so, the duration of the outbursts is just a few days. This is the case for all models, CEP parameters and AML/CV evolution schemes.

We emphasize here that the empirical relation used to compute the duration of the outbursts was derived from well-observed CVs in the solar neighbourhood (Smak, 1999). It includes only the normal outburst. This is, apparently, not the case for period-bouncers, where superoutbursts play a more important role and could increase the duration time significantly, since the outbursts can last for months. We discuss how this affects our results in more detail in Section 8.6.

7.4.3 Duty cycle

Finally, we present our main results regarding the behaviour of two important timescales associated with DNe. Thereafter, we will turn our attention to a description of the duty cycle distribution, which is illustrated in the right-hand column of Fig. 3.

Since the duration time is similar for all models and formulations, the main difference in the duty cycle comes from the recurrence time. The duty cycles of CVs formed in the MOCCA simulations are usually larger than $\sim$ 0.1 per cent (first quartile), and smaller than $\sim$ 1 per cent (third quartile), for both choices of CEP parameters. On the other hand, for the cCAML and eCAML formulations, most CVs have duty cycles smaller than $\sim$ 0.1 per cent and $\sim$ 0.001 per cent, respectively. This implies that GC CVs are even more difficult to detect upon adopting more realistic and accurate CV evolution schemes. This corroborates our previous result (Belloni et al., 2016), where we concluded the same thing based on the results of pure MOCCA simulations, which are out-dated with respect to their CV evolution schemes.

As for the recurrence times, CVs strongly influenced by dynamics usually have smaller duty cycles. This is because they usually have more massive WDs and smaller mass ratios. However, we caution that those CVs in the SDI group that have low-mass WDs can still have large duty cycles.

7.5 X-ray luminosity during quiescence

One important observational property of GC CVs is their X-ray luminosity. It is computed here as described in Belloni et al. (2016, see Section 3.3.4). Most of the predicted CVs in our models have X-ray luminosities between $\sim 10^{29}$ and $\sim 10^{30}$ erg s*-1*. However, a few of them have $L_{X}$ below $10^{29}$ erg s*-1*. These luminosities for our period-bouncers seem to agree with the values found for GW Lib and WZ Sge, namely $0.05^{+0.10}_{-0.02}$ $\times 10^{30}$ and $0.7^{+0.3}_{-0.1}$ $\times 10^{30}$ erg s*-1* [2-10 keV], respectively (Byckling et al., 2010). These systems are likely period-bouncers or CVs close to the period minimum, which makes them ideal for this comparison.

Even though our results seem to agree with the observed X-ray fluxes for CVs close to the period minimum, recent investigations have found a correlation between the duty cycle and the X-ray luminosity for X-ray bright DNe in the solar neighborhood (Britt et al., 2015). Our simulated CVs show a different behaviour than this, most likely because the authors primarily have bright X-ray DNe in their sample ($L_{X}>3\times 10^{30}$ erg s*-1*), whereas our simulated CVs are mostly period-bouncers or CVs close to the period minimum, with smaller X-ray luminosities. This makes a more rigorous comparison difficult.

7.6 Magnitude during quiescence

Another important observational property of GC CVs is their magnitudes. Once a potential candidate is identified via its X-ray luminosity, a confirmation via other techniques (such as identifying the optical counterparts) are required in order to improve the confidence level that the source is in fact a CV candidate222As in the Galactic field, only spectroscopy can confirm that a CV candidate is indeed a CV (Knigge, 2012). However, since GCs are crowded fields, spectroscopy might be a challenge. The use of a combination of different techniques (H$\alpha$ and FUV imaging, X-ray, colour, late and negative superhumps, etc.) will therefore be needed to confirm the GC CV candidates, especially for DNe..

In general, for period-bouncers, the optical flux is dominated by the WD, with little (if any) contribution from the hot spot corresponding to where the stream of matter from the donor strikes the accretion disk. This is especially true for the population of faint CVs found in the GC NGC 6397 (Cohn et al., 2010). As claimed by these authors, the WD should be heated by the accretion process overall, since isolated massive WDs are faint enough to avoid detection (i.e., they have efficient cooling).

In our simulations, we estimate the CV magnitude as described in Belloni et al. (2016). No additional heating is added to the WDs. For the cCAML and eCAML schemes, we do not have information about the WD magnitude or temperature. Thus, we present our results only for the MOCCA formulation. Here, we have information about the WD magnitude and can compute the CV magnitude by adding the flux contributions from the four main components: WD, donor, hot spot and disk.

Present-day CV magnitudes are strongly dependent on the time since CV formation, as well as on the formation mechanism. Basically, the later the CV is formed (i.e., mass transfer starts), the brighter it is. Since dynamically formed CVs are more massive, those newly formed CVs are brighter than CVs newly formed from primordial binaries. Here, we define newly formed CVs as CVs formed after $\sim$ 10 Gyr of cluster evolution, i.e. at most 2 Gyr ago.

Among those CVs that form at intermediate times (after $\sim$ 1 Gyr of cluster evolution but before $\sim$ 10 Gyr), we see a transition in dynamically formed CVs; those formed a long time ago are currently fainter but those formed more recently are currently brighter than those CVs formed from primordial binaries. This transition occurs at $\sim$ 6 Gyr. The reason is twofold. On the one hand, this effect is partially due to the above-mentioned cooling efficiency of the WDs. But, on the other hand, the hot spot becomes the brightest CV component at this particular formation-time (and its luminosity is related to the donor mass, i.e. the higher the donor mass, the higher the mass transfer rate, and in turn the higher the hot spot luminosity). In other words, CVs formed at maximum $\sim$ 6 Gyr ago have their optical fluxes dominated by their hot spots and those formed before this time have their optical fluxes dominated by their WDs. This makes dynamical CVs more luminous (currently) if formed after $\sim$ 6 Gyr of cluster evolution and CVs formed from primordial binaries more luminous (currently) if formed before $\sim$ 6 Gyr of cluster evolution.

For CVs that formed a long time ago (before $\sim$ 1 Gyr of cluster evolution), distinguishing between the different formation channels is more difficult. This is because these CVs have more massive WDs irrespective of their formation channels, such that their WD cooling efficiencies and magnitudes are similar.

The description above is illustrated in Fig. 4, where only CVs formed using the MOCCA scheme are shown. These are separated according to our choices for the CEP parameters and our classifications for the influence of dynamics in CV formation. Note that the trends associated with each formation channel are independent of the CEP parameters.

We emphasize here that although we do not consider changes in the WD magnitude arising from the accretion process, we expect the same overall results were this effect to be included. This is because the above-described picture concentrates on the effects of dynamics in shaping present-day GC CV magnitudes and detection limits. This feature or transitional path is similar in all models, as originally observed in the GC NGC 6397 (Cohn et al., 2010).

7.7 Spatial distributions

Finally, we provide a few words regarding the CV spatial distribution, as a function of the CV brightness. Cohn et al. (2010) found that there is no strong evidence in favour of the radial distribution of main-sequence-turn-off stars being different from the faint CVs ($p$-value $\sim$ 0.04). However, for bright CVs, the evidence is sufficient ($p$-value $\sim$ 0.001) to make the claim that bright CVs are more centrally concentrated.

In our models, we also typically find this trend. However, it strongly depends on the following properties: (i) the source of energy in the host cluster, (ii) the host cluster evolution, (iii) the average mass in the host cluster core, (iv) the WD-MS binary and CV formation times, (v) the WD-MS binary and CV masses, and (vi) the formation channel (Belloni et al., 2016, see Section 5.7). We emphasize that the inclusion of the six additional models considered in this paper does not change the picture described in Belloni et al. (2016), in particular CVs are most likely to be found between the core and half-mass radii.

8 DISCUSSION

In the first paper of this series (Belloni et al., 2016), we presented six specific MOCCA models with a focus on the properties of their present-day CV populations. In the second paper in this series (Belloni et al., 2017), we concentrated instead on the properties of the progenitor and formation-age populations, and the age-dependence of CV properties. In this paper, we extend these analyses further by considering twice the number of models, including two possible combinations of the CEP parameters and two AML and CV evolution formulations. That is, our focus in this paper is on the influence of our adopted prescriptions for the binary evolution in determining our simulated CV properties, at different phases in the life of a CV (formation-age, present-day, etc.).

In what follows, we discuss the principal implications of our results with respect to the simulated GC CV properties. We also discuss the dependences of our results on specific assumptions and approaches.

8.1 Initial binary populations

In this work, as in our previous papers, we simulate clusters that follow two different IBPs (Section 3), namely the Kroupa and Standard IBPs. The main distributions associated with both IBPs are shown in Figures 6 and 7 in the Appendix A.

The Standard IBP is associated with an uniform mass ratio distribution, an uniform distribution in the logarithm of the semi-major axis $\log(a)$ or a log-normal semi major axis distribution, and a thermal eccentricity distribution. In GC investigations, the Standard IBP is usually adopted because it supplies hard binaries that work as an energy source in the cluster over the long-term evolution. We note that in population synthesis codes for Galactic studies, the Standard IBP is also adopted, even though observational results show slightly different sorts of distributions. Such observational features are obtained directly from the Kroupa IBP after some stimulated evolution (see below).

It is usually accepted that the Galactic field population formed from the dissolution of embedded star clusters (e.g. Lada & Lada, 2003) after the expulsion of the residual gas during the star formation process. In order to reconcile this picture with the observations of G and M-dwarf binaries in the Galactic field (Duquennoy & Mayor, 1991; Mazeh et al., 1992; Fischer & Marcy, 1992), Kroupa (1995) developed a theory in which star clusters experience some degree of dynamical processing, including disruption of binaries, before dissolution. This dynamical processing should ultimately reproduce the observed properties of G and M-dwarfs in the Galactic field. Kroupa (1995) also showed that a single dynamical operator corresponds to the mean embedded cluster from which the Galactic field G and M-dwarf binaries originate (Marks & Kroupa, 2012). The confirmation of this theory came later when data of all late-type binary systems near the Sun became available (Reid & Gizis, 1997).

The Kroupa IBP has been tested against both numerical simulations and observations, and has successfully explained the observational features of young clusters, associations, and even binaries in old GCs (e.g. Kroupa, 2011; Marks & Kroupa, 2012; Leigh et al., 2015, and references therein).

In our investigation, our results show good agreement with the observed Galactic CV WD mass distribution for the Kroupa models. However, the CV WD mass distribution is sensitive to AML prescriptions (Schreiber et al., 2016) and more complete investigations should be carried out taking this into account.

Collectively, all of this suggests that the Kroupa IBP might be a better choice to seed population synthesis codes (and GC simulations), compared to the Standard IBP. However, a more thorough investigation considering observational selection effects and realistic Galactic star-forming processes (e.g. spatial distribution and stimulated evolution) needs to be conducted, including different constraints on, for example, Galactic WD-MS PCEBs and CVs.

8.2 Influence of dynamics for different sets of initial cluster conditions

We saw in Section 6.3 that the dominant formation channel in a cluster depends on its initial conditions. In general, as expected, initially sparse clusters mainly form CVs similarly to the Galactic field, i.e. without any help from dynamics. On the other hand, when clusters are initially dense ($\rho_{c}\gtrsim 10^{4}$ M⊙ pc*-3*), there is always a balance between dynamical creation and destruction of CV progenitors. This balance is strongly related to the cluster soft-hard boundary and the adopted IBP, and clearly outlining this dynamical picture might be challenging. Two important correlations, however, are: the importance of the CEP decreases as the cluster density increases, which was also found by Ivanova et al. (2006). The second correlation is found between the efficiency of the destruction of primordial CV progenitors and the initial cluster density; the higher the initial density, the more primordial CV progenitors are destroyed.

Belloni et al. (2017) concluded that the CV formation rate could be similar in many clusters. As we saw in Section 5, we rule out this initial impression, which was an artifact of small-number statistics. From the new models, we found that the CV formation rate is not the same for all clusters, although CVs are continuously added to the population. After $\sim 1-2$ Gyr, the rate could still be nearly constant (with only the “constant of proportionality” differing in different clusters).

Unfortunately, due to the limited sample of models analyzed here, we cannot reliably identify correlations between the initial or present-day cluster properties and the simulated CV properties, including the formation scenario, fraction of bright CVs, etc. This can be extended to include possible correlations between the CV formation rate and the initial cluster density or any present-day cluster properties. In future work, we aim to identify more correlations, by analyzing a larger subset of simulations from the MOCCA-SURVEY database (Askar et al., 2017) with difference choices for the binary stellar evolution parameters.

Finally, for the Kroupa models, Belloni et al. (2017) found that the net effect of weak dynamical interactions on the WDI group is strong. This is because CVs in the WDI group had similar properties as those CVs strongly influenced of dynamics (SDI group). In this work, we find that the net effect is weak (as in the Standard models), since CVs in the WDI group have similar properties to those CVs in the BSE group (i.e., formed without any influence from dynamics). This initial impression came from the fact that no CVs were found in the BSE group in our previous work (i.e., with larger efficiencies for the CEP). This is not the case when more realistic values of the CEP parameters are adopted.

8.3 CV evolution scheme

Given that the CV evolution model implemented in the BSE code (and used in the MOCCA code) is outdated, we evolved separately all close WD-MS binaries formed in our simulations using the code described in (Zorotovic et al., 2016), which is much more up-to-date. In this section, we discuss possible future improvements that can be made to the BSE code for more realistic simulations.

Compared to more recent implementations, the main differences in the BSE prescriptions for CV evolution are:

the donor mass-radius relation adopted in BSE differs drastically from that derived from observations, especially in the regime of low-mass donors (M-dwarf and brown dwarf); 2. 2.

the normalization factors adopted in BSE for AML due to magnetic braking and gravitational radiation are different to those derived from the best-fit sequence by Knigge et al. (2011) adopted for cCAML and renormalized for eCAML in the code from Zorotovic et al. (2016); 3. 3.

the Roche-lobe radius definition used in BSE is that derived by Eggleton (1983), which represents the donor as a point source in the potential. In the state-of-the-art model, the donor Roche-lobe depends on the polytropic index of the donor; 4. 4.

consequential angular momentum loss (CAML) that is associated with mass transfer (e.g. due to nova eruptions) is not considered in BSE; 5. 5.

the predicted donor radius does not increase in BSE once mass transfer begins. However, according to the observations, CVs above the orbital period gap are $\sim$ 30% larger with respect to the radius they would have without mass transfer (for more details Knigge et al., 2011, see Section 5.2). The code from Zorotovic et al. (2016) uses the mass radius relation from Knigge et al. (2011) for CVs, which takes into account an expansion factor that allows us to reproduce the observed orbital period gap (which can not be reproduced with the BSE code). 6. 6.

the critical mass ratio above which dynamical mass transfer from a low-mass MS star occurs is fixed to 0.695 in BSE. A more accurate critical mass ratio estimate stems from equating the adiabatic mass-radius exponent to the donor’s Roche-lobe mass-radius exponent, which depends strongly on the adopted CAML scheme (see Schreiber et al., 2016, for more details).

In Fig. 5 we show both the donor mass-radius relations and the CV evolution in the donor mass-period plane. The red dots in the left-hand panel correspond to all present-day CVs simulated with BSE, while the dashed line is the donor mass-radius relation given by Eqs. 16 and 17 in Knigge et al. (2011) that is used for the cCAML/eCAML formulations. In the right-hand panel, we also included the donor mass versus orbital period relation for all present-day CVs simulated with the code from Zorotovic et al. (2016) (blue crosses). Notice how different the results from both models are, especially in the domain of extremely-low-mass donors ($\lesssim 0.1$ M⊙). In the BSE code, CVs evolve rapidly towards longer periods after reaching the minimum period, which is not consistent with the observations.

Ultimately, all of the above differences should be corrected in BSE, which would allow for more realistic predictions for the subsequent CV evolution. Upgrading the BSE code would be a service to the stellar dynamics community, since the BSE code is used in many star cluster evolution codes.

8.4 Common-envelope phase parameters

In this work, we simulate clusters with two choices for the CEP parameters (Section 3), namely $\alpha$ = 3.0 and $\alpha_{\rm rec}$ = 0.5 and $\alpha$ = 1.0 and $\alpha_{\rm rec}$ = 0.0.

One important, albeit alarming, result found by Belloni et al. (2016) is that models following the Kroupa IBP do not produce CVs formed through isolated binary evolution, when a high CEP efficiency is adopted. In this work, we show that for more realistic values of the CEP parameters, this apparent problem is resolved. This is, in turn, an indication that assuming high values for the CEP parameters combined with a Kroupa IBP is an inadequate combination. This is consistent with studies performed to date using the Standard IBP (Zorotovic et al., 2010; Toonen & Nelemans, 2013; Camacho et al., 2014).

An important effect that arises by decreasing the CEP parameters is that the number of CVs increases, even though the number of PCEBs decreases. The decrease in the number of PCEBs is caused by the enhancement of mergers during the CEP. On the other hand, as the reduction in orbital energy is stronger for lower CEP efficiencies, systems that survive the CEP with a lower efficiency emerge from this phase at shorter periods, i.e, they need less time to start mass transfer and become CVs. This is seen in both IBPs adopted here.

Finally, we note an increase in the number of potential millisecond pulsars with He WD companions. Interestingly, recent searches for near-ultraviolet counterparts to millisecond pulsars in the GC 47 Tuc have identified five such systems to date (Rivera-Sandoval et al., 2015).

8.5 Consequential angular momentum loss treatment

In order to quantify the effects of CAML, which is not included in the BSE code, we consider in this investigation two formulations for CAML, namely the classical CAML (cCAML) from King & Kolb (1995) and the empirical CAML (eCAML) from Schreiber et al. (2016), as described in Section 2.3.

The main effects on our results stem from the stability criterion for dynamical mass transfer, which strongly affects the numbers and properties of our simulated CVs.

Let us first compare the criterion in the BSE code with the cCAML approach. As explained by Schreiber et al. (2016), the stability criterion can be translated into a critical mass ratio $q_{0}$ above which mass transfer is unstable (usually, dynamically or thermally). In the BSE code, the critical mass ratio for dynamically stable mass transfer is fixed to $q_{0}=0.695$ for low-mass MS stars, while the models that incorporate CAML derive a value for $q_{0}$ that depends on the donor mass. When cCAML is assumed, $q_{0}\sim 0.91$ for donor masses below $\sim$ 0.4 M⊙. Thereafter, it is a monotonically increasing function of the donor mass. However, although the cCAML formulation produces more CVs than the MOCCA models, we see only minor differences upon comparing the WD mass distributions.

Next, we compare the cCAML and eCAML formulations. In the latter, $q_{0}$ is strongly affected by the donor mass; it is a monotonically increasing function of the donor mass for all mass ranges (Schreiber et al., 2016, see their Fig. 2). Due to the way it was developed, eCAML removes all CVs whose WDs have He cores from the population due to dynamical mass transfer (which leads to merges). The Kroupa models do not show strong variations upon adopting either the cCAML or eCAML schemes, except for the total absence of low-mass WDs and the smaller number of CVs with the eCAML scheme. This is because most CVs formed in the Kroupa models have C/O WDs (with a small contribution from O/Ne/Mg WDs). These CVs are (presumably) weakly affected by enhanced AML as is provided by the eCAML scheme. For the Standard models, on the other hand, the WD mass distribution is strongly affected by the eCAML approach.

Based on our results, we emphasize that, in order to have better predictions for GC CVs for comparisons to the observations, the more realistic cCAML or eCAML formulations should be used, along with more accurate CEP parameters. This should help to avoid neglecting the formation of particular CVs that might contribute significantly to the global population, with the overall effects depending on the choice of IBP.

8.6 Empirical relations and assumptions

Throughout this work, we have relied on empirical relations to compute important CV properties, such as the duration of the outburst and its recurrence time. The most important such assumption is the donor mass-radius relation used in the cCAML and eCAML formulations. This is precisely what is shown in Fig. 5. This relation was derived from observations of CV donors whose masses are $\gtrsim$ 0.05 M⊙. We saw in Section 7.1 that the majority of CVs formed using the cCAML and eCAML schemes have donor masses below this value, independent of the formation channel. This is, in principle, expected since GC CVs are much older than CVs in the solar neighborhood, and they have had enough time to evolve far beyond the period minimum, at which point $M_{2}(P_{\rm min})\sim 0.07$ M⊙. On the other hand, if the donor mass-radius relation for extremely-low-mass BDs is somehow different, many CVs in our sample might have extremely low mass transfer rates, which would make their detection next to impossible with current instrumentation.

The effects on the overall CV population due to extrapolations of the recurrence time (Patterson, 2011) and outburst duration (Smak, 1999) empirical relations should also be explored. Another way that might potentially change CV population properties is the extrapolation of the empirical relations for recurrence times and durations of outbursts, in to the range of extremely-low-mass donors. These expressions were derived from well-observed CVs in the Galaxy and could be different for CVs far away from the period minimum. These relations are important because they are used here to estimate CV duty cycles, which might be different if the empirical relations are different.

Most period-bouncers or period-bouncer progenitors undergo superoutbursts, with just a few (if any) normal outbursts occurring during supercycles. Consequently, we might be missing important information that could lead to changes in our duty cycle estimates for period-bouncers.

8.7 Comparisons to CV candidates in GCs

Although rigorous comparisons between our models and real GC CVs are difficult, due to strong degeneracies associated with the initial cluster conditions (Askar et al., 2017), general conclusions can still be drawn. These are listed below.

We showed in Sections 7.6 and 7.7 that bright CVs were mainly formed a few Gyr ago due to strong dynamical interactions (mainly exchanges). They tend to be more centrally concentrated, i.e. close to, although outside, the core. This is in agreement with what Cohn et al. (2010) found for their bright CVs. They suggested that CVs are born in (or close to) the core and then migrate out of the core due to dynamical interactions, as they age. In our simulations, we see that dynamical interactions are unlikely after CV formation (Belloni et al., 2017; Leigh et al., 2016). Of those CVs that form in (or close to) the core and are kicked out to the outskirts, they quickly sink back in to the inner parts due to mass segregation.

This was also found by Hong et al. (2017), who investigated 81 simulated clusters with different initial masses, half-mass radii, Galactocentric distances, and primordial binary fractions. The authors concluded that bright CVs (donor masses greater than 0.1 M⊙) are more centrally concentrated than faint CVs. Additionally, the faint CVs have similar radial distributions relative to the main-sequence turn-off stars.

More generally, bright CVs are formed close to the core, but are often kicked out of the core in the process, before finally migrating back to the core due to mass segregation. Faint CVs can form (but not always) through strong dynamical interactions, close (or not) to the core. However, these are highly evolved systems (much older than bright CVs) and, since they have masses that are similar to (or even smaller than) main-sequence turn-off stars, they should not be as centrally concentrated as bright CVs. What is important from an observational point of view is that most CVs (bright or faint) should be located somewhere between the core and half-mass radii.

Our results suggest that detection through outbursts is possible for only a rather small fraction of the total CV population. This is because their duty cycles are usually extremely small, especially using a more realistic CAML treatment and CV evolution model. It follows that the apparent lack of outbursts in GC CVs could just be a selection effect, in the sense that our knowledge is limited to a small population of very close CVs that are frequently observed in outburst, since these are the easiest systems to detect (e.g. Servillat et al., 2011; Knigge, 2012; Belloni et al., 2016).

Finally, as has already been pointed out by Knigge (2012), the natural path toward improving our understanding of GC CVs is a deep survey for DNe in GCs, which would guarantee the detection of at least a few WZ Sge systems. This would allow for a much more thorough comparison between the predictions of theory and observations, a crucial step toward disentangling the true nature of GC CVs.

9 CONCLUSIONS AND PERSPECTIVES

In terms of the WD mass distribution, we found that models that follow the Kroupa IBP show better agreement with the observed Galactic CVs than the those that assume the Standard IBP.

Contrary to what we concluded in Belloni et al. (2017), the CV formation rate seems to depend on our assumptions, especially during the first 1-2 Gyr of cluster evolution. Also, CVs formed through weak interactions (WDI group) seem to have similar properties to the ones formed without any interaction (BSE group), contrary to those CVs formed by strong interactions (SDI group).

Upon assuming more realistic (lower) values for the common envelope efficiencies, we found that more CVs (and less PCEBs) are produced, especially in the models that follow the Kroupa IBP. Including the empirical approach for CAML from Schreiber et al. (2016), CVs with low-mass (helium-core) WDs are not produced, which is consistent with the observations.

Despite all the uncertainties involved in simulating CVs in GCs, we can infer from our simulations that bright CVs in GCs are young and mainly formed due to exchanges. These tend to be concentrated towards the centre of the cluster, but outside the core, and our simulations show that strong dynamical interactions become rare after the CV is formed.

Due to the small duty cycles, especially when more realistic models for CV evolution are used, our simulations suggest that only a small fraction of the total CV population should be detectable through outbursts. This would explain the apparent lack of outbursts in GC CVs.

Future investigations will concentrate on the influence of the CEP parameters as well as the CAML formalisms on populations of WD-MS PCEBs, as well as CVs formed from the Kroupa IBP. For that purpose, we intend to carry out a detailed population synthesis study incorporating a realistic Galactic star-formation model (temporal and spatial) along with observational selection effects intended to further constrain stellar evolution parameters by direct comparisons to the observed properties of Galactic WD-MS PCEBs and CVs. Additionally, we aim to analyse more GC models within the MOCCA-SURVEY to look for additional correlations between GC and CV properties.

Acknowledgements

DB was supported by the CAPES foundation, Brazilian Ministry of Education through the grant BEX 13514/13-0 and by the National Science Centre through the grant UMO-2016/21/N/ST9/02938. MZ acknowledges financial support from FONDECYT (3130559). MRS thanks for support from FONDECYT (1141269). MG and AA were supported by the National Science Centre through the grant DEC-2012/07/B/ST9/04412. AA would also like to acknowledge support from the National Science Centre through the grant UMO-2015/17/N/ST9/02573 and partial support from Nicolaus Copernicus Astronomical Centre’s grant for young researchers.

Appendix A Initial Binary Populations

In this investigation, we simulated models that follow two distinct initial binary populations (IBPs), namely: the Kroupa IBP and the Standard IBP (see Section 3 for more details). In order to provide to the reader an easy way to recognize the differences between them, in Fig. 6 we plot the main distributions associated with both IBPs for all initial binaries, i.e. primary mass (top left-hand panel), mass ratio (top right-hand panel), period (bottom left-hand panel), and eccentricity (bottom right-hand panel).

Additionally, since most CVs formed via pure binary stellar evolution come from initial binaries with primary masses between $\sim 0.8$ and $\sim 7.0$ M⊙, periods shorter than $\sim 10^{5}$ d, and mass ratios smaller than 0.5, we plot in Fig. 7 the same distributions in Fig. 6, but for the above-mentioned regions in the parameter space. This allows us to compare both the IBP for the entire population with the IBP of potential CV progenitors.