Stellar & Nebular Emission - PEGASE+Cloudy

In these models, the nebular component is generated from Cloudy simulations on a large grid of metallicities and ionization parameters (see Nebular Emission Recipe).

The source models are from PEGASE, assuming a Kroupa et al. (2001) IMF. Since PEGASE models older than 30 Myr produce no Lyman continuum emission by definition, a nebular component is not included at ages older than 30 Myr.

Imports

import numpy as np
from scipy.interpolate import interp1d
from lightning.stellar import PEGASEModelA24 as StellarModel
from lightning.stellar import PEGASEBurstA24 as BurstModel
from lightning.sfh import PiecewiseConstSFH
import matplotlib.pyplot as plt
import matplotlib as mpl
plt.style.use('lightning.plots.style.lightning-serif')
%matplotlib inline

Initialize Model

Here we’ll initialize our sort of ‘default’ stellar population model, which is integrated over stellar age bins to be used with a piecewise-constant SFH. We’ll do it twice, with and without the nebular component, to compare.

wave_grid = np.logspace(np.log10(0.01),
                        np.log10(10),
                        200)
filter_labels = ['SDSS_u', 'SDSS_g', 'SDSS_r', 'SDSS_i', 'SDSS_z',
                 'MOIRCS_J', 'MOIRCS_H', 'MOIRCS_Ks',
                 'IRAC_CH1', 'IRAC_CH2', 'IRAC_CH3', 'IRAC_CH4']
redshift = 0.0
age = [0, 1e7, 1e8, 1e9, 5e9, 13.4e9]

stars_neb = StellarModel(filter_labels,
                         age=age,
                         redshift=redshift,
                         wave_grid=wave_grid,
                         nebular_effects=True)

stars_noneb = StellarModel(filter_labels,
                           age=age,
                           redshift=redshift,
                           wave_grid=wave_grid,
                           nebular_effects=False)

burst = BurstModel(filter_labels, redshift=redshift, wave_grid=wave_grid)

Simple Stellar Population Models

We’ve included the left-hand panel of this plot in basically every Lightning-related paper for years.

fig, axs = plt.subplots(1,2, figsize=(8,4))

Nmod = len(age) - 1

cm_neb = mpl.colormaps['Oranges']
colors_neb = cm_neb(np.linspace(0.2, 0.9, Nmod))

cm_noneb = mpl.colormaps['Greens']
colors_noneb = cm_noneb(np.linspace(0.2, 0.9, Nmod))

for i in np.arange(Nmod):

    finterp = interp1d(stars_neb.wave_grid_rest, stars_neb.Lnu_obs[i,-2,1,:])
    f1 = finterp(1)

    axs[0].plot(stars_neb.wave_grid_rest,
                stars_neb.nu_grid_obs * stars_neb.Lnu_obs[i,-2,1,:] / f1,
                color=colors_neb[i])

    finterp = interp1d(stars_noneb.wave_grid_rest, stars_noneb.Lnu_obs[i,-2,:])
    f1 = finterp(1)

    axs[1].plot(stars_noneb.wave_grid_rest,
                stars_noneb.nu_grid_obs * stars_noneb.Lnu_obs[i,-2,:] / f1,
                color=colors_noneb[i])

axs[0].set_xscale('log')
axs[0].set_yscale('log')
axs[0].set_ylim(1e9, 5e17)
axs[0].axvline(0.0912, color='dimgray', linestyle='--')

axs[0].set_xlabel(r'Rest-Frame Wavelength [$\rm \mu m$]')
axs[0].set_ylabel(r'$\nu L_\nu\ [\rm a.u.]$')
axs[0].set_title('With Nebula')

axs[1].set_xscale('log')
axs[1].set_yscale('log')
axs[1].set_ylim(1e9, 5e17)
axs[1].axvline(0.0912, color='dimgray', linestyle='--', label='Lyman limit')
axs[1].legend(loc='lower right')

axs[1].set_xlabel(r'Rest-Frame Wavelength [$\rm \mu m$]')
axs[1].set_title('Without Nebula')

