The MOSDEF Survey: [S iii] as a New Probe of Evolving Interstellar Medium Conditions^*

The rest-optical emission lines of star-forming galaxies provide valuable insight into the physical properties of the ionized gas in H ii  regions. Diagnostics of nebular metallicity (O/H), ionization parameter (U), and electron density have been calibrated at z ∼ 0. An additional important consideration is the shape of the stellar spectrum ionizing the gas, primarily determined by the Fe/H of massive stars. Due to tight relations between these properties, H ii  regions and local star-forming galaxies follow excitation sequences in emission-line ratio diagrams, including the [O iii]/Hβ versus [N ii]/Hα and [S ii]/Hα Baldwin, Phillips & Terlevich (BPT) diagrams (Baldwin et al. 1981; Veilleux & Osterbrock 1987). Recent large near-infrared spectroscopic surveys at z ∼ 2 have demonstrated that high-redshift star-forming galaxies follow a sequence that is systematically offset from the z ∼ 0 sequence in the [N ii] BPT diagram (e.g., Steidel et al. 2014; Sanders et al. 2016; Kashino et al. 2019; Shapley et al. 2019). This shift in the excitation sequence implies that the ionized interstellar medium (ISM) conditions are changing with redshift, and in particular that at least some of the relations between O/H, U, and stellar Fe/H evolve.

Determining which parameters evolve at fixed O/H and by how much has proven difficult. Such work often relies on photoionization models to understand how each line ratio changes when the relevant ISM conditions are varied (e.g., Kewley et al. 2013; Steidel et al. 2014; Sanders et al. 2016; Strom et al. 2017). A particular challenge lies in characterizing the ionization parameter and hardness of the ionizing spectrum, as neither is simply tied to an observable and the two are highly degenerate when O/H is unknown. Useful constraints cannot be obtained when only a few rest-optical lines are available (the case for most high-redshift galaxies), and degeneracies and large uncertainties remain even when five to six strong rest-optical lines are available (Strom et al. 2018). Only one "pure" ionization parameter diagnostic is available among the strongest optical lines ([O iii]/[O ii]) and it is strongly affected by dust reddening. A better understanding of evolving ISM conditions can be obtained by introducing additional line ratios as constraints that move beyond the strongest optical lines ([O ii], Hβ, [O iii], Hα, [N ii], [S ii]), most preferably an additional ionization parameter diagnostic ratio of a single element that is relatively unaffected by reddening.

In this Letter we present the first measurements of [S iii]λλ9069,9531 for a representative sample at z > 1 and explore the utility of the S₃₂ ≡ [S iii]λλ9069,9531/[S ii]λλ6716,6731 ratio for constraining evolving ISM conditions. This analysis uses observations of star-forming galaxies at z ∼ 1.5 from the MOSFIRE Deep Evolution Field survey (MOSDEF; Kriek et al. 2015) in combination with galaxy and H ii  region samples from the local universe. In Section 2, we present the [S iii] detections for individual high-redshift galaxies and composite spectra, and describe the low-redshift comparison samples. We investigate the evolution of S₃₂ at fixed stellar mass (M_*) in Section 3. Finally, we present excitation sequences of [S iii] and [S ii] ratios in Section 4 and discuss implications for the evolving physical conditions in H ii  regions. Emission-line wavelengths are given in air. We assume a ΛCDM cosmology with H₀ = 70 km s⁻¹ Mpc⁻¹, Ω_m = 0.3, and Ω_Λ = 0.7.

2.1. The z ∼ 1.5 MOSDEF Sample

