Markarian 848 – Closing topics

I’m going to close out my analysis of Mrk 848 for now with three topics. First, dust. Like most SED fitting codes mine produces an estimate of the internal attenuation, which I parameterize with τV, the optical depth at V assuming a conventional Calzetti attenuation curve. Before getting into a discussion for context here is a map of the posterior mean estimate for the higher S/N target binning of the data. For reference isophotes of the synthesized r band surface brightness taken from the MaNGA data cube are superimposed:

mrk848_tauv_map
Map of posterior mean of τV from single component dust model fits with Calzetti attenuation

This compares reasonably well with my visual impression of the dust distribution. Both nuclei have very large dust optical depths with a gradual decline outside, while the northern tidal tail has relatively little attenuation.

The paper by Yuan et al. that I looked at last time devoted most of its space to different ways of modeling dust attenuation, ultimately concluding that a two component dust model of the sort advocated by Charlot and Fall (2000) was needed to bring results of full spectral fitting using ppxf on the same MaNGA data as I’ve examined into reasonable agreement with broad band UV-IR SED fits.

There’s certainly some evidence in support of this. Here is a plot I’ve shown for some other systems of the estimated optical depth of the Balmer emission line emitting regions based on the observed vs. theoretical Balmer decrement (I’ve assumed an intrinsic Hα/Hβ ratio of 2.86 and a Calzetti attenuation relation) plotted against the optical depth estimated from the SFH models, which roughly estimates the amount of reddening needed to fit the SSP model spectra to the observed continuum. In some respects this is a more favorable system than some I’ve looked at because Hβ emission is at measurable levels throughout. On the other hand there is clear evidence that multiple ionization mechanisms are at work, so the assumption of a single canonical value of Hα/Hβ is likely too simple. This might be a partial cause of the scatter in the modeled relationship, but it’s encouraging that there is a strong positive correlation (for those who care, the correlation coefficient between the mean values is 0.8).

The solid line in the graph below is 1:1. The semi-transparent cloud of lines are the sampled relationships from a Bayesian errors in variables regression model. The mean (and marginal 1σ uncertainty) is \(\tau_{V, bd} = (0.94\pm 0.11) + (1.21 \pm 0.12) \tau_V\). So the estimated relationship is just a little steeper than 1:1 but with an offset of about 1, which is a little different from the Charlot & Fall model and from what Yuan et al. found, where the youngest stellar component has about 2-3 times the dust attenuation as the older stellar population. I’ve seen a similar not so steep relationship in every system I’ve looked at and don’t know why it differs from what is typically assumed. I may look into it some day.

τV estimated from Balmer decrement vs. τV from model fits. Straight line is 1:1 relation. Cloud of lines are from Bayesian errors in variables regression model.

I did have time to run some 2 dust component SFH models. This is a very simple extension of the single component models: a single optical depth is applied to all SSP spectra. A second component with the optical depth fixed at 1 larger than the bulk value is applied only to the youngest model spectra, which recall were taken from unevolved SSPs from the updated BC03 library. I’m just going to show the most basic results from the models for now in the form of maps of the SFR density and specific star formation rate. Compared to the same maps displayed at the end of the last post there is very little difference in spatial variation of these quantities. The main effect of adding more reddened young populations to the model is to replace some of the older light — this is the well known dust-age degeneracy. The average effect was to increase the stellar mass density (by ≈ 0.05 dex overall) while slightly decreasing the 100Myr average SFR (by ≈ 0.04 dex), leading to an average decrease in specific star formation rate of ≈ 0.09 dex. While there are some spatial variations in all of these quantities no qualitative conclusion changes very much.

mrk848_sigma_sfr_sfr_2dust_maps
Results from fits with 2 component dust models. (L) SFR density. (R) Specific SFR

Contrary to Yuan+ I don’t find a clear need for a 2 component dust model. Without trying to replicate their results I can’t say why exactly we disagree, but I think they erred in aggregating the MaNGA data to the lowest spatial resolution of the broad band photometric data they used, which was 5″ radius. There are clear variations in physical conditions on much smaller scales than this.

Second topic: the most widely accepted SFR indicator in visual wavelengths is Hα luminosity. Here is another plot I’ve displayed previously: a simple scatterplot of Hα luminosity density against the 100Myr averaged star formation rate density from the SFH models. Luminosity density is corrected for attenuation estimated from the Balmer decrement and for comparison the light gray points are the uncorrected values. Points are color coded by BPT class determined in the usual manner. The straight line is the Hα – SFR calibration of Moustakas et al. (2006), which in turn is taken from earlier work by Kennicutt.

Model SFR density vs. Hα luminosity density corrected for extinction estimated from Balmer decrement. Light colored points are uncorrected for extinction. Straight line is Hα-SFR calibration from Moustakas et al. (2006)

Keeping in mind that Hα emission tracks star formation on timescales of ∼10 Myr1to the extent that ionization is caused by hot young stars. There are evidently multiple ionizing sources in this system, but disentangling their effects seems hard. Note there’s no clear stratification by BPT class in this plot. this graph strongly supports the scenario I outlined in the last post. At the highest Hα luminosities the SFR-Hα trend nicely straddles the Kennicutt-Moustakas calibration, consistent with the finding that the central regions of the two galaxies have had ∼constant or slowly rising star formation rates in the recent past. At lower Hα luminosities the 100Myr average trends consistently above the calibration line, implying a recent fading of star formation.

The maps below add some detail, and here the perceptual uniformity of the viridis color palette really helps. If star formation exactly tracked Hα luminosity these two maps would look the same. Instead the northern tidal tail in particular and the small piece of the southern one within the IFU footprint are underluminous in Hα, again implying a recent fading of star formation in the peripheries.

(L) Hα luminosity density, corrected for extinction estimated by Balmer decrement. (R) SFR density (100 Myr average).

Final topic: the fit to the data, and in particular the emission lines. As I’ve mentioned previously I fit the stellar contribution and emission lines simultaneously, generally assuming separate single component gaussian velocity dispersions and a common system velocity offset. This works well for most galaxies, but for active galaxies or systems like this one with complex velocity profiles maybe not so much. In particular the northern nuclear region is known to have high velocity outflows in both ionized and neutral gas due presumably to supernova driven winds. I’m just going to look at the central fiber spectrum for now. I haven’t examined the fits in detail, but in general they get better outside the immediate region of the center. First, here is the fit to the data using my standard model. In the top panel the gray line, which mostly can’t be seen, is the observed spectrum. Blue are quantiles of the posterior mean fit — this is actually a ribbon, although its width is too thin to be discernable. The bottom panel are the residuals in standard deviations. Yes, they run as large as ±50σ, with conspicuous problems around all emission lines. There are also a number of usually weak emission lines that I don’t track that are present in this spectrum.

