The MOSDEF Survey: Broad Emission Lines at z = 1.4–3.8^*

In this work, we use the first 2 years of data from the MOSDEF survey (Kriek et al. 2015), where we obtained near-infrared spectra for ∼1500 high-redshift galaxies using the MOSFIRE instrument (McLean et al. 2012) on the W. M. Keck telescope. These spectra were collected over the course of 48.5 nights from 2012 to 2016, and enable measurements of the rest-frame optical nebular emission lines for galaxies at 1.37 < z < 3.80. The MOSDEF survey targets galaxies in the AEGIS, COSMOS, GOODS-N, GOODS-S, and UDS extragalactic legacy fields, which have extensive ancillary data, including Chandra, Spitzer, Herschel, HST, VLA, and ground-based optical/near-IR data.

One-dimensional spectra were extracted using custom IDL software called BMEP,⁹ as described in the appendix. BMEP was tested with output from the MOSDEF team's custom 2D reduction, the MOSFIRE Data Reduction Pipeline,¹⁰ and the 2D optical spectra from the Keck Low Resolution Imaging Spectrometer (LRIS, Oke et al. 1995; Rockosi et al. 2010). For the MOSDEF data, both optimally weighted and unweighted spectra were extracted for each object, and we use the optimally weighted spectra for this analysis. The optimal extraction algorithm follows Horne (1986) but is modified to be able to extract fractions of pixels (see the appendix). To determine the weighting profile, center, and width of each object, we fit a Gaussian to the profile of each object in each filter. The profile was determined by summing flux only at those wavelengths with high S/N in either the continuum or emission lines. Using high S/N areas of the spectra creates clean weighting profiles for the optimal extraction, as wavelengths with little or no signal are excluded.

Galaxies in the MOSDEF Survey are split into three redshift bins (1.37 < z < 1.70, 2.09 < z < 2.61, and 2.95 < z <3.80) that were each observed using a different filter and grating set in order to maximize efficiency of detecting multiple rest-optical emission lines of interest (see Kriek et al. 2015, for details). We combine the 1.37 < z < 1.70 and 2.09 < z <2.61 galaxies into a single sample (hereafter the z ∼ 2 sample) because these galaxies have coverage of Hα and [O iii]. There may be some evolution between the galaxies at these two redshift ranges, but without combining them, broad emission is extremely difficult to detect. The 2.95 < z < 3.80 galaxies do not have coverage of Hα and are stacked separately (hereafter the z ∼ 3 sample).

The parent data set contains 555 galaxies with measured redshifts, including 503 galaxies with [O iii] detections and 394 Hα detections. We create a sub-sample where we remove galaxies for which it may be difficult to accurately measure the broad emission or have AGN (discussed later). When cleaning the sample, we consider the Hα and [O iii] detections separately, except when considering galaxy-wide effects that are mergers and AGN presence.

• 1.  
  We remove 50 [O iii] and 35 Hα detections where the galaxy was an IR, X-ray, or both IR and X-ray detected AGN (see Coil et al. 2015 and Azadi et al. 2017). AGN are a possible source of outflows (Leung et al. 2017), and removing them allows us to isolate the effects of star formation on outflows.
• 2.  
  We remove 20 [O iii] and 26 Hα detections where the galaxy would be classified as an AGN based on z ∼ 0 optical line-ratio diagnostics (above the Kauffmann et al. 2003 line in the N2-BPT diagram). This is a conservative cutoff, considering galaxies at z ∼ 2 are offset compared to z ∼ 0 galaxies (Shapley et al. 2015).
• 3.  
  We remove 182 [O iii] detections where S/N < 10 for [O iii] and 122 Hα detections where S/N < 10 for Hα.
• 4.  
  We remove 22 [O iii] and 41 Hα detections where the [O iii] or Hα emission line was on or near bright skylines. These galaxies had a skyline on the edge of the line where the broad emission dominates. These galaxies had a nearby skyline, and the broad component fit was severely effected by residuals from the sky subtraction.
• 5.  
  We remove 14 [O iii] and 17 Hα detections where the galaxy appears to be undergoing a merger as indicated from the images or spectra of nearby objects. Mergers may have complicated kinematics and may not be well fit by our fitting method (described in Section 3.1). We determined which galaxies were mergers by inspecting both images and spectra by eye. If the galaxy was very misshapen or there was a nearby companion, we removed it. In the spectra, if there were two profiles that overlapped, then those galaxies were removed. If there was any evidence of a merger, the galaxy was removed.
• 6.  
  We remove 7 [O iii] and 8 Hα detections where the [O iii] or Hα emission is near the edge of wavelength coverage and the shape of the profile is difficult to determine. We removed galaxies where the center of the main emission line was less than 10 pixels away from the edge of the detector.
• 7.  
  We remove 5 [O iii] and 7 Hα detections, for which we measured a FWHM > 275 from a single Gaussian. Some galaxies have broad lines simply from rotation and velocity dispersion. To isolate the broad emission, we restrict the narrow emission to have a FWHM <275 km s⁻¹, as described in Section 3.1. Including galaxies with FWHM > 275 km s⁻¹ creates false positives because the narrow component does not properly fit the narrow emission.

The final sample has 216 unique galaxies with 203 [O iii] measurements and 138 Hα measurements. There are 125 galaxies with both an Hα and [O iii] detection. We create stacks (discussed in Section 3.2), and the stacks at z ∼ 2 have an additional restriction to contain galaxies with wavelength coverage of Hα, [O iii], [S ii], [N ii], and Hβ, which results in 113 galaxies in the z ∼ 2 stack. There are 60 galaxies in the z ∼ 3 stack.

Figure 1 shows histograms of redshift for the full sample of galaxies in blue, the remaining galaxies after the S/N rejection in red, and the the final sample in green. There is no clear bias against any particular redshift or mass; however, more galaxies with SFR < 10 ${M}_{\odot }$ yr⁻¹ have been removed. The right side of Figure 1 shows the SFR versus stellar mass diagram for our sample, along with the star-forming main sequence for MOSDEF galaxies in the 2.09 < z < 2.61 range from Shivaei et al. (2015). Because we removed galaxies that have an Hα or [O iii] S/N < 10, our sample may be incomplete for galaxies below the main sequence, at low stellar mass, or with high levels of dust.

Zoom In Zoom Out Reset image size

We estimate the physical parameters for our sample, including stellar mass, SFR, and age, by comparing the photometric SEDs with stellar population synthesis models (Conroy et al. 2009) using the stellar population fitting code FAST (Kriek et al. 2009). We assume a Chabrier (2003) initial mass function (IMF) and the dust reddening curve from Calzetti et al. (2000). We use spectroscopic redshifts from the MOSDEF survey and broadband and mediumband photometric catalogs assembled by the 3D-HST team (Skelton et al. 2014) spanning observed optical to mid-infrared wavelengths. We include a template error function to account for the mismatch in less constrained sections of the spectrum. For a full description of the stellar population modeling procedure, see Kriek et al. (2015).

When available, we derive SFRs based on the Hα emission line by correcting for Balmer absorption (using the SED) and dust extinction (using the Balmer decrement of Hα/Hβ), and then converting the Hα luminosity into a SFR (Kennicutt 1998), adjusted for a Chabrier (2003) IMF (see Shivaei et al. 2015 for more details). Because galaxies in the z ∼ 3.3 bin do not have coverage of Hα, we use SED fitting to determine their SFRs. For any galaxy with Hα but no coverage of Hβ, we also used the SED fitting to determine the SFR.

We aim to measure an underlying broad component of emission lines of galaxies by assuming that emission lines are composed of narrow and broad components with Gaussian profiles. This section describes the fitting process as well as constraints on parameters.

The Hα, [N ii], and [S ii] lines are in the same filter (H if at z ∼ 1.5 and K if at z ∼ 2.3), while [O iii] and Hβ fall into a different filter (J if at z ∼ 1.5, H if at z ∼ 2.3, and K if at z ∼ 3.3). We do not simultaneously fit lines that are in different filters because the spectral resolution is slightly different in each filter.¹¹ Additionally, the seeing may vary between different filters. We do not include [S ii] in fits of individual galaxies because it is too faint to measure the broad component; however, it is included in fits of stacks (see Section 3.2). For individual galaxies at 1.37 < z < 2.61, we fit Hα and [N ii] λλ6549, 6583 simultaneously. For individual galaxies at 1.37 < z < 3.8, we fit [O iii] λ5008, [O iii] λ4959, and Hβ simultaneously.

For each set of lines, we perform two preliminary fits and one final fit. The first preliminary fit uses a single Gaussian to fit each emission line using MPFIT, a non-linear least squares fitting code (Markwardt 2009). We use this single Gaussian fit to subtract off a linear continuum and normalize the data so the fitted flux density of the peak of the brightest line for each set of lines (Hα or [O iii]) is unity, respectively. Next, we fit the data again using MPFIT, but this time each emission line is fit with two Gaussians—one broad and one narrow. We use the resulting values and errors of this second fit as initial values for the final fit, which is done with a custom Markov Chain Monte Carlo (MCMC) code MPMCMCFUN.¹² How MCMC fitting works and the reasons for using it are described in the next two paragraphs.