We draw a sample of high-redshift galaxies with [S iii] measurements from the MOSDEF survey, a four-year program that obtained rest-frame optical spectra of ∼1500 galaxies at 1.37 ≤ z ≤ 3.80. A detailed description of the survey and data reduction can be found in Kriek et al. (2015). The MOSDEF flux calibration is robust between filters, such that ratios of lines from different filters (e.g., [S iii]/[S ii]) are biased less than 0.05 dex on average with a scatter of 0.07 dex (Kriek et al. 2015). We utilize emission-line measurements, stellar masses, reddening estimates, and star formation rates (SFRs) from the MOSDEF catalogs. Stellar masses and continuum reddening are estimated from emission-line corrected photometry via spectral energy distribution (SED) fitting using the code FAST (Kriek et al. 2009), assuming constant star formation histories, solar metallicity, a Chabrier (2003) initial mass function (IMF), and the Calzetti et al. (2000) attenuation curve. Nebular reddening is estimated using the Hα/Hβ ratio when both lines are detected at signal-to-noise ratio (S/N) ≥ 3, assuming an intrinsic ratio of 2.86 and the Cardelli et al. (1989) extinction curve. When Hβ is not detected, we infer E(B–V)_gas from continuum reddening under the assumption that E(B–V)_gas ≈ E(B–V)_stars, as found to be true on average at z ∼ 2 (Kashino et al. 2013; Reddy et al. 2015). SFRs are derived from reddening-corrected Hα luminosity using the Hao et al. (2011) calibration converted to a Chabrier (2003) IMF.

Galaxies were targeted in three redshift bins: 1.37 ≤ z ≤ 1.70, 2.09 ≤ z ≤ 2.61, and 2.95 ≤ z ≤ 3.80. Masks in each redshift bin were observed in multiple near-infrared filters in which the [O ii], Hβ, [O iii], Hα, [N ii], and [S ii] emission lines fall. Accordingly, z ∼ 1.5 masks were observed in Y, J, and H bands; z ∼ 2.3 masks in J, H, and K; and z ∼ 3.4 masks in H and K only. [S iii] lies significantly redward of [S ii], falling in the K-band at z ∼ 1.5 and redshifted out of the bands of atmospheric transmission at z > 2. Thus, most MOSDEF z ∼ 1.5 targets do not have [S iii] measurements because observations in the K-band are lacking. However, there are 49 MOSDEF star-forming galaxies at 1.25 ≤ z ≤ 1.65 with K-band observations because they were either filler targets (26) or serendipitous detections (23) on z > 2 masks. At least one [S iii] line is detected at S/N ≥ 3 for 10 individual galaxies in this sample, with both lines detected in four. The spectra of these 10 objects are presented in Figure 1. In addition to individual detections and limits, we stacked spectra to obtain average [S iii] measurements for the sample and include non-detections. We created multiple composite spectra following the procedure of Sanders et al. (2018), requiring Hα S/N ≥ 3 as well as the additional criteria described in Table 1 for each. In brief, individual spectra were converted to luminosity density, normalized by Hα luminosity, shifted onto a common rest-wavelength grid, and median stacked without weighting. This process ensures that high-SFR galaxies do not dominate the composite line ratios. The Stack1 composite spectrum with coverage of both [S iii] lines is displayed in the top row of Figure 1.

Zoom In Zoom Out Reset image size

Table 1.  Description of the z ∼ 1.5 [S iii] Composite Spectra

Name	Selection Criteria	Ngala	$\langle \mathrm{log}(\tfrac{{M}_{* }}{{M}_{\odot }})\rangle$ b
Stack1	[S iii]λ9069 and [S iii]λ9531 coverage	28	9.78
Stack2	[S iii]λ9069 coverage and two bins in M*	35, 11c	9.54, 10.21c
Stack3	[S iii]λ9069, [O iii]λ5007, and Hβ coverage	29	9.85

Notes.

^aThe number of galaxies included in each composite. ^bThe median stellar mass of the galaxies in each composite. ^cThe low-mass bin is listed first, followed by the high-mass bin. The two mass bins, split at 10^9.94 M_⊙, were populated such that equivalent S/N is obtained on [S iii]λ9069 in both composite spectra.

Download table as:  ASCIITypeset image

2.2. The $z\sim 0$ Comparison Samples

2.3. The [S iii]λ9531/λ9069 Ratio

