bayesian modeling – Page 2

2022-04-26

Adding emission lines to non-parametric kinematic models

This is still experimental and I’m not sure how much farther I’ll pursue it, but I tried a straightforward way to add emission lines to non-parametric line of sight velocity distribution (LOSVD) models. The idea is simple enough: model the line profiles directly using Stan’s simplex data type with each modeled line represented by a vector of mostly zeroes and with the simplex centered on the line’s rest wavelength. Although not essential I’m assuming I will have estimates of redshift offsets for each fiber or spaxel in a MaNGA data file (RSS or cube), so any additional offsets should be small. I’ve chosen to ignore the fact that the discretized line profiles will differ between lines because their central wavelengths will fall at different points within their assigned wavelength bins. Also, different lines could arise from kinematically distinct regions, which is not uncommon in galaxies with broad line AGNs. The obvious solution to this is to allow multiple line profiles. For these initial exercises I’m using a single line profile for all modeled lines (I fit 18 at present). As I’ve done since I started these modeling exercises I am fitting emission lines and stellar contributions simultaneously, with the stellar part represented by a small set of eigenspectra derived from my usual EMILES based library.

Below the “fold” I’ve included the Stan code in it’s current (but certainly not final) form. About half of the code for modeling the stellar LOSVD, is adopted from the original version that I wrote about last year. The emission line model portion takes advantage of an odd feature of Stan, namely the ability to store a matrix in sparse form and perform one specific operation — matrix multiplication with a vector. I still haven’t figured out the particular matrix representation used, so I just create a dummy matrix for the emission lines in the transformed data section and extract the two vectors describing the positions of the non-zero elements of the matrix. In the model section the simplex vector representing the line profiles is repeated as many times as there are lines to fill in the non-zeros.

Another thing to note is that Stan doesn’t know how to work with missing data. In general there will be gaps in spectra that were masked for some reason, while the input templates must be complete over the covered spectrum (plus a few extra at each end for the convolution). This requires a bit of housekeeping that’s mostly done in the R code that sets up the input data.

In the models I ran last year I ignored dust reddening since I didn’t expect it to be significant in the passively evolving Coma cluster galaxies I tested them on. In general reddening isn’t ignorable though and there’s the potential for template mismatch without some allowance for it. For now I inserted the function I use for “modified Calzetti” attenuation but it’s not actually used. In the one test I tried including attenuation in the model significantly increased execution time. I will probably look at using a multiplicative polynomial for the same purpose. A final thing to note is there are no explicit priors for the two simplex vectors that represent the stellar velocity convolution kernel and emission line profile, which means that the priors default to a maximally diffuse dirichlet distribution. This turns out to be an important issue that I will discuss further below.

I’ve run sets of models for a small sample of galaxies so far. so lets look at some results. For these runs I used 500 warmup and 500 post-warmup iterations per chain with 4 parallel chains for a total post-warmup sample size of 2000. The stellar templates are represented by 6 eigenspectra created from singular value decompositions of my standard EMILES based SSP library. I currently fit 18 emission lines — Balmer lines from Hα through Hζ and a selection of the stronger forbidden lines. Model runs typically take 2-3 minutes for sampling on my old 4th generation Intel core I7 based PC. This can undoubtedly be reduced by factors of at least several with multithreading.

In this post I’m going to look in some detail at results for Markarian 888, which was the main topic of my last post. Recall this is one of Schawinski’s “blue early type galaxies” that turns out to have obvious spiral-like structure in its inner region and clear evidence for a relatively recent merger.

First, the graph below shows the central fiber spectrum, and below that the results of a model run for the stellar velocity convolution kernel. From those I can calculate velocity offsets¹

\( \delta v = 69 \times \sum\limits_{k = -\lfloor K/2 \rfloor}^{\lfloor K/2 \rfloor} k p_i \)

and velocity dispersions

\(
v_{disp} = 69 \times \sqrt{\sum\limits_{k = -\lfloor K/2 \rfloor}^{\lfloor K/2 \rfloor} k^2 p_i – \delta v^2}
\)

for each draw from the posteriors, and these are shown as histograms in the middle and right graphs. This was a high signal to noise spectrum with prominent emission and absorption lines, a favorable situation for this modeling exercise and in fact the results look very promising with an especially well determined distribution for the emission lines. In summary, the posterior mean stellar velocity offset was 25 km/sec. with a 95% credible interval of (19.8, 31.3), while the corresponding values for emission lines were 24.3 and (22.8, 25.9), so the credible interval for emission lines lies entirely within that for stars.

The corresponding values for velocity dispersions on the other hand differ quite a lot: 141 km/sec. for stars with a 95% credible interval of (133.6, 150.3) versus 109 km/sec and (106, 111) for gas. I’ll say some more about this below, but I think this is already a sign that the second moments (at least) of these velocity distributions need to be treated with caution.