mrk848_fit_central_spec
Fit to central fiber spectrum; model with single gaussian velocity distributions.

I have a solution for cases like this which I call partially parametric. I assume the usual Gauss-Hermite form for the emission lines (as in, for example, ppxf) while the stellar velocity distribution is modeled with a convolution kernel2I think I’ve discussed this previously but I’m too lazy to check right now. If I haven’t I’ll post about it someday. Unfortunately the Stan implementation of this model takes at least an order of magnitude longer to execute than my standard one, which makes its wholesale use prohibitively expensive. It does materially improve the fit to this spectrum although there are still problems with the stronger emission lines. Let’s zoom in on a few crucial regions of the spectrum:

Zoomed in fit to central fiber spectrum using “partially parametric velocity distribution” model. Grey: observed flux. Blue: model.

The two things that are evident here are the clear sign of outflow in the forbidden emission lines, particularly [O III] and [N II], while the Balmer lines are relatively more symmetrical as are the [S II] doublet at 6717, 6730Å. The existence of rather significant residuals is likely because emission is coming from at least two physically distinct regions while the fit to the data is mostly driven by Hα, which as usual is the strongest emission line. The fit captures the emission line cores in the high order Balmer lines rather well and also the absorption lines on the blue side of the 4000Å break except for the region around the [Ne III] line at 3869Å.

I’m mostly interested in star formation histories, and it’s important to see what differences are present. Here is a comparison of three models: my standard one, the two dust component model, and the partially parametric velocity dispersion model:

mrk848_centralsfr3ways
Detailed star formation history models for the northern nucleus using 3 different models.

In fact the differences are small and not clearly outside the range of normal MCMC variability. The two dust model slightly increases the contribution of the youngest stellar component at the expense of slightly older contributors. All three have the presumably artifactual uptick in SFR at 4Gyr and very similar estimated histories for ages >1 Gyr.

I still envision a number of future model refinements. The current version of the official data analysis pipeline tracks several more emission lines than I do at present and has updated wavelengths that may be more accurate than the ones from the SDSS spectro pipeline. It might be useful to allow at least two distinct emission line velocity distributions, with for example one for recombination lines and one for forbidden. Unfortunately the computational expense of this sort of generalization at present is prohibitive.

I’m not impressed with the two dust model that I tried, but there may still be useful refinements to the attenuation model to be made. A more flexible form of the Calzetti relation might be useful for example3there is recent relevant literature on this topic that I’m again too lazy to look up.

My initial impression of this system was that it was a clear false positive that was selected mostly because of a spurious BPT classification. On further reflection with MaNGA data available it’s not so clear. A slight surprise is the strong Balmer absorption virtually throughout the system with evidence for a recent shut down of star formation in the tidal tails. A popular scenario for the formation of K+A galaxies through major mergers is that they experience a centrally concentrated starburst after coalescence which, once the dust clears and assuming that some feedback mechanism shuts off star formation leads to a period of up to a Gyr or so with a classic K+A signature4I previously cited Bekki et al. 2005, who examine this scenario in considerable detail.Capturing a merger in the instant before final coalescence provides important clues about this process.

To the best of my knowledge there have been no attempts at dynamical modeling of this particular system. There is now reasonably good kinematic information for the region covered by the MaNGA IFU, and there is good photometric data from both HST and several imaging surveys. Together these make detailed dynamical modeling technically feasible. It would be interesting if star formation histories could further constrain such models. Being aware of the multiple “degeneracies” between stellar age and other physical properties I’m not highly confident, but it seems provocative that we can apparently identify distinct stages in the evolutionary history of this system.

Markarian 848 – detailed star formation histories

In this post I’m going to take a quick look at detailed, spatially resolved star formation histories for this galaxy pair and briefly compare to some of the recent literature. Before I start though, here is a reprocessed version of the blog’s cover picture with a different crop and some Photoshop curve and level adjustments to brighten the tidal tails a bit. Also I cleaned up some more of the cosmic ray hits.

mrk_848_hst_crop_square
Markarian 848 – full resolution crop with level curve adjustment HST ACS/WFC F435W/F814W/F160 false color image

The SFH models I’m going to discuss were based on MaNGA data binned to a more conservative signal to noise target than in the last post. I set a target S/N of 25 — Cappellari’s Voronoi binning algorithm has a hard coded acceptance threshold of 80% of the target S/N, which results in all but a few bins having an actual S/N of at least 20. This produced a binned data set with 63 bins. The map below plots the modeled (posterior mean log) stellar mass density, with clear local peaks in the bins covering the positions of the two nuclei. The bins are numbered in order of increasing distance of the weighted centroids from the IFU center, which corresponds with the position of the northern nucleus. The center of bin 1 is slightly offset from the northern nucleus by about 2/3″. For reference the angular scale at the system redshift (z ≈ 0.040) is 0.8 kpc/” and the area covered by each fiber is 2 kpc2.

Although there’s nothing in the Voronoi binning algorithm that guarantees this it did a nice job of partitioning the data into physically distinct regions, with the two nuclei, the area around the northern nucleus, and the bridge between them all binned to single fiber spectra. The tidal tails are sliced into several pieces while the very low surface brightness regions to the NE and SW get a few bins.

mrk848_stmass_map
Stellar mass density and bin numbers ordered in increasing distance from northern nucleus

I modeled star formation histories with my largest subset of the EMILES library consisting of SSP model spectra for 54 time bins and 4 metallicity bins. As in previous exercises I fit emission lines and the stellar contribution simultaneously, which is riskier than usual in this system because some regions, especially near the northern nucleus, have complex velocity structures producing emission line profiles that are far from single component gaussians. I’ll look at this in more detail in a future post, but for now let’s just take these SFH models at face value and see where they lead. As I’ve done in several previous posts what’s plotted below are star formation rates in solar masses/year over cosmic time, counting backwards from “now” (now being when the light left the galaxies about 500Myr ago). This time I’ve used log scaling for both axes and for the first time in these posts I’ve displayed results for every bin ordered by increasing distance from the IFU center. The first two rows cover the central ~few kpc region surrounding the northern nucleus. The southern nucleus covered by bin 44 is in row 8, second from left. Both the time and SFR scales are the same for all plots. The black lines are median posterior marginal estimates and the ribbons are 95% posterior confidence intervals.

