Initial Conditions

Running CMC comes in two distinct stages: generating the initial conditions, and then running CMC on those initial conditions. We briefly cover here how to create initial conditions and how to run and restart CMC.

Point-mass Initial Conditions

Generating initial conditions is the first step in running CMC. Since Version 3.4, COSMIC has been able to create CMC initial conditions, using a combination of the native binary samplers and dynamical samplers for drawing positions and velocities of clusters in virial and hydrostatic equilibrium. COSMIC can currently generate clusters from Plummer, Elson, or King profiles, as detailed below.

The examples below are also found in the examples folder of the main GitHub repository. Note that there are two distinct options for CMC clusters: cmc_point_mass (which we use in the Plummer and Elson examples below) will produce a cluster of point-mass particles: no stellar radii, mass functions, or binaries, typical for purely dynamics problems. You can see it’s options with:

In [1]: from cosmic.sample.sampler import cmc

In [2]: help(cmc.get_cmc_point_mass_sampler)
Help on function get_cmc_point_mass_sampler in module cosmic.sample.sampler.cmc:

get_cmc_point_mass_sampler(cluster_profile, size, **kwargs)
    Generates an CMC cluster model according  
    
    Parameters
    ----------
    cluster_profile : `str`
        Model to use for the cluster profile (i.e. sampling of the placement of objects in the cluster and their velocity within the cluster)
        Options include:
        'plummer' : Standard Plummer sphere.
            Additional parameters:
            'r_max' : `float`
                the maximum radius (in virial radii) to sample the clsuter
        'elson' : EFF 1987 profile.  Generalization of Plummer that better fits young massive clusters
            Additional parameters:
            'gamma' : `float`
                steepness paramter for Elson profile; note that gamma=4 is same is Plummer
            'r_max' : `float`
                the maximum radius (in virial radii) to sample the clsuter
        'king' : King profile 
            'w_0' : `float`
                King concentration parameter
            'r_max' : `float`
                the maximum radius (in virial radii) to sample the clsuter
    
    size : `int`
        Size of the population to sample
    
    Returns
    -------
    Singles: `pandas.DataFrame`
        DataFrame of Single objects in the format of the InitialCMCTable
    Binaries: `pandas.DataFrame`
        DataFrame of Single objects in the format of the InitialCMCTable

The more general cmc sampler (which we use in the King profile example below) has all the additional features for generating stellar binaries that COSMIC has. You can see it’s (significantly more diverse) options with

In [3]: from cosmic.sample.sampler import cmc

In [4]: help(cmc.get_cmc_sampler)
Help on function get_cmc_sampler in module cosmic.sample.sampler.cmc:

