Extragalactic

I’m going to end this series for now with a more detailed look at the spatially resolved star formation history of this galaxy and make a highly selective comparison to some recent theoretical and empirical work on mergers. This continues from part 1 and part 2.

My main theoretical sources for the following discussion include a paper titled “The fate of the Antennae galaxies” by Lahén et al. (2018) and a similar paper also simulating an Antennae analog by Renaud, Bournaud and Duc (2015). I didn’t pick these because I necessarily think this galaxy is an evolved analog of the Antennae; in fact I think there are important differences in the likely precursors. The main reason is these studies are among the first to run very high resolution simulations through and beyond coalescence. Some other significant papers include high resolution simulations with detailed treatment of feedback by Hopkins et al. (2013),  Di Matteo et al. (2008), who studied a large sample of merger and flyby simulations, Bekki et al. (2005) who specifically looked at the formation of K+A galaxies through mergers, and Ji et al. (2014) who studied the lifetime of merger features using a suite of simulations. This doesn’t begin to scratch the surface of this literature. Merger simulations are very popular!

Some recent observational papers that examine individual post-merger systems in more or less detail include Weaver et al. (2018) and Li et al. (2018). I will return to the latter paper, which was based on MaNGA data, in a later post. Other observational work that’s more statistical in nature includes Ellison et al. (2015), Ellison et al. (2013), and earlier papers in the same series.

As I’ve mentioned several times already the model spectra I use mostly come from the 2017 EMILES extension of the MILES SSP library with BaSTI isochrones and Kroupa IMF, supplemented with unevolved model spectra from  the 2013 update of BC03 models. The results reported in the previous two posts used a small subset of the library — just 3 metallicity bins and 27 time bins, with every other time bin in the full library starting with the second excluded. For my “full” MILES subset I use 4 metallicity bins with \([Z/Z_{\odot}] \in \{-0.66, -0.25, 0.06, 0.40\}\) and all 54 time bins (the lowest metallicity bin is dropped in the reduced subset). Execution time of the Stan SFH models is roughly proportional to the number of parameters, so everything else equal it takes about 2.5 times longer to run a single model with the full set.

For this exercise I binned the MaNGA spectra with a considerably higher target SNR than usual, ending up with 53 coadded spectra out of the original 183 fiber/pointing combinations. One of the mods I made to Cappellari’s Voronoi binning code was to drop a “roundness” check. The bins in this case still end up with reasonably compact shapes since the surface brightness within the IFU footprint is fairly symmetrical. All but one of the spectra in the inner 2 kpc. ended up unbinned, so the spatial resolution in that critical region is unchanged.

To directly compare the results from the current round  of models with the previous lower time resolution runs here are the modeled mass growth histories summed over the entire IFU footprint (blue is the EMILES subset):

Note: MaNGA’s dithering strategy results in considerable overlap in fiber positions in order to obtain a 100% filling fraction. That means the area in fibers exceeds the area of the IFU footprint by about 60%. So, sums of quantities like masses or star formation rates need to be adjusted downwards by about 0.2 dex.

The only real difference we see in the higher time resolution runs is that the initial, strongest starburst is a few percent weaker, and there are some subtle differences in late time behavior. The first large starburst begins at 1.25 Gyr in both sets of runs, and unfortunately this bin is 250 Myr in width in the full resolution set, vs 350 Myr in the subset, so we don’t really get significantly finer grained estimates of the SFH in the initial burst.