mrk848_sfh_bin20
Star formation histories for the Voronoi binned regions shown in the previous map. Numbers at the top correspond to the bin numbers in the map.

Browsing through these, all regions show ∼constant or gradually decreasing star formation up to about 1 Gyr ago1You may recall I’ve previously noted there is always an uptick in star formation rate at 4 Gyr in my EMILES based models, and that’s seen in all of these as well. This must be spurious, but I still don’t know the cause.. This of course is completely normal for spiral galaxies.

Most regions covered by the IFU began to show accelerated and in some areas bursts of star formation at ∼1 Gyr. In more than half of the bins the maximum star formation rate occurred around 40-60 Myr ago, with a decline or in some areas cessation of star formation more recently. In the central few kpc around the northern nucleus on the other hand star formation accelerated rapidly beginning ~100 Myr ago and continues at high levels. The southern nucleus has a rather different estimated recent star formation history, with no (visible) starburst and instead gradually increasing SFR to a recent maximum. Ongoing star formation at measurable levels is much more localized in the southern central regions and weaker by factors of several than the northern companion.

Here’s a map that I thought would be too noisy to be informative, but turns out to be rather interesting. This shows the spatial distribution of the lookback time to the epoch of maximum star formation rate estimated by the marginal posterior mean. The units displayed are log(age) in Myr; recall that I added unevolved SSP models from the 2013 update of the BC03 models to the BaSTI based SSP library, assigning them an age of 10Myr, so a value of 1 here basically means ≈ now.

mrk848_lbt_to_maxsfr
Look back time to epoch of maximum star formation rate as estimated by marginal posterior mean. Units are log(age in Myr).

To summarize, there have been three phases in the star formation history of this system: a long period of normal disk galaxy evolution; next beginning about 1 Gyr ago widespread acceleration of the SFR with localized bursts; and now, within the past 50-100 Myr a rapid increase in star formation that’s centrally concentrated in the two nuclei (but mostly the northern one) while the peripheral regions have had suppressed activity.

I previously reviewed some of the recent literature on numerical studies of galaxy mergers. The high resolution simulations of Hopkins et al. (2013) and the movies available online at http://www.tapir.caltech.edu/~phopkins/Site/Movies_sbw_mgr.html seem especially relevant, particularly their Milky Way and Sbc analogs. They predict a general but clumpy enhancement of star formation at around the time of first pericenter passage; a stronger, centrally concentrated starburst at the time of second pericenter passage, with a shorter period of separation before coalescence. Surprisingly perhaps, the merger timescales for both their MW and Sbc analogs are very similar to my SFH models, with ∼ 1 Gyr between first and second perigalactic passage and another few 10’s of Myr to final coalescence.

I’m going to wrap up with maps of the (posterior mean estimates) of star formation rate density (as \(\log10(\mathsf{M_\odot/yr/kpc^2})\)) and specific star formation rate, which has units log10(yr-1). These are 100Myr averages and recall are based solely on the SSP template contributions.

mrk848_sfr_ssfr
(L) Star formation rate density (R) Specific star formation rate

Several recent papers have made quantitative estimates of star formation rates in this system, which I’ll briefly review. Detailed comparisons are somewhat difficult because both the timescales probed by different SFR calibrators differ and the spatial extent considered in the studies varies, so I’ll just summarize reported results and compare as best as I can.

Yuan et al. (2018) modeled broad band UV-IR SED’s taking data from GALEX, SDSS, and Spitzer; and also did full spectrum fitting using ppxf on the same MaNGA data I’ve been examining. They divided the data into 5″ (≈ 4 kpc) radius regions covering the tidal tails and each nucleus. The broadband data were fit with parametric SFH models with a pair of exponential bursts (one fixed at 13Gyr, the other allowed to vary). From the parametric models they estimated the SFR in the northern and southern nuclei as 76 and 11 \(\mathsf{M_\odot/yr}\) (the exact interpretation of this value is unclear to me). For comparison I get values of 33 and 13 \(\mathsf{M_\odot/yr}\) for regions within 4 kpc of the two nuclei by calculating the average (posterior mean) star forming density and multiplying by π*42. They also calculated 100Myr average SFRs from full spectrum fitting with several different dust models and also from the Hα luminosity, deriving estimates that varied by an order of magnitude or more. Qualitatively we reach similar conclusions that the tails had earlier starbursts and are now forming stars at lower rates than their peak, and also that the northern nucleus has recent star formation a factor of several higher than the southern.

Cluver et al. (2017) were primarily trying to derive SFR calibrations for WISE W3/W4 filters using samples of normal star forming and LIRG/ULIRGs (including this system) with Spitzer and Herschel observations. Oddly, although they tabulate IR luminosities for the entire sample they don’t tabulate SFR estimates. But plugging into their equations 5 and 7 I get star formation rate estimates of just about 100 \(\mathsf{M_\odot/yr}\). These are global values for the entire system. For comparison I get a summed SFR for my models of ≈ 45 \(\mathsf{M_\odot/yr}\) (after making a 0.2 dex adjustment for fiber overlap).

Tsai and Hwang (2015) also used IR data from Spitzer and conveniently present results for quantities I track using the same units. Their estimates of the (log) stellar mass density, SFR density, and specific star formation rate in the central kpc of (presumably) the northern nucleus were 9.85±0.09, 0.43±0.00, and -9.39±0.09 respectively. For comparison in the fiber that covers the northern nucleus my models return 9.41±0.02, 0.57±0.04, and -8.84±0.04. For the 6 fibers closest to the nucleus the average stellar mass density drops to 9.17 and SFR density to about 0.39. So, our SFR estimates are surprisingly close while my mass density estimate is as much as 0.7 dex lower.

Finally, Vega et al. (2008) performed starburst model fits to broadband IR-radio data. They estimated a burst age of about 60Myr for this system, with average SFR over the burst of 225 \(\mathsf{M_\odot/yr}\) and current (last 10 Myr) SFR 87 \(\mathsf{M_\odot/yr}\) in a starforming region of radius 0.27kpc. Their model has optical depth of 33 at 1 μm, which would make their putative starburst completely invisible at optical wavelengths. Their calculations assumed a Salpeter IMF, which would add ≈0.2 dex to stellar masses and star formation rates compared to the Kroupa IMF used in my models.

