A little more on star formation history priors

After my not so insightful realization recounted last time that my attempt to modify the prior on star formation histories wasn’t actually doing anything I thought a little further about how to specify one. Gaussians are always popular choices for priors, so why not give them a try? For a first cut I added the following lines to the “transformed data” section of the Stan model:

  vector[nt*nz] norm_sfr;
  
  norm_sfr = (dT .* norm_st)/sum(dT .* norm_st) ;

Then I added a scale parameter sd_sf to the parameters section, and finally in the model section:

    sd_sf ~ normal(0., 1.);
    b_st_s ~ normal(norm_sfr, sd_sf);

Even though the prior allows values of the stellar contributions that are infeasible given the simplex declaration this causes no technical problems: the models sample without complaint and parameters satisfy all constraints. Execution times are comparable to the original model formulation and convergence diagnostics are OK. But, the model runs had some unexpected features. I did a set of model runs for a single MaNGA galaxy from the post-starburst ancillary sample — mangaid 1-201936 (plateifu 10220-3703). Model star formation histories compared to the original model are shown below for the 58 binned spectra:

MaNGA plateifu 10220-3703 (mangaid 1-201936). Model star formation histories with two different SFH priors.

What’s striking here is that several of the spectra in the low S/N outskirts of the galaxy have nearly constant star formation rates with very little sample variation. In other words the models are basically returning the priors. The cause of this behavior isn’t quite clear. Relatively low signal to noise seems to be necessary, but not sufficient since similarly noisy spectra have essentially the same SFH’s as the original model formulation. It also isn’t due to convergence failure because much longer runs with more adaptation iterations show the same behavior. It is possible perhaps that the posteriors are significantly multimodal and Stan is preferentially falling into one of them. Notable also is that the fits to the data measured by log-likelihood are virtually identical even for the runs with the anomalous SFH’s. At the very least this tells us that uncertainties in quantities derived from the models are considerably larger than within model run variations — of course I have always believed this and said so a number of times.

After trying several variations on this theme that either had none at all or undesirable effects on sampling, and after some additional consideration I think that, given the model parametrization, the uniform on the simplex prior for stellar contributions is actually the one I want. That leaves the question of what, if anything, to do about the abrupt jumps in model star formation rates.

One possibility is simply to redefine the endpoints of the age bins to be, say, halfway between nominal SSP ages instead of at the model ages as is my current practice. In the case of the EMILES library this would mean for example that the 3.75, 4, 4.5 Gyr bins would have widths of 0.25, 0.375, 0.5 Gyr instead of the present 0.25, 0.25, 0.5. This involves no change to the actual model runs at all, so most quantities derived from the models are unchanged.

Another solution is to adopt a library with a more uniform age progression. One with approximately equal increments in log age seems preferred. As yet there have been no published updates to the MaStar based SSP libraries mentioned last time, so I’m waiting for them, while still considering generating my own.

I’m going to briefly return to BPASS based models. After that I’m not sure.

Priors on stellar contributions

One thing I’ve wondered about for a while is the extent that priors on SSP model contributions affect modeled star formation histories. Previous experience suggests not much at all unless the prior is highly constraining. To be a little more specific in my current workflow I normalize the SSP model spectra to have average flux values of 1 in (approximately) the V band, and also adjust the galaxy fluxes in the same way. In the Stan model the stellar contribution parameters are declared as a simplex, that is a vector with non-zero elements summing to 1. That makes the parameter values the fractional contributions to the (unreddened) galaxy flux at V. My current working code doesn’t provide an explicit prior for the stellar parameters, but they have an implicit proper prior of a uniform distribution on the appropriate dimensional simplex. More technically the prior is a Dirichlet distribution with all concentration parameters equal to 1. Note the marginal distributions aren’t uniform, but they are all the same. One implication of that is that a typical draw from the prior1Stan doesn’t sample from the prior even for initialization, but of course the prior influences the posterior through Bayes’ rule. will have jumps in the star formation rate at exactly the times where the width of the age bins jump, a problem that I’ve noted several times before.

A possible solution to this problem is simply to alter the prior to encourage smoothly varying star formation rate rather than smoothly varying light contributions. It turns out I’ve been feeding my Stan code the data I need to do that: the initial mass in a given model SSP is inversely proportional to the normalization factor applied to the spectrum, and the star formation rate is just the mass divided by the time interval assigned to the SSP. Both of those quantities are passed as data to the Stan model even though they weren’t used in any way previously. I added just 3 lines to the code to change the prior on the stellar contributions. In the “transfored data” section

norm_sfr = 1.0 ./ (dT .* norm_st);
norm_sfr = norm_sfr / mean(norm_sfr);

defines the reweighting of the prior, which is written in the model section as

    target += dirichlet_lupdf((b_st_s .* norm_sfr)/sum(b_st_s .* norm_sfr)|rep_vector(1., nt*nz)); 

If I did this right the prior is completely agnostic to any star formation history, whereas the previous implicit prior was completely agnostic to any run of light contributions.

The modified code compiled without complaint and there’s no discernible difference in either execution time or convergence diagnostics.

I’ve only run this on one set of galaxy data, from MaNGA plateifu 8565-3703 (mangaid 1-92735). This is one of Schawinski’s “blue early type galaxies,” chosen mostly because the data were binned to just 34 spectra so a complete set of model runs only took a few hours. And, as seen below the other things that don’t change are the modeled star formation histories.

I need to do some more validation exercises, but it appears my long ago conclusion that the choice of prior on the star formation history has little effect was correct. The data dominates the model outcome through the likelihood.

Update

Since I hit publish a few days ago I realized two things. First, the weights I applied to the stellar contributions to “encourage” a smoothly varying star formation rate were inverted. What I should have used was:

transformed data {
  vector[nt*nz] norm_sfr = (dT .* norm_st) / mean(dT .* norm_st);
}

The second thing I realized is this makes no difference at all given the form of the prior. The transformation simply maps a point on the simplex to another point on the simplex that has exactly the same probability density since by construction the prior is uniform on the simplex.

So, it should have been expected, and it’s a good thing that it in fact happened that the model runs produced the same results. What happens if I use the following prior, with the weights supplying the parameters for the Dirichlet?

target += dirichlet_lupdf(b_st_s | norm_sfr);

This failed to sample almost always. I’m not entirely sure why, but I suspect this turns a problem with relatively simple geometry at least with regards to the prior to one with a complex and troublesome geometry.

This little experiment actually told me nothing about the effect of priors on model star formation histories. Two of the priors are actually the same, and the third fails for reasons that aren’t completely clear. I may experiment with different forms of prior. I’m still, of course, looking for a new SSP model library.

MaNGA in M31 – final word for now and some pitfalls

When I was doing my initial fits to the M31 MaNGA spectra I noted two that I initially thought were contaminated by foreground stars, and therefore I masked them to prevent further analysis. One of the two, in MaNGA plateifu 9677-12701 (mangaid 52-8) is a certain foreground star and won’t be discussed further. The other, in plateifu 9678-12703 (mangaid 52-23), turns out to be a luminous red supergiant that’s a genuine resident of M31. This is confirmed by two nearly contemporaneous catalogs of M31 red supergiants: the one by Ren et al. (2021) that I noted previously and Massey et al. (2021), which I stumbled upon more recently.

In fact there are two comparably bright red supergiants in this IFU. One that’s about 9″ north of the masked one probably should have been masked by whatever criteria I used, but it’s likely I failed to notice the fit to the data since I don’t have the patience to look at every spectrum and data fit that pops up. So, here is the spectrum, displayed in (negative) magnitudes with arbitrary zero point. The blue spectrum is the closest match in my SSP library, a 10 Myr old population with the highest metallicity I used (2.5 Z. This is one of the theoretical spectra from PyPopstar). This sorta looks right except it’s much too blue. The solution to that is, of course, to add some reddening through dust attenuation.

plateifu 9678-12703 (M31 10 kpc ring) Spectrum contaminated with red supergiant and closest match SSP model spectrum

The maximum likelihood (non-negative weighted least squares) fit did just that, with only a single stellar contributor and a very high dust attenuation of τV = 3.0. This still doesn’t quite work: the residuals are rather strongly sloped in the blue and the details of the absorption features in the red aren’t quite right.

plateifu 9678-12703 – NNLS fit to spectrum contaminated with red supergiant

I still use a Calzetti attenuation relation in my NNLS fits. The Bayesian fits using Stan have the more flexible attenuation prescription that I described back in this post, and that helped considerably with the continuum as seen in the plot below. The absorption features in the red still aren’t fit well. The model has an even more extreme attenuation estimate with a much “grayer” than Calzetti slope, with τV = 4.38 ± 0.05 and δ = -0.33 ± 0.011see the link above for the meaning of these parameters.

9678-12703 – posterior predictive fit to spectrum

The model star formation history (displayed as a mass growth history below) isn’t completely implausible. The presence of a very luminous evolved star indicates the region is at least some Myr old, and a rapid onset and decline of star formation is typical for star forming regions in mature spirals. The recent episode of star formation added about 7% to the present day stellar mass, while at least 60% was in place by 8 Gyr ago (per the model).

plateifu 9678-12703 (M31 10 kpc ring) Model mass growth history for a region containing a bright red supergiant

Nevertheless I consider the model results to be highly suspect, mainly based on the very large optical depth. Massey estimates the attenuation for the red supergiant to be AV ≈ 1.18, although this is apparently based on a formula rather than a direct empirical estimate. But another piece of evidence that’s close to a smoking gun is an estimate based on the Balmer decrement of emission. Despite the lack of apparent ionizing sources outside the bright H II regions in the west there is widespread diffuse emission in this region with star-forming like line ratios. The estimated optical depth derived from the Balmer decrement for this region is τV, bd = 1.47 ± 0.23 (1σ), reasonably consistent with the Massey estimate and with the values derived for the rest of the IFU.