MPMCMCFUN can fit any model to data using the Metropolis-Hastings algorithm (Metropolis et al. 1953; Hastings 1970) with a similar syntax to MPFIT. The best fit from MPMCMCFUN is the single model with the highest likelihood. For any implementation of the Metropolis-Hastings algorithm (including MPMCMCFUN), each parameter has a distribution of values from the Markov chain, and the error in each parameter can be calculated from this distribution. One way to do this is to fit a Gaussian function to the histogram of values, which gives a best-fit value as the centroid and an error from the σ of the fit. However, some parameters may have distributions that are non-Gaussian. For these distributions, the 1, 2, and 3σ error bars can be calculated by determining where 68%, 95%, and 99.7% of the distribution lies. If the distribution of values is very non-Gaussian (i.e., double peaked), it is best to show the entire distribution (also called the posterior); however, we do not need to do that for these models. Additionally, the correlation between parameters can be analyzed with a contour plot between two parameters.

The final MCMC fit is necessary because it offers a better characterization of errors and 3σ upper limits. In a least squares algorithm, only the error bar and covariance matrix are used to parameterize the error bars. This is especially relevant in the fits we do in this work because we are fitting one emission line with two Gaussian components and the correlation between parameters is significant. Additionally, when using MPFIT, if the best value for a parameter is the same as the limit set for that parameter, the error bar output is zero. This makes it impossible to calculate an upper limit. For this work, a 3σ upper limit is defined as the value where 99.85% of values in the posterior lie below. We consider a "significant" detection of the broad component if the lower 3σ limit on the broad component amplitude is larger than zero. The errors in the second preliminary fit are used as the parameter jump amplitudes (also called σ in step 1 of the preceding paragraph) for the final fit. For our data we used 100,000 iterations for initial testing and 400,000 iterations for the final fit. For the rest of this work, we use the subscripts S, B, and N to distinguish parameters for the single, broad, and narrow components, respectively.

When fitting multiple emission line components, we constrain the FWHM_B, FWHM_N, broad component shift (Δv), constant background, and the narrow component redshift to be the same for each line. This leaves each single emission line with two free parameters, broad amplitude (A_B), and narrow amplitude (A_N). Two exceptions to this are the [N ii] λ6585/[N ii] λ6550 and [O iii] λ5008/[O iii] λ4959 flux ratios, which are set to 2.93 and 2.98, respectively, according to atomic physics (Osterbrock 1989). Therefore, each fit has five free parameters shared by each line (FWHM_B, FWHM_N, Δv, narrow component redshift, and constant background) and two free parameters for each line (A_N and A_B).

The resulting best-fit parameters are likely to depend on the chosen limits. For instance, not placing a minimum on the FWHM_N can result in an unphysically narrow emission line. Also, not placing a minimum on the FWHM_B can result in the broad component not being representative of broadened emission. Therefore, we place physically motivated restrictions on all free parameters. For individual galaxies, we restrict the FWHM_N so that it cannot be lower than the average FWHM of skylines in that particular filter and mask. For galaxies smaller than the slit width, it is possible that the FWHM of the narrow component is smaller than that of skylines. In most cases, the seeing is not smaller than the width of the slit (07). Emission lines are also broadened by the spatial extent of the galaxy and the velocity distribution therein. Therefore, we do not expect FWHM_N to be much lower than the width of the skylines. We also restrict the FWHM_N to be less than 275 km s⁻¹. For this sample, we have removed galaxies where FWHM_S is larger than 275 km s⁻¹ (see Section 2.1).

In order to properly study outflows, we must be certain that the broad components measure a kinematically distinct feature from the rotation of the host galaxy. In other words, the broad flux must not be an artifact from a better fit to the narrow emission by using two Gaussian components. Therefore, we restrict the minimum FWHM_B to be a larger value than could be reasonably fit using only a single Gaussian component. Accordingly, we set the minimum FWHM_B to be 300 km s⁻¹, which provides some separation in the velocities of the narrow and broad components. Typical FWHMs for ionized outflows from star-forming galaxies are 300–600 km s⁻¹ (Genzel et al. 2011; Newman et al. 2012; Wood et al. 2015). Some studies have measured galactic outflow speeds >1000 km s⁻¹, but these are typically associated with AGN (Shapiro et al. 2009; Förster Schreiber et al. 2014; Genzel et al. 2014b). Since we have removed known AGN, we do not expect outflows of such high velocity. The upper limit of FWHM_B is set to 850 km s⁻¹, which is the typical maximum velocity deduced by the blueshifted interstellar absorption lines in the rest-frame UV of z ∼ 2 galaxies (Steidel et al. 2010).

When fitting with a single Gaussian component, the great majority of galaxies in this sample (95%) are under 275 km s⁻¹. We chose an upper limit of 275 km s⁻¹ for the narrow component based on this. Note the outflow velocities we are trying to measure, which are 300 to 600 km s⁻¹ (Genzel et al. 2011; Newman et al. 2012; Wood et al. 2015). With the separation of the narrow component and broad component of 25 km s⁻¹, we are certain that the broad component is actually detecting broad flux and is not simply a better fit to a single Gaussian. None of the galaxies with 3σ broad component detections run into the upper limit of the narrow component and are generally between 150 and 250 km s⁻¹. The limits of FWHM_N < 275 km s⁻¹ and FWHM_B > 300 km s⁻¹ are appropriate for our sample and have the downside of being insensitive to broad emission below 300 km s⁻¹. Other samples and fitting techniques might use a different definition of broad and narrow components.

We tested several other fitting restrictions, such as changing the limit on FWHM_N, changing the limit on FWHM_B, and setting the limit on FWHM_N and FWHM_B based on FWHM_S. Ultimately we decided on the threshold provided in the preceeding paragraphs because it is rather conservative in that it is more likely to not detect real broad emission rather than falsely detect broad emission when there is none.

When the amplitude of the broad component is consistent with zero, the center and width become unconstrained. Therefore it is important to restrict these parameters so they do not stray to unrealistic values. The centroid of the broad component is allowed to be anywhere within ±100 km s⁻¹ of the expected value. The broad component shift is the same for each line in the fit. No objects that had significant detections of the broad component ran into this limit. Other studies typically find shifts of <100 km s⁻¹ (Newman et al. 2012; Wood et al. 2015).

The amplitudes of the narrow Hα and [O iii] components are constrained to be between 0.2 and 1.05 (the peak of these lines were normalized to unity from the single Gaussian fit), and the broad component amplitude is constrained to −0.3 and 0.8. Since some galaxies may not show any signs of broad emission, the best value for the FWHM_B could be at or near zero for these galaxies. In these cases, it is still useful to put limits on the broad emission. Therefore, we allow the FWHM_B to be negative to fully sample the parameter space and set proper upper limits on the FWHM_B. In cases where the best value for the FWHM_B is less than zero, we interpret this galaxy as having no significant broad emission but still show the upper limit. For [N ii], [S ii], and Hβ we scale the restrictions to the relative peak of each line. All of the constraints are listed in Table 1.

Figure 2 shows fits for four galaxies that show evidence for a broad component. It is clear that a single Gaussian does not fit some of lines shown in these galaxies well, as evidenced by the "wave" pattern that is present in the residuals. The pattern shows the single fits underestimate flux at the peak, overestimate the wings, and underestimate the base. This pattern is particularly evident in the [O iii] lines of COSMOS-12015 and COSMOS-13015, and in the Hα lines of GOODS-N-12024 and GOODS-N-7231.

Zoom In Zoom Out Reset image size

The broad flux only dominates a small fraction of the line at high velocities. If the broad flux component is not approximately Gaussian, then the fitted broad flux might be different from what we measure. Other studies that analyze galaxies with higher signal to noise find that the broad emission is typically well fit by a Gaussian (Shapiro et al. 2009; Genzel et al. 2011; Newman et al. 2012).

The speed of outflows increases as a function of SFR and galaxy stellar mass, which is seen in observations (Martin 2005; Weiner et al. 2009) and simulations (Christensen et al. 2016; Muratov et al. 2015). At low outflow speeds (FWHM_B < 275 km s⁻¹), the emission from the outflow may be indistinguishable from the emission from the H ii regions. To quantify the detectability of low velocity outflows, we created two tests using simulated spectra where we can control the BFR and FWHM_B of galaxies used to make stacks. We can then measure the BFR and FWHM_B of the stacks and check if they are representative of the input galaxy parameters.

