The Data Analysis Pipeline for the SDSS-IV MaNGA IFU Galaxy Survey: Emission-line Modeling

The MaNGA survey is one of the three key components of SDSS-IV (Blanton et al. 2017) and aims to obtain IFS data for a representative sample of 10,000 galaxies in the redshift range 0.01 < z < 0.15 by 2020. The MaNGA instrument operates on the SDSS 2.5 m telescope at Apache Point Observatory (Gunn et al. 2006) and consists of a set of 17 hexagonal fiber bundles of different sizes, plus a set of minibundles and sky fibers used for flux calibration and sky subtraction, respectively (Drory et al. 2015; Law et al. 2015; Yan et al. 2016a, 2016b). All fibers are fed into the dual-beam BOSS spectrographs covering the wavelength range from 3600 to 10300 Å with a spectral resolution R ∼ 2000 (Smee et al. 2013).

MaNGA galaxies are selected from an extended version of the NASA-Sloan Atlas (NSA) and are observed out to 1.5 R_e (primary sample, comprising two-thirds of the total sample) or 2.5 R_e (secondary sample, comprising one-third of the total sample). Targets are selected to be representative of the overall galaxy population at each stellar mass in the range $9.0\lt \mathrm{log}({M}_{\star }/{M}_{\odot })\,\lt 11.0$. In practice, the absolute i-band magnitude is used for sample selection to avoid the systematic uncertainty intrinsic in deriving stellar masses (Wake et al. 2017).

The starting point of this paper is the datacubes produced for DR15 by the MaNGA data-reduction pipeline (DRP; Law et al. 2016), with additional modifications described in Aguado et al. (2019). The MaNGA DAP takes as input the reduced MaNGA datacubes generated with logarithmic wavelength sampling.

An overview of the DAP workflow is presented in Westfall et al. (2019). Here, we summarize the aspects of the algorithm that are most relevant to the derivation of emission-line properties, in order to motivate the tests performed in the rest of the paper. A graphical overview of the relevant components of the DAP workflow is presented in Figure 1 and discussed in the following sections.

Zoom In Zoom Out Reset image size

The DAP currently performs two full-spectrum fits, both using pPXF. The first fit primarily determines the stellar kinematics, and the second models the emission lines. Two example spectra fitted by the MaNGA DAP can be inspected in Figure 2. For DR15, the stellar kinematics are determined only for spectra binned to g-band S/N > 10, as produced by applying the Voronoi binning algorithm implemented in the python language¹⁶ by Cappellari & Copin (2003). The subsequent emission-line modeling is done for both the binned spectra and after deconstructing the bins into the individual spaxels; a full description of the emission-line-fitting module of the DAP is provided in Section 8 of Westfall et al. (2019). Because it is most relevant to emission-line science, we focus on the results provided by the latter approach, which we term the “hybrid” binning scheme (see below).

Zoom In Zoom Out Reset image size

2.2.1. Stellar-kinematics Fit (First Fitting Stage)

2.2.2. Hybrid Binning Approach

2.2.3. Simultaneous Fit of Gas and Stars (Second Fitting Stage)

To simultaneously optimize the fit to the stellar continuum and the emission lines, we use pPXF with Gaussian emission-line templates associated to kinematic parameters that are independent of those used for the stellar templates (Cappellari 2017). Simultaneous fits of emission-line and continuum templates have been introduced and recommended in previous work (Sarzi et al. 2006; Oh et al. 2011) as a way to minimize the bias resulting from masking of the wings of the Balmer absorption profiles and other stellar features on the recovered emission-line fluxes. This approach allows one to enforce the physical constraint that emission lines cannot be negative, while optimizing the fit to the stellar continuum. pPXF adopted the same idea for the gas fitting, but implemented this strategy in a different way, as described in Cappellari (2017). An example of the results of this fitting algorithm can be seen in Figure 2 for both a highly star-forming and an early-type galaxy. We assess the difference in the best-fit models derived by simultaneous versus subsequent fitting of the emission lines in Section 5.1.

The 22 emission lines fit by the DAP in DR15 are presented in Table 1. The flux ratios of doublets are fixed when such ratios are determined by atomic physics (Table 1; Osterbrock & Ferland 2006). For DR15 we tie the velocity of all the fitted emission lines, but do not tie the velocity dispersions with the exception of most of the doublets. Specifically, the line doublets with fixed flux ratios and the [O ii]λλ 3727,29 doublet have their velocity dispersions tied between the doublet lines, but each doublet and all other lines have independent velocity dispersions.¹⁷ While tying the kinematics of different emission lines may prove advantageous in recovering the flux of weak lines, we do not tie all velocity dispersions in DR15 to allow for modest inaccuracies (generally a few percent) in the wavelength-dependent line-spread function (LSF) determined by the MaNGA DRP; further discussion of the current line-tying strategy is presented in Section 5.2. Given the limited spectral range of the adopted MILES-HC stellar library, lines redder than ∼7400 Å are not fit for DR15. In future data releases, however, we plan to use templates with a larger wavelength range to allow for a determination of the continuum under the emission lines redder than 7400 Å (see Section 4).

Table 1.  Wavelengths and Ionization Potential of the Relevant Ion for Each Emission Line Fit for DR15, Subdivided into the Groups Defined to Study the Different Tying Schemes Described in Section 5.3

Line Name	Wavelength (Vacuum) (Å)	DAP String Name	Ionization Potential (eV)	Fixed Ratio
Hydrogen Balmer Lines
Hθ (H10)	3798.983	Hthe-3798	13.60	no
Hη (H9)	3836.479	Heta-3836	13.60	no
Hζ (H8)	3890.158	Hzet-3890	13.60	no
H (H7)	3971.202	Heps-3971	13.60	no
Hδ	4102.899	Hdel-4102	13.60	no
Hγ	4341.691	Hgam-4341	13.60	no
Hβ	4862.691	Hb-4862	13.60	no
Hα	6564.632	Ha-6564	13.60	no
Low-ionization lines
[O ii]λ3727	3727.092	OII-3727	13.61	no
[O ii]λ3729	3729.875	OII-3729	13.61	no
[O i]λ6300	6302.04	OI-6302	0.0	no
[O i]λ6364	6365.535	OI-6365	0.0	0.328 [O i]λ6300
[N ii]λ6548	6549.86	NII-6549	14.53	0.327 [N ii]λ6584
[N ii]λ6584	6585.271	NII-6585	14.53	no
[S ii]λ6717	6718.294	SII-6718	10.36	no
[S ii]λ6731	6732.674	SII-6732	10.36	no
High-ionization lines
[Ne iii]λ3869	3869.86	NeIII-3869	40.96	no
[Ne iii]λ3968	3968.59	NeIII-3968	40.96	no
He iiλ4687	4687.015	HeII-4687	54.41	no
[O iii]λ4959	4960.295	OIII-4960	35.12	0.340 [O iii]λ5007
[O iii]λ5007	5008.240	OIII-5008	35.12	no
He iλ5876	5877.243	HeI-5877	24.58	no

Note. Ritz wavelengths in vacuum are taken from the National Institute of Standards and Technology (NIST; http://physics.nist.gov/PhysRefData/ASD/Html/help.html). The DAP string name is reported in the header of the MAPS files for the emission-line extensions and allows the users to associate each map with the correct line (see Section 6.1). Ionization potentials are taken from Draine (2011). Lines redder than ∼7400 Å, corresponding to the red cutoff of the MILES-HC stellar templates, are not fit in DR15.

Download table as:  ASCIITypeset image

During the fitting procedure, the DAP does not constrain the emission-line fluxes to follow a specific attenuation law resulting from dust present in the host galaxy (as generally done in GANDALF, e.g., Oh et al. 2011, and optionally available in pPXF). All line fluxes are, however, corrected for Galactic foreground extinction using the maps of Schlegel et al. (1998) and the reddening law of O’Donnell (1994). Users comparing the DAP output with the output generated by other spectral-fitting pipelines may need to take this factor into account.

An eighth-order multiplicative Legendre polynomial is used in DR15 to match the overall spectral shape of the data, which can deviate from that of the models both because of dust extinction and small inaccuracies in the spectrophotometric calibration. A physically motivated extinction law may be used in the MaNGA DAP instead of multiplicative polynomials, but this was found to produce worse fits to the stellar continuum, especially at the blue end of the spectrum (see Section 5.1). Additive polynomials are not, and indeed ought not be, used in this stage as they modify the depth of stellar absorption lines in the templates, therefore potentially leading to degeneracies with emission-line strengths. This point is further elucidated in Section 5.1. Neither polynomials nor extinction corrections are applied to the emission-line templates by the fit, but only to the stellar continuum.

2.2.4. Nonparametric Emission-line Properties and EW

A full description of the DAP output data model is provided by Westfall et al. (2019), particularly in Sections 2 and 11 and Appendix A.

Emission-line fluxes (from both Gaussian fitting and the moments analysis), velocities, velocity dispersion, and associated errors and masks are consolidated into the main DAP output file, the MAPS file, for each analyzed datacube. The MAPS file is a multiextension fits file, where each extension provides a set of 2D maps of DAP measurements.¹⁸

The best-fit continuum and emission-line models are given as extensions in the DAP model LOGCUBE file. Most extensions in this file are three-dimensional datacubes, presented on the same world coordinate frame as the input MaNGA datacube. The model LOGCUBE provides the results of both the continuum-only fit used to determine the stellar kinematics and the combined fit used to simultaneously model the continuum and emission.¹⁹ We highlight here, however, that the continuum from the stellar-kinematics fit should not be used to recompute emission-line parameters.

While the vast majority of spaxels are successfully fit by the MaNGA DAP, users should be aware of some known failure modes, discussed in Section 10.2.2 of Westfall et al. (2019). In the context of line emission, it is particularly important to note that broad-line active galactic nuclei (AGNs) are generally not well fit. All current runs of the DAP assume a single Gaussian component per emission line, meaning that it is not currently possible to recover the broad and narrow components present in Type I AGNs. This represents a notable limitation of the current DAP, and we therefore recommend that users interested in Type I AGNs—roughly 1% of MaNGA galaxies (Sánchez et al. 2017)—perform their own spectral fitting. Further recommendations and descriptions of known bugs in the DR15 DAP run are presented in Section 6.1.

Defining line S/N and A/N. Several fit quality measures have been employed in the literature to describe the reliability of measured parameters for an emission line. The most-used characterization of fit quality is the fractional error on the recovered line flux ${\mathrm{Err}}_{\mathrm{Flux}}/\mathrm{Flux}$, which is often quoted in terms of the “signal-to-noise ratio” of a line, defined as

$$\begin{eqnarray}&&{\rm{S}}/{{\rm{N}}}_{\mathrm{line}}\equiv \mathrm{Flux}/{\mathrm{Err}}_{\mathrm{Flux}}.\end{eqnarray} \tag{ 1 }$$

A more empirical way to assess line detection relies on quantifying how much the line protrudes above the noise level in the spectrum. This is usually measured by the amplitude-over-noise ratio,

$$\begin{eqnarray}&&{\rm{A}}/{{\rm{N}}}_{\mathrm{line}}\equiv \mathrm{Amplitude}/\mathrm{rms},\end{eqnarray} \tag{ 2 }$$

where the amplitude refers to the best-fit Gaussian amplitude and the rms is calculated from the residuals between the data and the model in small regions on either side of the line.

Line S/N as a good measure of fit quality. In Figure 3 we plot the S/N_line obtained from the DAP MAPS file versus the A/N_line obtained using the DAP fit residuals around the position of the line for a sample of 300 random galaxies in DR15 (∼4 × 10⁴ spaxels). The rms is computed as the mean rms in side bands blueward and redward of the position of each line.²⁰

Zoom In Zoom Out Reset image size

The figure highlights the tight linear relation in log space between the two quantities across almost 3 dex in S/N_line for strong lines spanning a large fraction of the MaNGA wavelength range. This tight relation, which shows an increase in scatter only for S/N_line < 3, implies that the S/N computed by the DAP is equivalent to the more empirical A/N to very good accuracy. This scaling is indeed expected because most emission lines in MaNGA are unresolved, and their velocity dispersion is roughly comparable to one pixel in the spectral direction.

To check whether the fit residuals at the position of different emission lines are comparable with the error spectrum, we computed the χ² per degree of freedom (dof) of the Gaussian fit in 15 pixel windows around the fitted position of each line center. In the right panel of Figure 3, we show the median relation between the ${\chi }^{2}/\mathrm{dof}$ for each line and S/N_line. All strong lines considered follow a similar relation, except [O ii]λ3727, which suffers from worse χ² at fixed S/N_line, possibly as a result of the difficulty in correctly fitting the unresolved doublet. For S/N_line < 30, the other strong lines show a roughly constant ${\chi }^{2}/\mathrm{dof}\,\sim \,0.8$. At higher S/N_line, the ${\chi }^{2}/\mathrm{dof}$ increases sharply, up to three orders of magnitude. We interpret this as “template mismatch,” in the sense that our Gaussian model represents an increasingly worse representation of the observed line profiles at high S/N. In this regime, the χ² could be lowered by fitting each emission line with more than one Gaussian component (e.g., Gallagher et al. 2019), but this goes beyond the scope of the current MaNGA DAP. We note a similar behavior in the ${\chi }^{2}/\mathrm{dof}$ over the full spectrum in both full-spectrum-fitting modules of the DAP, as shown in Figure 19 of Westfall et al. (2019); however, the discrepancies between model and data are orders of magnitude stronger in these small windows near each line.

We note that this increase in ${\chi }^{2}/\mathrm{dof}$ does not mean that the fluxes obtained via Gaussian fitting are unreliable. In fact, we have checked that the fluxes obtained for Gaussian fitting agree exceedingly well with line fluxes obtained by simply summing the flux around the position of the line (SFLUX extension of the DAP MAPS files), as can be seen in Figure 4. For S/N < 3, we start to see a discrepancy between the two flux measurements, with Gaussian fluxes being higher on average. This effect is partially due to the fact that we fit a Gaussian with positive amplitude, while we allow negative summed fluxes. Exclusion of the negative summed fluxes from the comparison improves the median agreements at low S/N (not shown).

Zoom In Zoom Out Reset image size

In light of this discussion, we conclude that emission lines are statistically well fit by a single Gaussian model at low S/N. The increasing ${\chi }^{2}/\mathrm{dof}$ does not imply that line fluxes from Gaussian fitting are unreliable in this regime, as can be demonstrated by a comparison with the nonparametric summed fluxes. Since S/N_line correlates very well with the empirically derived A/N, we suggest that, despite its simplicity, S/N_line is an excellent metric for the uncertainty in the fit. Henceforth, in this paper, S/N_line is used instead of A/N.

The typical line S/N of MaNGA data. To conclude this section, we present in Figure 5 the S/N distribution of the same sample of spaxels used in Figure 3, which is representative of the line S/N distribution in the MaNGA data. Only spaxels with S/N > 0 are plotted; that is, we do not plot the large number of spaxels that have no detected line emission according to our fitting procedure. Colored lines show the S/N distribution in three different radial bins. The radial variation of these S/N distributions highlights the decrease in S/N even in the strong nebular lines for $R\gt 1.5\ {R}_{e}$. We note that the MaNGA sample includes both star-forming and passive galaxies that are characterized by low-S/N line emission. In particular, the bimodality between star-forming and low-ionization emission-line regions (LIERs; Belfiore et al. 2016) is evident as a bimodality in the S/N of the Balmer lines, especially at small galactocentric radii.

Zoom In Zoom Out Reset image size

In order to test the the presence of possible systematic errors in the recovery of emission-line parameters and the statistical correctness of the errors produced by the DAP, we have carried out a set of idealized recovery simulations. Four test galaxies were selected to span a wide range of stellar-continuum and emission-line properties (two star-forming blue galaxies and two red LIER galaxies). Considering all four galaxies, our mock data set consists of ∼5000 spaxels with S/N > 1 in Hα.

The MaNGA datacube for each galaxy was fit using the DR15 version of the DAP, and the best-fit model cube, including both continuum and emission lines, was used as a template for generating “mock” datacubes. For each spaxel in the model cube, Gaussian noise was added to the model spectrum, with a standard deviation given by the error vector in the input MaNGA data. Assuming the MaNGA DRP errors are accurate, this procedure generates mock cubes with the same noise level as the original data. Mock cubes with twice and half the noise level of the original data were also created.

All of the mock cubes were run through the MaNGA DAP in the same way as real MaNGA data. In particular, the same MILES-HC stellar templates were used to generate and then fit the mock cubes. In Section 4.4 we repeat this exercise using a different template set to fit the simulated data, and we discuss the effect of template mismatch.

In the ideal case, the fits to the mock datacubes would recover the input values of flux, velocity, and velocity dispersion for all emission lines with no bias, and the (1σ) errors for these quantities reported by the DAP should be equal to the standard deviation of the residuals between the output and input values. In other words, we expect $\lt {q}_{\mathrm{in}}-{q}_{\mathrm{out}}\gt =0$ and $\mathrm{std}(({q}_{\mathrm{in}}-{q}_{\mathrm{out}})/{\mathrm{err}}_{\mathrm{out}})=1$, where angle brackets denote averaging, std is the standard deviation, ${q}_{\mathrm{in}}$ and q_out are input and output values for a physical quantity, and err_out is the DAP-provided error in q_out. In this section, we will refer to $({q}_{\mathrm{in}}-{q}_{\mathrm{out}})/{\mathrm{err}}_{\mathrm{out}}$ as the normalized residual for quantity q.

In Figure 6(a) we plot the offset between input and output fluxes (in dex) as a function of measured (output) S/N for six strong emission lines ([O ii]λ3727, Hβ, [O iii]λ5007, Hα, [N ii]λ6584, and [S ii]λ6731). The plot demonstrates the existence of a small positive bias in the recovered flux for S/N < 6, which then becomes a sizable decrease in the recovered flux for the lowest S/N levels (∼2–3). Overall, the ability of the code to recover the input fluxes is better than 0.05 dex (12%) for S/N  > 6.

Zoom In Zoom Out Reset image size

In Figures 6(b)–(d), we plot the normalized residuals for the flux, velocity, and velocity dispersion as a function of the output S/N. The solid colored lines correspond to the median values of the normalized residuals in logarithmic bins of S/N, while the dashed lines represent the 16th and 84th percentiles. In all panels, in the case of perfect recovery, the median lines would lie at zero, and, assuming Gaussian errors, the 16th and 84th percentiles would follow horizontal lines at ±1.

From Figure 6(b), we observe that fluxes of all of the emission lines considered are recovered with negligible bias down to S/N ∼ 1.5. More notably, the errors are also correctly estimated, because the 16th and 84th percentiles closely follow the ±1 lines in normalized residuals. For S/N < 1.5, the flux is systematically underestimated. At these low S/N, the distribution of normalized residuals also deviates from a Gaussian, showing a long tail at low normalized residuals. We note that, while different lines cover different ranges in S/N, the behaviors of different lines are remarkably similar.

In Figure 6(c) we show the normalized residuals for the emission-line velocities as a function of Hα S/N. Only the Hα line is plotted in this panel, because all of the emission lines are fit with the same velocity. The figure shows that the input emission-line velocities are recovered with no bias down to S/N ∼ 1.5. However, the formal errors calculated by the DAP are underestimated for S/N < 10. In particular, at S/N ∼ 2, the output error is a factor of ∼3 lower than expected. At high S/N, on the other hand, the output errors are consistent with the scatter in the normalized residuals. The source of this underestimation likely lies in the fact that the formal error provided by pPXF for the gas fluxes is computed, for computational efficiency, from the covariance matrix of the gas emission templates alone. This implies that the uncertainties currently ignore the covariance between the fluxes (which are linear parameters in the fit) and the gas kinematics (which are nonlinear parameters). Proper uncertainties could be computed via bootstrapping at the expense of a significantly larger computation time, or by recomputing the covariance matrix with respect to all variables at the best-fitting solution.

In order to quantify this deviation, we have fit the observed relation with a simple functional form:

$$\begin{eqnarray}\begin{array}{rcl}({v}_{\mathrm{in}}-{v}_{\mathrm{out}})/{\mathrm{err}}_{v\mathrm{out}} & = & (0.8\pm 0.1)+(0.49\pm 0.07)\\ & & \times \ \mathrm{log}{\left({\rm{S}}/{\rm{N}}\right)}^{-1}.\end{array}\end{eqnarray} \tag{ 3 }$$

The resulting fit provides a very good representation of the data and is shown in light red in Figure 6(b). This correction is not applied to DR15 DAP output and needs to be taken into account by the user. We anticipate that users interested in fitting detailed kinematic models to the emission-line velocity field may need to take this correction factor into account.

Figure 6(d) shows the normalized residuals for the velocity dispersions of the different emission lines considered. The velocity-dispersion trends are similar to those observed for the flux and are indeed their likely cause, because the flux is positively correlated with dispersion. We observe remarkably good agreement in both the median and 16th and 84th percentile values down to S/N ∼ 1.5. Below that value, the dispersion shows a larger tail of negative normalized residuals.

Overall, idealized recovery simulations with no template mismatch demonstrate that the values and errors of flux and velocity dispersion can be recovered accurately with negligible bias down to S/N ∼ 1.5. The velocities can also be recovered reliably down to low S/N, but their associated errors appear to be underestimated for S/N < 10.

In this section, we further analyze the error statistics for the emission-line measurements provided by the DAP by using repeat observations. In DR15, 56 galaxies have been observed more than once, mainly for the purpose of testing random and systematic errors (Westfall et al. 2019, Table 1).²¹ After processing these galaxies through the DAP, their MAPS files were transposed into the same world coordinate system, in order to account for small shifts in the integral-field unit (IFU) bundle positions and orientations between observations. The world coordinate system is derived by the MaNGA DRP by matching the MaNGA cubes to preexisting SDSS photometry in the advanced astrometry module (Law et al. 2016, Section 8), and it therefore takes into account small shifts and rotations of the IFU fiber bundles. In comparing repeat observations, we do not, however, take into account possible changes in the seeing conditions.

Similar to the procedure adopted to analyze the recovery simulations in the previous section, we calculate the normalized residuals as a function of S/N. Because for repeat observations we do not know the true value of any physical quantity, we define the normalized residual as $({q}_{1}-{q}_{2})/{\left({\mathrm{err}}_{1}^{2}+{\mathrm{err}}_{2}^{2}\right)}^{1/2}$, where 1 and 2 refer to a pair of repeat galaxies, and err is the error in quantity q. Considering all repeat galaxies, we obtain a sample of ∼5 × 10⁴ pairs of spaxels with two independent measurements with S/N > 1 for Hα.

In Figure 7(a) we show the normalized flux residual as a function of S/N of the first galaxy in the pair. Following the same graphical conventions as in Figure 6, the solid colored lines represent the median, while the dashed colored lines represent the 16th and 84th percentiles as a function of S/N. The median residual is found to be close to zero. The 16th and 84th percentiles, on the other hand, are found to be close to ±1 at S/N ∼ 2 and show a systematic deviation toward larger values at higher S/N. This deviation is particularly evident for S/N > 10.

Zoom In Zoom Out Reset image size

In Figures 7(b) and (c), we show the normalized residuals as a function of S/N for the velocity and velocity dispersion. For velocity, the errors are underestimated at all S/N, with the worst discrepancy at low S/N, similar to what was found in our study of idealized simulations in the previous section. Different from what was seen in the previous section, the errors also diverge from expectations at high S/N (S/N > 20–30), while they appear to be underestimated by a factor less than 2 in the range S/N = [6, 50]. The velocity dispersion appears to be much better behaved, with no evidence for large error underestimation until S/N > 100. Interestingly, for S/N = [1, 2], the errors appear to be overestimated.

Figure 7(d) shows the distribution of normalized residuals for six strong lines (see legend) in the S/N range [1–10]. The gray dashed line shows a normalized Gaussian of unit standard deviation, which represents the theoretical expectation in case of ideal error measurements. In Figure 7(d) we show as vertical black dashed lines the 16th and 84th percentiles of the observed distribution for Hα (1.22 and −1.35, respectively). We note that, although the data presents slightly non-Gaussian tails, it is well fit to first order by a Gaussian with standard deviation 1.25 (fit not shown).

Figure 7(e) is the same as Figure 7(d), but represents the S/N range [10–100]. As already evident in Figure 7(a), at this S/N level, either the errors are underestimated by a factor of 2–3, or some other systematic enters into the comparison of repeat observations.

Because this large error underestimation is not seen in the idealized recovery simulation, we consider possible systematic effects that could cause this. First, as already seen in Section 3.1, at S/N > 20–30, our Gaussian model may be insufficient to accurately fit the line profiles in real data, leading to higher normalized residuals and, possibly, underestimated errors. Second, regions of bright line emission tend to be clumpy, and the measured fluxes are therefore particularly sensitive to differences in point-spread function (PSF). In this case, the increased error in flux is due to intrinsic scatter in the amplitude, and not to larger errors in the recovered velocity dispersion, which would be in agreement with the findings from Figure 7(c). The same effect would be caused by small astrometric misalignments between repeat observations.

We gained some insight into these issues by visually inspecting difference and normalized residuals maps for different pairs of repeat observations. This exercise clearly revealed that the largest normalized residuals are indeed associated with bright and clumpy line emission. We have therefore performed a simple test to quantify the effect of astrometric offsets and PSF differences. One of the galaxies showing the largest differences in normalized residuals was selected as an example. For this galaxy, we considered the output of the MaNGA astrometry module, which matches the MaNGA IFU data to the underlying SDSS photometry in order to correct for small deviations of the rotation and centroid position of the MaNGA IFU ferrules in a given exposure due to the mechanical tolerance of the ferrule and rotational clocking pin holes. We artificially added random error to the best-fit astrometric solution, consistent with the uncertainty calculated by the astrometry module (typically about 025 and 01 for the rotational and translational components, respectively). A datacube was then produced, following the usual MaNGA reduction recipes, and fit using the DAP. We compared the output produced by this datacube with additional astrometric error to the reference datacube generated for DR15. At low S/N, the dispersion in line fluxes between the two datacubes is negligible, but it flares at high S/N in a fashion consistent with that observed in Figure 7(a). In particular, we find that astrometric errors consistent with those expected by registering MaNGA data to SDSS photometry are sufficient to explain the observed increase in the error budget in repeat observations.

In summary, recovery simulations demonstrate that the errors for flux and velocity dispersion of different emission lines behave in a statistically correct fashion down to S/N ∼ 1.5. Errors in the velocity are underestimated for S/N < 10, and the source of this discrepancy is not know at the time of writing. Equation (3) quantifies this underestimation and can be utilized to rescale the errors based on the outcome of the recovery simulation. In DR15 we leave it to the user to apply this correction if deemed necessary to their science goal.

These trends are largely confirmed by the analysis of repeat observations. Repeat observations, however, also show underestimation of the errors in the high-S/N regime. We have demonstrated that this trend can be entirely explained by small astrometric errors in individual exposures, which are consistent with the uncertainties derived by the MaNGA DRP astrometric registration routine. In light of this discussion, we leave it up to the user to consider whether adding this extra error contribution is advisable for their specific science goal.

In the following section, we briefly outline the characteristics of the spectral libraries that we will discuss and compare in this paper.

Hierarchically clustered MILES templates (MILES-HC). As discussed in Section 5 of Westfall et al. (2019), we have applied a hierarchical clustering algorithm to the MILES stellar library, which consists of 985 stars covering the wavelength range 3525–7500 Å (Sanchez-Blazquez et al. 2006; Falcón-Barroso et al. 2011). The clustering algorithm subdivides the stars in the MILES library into a number of groups that are defined to be maximally different from each other. Forty-nine such groups are generated, and a composite spectrum is obtained for each group as the average of the spectra of the contributing stars. These 49 spectra were visually inspected, and seven of them were removed because of artifacts or the presence of emission lines (in flaring late-type stars), leading to a total of 42 stellar templates. The resolution of the MILES library has been independently derived by Beifiori et al. (2011) and Falcón-Barroso et al. (2011) and is 2.54 Å (FWHM). This library is used in generating all of the DAP DR15 data products.

Maraston 2011 SSP models based on MILES (M11-MILES). These SSP are generated using the MILES stellar library by the Maraston stellar-population synthesis code (Maraston & Strömbäck 2011).²² A total of 110 models are used with ages ranging from 6.5 Myr to 15 Gyr, three metallicities (Z = 0.01, 0.02, 0.04), and a Salpeter IMF. The spectral resolution is the same as that of the MILES library (2.54 Å FWHM).

Vazdekis MIUSCAT SSP Models (MIUSCAT). This is a set of 72 SSP models generated according to Vazdekis & Ricciardelli (2012)²³ with a Salpeter IMF (unimodal IMF with slope = 1.3) and a set of 24 ages (0.0631, 0.0794, 0.1000, 0.1259, 0.1585, 0.1995, 0.2512, 0.3162, 0.3981, 0.5012, 0.6310, 0.7943, 1.0000, 1.2589, 1.5849, 1.9953, 2.5119, 3.1623, 3.9811, 5.0119, 6.3096, 7.9433, 10.0000, 12.5893 Gyr) and three metallicities ([M/H] = −0.4, 0.0, 0.22). MIUSCAT extends the wavelength range of MILES to cover the full range 3465–9469 Å by making use of the near-infrared CaT library of Cenarro et al. (2001). Stellar spectra from the Indo-US stellar library are used to fill in the gap left between the MILES and CaT spectral ranges and also to extend toward the blue and red the wavelength coverage of the MILES and CaT libraries, respectively. The spectral resolution is the same as the MILES library, as the higher-resolution CaT and Indo-US libraries are convolved to the MILES spectral resolution. The MaNGA VAC generated by the Pipe3D team also uses a subset of MIUSCAT templates for spectral fitting, although the exact set of templates differs from the ones described above.

Bruzual and Charlot SSP models based on STELIB (BC03). This is the set of 40 Bruzual & Charlot (2003) models used in the MPA-JHU catalog for the SDSSDR4.²⁴ The SSP templates cover a range of 10 ages (0.005, 0.025, 0.101, 0.286, 0.640, 0.904, 1.434, 2.500, 5.000, 10.000 Gyr) and four metallicities (Z = 0.008, 0.004, 0.02, 0.05) using a Chabrier IMF. The SSP models are based on the STELIB stellar library (Le Borgne et al. 2003) and have nominal spectral resolution of 2.3 Å. These SSP models are different from the BC03 SSP based on MILES used to produce the MPA-JHU value-added catalog for later data releases, such as the latest DR8 version.²⁵

In this section, we investigate the consistency between the MILES-HC stellar library used in DR15 and the SSP models described in the subsection above. SSP models prescribe parameterized combinations or interpolations between elements of stellar libraries based on physical models of stellar interiors, atmospheres, and their evolution (isochrones). To the extent that these physical models are incomplete or simply inaccurate, one may worry about whether SSP and stellar spectra are consistent in both overall shape and in the details of their absorption features. SSP models based on observed stellar libraries mitigate, for example, incomplete line lists in synthetic spectra, but suffer problems of library incompleteness and stellar misclassification. Library incompleteness is a general concern for all continuum-fitting methods, and it is the specific concern motivating the work in this section.

In particular, the lack of hot (O-type) stars in the MILES library potentially undermines our ability to fit young stellar populations with MILES-HC. The hierarchical clustering performed on the spectra may further dilute the blue continua of the few B stars present in MILES. While O-star spectra are largely featureless, their very blue continua cannot be reproduced by a linear sum of other stellar spectra. A common solution to this problem, within the pPXF framework, is to allow the inclusion of additive or multiplicative polynomials in the fit. These polynomials are generally of low order to avoid any effect of the polynomial on the spectral features of individual absorption lines.

However, hot stars do have some critical features, notably in their hydrogen and helium lines, that have distinct and systematic changes with temperature. These changes in stellar absorption lines include the EW, the core width (due to rotation and winds), and the relative strength and shape of the wings (due to the Stark effect). None of these are modified by multiplicative polynomials, and only the EW (not the shape) can be modulated by additive polynomials. Further, because of the different and nonlinear temperature dependence of these features, it is not possible to accurately simulate spectra containing O stars with a library that does not include these stars. For a library with such a deficiency, we would expect to see Balmer lines that are too narrow and a deficiency of He absorption. Because this will lead to systematics in the continuum model at the location of key emission lines, the amplitude of such systematics is important to assess.

We therefore test the ability of the MILES-HC library to reproduce different stellar populations by fitting the BC03 SSP models of different ages (and solar metallicity) with MILES-HC spectra. The fit has been performed in three ways: (1) with no polynomials, (2) with an eighth-order additive polynomial, and (3) with an eighth-order multiplicative polynomial. Aside from the polynomial type, we perform this fit in the same way as the first (i.e., the stellar kinematics) fitting stage of the DAP. In each case, we compute the residual between the best-fit model and the input SSP and compute the resulting rms over the wavelength range 3700–7400 Å.

The fractional rms (i.e., the standard deviation of output–input/input) calculated over ∼40 Å windows is shown in Figure 8 for BC03 templates of different ages fitted with no polynomials.

Zoom In Zoom Out Reset image size

The figure demonstrates that the overall spectral shape of the BC03 SSPs are well fit by MILES-HC, with median residual rms values of 10^−2.0. The largest rms values are seen both at the blue end of the spectrum at very young ages and in very localized wavelength regions, generally corresponding to notable absorption lines. Balmer series lines are particularly problematic and increasingly so at young ages. Metal lines (such as Mgb and NaD), on the other hand, are fit worse at older ages. The inclusion of polynomials does not lead to an overall improvement of the fit quality, although it does have an effect on the fit around the Balmer and helium lines.

Figure 9 compares the fits of the MILES-HC library to the BC03 SSPs for ages of 5 and 25 Myr with and without polynomials. This figure is worth careful scrutiny. Inspection reveals MILES-HC underpredicts the hydrogen line depths, increasingly for lower-order lines. This is mostly ameliorated by either additive or multiplicative polynomials, which do a good job of matching the wings but fail to match the core. The same relative statements are true for the 25 Myr SSP, but the amplitudes of the differences are decreased, that is, the MILES-HC fit is substantially better on its own without polynomials. In contrast, the MILES-HC fits overpredict He i and underpredict He ii lines for BC03 for both ages. Polynomials do little to help remedy the mismatch in EW and often degrade the quality of the fit.

Zoom In Zoom Out Reset image size

We can interpret the overpredicting of the He i and the underpredicting of He ii lines as being due to the lack of very hot O stars in MILES-HC. B stars are the hottest stars in the library, and they do not have He ii lines. If MILES-HC lacks templates with strong He ii lines, then the fit will use more B stars to compensate for the spectral shape and end up overfitting He i, while still not producing the He ii features.

A detailed accounting of the stellar templates (and their weights) that go into the specific SSPs would be one way to make progress on this question, but this critical “deconstruction” of SSPs is beyond the scope of this paper.

We conclude that even with additive or multiplicative polynomials, MILES-HC is likely to have small systematic residuals in the cores of the hydrogen lines that lead to emission-line overestimates at the very youngest ages. The systematics for helium lines are more significant and varied, and in some cases are minimized without including polynomials. Polynomials do not significantly improve the overall match of the spectral between SSPs and stellar templates, which is excellent except for the youngest ages for λ < 4000 Å.

In order to explicitly test the effect on measured emission-line fluxes caused by the use of different stellar-continuum models, we fit a subset of DR15 MaNGA datacubes with the three sets of SSPs discussed in Section 4.1, in addition to the MILES-HC fit performed in DR15. In particular, we used SSPs to perform the second fitting stage in the DAP, but not for the extraction of the stellar kinematics, for which stars are generally recommended. When necessary, the difference between the intrinsic spectral resolution of the MILES stars and that of the SSP templates used for the second fitting stage has been taken into account.

The SSP fits were carried out for a sample of 15 galaxies, evenly sampling the NUV − r versus $\mathrm{log}({M}_{\star }/{M}_{\odot })$ plane, in order to have access to a wide variety of stellar populations. We only considered galaxies with line emission (including extended LIER galaxies on the red sequence). Considering the entire galaxy subsample, we obtain ∼2 × 10⁴ spaxels with Hα S/N >1.

In Figure 10 we compare the emission-line fluxes obtained using MILES-HC for both fitting stages (i.e., the DR15 data products) with the fluxes obtained after switching to an SSP template for the second fitting stage. The flux ratios are presented as a function of line S/N for different strong lines. The dashed gray lines represent the level at which the flux difference is comparable to the random error.

Zoom In Zoom Out Reset image size

Figure 10 shows a number of interesting features. For the Balmer lines (Hα and Hβ are shown in the figure), different templates give systematically different line fluxes because of the different best-fit stellar Balmer absorptions, especially at low S/N. MIUSCAT prefers deeper Balmer absorption, leading to larger Balmer line fluxes than MILES-HC. BC03 and M11-MILES lead to better agreement with DR15 for Hα but show significant differences in Hβ. It is interesting to note that the systematic discrepancies are substantial (up to 0.1 dex at S/N = 2 between DR15 and MIUSCAT). They are, however, smaller than the random errors at low S/N, while they become comparable to (or larger than) the random error at high S/N. We have checked the behavior of the flux differences as a function of EW of the lines, and the resulting plot is very similar to Figure 10, especially for the Balmer lines, because S/N largely tracks the EW.

On the other hand, for metal lines, such as [O iii]λ5007, [S ii]λ6732, and [O ii]λλ 3727,29, the discrepancies between fluxes obtained with different templates are less extreme and do not correlate as well with line S/N. [N ii]λ6585 and [S ii]λ6732 stand out from the other metal lines for showing comparatively larger discrepancies. In Figure 11 we show some example fits to the spectral regions around Hα and the [N ii] doublet for the central spaxel of four galaxies (low-mass blue, 7815–6101; massive blue, 8138–12704; low-mass red, 8329–1901; high-mass red, 8258–6102) spanning a range of properties in the NUV − r versus $\mathrm{log}({M}_{\star }/{M}_{\odot })$ plane. These example fits highlight the previously discussed differences in the core of the Hα line, but also the resulting effect on the nearby [N ii] lines, which are the outer edges of the Balmer absorption wings. For example, the BC03 templates generate best-fit models that have substantially different line wings from those of other template sets, therefore affecting the [N ii] flux in addition to Hα.

Zoom In Zoom Out Reset image size

Although these flux discrepancies are smaller than the random error, they are systematic and behave differently for the different lines considered, therefore leading to biases in the derived line ratios. In Figure 12 we show the differences in dex for several line ratios and other derived quantities between the cases fit with SSPs and DR15 as a function of EW(Hα). In the first row, we show the Balmer decrement (Hα/Hβ) and two classical BPT (Baldwin–Phillips–Terlevich; Baldwin et al. 1981) line ratios ([N ii]/Hα and [O iii]/Hβ). At low S/N, the Balmer decrement measured with MIUSCAT and M11-MILES differs substantially from that inferred in DR15 or using BC03. Estimating E(B − V) using an intrinsic ratio ${\rm{H}}\alpha /{\rm{H}}\beta$ = 2.86 and a Calzetti (2001) extinction curve, we see that deviations up to 0.1 dex are evident at EW(Hα) ∼ 2 Å, growing worse at even lower EW. Regarding the BPT line ratios, [N ii]/Hα is relatively unaffected by template choice, possibly because the vicinity of the two lines means they are affected by the best-fit continuum shape in a correlated way. [O iii]/Hβ, on the other hand, displays significant differences for low-EW lines (there is a 0.2 dex difference between MIUSCAT and DR15 at EW(Hα) ∼ 2 Å). These biases will have a measurable effect on the BPT diagram positions of low-EW regions, which tend to be associated with LIER emission and diffuse ionized gas.

Zoom In Zoom Out Reset image size

In the bottom row of Figure 12, we also show two metallicity-sensitive indices often employed in the literature, O3N2 = ([O iii]λ5007/Hβ)/([N ii]λ6583/Hα) (Pettini & Pagel 2004) and R23 = ([O ii]λλ 3727,29 + [O iii]λλ 4959,5007)/Hβ (Pagel et al. 1979). While O3N2 is relatively insensitive to dust extinction, we correct the measured line fluxes for extinction when computing R23. It is evident from the figure that discrepancies larger than 1/10th of a dex are present at low EW for both indicators. A cut on EW(Hα) > 6 is sometimes performed in studies of ISM metallicity in order to minimize the contamination from gas not directly associated with H ii regions (Sánchez et al. 2014). Here we show this threshold as a dashed black line for these two indicators in Figure 12, demonstrating for larger EWs that the systematic effects from the choice of continuum templates are nonnegligible. This exercise demonstrates that care is needed when comparing results from IFU surveys calculating emission-line fluxes with different underlying stellar or SSP models.

In the Appendix we perform a similar comparison on the line fluxes measured by the DAP and the Pipe3D VAC for DR15. The DAP and Pipe3D differ in many fundamental aspects besides the choice of continuum templates, so it is more difficult to attribute discrepancies to just one factor. However, for some line ratios (like Hα/Hβ and [O iii]/Hβ), we find discrepancies between the results of the two pipelines that are comparable, at least qualitatively, with those between the DAP DR15 run and the DAP run utilizing MIUSCAT templates. These differences may therefore be attributed at least partly to the different choice of continuum templates in the two pipelines.

However, while the results in this section show no strong effect of template choices on the [N ii]/Hα, we do find significant discrepancies for this line ratio between the DAP and Pipe3D. This fact is discussed further in Appendix A.2.

The aim of this section is not to select the “correct” set of templates, but simply to quantify the effect of different templates on the resulting emission-line fluxes. In general, we cannot determine which set of templates is the correct one for our galaxy data, so we devise an artificial exercise similar to the recovery simulation presented in Section 3.2 to study the effect of using the “wrong templates” in the presence of noise.

In particular, we take the best-fit model from the previous section based on a set of SSP templates and add noise in the same way as was done in Section 3.2. Here we discuss the results of using mock datacubes generated using the MIUSCAT best-fit models, which are then fit using the standard DR15 approach (i.e., using the MILES-HC library). In light of the results of the previous section, we expect the recovered Balmer line fluxes to be lower than the input ones on average, given the preference for MIUSCAT to fit deeper Balmer absorption features.

Figure 13 demonstrates this effect. The recovered Hα and Hβ fluxes are indeed systematically lower than the input ones. It is interesting to note, however, that while the median flux is systematically biased, the 16th and 84th percentiles still approximately correspond to ±1 with respect to the median, indicating that template mismatch does not dramatically affect the statistical validity of the emission-line flux errors. At low S/N, we also observe a bias in the recovered fluxes, in the sense that the flux tends to be underestimated, as already noted in Section 3.2 (see Figure 6(a)).

Zoom In Zoom Out Reset image size

In this section, we assess the role and importance of multiplicative polynomials in the second fitting stage of the DAP (the simultaneous fitting of continuum and emission lines). A complementary discussion of the role of additive polynomials during the first fitting stage of the DAP (stellar kinematics) is presented in Section 7.3.3 of Westfall et al. (2019). We remind the reader that the multiplicative polynomials are only applied to the stellar-continuum templates and not to the emission-line (Gaussian) templates. In this sense, their effect on the line fluxes is only indirect. We nonetheless address this issue here as a check on the quality of our continuum model and for its relation to the overall flux calibration of the MaNGA survey.

5.1.1. Role of Polynomials

5.1.2. Typical Shape of the Polynomial Correction in DR15

In this section, therefore, we start by looking at the typical shapes of the multiplicative polynomials used in the second fitting stage for the DR15 DAP run. To do so, we selected a random sample of 100 DR15 galaxies and reconstructed the multiplicative polynomials used in each of their spaxels. In Figure 14 we show the shape of the median multiplicative correction applied as a function of rest-frame wavelength. The sample of galaxies is subdivided into two mass bins (red for $\mathrm{log}({M}_{\star }/{M}_{\odot })\gt 10;$ blue for $\mathrm{log}({M}_{\star }/{M}_{\odot })\lt 10$) and the shaded areas correspond to the 16th and 84th percentiles of the distribution. We also show in black the expected multiplicative correction for a Calzetti extinction curve and two values of E(B − V).²⁶ The extinction curves are scaled arbitrarily to the median of the polynomial corrections to highlight the similarity in relative shape.

Zoom In Zoom Out Reset image size

The trends observed in the figure can be qualitatively interpreted as follows. The massive bin contains a larger number of passive galaxies, which are largely devoid of gas and thus suffer lower extinction. The shapes of the polynomials are consistent with the values of E(B − V) measured for the continuum by full spectral fitting in the outskirts of MaNGA galaxies (Goddard et al. 2017).

5.1.3. Deviations from Smooth Polynomial Shapes and Consequences for Flux Calibration

The key features in Figure 14 are the deviations from the expected smooth extinction curves, namely the upturn in the mean correction at the blue end of the MaNGA wavelength range and a similar downturn redder than 7000 Å. We determined that these deviations are not due to imperfections in the MaNGA flux calibration for the following reasons.

First, the upturn in the blue occurs at different observed wavelengths for galaxies at different redshifts. For example, if one considers massive galaxies in the MaNGA primary and secondary samples, which are selected in the same fashion but separated by a small redshift interval, the blue upturn moves to longer observed wavelength in the secondary sample. If the upturn was due to imperfections in the MaNGA flux calibration derived from standard-star spectra, it would always appear at the same observed wavelength. Second, the downturn observed redder than 7000 Å occurs in the middle of the MaNGA spectral coverage (but at the edge of the spectral coverage of MILES-HC) and is therefore more likely to be originating from the MILES-HC than the MaNGA data.

In order to test whether imperfect relative flux calibration of the MILES-HC library is responsible for the deviations observed in Figure 14, we selected one test galaxy (a massive red galaxy, 8258–6102, with good S/N throughout) and examined the stacked polynomial shapes obtained after fitting the galaxy with different template sets. The results are presented in Figure 15. It is interesting to note that between 4000 and 7000 Å, MILES-HC (labeled DR15), M11-MILES, MIUSCAT (all based on MILES stars), and BC03 (based on STELIB) agree to better than 10%. Bluer than 4000 Å, BC03 presents a downturn, while both M11-MILES and MILES-HC show an upturn. MIUSCAT, on the other hand, gives rise to a flattening. It should also be noted that the downturn at 7000 Å is present both in MILES-HC and M11-MILES, pointing toward a problem with the MILES stars.

Zoom In Zoom Out Reset image size

To check the behavior of the code at the edges of the wavelength range, we changed the degree of multiplicative polynomials from eight to seven (i.e., from even to odd parity). If the behavior of the polynomials at the edges was entirely dictated by the fit within the central wavelength region, we would expect that a change of parity would lead to a change in symmetry of the recovered multiplicative correction, which is however not seen in Figure 14. We concluded, therefore, that the offsets seen at the edges of the fitted MaNGA wavelength range are likely to be real.

It is possible that MILES-HC suffers from the lack of hot stars, including blue horizontal branch stars. These stellar types may have been accounted for differently by different SSP models, generating the discrepancy observed between BC03, MIUSCAT, and M11-MILES at the blue edge of the optical wavelength range.

Interestingly, the SSP templates that extend redder than 9000 Å show the need for an upward correction to match the MaNGA data. The presence of this red upturn has been identified via visual inspection in some of the MaNGA spectra. Since this spectral range is not fit in DR15, we postpone further study of this potential systematic effect.

The MaNGA DAP implements simultaneous fitting of emission-line and continuum templates following the recommendation from previous work (Sarzi et al. 2005; Oh et al. 2011). Sarzi et al. (2005), in particular, demonstrated the advantage of this algorithmic choice when dealing with the limited wavelength range of the SAURON data, where the emission lines lie in close vicinity to the key metallicity- and age-sensitive features. In Section 4.3, however, we have demonstrated that, even by performing simultaneous fitting, residual degeneracies between Balmer absorption and line emission are still present, leading to noticeably different best-fit models when using different template libraries. In this section, therefore, we perform some illustrative tests to evaluate the impact of the masking on the recovered best-fit continuum under the Balmer lines.

We first examine the difference between the best-fit continua obtained by the first and second fitting stages in the DAP for the central spaxels in four test galaxies (Figure 16, same galaxies as in Figure 11). A difference between the two best-fit stellar-continuum models in this comparison may be due to either the effect of masking (emission-line regions are masked in the first fit but not in the second) or the difference in the use of polynomials (additive polynomials in the first fit and multiplicative polynomials in the second fit). In Figure 16 we plot the difference between the two best-fit models normalized by the input spectrum. At the wavelength of a specific emission line, this can be interpreted as the fractional error in the amplitude of the line introduced by these different choices in continuum fitting.

Zoom In Zoom Out Reset image size

The largest deviations are seen in regions corresponding to strong absorption lines (like the NaD doublet, evident in all four examples except the low-mass blue galaxy) and the Balmer absorption lines. However, in the case of Balmer lines, the differences can be both positive and negative. We attribute this behavior to the different implementation of polynomials in the two fitting stages. In the case of the low-mass blue galaxy, where the most prominent absorption features are the Balmer lines, additive polynomials lead to shallower absorption-line profiles, which result in positive residuals at the positions of the Balmer lines in Figure 16. For the higher-mass galaxies, other metal absorption lines dominate the spectrum and therefore determine the shape of the additive polynomials, causing positive, null, or negative residuals around the Balmer lines. We note that the deviations observed in the low-mass blue galaxy are much larger than those observed in the red galaxies, with significant changes already seen at Hβ (∼3%) and increasing to 13% at H.

A cleaner test to isolate the effect of masking is to apply the same type of polynomials to both the first and second fitting stages. We have therefore repeated the exercise just described by using eighth-order multiplicative polynomials for both fitting stages. The resulting normalized flux differences are shown in Figure 17 and look substantially different from Figure 16. Now the differences around metal absorption lines are reduced, and the Balmer lines correspond to the largest residuals (of the order of ∼1%–2%). These are again stronger for the high-order Balmer lines (in particular Hδ and Hγ). Integrating over the line profile, we see the systematic differences in Hα flux are less than 2%, which is negligible in most cases when compared to the discrepancies caused by changes in the template library. However, the changes in measured fluxes are more substantial for the high-order lines.

Zoom In Zoom Out Reset image size

There is an interesting difference between the young spectrum of the low-mass blue galaxy, which displays deeper Balmer absorption in the second (unmasked) fit, and the other spectra characterized by older stellar populations, which show shallower absorption in the second fit. The reasons for this difference must be related to how the inclusion of the masked regions affects the best-fit template mix. In the future, it would be of interest to repeat the same exercise in the context of stellar-population synthesis and assess the effect of masking on the recovery of stellar-population parameters, which is likely more significant than the effect on the emission lines.

When several emission lines are fit across a large wavelength range, whether or not to tie the kinematic parameters for different lines becomes a debatable problem. In general, tying the velocities and velocity dispersions of different lines is an advantage in the low-S/N regime, where stronger lines contribute much more to the overall χ², therefore effectively determining the kinematic parameters of weaker ones. Several reason exist, however, to be skeptical of tying kinematic parameters. First, without an accurate knowledge of the LSF and its change with wavelength, it is not possible to correctly fix the astrophysical velocity dispersions of widely separated lines. Likewise, small errors in the wavelength calibration can induce problems when fitting all emission lines with a common velocity. Finally, there are astrophysical reasons to expect emission lines emitted by different ionic species in different ionization stages to have different kinematics.

In this section, we test the effect of making different assumptions regarding the tying of kinematic parameters. We considered the sample of 15 galaxies described in Section 4.3, selected to evenly sample the NUV − r versus $\mathrm{log}({M}_{\star }/{M}_{\odot })$ plane. We consider the schemes described below.

All parameters free (all free). In this scheme, all velocities, dispersions, and amplitudes of the different emission lines are fit individually as free parameters.

Tie velocities (DR15). In this scheme, the velocities of all lines are tied together, while the velocity dispersions are fit independently. This tying scheme may be beneficial when uncertainties in the LSF prevent the tying of the astrophysical dispersions and is the scheme adopted in the DR15 run.

Tie velocities in three groups (v groups). In this scheme, we define three groups of emission lines:

1. 1.  
   Balmer lines: Hα, Hβ, Hγ, Hδ, H, Hζ, Hη, Hθ.
2. 2.  
   Low-ionization lines: [O ii]λλ 3727,29, [O i]λλ 6300,64, [N ii]λλ 6548,84, [S ii]λλ 6717,31.
3. 3.  
   High-ionization lines: [Ne iii]λλ 3869,3968, He iiλ4687, [O iii]λλ 4959,5007, He iλ5876.

The velocities of different lines are tied within the same group. None of the velocity dispersions are tied. This scheme is a variant of the DR15 run that allows for different astrophysical velocities for lines of different species and ionization states.

Tie velocities and dispersions in three groups (v + σ groups). The scheme is the same as the previous one, but we tie both velocity and velocity dispersions for the lines in the same group.

All parameters tied (all tied). In this scheme we tie the velocity and velocity dispersions of all lines together.

In all schemes, except the first one (all free), we impose the flux ratios set by atomic physics when fitting the line doublets of [O i]λλ 6300,64, [N ii]λλ 6548,84, and [O iii]λλ 4959,5007; see Table 1.

In Figure 18 we show the flux ratio (in dex) between each of the test runs and DR15, as a function of line S/N for different emission lines. The thick gray dashed lines correspond to the deviations expected given the random errors in the flux measurements.

Zoom In Zoom Out Reset image size

For S/N < 2, the all-free case gives larger fluxes than DR15 for all lines considered. Figure 19 shows a similar plot but for the velocity dispersions, demonstrating that the deviations toward larger fluxes are accompanied by lower sigmas. This may be due to the fact that, when velocities are not tied, the algorithm may be fitting noise spikes at low S/N. These spikes tend to have a width of one pixel, leading to smaller dispersions and higher amplitudes.

Zoom In Zoom Out Reset image size

The “v groups” case, where velocities are allowed to vary within emission line groups, is indistinguishable from DR15 in terms of fluxes and sigmas.

Finally, the “v + σ groups” and “all tied” test runs produce marginally lower sigmas and fluxes than DR15.

In Figure 20 we show the velocity difference between DR15 and different test runs for Hα [N ii]λλ 6584 and [O iii]λλ 5007. It is evident that there is no systematic velocity shift if velocities are tied in groups rather than all together. If the velocities are not tied (or tied in groups), the scatter in the resulting velocities with respect to DR15 is comparable to the random error (estimated in DR15) for the metal (both high- and low-ionization) lines, but smaller for the Balmer lines. Overall, the consistency of the velocities determined by following different tying prescriptions validates the DR15 approach and demonstrates that the MaNGA data do not suffer from any detectable systematic in the wavelength calibration (see also Figure 19 of Law et al. 2016), which would necessarily invalidate some of our tying schemes.

Zoom In Zoom Out Reset image size

Finally, in Figure 21 we show the ratio (in dex) between the Hα velocity dispersion and that of other Balmer (in black) and metal (in blue) emission lines as a function of wavelength, after subtracting in quadrature the DR15 estimate of the instrumental velocity dispersion. Only spaxels with S/N > 10 on the specific line are considered. The figure demonstrates that in DR15 the velocity dispersions of different emission lines within the MILES wavelength range are in good agreement with each other, with no significant wavelength-dependent systematic. The largest discrepancies are found in the blue end of the wavelength range, where both [O ii] and Hδ are larger than Hα by 0.07 dex (∼17%) on average.

Zoom In Zoom Out Reset image size

Here we briefly summarize our recommendations for usage of the DAP DR15 output regarding emission lines.

Users whose science goal would benefit from the best spatial resolution afforded by MaNGA should use the output from the hybrid binning scheme contained in the HYB10-GAU-MILESHC directory. If, on the other hand, one requires the emission-line properties to be computed on the same (Voronoi) binning scheme as the continuum, the output in VOR10-GAU-MILESHC should be used instead.

All of the emission-line maps are included as extensions in the DAP MAPS files. Each extension corresponds to a 3D array, where two dimensions correspond to the on-sky spatial pixels and the third dimension allows the user to choose a specific emission line. These extensions can be opened with a standard fits viewer (e.g., DS9,²⁷ QFitsView²⁸ ) and appear as a datacube. The correspondence between the index in the MAPS file extension and the line name is given in the extension header. A practical example of how to perform this association in an automatic fashion in python is given below.

We strongly encourage users to take close consideration of the masks provided in the MAPS file. In the python programming language, masked arrays, as implemented in the numpy package, are particularly suitable for the task of manipulating data with associated mask information. Masked arrays allow one to perform arithmetic and other simple operations (like taking a median) while automatically ignoring the masked pixels. The code example below provides an example of how to obtain the Hα flux map from the DAP MAPS file and encode it into a numpy masked array.

Users interested in a comprehensive software framework in the python language to access and manipulate the MaNGA data are encouraged to use our purpose-built Marvin package (Cherinka et al. 2019). Visualization tools for the DAP data products are also directly available on the Marvin web interface at https://dr15.sdss.org/marvin/.

Finally, we warn users of the following bugs that were discovered in the DR15 DAP output. These problems have already been resolved at the software level, and the DAP output will be corrected in the next data release.

1. 1.  
   The [O ii]λλ 3727,29 velocity-dispersion masks and errors are incorrect in DR15. In particular, the vast majority of spaxels in the velocity-dispersion maps of [O ii] are masked in DR15 because of a bug. The maps of the velocity dispersions themselves are, however, correct. We therefore recommend ignoring the mask extensions for the [O ii]λλ 3727,29 line velocity dispersions and applying the Hα velocity dispersion masks instead. Unfortunately, the same bug caused the inverse variance extensions for the [O ii] lines to be filled with zeros. This has been fixed for future releases, but for DR15, uncertainties for the [O ii] velocity dispersions are not available.
2. 2.  
   The Hζ line parameters are unreliable in DR15. This is due to a blend with the nearby He i line at 3889.749 Å (vacuum) that was not included in the line list for DR15, but will be included in future releases. We recommend against use of this line in DR15.

There are many ways to improve on DR15 with respect to the emission-line properties. First, we would like to transition away from MILES-HC toward a stellar library derived from MaStar spectra (Yan et al. 2018). MaStar is a new library of stellar spectra observed with the MaNGA instrument suite at APO and offers several advantages over MILES, in terms of its carefully controlled flux calibration, wider coverage in terms of stellar parameters, wider wavelength range, and LSF similar to the MaNGA galaxy spectra.

The generation of a hierarchically clustered set of MaStar spectra for kinematics extraction and the production of a new generation of SSP templates based on MaStar spectra are currently being pursued by the MaNGA team. We anticipate incorporating a “MaStar-HC” library or a set of MaStar-based SSP templates for the simultaneous fitting of continuum and emission lines in the second fitting stage of the DAP. The use of continuum templates with wider wavelength coverage would also allow us to provide accurate modeling of the continuum for the [S iii]λ8831,9071 9533 lines in the near-infrared, which are sensitive tracers of the ionization parameter of the ISM (Kewley & Dopita 2002).

A detailed characterization of the MaNGA LSF is currently underway and will be described in a forthcoming publication. Nonetheless, the results presented in Section 5.3, and Figure 21 in particular, demonstrate that the current estimate of the MaNGA LSF is sufficiently accurate to move toward a tying scheme where both velocity and velocity dispersion are tied for physically motivated groups of lines. We anticipate that the next MaNGA data release will adopt either the “all tied” or the “v + σ groups” tying approach described in Section 5.3.

Further improvements may include a treatment of multi-Gaussian kinematic components and broad emission lines in AGNs. While automatically detecting very broad emission lines (σ > 1000 km s⁻¹) is a relatively simple task, characterizing the significance of line asymmetries or additional kinematic components for millions of low-S/N, medium-resolution MaNGA spaxels represents a significant task, likely to remain outside the scope of a general-purpose data analysis pipeline like the DAP. Gallagher et al. (2019) describe the first attempt at this type of analysis on the MaNGA data.

Given the timescale of the SDSS-IV project, which is scheduled to terminate data collection in 2020, we do not expect further substantial project-led developments of the MaNGA DAP. Members of the astronomical community interested in adapting the MaNGA DAP to their own specific data format and scientific interests are welcome to make use of the DAP source code, which is publicly released on GitHub.²⁹

We first compare the fluxes between Pipe3D and the DAP for a set of strong lines ([O ii]λλ 3727,3729, Hβ, [O iii]λ 5007, Hα, [N ii]λ 6584, and [S ii]λ 6731). We perform this comparison using a large sample of 3  ×  10⁶ spaxels from 4565 datacubes. We compare the ratio of the fluxes from each suite pipeline as a function of the DAP S/N in Figure 22. At high S/N, the line fluxes agree to within 0.02 dex, or ∼5%. For S/N < 10, the DAP consistently estimates emission lines to be brighter than Pipe3D. For Hα the median of this effect is of order 0.05 dex (12%) for 3 < S/N < 10 and as much as 0.12 dex (28%) for [N ii]λ 6584 in the same S/N range. Hβ behaves differently from the other lines tested, as the DAP reports lower average fluxes compared to Pipe3D at low S/N. For the Balmer lines, this may be related to the differences in the stellar-absorption-line fits between the two pipelines; however, this explanation seems less likely for the forbidden lines. We note that for weak lines, the DAP stipulation that emission-line fluxes must be positive will introduce a positive bias, which may partially explain the upward skew seen in the distributions for low-S/N lines in Figure 22.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The small systematic offsets observed in the line fluxes can be enhanced when one computes some commonly used line ratios. Given the importance of establishing the consistency of these higher-level measurements, we compare the values of $E(B-V)$, [N ii]/Hα, and [O iii]/Hβ obtained by DAP and Pipe3D as a function of EW(Hα) (see Figure 23).

For EW(Hα) > 20 Å, the line ratios considered are in reasonable agreement between the two pipelines, although both for $E(B-V)$ and [N ii]/Hα a roughly constant systematic offset (∼0.02 dex) is observed even at high EW. The agreement considerably worsens at low EW. Two main differences exist between the methodologies of Pipe3D and the DAP that might influence the fluxes of the emission lines. These two pipelines use different template libraries to fit the stellar continua. This difference can result in differences in line fluxes, particularly at low EW. Since the [N ii]λ 6584 line is adjacent to ${\rm{H}}\alpha$, the features in the wings of the absorption line can impact the measured flux of the [N ii]λ 6584 line, particularly in spectra with high stellar velocity dispersion.

The reasons for this behavior are explored in the left-hand panel of Figure 24, where we show four example spectra and the stellar-continuum fits performed by the DAP and Pipe3D. These spectra were chosen to cover a range of [N ii] EW and σ_⋆ values. In these examples, we can see that differences in the adopted continuum model can have an effect on the flux in the emission lines, even in the spectral region around the forbidden lines. For the highest stellar velocity dispersions, the Hα absorption can spread underneath the [N ii] lines, leading to an overestimate of the [N ii] flux.

Zoom In Zoom Out Reset image size

The potential effect on a specific scientific result of the discrepancy between the DAP and Pipe3D fluxes is illustrated in Figure 24, right panel, where we show the resolved mass–metallicity relation obtained from the two pipelines using the O3N2 metallicity calibrator and the Pettini & Pagel (2004) calibration. We consider only spaxels classified as star-forming using the Kewley et al. (2001) line in the [S ii] BPT diagram. All spaxels with S/N < 3 for the emission lines required for the O3N2 diagnostic are excluded. The stellar mass surface density values for each spaxel are taken from the Pipe3D VAC. These values are corrected for the effect of dust extinction, but not for the potential effect of galaxy inclination. The shape of the resulting resolved mass–metallicity relation agrees well with previous studies based on a smaller sample of MaNGA galaxies (Barrera-Ballesteros et al. 2016). We note, moreover, that the determinations obtained from both pipelines are in very good agreement with regards to the shape of the relation, despite a systematic shift in the median metallicity. This shift is of the order of 0.03 dex at low metallicity, decreasing to 0.01 dex at high metallicity. Since the systematics associated with metallicity calibrations are larger than this offset, we consider that the choice of pipeline does not significantly affect this science case.

The comparison performed in this appendix highlights the difficulties in measuring the strengths of emission lines in galaxy spectra. Systematic differences in the derived fluxes for strong, high-EW lines are within a few percent for the majority of the lines examined here. We have checked that the effect on the resolved mass–metallicity relation is small (0.03 dex at most) but clearly systematic. However, caution should be taken not to overinterpret fluxes and flux ratios of emission lines, particularly in the regime of low line EW and for weak lines that are likely to be severely affected by the quality of the continuum subtraction.