How widespread a problem is this? Below I plot the Balmer decrement derived optical depths against the stellar based estimates for all spectra in M31 MaNGA with star forming emission line ratios, about 11% of the entire sample. The 5 most extreme outliers are in this IFU in the regions surrounding the two bright red supergiants (the masked spectrum would also be in this region of the plot). The same 5 regions are also extreme outliers in the SFR vs. stellar mass and SFR vs. Hα plots that I showed early on. So, even though there are many cataloged supergiants in the study region these two appear to be uniquely bright and to have had the largest impact on model results.

M31 MaNGA – Optical depth off attenuation estimated from Balmer decrement vs. model values of τV

Of course there are hotter bright stars in the study area and these could affect results in different and possibly unexpected ways. For example the outermost IFU contains one bright star that GAIA estimates has a surface temperature of 5500 C, which would make it a G supergiant if it’s in M31. I noted in the last post that the model star formation history for that region looks like a post starburst with an age around 800 Myr. This is, I think, several galactic rotation periods, and stars born that long ago should have dispersed by now unless they’re gravitationally bound. There’s no sign of a star cluster there nor is there a cataloged one nearby, so it seems likely to me that the “starburst” is an artifact. As I noted in the last post though the fit to the data is quite good.

These examples illustrate an issue that’s fairly well known, which is that using simple stellar populations as building blocks of low mass stellar systems are potentially affected by so-called “stochastic” effects, which simply means that the distribution of stellar masses can vary randomly from what’s assumed in the SSP models. Specifically, in M31 there are individual stars luminous enough to affect spectra. One possible solution might be to add some stellar spectra to the library. I might give that a try some day.

I’m going away and won’t be writing for a while. I’m hoping to acquire or build a SSP model library based on SDSS MaStar spectra yet this year. This is a much larger collection of stellar spectra than has been previously available and it has the advantage of having the same flux calibration and (approximately) spectral resolution as the SDSS galaxy spectra. I also plan to return to my study of post-starburst galaxies.

The MaNGA M31 ancillary program – model details (part 1?)

After a fairly long break I want to get back to M31 and MaNGA for one, or perhaps several posts and take a more detailed look at my model results. I still haven’t decided where I’m going to take this investigation. I may examine every IFU or just the ones that I found most interesting, and I’m not sure which of the many quantities that I estimate I’ll discuss. Besides my models I’ve retrieved a number of catalogs of interesting objects using Aladin. These include in particular H II regions (Azimlu et al. 2011), OB associations (Magnier et al. 1993), and red supergiants (Ren et al. 2021). All of these are products of recent or ongoing star formation. There are of course a huge number of catalogs of just about every type of astronomical object found in galaxies, and I may examine some more depending on what interests me.

For orientation here’s a screencap of the Legacy Survey sky browser’s false color GALEX image of the northern half of M31 with the IFU positions overlaid and labelled with MaNGA’s plateifu identifiers. As a reminder these are all located within the PHAT survey footprint and specifically within the region for which star formation histories were estimated by Williams et al. (2017).

lg_manga
Screen capture of Legacy Survey Galex image of M31 with MaNGA IFU overlay

Before getting to individual IFU’s here is one more set of IFU-wide results. The following three graphs are model mass growth histories in units of present day solar mass per kiloparsec2. These are uncorrected for projection effects.

There are a couple interesting points here. There’s a clear stratification of mass density with projected radius, with about a factor 30 decline from the innermost to outermost IFU. This is in fairly good agreement with Williams’ estimate in their Figure 14.

The other thing to note is that all regions had most (> 55%) of their stellar mass in place by 8 Gyr ago and 92-99% in place by 1 Gyr ago. The largest fraction of recent star formation is in the IFU 9678-12703, which is very close to the region with the highest SFR in this half of the galaxy. There is also a trend towards later mass build up with increasing radius, which is completely consistent with the “inside-out” growth paradigm. The outermost IFU, 9678-12701 at about 16kpc radius has formed about 5% of its present day stellar mass in the past Gyr.

As I said in the previous post I don’t see clear evidence for a widespread burst of star formation that’s widely believed to have occurred around 2-4 Gyr ago. A confounding factor in my models is that they invariably show jumps in SFR at times when the interval between SSP model ages change and the two oldest of these occur at 1 and 4 Gyr, so this produces a possibly spurious period of apparently accelerated star formation. I hope to find (or perhaps produce) a set of SSP models with a better age distribution this year.

mghden_innerdisk
Growth of stellar mass density – inner disk M31 MaNGA IFU’s
Growth of stellar mass density – M31 MaNGA IFU’s in 10 kpc ring
Growth of stellar mass density – outer disk M31 MaNGA IFU’s

I think I’m going to hit publish now and resume with inner disk IFU’s next time.

The MaNGA M31 ancillary program – comparing star formation history models