In the low velocity test (Test 1), the input FWHM_B increase from 100 to 800 km s⁻¹ between stellar masses of 10⁹ M_* and 10¹¹ M_*. For Test 1, galaxies below 10¹⁰ M_* have velocities below 300 km s⁻¹. In the high velocity test (Test 2), the input FWHM_B increase from 275 to 550 km s⁻¹ as a function of stellar mass. Both tests use the same distribution of BFRs and a FWHM_N of 200 km s⁻¹. The slope of the BFR versus M_* for the tests is the same as the slope of η versus M_* form FIRE simulations (Equation (8) from Muratov et al. 2015). We normalize this relation using simplified version of our Equation (2), η = C∗BFR and determine constant C using our highest mass stack (BFR = 0.68 at log M_* = 10.4), where the effects of low velocity outflows should be minimal. We use the distribution of SFRs and stellar masses of the sample (shown in Figure 1) to make stacks by mass and by SFR.

The results of these tests are shown in Figure 3. In Test 1, stacks underestimated the BFR for galaxies in the low and medium mass stacks. The stacks in Test 1 show an increase in BFR as a function of mass despite the input galaxies having a decrease with mass. The low and medium mass stacks contain 100% and 55% galaxies with FWHM_B < 275 km s⁻¹, respectively. The FWHM measured for these stacks is too large compared to the input galaxy values, although this is expected because we do not allow the FWHM_B parameter to go below 300 km s⁻¹ (see Section 3.4). The stacks in the test show that we could underestimate the extent of broad flux when the FWHM_B is <275 km s⁻¹ for a large fraction of galaxies within the stack. We cannot measure broad flux at low velocities with the measurements made in this work, but this is possible with the data in the FIRE simulation.

Zoom In Zoom Out Reset image size

The broad component is difficult to separate from the narrow emission. The galaxies in Figure 2 were chosen because they show the strongest evidence for broad emission. The faint, high velocity wings of the broad component are difficult to distinguish from the noise for the majority of the individual galaxies. In order to achieve a higher S/N, we create stacks of galaxies in bins of stellar mass, such that each stack has approximately the same number of galaxies. To create each stack, we interpolate the flux for each galaxy to a common rest-frame wavelength grid, subtract off any continuum, convert each spectrum from flux density to luminosity density, divide by the total luminosity of either Hα or [O iii] depending on the wavelength coverage of the stack, and then sum each spectrum with no weighting. To avoid adding significant noise from skyline subtraction residuals, we remove pixels associated with skylines where the error spectrum is above 1.5× the median error. The error in the stacked spectrum is calculated by making the stack 200 times but using input spectra with added Gaussian noise according to the associated error spectrum for each individual object; the error is calculated by taking the standard deviation of the 200 stacks at each wavelength. The z ∼ 2 stacks are shown in Figure 4.

Zoom In Zoom Out Reset image size

With this sample, it is also possible to create stacks in bins of SFR, sSFR, and Σ_SFR. These stacks are not truly independent from the stacks by mass because these properties correlate with mass. We chose to use mass for the primary analysis in this work because, of the physical parameters we considered, mass is the only parameter that is not estimated by using the Hα emission line, which may be influenced by broad emission.

The line-fitting process for stacks is the same as described in Section 3.1, with some exceptions. The doublet [S ii] λλ6718, 6733 is included when fitting Hα and [N ii]. For the lower limit of FWHM_N, we use the average skyline FWHM for each galaxy, which is 80 km s⁻¹. The fit to each stack is plotted in Figure 4. This figure shows Hα and [N ii] for the stacks at 1.37 < z < 2.61. For both Hα and [O iii] in all stacks, the amplitude of the broad component is significant at the >3σ confidence level.

The broad component dominates the line only at the highest velocities, which is also where the S/N is the lowest. Here we test the fitting process to show that measured broad line parameters are consistent with simulated input parameters.

For this test, we take a single Gaussian, add noise, and fit the Gaussian using the method described in Section 3.1. Since this idealized Gaussian has no actual broad component, any broad component that we measure is a false positive. We performed this test on 200 simulated emission lines, each with 10 different FWHMs between 75 and 275 km s⁻¹, which span the range of measured narrow components from the MOSDEF sample. We used the same resolution and wavelength as for an Hα line for a z ∼ 2.3 galaxy. We add a constant amount of noise to each spectrum such that the S/N of the Hα line ranges between 10 and 300. The assumption of constant noise is an appropriate approximation of the error for a single emission line in the actual spectra because we are limited by the bright sky rather than Poisson noise from the object. The noise does increase as a function of wavelength, particularly in the K band, but this test was only on a single emission line and the error does not change much over the span of a single emission line.

We find that only 11%, 1%, and 0.05% of simulated galaxies have a false positive of 1σ, 2σ, or 3σ, respectively. These are lower than the expected rates based on Gaussian statistics, which are 16%, 2%, and 0.1%. The slightly lower values of the test result are because we allow the narrow peak to exceed 1.0 (the max is 1.05). The average broad component is slightly less than 0, which creates an offset in the number of 1σ, 2σ, or 3σ detections. Another result of this test is that the fraction of false positives did not change as a function of width or noise added.

In addition to false positives in individual spectra, there is a possibility of creating an artificial broad component when making the stacks. We performed several tests to ensure that in creating the stacks, we did not also create an artificial broad component. The first test takes 50 Gaussians of random FWHM between 75 and 250 km s⁻¹, adds noise, and creates a stack as described in Section 3.2. Each added Gaussian has a constant amount of noise across each wavelength element that is similar to the level in actual spectra, except for skylines, which are not included in these simulated stacks. We found no evidence of introducing false positives when creating stacks. We repeat this test and add a random shift between ±1 Å to the centroid of the Gaussian, which simulates imperfect redshift estimates. The average redshift error for this sample is 6 × 10⁻⁵, which is ∼0.13 Å at these redshifts. These stacks also failed to produce false positives.

These tests have shown that we do not expect false positives to be an issue when using the fitting method described in Section 3.1. We have also shown that creating stacks of galaxies does not introduce a broad emission signature.

Previously, we stated that we constrained the FWHM_B, FWHM_N, broad component shift (Δv), and the narrow component redshift to be the same for each line. Unless there is some extreme metallicity gradient from the gas emitting these components, we believe this is a reasonable assumption. The kinematics of [N ii] is very likely to match the kinematics of the Hα emission, and the kinematics of Hβ is very likely to match the kinematics of the [O iii] emission. It is not clear if other works have made these same assumptions about the widths of the emission lines (Newman et al. 2012, 2014).

These assumptions have an impact on the significance of our fits in the weaker lines [N ii] and Hβ. For the stronger lines (Hα and [O iii]), the single Gaussian fit does not adequately explain the flux profile. This is evidenced by the "wave" pattern in the residuals (see Figure 2). However, for the weaker lines, no such wave pattern exists from the single Gaussian fit, indicating that this fit is a sufficient fit to those data. In every case where there is a significant broad component (in stacks or individual galaxies), FWHM_N < FWHM_S by 25 to 50 km s⁻¹. So, even though the single Gaussian model is sufficient to fit the emission of the weak line, the narrow component of the narrow+broad fit is typically a poor fit because FWHM_N < FWHM_S. The narrow component alone does not do a good job of fitting the flux of the weak line, which makes the broad component of the weaker lines significant in most stacks.

If we did not constrain FWHM_B, FWHM_N, Δv, and the narrow component redshift to be the same for each line, then for [N ii] and Hβ in the stacks, the FWHM_N would be equal to FWHM_S. In this case, the broad component would not be significant for [N ii] and Hβ because the narrow component would be a sufficient fit to the emission line. Since we have a physical motivation for locking the velocities of each line to one another, we assert that the broad emission in the weaker lines is significant even if the single Gaussian fit is a sufficient fit to the data.

After fitting each galaxy and stack, we parameterize the broad emission we measured as a broad to narrow flux ratio (broad flux/narrow flux, BFR). We chose this parameterization because other studies have used this and using the same parameterization allows for easy comparison (e.g., Newman et al. 2012). The BFR is also used to estimate the mass loading factor (Section 5.3). The other natural parameterization, broad flux to total flux, can be calculated as $(\mathrm{broad}\ \mathrm{flux}/\mathrm{total}\ \mathrm{flux})=1/(1+{\mathrm{BFR}}^{-1})$.

The left side of Figure 5 shows the BFR measured from the Hα line as a function of mass. For individual galaxies, there are 10 detections with >3σ significance out of 138 galaxies (7%). For the stacks and the galaxies with detections, the broad flux accounts for 10%–70% of the total flux in nebular emission lines. The MOSDEF measurements for BFR are consistent with the measurements from Newman et al. (2012), who did a similar analysis for galaxies at the same mass range. The details of the fits are in Table 2. The small differences in the number of significant detections in this study and Leung et al. (2017) for the same sample can be attributed to slight differences in codes used to fit the data.

