Many astronomical investigations require simulating populations of stars, and isochrones contains some utilities to help enable this. Given population distributions of the quantities required to simulate individual stars, a StarPopulation object can be defined and used to generate sample populations following this distribution. Binary stars, ubiquitous as they are, are necessarily built into this framework, so the parameters needed to simulate an individual stellar observation are the following:
where $M_A, M_B$ are the primary and (if present) secondary masses, $T$ is age, $[Fe/H]$ is the metalicity, $d$ is distance, and $A_V$ is the $V$-band extinction, quantifying the effect of dust along the line of sight. Generating a population of such stars then requires sampling from distributions of each of the above quantities. A StarPopulation takes metallicity, distance, and extinction distributions as arguments, and samples from each of those distributions when generating a sample population.
Sampling primary/secondary masses and ages is a bit less straightforward. For $M_A, M_B$, isochrones parametrizes the distribution with a primary initial mass function (IMF), binary fraction $f_B$, and mass-ratio ($q = M_B/M_A$) distribution $p(q) \propto q^\gamma$. The age distribution of stars in a population is often described as a "star-formation history" (SFH)---sampling a population with a given SFH is the same as treating the SFH as the probability distribution function of stellar age, sampling ages from this distribution, and then truncating any stars that have reached the end of their evolution. Practically, this truncation happens by rejection sampling: evaluating the ModelGridInterpolator at each sampled set of parameters, and rejecting samples for which the interpolator returns np.nan values for the observed stellar properties (which will happen when trying to interpolate out-of-bounds, which happens when a star is requested beyond the end of its lifetime).
In [1]:
from scipy.stats import uniform, norm
from isochrones import get_ichrone
from isochrones.priors import GaussianPrior, SalpeterPrior, DistancePrior, FlatPrior
from isochrones.populations import StarFormationHistory, StarPopulation
# Initialize interpolator
mist = get_ichrone('mist')
# Initialize distributions
# Ingredients required to generate primary & secondary masses
imf = SalpeterPrior(bounds=(1, 10)) # minimum 1 Msun
fB = 0.4
gamma = 0.3
# SFH distribution takes a scipy stats distribution, of age in Gyr
sfh = StarFormationHistory(dist=uniform(0, 10))
# The following are all isochrones.priors.Prior objects,
# or anything with a .sample(N) method
feh = GaussianPrior(-0.2, 0.2)
distance = DistancePrior(max_distance=3000)
AV = FlatPrior(bounds=[0, 1])
pop = StarPopulation(mist, imf=imf, fB=fB, gamma=gamma, sfh=sfh, feh=feh, distance=distance, AV=AV)
Once the object is created, it can be used to generate a population of stars.
In [2]:
df = pop.generate(1000)
df.head()
Out[2]:
Note that this operation is not nearly as fast as directly interpolating an isochrone or evolution track grid (since generating properites given mass, age, and metallicity necessarily involves solving for EEP first):
In [3]:
%timeit pop.generate(1000)
Also, this can be made much faster if you loosen the requirement on getting exactly a particularly desired number of stars (as part of the generating algorithm involves replacing stars that come out as nan until no nans are left):
In [4]:
print(len(pop.generate(1000, exact_N=False)))
%timeit pop.generate(1000, exact_N=False)
The full column list of this table of simulated stars is the following:
In [5]:
', '.join(df.columns)
Out[5]:
All quantities with a tag _0 refer to the primary star; all quantities with _1 refer to the secondary. Columns ending in just _mag represent the combined magnitude of both primary and secondary component. Let's look the Gaia color-magnitude diagram for this simulated population. Note also the A_[x] columns, which give the specific extinction per band for each system (and for the individual components of the binary).
In [6]:
import holoviews as hv
hv.extension('bokeh')
import hvplot.pandas
def hr_plot(df):
df['BpRp'] = df.BP_mag - df.RP_mag
hr = df.hvplot.scatter('BpRp', 'G_mag',
hover_cols=['mass_0', 'mass_1', 'age_0', 'AV_0'],
color='feh_0')
return hr.options(height=400, width=500, invert_yaxis=True)
hr_plot(df)
Out[6]:
There is also a simple utility function that can "deredden" a generated population dataframe (e.g., recover the true intrinsic magnitudes of each star in the absence of dust), by subtracting off the A_x extinction values from the magnitudes, and setting all extinctions to zero. Let's use this to deredden the above hr diagram:
In [7]:
from isochrones.populations import deredden
dereddened = deredden(df)
hr_plot(df).options(size=3, alpha=0.2, color='red') * hr_plot(dereddened).options(alpha=0.2, color='black', size=3)
Out[7]:
See how the dust (reddened points) moves each star down (fainter) and to the right (redder).
In [8]:
mass_A = 1.0
mass_B = [0.8, 0.6, 0.4, 0.2]
age, feh, distance, AV = (9.6, 0.02, 100, 0.1)
mist.generate_binary(mass_A, mass_B, age, feh, distance=distance, AV=AV)[['G_mag', 'BP_mag', 'RP_mag']]
Out[8]: