Download this as a Jupyter notebook!

Spatially-Resolved Mass-Metallicity Relation

We’re going to construct the spatially-resolved mass-metallicity relation (MZR) for a MaNGA galaxy, where mass refers to stellar mass and metallicity refers to gas-phase oxygen abundance.

Roadmap

1. Compute metallicity.

2. Select spaxels that are

3. star-forming,

4. not flagged as “bad data,” and

5. above a signal-to-noise ratio threshold.

6. Compute stellar mass surface density.

7. Plot metallicity as a function of stellar mass surface density.

import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

import os
from os.path import join
path_notebooks = os.path.abspath('.')
path_data = join(path_notebooks, 'data')

Load Maps for Galaxy

Import the Marvin Maps class from marvin.tools.maps and initialize a Maps object for the galaxy 8077-6104.

from marvin.tools.maps import Maps

# REMOVE FROM NOTEBOOK
filename = '/Users/andrews/hacks/galaxies-mzr/data/manga-8077-6104-MAPS-SPX-GAU-MILESHC.fits.gz'
maps = Maps(filename=filename)

# maps = Maps('8077-6104')

INFO: No release version set. Setting default to MPL-6

Measure Metallicity

Pettini & Pagel (2004) N2 metallicity calibration

We are going to use the N2 metallicity calibration (their Equation 1) from Pettini & Pagel (2004):

12 + log(O/H) = 8.90 + 0.57 \(\times\) log( \(\frac{F([NII])}{F(H\alpha)}\) ).

One of the benefits of this calibration is that the required lines are very close in wavelength, so the reddening correction is negligible.

Get [NII] 6585 and Halpha flux maps from the Marvin Maps object. Note: MaNGA (and Marvin) use the wavelengths of lines in vaccuum, whereas they are usually reported in air, hence the slight offsets.

nii = maps.emline_gflux_nii_6585
ha = maps.emline_gflux_ha_6564

Calculate the necessary line ratio.

Marvin can do map arithmetic, which propagates the inverse variances and masks, so you can just do +, -, *, /, and ** operations as normal. (Note: taking the log of a Marvin Map will work for the values but the inverse variance propagation does not correctly propagate the inverse variance yet.)

n2 = nii / ha
logn2 = np.log10(n2)

Finally, calculate the metallicity.

oh = 8.90 + 0.57 * logn2

Select Spaxels

Using the BPT Diagram to select star-forming spaxels

Metallicity indicators only work for star-forming spaxels, so we need a way to select only these spaxels.

The classic diagnostic diagram for classify the emission from galaxies (or galactic sub-regions) as star-forming or non-star-forming (i.e., from active galactic nuclei (AGN) or evolved stars) was originally proposed in Baldwin, Phillips, & Terlevich (1981) and is known as the BPT diagram.

The BPT diagram uses ratios of emission lines to separate thermal and non-thermal emission.

The classic BPT diagram uses [OIII]5007 / Hbeta vs. [NII]6583 / Halpha, but there are several versions of the BPT diagram that use different lines ratios.

BPT Diagrams with Marvin

Let’s use Marvin’s maps.get_bpt() method to make BPT diagrams for this galaxy.

red line: maximal starbust (Kewley et al 2001) – everything to the right is non-star-forming.

dashed black line: conservative star-forming cut (Kauffmann et al. 2003) – everything to the left is star-forming.

Line ratios that fall in between these two lines are designated “Composite” with contributions from both star-forming and non-star-forming emission.

blue line: separates non-star-forming spaxels into Seyferts and LINERs.

Seyferts are a type of AGNs.

LINERs (Low Ionization Nuclear Emission Regions) are not always nuclear (LIER is a better acronym) and not always AGN (oftern hot evolved stars).

Sometimes these diagnostic diagrams disagree with each other, hence the “Ambiguous” designation.

Try using maps.get_bpt? to read the documentation on how to use this function.

masks_bpt, __, __ = maps.get_bpt()

The BPT masks are dictionaries of dictionaries of a boolean (True/False) arrays. We are interested in the spaxels that are classified as star-forming in all three BPT diagrams are designated as True, which is designated with the global key. Print this mask.

masks_bpt['sf']['global']