In this post I’m going to compare IFU wide star formation histories from my models to those of Williams et al. (2017) in the nearest 83″ by 83″ PHAT tile to each MaNGA IFU in the study. I picked the Williams paper for comparison mostly because it’s possible to! They give a complete tabulation of model results for all regions and all 4 sets of isochrones that they used, and these are available through the Vizier service. Specifically I used their Table 2, which provides star formation rate densities summed over all metallicities. Since the SSP model spectra I use are based on BaSTI isochrones I initially compared to their BaSTI based models. One problem with the Williams comparison is the authors had a very wide youngest time bin of 300 Myr, which is where my models should generally have the highest precision (I make no strong claims about accuracy). It would be nice to do a similar comparison to the earlier companion paper on recent star formation by Lewis et al., which gives a much finer grained view of the last ~half Gyr, but unfortunately there is no published tabulation of their results.

At the other end of the timeline the oldest bin is also very wide, from 8 to 14 Gyr lookback time. This isn’t a surprise: the limits for reliable photometry of individual stars were rather shallow, no fainter than m = 28 or Mg ≈ 3.6 according to Lewis. This is brighter than the main sequence turnoff at 8 Gyr, so any information about the truly ancient star formation history is coming from giant branch stars which have very similar evolutionary tracks at old ages1Checked by downloading a few isochrones from the BaSTI website.

At the end of my last post I mentioned the necessity to correct densities for the rather large inclination of M31’s disk. It turns out though that I reproduce Williams’ Table 3 from their Table 2 if the densities are uncorrected. Their tabulated SFR densities are in units of 10-4 Marcmin-2/year. One arcminute at their adopted M31 distance is about 0.227 kpc, so to convert to star formation rates per kpc2 the values in table 2 are multiplied by 19.8 × 10-4. From my models I sum the star formation rates over all modeled spectra in each IFU and divide by the total area in fibers, with each fiber covering a projected area of 42.78 pc2. Note that I do not try to analyze a single composite spectrum summed over the entire IFU since dust attenuation is quite patchy.

The graphs below overlay my modeled star formation rate densities on those of Williams in the tile with the nearest center to that of each IFU. The ribbons indicate the nominal 95% credible limits of SFR. These are certainly wildly optimistic. Table 2 of the paper includes uncertainty estimates, which I chose not to include. SFR densities are linearly scaled with different limits for each plot. The time scale is logarithmic. Something in between linear and logarithmic would seem more appropriate since this perhaps gives too much space to very recent times, but I haven’t found a suitable scaling.

These are ordered into three groups by location: the inner disk which is everything inside the 10 kpc ring, the 10 kpc ring, and the outer disk which is everything outside the 10 kpc ring. There’s some ambiguity about the locations of the plateifu’s 9678-9101 and 9678-12704. The first of these is about 12 kpc from the nucleus in what could be either a wide section of the 10 kpc ring or a separate structure. 9678-12704 appears in projection to be just outside the 10 kpc ring but it may be considerably farther out in a segment of spiral arm at ~15 kpc.

Commentary will continue after the graphs. I will discuss the individual IFU’s in more detail in later post(s).

Inner Disk

sfh_innerdisk
M31 MaNGA ancillary program – my star formation history models summed over each IFU compared to nearest tile in Williams et al. (2017) with BaSTI isochrones. Inner disk IFU’s.

10 Kpc ring

sfh_10kpc
M31 MaNGA ancillary program – my star formation history models summed over each IFU compared to nearest tile in Williams et al. (2017) with BaSTI isochrones. 10 kpc ring IFU’s.

Outer disk

sfh_outerdisk
M31 MaNGA ancillary program – my star formation history models summed over each IFU compared to nearest tile in Williams et al. (2017) with BaSTI isochrones. Outer disk IFU’s.

The first two things I noticed were that star formation in every region declines monotonically from very early times to at least 4 Gyr ago. It also starts out lower than in the PHAT team’s models. Because of this early time mass deficit all of my models have smaller current day stellar mass densities by varying amounts. I don’t really have a pat explanation for this. Some authors have posited a “dazzle effect”2I’m going to discuss this a little further at the end of the post where recent star formation obscures the contribution of old populations. It’s certainly likely that this occurs, but if these Bayesian models are behaving as I hoped this lack of information should manifest as larger uncertainties rather than a systematic bias. Well, my hope could be wrong. On the other hand I don’t see strong evidence in these models for such an effect. From my eyeball analysis I don’t see an obvious correlation between present day star formation and the size of the early time deficit.

Another possibility is a systematic difference in the amount or shape of attenuation between my models and theirs. There is another well known “degeneracy” between stellar age and attenuation in SFH modeling, but I haven’t yet investigated whether this could be occurring here.

The PHAT models have a very long interval from 8 to ~2-3 Gyr ago with very little star formation. Some authors find evidence for a large increase from about 2-4 Gyr which is usually attributed to a merger or perhaps close encounter with M33. This isn’t seen in the BaSTI based models but there is a large more recent burst from about 1-2 Gyr lookback time. My models see neither a cessation of star formation nor a particularly large burst at intermediate ages. As I’ve noted before my models “want” to have smoothly time varying light3and therefore mass contributions and this might make a modest burst at moderately large ages difficult to discern. Another confounding factor arises from the abrupt changes in age intervals (at 0.1, 0.5, 1, and 4 Gyr) which results in the sawtooth pattern in SFR that’s obvious in every plot above.

