Build a galaxy model from an interpolated SFH

GalaPy provides a non-parametric, interpolated, step-wise SFH model with derived components (like dust/gas mass and metallicity) treated as free parameters. This model is designed to predict the emission from galaxies for which the stellar mass growth history is available (e.g. obtained from hydro-dynamical simulations or with semi-analytical models) or to test the behaviour of exotic and arbitrarily complex SFH shapes.

Load dataset

The dataset at hand might look something like

import numpy
tau = numpy.array([
    0.0000000e+00, 2.0674900e+08, 4.2752600e+08, 6.6250900e+08,
    9.1218400e+08, 1.1766550e+09, 1.4557920e+09, 1.7499790e+09,
    2.0585520e+09, 2.3815860e+09, 2.7187860e+09, 3.0688570e+09,
    3.4325970e+09, 3.8068280e+09, 4.1934690e+09, 4.5884250e+09,
    4.9938270e+09, 5.4046790e+09, 5.8241890e+09, 6.2453270e+09,
    6.6728250e+09, 7.0994560e+09, 7.5283510e+09, 7.9562340e+09,
    8.3802740e+09, 8.8046430e+09, 9.2189360e+09, 9.6296370e+09,
    1.0036022e+10, 1.0428195e+10, 1.0814397e+10, 1.1194359e+10,
    1.1555720e+10
])
sfr = numpy.array([
    0.        , 0.        , 0.00258615, 0.04669582, 0.21555495,
    0.56162186, 0.60869537, 0.55047424, 0.36031803, 0.27614234,
    0.20609573, 0.16298589, 0.2459535 , 0.47935471, 0.57768621,
    0.55990713, 0.52003312, 0.48266521, 0.45552792, 0.43251383,
    0.40347906, 0.36103059, 0.28409527, 0.22353781, 0.17772765,
    0.15988216, 0.1646592 , 0.17408446, 0.19488885, 0.23297054,
    0.2524771 , 0.25091753, 0.23612502
])

where tau is an array of sampled ages (in units of \([\text{yr}]\)) in the evolution history of the galaxy and sfr is the star formation rate, \(\psi(\tau)\) (in units \([M_\odot/\text{yr}]\)) corresponding to each value in ŧau.

This particular data-set has been extracted from a catalogue of galaxies simulated with the semi-analytical model L-GALAXIES2020 in Parente et al. (2023) in which the autors include new features to the code, namely: (i) a state-of-the-art dust model that adopts the two-size approximation and (ii) a new disc instability criterion that triggers bulge and central black hole growth.

Let’s fix some further properties of the object simulated by the SAM:

# redshift of the source
zz = 1.e-4
# current age of the object
age = tau.max()
# average absolute metallicity
Zgxy = 6.74e-3
# dust mass
Mdust = 5.e+6

Build galaxy model

From the information above we can build a galaxy model with GalaPy, let’s first import the class wrapping up all the physical components necessary:

from galapy.Galaxy import GXY

We choose a parsec22.NTL SSP library that along with the emission from stellar atmospheres also includes (i) non-thermal emission due to synchrotron from super-novae and, (ii) nebular thermal emission including lines (see Choose the SSP library for further details).

We set the 'interpolated' model by passing the relevant keyword to the sfh argument dictionary and also set some further properties as derived from the SAM and listed above:

gxy = GXY(
    age = age, redshift = zz,
    csp = {'ssp_lib':'parsec22.NTL'}, # set the SSP library
    sfh = {
        'model':'interpolated',       # choose the interpolated model
        'tau':tau, 'sfr':sfr,         # pass the dataset simulated with SAM
        'Zgxy':Zgxy, 'Mdust':Mdust    # pass the further properties from the SAM
    },
)