Overall I find it encouraging that my model SFR estimates are no worse than factors of several lower than what are obtained from IR data — if the Vega+ estimate of the dust optical depth is correct most of the star formation is well hidden. Next time I plan to look at Hα emission and tie up other loose ends. If time permits before I have to drop this for a while I may look at two component dust models.

Markarian 848

This galaxy1also known as VV705, IZw 107, IRAS F15163+4255, among others has been my feature image since I started this blog. Why is that besides that it’s kind of cool looking? As I’ve mentioned before I took a shot a few years ago at selecting a sample of “transitional” galaxies from the SDSS spectroscopic sample, that is ones that may be in the process of rapidly shutting down star formation2See for example Alatalo et al. (2017) for recent usage of this terminology.. I based the selection on a combination of strong Hδ absorption and either weak emission lines or line ratios other than pure starforming, using measurements from the MPA-JHU pipeline. This galaxy pair made the sample based on the spectrum centered on the northern nucleus (there is also a spectrum of the southern nucleus from the BOSS spectrograph, but no MPA pipeline measurements). Well now, these galaxies are certainly transitioning to something, but they’re probably not shutting down star formation just yet. Simulations of gas rich mergers generally predict a starburst that peaks around the time of final coalescence. There is also significant current star formation, as high as 100-200 \(M_\odot/yr\) per various literature estimates, although it is mostly hidden. So on the face of it at least this appears to be a false positive. This was also an early MaNGA target, and one of a small number where both nuclei of an ongoing merger are covered by a single IFU:

Markarian 848. SDSS image thumbnail with MaNGA IFU footprint.

Today I’m going to look at a few results of my analysis of the IFU data that aren’t too model dependent to get some insight into why this system was selected. As usual I’m looking at the stacked RSS data rather than the data cubes, and for this exercise I Voronoi binned the spectra to a very low target SNR of 6.25. This leaves most of the fibers unbinned in the vicinity of the two nuclei. First, here is a map of BPT classification based on the [O III]/Hβ vs. [N II]/Hα diagnostic as well as scatterplots of several BPT diagnostics. Unfortunately the software I use for visualizing binned maps extends the bins out to the edge of the bounding rectangle rather than the hexagonally shaped edge of the data, which is outlined in red. Also note that different color coding is used for the scatter plots than the map. The contour lines in the map are arbitrarily spaced isophotes of the synthesized R band image supplied with the data cube.

Map of BPT class determined from [O III]/Hβ vs. [N II}/Hα diagnostic diagram, and BPT diagnostics for [N II], [S II] and [O I]. Curves are boundary lines form Kewley et al. (2006).

What’s immediately obvious is that most of the area covered by the IFU including both nuclei fall in the so-called “composite” region of the [N II]/Hα BPT diagram. This gets me back to something I’ve complained about previously. There was never a clear physical justification for the composite designation (which recall was first proposed in Kauffmann et al. 2003), and the upper demarcation line between “pure” AGN and composite systems as shown in the graph at upper right was especially questionable. It’s now known if it wasn’t at the turn of the century (which I doubt is the case) that a number of ionization mechanisms can produce line ratios that fall generally in the composite/LINER regions of BPT diagnostics. Shocks in particular are important in ongoing mergers. High velocity outflow of both ionized and neutral gas have been observed in the northern nucleus by Rupke and Veillux (2013), which they attributed to supernova driven winds.

The evidence for AGN in either nucleus is somewhat ambiguous. Fu et al. (2018) called this system a binary AGN, but that was based on the “composite” BPT line ratios from the same MaNGA data as we are examining (their map, by the way, is nearly identical to mine; see also Yuan et al. 2018). By contrast Vega et al. 2008 were able to fit the entire SED from NIR to radio frequencies with a pure starburst model and no AGN contribution at all, while more recently Dietrich et al. 2018 estimated the AGN fraction to be around 0.25 from NIR to FIR SED fitting. A similar conclusion that both nuclei contain both a starburst and AGN component was reached by Vardoulaki et al. 2015 based on radio and NIR data. One thing I haven’t seen commented on in the MaNGA data that possibly supports the idea that the southern nucleus harbors an AGN is that the regions with unambiguous AGN-like optical emission line ratios are fairly symmetrically located on either side of the southern nucleus and separated from it by ∼1-2 kpc. This could perhaps indicate the AGN is obscured from our view but visible from other angles.

There are also several areas with starforming emission line ratios just to the north and east of the northern nucleus and scattered along the northern tidal tail (the southern tail is largely outside the IFU footprint). In the cutout below taken from a false color composite of the F435W+F814W ACS images several bright star clusters can be seen just outside the more heavily obscured nuclear region, and these are likely sources of ionizing photons.

mrk_848_hst_crop
Markarian 848 nuclei. Cutout from HST ACS F814W+F435W images.

Finally turning to the other component of the selection criteria, here is a map of the (pseudo) Lick HδA index and a plot of the widely used HδA vs Dn(4000) diagnostic. It’s a little hard to see a clear pattern in the map because this is a rather noisy index, but strong Balmer absorption is seen pretty much throughout, with the highest values outside the two nuclei and especially along the northern tidal tail.

Location in the Hδ – D4000 plane doesn’t uniquely constrain the star formation history, but the contour plot taken from a largish sample of SDSS spectra from the NGP is clearly bimodal, with mostly starforming galaxies at upper left and passively evolving ones at lower right, with a long “green valley” in between. Simple models of post-starburst galaxies will loop upwards and to the right in this plane as the starburst ages before fading towards the green valley. This is exactly where most of points in this diagram lie, which certainly suggests an interesting recent star formation history.

mrk848_hda_d4000
(L) Map of HδA index. (R) HδA vs Dn(4000) index. Contours are for a sample of SDSS spectra from the north galactic pole.

I’m going to end with a bit of speculation. In simulations of gas rich major mergers the progenitors generally complete 2 or 3 orbits before final coalescence, with some enhancement of star formation during the perigalactic passages and perhaps some ebbing in between. This process plays out over hundreds of Myr to some Gyr. What I think we are seeing now is the 2nd or third encounter of this pair, with the previous encounters having left an imprint in the star formation history.