Not all of our targets are detected in both [S iii] lines. Accordingly, we need to convert the line flux of one line to the total doublet flux. This approach is possible because the ratio [S iii]λ9531/λ9069 is fixed to a value of ≈2.5 according to the transition probabilities (Osterbrock & Ferland 2006; Tayal et al. 2019). We show the [S iii]λ9531/λ9069 ratios for individual z ∼ 1.5 galaxies, the Stack1 z ∼ 1.5 composite, and the z ∼ 0 comparison samples in Figure 2. We find that Stack1 and three out of four z ∼ 1.5 galaxies are consistent with the theoretically expected ratio within 1σ. The single offset z ∼ 1.5 galaxy displays possible skyline contamination on the red wing of [S iii]λ9531 (Figure 1, fifth row). The z = 0 CHAOS H ii  regions also match the expected value on average. The median ratio of the z ∼ 0 MaNGA sample is 2.12, significantly lower than the expected value. This offset may indicate that Paschen- absorption at 9546 Å has not been fully accounted for. To avoid biasing the total [S iii] fluxes low, we only use [S iii]λ9069 for the MaNGA sample. For all samples and composites, we infer the total [S iii]λλ9069,9531 flux assuming that [S iii]λ9531/λ9069 = 2.5 when only one [S iii] line is detected.

Zoom In Zoom Out Reset image size

We now investigate the evolution of the emission-line ratio S₃₂ ≡ [S iii]λλ9069,9531/[S ii]λλ6716,6731. We begin by characterizing the global galaxy properties of the z ∼ 1.5 [S iii] sample. In the top panel of Figure 3, we show SFR versus M_* for the z ∼ 1.5 MOSDEF and z ∼ 0 MaNGA samples. The z ∼ 1.5 galaxies with [S iii] detections (red) and those with [S iii] coverage but no detections (pink) fall on the mean z ∼ 1.5 M_*–SFR relation described by the full MOSDEF z ∼ 1.5 star-forming sample (cyan). The [S iii] subset has a lower average M_* (∼10^9.7 M_⊙) than the full MOSDEF sample (∼10^10.0 M_⊙), but is not significantly biased in SFR at fixed M_*. Both individual [S iii]-detected galaxies and stacks including non-detections are representative of the typical z ∼ 1.5 star-forming population.

Zoom In Zoom Out Reset image size

The bottom panel of Figure 3 displays S₃₂ versus M_*. Individual z ∼ 1.5 galaxies span a wide range of S₃₂, with the highest S₃₂ values occurring in the low-mass half of the sample. Blue squares show stacks of z ∼ 1.5 spectra in two M_* bins (Stack2). At fixed M_*, [S iii]-detected galaxies lie at higher S₃₂ than the stacks due to preferential detection of high-excitation galaxies with stronger [S iii]. We find that S₃₂ decreases with increasing M_* for both the z ∼ 1.5 and z ∼ 0 samples, indicative of lower ionization parameter at higher M_* and higher metallicity. The anticorrelation between S₃₂ and M_* displays a similar slope at z ∼ 1.5 and z ∼ 0, but the z ∼ 1.5 galaxies are offset 0.25 dex higher in S₃₂ at fixed M_*.

It has been shown that [S ii] is significantly enhanced in z ∼ 0 integrated galaxy spectra due to diffuse ionized gas (DIG) emission, which is expected to be negligible at high redshifts (Zhang et al. 2017; Sanders et al. 2017; Shapley et al. 2019). DIG contamination thus biases redshift evolution comparisons. We correct for DIG contamination of [S ii] in the MaNGA line ratios according to the prescription of Sanders et al. (2017), removing the DIG contribution and yielding the contribution from H ii  regions only. The effect of DIG on galaxy-integrated [S iii] emission has not yet been investigated. However, DIG contribution to [S iii] is not expected to be strong because DIG primarily enhances low-ionization species. After correcting the MaNGA sample for DIG contamination of [S ii], the relation between S₃₂ and M_* appears to be nearly the same at z ∼ 0 (black line) and z ∼ 1.5 (blue squares).