The galaxy model above is sufficient to have a simulated emission from the object (and of course further physical properties can be tuned via the free-parameters of the model (e.g. properties of dust), for the full list of tunable free-parameters the user can refer to the following page in this docs: Free parameters of the galaxy model

Check the interpolation goes as expected

Let’s check what we are modelling by comparing the input SFH with the interpolated value returned by gxy.sfh. First, we can import the matplotlib.pyplot formatted version present in the galapy.analysis.plot submodule and some other function for cosmetics:

from galapy.analysis.plot import plt, format_axes_ticks

we also want to model the SFH on a time-grid thinner than the original one, but still containing the original values, we can build one as follows:

tt = numpy.sort(numpy.append(
    tau, # original array
    numpy.logspace( 8.0, numpy.log10(age), 256 ) # log-spaced grid of times
))

The tt array goes from \(10^8\ \text{yr}\) up to the age of the object.

# make a grid of subplots
fig, axs = plt.subplots(2,1, sharex=True,
                        gridspec_kw={'hspace':0.0},
                        tight_layout=True)

#########################################
# first sub-plot

# fix logarithmic scale for the axes and y-axis label
axs[0].set(
    #xscale='log',
    yscale='log',
    ylabel='$\\Psi(\\tau)\\ [M_\\odot\\;\\mathrm{yr}^{-1}]$'
)

# plot the interpolated model:
axs[0].plot(tt, gxy.sfh(tt),
            ls='--', color='grey') # cosmetics

# plot the SAM input values
axs[0].plot(tau, sfr,
            ls='none', marker='o', # from here just cosmetics
            markerfacecolor='white',
            markersize=6,
            markeredgewidth=1.75)

#########################################
# second sub-plot

# fix logarithmic scale for the axes and x- and y-axis labels
axs[1].set(
    #xscale='log',
    yscale='log',
    xlabel='$\\tau\\ [\\mathrm{yr}]$',
    ylabel='$M_\\star(\\tau)\\ [M_\\odot]$'
)

# also plot the evolution of the stellar mass
axs[1].plot(tt, gxy.sfh.Mstar(tt, 1000), color='tab:red')

[<matplotlib.lines.Line2D at 0x7fc79a0b9290>]

We can compute the stellar mass of the object at current age with

print( f'{gxy.sfh.Mstar(age):.3e} Msol' )

1.917e+09 Msol

Compute emission

We can use the gxy object to simulate a grid of fluxes received from the mock source

# array of rest-frame wavelengths
wave = gxy.wl()
# array of fluxes in mJy
flux = gxy.get_SED()

By calling the gxy.get_SED() method, also the contribution from the different components of the galaxy model is computed and stored internally, it can be accessed from the class attribute gxy.components which, though, returns luminosities, or by calling the dedicated method to convert it into fluxes:

components = gxy.components_to_flux()

Note that the components object is a dictionary where the keys point out the component short name and the values are arrays of fluxes in \(mJy\) computed on the same grid of wavelengths returned by the gxy.wl() method above

list(components.keys())

['stellar', 'extinct', 'MC', 'DD']

GalaPy provides some functions for formatted plotting of the SED, they are made available in the galapy.analysis.plot sub-module:

from galapy.analysis.plot import sed_layout, sed_components, sed_flux

We can therefore plot the overall SED and components by calling

fig, ax = plt.subplots(1,1,figsize=(6,3), constrained_layout=True)

# prepares the plot axes
ax = sed_layout(gxy.redshift, frame='rest', ax = ax, xlim=(50., 8.e+9), ylim=(1.e-6,1.e+10))

# plot the different components
_ = sed_components(
    wave, components,
    redshift=gxy.redshift, frame='rest',
    ax=ax
)

# plot the total flux:
_ = sed_flux(
    wave, flux,
    redshift=gxy.redshift, frame='rest',
    ax=ax
)

# plot a legend
ax.legend(ncols=2)

<matplotlib.legend.Legend at 0x7fc796d4a510>

Compute photometry

We might want to check what the emission from the SED model above looks like when transmitted through a set of bandpass filters.

First let’s select some filters from the database, let’s say we want the SDSS filters, Herschel’s PACS filters and ALMA bands we can extract the correct names and list them by calling

from galapy.PhotometricSystem import list_filters
bands = list_filters('SDSS')+list_filters('Herschel.PACS')+list_filters('ALMA')
print(bands)

['SDSS.r', 'SDSS.g', 'SDSS.i', 'SDSS.u', 'SDSS.z', 'Herschel.PACS.green', 'Herschel.PACS.blue', 'Herschel.PACS.red', 'ALMA.B7', 'ALMA.B3', 'ALMA.B6', 'ALMA.B8']

We can then build a photometric-system object by instantiating

from galapy.PhotometricSystem import PMS
pms = PMS(*bands)

The pivot wavelength associated with each band is contained in the lpivot attribute:

lpiv = pms.lpiv
lpiv

array([3.55652397e+03, 4.70249528e+03, 6.17557888e+03, 7.48997685e+03,
       8.94670956e+03, 7.14985538e+05, 1.01855375e+06, 1.64412822e+06,
       6.96722749e+06, 9.46015413e+06, 1.24698305e+07, 3.08925524e+07])

which can be plotted using

from galapy.analysis.plot import photometric_system as pms_plot
fig, ax = plt.subplots(1,1,figsize=(12,3), constrained_layout=True)
_ = pms_plot(pms, ax=ax)

To see the emission in bands we will have to build a galapy.Galaxy.PhotoGXY model instead of the galapy.Galaxy.GXY model we have built above.

from galapy.Galaxy import PhotoGXY

Note that PhotoGXY inherits from GXY, therefore all the functionalities of the gxy object built in previous sections of this tutorial.

Equivalently to what done in the previous sections, we can build a photo-galaxy object by calling

pgxy = PhotoGXY(
    pms = pms, # <--- this is the only argument different from the base GXY object
    age = age, redshift = zz,
    csp = {'ssp_lib':'parsec22.NTL.refined'}, # set the SSP library (in its 'refined' version)
    sfh = {
        'model':'interpolated',       # choose the interpolated model
        'tau':tau, 'sfr':sfr,         # pass the dataset simulated with SAM
        'Zgxy':Zgxy, 'Mdust':Mdust    # pass the further properties from the SAM
    },
)

Note that, since the ALMA bands act on a region of the wavelength space that is otherwise under-sampled, we have also changed the SSP library, using in this case the refined version of the same model

We can compute the transmitted flux in the different bands with

pflux = pgxy.photoSED()

and we can plot it together with the full SED flux, being carefull to use observed wavelength this time, and not the rest-frame one (this can be done by calling the function returning the wavelength grid with the argument pgxy.wl(obs=True) and setting the frame of the plotting functions from 'rest' to either 'obs' or 'both' (as done in the plot below)

fig, ax = plt.subplots(1,1,figsize=(6,3), constrained_layout=True)

# prepares the plot axes
ax = sed_layout(pgxy.redshift, frame='both', ax = ax, xlim=(50., 8.e+9), ylim=(1.e-6,1.e+10))

# plot the total flux:
_ = sed_flux(
    pgxy.wl(obs=True), pgxy.get_SED(),
    redshift=pgxy.redshift, frame='both',
    ax=ax
)

# plot the band fluxes here
plt.scatter(
    lpiv, pflux,
    color=plt.cm.plasma(numpy.linspace(0.1,0.9,lpiv.size)),
    zorder=2 # otherwise pyplot would put it automatically behind the solid line
)

# plot a legend
ax.legend(ncols=2)

<matplotlib.legend.Legend at 0x7fc796780ad0>