Text(0.5, 1.0, 'Without Nebula')

Where darker shades represent older ages. Note that nebular emission lines are only present for the two models with ages 100 Myr. The effect of free-free nebular continuum emission can be seen by comparing the lightest yellow curve on the left with the lightest green curve on the right. Note also that the Lyman continuum emission is not completely attenuated by the nebula (and old stellar populations with hot stripped stars produce low-level Lyman continuum emission). However, our dust models are opaque to the Lyman continuum by definition.

Note that the PEGASEBurstA24 model class provides a cleaner interface to the single-age population models, and can actually be used to fit SEDs. We’ll use it now to look at the variation of the spectrum with \(Z\) and \(\log U\) for a 1 Myr old population:

params_Z = np.tile([6.0, 6.0, 0, -2.0], (len(burst.Zmet), 1))
params_Z[:,2] = burst.Zmet

params_logU = np.tile([6.0, 6.0, 0.02, 0.0], (len(burst.logU), 1))
params_logU[:,3] = burst.logU

lnu_burst_Z,_,_ = burst.get_model_lnu_hires(params_Z)
lnu_burst_logU,_,_ = burst.get_model_lnu_hires(params_logU)

fig, axs = plt.subplots(1,2, figsize=(8,4))

Nmod = len(burst.Zmet)
cm_neb = mpl.colormaps['Oranges']
colors_Z = cm_neb(np.linspace(0.2, 0.9, Nmod))

for i in range(Nmod):

    axs[0].plot(burst.wave_grid_obs,
                burst.nu_grid_obs * lnu_burst_Z[i,:],
                color=colors_Z[i])

axs[0].set_xscale('log')
axs[0].set_yscale('log')
axs[0].set_xlim(0.01, 10)
axs[0].set_ylim(1e4, 5e10)
axs[0].set_title(r'Varying $Z$')
axs[0].set_ylabel(r'$\nu L_\nu\ [\rm L_{\odot}]$')
axs[0].set_xlabel(r'Rest-Frame Wavelength [$\rm \mu m$]')

Nmod = len(burst.logU)
cm_noneb = mpl.colormaps['Greens']
colors_logU = cm_noneb(np.linspace(0.2, 0.9, Nmod))

for i in range(Nmod):

    axs[1].plot(burst.wave_grid_obs,
                burst.nu_grid_obs * lnu_burst_logU[i,:],
                color=colors_logU[i])

axs[1].set_xscale('log')
axs[1].set_yscale('log')
axs[1].set_xlim(0.01, 10)
axs[1].set_ylim(1e4, 5e10)
axs[1].set_title(r'Varying $\log U$')
axs[1].set_xlabel(r'Rest-Frame Wavelength [$\rm \mu m$]')

[[ 6.00000000e+00  1.00000000e+06  6.47187320e-04 -2.00000000e+00]
 [ 6.00000000e+00  1.00000000e+06  8.14760564e-04 -2.00000000e+00]
 [ 6.00000000e+00  1.00000000e+06  1.02572278e-03 -2.00000000e+00]
 [ 6.00000000e+00  1.00000000e+06  1.29130847e-03 -2.00000000e+00]
 [ 6.00000000e+00  1.00000000e+06  1.62566105e-03 -2.00000000e+00]
 [ 6.00000000e+00  1.00000000e+06  2.04658600e-03 -2.00000000e+00]
 [ 6.00000000e+00  1.00000000e+06  2.57649913e-03 -2.00000000e+00]
 [ 6.00000000e+00  1.00000000e+06  3.24362023e-03 -2.00000000e+00]
 [ 6.00000000e+00  1.00000000e+06  4.08347593e-03 -2.00000000e+00]
 [ 6.00000000e+00  1.00000000e+06  5.14079162e-03 -2.00000000e+00]
 [ 6.00000000e+00  1.00000000e+06  6.47187320e-03 -2.00000000e+00]
 [ 6.00000000e+00  1.00000000e+06  8.14760564e-03 -2.00000000e+00]
 [ 6.00000000e+00  1.00000000e+06  1.02572278e-02 -2.00000000e+00]
 [ 6.00000000e+00  1.00000000e+06  1.29130847e-02 -2.00000000e+00]
 [ 6.00000000e+00  1.00000000e+06  1.62566105e-02 -2.00000000e+00]
 [ 6.00000000e+00  1.00000000e+06  2.04658600e-02 -2.00000000e+00]
 [ 6.00000000e+00  1.00000000e+06  2.57649913e-02 -2.00000000e+00]]