array([[False, False, False, ..., False, False, False],
       [False, False, False, ..., False, False, False],
       [False, False, False, ..., False, False, False],
       ...,
       [False, False, False, ..., False, False, False],
       [False, False, False, ..., False, False, False],
       [False, False, False, ..., False, False, False]])

Masks

MaNGA (and SDSS generally) use bitmasks to communicate data quality.

Marvin has built-in methods to convert from the bitmasks integer values to individual bits or labels and to create new masks by specifying a set of labels.

Show the mask schema with n2.pixmask.schema.

n2.pixmask.schema

	bit	label	description
0	0	NOCOV	No coverage in this spaxel
1	1	LOWCOV	Low coverage in this spaxel
2	2	DEADFIBER	Major contributing fiber is dead
3	3	FORESTAR	Foreground star
4	4	NOVALUE	Spaxel was not fit because it did not meet sel...
5	5	UNRELIABLE	Value is deemed unreliable; see TRM for defini...
6	6	MATHERROR	Mathematical error in computing value
7	7	FITFAILED	Attempted fit for property failed
8	8	NEARBOUND	Fitted value is too near an imposed boundary; ...
9	9	NOCORRECTION	Appropriate correction not available
10	10	MULTICOMP	Multi-component velocity features present
11	30	DONOTUSE	Do not use this spaxel for science

Select non-star-forming spaxels (from the BPT mask) and set their mask value to the DAP’s DONOTUSE value with the n2.pixmask.labels_to_value() method. Note that we are selecting spaxels that we want from the BPT mask (i.e., True is a spaxel to keep), whereas we are using the pixmask to select spaxels that we want to exclude (i.e., True is a spaxel to ignore).

mask_non_sf = ~masks_bpt['sf']['global'] * n2.pixmask.labels_to_value('DONOTUSE')

Select spaxels classified by the DAP as bad data according to the masks for spaxels with no IFU coverage, with unreliable measurements, or otherwise unfit for science. Use the n2.pixmask.get_mask method.

mask_bad_data = n2.pixmask.get_mask(['NOCOV', 'UNRELIABLE', 'DONOTUSE'])

Select spaxels with signal-to-noise ratios (SNRs) > 3 on both [NII] 6585 and Halpha.

ha.ivar = inverse variance = \(\frac{1}{\sigma^2}\), where \(\sigma\) is the error.

min_snr = 3.
mask_nii_low_snr = (np.abs(nii.value * np.sqrt(nii.ivar)) < min_snr)
mask_ha_low_snr = (np.abs(ha.value * np.sqrt(ha.ivar)) < min_snr)

Do a bitwise (binary) OR to create a master mask of spaxels to ignore.

mask = mask_non_sf | mask_bad_data | mask_nii_low_snr | mask_ha_low_snr

Plot the Metallicity Map

Plot the map of metallicity using the plot() method from your Marvin Map metallicity object. Also, mask undesirable spaxels and label the colorbar.

Note: solar metallicity is about 8.7.

fig, ax = oh.plot(mask=mask, cblabel='12+log(O/H)')

Compute Stellar Mass Surface Density

1. Read in spaxel stellar mass measurements from the Firefly spectral fitting catalog (Goddard et al. 2017).

• The Firefly stellar population fitting results file is large (1.8 GB), so we have extracted only the measurements needed for our purposes DOWNLOAD CSV FILE.

• Summary of the Firefly Value Added Catalog

• Datamodel of the Firefly Value Added Catalog

2. Convert spaxel angular size to a physical scale in pc.

3. Divide stellar mass by area to get stellar surface mass density.

Read in stellar masses

Use pandas to read in the csv file with stellar masses.

import pandas as pd
mstar = pd.read_csv(join(path_data, 'manga-{}_mstar.csv'.format(maps.plateifu)))

Plot stellar mass map using ax.imshow(). MaNGA maps are oriented such that you want to specify origin='lower'. Also include a labelled colorbar.

fig, ax = plt.subplots()
p = ax.imshow(mstar, origin='lower')
ax.set_xlabel('spaxel')
ax.set_ylabel('spaxel')
cb = fig.colorbar(p)
cb.set_label('log(Mstar) [M$_\odot$]')

Calculate physical size of a spaxel