We can get a rough idea of the nature of the precursors from this graph. The present day stellar mass is \(4 \times 10^{10} M_{\odot}\), of which just about a third was formed in or after the initial starburst. The present day combined stellar mass of the precursors was just over \(2.6 \times 10^{10} M_{\odot}\) just before the burst — the then stellar mass would have been somewhat larger since this mass growth estimate includes mass loss to stellar evolution. Also, there’s considerable light and therefore stellar mass outside the IFU footprint — the drpcat value for the effective radius is 5.5″ while the Petrosian radius per the SDSS photo pipeline is 12.5″. The IFU footprint is about 10″ radius, or about 1.8 effective radii. Taking a guess that we’re sampling about 80% of the stellar mass the values above need to be adjusted upwards by 25%. Assuming the precursors were equal mass and rounding up their pre-merger stellar masses would be around \(1.65  \times 10^{10} M_{\odot}\). Guessing the current total stellar mass to be \(5 \times 10^{10} M_{\odot}\) then \(1.7 \times 10^{10} M_{\odot}\) was added by the merger. The amount of gas turned into stars was actually higher — about \(2.4 \times 10^{10} M_{\odot}\), but some, maybe most, of the gas lost to stellar evolution might have been recycled. But, I’ll be conservative and adopt the higher value. If that was divided equally between the progenitors they would have had gas masses around \(1.2 \times 10^{10} M_{\odot}\) just before the merger, or just over 40% of the baryonic mass. While high that’s not extraordinarily so for a late type spiral in that stellar mass range (see for example Ellison et al. 2015, figures 6-7).

These values are rather different from the putative Antennae analogs in the two studies linked above, which are in turn rather different from each other. Lahén et al. assign equal stellar masses of \(4.7 \times 10^{10} M_{\odot}\) to each progenitor, equal gas masses of \(0.8 \times 10^{10} M_{\odot}\) (14.5% of the total baryonic mass), and a baryonic to total mass ratio of 0.1. The simulations of Reynaud et al. have stellar masses of \(6.5 \times 10^{10} M_{\odot}\) apiece, gas masses of just \(0.65 ~\mathrm{and}~ 0.5 \times 10^{10} M_{\odot}\), and dark matter halos of \(23.8 \times 10^{10} M_{\odot}\) (for a baryonic to total mass ratio of about 23%). The merger timescales are rather different too: 180 Myr from first pericenter passage to coalescence for the Reynaud et al. simulations vs. ~600 Myr for Lahén et al. (the latter seems closer to the observational evidence; for example Whitmore et al. 1999 found clusters with 3 distinct age ranges up to ~500Myr, with the oldest argued to have formed in the first pericenter passage).

The qualitative features of the merger progress are similar however, and some are fairly generically observed in simulations. There are two pericenter passages, with a second period of separation followed shortly by coalescence at the third pericenter. Star formation is widespread but clumpy in the first flyby, with a large but short starburst and slow decline. A stronger, shorter, and more centrally concentrated starburst occurs around the time of the second pericenter passage through coalescence. The Lahén simulations follow the merger remnant for another Gyr — they show an exponentially declining SFR from an elevated level immediately after coalescence. Neither of these sets of simulations attempt to model AGN feedback, which would presumably cause more rapid quenching. Notably, the Lahén merger remnant develops a kinematically decoupled core although overall the remnant is a fast rotator. Unless the rotation axis is pointing right at us this galaxy is a slow rotator except for the core.