get_cmc_sampler(cluster_profile, primary_model, ecc_model, porb_model, binfrac_model, met, size, **kwargs)
    Generates an initial binary sample according to user specified models
    
    Parameters
    ----------
    cluster_profile : `str`
        Model to use for the cluster profile (i.e. sampling of the placement of objects in the cluster and their velocity within the cluster)
        Options include:
        'plummer' : Standard Plummer sphere.
            Additional parameters:
            'r_max' : `float`
                the maximum radius (in virial radii) to sample the clsuter
        'elson' : EFF 1987 profile.  Generalization of Plummer that better fits young massive clusters
            Additional parameters:
            'gamma' : `float`
                steepness paramter for Elson profile; note that gamma=4 is same is Plummer
            'r_max' : `float`
                the maximum radius (in virial radii) to sample the clsuter
        'king' : King profile 
            'w_0' : `float`
                King concentration parameter
            'r_max' : `float`
                the maximum radius (in virial radii) to sample the clsuter
    
    primary_model : `str`
        Model to sample primary mass; choices include: kroupa93, kroupa01, salpeter55, custom
        if 'custom' is selected, must also pass arguemts:
        alphas : `array`
            list of power law indicies
        mcuts : `array`
            breaks in the power laws.
        e.g. alphas=[-1.3,-2.3,-2.3],mcuts=[0.08,0.5,1.0,150.] reproduces standard Kroupa2001 IMF
    
    ecc_model : `str`
        Model to sample eccentricity; choices include: thermal, uniform, sana12
    
    porb_model : `str`
        Model to sample orbital period; choices include: log_uniform, sana12
    
    msort : `float`
        Stars with M>msort can have different pairing and sampling of companions
    
    pair : `float`
        Sets the pairing of stars M>msort only with stars with M>msort
    
    binfrac_model : `str or float`
        Model for binary fraction; choices include: vanHaaften or a fraction where 1.0 is 100% binaries
    
    binfrac_model_msort : `str or float`
        Same as binfrac_model for M>msort
    
    qmin : `float`
        kwarg which sets the minimum mass ratio for sampling the secondary
        where the mass ratio distribution is flat in q
        if q > 0, qmin sets the minimum mass ratio
        q = -1, this limits the minimum mass ratio to be set such that
        the pre-MS lifetime of the secondary is not longer than the full
        lifetime of the primary if it were to evolve as a single star
    
    m2_min : `float`
        kwarg which sets the minimum secondary mass for sampling
        the secondary as uniform in mass_2 between m2_min and mass_1
    
    qmin_msort : `float`
        Same as qmin for M>msort; only applies if qmin is supplied
    
    met : `float`
        Sets the metallicity of the binary population where solar metallicity is zsun 
    
    size : `int`
        Size of the population to sample
    
    zsun : `float`
        optional kwarg for setting effective radii, default is 0.02
    
    Optional Parameters
    -------------------
    virial_radius : `float`
        the initial virial radius of the cluster, in parsecs
        Default -- 1 pc
    
    tidal_radius : `float`
        the initial tidal radius of the cluster, in units of the virial_radius
        Default -- 1e6 rvir
    
    seed : `float`
        seed to the random number generator, for reproducability
    
    Returns
    -------
    Singles: `pandas.DataFrame`
        DataFrame of Single objects in the format of the InitialCMCTable
    
    Binaries: `pandas.DataFrame`
        DataFrame of Single objects in the format of the InitialCMCTable

Plummer Sphere

One of the most common initial conditions for a star cluster are those of Plummer (1911). They are not the most representative of star clusters in the local universe (see the King profile below), their popularity endures because their main quantities (such as enclosed mass, density, velocity dispersion, etc.) can be expressed analytically.

To generate a Plummer profile with 10,000 particles (here we only consider point-mass particles, as opposed to real stars), we can use the following COSMIC functions. Note that this returns two HDF5 tables, containing all the relevant parameters for the single stars and binaries. The two tables are linked, such that the Singles table contains all the single stars and binary centers of mass, as well as indices pointing to the relevant binaries in Binaries.

In [5]: from cosmic.sample import InitialCMCTable

In [6]: Singles, Binaries = InitialCMCTable.sampler('cmc_point_mass', cluster_profile='plummer', size=10000, r_max=100)

In [7]: print(Singles)
           id    k       m  Reff          r        vr        vt  binind
0         1.0  0.0  0.0001   0.0   0.029964  0.483857  0.856766     0.0
1         2.0  0.0  0.0001   0.0   0.033492 -0.461177  0.664207     0.0
2         3.0  0.0  0.0001   0.0   0.034694 -0.638679  0.997434     0.0
3         4.0  0.0  0.0001   0.0   0.047434 -0.468556  0.808870     0.0
4         5.0  0.0  0.0001   0.0   0.049135 -0.601994  0.748868     0.0
...       ...  ...     ...   ...        ...       ...       ...     ...
9995   9996.0  0.0  0.0001   0.0  29.373008  0.071705  0.206904     0.0
9996   9997.0  0.0  0.0001   0.0  37.808575  0.000110  0.067377     0.0
9997   9998.0  0.0  0.0001   0.0  45.144674  0.121053  0.032027     0.0
9998   9999.0  0.0  0.0001   0.0  52.101847 -0.011594  0.122011     0.0
9999  10000.0  0.0  0.0001   0.0  61.287149  0.051053  0.064878     0.0