Many strong-line ratios (e.g., [O iii]/Hβ, [O iii]/[O ii]) display significant evolution toward higher excitation at fixed M_* even after accounting for z ∼ 0 DIG (Sanders et al. 2016, 2018), thought to primarily reflect evolution toward lower metallicity at fixed M_* with increasing redshift. We now examine whether the observed lack of significant S₃₂ evolution after DIG correction is consistent with this picture. Mass–metallicity relation studies at both z ∼ 0 and z ∼ 1.5 find Δlog(O/H) ∼ 0.6 dex over a decade in M_* at fixed redshift (Zahid et al. 2014; Kashino et al. 2017), while Figure 3 shows S₃₂ changes by ∼0.2 dex over the same mass interval. Accordingly, S₃₂ has a very weak dependence on metallicity (${\rm{\Delta }}\mathrm{log}({S}_{32})\approx 0.33\times {\rm{\Delta }}\mathrm{log}({\rm{O/H}})$). Indeed, S₃₂ displays a much weaker dependence on metallicity than [O iii]/[O ii] in photoionization models (Kewley et al. 2019). From z ∼ 0 to z ∼ 1.5, O/H decreases by ∼0.2 dex at fixed M_* corresponding to an expected increase in S₃₂ of only ∼0.07 dex. This difference is smaller than the Stack2 error bars. We thus find that the lack of significant S₃₂ evolution at fixed M_* after DIG correction is fully consistent with the evolution of the mass–metallicity relation.

We now turn to the evolution of excitation sequences in emission-line ratio diagrams and the changing the physical conditions of ionized gas in H ii  regions with redshift. In Figure 4, we show [O iii]/Hβ versus [S ii]/Hα (left column) and S₃₂ (right column). The top row displays empirical data sets. In the top-left panel, the z ∼ 1.5 sample displays larger [S ii]/Hα at fixed [O iii]/Hβ than z = 0 H ii  regions on average, but is offset toward lower [S ii]/Hα compared to the z ∼ 0 MaNGA galaxies. The severe offset between the H ii  regions and z ∼ 0 galaxies demonstrates the strong influence of DIG emission on [S ii] (Sanders et al. 2017; Shapley et al. 2019). Correcting [S ii] and [O iii] emission for DIG contamination following Sanders et al. (2017) yields the black line that closely matches the median sequence of the H ii  regions. In the top-right panel, the z ∼ 1.5 sample lies on the sequence described by z ∼ 0 MaNGA galaxies, but is offset from the z = 0 H ii  region sequence toward lower S₃₂ at fixed [O iii]/Hβ. Once again, the influence of DIG emission biases the comparison to low-redshift galaxies.

Zoom In Zoom Out Reset image size

To interpret these offsets, we employ the set of photoionization models described in Sanders et al. (2020) to identify the qualitative shift in these line ratios when varying H ii  region physical conditions. These models were run using Cloudy (Ferland et al. 2017) with BPASS v2.2.1 binary models (Stanway & Eldridge 2018) as the input radiation field, where the stellar metallicity (Z_* = Fe/H) is a proxy for the hardness of the ionizing spectrum. The grids span 0.05–1.5 Z_⊙ in nebular metallicity and log U = −2.5 to −3.5, and are color-coded by Z_*. Lines of constant ionization parameter are solid and lines of constant O/H are dotted. The models adopt an electron density of n_e = 250 cm⁻³, typical at z ∼ 2 (Sanders et al. 2016). Varying n_e between 25 cm⁻³ (typical at z ∼ 0) and 250 cm⁻³ changes [O iii]/Hβ and [S ii]/Hα by ≲0.05 dex and does not change S₃₂, and does not affect our results.

The grids fail to fully overlap the z ∼ 1.5 sample in the bottom-right panel of Figure 4. This failure is likely due to a known problem wherein photoionization models underproduce [S ii] relative to other lines (see Kewley et al. 2019, and references therein), preventing us from placing quantitative constraints on the gas properties. However, the qualitative direction of line ratio shifts should be robust even in the face of systematic [S ii] underestimation. We thus proceed by comparing the observed offsets between the z ∼ 0 and z ∼ 1.5 samples in Figure 4 with the qualitative shifts in model line ratios when varying ionization parameter and the hardness of the ionizing spectrum.

