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.