At ages younger than 1 Gyr there’s generally good agreement about the course of star formation up until the youngest age bin of width 300Myr in the PHAT models. My models have anywhere from slightly to dramatically higher SFR densities averaged over the most recent age bin. I suspect this is because many of the IFU positions were chosen to be in regions with active star formation. In particular the plateifu 9678-12703 (mangaid 52-23) is very close to the region in the 10 kpc ring with the highest density of ongoing star formation in the northern half of the disk.

I plan to discuss the individual IFUs in more detail in a later post. Below the fold are some more graphics: mass growth histories and SFR densities compared to the PADOVA isochrone based models.

Continue reading “The MaNGA M31 ancillary program – comparing star formation history models”

The MaNGA M31 ancillary program – preliminaries

One of the ancillary programs (with principal investigator Julianne Dalcanton) in the final MaNGA release targeted 18 fields in the disk of the Andromeda galaxy M31. The targets were selected from within the footprint of the “Panchromatic Hubble Andromeda Treasury,” aka PHAT1not my coinage., also with PI Dalcanton. The initial PHAT survey description was in Dalcanton et al. (2012) and was followed by a lengthy series of papers. Especially relevant for this discussion are two papers describing estimates of the recent and ancient star formation histories of the disk outside the area dominated by bulge light: Lewis et al. (2015), “The Panchromatic Hubble Andromeda Treasury. XI. The Spatially Resolved Recent Star Formation History of M31” and Williams et al. (2017), “PHAT. XIX. The Ancient Star Formation History of the M31 Disk.” For reference here is a mosaic of HST images in the F475W filter with the IFU locations overlaid:

m31manga_phathst
Mosaic of HST F475W images of PHAT study region with M31 MaNGA IFU positions overlaid

Zooming out to show the whole disk here they are overlaid on a false color FUV+NUV image from GALEX, which gives a pretty good picture of where stars are actually forming:

m31manga_galex
GALEX false color NUV+FUV image of M31 with MaNGA IFU positions overlaid – screencap from Aladin

This data set provides an excellent opportunity to compare my SFH modeling code to a completely different, more direct, method of inferring star formation histories namely counting resolved stars in color magnitude diagrams. I recently completed model runs for all 18 IFU’s with the same Voronoi binning of stacked RSS spectra, the same modeling code and SSP model spectra as I’ve used for a while now.

There’s no redshift listed in the DRP catalog; NED gives a heliocentric redshift of -0.001, but for purposes of calculating intrinsic quantities I need the “Hubble flow” redshift. I adopted a distance of 761 (± 11) kpc or distance modulus of (m-M)0 = 24.407 from Li et al. (2021), which is the most recent and according to the authors most precise determination to date. With my adopted Hubble constant of H0 = 70 km/sec/Mpc this makes the Hubble flow recession velocity 53.27 km/sec or zdist = 0.0001777. The angular scale is 3.69 pc/arc-second. This distance estimate is a few percent smaller than the PHAT team authors and most other recent literature I reviewed, but fortunately most other sources of uncertainty are much larger.

An issue I noticed early on was the modelled values of the optical depth of attenuation were right at 0 for almost all spectra with only a few much larger exceptions. A quick check of the metadata showed that the values adopted for the foreground galactic extinction almost certainly were taken from the SFD dust maps which faithfully capture the intrinsic dust content of M31 albeit at rather low resolution. These hugely overestimate the actual foreground galactic extinction and that has multiple undesirable consequences. So, I assigned a single extinction value of E(B-V) = 0.055 to all IFU’s, consistent with the NED value of AV = 0.17 mag. The preliminary runs were redone with the newly adopted extinction value.

After binning to a minimum mean SNR of 5 there were 2,624 spectra in the 18 IFUs, of which I ran models for 2,621. Three spectra had apparent foreground stars, although one of those might actually be a red supergiant in M31. The fibers are basically sampling star cluster size and stellar mass regions so a single extremely luminous star could potentially affect a spectrum.

I’m only going to show a few summary results for the entire sample in this post. My goal is to do a more detailed quantitative comparison to (at least) the SFH models of Wilson, for which there are extensive results tabulated. There are of course many catalogs of interesting objects within M31, and I plan to look at some of them.

First, here is a plot of the (100 Myr averaged) star formation rate density against stellar mass density, color coded by BPT diagnostic. The solid line is my estimate of the “spatially resolved star forming main sequence” based on a small sample of non-barred spiral galaxies. The dashed line is the estimate of Bluck et al. (2020), which I commented previously appears to mark approximately the location of the green valley at least with regard to my models. A striking feature of this plot is the apparent stratification into at least three distinct groups that can be interpreted as starforming, quiescently evolving, and passively evolving. I suspect this observed stratification is just the result of hand picking a small number of “interesting” regions. Most or perhaps all of the points in the passively evolving group are in the IFU closest to the bulge, while most of those along and above the SFMS lie near the most vigorously star forming regions in the PHAT footprint. Especially noteworthy are 5 outliers that are well above any others in the plot in terms of SFR density. These are all in the same IFU (plateifu 9678-12703) which is located within the largest star forming region in that quadrant of the “10 kpc ring.”