Zoom In Zoom Out Reset image size

Table 2.  Measurements from the Stack Fits

Parametera	Avg.b	Rangec	BFRd	BFRmine	BFRmaxf	FWHMBg	ΔVBh	ηi
(Hα)-z ∼ 2	9.56	8.97–9.80	0.15 $\displaystyle \genfrac{}{}{0em}{}{+0.071}{-0.041}$	0.056	0.35	300 ± 200	−9.0 ± 20	0.26 $\displaystyle \genfrac{}{}{0em}{}{+0.15}{-0.086}$
(Hα)-z ∼ 2	9.95	9.82–10.1	0.40 $\displaystyle \genfrac{}{}{0em}{}{+0.092}{-0.15}$	0.15	0.70	300 ± 60	16. ± 10	0.64 $\displaystyle \genfrac{}{}{0em}{}{+0.21}{-0.28}$
(Hα)-z ∼ 2	10.3	10.1–10.7	0.68 $\displaystyle \genfrac{}{}{0em}{}{+0.15}{-0.18}$	0.35	1.2	340 ± 30	−10. ± 6	1.4 $\displaystyle \genfrac{}{}{0em}{}{+0.41}{-0.42}$
([O iii])-z ∼ 2	9.53	8.97–9.80	0.11 ± 0.028	0.029	0.20	390 ± 100	−63 ± 40	⋯
([O iii])-z ∼ 2	9.92	9.82–10.1	0.56 ± 0.12	0.21	0.92	300 ± 20	−18 ± 8	⋯
([O iii])-z ∼ 2	10.4	10.1–10.7	0.58 ± 0.17	0.060	1.1	300 ± 30	8.5 ± 9	⋯
([O iii])-z ∼ 3	9.5	9.2–9.8	0.75 ± 0.14	0.32	1.2	340 ± 30	−13 ± 9	⋯
([O iii])-z ∼ 3	10.0	9.82–10.1	0.37 ± 0.16	0.08	0.85	300 ± 100	27 ± 20	⋯
([O iii])-z ∼ 3	10.5	10.1–11.0	0.71 ± 0.37	0.1	1.8	320 ± 80	25 ± 20	⋯

Notes.

^aParameter by which the stack was created. ^bThe geometric mean of the galaxies included. ^cMass range of galaxies included. ^dBroad flux fraction of the stack and 1σ errors. ^eThe 3σ lower limit on the BFR. ^fThe 3σ upper limit on the BFR. ^gThe FWHM of the broad component (km s⁻¹) and the 1σ error. ^hThe velocity offset between the peaks of the broad and narrow components (km s⁻¹). A negative value indicates a blueshift. The 1σ error in the velocity offset is included. ⁱThe mass loading factor (see Section 5.1).

Download table as:  ASCIITypeset image

The right side of Figure 5 shows the BFR measured from the [O iii] lines as a function of mass. For [O iii] there are 21 detections with >3σ significance out of 201 galaxies (10%). For the stacks and galaxies with detections, the broad flux accounts for 20% to 50% of the total flux in nebular emission lines. For the z ∼ 3 stacks, the broad component comprises 30% to 60% of the flux and the BFR is slightly higher than in the z ∼ 2 stacks on average.

The stacks in Figure 5 show an apparent correlation between the BFR and mass. Measuring a broad component of the emission line relative to the narrow component in the lowest mass stack is difficult because most of the broad emission may be at FWHM < 275 km s⁻¹, which is hard to measure due to our assumptions in the fitting process (as discussed in Section 3.5). Non-measurement of the broad fluxes for the lowest mass galaxies may introduce a bias in the BFR versus stellar mass relation. Additionally, the [O iii] broad emission at z ∼ 2 does not show any increase above 7 × 10⁹ ${M}_{\odot }$, and the z ∼ 3 stacks show no change with mass. There is also no correlation between the detected broad emission in individual galaxies and mass. For these reasons, we cannot confirm a correlation between the BFR and mass.

Figure 6 shows the BFR measured from the Hα line as a function of S/N. Here, we define the S/N as the fitted flux by a single Gaussian divided by the error in the flux for either Hα or [O iii] λ5007. Galaxies with detections have higher S/N. For Hα, 66% of galaxies with S/N > 70 have broad component detections, but only 1.6% of galaxies with S/N < 70 have detections. For [O iii], 32% of galaxies with S/N > 45 have broad component detections, but only 5% of galaxies with S/N < 45 have detections. These two thresholds were chosen by eye to emphasize the S/N dependence on detecting the broad component. It is easier to detect broad emission in [O iii] than in Hα because we include the [O iii] λ4959 emission when fitting [O iii] λ5007. Using both lines in the fit provides a better constraint on the shape of the broad and narrow emission profiles than using only one line.

Zoom In Zoom Out Reset image size

The dependence of the detection of the broad flux on S/N implies that the 10% detection rate is a lower limit. Because outflows are supposedly ubiquitous at z ∼ 2, we would likely see more broad component detections with deeper data. A dependence on S/N and a broad component detection was also seen in Leung et al. (2017) for AGN in the MOSDEF sample.

As described in Section 3.1, we fit narrow and broad components to the [O iii], Hβ, Hα, [N ii], and [S ii] emission lines in stacked spectra. From this analysis, we are able to calculate the [N ii]/Hα, [S ii]/Hα, and [O iii]/Hβ ratios, and place each component on the N2-BPT and S2-BPT diagrams. Figure 7 shows the N2-BPT and S2-BPT diagram for the low, medium, and high mass stacks. The mass range of each stack is in the first three rows of Table 2. We do not include individual galaxies here because there were not enough 3σ detections of the broad components of Hβ, [N ii], and [S ii], and we could not create robust line ratios. The blue dashed line is measured from Kewley et al. (2013) for local galaxies. The orange dashed line is measured from Shapley et al. (2015) for z ∼ 2.3 galaxies in the MOSDEF survey. The dashed black line is from Kauffmann et al. (2003), and separates star-forming galaxies and AGN in the local universe. The dotted black line is from Kewley et al. (2001) and is the "maximum starburst" line where galaxies containing AGN lie above this line.

Zoom In Zoom Out Reset image size

For individual galaxies, Balmer emission line fluxes can be corrected for underlying stellar absorption based on the equivalent widths of stellar Balmer features as estimated from the stellar population synthesis model fit to the SED of each galaxy (Reddy et al. 2015). For each stack, we estimate the Balmer absorption by calculating the average absorption for each galaxy in the stack. This gives us an estimate for the total fraction of flux that was absorbed but no information about the shape. Without knowing the exact shape/width of the absorption feature, we do not know how much of the correction should be applied to the narrow feature and how much should be applied to the broad feature. Therefore, we calculate the Balmer absorption correction assuming the broad component is affected by 0%, 33%, 66%, and 100% of the Balmer absorption and the narrow component is affected by 100%, 66%, 33%, and 0% respectively. This gives a general idea of the most the Balmer absorption could affect each line ratio. In Figure 7, the 0% and 100% absorption cases correspond to the hollow point and the solid point furthest from the hollow point, respectively. For the single Gaussian fits (square points), there is only one solid point because the Balmer emission is not split between narrow and broad components and the magnitude of the Balmer absorption correction is unambiguous. The shape of the Balmer absorption may also affect the fits in a manner that is difficult to predict.

In Figure 7, the ratios from the single Gaussian fits (squares) are lower than results from previous MOSDEF studies (Shapley et al. 2015). This can be explained by the fact that we required a S/N > 10 for the Hα and [O iii] lines. This requirement preferentially removed lower mass galaxies, which are typically more offset from the local relation (Shapley et al. 2015).

The narrow component ratios (triangles) tend to lie more toward the z ∼ 0 relationship (blue dashed line, Kewley et al. 2013) than other measurements at z ∼ 2.3. The broad components of the narrow+broad fits (diamonds) lie in the composite region or above the Kewley et al. (2001) line. The Balmer correction is large for the [O iii]/Hβ ratio, and this makes it difficult to conclude if the broad components have higher [O iii]/Hβ than their narrow counterparts for the low and high mass stacks. The broad component of the medium mass stack lies in the composite region even after Balmer correction. The Balmer correction is smaller for the [N ii]/Hα ratio, and it is clear that the broad component for the highest mass stack has higher [N ii]/Hα than the narrow components, even after Balmer absorption correction. The differences between the broad and narrow components vary between mass bins. This is likely because the fraction of flux in the broad component increases as a function of mass (Figure 5).

The right side of Figure 7 shows the S2-BPT diagram for each stack and for each component. The [S ii] line typically has less flux than the [N ii] line, making measuring the broad component more difficult. We are only able to place 1σ limits on the broad components of the [S ii] line in stacks. Nevertheless, these limits are consistent with a higher [S ii]/Hα ratio for the broad components.