[10000 rows x 8 columns]

Then we can save the initial conditions to an HDF5 file (for use in CMC) with:

In [8]: InitialCMCTable.write(Singles, Binaries, filename="plummer.hdf5")

InitialCMCTable.write automatically scales the cluster to Hénon units, where \(G = 1\), the cluster mass is \(M_{c}=1\), the cluster energy is \(E=-0.25\), and the virial radius is \(r_v \equiv G M_c^2 / 2|U| = 1\). These are the internal code units used by CMC.

Warning

InitialCMCTable.write scales the Singles and Binaries tables in place, so if you want to make any changes to the initial conditions after writing them, you’ll have to either create copies of the tables or generate new ones.

We can check that the Plummer function reproduces what we would expect from analytic predictions. The enclosed mass a Plummer sphere is given by

\[M(r) = M_{\rm total}\left(1 + \frac{a^2}{r^2}\right)^{-3/2}\]

where \(a\) is an arbitrary scale factor (which is \(3\pi / 16\) when the virial radius is normalized to 1). If we compare the mass-weighted cumulative radii of our Singles Pandas table to the analytic results, we can see:

In [9]: import numpy as np

In [10]: import matplotlib.pyplot as plt

In [11]: r_grid = np.logspace(-1.5,2,100)

In [12]: m_enc = (1 + 1/r_grid**2)**-1.5

In [13]: rv = 16/(3*np.pi) # virial radius for a Plummer sphere

In [14]: plt.plot(r_grid/rv,m_enc,lw=3);

In [15]: plt.hist(Singles.r,weights=Singles.m,cumulative=True,bins=r_grid);

In [16]: plt.xscale('log')

In [17]: plt.xlabel("Radii",fontsize=15);

In [18]: plt.ylabel(r"$M (< r) / M_{\rm total}$",fontsize=15);

In [19]: plt.legend(("Plummer","COSMIC Samples"),fontsize=14);
../_images/plot_plummer.png

Elson Profile

The Elson (1987) profile is a generalization of the Plummer profile that has been shown to better represent young massive clusters in the local universe. The density at a radius \(r\) is given by

\[\rho(r) = \rho_{0}\left(1 + \frac{a^2}{r^2}\right)^{-\frac{\gamma + 1}{2}}\]

Note that \(\gamma = 4\) gives a Plummer profile (the above code actually just calls the Elson profile generator with \(\gamma=4\)), though most young clusters are better fit with \(\gamma\sim2-3\). The enclosed mass is correspondingly more complicated:

\[M(r) = \frac{4 \pi \rho_0}{3} r^3 \,_2F_1\left(\frac{3}{2},\frac{\gamma + 1}{2} ; \frac{5}{2} ; -\frac{r^2}{a^2}\right)\]

Where \(\,_2F_1\) is the ordinary hypergeometric function.

Unlike both the Plummer and King profiles, the distribution function for the Elson profile cannot be written analytically. To generate the initial conditions, we directly integrate the density and potential functions to numerically compute \(f(E)\), and draw our velocity samples from that (see appendix B of Grudic et al., 2018). This produces a handful of warnings in the SciPy integrators, but the profiles that it generates are correct.

To generate an Elson profile with \(\gamma=2.5\), we can use

In [20]: Singles, Binaries = InitialCMCTable.sampler('cmc_point_mass', cluster_profile='elson', gamma=2.5, size=10000, r_max=100)

Comparing with the theoretical calculation for the enclosed mass, we find similarly good agreement:

In [21]: from scipy.special import hyp2f1

In [22]: gamma = 2.5