mrk888_centralspec_losv — Markarian 888 – mangaid 9894-3703 Central fiber spectrum, modeled stellar and ionized gas losvd

After running the code for every binned spectrum in the IFU I get stellar and gas velocity maps as shown below. These look similar to each other and to the map in the last post based on my hybrid red shift offset fitting routine, although a closer look will show a systematic decoupling of gas and stellar velocities.

mrk888_st_em_losvd — Stellar and ionized gas velocity maps for Mrk 888 (MaNGA plateifu 9894-3703)

I often take a look at the official MaNGA data analysis pipeline results, usually through the “Marvin” interface. Generally my results look similar for any given quantity, but it’s hard to tell how closely since we’re looking at different things (binned RSS spectra vs. cubes, for one) and I’m really just doing a visual comparison. As more or less a one off exercise I decided to compare velocity fields from the cubes. First, here are stellar and Hα maps from the DAP. These are from the “HYB10-MILESHC-MASTARSSP” sequence, which refers to the binning strategy and template inputs.

mrk888_vel_maps_per_dap — Mrk 888 – stellar and gas velocity fields derived from non-parametric LOSVD models

Next are my stellar and ionized gas velocity maps. Evidently these are considerably “noisier” in appearance, but similar overall. The noisiness may be due to different binning strategies: I modeled velocity distributions for every spaxel that met a minimal signal to noise requirement.

mrk888_vel_maps_per_moi — Mrk 888 – stellar and gas velocity fields derived from non-parametric LOSVD models

For a more quantitative comparison here are scatterplots of the mean stellar and gas velocity offsets from my runs against the MaNGA DAP. Both of these show slight tilts, that is the slopes are less than one, and small offsets. I’m tempted to suspect systematic errors somewhere, but I haven’t found any yet. And I have no hypothesis about the decidedly non-random behavior of the gas velocity plots. So there are still questions to be answered, but the magnitude of differences is no more than about ±1 wavelength bin (69 km/sec). For my main purposes this is good enough.

mrk888_vel_comp — Mrk 888 – scatter plot of mean stellar and gas velocities – my models vs. DAP

One more pair of graphs for this galaxy. In the initial maximum likelihood fit to a spectrum I estimate stellar and gas velocity dispersions with a single component gaussian. The graphs below plot the velocity dispersions from the Bayesian models calculated by the formula above against the ML fit values. The stellar velocity dispersions behave pretty much the same as I noted before: there’s a weak positive correlation between the non-parametric models and the maximum likelihood estimates with the former generally offset to higher values. The gas velocity dispersion estimates show an even larger discrepancy along with an apparently nonlinear trend.

mrk_888_vdisp — Mrk 888 – velocity dispersions from Bayesian nonparametric losvd models vs. estimates from maximum likelihood fitting

More briefly now, since this got longer than I wanted. I also ran models for the “Zoogems” target NGC 810, which I noted two posts ago has interesting kinematics, in particular a rapidly rotating disk that isn’t evident in the stellar velocity map from the DAP. From the maps below we see that the rapidly rotating component is the ionized gas, which has a mean rotation velocity as much as 150 km/sec larger than the stellar component.

ngc_810_st_em — NGC 810 – mean stellar and gas velocity maps

ngc_810_em-st — NGC 810 – difference between mean gas and stellar velocities