I’ve done SFH modeling for this binning of the data, and also for data binned to higher SNR and modeled with higher time resolution models. Next post I’ll look at these in more detail.

Revisiting the Baryonic Tully-Fisher relation… – Part 3

This post has been languishing in draft form for well over a month thanks to travel and me losing interest in the subject. I’m going to try to get it out of the way quickly and move on.

Last time I noted the presence of apparent outliers in the relationship between stellar mass and rotation velocity and pointed out that most of them are due to model failures of various sorts rather than “cosmic variance.” That would seem to suggest the need for some sample refinement, and the question then becomes how to trim the sample in a way that’s reproducible.

An obvious comment is that all of the outliers fall below the general trend and (less obviously perhaps) most have very large posterior uncertainties as well. This suggests a simple selection criterion: remove the measurements which have a small ratio of posterior mean to posterior standard deviation of rotation velocity. Using the asymptotic circular velocity v_c in the atan mean function and setting the threshold to 3 standard deviations the sample members that are selected for removal are circled in red below. This is certainly effective at removing outliers but it’s a little too indiscriminate — a number of points that closely follow the trend are selected for removal and in particular 19 out of 52 points with stellar masses less than \(10^{9.5} M_\odot\) are selected. But, let’s look at the results for this trimmed sample.

lgm_logvc_circled_bad_2ways
Posterior distribution of asymptotic velocity `v_c` vs stellar mass. Circled points have posterior mean(v_c)/sd(v_c) < 3.

Again I model the joint relationship between mass and circular velocity using my Stan implementation of Bovy, Hogg, and Roweis’s “Extreme deconvolution” with the simplification of assuming gaussian errors in both variables. The results are shown below for both circular velocity fiducials. Recall from my previous post on this subject the dotted red ellipse is a 95% joint confidence interval for the intrinsic relationship while the outer blue one is a 95% confidence ellipse for repeated measurements. Compared to the first time I performed this exercise the former ellipse is “fatter,” indicating more “cosmic variance” than was inferred from the earlier model. I attribute this to a better and more flexible model. Notice also the confidence region for repeated measurements is tighter than previously, reflecting tighter error bars for model posteriors.

tf_subset1
Joint distribution of stellar mass and velocity by “Extreme deconvolution.” Inner ellipse: 95% joint confidence region for the intrinsic relationship. Outer ellipse: 95% confidence ellipse for new data. Top: Asymptotic circular velocity v_c. Bottom: Circular velocity at 1.5 r_eff.

Now here is something of a surprise: below are the model results for the full sample compared to the trimmed one. The red and yellow ellipses are the estimated intrinsic relations using the full and trimmed samples, while green and blue are for repeated measurements. The estimated intrinsic relationships are nearly identical despite the many outliers. So, even though this model wasn’t formulated to be “robust” as the term is usually understood in statistics in practice it is, at least as regards to the important inferences in this application.

tf_alldr15
Joint distribution of stellar mass and velocity by “Extreme deconvolution” (complete sample). Inner ellipse: 95% joint confidence region for the intrinsic relationship. Outer ellipse: 95% confidence ellipse for new data. Top: Asymptotic circular velocity v_c. Bottom: Circular velocity at 1.5 r_eff.

Finally the slope, that is the exponent in the stellar mass Tully-Fisher relationship \(M^* \sim V_{rot}^\gamma\) is estimated as the (inverse of) slope of the major axis of the inner ellipses in the above plots. The posterior mean and 95% marginal confidence intervals for the two velocity measures and both samples are:

v_c (subset 1) \(4.81^{+0.28}_{-0.25}\)

v_c (all) \(4.81^{+0.28}_{-0.25}\)

v_r (subset 1) \(4.36^{+0.23}_{-0.20}\)

v_r (all) \(4.33^{+0.23}_{-0.21}\)

Does this suggest some tension with the value of 4 determined by McGaugh et al. (2000)? Not necessarily. For one thing this is properly an estimate of the stellar mass – velocity relationship, not the baryonic one. Generally lower stellar mass galaxies will have higher gas fractions than high stellar mass ones, so a proper accounting for that would shift the slope towards lower values. Also, and as can be seen here, both the choice of fiducial velocity and analysis method matter. This has been discussed recently in some detail by Lelli et al. (2019)1These two papers have two authors in common..

Next time, back to star formation history modeling.

Revisiting the Baryonic Tully-Fisher relation… – Part 2

Last time I left off with the remark that while most of the sample of disk galaxies clearly exhibits a tight relationship between circular velocity and stellar mass, there are some apparent outliers as well. While some “cosmic variance” is expected most of the apparent outliers are due to model failures, which have several possible causes:

  1. Violation of the physical assumptions of the model, namely that the stars and gas are rotating (together) in the plane of a thin disk that’s moderately inclined to our line of sight (see my original post on this topic).
  2. Errors in the photometry. I use two photometric quantities (specifically nsa_elpetro_ba and nsa_elpetro_phi) from the MaNGA DRPALL catalog to set priors for the kinematic parameters cos_i and phi(cosine of the disk inclination and angle to receding side) and also to initialize the sampler. Since proper and in practice fairly informative priors are required for these parameters errors here will propagate into the models, sometimes in ways that are fatal. I’ll look in more detail at some examples below.
  3. Bad velocity data.
  4. Not enough data.
  5. Sampler convergence failures with no obvious cause.

The first two bullet points are closely related: most of the failures to satisfy the physical assumptions are directly related to errors in the photometric decompositions. One fairly common failure was galaxies that were too nearly face on to obtain reliable rotation curves. As an example here is the galaxy with the lowest estimated rotation velocity in the sample, mangaid 1-135054 (plateifu 8550-12703):

8550-12703_vf_vrot
Mangaid 1-135054 (plateifu 8550-12703). (L) Measured velocity field and (R) posterior predictive estimate of circular velocity with 95% confidence band.

Besides showing no sign of rotation the velocity field hints at possible large scale, low velocity outflow in the central region. There are also a few apparent outliers, although these had little effect on model results. Fortunately the model output gives us plenty of clues that the results shouldn’t be trusted. The median circular velocity estimate is unrealistically low with very large posterior uncertainty (above right), while the posterior marginal density for cos_ihas a mode near 1 and also very large uncertainty (below).

8550-12703_post_cosi
Mangaid 1-135054 (plateifu 8550-12703). Posterior distribution of cosine of disk inclination.