m31manga_sigmam_sigmasfr
100 Myr average star formation rate density vs. stellar mass density for 2621 binned spectra in M31 disk. Solid and dashed lines are my and Bluck’s central estimates of the “spatially resolved star forming main sequence.

Next are plots of star formation rate density against Hα luminosity density. The left panel is for all spectra color coded by BPT diagnostic, with Hα adjusted by the modeled amount of stellar attenuation. The right panel shows regions with star forming BPT diagnostics only, with Hα corrected by the observed Balmer decrement. The solid line in both panels is Calzetti’s calibration of the Hα – SFR relationship. The relationships plotted here are consistent with what I’ve seen in other MaNGA samples and with published values, which is encouraging.

m31manga_sigmaha_sigmasfr
Star formation rate density vs. Hα luminosity density for 2621 binned spectra in M31 disk. (L) Emission corrected for modeled stellar attenuation. (R) For regions with star forming emission line ratios only: emission corrected from estimated Balmer decrement.

The obvious point of comparison to my models are the detailed star formation histories in the two PHAT papers mentioned at the top. Unfortunately there is no detailed tabulation of model results in the paper by Lewis et al. The paper by Williams et al. has extensive tables, but there are still a few obstacles to detailed comparisons which I will discuss next time.

A few more items from my handwritten notes that I want to get in pixels. I have never previously tried to correct surface densities for inclination in disk galaxies, but for comparison purposes and because of the large inclination of M31’s disk I need to do so here. I adopted an inclination angle of 77°, so a 1″ radius fiber covers a 3.69 x 16.4 pc (semi major and minor axes) elliptical region, or 190 pc2. Densities need to be adjusted downward by a factor 4.45 or -0.648 dex2This adjustment was not made in the plots above. Since these are plots of densities against densities all points would just shift downwards along lines of slope one..

In order to achieve 100% coverage of the IFU footprints the exposures were dithered to three different positions with overlapping fiber positions. Comparing the area in fibers to the area in spaxels in the cubes the overfilling factor averages 0.217 dex or 65%. The total area in all cubes is 10,731 arcsec2, or a deprojected area of 0.65 kpc2. The most distant IFU from the nucleus is at a projected radius of about 16 kpc. A simple extrapolation to the ≈800 kpc2 area of the disk within that radius is probably unsafe.

One final map to anticipate the next post(s). Wilson provides tables of model star formation rates for 16 age bins, 826 regions, and 4 different sources of isochrones including the same BASTI isochrones I use. The complete data set is available through Vizier. In the plot below I created a map of the recent star formation rate density interpolated to nominal 10″ resolution from their Table 2 models with BASTI isochrones. This should be compared to their Figure 16 (they use logarithmic scaling).

m31manga_wilsonsfr
Current (300 Myr average) star formation rate density in the PHAT footprint per models of Wilson et al. (2017) with positions of MaNGA IFUs overlaid.

Journal notes: Haines et al. (2015), “Testing the modern merger hypothesis…”

While browsing through the ADS listing of papers that cite Schawinski’s paper that I’ve been discussing for a while I came across this one by Haines et al. with the full title “Testing the modern merger hypothesis via the assembly of massive blue elliptical galaxies in the local Universe”. Besides being on the same theme of searching for post-starburst or “transitional” galaxies in the local universe that I’ve been pursuing for some time the paper was interesting because it made use of IFU based spectroscopic data that predates MaNGA. As it happens 4 of the 12 galaxies have observations in the final MaNGA release, providing an excellent opportunity to compare results from completely independent data sets.

The “modern merger hypothesis” that the authors tested relates to a topic I’ve discussed before, which is that N-body simulations show that strong, centrally concentrated starbursts are a possible outcome of major gas rich galaxy mergers around the time of coalescence. If some feedback process (an AGN or supernovae) rapidly quenches star formation there will ensue a period of time when the galaxy will be recognizable as post-starburst.

In a series of long and rather difficult (and influential judging by the number of citations) Hopkins and collaborators (2006, 2008a, 2008b) have made a case that major gas rich mergers with accompanying starbursts are in fact the major pathway to the formation of modern elliptical galaxies. They claim that their merger hypothesis accounts for a variety of phenomena, including the growth and evolution of supermassive black holes and quasars.

The specific aspect of the merger hypothesis this study tried to address was the prevalence of strong centrally concentrated starbursts in a sample of ellipticals in the process of forming as evidenced by visible disturbances consistent with recent mergers. The main tool they used was a suite of simple star formation history models with exponentially decaying star formation rate with single (also exponentially decaying) bursts on top of varying ages and decay time scales. They used these to predict just two quantities: Balmer absorption line strength measured by the average of the Lick HδA and HγA indexes, and the 4000Å break strength index Dn4000. For reference here is a screen grab of their model trajectories:

Predected trajectories in the Hδ – Dn4000 plane per Haines et al. (2015). Clipped from the electronic journal paper.