In addition to the fits, we calculated line ratios by integrating the flux in several velocity bins with respect to the centroid of the line. This measurement provides a non-parametric estimation of the line ratios that is independent of any model and is shown in the appendix. We found that the higher velocity bins generally had higher [N ii]/Hα and [O iii]/Hβ ratios, which is consistent with what we measure with the fits.

Emission line ratios from shocks differ from ratios in photoionized gas. Shock-heated gas can become ionized by high-energy photons from the shock or excited by collisions. Emission line ratios shift in the presence of shocks and the magnitude and direction of the shift depends on the metallicity, electron density, magnetic field, and shock velocity (Allen et al. 2008). Shocked emission tends to have higher [N ii]/Hα and [S ii]/Hα ratios relative to what is produced in photoionized H ii regions (Allen et al. 2008). Since the broad components in Figure 7 have higher [N ii]/Hα ratios than the narrow components or single Gaussian fits, this may indicate the presence of shocks. In this section, we investigate if the broad emission can be explained by shocks by creating the N2-BPT and S2-BPT diagrams using data from the shock models by Allen et al. (2008) ¹³ and comparing these models to the broad emission line ratios.

The shock models simulated emission line ratios for shocked gas, the precursor to the shock, and a shock+precursor, which combines the shock and precursor components. The precursor is material that is photoionized by the shock but not directly shocked itself. Because we do not spatially resolve the emission from these galaxies, we are unable to separate the different components of the shock. Therefore, we compare our measurements of the broad emission to the combined shock+precursor ratios.

Allen et al. (2008) measured shock+precursor emission line ratios for two sets of models: one at a fixed electron density with varying metallicity (n_e = 1 cm⁻³ at log(O/H)+12 of 8.03, 8.35, 8.44, and 8.93), and another at fixed metallicity with varying electron density (log(O/H)+12 = 8.93 at n_e = 1, 10, 100, 1000 cm⁻³). We restrict the models shown to those that have a magnetic field strength at pressure equipartition. The shock velocity for the models range from 100 to 1000 km s⁻¹, but we only show shock velocities of 200–500 km s⁻¹ based on the velocities measured in Table 2.

In Figure 8, we show the shocked models for the N2-BPT and S2-BPT diagrams. The top row shows the effect of changing metallicity on shocked diagnostic ratios, and the bottom row shows the effect of changing density. The galaxies in this sample (with S/N of [N ii] > 3) have a median metallicity of log(O/H)+12 = 8.43, with 80% of galaxies between 8.27 < log(O/H)+12 < 8.59 calculated using the [N ii]/Hα ratio as in Sanders et al. (2015). The electron density of the MOSDEF galaxy sample at 2.09 < z < 2.61 is ${290}_{-169}^{+88}$ cm⁻³ (Sanders et al. 2016). This was calculated using the entire [S ii] line and the electron density of the material causing the broad emission may be different. Newman et al. (2012) measured a density of ${10}_{-50}^{+590}$ cm⁻³ from a stack of 14 galaxies, and this value is consistent with our assumption of 290 cm⁻³. Since none of the simulations span exactly the range of metallicities and densities of the MOSDEF galaxies, we are forced to extrapolate between the effects of metallicity and density. We highlight the point that is the best match to the metallicity, electron density, and shock velocity in green. This green point corresponds to the shock model which has (v =300 km s⁻¹, log(O/H)+12 = 8.44, n_e = 1cm⁻³) in the top row and (v = 300 km s⁻¹, log(O/H)+12 = 8.93, n_e = 100 cm⁻³) in the bottom row.

Zoom In Zoom Out Reset image size

In Figure 8, there is a strong metallicity dependence on the [N ii]/Hα ratio, and there is almost no change as a function of electron density except at the highest density (n_e = 1000 cm⁻³) where [N ii]/Hα decreases. The broad lines measured from stacks are consistent with the [log(O/H)+12 = 8.44, n_e =1 cm⁻³] and [log(O/H)+12 = 8.35, n_e = 1 cm⁻³] points. The shock model that best matches the physical parameters (green point) is very near the broad emission line ratios. Therefore, it is feasible that the positions of the broad components in the N2-BPT diagram can be explained by shocks.

In the top row of Figure 8, the model that best matches the MOSDEF data has a higher [S ii]/Hα ratio than any of the broad emission. This may be due to the limitation that the models with varying metallicity have an electron density of 1 cm⁻³. The models in the bottom row show that [S ii]/Hα decreases as electron density increases. Since the electron density of the MOSDEF galaxy sample is ${290}_{-169}^{+88}$ cm⁻³, the models in the top row would likely shift to lower [S ii]/Hα ratios at higher electron densities. Therefore, it is possible that the positions of the broad components in the S2-BPT diagram can be explained by shocks.

From these line ratios, we conclude that it is possible that the broad emission is a result of shocked emission. This explanation does not rule out other sources of emission, such as AGN (discussed in Section 5.2) and photoionized outflows (discussed in Section 5.3).

When creating the sample presented in this work, we removed all X-ray, IR, and optically identified AGN because our goal is to study star formation driven outflows, and AGN are also known to drive outflows at z ∼ 2 (e.g., Förster Schreiber et al. 2014; Leung et al. 2017). However, there may be low-luminosity AGN that were not detected with these methods. There is an observational bias against identifying AGNs at all wavelengths in low-mass galaxies (Azadi et al. 2017). This bias may lead to some galaxies that host AGN being included in the sample. AGN typically have higher [N ii]/Hα and [O iii]/Hβ ratios because of harder ionization coming from the accretion disk, which is consistent with the line ratios we find in the broad component.

With integral field spectroscopy, it is possible to create spatially resolved line ratios (Wright et al. 2010; Newman et al. 2014), which can be used to determine if a galaxy in the "composite" region of the N2-BPT diagram has an AGN. Using spatially resolved emission line maps, Newman et al. (2014) found some galaxies that lie in the composite region of the BPT diagram host AGN. The cores of these galaxies lie in the AGN region, while the outer edges lie in the star-forming region. The high [N ii]/Hα and [O iii]/Hβ ratios of the core indicated the presence of an AGN that would not be detected in spectra of low spatial resolution.

Although we have removed all X-ray, IR, and optically identified AGN from this sample, we cannot completely rule out some contribution to the emission lines from low-luminosity black holes. We used a conservative cutoff for optically identified AGN candidates. If low-luminosity AGN mixing is a significant source of emission in the stacks, then the presence of AGN would have to be extremely widespread among z ∼ 2 galaxies that otherwise appear to be dominated by star formation alone (Coil et al. 2015). Since SFRs are higher at z ∼ 2 than at z ∼ 0, weak AGN would not be bright enough to significantly change line ratios (Coil et al. 2015). It seems more likely that shocks could be commonplace at z ∼ 2 due to the high SFRs, instead of AGN being ubiquitous. If there is any AGN contribution, it is likely small.

It is unlikely that we are detecting outflows driven primarily by AGN because AGN-driven outflows have more extreme kinematics compared to star formation driven outflows. AGN-driven outflows are typically 500 to 5000 km s⁻¹ (Förster Schreiber et al. 2014; Genzel et al. 2014b; Leung et al. 2017), which is much faster than typical outflow velocities from star-forming galaxies (300–550 km s⁻¹; Shapiro et al. 2009; Newman et al. 2012). Additionally, the velocity difference between the narrow and broad components in AGN-driven outflows is 100 to 500 km s⁻¹ (Leung et al. 2017), while star formation driven outflows typically have velocity offsets of <100 km s⁻¹ (Newman et al. 2012). Since the lower end of the velocity range of AGN outflows overlaps with the velocity range in outflows and broad components measured here, this argument does not completely exclude some AGN contribution.

One interpretation of the broad component is that it traces ionized outflowing materials (Bland & Tully 1988; Heckman et al. 1990; Gergeev 1992; Phillips 1993; Lehnert & Heckman 1996; Veilleux et al. 2001; Colina et al. 2004; Westmoquette et al. 2007, 2008; Shapiro et al. 2009; Newman et al. 2012; Rupke & Veilleux 2013; Genzel et al. 2014a; Feruglio et al. 2015; Leroy et al. 2015). In this section, we interpret the broad component as a photoionized outflow, calculate the mass loading factor η (outflow mass rate/SFR), and compare to other observations as well as simulations.

Using some assumptions about the outflow velocity, radius, temperature, and density, we can convert the BFR into an estimate of η. We adopt the outflow model from Genzel et al. (2011). This model assumes the broad component to be photoionized and the emission of Hα to be a result of case-B recombination.

The assumptions that the broad component is an outflow and that the broad component is shocked (as we assumed in the previous sub-section) are not mutually exclusive. The broad component could be a shocked outflow. There could also be a mix of photoionized and shocked gas (Newman et al. 2014). If we assume a 100% collisionally ionized outflow, the mass loading factor would be smaller by a factor of ∼2 (see the appendix of Genzel et al. 2011).