Zooming out a bit on SDSS imaging gives a likely explanation for the peculiar velocity field. The elliptical galaxy just to the NW has nearly the same redshift (the velocity difference is ∼75 km/sec) and is almost certainly interacting with our target.

MaNGA target and companion; credit SDSS

A key assumption in using photometric properties as proxies for kinematic quantities is that disk galaxies have intrinsically circular surface brightness profiles. This is never quite the case in practice and sometimes morphological features like strong bars can make this assumption catastrophically wrong. Here was perhaps the most extreme example in DR14:

mangaid 1-185287 (plateifu 8252-12704). SDSS thumbnail with IFU overlay
8252-12704_vf
mangaid 1-185287 (plateifu 8252-12704) Measured velocity field from stacked RSS spectra

The photometric major axis angle was estimated to be 98.4o, that is just south of east, while the position angle of the maximum recession velocity is close to due south. When I first examined this velocity field I had a hard time reconciling it with rotation in a thin disk. This was before I learned how to do Voronoi binning though. The image below shows the binned velocity field (with a target S/N of 6.25). This shows that the relative velocity does increase on a roughly north to south line, indicating that this is indeed a rotating disk galaxy.

8252-12704_vf_binned
mangaid 1-185287 (plateifu 8252-12704) Measured velocity field from binned RSS spectra. Black arrow indicates major axis position angle from photometry. Gray arrow: position angle of receding side from velocity model with prior guess of 180o

I mentioned in an earlier post that without a proper prior on the kinematic position angle phi these models are inherently multi-modal (in fact they are infinitely modal and therefore would have improper posteriors). The solution to that of course is to have a proper prior. But if the prior is seriously in error the posterior estimates for the components of the velocity will end up scrambled as can be seen in the top row of the graph below, which shows the posterior distributions of the circular and “expansion” velocities1Remember that the photometric position angle is determined modulo π while the direction to the maximum recession velocity is measured modulo 2π. We don’t care about a π radian error in the prior though because that just flips the signs of the velocity components, which causes no sampling issues and is trivially fixable in the generated quantities block. It’s smaller errors that cause problems.

The obvious solution to a bad prior is to correct it2Using the data you’re trying to model to establish a prior is, technically, cheating (the polite term is “Empirical Bayes”), but this seems a relatively benign use of the data., which is easy enough. The bottom row shows the results of re-running the model with a prior phicentered on 180o and the same data. Now both the circular and expansion velocity curves are at least plausible. The posterior mean of phiis ≈185o, which is very close to correct as can be seen in the binned velocity field shown above.

8252-12704_rot_exp_guessphi_guess180
mangaid 1-185287 (plateifu 8252-12704) Top row: model rotation (L) and expansion (R) velocities with prior for major axis angle taken from photometry (phi = 98.4). Bottom row: same but with prior phi = 180o.

One more example. Bars were the most common cause of misleading photometric decompositions, but not the only one. Some galaxies are just asymmetrical. Here is the one that had the largest offset between the photometric and kinematic position angles:

mangaid 1-201291 (plateifu 8145-6103) – SDSS thumbnail with IFU footprint overlay

And the velocity field (again this agrees well with the Marvin measurements of the data cube):

8145-6103_velocityfield
mangaid 1-201291 (plateifu 8145-6103). Velocity field from stacked RSS spectra with major axis angle from nsa_elpetro_phi

This time the velocity field looks unremarkable, but again because of the prior the estimated circular and expansion velocities are scrambled together. And once again also, changing the prior to be centered on the approximate actual position angle of the receding side produces reasonable estimates for both:

mangaid 1-201291 (plateifu 8145-6103).
mangaid 1-201291 (plateifu 8145-6103). Top: Posterior predictive distributions of circular and expansion velocity with prior on `phi` from photometry. Bottom: Same with prior centered on -10o

Next up, I’ll take a more holistic look at final sample selection, and maybe get to results.

Revisiting the Baryonic Tully-Fisher relation with DR15 data

As I mentioned in the previous two posts SDSS Data Release 15 went public back in December and a query for “normal” disk galaxies as judged by Galaxy Zoo 2 classifiers returned 588 hits. I’ve finally run the GP velocity models on all of the new data and made a second run on around 40 that were contaminated by foreground stars or neighboring galaxies. So far I haven’t found an alternative to selecting these by eye and doing the masking manually, so that’s an error prone process. The month+ long gap between postings by the way was due to travel — my computer wasn’t grinding away on these models for all that time. As I mentioned last post the sampling properties including execution time of the GP model with arctangent mean function are usually quite favorable using Stan. The median wall time for these runs was about a minute, with a range from 25 to 1600 seconds. All model runs used 500 warmup iterations and 500 post-warmup with 4 chains run in parallel. This is more than enough for inference.

Before I discuss the results I’ll show them. As I did for the first pass at this way back in July I retrieved stellar mass and uncertainty estimates made by the MPA-JHU group from CasJobs; all but a handful also have mass estimates from the Wisconsin group. I may look at those later but don’t anticipate any very significant differences.

There are now at least two plausible choices for reference circular velocities: the velocity at a fiducial radius, and again I choose 1.5 effective radii since the MaNGA IFUs are meant to cover out to that radius in the primary sample. The other obvious choice is the asymptotic velocity vc in the arctangent mean function. This seems in principle to be the better choice since it estimates the circular velocity in the flat part of the rotation curve, but it might be a considerable extrapolation from the actual data in some cases.

Both sets of results are shown below for all model runs that ran to completion (N=582). Plotted are the median, 2.5, and 97.5 percentiles (≈ ± 2σ) of posterior predictions for the log-circular velocity at 1.5reff (top graph) and the same quantiles of the posteriors of the parameter v_c (bottom graph). These are plotted against the median, 16, and 84 percentiles (≈ ± 1σ) of the total stellar mass estimates per the MPA group.

logvr_mstar_dr15
Estimated rotation velocity at 1.5 effective radii vs. stellar mass estimate from MPA-JHU models
logvc_mstar_dr15
Estimated asymptotic rotation velocity against stellar mass from MPA-JHU models. Vertical error bars mark the 2.5 and 97.5 percentiles of the model posteriors of (log) velocity in km/sec. Horizontal error bars mark the 16 and 84 percentiles model posteriors of (log) stellar mass.