Getting back to data for this galaxy I show some star formation history plots below, first for the inner 9 fibers (all within <1.25 kpc of the nucleus). Below that are a pair of plots of aggregate SFH for fibers within 2.5 kpc (left) of the nucleus and outside that radius (right). The time axis on these plots is logarithmic since I want to focus on the late time behavior (and also this is closer to the real time resolution). Note that in the inner region there are 3 broad peaks in star formation rate — the previously discussed one at 1-1.25Gyr, one centered around 500Myr, and a recent (\(\lesssim\)100Myr revival that peaks at the youngest BaSTI age of 30Myr. The relative and absolute sizes of the peaks are highly variable, indicating that the stars aren’t fully relaxed yet. The outer region covered by the IFU by contrast has just a single peak at 1-1.25Gyr, with a more or less monotonic decline thereafter.

If we believe the models the straightforward interpretation is the first pericenter passage happened about 1-1.25Gyr ago, with coalescence around 500Myr ago and continued infall or recycling driving the ongoing star formation. This is consistent with simulations, which can have merger timescales ranging anywhere from a few hundred Myr to several Gyr. The main problem with this interpretation is that in most simulations the highest peak SFR occurs around coalescence, with the integrated SF roughly equally divided between the first passage and during/after coalescence. In these models about 75% of the post-burst star formation occurs in the 1-1.25Gyr bin, with about 23% in the later burst. An alternative scenario then would be that the entire merger sequence occurred in the 250Myr window of the 1-1.25Gyr bin, and everything since is due to post-merger infall. Of course a third possibility, which I don’t dismiss, is that the details of the modeled episodes of star formation are just artifacts having no particular relationship to the real recent star formation history. The main argument I have against that is that all the spectra seem to tell a consistent story of multiple recent episodes of star formation of varying strength in space and time, and with clear radial trends.

Regardless of the details this galaxy nicely confirms several of the predictions of Bekki et al. (linked above). As we saw in the last post there is a positive color gradient and negative Balmer absorption gradient as they predict for dissipative major mergers, with a break in the color gradient trend that corresponds approximately with the boundary between multi-peaked and single peaked late time star formation histories.

The obvious morphological indicators of disturbance are provocative but don’t seem to be highly constraining. My recollection is that folklore guesstimates that tidal features have a visible lifetime of ~1 Gyr. Ji et al. (linked above) on the other hand found based on a suite of simulations that merger features remain visible at the 25 mag./arcsec\(^2\) somewhat longer, on average about 1.4Gyr. It may be significant that the disturbance is easily seen even in the SDSS Navigate imaging, especially the northern loop. The somewhat similar star forming elliptical SDSS J1420555.01+400715.7 discussed by Li et al. (which I will turn to later) required deeper imaging and some enhancement to show signs of disturbance. One thing I’m working on learning is galaxy profile fitting and photometric decomposition. I’d like to see if the single component Sersic fit that fit the mass density profile also works for the surface brightness (I suspect there is a steeper core). I’d also like to quantify the surface brightness in the tidal features. These are exercises for later.

I’m going to conclude with a few plots of quantities that I thought might benefit from the higher S/N in the binned spectra. First is the gas phase metallicity from the O3N2 index plotted against radius. Unfortunately emission line strengths are still too weak beyond about 2.5 kpc radius to say much about the outskirts, but we still see a fairly flat trend with radius up to ~1.5 kpc with a turnover to lower values that continues at least out to ~2.5kpc. This is broadly similar with the simulations of Lahén et al., who obtained a shallower metallicity gradient for their merger remnant than an undisturbed spiral galaxy.

Here I again plot the estimated SFR density against radius and the estimated SFR density against Hα luminosity density. The high luminosity end is essentially unchanged from the plot in the last post since these spectra are unbinned. Again, the points at the high luminosity end lie above the Hα-SFR calibration of Moustakas et al. by considerably more than the nominal length of the error bars, which points to a likely recent (< ~10Myr) fading of star formation. This is consistent with the detailed SFH models. At the other end of the luminosity scale there seems to be an excess of points to the right of the calibration line. This could indicate that the main ionization source in the outer regions of the IFU footprint is something other than massive young stars (presumably shocks or hot evolved stars).

Despite the “composite” emission line ratios in the central region there’s little evidence for an AGN. The galaxy is a compact FIRST radio source, but the size of the starforming region is right around the FIRST resolution. The raw WISE colors are W1-W2 = +0.07, W2-W3 = +2.956, which is well within the locus of spiral colors.

In part two of this series that I began here I’m going to look at various quantities derived from the modeled star formation histories and emission line fits. Later I’ll dive into some selected recent literature on major gas rich mergers and do some comparisons. I’m especially interested in whether the model star formation histories have anything to tell us about the chronology of the merger and the nature of the progenitors.

First, here is a representative sample of spectra. These are reduced to the local rest frame and fluxes are corrected for foreground galactic extinction but not any internal attenuation. The middle spectrum in the second row is from the central fiber/pointing combination (remember these are from the stacked RSS data). The peripheral spectra are roughly equally spaced along a ring about 7″ or 4kpc from the center.

The central spectrum does show moderately strong emission lines along with fairly strong Balmer absorption, which points to some ongoing star formation along with a significant intermediate age population. Emission lines are evidently lacking or very weak in the peripheral spectra (most of the spikes are terrestrial night sky lines or simply noise), indicating the galaxy is currently passively evolving at 4kpc radius. In more detail…

The most basic but also finest grain quantities I look at are posterior estimates of star formation histories and mass growth histories. These contain similar, but not quite the same information content. I show star formation histories as the “instantaneous” star formation rate vs. look back time defined as the total stellar mass born in an age bin divided by its width in years.

Mass growth histories estimate the cumulative fraction of the present day stellar mass. These incorporate a recipe for mass loss to stellar evolution and as such aren’t quite integrals of the star formation histories. By convention remnant masses are included in stellar mass, and there are recipes for this as well (in this case taken from the BaSTI web page). I display posterior marginal means and 95% confidence intervals for both of these against time.

I think I first encountered empirical estimates of mass growth histories in McDermid et al. (2015), although earlier examples exist in the literature (eg Panter et al. 2007). I may be the first to display them for an individual galaxy at the full time resolution of the input SSP library. Most researchers are reluctant to make detailed claims about star formation histories of individual galaxies based on SED modeling. I’m not quite so inhibited (not that I necessarily believe the models).

The most obvious feature in these is a burst of star formation that began ≈1.25 Gyr ago that was centrally concentrated but with enhanced star formation galaxy wide (or at least IFU footprint wide). In the central fiber the burst strength was as much as 60%, with the amplitude fading away with distance. Is this percentage plausible? First, just summing over all fibers the mass formed during and after the initial burst is about 1/3 of the present day stellar mass. This probably overestimates the burst contribution since the galaxy extent is considerably larger than the IFU and the outer regions of the progenitors likely experienced less enhanced star formation, but even 1/3 is consistent with the typical neutral hydrogen content of present day late type spirals. I will look at the merger simulation literature in more detail later, but simulations do predict that much of the gas is funneled into the central region of the merger remnant where it’s efficiently converted into stars, so this is broadly consistent with theory. But second, I have some evidence from simulated star formation histories that the presence of a late time burst destroys evidence of the pre-burst history, and this leads to models underestimating the early time mass growth. This could have an ≈0.1 dex effect on the total stellar mass estimate.

Most of the SFH models show more or less monotonically declining star formation rates after the initial burst, but the central fiber shows more complex late time behavior with several secondary peaks in star formation rate. It’s tempting to interpret these literally in hopes of constraining the timeline of the merger. In the next post I’ll look at this in more detail. One limitation of this set of models is the coarse time resolution. I only used half the available SSP ages from the EMILES/BaSTI library — the critical 1.25Gyr model in particular is 350Myr wide, which is a significant fraction of the expected timescale of the merger. Next time I’ll look at the results of some higher time resolution model runs.

I look at a large and still growing number of quantities derived from the models. For visualization it’s sometimes useful to create maps like the ones I showed in the first post on this galaxy. Many interesting properties of this particular galaxy are fairly radially symmetric with strong radial trends, so simple scatter plots are most effective.

I posted several graphs in a post on the current GZ Talk. Since those posts I’ve learned how to do Voronoi binning by translating Cappellari’s Python code to R. I also noticed that some of the EMILES spectra are truncated in the red within the range I was using for SED fitting, so I trimmed it to a wavelength range with strictly positive model fluxes (currently rest frame wavelengths 3465-8864 Å) and reran the models on the binned data. Binning has little effect on the results since out of 183 fiber/pointing combinations I ended up with 179 binned spectra. There are small changes due to the reduced wavelength range.

One of the more interesting plots I posted there was the estimated stellar mass density (in log solar masses/kpc^2) against radius (measured in effective radii; the effective radius from the drpall catalog is 5.48″ or about 3.2 kpc). The updated plot is below. The blue line is a single component Sersic model with Sersic index n=3.35 from a nonlinear least squares fit using the R function nls() (no Bayes this time). This is close enough to a deVaucouleurs profile (n=4) and definitely not disky (n=1). I will someday get around to comparing this to a photometric decomposition — they should be similar, although this fit is undoubtedly affected by the low spatial resolution of the IFU data.

The other plots I posted there included the star formation rate density (in \(M_\odot/yr/kpc^2\), averaged over an arbitrarily chosen 100 Myr) against radius and Hα luminosity density (uncorrected for internal attenuation) against radius:

The trends look rather similar, which leads to the obvious (SFR density vs. luminosity density):

This is one of the more remarkable and encouraging results of these modeling exercises — the star formation rate estimates come entirely from the stellar component fits and are for an order of magnitude longer time scale than is probed by Hα emission. Optical SED fitting is rarely considered suitable for estimating the recent star formation rate, yet we see a tight linear relationship between the model estimates and Hα luminosity, one of the primary SFR calibrators in the optical, that only begins to break down at the lowest luminosities. In the right pane Hα is corrected for attenuation using the Balmer decrement with an assumed intrinsic Hα/Hβ ratio of 2.86. The straight line in both graphs is the Hα-SFR calibration of Moustakas, Kennicutt, and Tremonti (2006) with a 0.2 dex intercept shift to account for different assumed IMFs. At the well constrained high luminosity end most of the points appear to lie above the relation, which could indicate that recent star formation has declined relative to the 100Myr average (or of course it could be a fluke).

Two other measures of star formation rate I track are specific star formation rate, that is SFR divided by stellar mass (which has units of inverse time), and relative star formation rate, SFR divided by the average SFR over cosmic time. Trends with radius (note that these and just about all other quantities I track are scaled logarithmically):

On longer time scales I track the stellar mass fraction in broad age bins (this has been the customary sort of quantity reported for some years). Here is the intermediate age mass fraction (0.1 < T ≤ 2.3 Gyr), which basically measures the burst strength for this galaxy:

Oddly, the peak burst strength is somewhat displaced from the photometric center.

For completeness, here are trends of the modeled optical depth of stellar attenuation and intrinsic g-i color (a long wavelength baseline color serves as a rough proxy for specific star formation rate):

This is very similar to the color trend calculated from the spectral cubes I showed in the first post.

Finally, in recent months I have begun to track several “strong line” gas phase metallicity indicators. Here is an oxygen abundance estimate from the O3N2 index as calibrated by Pettini and Pagel (2004). Astronomers present “heavy” element abundances in the peculiar form 12 + log(O/H) where O/H is the number ratio of Oxygen (or other element) to Hydrogen.

Unlike just about everything else there’s no clear abundance trend here, although the precision drops dramatically outside the central few kpc. For reference the solar Oxygen abundance (which is surprisingly unsettled) is around 8.7, so the gas phase metallicity is around or perhaps a little above solar.

So, to summarize, there was a powerful burst of star formation that began ~1Gyr ago that was clearly triggered by a gas rich major merger. The starburst was centrally concentrated but enhanced star formation occurred galaxy wide (perhaps not uniformly distributed). Star formation is ongoing, at least within the central few kiloparsecs, although there is some evidence that it is fading now and is certainly much lower than the peak.

The stellar mass distribution is elliptical-like although perhaps not quite fully relaxed to a new equilibrium. As we saw in the previous post there is no evidence of regular rotation, although there is a rotating kinematically decoupled core that coincides in position with the  region where the starburst was strongest.

Next post I’ll look at some of the recent merger literature, and also show some results from higher time resolution modeling.

 

As I mentioned in the first post I’ve been doing star formation history (SFH) modeling by full spectrum fitting of galaxy spectra for several years. I’m going to quickly review what I do as a prelude to the next several posts where I plan to look at my proto-K+A galaxy in more detail. I’ll probably return to modeling details in the future. I still haven’t reached the point of having version controlled code let alone “publishing” it, but that’s not far off (I hope).

Full spectrum fitting is actually pretty simple, especially since other people have already done the hard work. The basic observational data is a set of monochromatic fluxes \(g_{\lambda}\) with estimated uncertainty at each point \(\lambda\) (\(\lambda\) is just a discrete valued index here) \(\sigma_{\lambda}\) or precision (inverse variance) \(p_{\lambda} = 1/\sigma^2_{\lambda}\) . For now let’s assume that observed fluxes are independently gaussian distributed with known variances (yes, this is a simplification).

The second main ingredient is a set of template spectra \(T_{\lambda, n}, n = 1, \ldots, N\), which nowadays usually comprises model spectra from a library of simple stellar population (SSP) models. If the physical parameters of the templates span those of the galaxy being modeled its spectrum would be fit with some nonnegative combination of the template spectra:

\(g_{\lambda} = \sum_{n=1}^N \beta_n T_{\lambda, n} + \epsilon_{\lambda}, ~~\beta_n \ge 0~~ \forall n\)

 

and by assumption the errors \(\epsilon_{\lambda}\) are gaussian distributed with mean 0 and standard deviation \(\sigma_{\lambda}\).

The spectroscopy  we deal with is high enough resolution to be blurred by the constituents of the galaxy, so an additional ingredient is the “line of sight velocity distribution” (LOSVD), which must be convolved with the templates. Finally, it’s common to include additive or multiplicative functions to the model to allow for calibration variations or to model attenuation. I hopefully assume calibrations are well matched between observations and models and only attempt to model attenuation, usually with a Calzetti relation parametrized by optical depth. So, the full model looks like:

\(g_{\lambda} = \sum_{n=1}^N \beta_n  (\mathcal{L}(v_n) \ast T_{\lambda, n}) \exp(-\tau_{V, n} A_{\lambda}) + \epsilon_{\lambda}, ~~\beta_n \ge 0~~ \forall n\)

 

where the \(v_n\) are (possibly vector valued) parameters of the LOSVD and \(\tau_{V,n}\) is the dust optical depth at V. These are nominally indexed by n to indicate that there may be different values for different model constituents. In practice I only consider a few distinct components.

An important insight of Cappellari and Ensellem (2004) that was dismissed in two sentences is that conditional on the parameters of the LOSVD and any multiplicative factors the subproblem of solving for the coefficients \(\beta_n\) is a linear model in the parameters with nonnegativity constraints, for which there is a fast, deterministic least squares algorithm due originally to Lawson and Hanson (1974). This greatly simplifies the nonlinear optimization problem because whatever algorithm is chosen for that only needs to optimize over the parameters of the LOSVD and optical depth. I currently use the general purpose optimization R package Rsolnp along with nnls for SED fitting. This is all similar to Cappellari’s ppxf but less versatile; my main interest is in star formation histories and I’ve found a better way to estimate them in a fully Bayesian framework. These days I mostly use the maximum likelihood fits as inputs to the Bayesian model implemented in Stan.

I have used several SSP libraries, but the workhorses recently have been subsets of the EMILES library with BaSTI isochrones and Kroupa IMF, supplemented with some unevolved models from BC03. One thing I do that’s not unprecedented but still somewhat unusual is to fit emission lines simultaneously with the stellar library fits. I do this by creating fake emission line spectra, one for each line being fit (usually the lower order Balmer lines and selected forbidden lines).

The models that I’m likely to discuss for the foreseeable future use separate single component gaussian fits to the stellar and emission line velocity distributions. I assume the velocity offsets calculated in the first step of the analysis are accurate, so the mean of the distributions is set to 0. I also use a single component dust model, usually with Calzetti attenuation relation. I don’t try to tie the relative strengths of emission lines, so attenuation is only directly modeled for the stellar component. These assumptions are obviously too simple for galaxies with broad emission line components or interacting galaxies with overlapping stellar components. For now I plan to ignore those.

One interesting property of the NNLS algorithm for this application is that it always returns a parsimonious solution for the stellar component, with usually only a handful of contributors. An example of a solution from a fit to the central spectrum of the proto-K+A galaxy I introduced a few posts ago is shown below. The grey bars are the optimal contributions (proportional to the mass born in each age and metallicity bin) from the elements of the SSP library. This particular subset has 3 metallicity bins and 27 age bins ordered from low to high.

Only six (actually a seventh makes a negligible contribution) SSP model components are in the optimal fit, out of 81 in this particular library. Usually in statistics sparsity is considered a highly desirable property, especially in situations like this with a large number of explanatory variables. Considerable effort has been made to derive conditions where non-negative least squares, which has no tuning parameters, automatically produces sparse solutions (see for example Meinshausen 2013). However in this context sparsity is maybe not so desirable since we expect star formation to vary continuously over time in most cases rather than in a few short bursts as would be implied here: as a Bayesian my prior assigns positive probability density to star formation at all times.

But in the absence of a parametric functional form for star formation histories a fully Bayesian treatment of SED modeling is an extremely difficult problem. The main issue is that model spectra are highly correlated; if we could visualize the N dimensional likelihood landscape we’d see a long, curving, steep sided ridge near the maximum. I’ve tried a number of samplers over the years including various forms of Metropolis-Hastings and the astronomers’ favorite EMCEE and none of them are able to  achieve usable acceptance rates. A few years ago I discovered that the version of “Hamiltonian Monte Carlo” (informally called the No U-Turn Sampler, or NUTS) implemented in the Stan modeling language is capable of converging to stationary distributions in a tractable amount of time. In the plot below the red dots indicate the range of values taken by each SSP model contribution in the post-warmup sample from a Stan run. Now all time and metallicity bins are represented in the model in amounts that are at least astrophysically plausible. For example instead of a single old component as in the NNLS fit the pre-burst history indicates continuous, slowly declining star formation, which is reasonable for the likely spiral progenitors.

A trace plot, that is a sequential plot of each draw from the sampler, is shown for a few parameters below — these are the single old component in the maximum likelihood fit and the two in the surrounding age bins. Samples spread out rather quickly from the initialized values at the MLE, reaching in this case a stationary distribution about 2/3 of the way through the adaptation interval shown shaded in gray; I used a rather short warmup period of just 150 iterations for this run (the default is 1000, while I generally use 250 which seems adequate). There were 250 post warmup iterations for each of 4 chains in this run. No thinning is needed because the chains show very little autocorrelation, making the final sample size of 1000 adequate for inference. Most physical quantities related to the stellar contributions involve weighted sums of some or all coefficients, and these typically show little variation between runs.

Posterior marginal distributions often have very long tails as shown here. This can be a problem for MCMC samplers and Stan is quite aggressive in warning about potential convergence issues. The quantities that I track tend to have much more symmetrical distributions.

Fits to emission lines tend to behave differently. If they’re actually present the MLE will have positive contributions (the coefficients are proportional to emission line fluxes). On the other hand in a passively evolving galaxy with no emission lines the MLE estimate is typically sparse, that is coefficient estimates will tend to be 0 rather than spuriously small positive values.

Sample values from Stan tend to stay centered on the ML fits (Hα and the surrounding [N II] lines are shown here):

Marginal posterior distributions tend to be symmetrical, close to Gaussian, and uncorrelated with anything else (except for the [O II] 3727-3729Å doublet, which will have negatively correlated contributions).

One thing that’s noteworthy that I haven’t discussed yet is that usually all available metallicity bins make contributions at all ages. I haven’t been able to convince myself that this says anything about the actual distribution of stellar metallicities, let alone chemical enrichment histories. In recent months I’ve started looking at strong emission line metallicity indicators. These appear to give more reliable estimates of present day gas phase abundances, at least when emission lines are actually present.