The model assumes a spherical outflow with a constant velocity (see Steidel et al. 2010). The mass outflow rate, ${\dot{M}}_{\mathrm{out}}$, can be calculated as

$$\begin{eqnarray}&&{\dot{M}}_{\mathrm{out}}=\displaystyle \frac{1.36{m}_{{\rm{H}}}}{{\gamma }_{{\rm{H}}\alpha }{n}_{{\rm{e}}}}\left({L}_{{\rm{H}}\alpha }\displaystyle \frac{{F}_{\mathrm{broad}}}{{F}_{\mathrm{narrow}}}\right)\displaystyle \frac{{V}_{\mathrm{out}}}{{R}_{\mathrm{out}}},\end{eqnarray} \tag{ 1 }$$

where m_H is the atomic mass of hydrogen, V_out is the velocity of the outflow, R_out is the radius, γ_Hα is the Hα emissivity, n_e is the electron density in the outflow, L_Hα is the total extinction corrected Hα luminosity, and F_broad/F_narrow is the BFR.

We attempt to measure each component of Equation (1) from our stacks (as described later), but sometimes we do not have sufficient signal to do so. For physical parameters we cannot estimate, we adopt values from Newman et al. (2012; hereafter N12) who preformed a similar analysis on 27 star-forming galaxies at z ∼ 2. The sample from N12 has a similar mass range to our sample but has higher SFRs (∼90 ${M}_{\odot }$ yr⁻¹ on average), which may lead to physical differences. However, N12 is currently the most similar study to ours with measurements of the parameters in Equation (1), and we use their values when we are unable to measure them from our data.

The electron density for the outflow can be measured using the broad [S ii] λ6718/[S ii] λ6733 ratio (Osterbrock 1989; Newman et al. 2012; Sanders et al. 2016). We attempt to measure this ratio for the broad components for the stacks. The flux in the [S ii] lines is low, which results in a large measurement uncertainty in the ratio. We are unable to constrain the density using the broad component from this work. We adopt the value used by N12 of ${50}_{-50}^{+550}$ cm⁻³.

The term V_out/R_out is the inverse of the characteristic timescale of the outflow. In an attempt to measure the radius of the outflow, we made a stack of the 2D spectra and attempted to find broad flux in the spatial direction (e.g., Martin 2006; Leung et al. 2017). We were unable to measure a spatially extended component in the stacked spectrum. For R_out, we adopt the value of 3 kpc as measured by N12. This value is reasonable, given the angular size at this redshift is ∼8 kpc arcsec⁻¹ and our best seeing is ∼06. For V_out, we use the "maximum" velocity of the outflow defined as ${V}_{\max }\,=| {\rm{\Delta }}{v}_{{\rm{B}}}$–$2{\sigma }_{{\rm{B}}}|$ (Genzel et al. 2011; Wood et al. 2015). This value represents the velocity of the outflow if one assumes the outflow is spherically symmetric with a constant velocity.

We use an Hα emissivity of 3.56 × 10⁻²⁵ erg cm³ s⁻¹, which assumes an electron temperature of T_e = 10⁴ K.

If we use the Kennicutt (1998) relation between SFR and Hα luminosity corrected for a Kroupa IMF (SFR[${M}_{\odot }$ yr⁻¹] =7.9 × 10⁻⁴² L_Hα [erg s⁻¹]), we can divide Equation (1) by SFR and simplify to

$$\begin{eqnarray}&&\eta \approx 2.0\left(\displaystyle \frac{50\ {\mathrm{cm}}^{-3}}{{n}_{{\rm{e}}}}\right)\left(\displaystyle \frac{{V}_{\mathrm{out}}}{300\ \mathrm{km}\ {{\rm{s}}}^{-1}}\right)\left(\displaystyle \frac{3\ \mathrm{kpc}}{{R}_{\mathrm{out}}}\right)\left(\displaystyle \frac{{F}_{\mathrm{broad}}}{{F}_{\mathrm{narrow}}}\right).\end{eqnarray} \tag{ 2 }$$

Table 2 lists the η calculated which is between 0.26 and 1.4. The error calculation includes measurement uncertainties from the BFR and the FWHM_B. We do not include errors in the radius, electron density, and temperature assumed, and including these errors would increase the error on the mass loading factor by at least an order of magnitude.

The mass loading factor for the highest mass stack ($\eta ={1.4}_{-0.42}^{+0.41}$) is lower than the predicted value from FIRE at that mass (η = 2.6, Muratov et al. 2015). The FWHM_B measured for this stack is 340 ± 30 km s⁻¹ (Table 2), so the difficulty of detecting low velocity outflows should not be a factor. Using a smaller electron density or smaller radius in Equation (2) would bring these into better alignment; however, there is no evidence to justify such changes. An alternate explanation for why η is lower than what is measured in the FIRE simulations is that some fraction of the outflow is neutral and not visible as a broad Hα component. Outflows measured in the Hα line are a measure of the ionized component of the outflow, but outflows are multi-phase and have neutral, ionized, and dusty components (Wood et al. 2015; Martín-Fernández et al. 2016). An undetected neutral component may account for the factor of 2 difference between what is measured here and the FIRE simulations.

The errors and assumptions needed to calculate η make any comparison of our estimate to theoretical or other observational estimates difficult. Additionally, outflows measured in the Hα line are a measure of the ionized component of the outflow, but outflows are multi-phase and have neutral, ionized, and dusty components (Feruglio et al. 2015; Leroy et al. 2015; Wood et al. 2015).

The abbreviations introduced in this section are

$$\begin{eqnarray*}&&{\rm{N}}2=[{\rm{N}}\,{\rm{II}}]/{\rm{H}}\alpha \end{eqnarray*}$$

$$\begin{eqnarray*}&&{\rm{O}}32=[{\rm{O}}\,{\rm{III}}]\lambda \lambda 4959,5007/[{\rm{O}}\,{\rm{II}}]\lambda \lambda 3726,3729\end{eqnarray*}$$

$$\begin{eqnarray*}&&{\rm{O}}3{\rm{N}}2=([{\rm{O}}\,{\rm{III}}]\lambda \lambda 4959,5007)/([{\rm{N}}\,{\rm{II}}]/{\rm{H}}\alpha )\end{eqnarray*}$$

$$\begin{eqnarray*}&&{\rm{R}}23=([{\rm{O}}\,{\rm{III}}]\lambda \lambda 4959,5007+[{\rm{O}}\,{\rm{II}}]\lambda \lambda 3726,3729)/{\rm{H}}\beta ).\end{eqnarray*}$$

For the MOSDEF sample, Shapley et al. (2015) and Sanders et al. (2016) showed that the z ∼ 2 galaxies are offset from the z ∼ 0 galaxies in diagnostic diagrams that include [N ii]. Specifically, the galaxies were offset in the N2-BPT, O32 versus O3N2, and O32 versus N2 diagrams, and did not show any significant offset in the O32 versus R23 and S2-BPT diagrams. While there is much speculation, there is no definitive explanation for the offset in diagrams that include [N ii] (e.g., Masters et al. 2014, 2016; Steidel et al. 2014, 2016; Shapley et al. 2015; Sanders et al. 2016). In this section, we test if the offsets in the diagnostic diagrams could be caused only by adding the emission from shocks to the z ∼ 0 spectra. Using the distribution of galaxies at z ∼ 0 as an estimate of the intrinsic emission of galaxies at z ∼ 2 is merely an approximation because z ∼ 2 galaxies have different properties that affect their location on diagnostic diagrams, such as higher electron density (Sanders et al. 2016) and lower metallicity at fixed stellar mass and SFR (Sanders et al. 2018).

Figure 9 shows the O32 versus R23, O32 versus O3N2, and O32 versus N2 diagrams with data at z ∼ 0 from SDSS in black and at z ∼ 2 from Sanders et al. (2016). We overlay the same shock models shown in Figure 8. The shock models in Figure 9 show the same general trend as the z ∼ 2 data when compared to the SDSS data: no clear offset in O32 versus R23, a slight offset in O32 versus O3N2, and a large offset in O32 versus N2.

Zoom In Zoom Out Reset image size

The electron density affects the location of the shocked emission, and we consider two measurements of the electron density. The first is a value of ${290}_{-169}^{+88}$ cm⁻³ (Sanders et al. 2016) from the MOSDEF sample. The second is a value of ${50}_{-50}^{+550}$ cm⁻³ from N12. Both of these densities were measured using [S ii] λ6718/[S ii] λ6733, but the N12 values were for the broad components in their stacks, and the value from Sanders et al. (2016) is from the entire line. If we believe the shocks are within the galaxy, at the base of an outflow (as described in Wood et al. 2015), then we should choose 290 cm⁻³ from Sanders et al. (2016). If the shocks are originating in outflows kiloparsecs from the galaxy, then we should use the value of 50 cm⁻³ from N12. Because the models by Allen et al. (2008) only have electron densities of n_e = 1, 10, 100, 1000 cm⁻³, we pick the 100 cm⁻³ model, which is between the two measurements.