MaNGA’s maps (and data cubes) have a spaxel size of 0.5 arcsec. Let’s convert that into a physical scale for our galaxy.

spaxel_size = 0.5  # [arcsec]

# or programmatically:
# spaxel_size = float(maps.getCube().header['CD2_2']) * 3600

Get the redshift of the galaxy from the maps.nsa attribute.

redshift = maps.nsa['z']

We’ll use the small angle approximation to estimate the physical scale:

\(\theta = \mathrm{tan}^{-1}(\frac{d}{D}) \approx \frac{206,265 \, \mathrm{arcsec}}{1 \, \mathrm{radian}} \frac{d}{D}\),

where

\(\theta\) is the angular size of the object (in our case spaxel) in arcsec,

\(d\) is the diameter of the object (spaxel), and

\(D\) is the angular diameter distance.

The distance (via the Hubble Law — which is fairly accurate for low redshift objects) is

\(D \approx \frac{cz}{H_0}\),

where

\(c\) is the speed of light in km/s,

\(z\) is the redshift, and

\(H_0\) is the Hubble constant in km/s/Mpc.

Calculate \(D\).

c = 299792  # speed of light [km/s]
H0 = 70  # [km s^-1 Mpc^-1]
D = c * redshift / H0  # approx. distance to galaxy [Mpc]

Rearrange the small angle formula to solve for the scale (\(\frac{d}{\theta}\)) in pc / arcsec.

scale = 1 / 206265 * D * 1e6  # 1 radian = 206265 arcsec [pc / arcsec]

Now convert the spaxel size from arcsec to parsecs and calculate the area of a spaxel.

spaxel_area = (scale * spaxel_size)**2  # [pc^2]

Finally, we simply divide the stellar mass by the area to get the stellar mass surface density \(\Sigma_\star\) in units of \(\frac{M_\odot}{pc^2}\).

sigma_star = np.log10(10**mstar / spaxel_area)  # [Msun / pc^2]

Let’s plot metallicity as a function of \(\Sigma_\star\)! Remember to apply the mask. Also set the axis range to be [0, 4, 8, 8.8].

fig, ax = plt.subplots(figsize=(6, 6))
ax.scatter(sigma_star.values[mask == 0], oh.value[mask == 0], alpha=0.15)
ax.set_xlabel('log(Mstar) [M$_\odot$]')
ax.set_ylabel('12+log(O/H)')
ax.axis([0, 4, 8.0, 8.8])

[0, 4, 8.0, 8.8]

MaNGA Spatially-Resolved Mass-Metallicity Relation

We have constructed the spatially-resolved MZR for one galaxy, but we are interested in understanding the evolution of galaxies in general, so we want to repeat this exercise for many galaxies. In Barrera-Ballesteros et al. (2016), Jorge Barrera-Ballesteros (who gave a talk at Pitt in November 2017) did just this, and here is the analogous figure for 653 disk galaxies.

The best fit line from Barrera-Ballesteros et al. (2016) is given in the next cell.

# fitting formula
aa = 8.55
bb = 0.014
cc = 3.14
xx = np.linspace(1, 3, 1000)
yy = aa + bb * (xx - cc) * np.exp(-(xx - cc))

Remake the spatially-resolved MZR plot for our galaxy showing the he best fit line from Barrera-Ballesteros et al. (2016).

fig, ax = plt.subplots(figsize=(6, 6))
ax.scatter(sigma_star.values[mask == 0], oh.value[mask == 0], alpha=0.15)
ax.plot(xx, yy)
ax.set_xlabel('log(Mstar) [M$_\odot$]')
ax.set_ylabel('12+log(O/H)')
ax.axis([0, 4, 8.0, 8.8])

[0, 4, 8.0, 8.8]

The spaxels in our galaxy are typically above the best fit relation. Part of the offset may be due to systematic differences in the metallicity calibrator used, but the overal trend of flat metallicity as stellar mass surface densities decreases seems to be in tension with their best fit. It would be worth investigating this effect for more galaxies to understand if individual galaxies typically obey the best fit relation or whether they typically exhibit a flat trend in this space.

Ultimately, Barrera-Ballesteros et al. (2016) concluded that the spatially-resolved MZR is a scaled version of the global MZR (e.g., see Tremonti et al. 2004).