Finally, here are two more graphs of the kinematic models for the spectra with the largest positive and negative differences between stellar and gas velocity (these are the reddest and bluest patches in the map above. These show the spectra in the top row, followed by the posterior distribution of the stellar LOSVD and distributions of the first two moments, then the same for the emission line profiles.

ngc810_175_spec_los — NGC 810 – sample MaNGA spectrum, non-parametric stellar LOSVD and emission line profile

ngc810_90_spec_los — NGC 810 – sample MaNGA spectrum, non-parametric stellar LOSVD and emission line profile

The second set of graphs illustrates a significant problem: when there’s poor data — in this case no detectable emission lines — the models tend to return the prior, which for now is a maximally diffuse dirichlet. This will bias second moment estimates to the high side and can lead to spurious correlations, for example the increase of velocity dispersion with radius that I noted in the last post on this topic.

At the same time this spectrum is in the region of overlap between the main galaxy and its companion, and two peaks in the stellar velocity distribution are clearly seen at about the right velocity offset (~340 km/sec).

A short to do list:

Investigate informative priors.
Allow for dust reddening. I tried using a modified Calzetti attenuation relation in the code (note its presence in the code below, although it’s not used in the model), but it adds too much to execution time. A low order multiplicative polynomial might suffice.
Investigate systematics between my estimates and those in the DAP.

2022-01-152022-01-17

NGC 810 – interesting kinematics in a Zoogems and MaNGA target

The final release of SDSS MaNGA went public back in early December as promised, and I’ve spent the last month or so of my hobby time looking for manageable sized samples of interesting galaxies. One sample I looked at out of curiosity was the Zoogems target list, which is an HST gap filler imaging program with about 300 galaxies selected (mostly) by Galaxy Zoo volunteers. It turns out there are 11 targets with MaNGA data, of which 5 have been observed by HST so far.

As can be seen from the thumbnails above this is a pretty diverse lot, with several in progress mergers and merger remnants, some normal looking spirals at least two of which were from Masters’ red spirals sample, and 3 of Schawinski’s blue early type galaxies. Only one of those 3 has HST imaging so far (number 8 in the thumbnails above), although there are a surprising 24 blue ellipticals in the final MaNGA release out of 215 in Schawinski’s original sample.

Of the 5 Zoogems galaxies that have been observed so far the one that caught my eye as deserving an early look was number 3 in the top row, NGC 810, an apparent elliptical with an unusual dust lane that’s almost perfectly aligned with the minor axis. There are also hints of shells indicating a likely merger sometime in the past.

NGC 810 HST ACS, proposal id 15445, PI Keel

The MaNGA data, which is new in DR17, only covers the central part of the galaxy with the companion just photobombing the edge. A larger IFU would have been nice for this observation, but the data quality is better than average in terms of nominal signal to noise. I was able to use all 183 fiber/position combinations in the RSS file without binning.

NGC 810 – plateifu 9514-6101 – MaNGA IFU footprint

The first step in the analysis process after loading the data is to estimate redshift offsets from the system redshift for each spectrum, and from that it’s straightforward to calculate a velocity field, which in this galaxy looks like¹:

NGC 810 (plateifu 9514-6101) – losvd estimated from cube

It turns out the redshift assigned to this system was that of the companion galaxy, which was the only SDSS spectroscopic target in the immediate vicinity and is evidently blueshifted by ~350 km/sec from the target. What’s more interesting though² is the apparent rapidly rotating disk that’s aligned with and somewhat thicker than the dust lane. There may also be overall prolate rotation outside the disk although the presence of the companion makes it hard to tell based solely on visual inspection. In hopes of separating out multiple velocity components I returned to the non-parametric line of sight velocity distribution models that I wrote some posts about last year. Unlike the ones I practiced on previously this galaxy has non-negligible amounts of emission, at least in the central region, so I just temporarily masked the regions around the emission lines that I fit. That results in a pure stellar velocity distribution. The results were a bit surprising:

NGC 810 (plateifu 9514-6101) (L) velocity field estimated from RSS file (M) stellar velocity offsets (R) net stellar velocity

In the left pane above is the velocity field from the RSS data, with the system redshift adjusted to the IFU center. For the LOSVD models I set the adopted redshift of each spectrum to the system redshift plus the offset calculated previously. Now I had hoped to be able to cleanly separate the contribution of the companion from that of the main galaxy, which so far I haven’t been able to do. But what I did find that was unexpected is that the average stellar velocities in the disk partially offset the original measurements (middle pane), so the net stellar velocity field shows a much more slowly rotating stellar disk.

As I’ve written before I use a set of 15 eigenspectra from a principal components analysis of some tens of thousands of SDSS spectra that I performed some years ago for redshift offset estimation. Those galaxies were of all types and include systems both with and without significant emission. The redshift estimation routine just does straightforward template matching and returns a single value for the best fitting offset. Since the templates encode information about both emission and absorption lines that estimate could be most applicable to the ionized gas, stars, or some combination. In this case it’s possible emission lines were driving the original fits, implying the gas and stars in the disk are kinematically decoupled. I have not verified that though.

Another issue I noticed is that the stellar velocity field from the official MaNGA Data Analysis Pipeline looks rather different from mine, with barely a hint of a kinematically distinct disk. This wasn’t really evident in the Marvin webpage, which makes some really unfortunate choices for color palettes. So here is the same data rendered with a more nearly perceptually uniform rainbow palette³:

NGC 810 (plateifu 9514-6101) Stellar velocity field per MaNGA DAP

I decided to re-run my LOSVD modeling code on the RSS data, this time setting the redshift offset to 0 for each fiber, so this is now measuring velocities relative to the overall system velocity. I also used a larger convolution kernel (25 vs. 21 in the first set of runs). The map of the average velocity offsets is:

NGC 810 (plateifu 9514-6101) Stellar losvd from nonparametric model

Although not a perfect match this is somewhat closer to the DAP map. I suspect what’s happening here is that there really are at least two, and more likely 3 distinct kinematic components. I haven’t read the DAP release paper in a while and don’t know exactly how they estimate stellar velocities, but in any case their model just returns a single value for visualization purposes. To see the (possible) complexity of the actual data here are the results for a single fiber with the largest positive velocity offset in the map above. Again, I don’t know how much of the structure in the posterior distribution of the convolution kernel is real, but it’s evident there’s more complexity than is captured in the first two moments shown in the middle and right panes.

Besides the kinematic modeling I did run star formation history models on the full RSS data set using the same tools as in recent posts. I’m not going to discuss them in detail, but some summary maps are worth showing. In the graphic below these are, from top left, stellar mass density, Hα luminosity density uncorrected for attenuation, SFR density (as usual 100Myr average), and stellar dust attenuation.

There’s no sign of a disk in the stellar mass map, which faithfully follows the distribution of light. A disk is visible in Hα and the small amount of recent starformation is also confined to a disk and nuclear region. In the fourth pane I show the modeled stellar dust attenuation, mostly just to demonstrate that this component of the model does capture something of reality.

ngc810_model_summaries — NGC 810 – a sampling of quantities derived from star formation history models

Getting back briefly to the first paragraph, there are 8 post-starburst galaxies from the catalog of Melnick and de Propris and 24 from an ancillary program to observe post-starburst galaxies from various sources that was added for DR17. There’s just one galaxy from the former catalog in the latter set, so that makes 31 total, an easily manageable number. There are also 24 of Schawinski’s blue ellipticals. Of course there are many disk galaxies, far too many for me to look at.

2021-08-102021-08-10

KUG 0859+406 – unravelling the differences between 2018 and 2021 model runs

I did have my old data and model runs of course, in fact they were spread over several directories on two machines. I’m going to refer to it by this catalog designation, KUG standing for the “Kiso survey of Ultraviolet-excess Galaxies.” It’s also a low power radio source with catalog entries in both FIRST and NVSS, and of course it’s in MaNGA with plateifu ID 8440-6104 (mangaid 01-216976).

In my 2018 model runs, which were interesting enough to write 3 posts about, I found this galaxy had undergone an extraordinarily large burst of star formation that began ~1 Gyr ago with locally as much as 60% of the present day stellar mass born in the burst and something like 40% of the mass over the footprint of the IFU. In this years model runs the peak burst fraction was a considerably more modest ~25% and globally barely amounted to a slight enhancement of star formation. The starburst was also much more localized than in the earlier runs:

burst_strengths_compared — Fractional stellar mass in stars between 0.1 and 1.75 Gyr old in 2018 and 2021 model runs

So what happened? First, here is a comparison of modeled star formation histories for the innermost fiber, which got the largest injection of mass in the starburst.

central_sfh_compared — Model star formation histories for central fiber of MaNGA plateifu 8440-6104, 2018 and 2021 model runs

The obvious remark is the double peaked starburst noted back in 2018 (and discussed at some length) has been replaced with a single narrow peak with a slow ramp up and fast decay. The peak SFR is a little larger than before but the total mass in the burst is lower.

I’ve made several changes in model formulation since 2018, of which the most important in the current context is adopting the more flexible “modified Calzetti” attenuation relation that adds an additional slope parameter to the prescription. In the current year model runs a steeper than Calzetti relation is favored throughout the IFU footprint, particularly in the central region where the starburst was strongest:

map_delta — Map of modified Cal;zetti slope parameter δ — MaNGA plateifu 8440-6104

A smaller optical depth of attenuation is also favored throughout:

tauv_compared — Modeled optical depth of attenuation – 2021 runs vs. 2018 MaNGA plateifu 8440-6104

This has a couple predictable consequences. Steeper attenuation will favor an intrinsically bluer, hence younger population while a lower optical depth requires less light, and hence mass in the stellar population. I can test this directly by returning to a model with Calzetti attenuation, and here is the result for the central fiber (this model run is labeled 2021 (c) in the legend below):

mgh_compared — Mass growth histories – 2021 run 2021 run with Calzetti attenuation 2018 run Central fiber of MaNGA plateifu 8440-6104

So, an eyeball analysis suggests about 3/4 of the difference between the 2018 and 2021 runs is due to the modification to the attenuation relation. The other changes I’ve made to the models are to change the stellar contribution parameters from a non-negative vector to a simplex, and at the same time changing the way I rescale the data. In early runs the SSP model fluxes were scaled to make the maximum stellar contribution ≈ 1, while the current models scale both the galaxy and SSP fluxes to ≈ 1 in the neighborhood of V, making the individual stellar contributions approximately the fraction of light contributed. An additional scale factor parameter in the model is used to adjust the overall fit. Assuming I did this right this should have no effect on a deterministic maximum likelihood solution, but with MCMC who knows?

Although the fit to the data looks about the same between the model with and without the attenuation modification the summed log-likelihood is consistently about 1% higher for the modified Calzetti model with no overlap at all in the distribution of likelihood. This suggests the case for a steeper than Calzetti attenuation is a fairly robust result.

ppfits — “Posterior predictive” fits to galaxy flux data – modified Calzetti attenuation vs. Calzetti – central fiber of MaNGA plateifu 8440-6104

The galaxy flux data also changed a little bit. The early runs were on the DR14 release (version 2_1_2 of the MaNGA DRP) while the recent ones used the DR15 release (ver 2_4_3). Most of the calibration differences resemble random noise, but there is some curvature that systematically affects both the red and blue ends of the spectrum and could cause some change in the temperature distribution of the models:

measured_flux_compared — Difference in measured flux from DR14 to DR15 – central fiber of MaNGA plateifu 8440-6104

While the detailed star formation histories changed, quantities that aren’t too model dependent didn’t very much. One example is shown below. Also, the kinematic maps agree with the earlier ones in detail.

halpha_luminosity_compared — Hα luminosity density – 2021 runs vs. 2018 MaNGA plateifu 8440-6104

velocity_field_8440_6104 — Velocity field of MaNGA plateifu 8440-6104 from 2021 model runs. Map interpolated from RSS file spectra.

One input that hasn’t changed are the emiles SSP model spectra, although there have been some procedural changes in how I handle the modeling. Early on I often used a much smaller subset of SSP models with just 27 time bins and 3 metallicities for preliminary modeling, including my first models on the same binning of these data. I also routinely ran 250 warmup iterations with just 250 more per chain. My current standard practice is always to use the largest emiles subset with 216 SSP models in 54 time bins and 4 metallicities, and I generally run 750 post-warmup iterations per chain but still with only 250 warmup iterations. This is generally enough and if adaptation fails it is usually fairly obvious. The small sample size of the earlier runs mostly effects the precision of inferences rather than means.

To conclude for now, my speculation about whether it might be possible to say something about the timing of critical events in a merger from the model star formation history was too optimistic. On a positive note though both sets of model runs retrodict that coalescence occurred at a lookback time around 500Myr ago, which is consistent with the fact that tidal tails and other merger signatures are clearly visible even in SDSS imaging. Both sets of model runs also have that odd uptick in star formation at 30Myr in the central fiber. And while the difference in burst mass contributions is a little disconcerting the current runs are more consistent with the likely gas content of ordinary spiral galaxies.

This example illustrates another well known “degeneracy” among attenuation (and adopted attenuation relation), mass, and stellar age. Whether I’ve broken the degeneracy by adopting the more flexible attenuation prescription described some time ago remains to be validated.

2021-07-142021-07-14

Do the details of the line of sight velocity distribution matter for star formation history modeling? Probably not much.

I decided to try a set of models for one galaxy – NGC 4889 (with MaNGA plateifu 8479-12701), which had the highest overall velocity dispersion of the Coma sample I’ve been discussing in the last several posts. It also has some evidence for multiple kinematic components which isn’t too much of a surprise since it’s one of the central cD galaxies in Coma. The SSP model spectra fed to the SFH models were preconvolved with the element wise means of the LOSVD convolution kernels from the velocity distribution modeling exercise. Again, this is an expedient to avoid what could otherwise be prohibitively computationally expensive. The models I ran were the same as described back in this post — these ignore emission but do model dust attenuation with the usual modified Calzetti attenuation relation.

To get quickly to the results, here are model star formation histories compared to the previous runs that used the full model in its current form. Usually I like these plots of results from all spectra in an IFU, but in this one all 381 spectra met my S/N criterion, so the plot is pretty crowded. You really need to see it live on a 4K monitor to see the details.

NGC 4889 (MaNGA plateifu 8479-12701) Model star formation histories for all spectra, runs with non-parametric LOSVD vs. single gaussian stellar velocity dispersions

Well it’s pretty hard to see but differences in model SFH’s are mostly in the youngest age bins, which are very poorly constrained anyway in these presumably passively evolving galaxies. Here’s a closer look at a single model run that had the largest difference in estimated stellar mass density (more on this right below) of about 0.19 dex:

mgh_comparison_8479-12701 — NGC 4889 (MaNGA plateifu 8479-12701) Model mass growth histories for a single spectrum – runs with non-parametric LOSVD vs. single Gaussian stellar velocity dispersion

So, the difference in star formation histories was slower mass build up between about 12-5 Gyr look back times in the second run, which was responsible for the lower current day stellar mass density. How this resulted from the choice of LOSVD is not at all obvious.

Let’s look at a few summary results. First, the model stellar mass surface densities:

sigma_mstar_comparison_8479-12701 — NGC 4889 (MaNGA plateifu 8479-12701) Model Σ_M* – runs with non-parametric LOSVD vs. single Gaussian stellar velocity dispersion

These fall on an almost exactly one to one relation with a few hands full of outliers. Oddly these are mostly in the higher signal to noise area of the IFU (i.e. near the center).

Results for star formation rate density and specific star formation rate are even more consistent between runs, with essentially no differences larger than the nominal 1 σ error bars.

sigma_sfr_comparison_8479-12701 — NGC 4889 (MaNGA plateifu 8479-12701) Model Σ_sfr – runs with non-parametric LOSVD vs. single Gaussian stellar velocity dispersion

ssfr_comparison_8479-12701 — NGC 4889 (MaNGA plateifu 8479-12701) Model SSFR – runs with non-parametric LOSVD vs. single Gaussian stellar velocity dispersion

One problem I encountered was that I had to re-run some models either for technical reasons or because of obvious convergence failures. I suspect there could have been some convergence issues in both sets of runs and am slightly worried that could be the source of the few differences in summary measures seen. Oddly, there were almost no suspicious convergence diagnostics in either set of runs (once the latter were run to satisfactory conclusion), and Stan is quite aggressive about reporting possible convergence issues.

Anyway, modeling kinematics remains an interesting topic to me, but it seems somewhat decoupled from modeling star formation histories. Right now I’m waiting for the final SDSS data release to decide what projects I want to tackle.

I’m going to end with a couple of asides. First, I recognize that all of these error bars are overoptimistic, maybe by a lot. The main reason, I think, is that I treat the flux values as independent which they clearly are not¹, which effectively results in overestimating the sample size. One possible partial solution is to allow the flux uncertainties to vary from their nominal values by, for example, a factor > 1. This would involve adding as few as one parameter to the models, which is something I’ve actually tried in the past. I may relook at that.

One interesting feature of the previous two graphs is the rather obvious systematic trend with radius of both SFR density and specific star formation rate, as shown more directly below taken from the first set of model runs:

sfr_ssfr_d_8479-12701 — NGC 4889 (MaNGA plateifu 8479-12701) Σ_SFR and SSFR vs. distance from IFU center

Are these real trends? I don’t know, but I don’t see an obvious reason why they might be spurious features of the models. In normal star forming galaxies I encounter trends with radius in both directions and sometimes no trend at all.

As a final and related aside there was a paper by Sedgwick et al. that showed up on arxiv not long ago that presented estimates of star formation rates of early type galaxies from observations of core collapse supernovae carefully matched to host galaxies with high confidence morphological classifications. To oversimplify their conclusions they found that typically massive ellipticals might have specific star formation rates ∼ 10^-11 / yr, which is somewhat higher than usually supposed. As I mentioned in my last post my models will always have some contribution from young stars and I typically get central estimates of SSFR > ∼ 10^-11.5 even in galaxies with no hint of emission (as is the case with this Coma sample). This particular galaxy has a total stellar mass within the IFU of ∼ 10^11.5 M_⊙ , so it could be forming stars at a rate of ∼ 1 M_⊙ / yr.

Well, I think I have one more post to write before the SDSS DR17 release.

2021-05-27

A little more on nonparametric line of sight velocity distribution modeling – is regularization needed?

A while back I completed LOSVD models for 3348 spectra in the 33 Coma Cluster galaxies I chose for an initial sample. I want to look briefly at a couple summary statistics from the models. For each sample in the posterior I calculate first and second moments of the velocity distribution, and from those I calculate means and standard deviations of velocity offsets and velocity dispersions. Recall from the earlier posts that wavelengths are in the spectrum rest frame with peculiar velocities estimated by fitting a set of eigenspectra to the data. The first plot below is the mean velocity offset with errorbars ± 1 standard deviation versus the signal to noise in each spectrum. The horizontal line is the overall median of 5.25 km/sec. This is less than 1/10th the width of a wavelength bin and, I think, consistent with the absolute wavelength calibration accuracy claimed by the MaNGA team. From an eyeball analysis there doesn’t appear to be any trend in the mean with S/N, but the dispersion in estimates increases as the data quality gets worse. This seems unlikely to be a real problem. Two outliers cut off from the graph had mean velocity offsets of ≈ 200 km/sec, which is only 3 wavelength bins.

voff_v_snr — Mean velocity offset and standard deviation of mean vs. spectrum S/N for 33 passively evolving Coma Cluster galaxies. Solid horizontal line is median offset.

Possibly more concerning is that the mean velocity dispersion also shows a systematic trend with S/N:

vdisp_v_snr — Mean velocity dispersion and standard deviation of mean vs. spectrum S/N for 33 passively evolving Coma Cluster galaxies.

Of course there is a wide range of velocity dispersions since there’s a range of galaxy masses in this sample, but each galaxy shows a similar trend. This matters because, of course, S/N correlates strongly with distance from the nucleus which is usually the IFU center and this in turn leads to potentially spurious correlations with distance.

MaNGA spectrum S/N vs. distance from IFU center in kpc. — 33 passively evolving Coma Cluster galaxies

For example here are results for the cD galaxy NGC 4889 with MaNGA plateifu 8479-12701. This has the highest overall velocity dispersion in the sample by a fair margin as seen in the plots above. Plotting the estimated velocity dispersion against distance from the IFU center shows an apparent systematic increase¹:

NGC 4889 (MaNGA plateifu 8479-12701) – mean velocity dispersion vs. radius per nonparametric LOSVD fit.

Now this could be a real trend: my recollection is the outskirts of cD galaxies have essentially the velocity dispersion of the clusters they live in, which is taken as strong evidence that they grow by cannibalism. But the IFU only covers the inner part of this galaxy and in any case the confounding effect of the correlation between S/N and modeled velocity dispersion prevents any conclusion. Just for comparison here is the same relation using the single gaussian estimates from the preliminary maximum likelihood fits:

NGC 4889 (MaNGA plateifu 8479-12701) – velocity dispersion vs. radius per maximum likelihood fit.

So, if we’re interested in kinematics or dynamics some regularization with an informative prior seems advisable. What’s not so clear is whether this matters in SFH modeling, which as I’ve said many times is what most interests me. There’s one way to find out…

2021-01-182021-01-18

Arxiv notes: “BAYES-LOSVD: a bayesian framework for … line-of sight velocity distributions”

This paper by J. Falcón-Barroso and M. Martig that showed up on arxiv back in November interested me for a couple reasons. First, this was one of the first research papers I’ve seen to make use of the Stan language or indeed any implementation of Hamiltonian Monte Carlo. What was even more interesting was they were doing exactly something I experimented with a few years ago, namely estimating the elements of convolution kernels for non-parametric description of galaxy stellar kinematics. And in fact their Stan code posted on github has comment lines linking to a thread on discourse.mc-stan.org I initiated where I asked about doing discrete convolution by direct multiplication and summation.

The idea I was pursuing at the time was to use sets of eigenspectra from a principal components decomposition of SSP model spectra as the bases for fitting convolved model spectra, with the expectation being that considerable dimensionality reduction could be achieved by using just the leading PC’s, which would in turn speed up model fitting. In fact this was certainly the case, although at the time I perhaps overestimated the number of components needed ¹ and I found the execution time to be disappointingly large. Another thing I had in mind was that I wanted to constrain emission line kinematics as well, and I couldn’t quite see how to do that in the same modeling framework. So, I soon turned away from this approach and instead tried to constrain both kinematics and star formation histories using the full SSP library and what I dubbed a partially parametric approach to kinematic modeling, using a convolution kernel for the stellar component and Gauss-Hermite polynomials for emission lines. This works, more or less, but it’s prohibitively computationally intensive, taking an order of magnitude or more longer to run than the models with preconvolved spectral templates. Consequently I’ve run this model fewer than a handful of times and only on a small number of spectra with complicated kinematics.

Falcón-Barroso and Martig started with exactly the same idea that I was pursuing of radically reducing the dimensionality of input spectral templates through the use of principal components. They added a few refinements to their models. They include additive Legendre polynomials in their spectral fits, and they also “regularize” the convolution kernels with informative priors of various kinds. Adding polynomials is fairly common in the spectrum fitting industry, although I haven’t tried it myself. There is only a short segment in the red that is sometimes poorly fit with the EMILES library, and that would seem to be difficult to correct with low order polynomials. What does affect the continuum over the entire wavelength range and could potentially cause a severe template matching issue is dust reddening, and this would seem to be more suitable for correction with multiplicative polynomials or simply with the modified Calzetti relation that I’ve been using for some time.

I’m also a little ambivalent about regularizing the convolution kernel. I’m mostly interested in the kinematics for the sake of matching the effective SSP model spectral resolution to that of the spectra being modeled for star formation histores, and not so much in the kinematics as such.

I decided to take another look at my old code, which somewhat to my surprise I couldn’t find stored locally. Fortunately I had posted the complete code on that discourse.mc-stan thread, so I just copied it from there and made a few small modifications. I want to keep things simple for now, so I did not add any polynomial or attenuation corrections and I am not trying to model emission lines. I also didn’t incorporate priors to regularize the elements of the convolution kernel. Without an explicit prior there will be an implicit one of a maximally diffuse dirichlet. I happen to have a sample of MaNGA spectra that are well suited for a simple kinematic analysis — 33 Coma cluster galaxies from a sample of Smith, Lucey, and Carter (2012) that were selected to be passively evolving based on weak Hα emission, which virtually guarantees overall weak emission. Passively evolving early type galaxies also tend to be nearly dust free, making it reasonable not to model it.

I’ve already run SFH models for all of these galaxies, so for each spectrum in a galaxy I adjust the wavelengths to rest frame using the system redshift and redshift offsets that were calculated in the previous round. I then regrid the SSP library (as usual the 216 component EMILES subset) templates to the galaxy rest frame and extract the region where there is coverage of both the galaxy spectrum and the templates. I then standardize the spectra (that is subtract the means and divide by the standard deviations) and calculate the principal components using R’s svd() function. I again rescale the eigenvectors to unit variance and select the number I want to use in the models. I initially did this by choosing a threshold for the cumulative sum of the eigenvalues as a fraction of the sum of all. I think instead I should use the square of the eigenvalues since that will make the choice based on the fraction of variance in the selected eigenvectors. With the former criterion a threshold of 0.95 will result in selecting 8 eigenspectra, 0.975 results in 13 (perhaps ± 1). In fact as few as 2 or 3 might be adequate as noted by Falcón-Barroso and Martig. The picture below of three eigenspectra from one model run illustrates why. The first two look remarkably like real galaxy spectra, with the first looking like that of an old passively evolving population, while the second has very prominent Balmer absorption lines and a bluer continuum characteristic of a younger population. After that the eigenspectra become increasingly un-spectrum like, although familiar features are still recognizable but often with sign flips.

three_eigenspectra — First three eigenspectra from a principal components decomposition of an SSP model spectrum library

The number of elements in the convolution kernel is also chosen at run time. With logarithmic wavelength binning each wavelength bin has constant velocity width — for SDSS and MaNGA spectra of 69.1 km/sec. If a typical galaxy has velocity dispersion around 150 km/sec. and we want the kernel to be ≈ ±3 σ wide then we need 13 elements, which I chose as the default. Keep in mind that I’ve already estimated redshift offsets from the system value for each spectrum, so a priori I expected the mean velocity to be near 0. If peculiar velocities have to be estimated the convolution kernel would need to be much larger. Also in giant ellipticals and other fairly rare cases the velocity dispersion could be much larger as we’ll see below.

Lets look at a few examples. I’ve only run models so far for three galaxies with a total of 434 spectra. One of these (SDSS J125643.51+270205.1) appears disky and is a rapid rotator. The most interesting kinematically is NGC 4889, a cD galaxy at the center of the cluster. The third is SDSS J125831.59+274024.5, a fairly typical elliptical and slow rotator. These plots contain 3 panels. The left shows the marginal posterior distribution of the convolution kernels, displayed as “violin” plots. The second and third are histograms of the mean and standard deviation of the velocities. For now I’m just going to show results for the central spectrum in each galaxy. These are fairly representative although the distributions in each velocity bin get more diffuse farther out as the signal/noise of the spectra decrease. I set the kernel size to the default for the two lower mass galaxies. The cD has much higher velocity dispersion of as much as ~450 km/sec within the area of the IFU, so I selected a kernel size of 39 (± 1320 km/sec) for those model runs. The kernel appears to be fairly symmetrical for the rotating galaxy (top panes), while the two ellipticals show some evidence of multiple kinematic components. Whether these are real or statistical noise remains to be seen.

centrallosvd_8931_1902 — Distribution of kernel, velocity offset, and velocity dispersion. Central spectrum of MaNGA plateifu 8931-1902

centrallosvd_8933_6101 — Distribution of kernel, velocity offset, and velocity dispersion. Central spectrum of MaNGA plateifu 8933-6101

centrallosvd_8479_12701 — Distribution of kernel, velocity offset, and velocity dispersion. Central MaNGA spectrum of NGC 4889 (plateifu 8479-12701)

There are a couple surprises: the velocity dispersions summarized as the standard deviations of the weighted velocities are at best weakly correlated and generally larger than the values estimated in the preliminary fits that I perform as part of the SFH modeling process. More concerning perhaps is the mean velocity offset averages around -50 km/sec. This is fairly consistent across all three galaxies examined so far. Although this is less than one wavelength bin it is larger than the estimated uncertainties and therefore needs explanation.

hist_velocity_offset — Distribution of mean velocity offsets from estimated redshifts for 434 spectra in 3 MaNGA galaxies.

Some tasks ahead: Figure out this systematic. Look at regularization of the kernel. Run models for the remaining 30 galaxies in the sample. Look at the effect on SFH models.