To add the flux from shocks to local data, we select two shock models that are closest to the mean electron density, metallicity, and shock velocity of the MOSDEF sample (${290}_{-169}^{+88}$ cm⁻³ Sanders et al. 2016, log(O/H) +12 = 8.43, Sanders et al. 2015, and shock velocity of 300 km s⁻¹). These models are shown as green points. We then add the SDSS distribution and shocked data point together, assuming a BFR of 0.4 (which is the average BFR for the stacks by Hα and corresponds to 29% of the total flux being shocked) and plot the result as a green line. These SDSS+shocks models represent what local galaxies would look like with the addition of the best-fitting shock model from Figure 8.

The SDSS+shock data in Figure 9 generally show good agreement with the z ∼ 2 galaxies: no clear offset in O32 versus R23, a slight offset in O32 versus O3N2, and a large offset in O32 versus N2. The O32 versus R23 and O32 versus N2 diagrams for the n_e = 100 cm⁻³ SDSS+shock models (bottom row) are higher than expected but could be explained because these models are at a metallicity of log(O/H)+12 = 8.93, which is higher than the average z ∼ 2 galaxy metallicity. Generally, the N2 ratio decreases with decreasing metallicity and the R23 ratio increases with decreasing metallicity. A decrease in the metallicity of these two models would bring them into better agreement with the z ∼ 2 galaxies.

These models assume a single shock velocity, electron density, and metallicity for the whole sample. Galaxies at z ∼ 2 have a wide range of metallicities (Sanders et al. 2015), electron densities (Sanders et al. 2016), and shock velocities (Newman et al. 2014, Table 2). It is likely that the broad components occupy some region of the N2-BPT diagram based on these physical properties, not just a single point. We can constrain some parameters by looking at which shock models are unreasonable compared to the z ∼ 2 data in Figure 9. At velocities above 400 km s⁻¹, the SDSS+shock models would be significantly higher than the z ∼ 2 data in the O32 versus N2 diagram, and at velocities below 250 km s⁻¹, the SDSS+shock models would be lower than the z ∼ 2 data in the O32 versus N2 and O32 versus O3N2 diagrams. Changing the metallicity of the shock models (top row) does not result in large shifts in line ratio space, except in the O32 versus N2 diagram. All of the metallicities would fit the data, except for log(O/H)+12 = 8.03, which would be too low in the O32 versus N2 diagram. There is little change with electron density in Figure 9 (bottom row), but there is a large dependence in the S2-BPT diagram (Figure 8). We can only place limits on the broad [S ii]/Hα, and the 3σ limits are consistent with the lowest electron density models. Therefore, we conclude that all of the electron densities in the SDSS+shock models would match the z ∼ 0 to z ∼ 2 offsets.

For the O32 versus R23, O32 versus O3N2, and O32 versus N2 diagrams, adding shocks to SDSS data at z ∼ 0 could shift the line ratios toward the values measured at z ∼ 2. Shock velocities of 250 to 400 km s⁻¹ with metallicities ranging from log(O/H)+12 = 8.35–8.93 and electron densities between n_e = 1–1000 cm⁻³ are plausible. This does not prove that the offset is caused by shocks—only that they are a possibility. Our results do not rule out contribution from AGN as a driver of these offsets, as AGN occupy similar regions of diagnostic line-ratio diagrams as shocks.

A great deal of study has been done on the cause of the offset of z ∼ 2 galaxies from the z ∼ 0 galaxies in the N2-BPT diagram (Shapley et al. 2005, 2015; Erb et al. 2006; Liu et al. 2008; Steidel et al. 2014; Masters et al. 2016; Strom et al. 2017). Here, we calculate the same shock+SDSS models for the N2-BPT diagram and calculate where the broad component should lie if it is due to shocks.

Figure 10 shows the SDSS+shock model along with the location of the z ∼ 2 galaxies from Shapley et al. (2015). The SDSS+shock model lines show generally good agreement with the line from Shapley et al. (2015). This implies that adding shocks to the spectra of local galaxies could, in part, explain the offset of the z ∼ 2 galaxies. This conclusion comes with the caveat that we are unable to completely rule out some contribution from AGN in our stacks. Since AGN have similar line ratios to shocks, the same argument holds that low-luminosity AGN (instead of shocks) could explain the offset of the z ∼ 2 galaxies.

Zoom In Zoom Out Reset image size

As shown in Figure 8, the shocked [S ii]/Hα ratios that best match the properties of the MOSDEF galaxies are offset to higher values than the SDSS data. If shocks are the cause of the offset in the N2-BPT diagram, then one might also expect an offset in the S2-BPT diagram as well. However, there is no measured offset between the SDSS and the z ∼ 2 data in the S2-BPT diagram (Shapley et al. 2015). We have two possible explanations: the electron density dependence on the shocked [S ii]/Hα ratio and contribution from diffuse ionized gas.

The electron density of the MOSDEF sample is ${290}_{-169}^{+88}$ cm⁻³, and the shocked line ratios for that particular density would lie between the 100 and 1000 cm⁻³ shock models. The SDSS galaxies also lie between the 100 and 1000 cm⁻³ shock models (see Figure 8). It is possible that the shocked [S ii]/Hα ratio for and electron density of 290 cm⁻³ lies close to the SDSS distribution. If this is the case, including the shocked emission would not change the [S ii]/Hα ratios of the z ∼ 0 galaxies much. The small difference between the shocked [S ii]/Hα ratio and the photoionized [S ii]/Hα ratio could explain the lack of an offset in the S2-BPT diagram.

Another explanation for no offset in the S2-BPT diagram could be because of less contribution from diffuse ionized gas at z ∼ 2 compared to z ∼ 0. The fraction of Hα emission from diffuse ionized gas decreases as Σ_Hα increases (Oey et al. 2007). As emission from diffuse ionized gas decreases, the [S ii]/Hα ratio decreases while [O iii]/Hβ stays the same (Sanders et al. 2017; Zhang et al. 2017). Since galaxies at z ∼ 2 have higher SFR (Madau et al. 1998; Reddy & Steidel 2009; Madau & Dickinson 2014) and are more compact (e.g., Shen et al. 2003; Barden et al. 2005; Trujillo & Pohlen 2005), they have higher Σ_Hα, which implies they will have less contribution from diffuse ionized gas if these local relations hold at z ∼ 2. Local galaxies with high Σ_SFR do lie at lower [S ii]/Hα at a given [O iii]/Hβ than those with low Σ_SFR on the S2-BPT diagram (Masters et al. 2016). If z ∼ 2 galaxies follow these same trends, then we should expect a lower [S ii]/Hα ratio at a given [O iii]/Hβ ratio. Therefore, the lack of an offset in the S2-BPT diagram could be because less contribution from diffuse ionized gas causes the narrow emission to lie at lower [S ii]/Hα than average z ∼ 0 galaxies, while the broad emission is higher [S ii]/Hα due to shocks. In this scenario, these effects compete with each other and ultimately cancel each other out, resulting in no net offset in the S2-BPT diagram.

The presence of shocks can also affect measurements made from single emission lines, such as SFR from Hα. If the broad emission should be removed when calculating SFR, then not doing so would overpredict the SFR. Given the measured BFRs, SFRs would be overpredicted by 15%, 40%, and 68%, respectively, in our low, medium, and high stellar mass stacks. However, given the large number of systematic uncertainties in these calculations (extinction curves, nebular versus continuum extinction, IMFs, and star formation histories), a 15% to 70% offset could go undetected. Furthermore, despite only detecting a BFR of 0.15 in the lowest mass bin, the contribution from broad emission may be higher because of our inability to detect broad emission with FWHM <300 km s⁻¹.

Reduced 2D spectra have two dimensions: spatial and wavelength. The overall goal of BMEP is to optimally extract 1D integrated spectra, which sums the flux over the spatial dimension and leaves the user with flux versus wavelength. The first step in extracting spectra is to find the primary object. BMEP draws a line over the object's expected position, which helps distinguish the primary object from serendipitous objects. Next, the user interactively bins the data in the wavelength direction to create a flux profile. Finally, the user fits a Gaussian to the profile. The center and width of this Gaussian determine the spatial region to sum, as well as the weighting profile for an optimal extraction (Horne 1986). In some objects with high S/N, the profile is non-Gaussian, and the user can choose to weight by the actual profile instead of the Gaussian fit.