In [23]: def m_enc(gamma,r,rho_0):
   ....:     return 4*np.pi*rho_0/3 * r**3 * hyp2f1(1.5,(gamma+1)/2,2.5,-r**2)
   ....: 

In [24]: rv = 6. ## virial radius for a gamma=2.5 Elson profile

In [25]: rho_0 = 1/m_enc(gamma,100*rv,1)

In [26]: plt.plot(r_grid/rv,m_enc(gamma,r_grid,rho_0),lw=3); # note we scale by rv, rather than set scale factor

In [27]: plt.hist(Singles.r,weights=Singles.m,cumulative=True,bins=r_grid);

In [28]: plt.xscale('log')

In [29]: plt.xlabel("Radii",fontsize=15);

In [30]: plt.ylabel(r"$M (< r) / M_{\rm total}$",fontsize=15);

In [31]: plt.legend(("Elson ($\gamma=2.5$)","COSMIC Samples"),fontsize=14);
../_images/plot_elson.png

King Profile

An idealized cluster in thermodynamic equilibrium could be described as an isothermal sphere, where the velocities of stars resembled a Maxwell-Boltzmann distribution. But the isothermal sphere has infinite mass, and in any realistic star cluster, the distribution of stars should go to zero near the tidal boundary. The King (1966) profile accomplishes this by sampling from a lowered isothermal distribution

\[f(E) = \begin{cases} \rho_0 (2\pi\sigma^2)^{-3/2}(e^{E/\sigma^2})& E > 0;\\ 0& E \leq 0 \end{cases}\]

The King initial conditions can be created with COSMIC using:

In [32]: Singles, Binaries = InitialCMCTable.sampler('cmc_point_mass', cluster_profile='king', w_0=6, size=10000, r_max=100)

The analytic form of \(M(<r)\) cannot be written down for a King profile, but we can solve the ODE directly (this is done when generating the samples)

In [33]: from cosmic.sample.cmc import king

In [34]: radii,rho,phi,m_enc = king.integrate_king_profile(6)

In [35]: rho /= m_enc[-1]

In [36]: m_enc /= m_enc[-1]

In [37]: rv = king.virial_radius_numerical(radii, rho, m_enc) # just compute R_v numerically

In [38]: plt.plot(radii/rv,m_enc,lw=3);

In [39]: plt.hist(Singles.r,weights=Singles.m,cumulative=True,bins=r_grid);

In [40]: plt.xscale('log')

In [41]: plt.xlim(0.05,10);

In [42]: plt.xlabel("Radii",fontsize=15);

In [43]: plt.ylabel(r"$M (< r) / M_{\rm total}$",fontsize=15);

In [44]: plt.legend(("King ($w_0=6$)","COSMIC Samples"),fontsize=14);
../_images/plot_king.png

Realistic Initial Conditions

So far the above examples have only used the cmc_point_mass sampler. To generate realistic initial conditions, with stellar masses and binaries, we want to use the cmc sampler instead. This enables all the additional options found in the independent population sampler that COSMIC uses (see here for more details).

To generate the above King profile, but with all the additional stellar physics (initial mass function, binaries, etc.) we would use

In [45]: Singles, Binaries = InitialCMCTable.sampler('cmc', binfrac_model=0.1, primary_model='kroupa01',
   ....:                                             ecc_model='thermal', porb_model='log_uniform', qmin=0.1, m2_min=0.08,
   ....:                                             cluster_profile='king', met=0.00017, zsun=0.017, size=100000,w_0=6,
   ....:                                             seed=12345,virial_radius=1,tidal_radius=1e6)
   ....: 

This example is also found in the examples folder in the CMC repository.

Warning

At this step the metallicity (met) and value of the solar metallicity (zsun) must be specified to generate the correct initial radii for stars. If you have a very compact cluster and the value of zsun differs between here and your ini file, you may produce binaries which merge during cluster initialization, which breaks the code.