This is a pretty standard calculation variations of which have been performed for decades, and this graph looks much like others I have seen in the literature. A fairly basic problem with it though is that position in the Balmer – D4000 plane doesn’t uniquely constrain even the recent stellar evolution. In astronomers’ parlance there is a “degeneracy”1the term refers to a situation in which multiple combinations of some parameters of interest produce effectively equivalent values of some observable(s), or of course the converse. The best known example is the “age-metallicity degeneracy,” which refers to the fact that an old metal poor population looks like a younger metal rich one in several respects such as broad band colors. between burst strength (if any) and burst age. This is a well known problem with the Balmer line strength index that was already recognized by Worthey and Ottaviani (1997), who developed these indexes. Adding a second index in the form of the 4000Å break strength doesn’t break the degeneracy: there are regions of the plane where bursting and non-bursting populations overlap, as can be seen clearly in the graphic above. This is actually a problem for any attempt to identify post-starburst galaxies. After correcting for emission most ordinary starforming galaxies have strong Balmer absorption lines, so using that index alone will certainly produce many false positives. On the other hand selection criteria like those used by Goto and many others before and after — selecting for both strong Balmer absorption and weak emission — will capture only a small interval in post-starburst galaxies’ life cycles.

hd_d4000_bigsample
Hδ line strength vs. 4000Å break index for a large (~380K) sample of SDSS galaxy spectra. Measurements from the MPA-JHU analysis pipeline downloaded from SDSS Skyserver

Let’s get to results. Some basic details of the sample are in the table below. Morphological classifications are from McIntosh et al. (2014) as given in this paper. The abbreviations are SPM: spherical post merger; pE: peculiar Elliptical. The two marked pE/SPM didn’t have a strong consensus among several professional classifiers. I list them in order of my own visual impression of degree of disturbance. I also list redshifts taken from the MaNGA catalog and Petrosian colors.

NED nameNYU IDmangaidplateifuMorphzu-rg-i
NGC 39215410441-61744510510-6103SPM0.0191.970.86
MRK 3857194861-6049708940-6102pE/SPM0.0281.430.63
MRK 3661009171-6033097993-1902pE/SPM0.0271.590.79
NGC 1149223181-371558154-6103pE0.0292.291.11
Columns: (1) Common catalog designation (NED name). (2) NYU VAC ID. (3) MaNGA mangaid. (4) MaNGA plateifu. (5) Morphology (see text). (6) redshift from MaNGA DRP catalog. (7-8) Petrosian u-r and g-i colors from NYU VAC via the MaNGA DRP catalog.

The main prediction of the merger with accompanying centrally concentrated starburst hypothesis the paper tests is that the Balmer absorption index should be large and have a negative gradient with radius while the 4000Å break strength should be low with a positive gradient. The authors concluded that only one member of their sample — nyu541044 — clearly falls in the post-starburst region (marked as region 4 in the graph above) of the <Hδ, Hγ> – Dn4000 plane. The two pE/PM galaxies, both of which are in my sample, lie in the starforming region 1. They inferred from this that these galaxies are undergoing at most a weak burst. I’m going to mildly disagree with that conclusion.

Screenshot from 2022-07-07 15-23-36
Measured values for the specified indexes from Haines et al. (2015). Clipped from the electronic journal paper.

I have calculated the pseudo Lick index HδA and Dn4000 as part of my analysis “pipeline” since I started this hobby. I actually make these measurements in the initial maximum likelihood fitting step since they don’t depend on modeling except for small (usually) emission corrections. I don’t calculate an Hγ index, but its theoretical behavior is similar to Hδ. I’m trying here just to verify the approximate magnitude and radial trends of the chosen indexes. The two IFUs used in the Haines study had larger spatial coverage than these MaNGA observations (but much smaller wavelength coverage, which will become important). Instead of their strategy of binning in annuli I used my usual Voronoi binning strategy with a minimum target S/N. There were some oddities in the NYU estimates of effective radii so I chose to use distances from the IFU center in kpc for these plots. The distances assigned to the multiply binned spectra are the same as Cappelari’s published code produces; for single fiber spectra it’s just the position of the fiber center.

My measurements agree reasonably well with those of Haines et al. All three of the most disturbed galaxies have central Hδ indexes > 5Å with NGC 3921 (plateifu 10510-6103, nyu541044) having a larger central value and steeper gradient in the inner few kpc than the two pE/SPM galaxies. The fourth galaxy shows no obvious trend in either index with radius2The next several plots show trend lines for each galaxy computed by fitting simple loess curves to the data using the default parameters in ggplot2. These, and especially the confidence bands included in the plots, should not be taken seriously!. The central values where the S/N is highest are in good agreement.

Lets turn to the results of star formation history models, which I ran on all 4 data sets. First, here are 100Myr averaged star formation rate density and specific star formation rate versus distance:

Star formation rate density vs. distance from IFU center (kpc) for 4 disturbed early type galaxies.
Specific star formation rate density vs. distance from IFU center (kpc) for 4 disturbed early type galaxies.