Although the above process sounds simple, it is difficult to determine which wavelengths to bin to create the highest S/N spatial profile. Many galaxies at high redshift have bright emission lines compared to their continuum (Stark et al. 2013). Summing all wavelengths results in a spatial profile dominated by noise, which makes finding the center and width of the object impossible. The benefit of using BMEP is the ability to interactively create the best profile by adding or removing wavelengths when creating the spatial profile. Additionally, some galaxies do have continuum, but summing all wavelengths results in a noisy spatial profile because skylines would be included. The user can enable a "continuum mode," which does not include skylines in the spatial profile by removing wavelengths where the variance is higher than the median variance.

After extraction, the spectrum is plotted and can be inspected. The locations where the user clicked are also drawn in red on the plot. In noisy spectra, this allows the user to easily find emission lines in the 1D spectra. Once an emission line is found, the user can fit the line to a Gaussian, input which line it is, and calculate a redshift. All emission lines and calculated redshifts are saved in a catalog. A separate program consolidates all the lines fit for each object and calculates the most likely redshift.

In cases where an object had no obvious emission lines or continuum in the 2D spectrum, a "blind" extraction was performed. For objects with no signal in any bands, the blind extraction uses the expected position of the object calculated from the mask file and uses the same extraction width as the star's width in each filter. For objects with signal in one or more bands, the blind extraction used the average extraction widths and centers from filters where there was signal from the object. The widths from each filter are corrected for seeing differences. These blind spectra allow us to put upper limits on emission lines for spectra if we know the redshift. Currently, BMEP is only able to read in MOSFIRE mask files for the blind extraction and would need to be modified to be able to read in mask files of a different format.

The optimal extraction used in BMEP is based on Horne (1986). While testing the software, we compared extractions of a bright object done by several different users. Some extractions differed in extracted flux by 2% to 4% at all wavelengths for some users. We traced the cause of this to rounding differences between two extraction profiles. The extraction width is determined by a Gaussian fit to the profile and in some cases, the extracted widths would be 1 pixel different because we rounded the extraction to the nearest whole pixel. The optimal extraction of Horne (1986) sums over an integer number of pixels in the spatial direction to calculate the flux at each wavelength. To fix this, we created a sub-pixel optimal extraction algorithm. This algorithm extracts the central pixels exactly the same as in Horne (1986) but adds a fraction of a pixel at each end.

We base the sub-pixel algorithm on Equation (8) from Horne (1986). However, it is simplified for MOSFIRE reduction because there is no sky subtraction or cosmic ray rejection needed as these are done during the 2D reduction. The equations from Horne (1986) with these simplifications are

$$\begin{eqnarray*}&&{x}_{b}^{{\prime} }={\boldsymbol{R}}(c-w)\end{eqnarray*}$$

$$\begin{eqnarray*}&&{x}_{t}^{{\prime} }={\boldsymbol{R}}(c+w)\end{eqnarray*}$$

$$\begin{eqnarray}&&\displaystyle \sum _{{x}_{b}^{{\prime} }}^{{x}_{t}^{{\prime} }}D={F}_{\mathrm{box}}^{{\prime} }\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&\displaystyle \sum _{{x}_{b}^{{\prime} }}^{{x}_{t}^{{\prime} }}V={V}_{\mathrm{box}}^{{\prime} }\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&\displaystyle \frac{{\displaystyle \sum }_{{x}_{b}^{{\prime} }}^{{x}_{t}^{{\prime} }}{PD}/V}{{\displaystyle \sum }_{{x}_{b}^{{\prime} }}^{{x}_{t}^{{\prime} }}{P}^{2}/V}={F}_{\mathrm{opt}}^{{\prime} }\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&\displaystyle \frac{{\displaystyle \sum }_{{x}_{b}^{{\prime} }}^{{x}_{t}^{{\prime} }}P}{{\displaystyle \sum }_{{x}_{b}^{{\prime} }}^{{x}_{t}^{{\prime} }}{P}^{2}/V}={V}_{\mathrm{opt}}^{{\prime} }.\end{eqnarray} \tag{ 6 }$$

Unnumbered equations are defining relationships or variables. Bold letters indicate functions. R is the round function, c is the center of the object from the Gaussian fit to the profile, w is half the width to extract, D is the the flux in one pixel, V is the variance in one pixel, P is the weighting profile, x_b is the pixel at the bottom of the profile, and x_t is the pixel at the top of the profile. The weighting profile comes from the Gaussian fit to the spatial profile. ${F}_{\mathrm{box}}^{{\prime} }$ is the boxcar flux, ${V}_{\mathrm{box}}^{{\prime} }$ is the boxcar variance, ${F}_{\mathrm{opt}}^{{\prime} }$ is the optimal flux, and ${V}_{\mathrm{opt}}^{{\prime} }$ is the optimal variance for the Horne (1986) algorithm that does not include sub-pixel corrections. We remove the wavelength subscript for simplification.

We extend this equation to extract a fraction of each pixel. The central region of extraction is extracted the same as in Horne (1986), then a fraction of the outer pixels are added to this flux. First we calculate the range where the flux is extracted in the same manner as Equations (3)–(6). This is between ${xb}^{\prime}$ and ${xt}^{\prime}$, which are calculated as follows:

$$\begin{eqnarray*}&&{x}_{b}={\boldsymbol{L}}(c-w+1)\end{eqnarray*}$$

$$\begin{eqnarray*}&&{x}_{t}={\boldsymbol{L}}(c+w).\end{eqnarray*}$$

L is the "Floor" function. Next, we calculate the "remainder" from the bottom (R_b) and the top (R_t):

$$\begin{eqnarray*}&&{R}_{b}=1-({x}_{b}-{x}_{b}^{{\prime} })\end{eqnarray*}$$

$$\begin{eqnarray*}&&{R}_{t}={x}_{t}-{x}_{t}^{{\prime} }.\end{eqnarray*}$$

Now we calculate the weighting for sub-pixel extraction on the edges:

$$\begin{eqnarray*}&&{P}_{t}=P({xt}+1){R}_{t}\end{eqnarray*}$$

$$\begin{eqnarray*}&&{P}_{b}=P({xb}-1){R}_{b}.\end{eqnarray*}$$

The boxcar extraction for the sub-pixel algorithm is:

$$\begin{eqnarray*}&&B=D({xb}-1){R}_{b}\end{eqnarray*}$$

$$\begin{eqnarray*}&&T=D({xt}+1){R}_{b}\end{eqnarray*}$$

$$\begin{eqnarray}&&\displaystyle \sum _{{x}_{b}}^{{x}_{t}}D+B+T={F}_{\mathrm{box}}\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&\displaystyle \sum _{{x}_{b}}^{{x}_{t}}V+{V}_{b}{R}_{b}+{V}_{t}{R}_{t}={V}_{\mathrm{box}},\end{eqnarray} \tag{ 8 }$$

where B and T are the flux from the bottom and top pixels to be added to the central region. For the optimal extraction, this extra flux is

$$\begin{eqnarray*}&&B=\displaystyle \frac{{P}_{b}D({x}_{b}-1){R}_{b}}{V({x}_{b}-1)}\end{eqnarray*}$$

$$\begin{eqnarray*}&&T=\displaystyle \frac{{P}_{t}D({x}_{t}+1){R}_{t}}{V({x}_{t}+1))}.\end{eqnarray*}$$

Calculate sub-pixel flux and variance:

$$\begin{eqnarray}&&\displaystyle \frac{\left({\displaystyle \sum }_{{x}_{b}}^{{x}_{t}}{PD}/V\right)+B+T}{{\displaystyle \sum }_{{x}_{b}}^{{x}_{t}}{P}^{2}/V+{P}_{b}^{2}/{V}_{b}+{P}_{t}^{2}/{V}_{t}}={F}_{\mathrm{opt}}\end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\left({\displaystyle \sum }_{{x}_{b}}^{{x}_{t}}P\right)+{P}_{b}+{P}_{t}}{{\displaystyle \sum }_{{x}_{b}}^{{x}_{t}}{P}^{2}/V+{P}_{b}^{2}/{V}_{b}+{P}_{t}^{2}/{V}_{t}}={V}_{\mathrm{opt}}.\end{eqnarray} \tag{ 10 }$$

If the spatial range to extract are integers, then R_b = 1, R_t = 0, ${x}_{b}^{{\prime} }={x}_{b}+1$, and ${x}_{t}^{{\prime} }={x}_{t}$, and one can recover the original equations from Horne (1986).

Figure 11 shows a comparison of the sub-pixel optimal versus normal optimal. This figure was produced by selecting a flat, featureless section of a star that has no skylines. Within this region, we calculated the average flux and S/N, and then we varied the extracted width. As one might expect, the Horne extraction has steps where the width rounds to the next pixel and the sub-pixel extraction is smooth. Though the jumps in flux are severe when the width is small, the steps flatten out as width increases. At the width where we extract (2× FWHM), the jumps between steps are quite small, at worst around 4%. However, because we use a standard star to calculate the absolute flux, this can cause every flux for a mask to be 4% different when two different people extract a mask only, due to rounding.

Zoom In Zoom Out Reset image size