Evidently most of the sample follows a tight linear relationship with either measure of circular velocity, but there are some apparent outliers as well. I’m feeling a bit blocked right now, so I’ll end the post here. Next time I’ll look at some of the causes of model failure, what to do about them, and get to the results.

Yet more on rotation curve modeling — why the mean function matters

When I first began modeling disk galaxy rotation curves using low order polynomials for the circular velocity I noticed two rather frequent systematics in the model residuals:

  1. Lobe like areas symmetrically located around the nucleus with approximately equal and opposite signs. Sometimes these are co-located with bar ends but a bar is not always obvious.
  2. A contrast of a few 10’s of kilometers/sec between spiral arms and interarm regions. This is rather common in grand design spirals.

Here’s a particularly dramatic example of symmetrical lobes in mangaid 1-339041 (plateifu 8138-12704), IAU name SDSS J074637.70+444725.8. First, here are the measured line of sight velocities for the fiber spectra:

vf_rss_8138-12704
(L) Velocity field measured from stacked RSS file (R) Interpolated velocity field mangaid 1-339041 (plateifu 8138-12704)

The left plot shows the actual measurements from the stacked RSS file. The right is just an interpolated version of the left. Since value added data is now available it’s worth comparing this to output from “Marvin“. For reference here is the Hα velocity map:

Hα velocity map from MaNGA DAP

It’s hard to tell in any detail, but these look similar enough and the stellar velocity field as measured by the DAP also looks similar.

Next, here are the mean residuals from the posterior predictive fits shown as interpolated maps derived from the fits at the observed positions. As promised the left hand map from the low order polynomial fit has prominent lobes situated on either side of the nucleus and a more subtle contrast between spiral arms and interarm regions. The right hand map from the GP model appears to be largely free of systematic patterns. Why the difference?

vf_res_2ests_8138-12704
Mean residuals from models: (L) Polynomial rotation curve model (R) GP model with atan mean function

In this case the arctangent mean function I introduced in the last post worked very well, with the estimated circular velocity rising quickly to an asymptotic value of ∼300km/s. The low order polynomial representation is necessarily constrained by the possible shapes of a low order polynomial (in this case cubic), resulting in a shallower initial slope and a first local maximum farther out than in the GP model. The lobed residuals in the polynomial model are therefore seen to be due to an inner disk that’s rotating more rapidly than can be modeled (and not due to a kinematically distinct component or to streaming material).

rot_curves_2ests_8138-12704
Rotation and “expansion” velocity curves (T) polynomial model (B) GP model with atan mean function

As a brief morphological note, GZ2 classifiers thought this was a normal looking disk galaxy by an overwhelming majority. It’s hard to say they were wrong based on the SDSS imaging, but the deeper and wider legacy survey thumbnail clearly shows the outer disk to be disturbed, presumably by the edge on disk galaxy to the north — the relative velocities are ∼340km/sec, so they are likely in close proximity.

Plateifu 8138-12704 IFU footprint

As time allows I may take a closer look at model diagnostics available in Stan or some more examples. Longer term I plan to take another look at the Baryonic Tully-Fisher relationship for the larger sample available in DR15.

Rotation curves revisited; SDSS DR15; miscellany

I haven’t had a lot of time for this hobby lately and won’t for another month or so, but I haven’t been completely inactive. When I first started modeling disk galaxy rotation curves with Gaussian processes I used a low order polynomial for the mean function. This was an admitted hack, and not a very good one at that. It’s not uncommon in the rotation curve modeling industry to adopt a functional form for the circular velocity. Sometimes it’s physically motivated (usually based on a bulge/disk/dark matter halo model), sometimes it’s purely phenomenological. With a bit of digging in the literature I found a particularly simple form in Courteau (1997):

\(v(r) = v_{0} + \frac{2}{\pi} v_{c} \arctan(r/r_{t})\)

Making this the mean function for the GP model involves just a few simple changes to the Stan code I outlined in that post. We only need two parameters for the curve, and the parameter block now looks like:

parameters {       

real phi;

real<lower=0., upper=1.> cos_i; // cosine disk inclination

real x_c; //kinematic centers

real y_c;

real v_sys; //system velocity offset (should be small)

real v_c;

real<lower=0> r_t;

real<lower=0.> sigma_los;

real<lower=0.> alpha;

real<lower=0.> rho;

}

The key line in the model block is:

  v ~ multi_normal_cholesky(v_sys + sin_i * (2./pi() * v_c * atan(r/r_t) .* yhat ./ r), L_cov);       

And then in the generated quantities block are posterior predictions for new data at the observed points and in concentric rings:

  // model and residuals at the original points       



v_rot = 2./pi() * v_c * atan(r/r_t);

vf_model = gp_pred_rng(xyhat, v-v_sys-sin_i * (v_rot .* yhat ./ r), xyhat, alpha, rho, sigma_los);

vf_model = vf_model + v_sys + sin_i * (v_rot .* yhat ./r);

vf_res = v - vf_model;



// model in rings



v_gp = gp_pred_rng(xy_pred, v-v_sys-sin_i * (v_rot .* yhat ./ r), xyhat, alpha, rho, sigma_los);

v_pred = (v_gp + sin_i * 2./pi() * v_c * atan(r_pred/r_t) .* y_pred ./ r_pred) *v_norm/sin_i;

vrot_pred[1] = 0.;

vexp_pred[1] = 0.;

The full code is available in my github repository at https://github.com/mlpeck/vrot_stanmodels.

This model turned out to have very nice sampling properties using Stan’s version of NUTS, and convergence is fast enough in terms of both number of iterations and walltime to make it easily feasible to model a large sample.

Since I’ve been using this as a running example here are posterior circular and expansion velocity estimates for plateifu 8942-12702 (mangaid 1-218280):

gp_rot_exp_8942-12702
Circular and Expansion velocities, mangaid 1-218280 estimated by GP model with atan mean function

These curves are shaped just slightly differently from the ones I posted in July, but the confidence bounds are somewhat tighter. Whether this is actually due to better inference remains to be seen.

Second topic: SDSS Data release 15 went public in December as promised. I haven’t had time to dig into it much, but I did rerun the query I used to get a sample of high confidence disk galaxies based on GZ2 classifications. Run in the DR15 context that query got 588 hits, an increase of 229 from DR14 and still about 13% of the entire MaNGA galaxy sample excluding ancillary targets that aren’t also part of the main samples.