[[ 6.0e+00  1.0e+06  2.0e-02 -1.5e+00]
 [ 6.0e+00  1.0e+06  2.0e-02 -2.0e+00]
 [ 6.0e+00  1.0e+06  2.0e-02 -2.5e+00]
 [ 6.0e+00  1.0e+06  2.0e-02 -3.0e+00]
 [ 6.0e+00  1.0e+06  2.0e-02 -3.5e+00]
 [ 6.0e+00  1.0e+06  2.0e-02 -4.0e+00]]

Text(0.5, 0, 'Rest-Frame Wavelength [$\\rm \\mu m$]')

The models with lower metallicity (lighter colors) have more significant ionizing fluxes, while the models with larger ionization parameters have more luminous strong line emission and slightly larger ionizing escape from the nebula.

Composite Stellar Population Models

To construct composite stellar populations we must of course assume a SFH for the population. Since we binned our simple stellar populations, we must use the PiecewiseConstantSFH model.

sfh = PiecewiseConstSFH(age)

fig, axs = plt.subplots(1,2, figsize=(8,4))

coeffs= np.array([[1,1,0,0,0],
                  [0,0,1,1,1]])

Z = np.array([0.02, 0.02])
logU = np.array([-2.0, -2.0])

params = np.stack([Z, logU], axis=-1)

Nmod = coeffs.shape[0]

cm_neb = mpl.colormaps['Oranges']
colors_neb = cm_neb(np.linspace(0.2, 0.9, Nmod))

cm_noneb = mpl.colormaps['Greens']
colors_noneb = cm_noneb(np.linspace(0.2, 0.9, Nmod))

lnu_hires_neb,_,_ = stars_neb.get_model_lnu_hires(sfh, coeffs, params)
lnu_hires_noneb,_,_ = stars_noneb.get_model_lnu_hires(sfh, coeffs, Z)

for i in np.arange(Nmod):

    axs[0].plot(stars_neb.wave_grid_rest,
                stars_neb.nu_grid_obs * lnu_hires_neb[i,:],
                color=colors_neb[i])

    axs[1].plot(stars_noneb.wave_grid_rest,
                stars_noneb.nu_grid_obs * lnu_hires_noneb[i,:],
                color=colors_noneb[i])

axs[0].set_xscale('log')
axs[0].set_yscale('log')
axs[0].set_ylim(2e6, 1e11)

axs[0].set_xlabel(r'Rest-Frame Wavelength [$\rm \mu m$]')
axs[0].set_ylabel(r'$\nu L_\nu\ [\rm L_{\odot}]$')
axs[0].set_title('With Nebula')

axs[1].set_xscale('log')
axs[1].set_yscale('log')
axs[1].set_ylim(2e6, 1e11)
axs[1].axvline(0.0912, color='dimgray', linestyle='--', label='Lyman limit', zorder=-1)
axs[1].legend(loc='lower right')

axs[1].set_xlabel(r'Rest-Frame Wavelength [$\rm \mu m$]')
axs[1].set_title('Without Nebula')

Text(0.5, 1.0, 'Without Nebula')

Here the lighter colored populations are star forming, and the darker ones are quiescent. Note that, as implied above, even if you were to use stellar populations without nebular extinction for your SED fitting in Lightning, the resulting galaxy would have no Lyman continuum leakage, as our ISM attenuation models are defined to be opaque to Lyman continuum radiation.