Stellar Mass-Luminosity Bias

Suppose that stars have a mass-luminosity relationship such that \(L \propto M^3\). If we have a flux-limited survey, it will bias us towards observing more massive stars which are not representative of the full mass distribution. Let’s see how to set this up in popsynth.

Setup the problem

First, we will assume that we have some initial mass function (IMF) for our stars that describes how there masses are distributed. For simplicity, we will assume that this IMF is just a log-normal distribution:

[1]:

%matplotlib inline import matplotlib.pyplot as plt from jupyterthemes import jtplot jtplot.style(context="notebook", fscale=1, grid=False) purple = "#B833FF" yellow = "#F6EF5B" import networkx as nx import numpy as np import warnings warnings.simplefilter("ignore") import popsynth # create a sampler for mass # we do not directly observe the mass as it is a latent quantity initial_mass_function = popsynth.LogNormalAuxSampler(name="mass", observed = False)

We now assume the dependent variable is the luminosity, so we need a DerivedLumAuxSampler that generates luminosities given a mass:

[2]:
class MassLuminosityRelation(popsynth.DerivedLumAuxSampler):
    _auxiliary_sampler_name = "MassLuminosityRelation"

    def __init__(self, mu=0.0, tau=1.0, sigma=1):
        # this time set observed=True
        super(MassLuminosityRelation, self).__init__("mass_lum_relation", uses_distance=False)

    def true_sampler(self, size):

            # the secondary quantity is mass

        mass = self._secondary_samplers["mass"].true_values

            # we will store the log of mass cubed
        self._true_values = 3 * np.log10(mass)

    def compute_luminosity(self):
        # compute the luminosity
            # from the relation
        return np.power(10., self._true_values)


luminosity = MassLuminosityRelation()


Now we can put everything together. First, we need to assign mass as a secondary quantity to the luminosity

[3]:
luminosity.set_secondary_sampler(initial_mass_function)

Finally, we will use a simple spherical geometry to hold our stars. We will also put a hard flux limit on our survey to simulate a flux-limited catalog.

[4]:
pop_gen = popsynth.populations.SphericalPopulation(1, r_max=10)

# create the flux selection

flux_selector = popsynth.HardFluxSelection()
flux_selector.boundary = 1e-2
pop_gen.set_flux_selection(flux_selector)

# now add the luminisity sampler

pop_gen.add_observed_quantity(luminosity)

Now let’s draw our survey.

[5]:

pop = pop_gen.draw_survey(flux_sigma=0.5)

We can now look at the distribution of the masses:

[6]:
fig, ax = plt.subplots()

bins = np.linspace(0,20,50)

ax.hist(pop.mass,
        bins=bins,
        label='all',
        color=purple,
        histtype="step",
        lw=2

       )

ax.hist(pop.mass[pop.selection],
        bins=bins,
        label='selected',
        color=yellow,
        histtype="step",
        lw=2
       )

ax.set_xlabel('stellar mass')
ax.legend()

[6]:
<matplotlib.legend.Legend at 0x7feb7ce83250>
../_images/notebooks_stellar_mass_11_1.png

We can see that indeed our selected masses are biased towards higher values.

Let’s look in the mass-luminostiy plane:

[7]:
fig, ax = plt.subplots()

bins = np.linspace(0,20,50)

ax.scatter(pop.mass[~pop.selection],
           pop.luminosities_latent[~pop.selection],
           label='not selected',
           color=purple,
           alpha=0.5,
           s=10



       )

ax.scatter(pop.mass[pop.selection],
           pop.luminosities_latent[pop.selection],
           label='selected',
           color=yellow,
           alpha=0.5,
           s=5
          )



ax.set_xlabel('stellar mass')
ax.set_ylabel('luminosity')
ax.set_yscale('log')


ax.legend()


[7]:
<matplotlib.legend.Legend at 0x7feb79aa1eb0>
../_images/notebooks_stellar_mass_13_1.png