We first consider the scenario where the ionizing spectrum varies with redshift at fixed O/H (i.e., moving between colors at a fixed grid vertex), We find that a harder ionizing spectrum (lower Z_*) leads to higher [O iii]/Hβ and [S ii]/Hα while leaving S₃₂ unchanged. As a result, excitation sequences shift toward higher [S ii]/Hα and lower S₃₂ at fixed [O iii]/Hβ, in agreement with the observed offsets between z ∼ 1.5 galaxies and z = 0 H ii  regions or DIG-corrected z ∼ 0 MaNGA galaxies. We conclude that the shift in line-ratio excitation sequences between z ∼ 0 and z ∼ 1.5 is primarily driven by a harder ionizing spectrum at fixed nebular metallicity and does not require significant changes to U at fixed O/H.

If we instead consider varying the ionization parameter while keeping all other parameters fixed, increasing U (i.e., moving along dotted lines of a single color) leads to an increase in [O iii]/Hβ and S₃₂, and a decrease in [S ii]/Hα. The net effect is to shift galaxies toward lower [S ii]/Hα at fixed [O iii]/Hβ, and along the [O iii]/Hβ−S₃₂ sequence producing no significant offset in S₃₂ at fixed [O iii]/Hβ as lines of constant metallicity run roughly parallel to the full empirical sequences in the lower-right panel. Thus, larger U at fixed O/H can account for the offset (or lack thereof) between the z ∼ 1.5 sample and the z ∼ 0 MaNGA sample without DIG correction (dashed black lines in the lower panels). This conclusion was reached by past studies based on the position of z ∼ 1.6 star-forming galaxies in the [S ii] BPT diagram relative to a z ∼ 0 SDSS sample in which DIG was not accounted for (Kashino et al. 2017, 2019). However, because DIG emission is expected to be negligible in the highly star-forming compact galaxies at high redshifts (Shapley et al. 2019), a fair comparison is either to the H ii  region sample or to DIG-corrected integrated galaxy spectra, from which the z ∼ 1.5 galaxies are offset toward higher [S ii]/Hα and lower S₃₂ at fixed [O iii]/Hβ. Higher U at fixed O/H fails to account for these offsets.

Our results thus favor a harder ionizing spectrum at fixed O/H with increasing redshift, in agreement with recent work at z ∼ 2 based on deep rest-ultraviolet continuum spectroscopy (Steidel et al. 2016) and electron temperature metallicities (Sanders et al. 2020). However, the measurements utilized in this work can be acquired for many galaxies with a significantly smaller observational investment. Increasing the sample of z > 1 galaxies with [S iii] detections thus presents a viable path forward to an understanding of ISM conditions in individual high-redshift galaxies spanning a wide range in M_*, SFR, and metallicity. Fully leveraging new [S iii] observations requires more realistic photoionization models to turn qualitative conclusions into quantitative constraints. Photoionization modeling will be particularly discriminating in cases where both S₃₂ and [O iii]/[O ii] are available, where the use of two independent ionization parameter diagnostics simultaneously can break degeneracies between the ionizing spectrum and U. In the next few years, obtaining high-redshift data sets that expand beyond the strongest rest-optical emission lines should be a priority. Such observations are necessary to develop the tools required to interpret the wealth of information from the NIRSpec instrument on the James Webb Space Telescope, which will have sufficient sensitivity and wavelength coverage to measure [S iii] and numerous other weak lines up to z ∼ 5.5.

We acknowledge support from NSF AAG grants AST-1312780, 1312547, 1312764, and 1313171, grants AR-13907 and GO-15077 provided by NASA through the Space Telescope Science Institute, and grant NNX16AF54G from the NASA ADAP program. We also acknowledge a NASA contract supporting the WFIRST Extragalactic Potential Observations (EXPO) Science Investigation Team (15-WFIRST15-0004), administered by GSFC. We wish to extend special thanks to those of Hawaiʻian ancestry on whose sacred mountain we are privileged to be guests. Without their generous hospitality, the work presented herein would not have been possible.