I downloaded the RSS files and have done a first pass at rotation curve modeling of the 229 new galaxies. Unfortunately I won’t have time to actually analyze what I have until, most likely, March. To conclude for now here is the kinematically most interesting example I saw from the new to DR15 sample. GZ2 classifiers called this a normal disk galaxy by a large majority


SDSS J114526.10+495244.5
legacy survey cutout

and were about evenly divided on whether it has a bar. And here is its velocity field measured from the stacked RSS file on the left, with the modeled velocity field on the right. The left graph is interpolated to the same resolution as the data cubes.

vf_8989-9101
Interpolated velocity field and GP based model, mangaid 1-174947 (plateifu 8989-9101)

This is an incomplete map because of the signal/noise threshold I set, but it’s apparent that the outer disk is counter-rotating relative to the inner disk. This is confirmed for the stellar velocity field in the SDSS DAP, while the ionized gas rotates in the same sense as the outer disk throughout. Keep in mind that I measure blended velocities that in any given case could be dominated either by stellar light or ionized gas — in this galaxy the emission lines are very weak, so I measured the stellar velocity field. The GP based model does a good job of fitting the data, but the rotation curve isn’t exactly the standard shape:

vrot_vexp_8989-9101
Modeled circular and expansion velocity

One thing this demonstrates is the versatility of the gaussian process model with a simple mean function (even if it’s wrong).

A random pair of passively evolving galaxies

This post will be short. Eventually I’d like to add a reproducible real data example or two to the documentation of the software I just officially published. Even though I’m mostly working with MaNGA data these days it may not be practical to include an example since the smallest RSS file is ∼20MB and it can take days to fully process a single data set. The single fiber SDSS spectra on the other hand are only ~170-210KB in “lite” form. While I’m hunting around for suitable examples I’ll post a few results here. For today’s post a couple of randomly selected early type galaxies from my north galactic pole sample:

J1218+2802
SDSS J121810.72+280230.4

 

J1221+2221
SDSS J122140.20+222108.3

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

These don’t seem to be in any historical catalogs, so I’ll just call them by their SDSS designations of J121810.72+280230.4 and J122140.20+222108.3, or J1218+2802 and J1221+2221 for shorthand. The first appears to me to be a garden variety elliptical; the second might be an S0 — neither has a morphological classification listed in NED (a majority of GZ classifiers in both GZ1 and GZ2 called the second an elliptical). Galaxy 2 has a redshift of z=0.023, which puts it in or near the Coma cluster. Galaxy 1 is well in the background at z=0.081.

I modeled the star formation histories with my standard tools and the “full” EMILES ssp library (that is all 54 time bins in 4 metallicity bins). The Stan models both completed without warnings (one required some adjustment of control parameters). Here are some basic results. First, star formation histories, displayed as star formation rate against lookback time. I find it can be useful to use linear or logarithmic scales on either or both axes. Here I show SFR scaled both ways:

sfh_twogalaxies
Star formation histories for 2 passively evolving galaxies

So-called non-parametric models for star formation histories are hoped for gold standards for SFH modeling, and here we see part of the reason: these have rather different estimated mass assembly histories despite similar looking spectra (see below or follow the links), and neither estimate is particularly close to commonly chosen functional forms.

There are some oddities to beware of however. For example literally 100% of the models I’ve run with any variation of the EMILES library show an uptick in mean star formation rate at 4Gyr, independent of the absolute estimated SFR. This can’t be real obviously enough, but I haven’t yet found an error or other specific reason for it. There may be other persistent zigs and zags in SFR at specific times, and these can vary with the SSP library. Overall I assign no particular significance to small variations in mean estimated star formation rates, especially at early times. I’m somewhat more inclined to take seriously late time variations (say, most recent Gyr), which may be wishful thinking.

 

The data aren’t inconsistent with a small amount of ongoing star formation (estimates of specific star formation rate averaged over 100Myr are shown — these have units of yr-1). SSFR’s in a range ~10-11 yr-1 are about an order of magnitude lower than an ordinary starforming galaxy in the local universe.


ssfr_twogalaxies
froPosterior distribution of specific star formation rate (100 Myr timescale)

Fits to the data are comparable, and satisfactory over most of the wavelength range of the spectra. A recent paper on arxiv (article id 1811.09799) mentioned some fitting issues with EMILES SSP libraries in the red. This can be seen here in the range ~6750-7500Å (rest frame). I don’t have an explanation either. It might be significant that the MILES stellar library spectra extend only to ~7300Å, with other data sources grafted on to the red extension, so there could be a calibration problem. On the other hand the prominent dip in this range is at least partly due to TiO absorption, so this could point to metal abundance variations not captured by EMILES.


 

Finally, I find it useful to model attenuation, even in ETGs that are often considered dust free. I just use a Calzetti attenuation curve in these models even though it was based on starburst galaxy SEDs and therefore may not be appropriate. This is a topic for possible future model refinement.


tauv_twogalaxies
Posterior estimates of dust attenuation

Software officially published and version controlled

It took longer than I had hoped, but I finally have the R and Stan software that I use for star formation history modeling under version control and uploaded to github. There are 3 R packages, two of which are newly available:

    • spmutils: The actual code for maximum likelihood and Bayesian modeling of star formation histories. There are also tools for reading FITS files from SDSS single fiber spectra and from MaNGA, visualization, and analysis.
    • cosmo: A simple, lightweight cosmological distance calculator.
    • sfdmap: A medium resolution version of the all sky dust map of Schlegel, Finkbeiner and Davis (SFD) and a pair of functions for rapid access.

So far only cosmo is documented to minimal R standards. It’s also the most generally useful, at least if you need to calculate cosmological distances. sfdmap has skeleton help files; spmutils not even that. I plan to correct that, but probably not in the next several weeks.

The R packages are pure R code and can be installed directly from github using devtools:

cosmo: devtools::install_github("mlpeck/cosmo")

sfdmap: devtools::install_github("mlpeck/sfdmap")

spmutils: devtools::install_github("mlpeck/spmutils/spmutils")

The Stan files for SFH modeling that I’m currently working with are in spmutils/stan. I keep the Stan files in the directory ~/spmcode on my local machines, and that’s the default location in the functions that call them.

This code mushroomed over a period of several years without much in the way of design, so it’s not much use without some documentation of typical workflow. I will get around to that sometime soon (I hope).