Basic Properties and Star Formation History for NGC 337

The following example is intended to give an idea of the typical use-case of Lightning for a “normal” galaxy. We use a relatively simple, general model to fit the SED of NGC 337.

Data

We utilize the UV-to-submm photometry for NGC 337 from Table 2 in Dale et al. (2017), correcting it for Galactic extinction using the \(E(B-V)\) and \(A_\lambda / A_V\) given in their Tables 1 and 2.

The relevant data file can be found in examples/basic_NGC_337/ as ngc337_dale17_photometry.fits.

Configuration

Since this example is intended to be basic and general, we’ll stick close to the default configuration.

We’ll set the output filename:

OUTPUT_FILENAME: 'ngc337_lightning_output'     ,$

and increase the model uncertainty to \(5\%\):

MODEL_UNC: 0.05d0                             ,$

Since we have quite a lot of high-quality data, we’ll use the modified Calzetti Dust Attenuation curve for more flexibility:

ATTEN_CURVE: 'CALZETTI_MOD'                   ,$

but we’ll leave its priors and parameters at their default settings.

Everything else we want is already enabled by default: the restricted Draine and Li (2007) Dust Emission, the Affine-Invariant MCMC sampler, and the automatic Results and Post-Processing of the MCMC chains.

Running Lightning

Note

The IDL code snippets below are also available in batch file format, as examples/basic_NGC_337/NGC337_batch.pro.

To run Lightning, we now need only open an IDL session and do

restore, !lightning_dir + 'lightning.sav'

cd, !lightning_dir + 'examples/basic_NGC_337/'
lightning, 'ngc337_dale17_photometry.fits'

Analysis

We first load the results and check that the MCMC sampler converged and reached a statistically acceptable fit:

print, '//Convergence for NGC 337//'
print, 'Mean acceptance fraction: ', strtrim(mean(res[0].acceptance_frac), 2)
print, 'Convergence flag: ', strtrim(res[0].convergence_flag, 2)
print, 'Short chain flag: ', strtrim(res[0].short_chain_flag, 2)
print, 'Number of "stranded" walkers: ', strtrim(total(res[0].stranded_flag), 2)
print, 'PPC p-value: ', strtrim(res[0].pvalue, 2)
print, ''

//Convergence for NGC 337//
Mean acceptance fraction: 0.33788267
Convergence flag: 0
Short chain flag: 0
Number of "stranded" walkers: 1.00000
PPC p-value: 0.049500000

We can see that the acceptance fraction is good (more than \(20\%\) , less than \(50\%\) ) and that Lightning is confident that the MCMC sampler converged to a solution. We have the requested 1000 samples of our posterior distributions, and only 1 walker was “stranded.” The p-value is low but not extremely low, so we’ll inspect the best fit to see how it looks. The function NGC337_sed_plot.pro has been provided to help us do this:

p = NGC337_sed_plot(res)

The plot looks like so:

We can see from the SED plot and the residuals that the fit is actually quite good: it represents the data pretty well, and any under-fitting leading to the low p-value seems to be coming from one or two bands.

Now, one of the typical goals of SED fitting is constraint of the stellar mass and star formation rate. We’ll calculate these quantities (the stellar mass is actually already calculated in post-processing) and their posterior percentiles:

mstar_percentiles = percentile(res.mstar, [0.16, 0.50, 0.84])
SFR = 0.1 * res.psi[0,*] + 0.9 * res.psi[1,*]
SFR_percentiles = percentile(SFR, [0.16, 0.50, 0.84])

where we’re looking specifically at the star formation rate average over the last 100 Myr (i.e., the first two steps of our SFH). Then we can print them out

print, '//Basic properties for NGC 337//'
upstr = string(alog10(mstar_percentiles[2]) - alog10(mstar_percentiles[1]), format='(" (+",F5.2,")")')
downstr = string(alog10(mstar_percentiles[1]) - alog10(mstar_percentiles[0]), format='(" (-",F5.2,")")')
midstr = string(alog10(mstar_percentiles[1]), format='(F5.2)')
print, 'log(Mstar / Msun) = ', midstr, upstr, downstr

upstr = string(SFR_percentiles[2] - SFR_percentiles[1], format='(" (+",F5.2,")")')
downstr = string(SFR_percentiles[1] - SFR_percentiles[0], format='(" (-",F5.2,")")')
midstr = string(SFR_percentiles[1], format='(F5.2)')
print, 'SFR / [Msun yr-1] = ', midstr, upstr, downstr

and see

//Basic properties for NGC 337//
log(Mstar / Msun) =  9.55 (+ 0.07) (- 0.06)
SFR / [Msun yr-1] =  1.59 (+ 1.41) (- 1.00)

where the upper and lower uncertainties (from the 16th and 84th percentiles of the posterior) are given in parentheses. We can visualize the posteriors on these two quantities and the SFH with the NGC337_mstar_sfr_sfh.pro function:

p = NGC337_mstar_sfr_sfh_plots(res)

which produces the following graphic:

In this visualization, it’s easier to see that the SFR over the last 100 Myr is consistent with 0 – there is relatively little recent star formation. The SFH plot backs this up: while the overall shape of the SFH is uncertain (as we typically expect from SED fitting), it’s consistent with a declining burst of star formation between 100 Myr and 1 Gyr ago.

Of course, these are not the only properties of interest. Even the simple model in this example has a 10-dimensional parameter space, which gives us information mainly about the dust in the galaxy. For the final part of this example we’ll make the traditional corner plot of the posterior distribution to see the constraints on the model parameters.

parameter_arr = transpose([[reform(res.psi[0,*])], $
                           [reform(res.psi[1,*])], $
                           [reform(res.psi[2,*])], $
                           [reform(res.psi[3,*])], $
                           [reform(res.psi[4,*])], $
                           [res.tauv_diff], $
                           [res.delta], $
                           [res.umin], $
                           [res.gamma], $
                           [res.qpah]])

; Labels for each parameter
labels = ['$\psi_1$', '$\psi_2$', '$\psi_3$', '$\psi_4$', '$\psi_5$',$
          '$\tau_V$', '$\delta$',$
          '$U_{rm min}$', '$\gamma$', '$q_{\rm PAH}$']

; Use the corner_plot function from the lightning-visualization repository
; to make our 10D corner plot
p = window(dimension=[900, 900])
p = corner_plot(parameter_arr, labels, $
                /normalize, $
                /show_median, $
                contour_levels=[0.68, 0.95], $
                contour_smooth=3, $
                tickinterval=[1,1,1,0.75,0.25,0.1,0.5,1.0,0.01,0.005], $
                distribution_thick=2, $
                contour_thick=2, $
                font_size=8, $
                position=[0.1, 0.1, 0.9, 0.9], $
                /current)

This produces the following corner plot:

We can see that the galaxy is relatively unobscured (\(\tau_V < 0.3\) ), and that the dust attenuation and emission parameters are in line with expectations from local galaxies, with \(U_{\rm min} < 5\) and \(\gamma \approx 0.01\) , as in e.g. Dale et al. (2017).