Three of these galaxies are clearly experiencing centrally concentrated episodes of star formation, and two are at or near starburst levels in specific star formation rate near their centers. As seen below two of these straddle my estimate of the “spatially resolved star forming main sequence” while the one presumed post-starburst galaxy reaches it in the central region.

mstar_sfr_4spm
Star formation rate density versus stellar mass density for 4 disturbed early type galaxies

As I’ve shown several times before there’s a reasonably tight linear relationship between modeled star formation rate and Hα luminosity density. The plot shows Hα luminosity density corrected for modeled stellar redenning, which certainly underestimates attenuation in emission regions. The modeled star formation rates are consistently above the Kennicut relation shown as the straight line as I’ve seen in every sample I’ve looked at.

Star formation rate density vs. Hα luminosity density for 4 disturbed early type galaxies

Finally, lets take a look at detailed star formation histories. Instead of my usual practice of plotting them all in a grid here I just display 2 binned star formation histories. One comprises the innermost 7 bins, which since the fibers are arranged in a hexagonal grid should form a regular hexagon around the IFU center. These range in “radius” from about 0.75 to 1.1 kpc in these four galaxies. The second is for an “annulus” in approximately the outer kpc of each IFU. The extent of the IFU footprints ranges from 3.1 to 5.9 kpc. I calculate these by summing the contributions in each SFH model contributing to the bins, not by running new models for binned spectra. Since the dithered fiber positions overlaps this overestimates the total mass in each bin, but I care about the shape and timing of events rather than the absolute values of star formation rate estimates.

The next 4 plots display the results. Lookback time is logarithmically scaled with the same range and ticks for each SFH. Vertical scales are linear and differ for each graph. The graphs are in the same order as the basic information table above. As I’ve written before these models “want” to have smoothly varying mass per time bin which has the unfortunate effect of producing jumps in the apparent SFR when the bin widths change. In the BaSTI isochrone based SSP models these occur at 100 Myr, 1 Gyr, and 4 Gyr and can sometimes be quite prominent.

With caveats out of the way the one clear post-starburst in the sample had (per the model) a powerful and short starburst at ≈300 Myr lookback time, with a small amount continuing to the present (this can’t be seen at the scale of the graph, but ongoing star formation is ~1 M/yr). The total mass contribution from the burst and subsequent star formation is around 15%.

The two apparent ongoing starbursts have later bursts of star formation that are slightly weaker in terms of total mass contribution and peak star formation rate, but still quite significant. All three of the starburst/post-starburst galaxies appear to have had two major waves of late time (last ~2 Gyr or less) star formation. As I’ve written before in merger simulations the progenitors usually complete a few orbits before coalescence, with some enhanced star formation around each perigalactic passage. I hesitate to take these models that literally.

Turning finally to the last and least disturbed galaxy, NGC 1149, despite the bursty appearance of the SFH there’s no evidence for a major starburst in the cosmologically recent past. Whether an older starburst can be detected in this kind of modeling approach needs investigating.

One last set of graphs that may be useful. These show cumulative star formation histories — basically the cumulative sum of mass contributions starting from the oldest time bin. This is similar to a mass growth history which is a popular visualization. In my calculation of the latter the contributions are to the present day stellar mass, so an allowance for mass loss and remnant mass is made3these come from the source of the SSP models and are themselves models. Probably they are somewhat better than guesses. These things are basically black boxes to users.. The graphs are for the central regions only. Note the major virtue of these is that the contributions of major episodes of star formation can be estimated at a glance.

Cumulative star formation histories for central regions of 4 disturbed early type galaxies

To wrap up this part of the post 3 of these galaxies are compatible with the “modern merger hypothesis,” that is they have experienced centrally concentrated but spatially wide spread starbursts. The reason two of them don’t have post-starburst characteristics in the Hδ – D4000 plane is their starbursts are still underway. The current burst of star formation contributes about 5-10% of the mass in the central regions of these two. How much more is available is unknown (at least to me until I get around to finding out if there are HI mass estimates available).

Future plans: I’ve completed model runs on the 24 “post-starburst” galaxies in the MaNGA ancillary program dedicated to them. I may have something to say about them. I also may have something to say about one of the Zoogems targets that I had a small part in selecting.

Continue reading “Journal notes: Haines et al. (2015), “Testing the modern merger hypothesis…””

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 offsets1these formulas are approximate but close enough for the present purpose.

\( \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:

  1. Investigate informative priors.
  2. 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.
  3. Investigate systematics between my estimates and those in the DAP.

Continue reading “Adding emission lines to non-parametric kinematic models”

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.

thumbnails
SDSS thumbnail images of Zoogems targets with MaNGA data

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 like1this is actually from the data cube:

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 though2interesting enough that I made a couple posts on the Galaxy Zoo talk forum, which I rarely do anymore. 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 palette3I know data visualization experts frown on rainbows, but I think they’re OK for things like velocities or redshifts:

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.

NGC 810 – sample nonparametric losvd

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.

Continue reading “NGC 810 – interesting kinematics in a Zoogems and MaNGA target”

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.