CLEAR: The Ionization and Chemical-enrichment Properties of Galaxies at 1.1 < z < 2.3

As a low-redshift comparison sample, we make use of data from the SDSS Data Release 14 (DR14; Abolfathi et al. 2018), which includes emission-line fluxes and value-added catalogs. This catalog includes emission-line fluxes corrected for Balmer absorption and dust attenuation for SDSS III (including a reanalysis of galaxies from SDSS II; Thomas et al. 2013). We opt to use the stellar masses derived in the value-added catalog of Chen et al. (2012) using the Bruzual & Charlot (2003) stellar populations (and a Kroupa initial mass function, IMF) as these more closely match those derived for our CLEAR sample. For consistency in the comparison, we rederive the gas-phase oxygen abundances and ionization parameters of the SDSS galaxies using the same emission lines ([O II], [O III], Hβ) and method applied to the CLEAR sample (discussed below, Section 4.2).

The CLEAR HST/WFC3 pointings all lie within the CANDELS (Grogin et al. 2011; Koekemoer et al. 2011) GOODS-N and GOODS-S fields. The foundation of the CLEAR photometric catalog is the 3D-HST catalog from Skelton et al. (2014), which provides multiwavelength catalogs with photometric coverage from 0.3 to 8 μm. We have added to these HST F098M and/or F105W (i.e., Y-band) imaging as described in Estrada-Carpenter et al. (2019; see also R. Simons et al., in preparation). We then re-derived photometric redshifts and rest-frame colors (U − V and V − J) and derived stellar masses using an updated version of EAZY (Brammer et al. 2008; Brammer 2021). We refer to this catalog as 3D-HST+. We use the 3D-HST+ catalog for preliminary selection of the samples used here. We subsequently performed more sophisticated fits to the spectral energy distributions (SEDs) including both the 3D-HST+ catalog broadband photometry and WFC3 G102 and G141 grism spectra, and use the quantities derived from these latter fits for the analysis here (see Section 2.5).

The CLEAR program provides deep WFC3/G102 slitless spectroscopy in 12 pointings in the GOODS-N and GOODS-S fields. These data use observations of 10- or 12-orbit depth with WFC3/G102, which observe wavelengths 0.80–1.15 μm with R ∼ 210. We combined these data with all other available G102 data that overlap the CLEAR fields, including those data from programs GO-13420 (PI: Barro; see Barro et al. 2019), GO/DD-11359 (PI: O'Connell; see Straughn et al. 2011) and GO-13779 (PI: Malhotra; see Pirzkal et al. 2018). These data are described fully in a forthcoming paper (see Estrada-Carpenter et al. 2020; Simons et al. 2021, and R. Simons et al., in preparation).

We augment the CLEAR data with HST WFC3 slitless spectroscopy with the G141 grism from 3D-HST (Momcheva et al. 2016) that cover the CLEAR fields. The G141 data cover observed wavelengths 1.08–1.70 μm with R ∼ 130 and achieve flux limits for emission lines of 2.1 × 10⁻¹⁷ erg s⁻¹ cm⁻² (3σ for point sources; Momcheva et al. 2016).

We processed both the CLEAR and all ancillary WFC3 grism data in the CLEAR fields (see, Simons et al. 2021 and R. Simons et al., in preparation) using the grism redshift line and analysis software Grizli (Brammer 2022). The full process is described elsewhere (Simons et al. 2021, and see also Estrada-Carpenter et al. 2019, 2020; Matharu et al. 2022). In brief, we first reprocess the WFC3 G102 data, applying steps to correct for variable backgrounds and the flat-field, and we perform a sky-subtraction using the "Master Sky" provided in Brammer et al. (2015). We derive relative astrometric corrections to the processed data by aligning to the WFC3 F140W mosaic from Skelton et al. (2014).

We then use Grizli to model the G102 and G141 spectra of each object using the F105W and F140W direct images, respectively, and a coarse model fit to each galaxy's SED. We correct for galaxy contamination by subtracting the models for the spectra from nearby objects. This process is iterative. On the first pass, we model the spectra of all objects with m_F105W < 25 AB mag. We repeat the steps above. On the second pass, we apply a finer model correction to all objects with m_F105W < 24 AB mag. The adopted magnitudes for these steps are similar to those applied in the processing of the 3D-HST data (Brammer et al. 2012; Momcheva et al. 2016). Because the CLEAR G102 data are similar in depth to the 3D-HST G141 data, we achieve similar results (and visual inspection of the spectra and their residuals shows this accounts for the majority of contamination).

Finally, we use Grizli to extract two-dimensional (2D) and one-dimensional (1D) spectra for all galaxies in the 3D-HST+ catalog that fall in the CLEAR fields with brightness, m_F105W ≤ 25 AB mag, including the corrections for contamination described above. We used Grizli to measure spectroscopic redshifts, emission-line fluxes, and stellar-population parameters from the spectral fits to the continua and emission lines from the G102 and G141 data and the available multiwavelength broadband photometry from the 3D-HST+ catalog (see Section 2.2). For this process, Grizli uses a set of template basis functions derived from the Flexible Stellar Populations Synthesis models (FSPS; Conroy & Gunn 2010a) that include a range of stellar populations and nebular emission lines. Grizli integrates each model with the transmission functions of the broadband filters (including the system throughput of the telescope and detectors), and projects each stellar population model to match the G102 and G141 2D spectral "beams" using the observed direct image (F105W for G102; F140W for G141), matching the role angle (ORIENT) of HST and object morphology as closely as possible. This approach is required to model the unique morphological broadening of the spectral resolution (required for slitless spectroscopy). Grizli performs a nonnegative linear combination of the template spectra and determines a redshift through χ² minimization and a marginalization over redshift. Grizli fits emission-line fluxes using the best-fit redshift. It subtracts the continua (correcting for absorption features, e.g., from Balmer lines) using the best-fit stellar population model. The depth of the G102 data varies slightly in some of the fields (which contain different numbers of orbits from ancillary data), and the sensitivity depends somewhat on wavelength. Nevertheless, the bulk of the data (assuming the nominal 12-orbit depth) are sensitive to emission-line fluxes for point sources of ≈ 2 × 10¹⁷ erg s⁻¹ cm⁻² (3σ), comparable to the G141 data (R. Simons et al. 2022, in preparation).

Here, we use the CLEAR v3.0 catalogs, which are an internal team release. These include emission-line fluxes, spectroscopic redshifts, and other derived quantities and their respective uncertainties for 6048 objects from Grizli run on the combination of the G102 and G141 grism data and broadband photometry using the 3DHST+ catalogs. Of these galaxies, 4707 galaxies have coverage with both G102 and G141. These will be described fully in the forthcoming paper on the data release (R. Simons et al. 2022, in preparation) and have been discussed in other papers using these data (e.g., Estrada-Carpenter et al. 2019, 2020; Simons et al. 2021; Jung et al. 2022; Backhaus et al. 2022; Cleri et al. 2022; Matharu et al. 2022).

Here we use the CLEAR spectroscopy to study galaxies with coverage of strong emission lines that are tracers of the gas-phase oxygen abundance and nebular ionization, namely [O II] λ λ3727, 3729, Hβ λ4862, and [O III]λ λ4959, 5008 (e.g Maiolino et al. 2008; Sanders et al. 2015, 2021; Curti et al. 2017; Strom et al. 2018; Kewley et al. 2019a; Maiolino & Mannucci 2019; Henry et al. 2021, and many others; see references therein). Our G102 and G141 data cover all of these lines for galaxies at redshifts 1.1 < z < 2.3. In addition, for galaxies in the redshift range 1.1 < z < 1.6, our data include coverage of Hα λ6564 + [N II] λ6548, 6584 (which are blended at the grism data, see below). The spectra also provide coverage of [Ne III] λ3869, which is not detected in the majority of galaxies (but see Backhaus et al. 2022), but is observed in galaxy stacks (see below in Section 3).

We selected galaxies from CLEAR for this study using the following criteria:

• 1.  
  Spectroscopic redshift derived from the grism data in the range, 1.1 < z_grism < 2.3. This redshift range ensures that all three of the lines [O II], [O III], and Hβ are all contained by the G102 and/or G141 data.
• 2.  
  Detection of [O II], Hβ, or [O III] with signal-to-noise ratio (S/N) ≥ 3 in at least one line in the total (combined) 1D spectra.
• 3.  
  Galaxies are undetected in X-ray catalogs based on the Chandra X-ray catalogs for Chandra Deep Field–North (CDF)–N (Xue et al. 2011) and CDF–S (Luo et al. 2017); we rejected objects within 1'' of sources flagged as Type = AGN. This step excludes strong active galactic nuclei (AGNs), and removes 5% of sources in the GOODS-N and 6% of sources in the GOODS-S 3DHST+ parent catalogs. As an additional test, we checked if any additional objects in our sample are flagged as potential AGNs using the "mass-excitation" (MEx) diagnostic of Juneau et al. (2014), modified to account for redshift (Coil et al. 2015; Henry et al. 2021). This removed no additional objects (which we interpret as evidence that all candidate AGNs in our sample are identified as such in the ultra-deep CDF–N and CDF–S X-ray data).

In addition, we remind the reader that all galaxies have m_F105W ≤ 25 AB mag, as they are drawn from our CLEAR 3DHST+ catalog (see Section 2.3).

The selection produces a sample containing 196 galaxies. Figure 1 shows their redshift distribution. The redshifts span 1.1 < z < 2.3 with a median of 1.5. The distribution is highly peaked at the first redshift bin, with z ∼ 1.25. This is largely a result of galaxies in GOODS-N (GN), which makes up a larger number of sources (120 galaxies) in our sample compared to GOODS-S (GS; 76 galaxies). We consider two bins in redshift, each containing roughly 50% of the sample, with one bin defined with 1.1 < z < 1.5 (median z = 1.3) and the other with 1.5 < z < 2.3 (median z = 1.8). This allows us to test for redshift evolution in the properties of the sample.

Zoom In Zoom Out Reset image size

Figure 2 shows the stellar-mass–SFR distribution for our CLEAR sample of 1.1 < z < 2.3 galaxies using the selection criteria above, compared to the 3D–HST+ parent sample in the same redshift range. The stellar-mass distribution of the CLEAR sample is consistent with the 3D–HST+ sample when we restrict the stellar-mass range to $9.2\lt \mathrm{log}{M}_{* }/{M}_{\odot }\lt 10.2$ (where both samples are reasonably complete). The SFR distributions show that the CLEAR sample here is biased toward higher SFRs, by 0.18 dex (a factor of 1.5) compared to the 3D–HST+ parent sample. The bias can be explained as a result of the emission-line selection: we require galaxies to have S/N > 3 in Hβ, [O II], and/or [O III]. The CLEAR line-flux detection limit is 2 × 10⁻¹⁷ erg s⁻¹ cm⁻² (3σ), which for Hβ corresponds to SFR ≃ 3–7 M_⊙ yr⁻¹ (with no dust attenuation) at z = 1.5 − 2.0 (assuming the calibration of Kennicutt 1998 for a Chabrier IMF). Comparing this to Figure 2 we see that this effectively limits our study to objects with higher SFRs than the median. This bias in SFR is similar to other studies of emission-line selected studies of galaxies (cf., Shivaei et al. 2015; Sanders et al. 2018). We expect this bias to have only a minor impact on our results as previous studies have shown that a change in SFR of 1 dex corresponds to a change in metallicity of ≃ 0.3 dex (e.g., Henry et al. 2021). Based on this argument, the bias in SFR between the emission-line-selected sample and the parent sample, ${\rm{\Delta }}(\mathrm{log}\,\mathrm{SFR})\simeq 0.18$ dex (Figure 2), corresponds to ${\rm{\Delta }}(\mathrm{log}Z)=0.05$ dex.

Zoom In Zoom Out Reset image size

In Appendix A we also consider a subset of 87 galaxies from this sample with 1.2 < z < 1.5 for which Hα+[N II] are covered by the data. These galaxies have a median redshift z = 1.30. We use this subsample to test how incorporating additional lines impacts the constraints on the gas-phase metallicity and ionization (cf. Henry et al. 2021).

In what follows, we compare the galaxy emission-line properties (including derived quantities such as gas-phase metallicity and ionization parameter) to galaxy stellar population parameters, including stellar masses, SFRs, and sSFRs. To derive these latter quantities, we use a custom-designed method that fits stellar population synthesis models to the broadband photometry (from our 3D-HST+ catalog; see Section 2.2) and the WFC3 G102 and G141 1D spectra (see Section 2.3). The method is discussed in detail elsewhere (Estrada-Carpenter et al. 2020; V. Estrada-Carpenter et al. 2022, in preparation), and we summarize it here.

We use the FSPS models (Conroy & Gunn 2010b) with a Kroupa (2001) IMF. We fit a total of 23 parameters, including metallicity (of the stellar population, Z_*), age, dust attenuation (A_V , assuming the Calzetti et al. 2000 model), and a flexible star formation history (allowing for 10 bins of SFR dynamically spaced in time, following Leja et al. 2019). We also include eight additional nuisance parameters to allow offsets in the normalization/calibration between the spectra and the photometry, and to allow for correlated noise between spectral data points (see Estrada-Carpenter et al. 2020, and V. Estrada-Carpenter et al., in preparation for more details). We then fit to the broadband data and grism spectroscopy (for this modeling, we currently exclude regions of strong-line emission) using a Bayesian formalism with a nested sampling to predict the posteriors. We marginalize the posterior probability distribution functions to derive constrains on the stellar population parameters. We find that excluding the emission lines from the SED fitting causes a small bias in the SFRs (and the specific SFRs) of the galaxies in our sample. We compared the SED-derived SFRs for objects in our sample to SFRs estimated from dust-corrected Hβ emission measured in the grism data (assuming Case-B recombination and the calibration of Kennicutt & Evans 2012). For galaxies with SFRs ≲5 M_⊙ yr⁻¹, the Hβ-derived SFRs estimated are higher by about 0.25 dex (for galaxies with higher SFRs, the bias is negligible). Because we use SED-derived specific SFRs, we ensure we are not biased toward galaxies with strong emission lines only. Furthermore, this potential bias in SFR (and specific SFR) has a negligible impact on our conclusions related to the specific SFR, as this bias is smaller than the trends seen in the data and remains present if we replace the specific SFR with alternative measures (such as Hβ EW; see Reddy et al. 2018, and Sections 5.3 and 6.4).

Using this method, we fit the stellar population parameters for all of the galaxies in our sample. Here we focus on the stellar population constraints for stellar mass (M_*), SFR, and dust attenuation (A_V ) for the study here. We will present results derived from these parameters elsewhere (V. Estrada-Carpenter et al. 2022, in preparation). Figure 3 shows results from the fitting for two galaxies in our sample. The figure includes the best-fit SED (with parameters that maximize the likelihood) along with the broadband photometry and grism data. The figure also shows the marginalized posteriors for the three parameters above. We take the mode and highest density interval (HDI; Bailer-Jones et al. 2018) from the posteriors as the measurement and inter-68th percentile range (e.g., the inter–16th-to-84th percentile range) for each parameter, respectively. In what follows, we refer to these as "SED-derived" values, as they were derived from fitting models to the galaxy SEDs.

Zoom In Zoom Out Reset image size

To facilitate with the interpretation of the HST spectra, we constructed stacked spectra for galaxies in our sample in bins of stellar mass. We first divided the galaxies into subsamples of stellar mass, $9.0\lt \mathrm{log}{M}_{* }/{M}_{\odot }\lt 9.5$, $9.5\lt \mathrm{log}{M}_{* }/{M}_{\odot }\lt 10$, $10\lt \mathrm{log}{M}_{* }/{M}_{\odot }\lt 10.5$, and $10.5\lt \mathrm{log}{M}_{* }/{M}_{\odot }\lt 11.5$. To create the stacks, we corrected the spectra for dust attenuation assuming the Calzetti et al. (2000) model and the A(V) values derived from the SED fits (see Section 2.5). We then shifted all 1D spectra for the galaxies to the rest-frame using the measured redshift from the grism data. We linearly interpolated the data to a wavelength grid over 2500–7500 Å at a resolution of δ λ = 5 Å. We normalized each galaxy in the rest-wavelength range 4300–4500 Å (a window that avoids strong emission features, following Zahid et al. 2017) and coadded the spectra, weighting by the inverse variance of the flux density. We then divided the spectra by the total weights to obtain a mean spectrum. We also created a weighted sum of the variance of the flux density in the same way to study the variation of the spectra among galaxies in the sample.

Figure 5 shows the stacked spectra of the CLEAR galaxies in the bins of stellar mass. The spectra show common features, most prominently strong emission from [O II], [O III], Hβ, and Hα. In addition, weaker lines are also evident, including [Ne iii], Hγ, [S ii], and [O i]. The shaded region of the stacked spectra in the figure shows the scatter of the population in each stack (i.e., this is not the uncertainty on the mean). The shading indicates the scatter is generally larger at shorter wavelengths, which we attribute to variations in star formation histories (although some of this may be caused by the lower sensitivity of the G102 grism at bluer wavelengths).

Zoom In Zoom Out Reset image size

In general, the strength of the spectral features increases with decreasing stellar mass. In particular, it is in the lowest-mass galaxies ($9.0\lt \mathrm{log}{M}_{* }/{M}_{\odot }\lt 9.5$) where the strongest emission is seen, and where weaker lines such as [Ne iii] and Hγ become prominent. The relative strength of [O III]/[O II] also increases with decreasing stellar mass. This is an indication of increasing gas ionization and/or decreasing gas-phase oxygen abundance. We will explore this quantitatively below.

Because of this sizable variation in the spectral properties of the galaxies even at fixed stellar mass, we opt to study the individual spectra in most of the analysis that follows. However, we also use emission-line ratios and EW measurements from the stacks to interpret the average evolution of galaxies as a function of stellar mass and gas-phase metallicity. This complements work that analyzes the average properties of galaxies from stacked HST WFC3 grism spectra (including, e.g., Henry et al. 2021, which in part uses stacks that include the same data used here).

To measure the emission-line fluxes from stacked spectra in Figure 5, we adapted the penalized pixel-fitting method (pPXF; Cappellari 2017) for the CLEAR data. The primary difference between our use of pPXF and other data sets is that the CLEAR WFC3 grism data are lower resolution (R ∼ 100–200) than other spectroscopic studies (see, Cappellari 2017). Nevertheless, pPXF fits simultaneously the stellar components and nebular emission (i.e., correcting the nebular emission for stellar absorption), fitting for the line width (which is important here as our stacks include individual spectra with different spectral resolution owing to the morphological broadening). And, pPXF can separate the Balmer emission from the metal emission features. For our purposes, we added to pPXF the [Ne III] λ3868 emission line as this is prominent in our data (Figure 5). We then ran pPXF on the stacked data in Figure 5 (using stacks where the input spectra have been corrected for dust attenuation). We ran pPXF in two modes: one where each Balmer emission line is fit separately and one where we tie the Balmer emission lines to their theoretical Case-B values (e.g., Osterbrock 1989), and we use the latter for cases where Hγ is too weak to be visible in the spectra. We report the results in Table 1. Because the flux density in the stacked spectra have been normalized, we report emission-line flux ratios and EWs of prominent features. We use these results in the Section 6.4 to interpret the bulk trends between emission-line ratios and galaxy properties.

Table 1. Relative Emission-line Flux Ratios and Line Equivalent Widths from Stacked Spectra of CLEAR Galaxies at 1.1 < z < 2.3

	Number					EW Hβ	EW Hγ	EW Hc
Sample	of Galaxies	$\langle 12+\mathrm{log}({\rm{O}}/{\rm{H}})\rangle$ a	log R23	log O32	log [Ne III]/[O II]	(Å)	(Å)	(Å)
Stacked Spectra of Galaxies in bins of Stellar Mass
$\mathrm{log}{M}_{* }/{M}_{\odot }\gt 10.5$	9	8.83 ± 0.05	0.55 ± 0.06	−0.40 ± 0.12	<−1.33	21.9 ± 1.0	4.3 ± 1.1 b	2.5 ± 1.1 b
$\mathrm{log}{M}_{* }/{M}_{\odot }=[10,10.5)$	36	8.65 ± 0.03	0.70 ± 0.06	−0.15 ± 0.11	−1.18 ± 0.12	23.6 ± 1.2	9.7 ± 1.0 b	4.7 ± 1.1 b
$\mathrm{log}{M}_{* }/{M}_{\odot }=[9.5,10)$	84	8.55 ± 0.02	0.80 ± 0.03	0.07 ± 0.11	−1.03 ± 0.17	44.2 ± 1.9	16.3 ± 2.0	4.7 ± 1.7
$\mathrm{log}{M}_{* }/{M}_{\odot }=[9.0,9.5)$	66	8.42 ± 0.03	0.92 ± 0.03	0.27 ± 0.11	−0.87 ± 0.19	49.6 ± 3.8	18.5 ± 3.2	8.0 ± 3.5
Stacked spectra of galaxies with high- and low-ionization parameters, all with $\mathrm{log}{M}_{* }/{M}_{\odot }=[9.3,9.7]$ and $12+\mathrm{log}({\rm{O}}/{\rm{H}})\gt 8.3$
High ionization,	31	8.49 ± 0.04	0.89 ± 0.01	0.35 ± 0.11	−0.96 ± 0.11	61.3 ± 0.6	23.7 ± 0.5	8.2 ± 0.5
$\mathrm{log}q\gt 7.8$
Low ionization,	32	8.52 ± 0.04	0.86 ± 0.03	−0.03 ± 0.11	−1.12 ± 0.11	35.4 ± 0.5	13.0 ± 0.4	4.4 ± 0.5
$\mathrm{log}q\lt 7.8$

Notes. All line ratios and emission-line EWs are measured from the stacked spectra using pPXF as described in the text (Sections 3.2 and 6.4). EWs are measured in the rest frame.

^a Mean metallicity, and uncertainty on the mean, derived from the measurements of the individual galaxies in the sample (see the text). ^b Hγ and/or H weakly detected; flux ratios forced to their theoretical values in model fitting, Hγ/Hβ = 0.468 and H/Hβ = 0.159 (Osterbrock 1989). ^c Blended with [Ne iii] 3968.

Download table as:  ASCIITypeset image

Figure 6 compares the distribution of R₂₃ and O₃₂ for galaxies in the CLEAR sample to those from SDSS DR14 (Thomas et al. 2013). For both SDSS and CLEAR galaxies, the emission lines have been corrected for dust extinction. For SDSS, we used the values provided by Thomas et al. (2013). For our CLEAR galaxies, we corrected the line ratios using dust attenuation estimates from our SED fitting to the grism spectra and photometry (see Section 2.5) as most galaxies in our sample do not have measurements of multiple Balmer emission lines. We also assume the Calzetti et al. (2000) attenuation law and that the attenuation in the nebular gas is the same as for the stellar continuum (see Reddy et al. 2015).

Zoom In Zoom Out Reset image size

The CLEAR galaxies at 1.1 < z < 2.3 lie at the upper end of the R₂₃–O₃₂ distribution defined by SDSS (the latter is shown using a kernel density estimator, KDE). This is similar to the findings of other studies of high-redshift galaxies (e.g., Sanders et al. 2016; Strom et al. 2018; Runco et al. 2021), which interpret the data as an increase in ionization parameter, harder ionizing spectrum, decrease in metallicity, and possibly elevated nitrogen abundances and higher α/Fe abundance ratios. The CLEAR galaxies support many of these assertions, and we discuss these further below (see, Section 6.4). ¹⁶

Figure 6 also compares the line ratios to photoionization models. The models include the Dopita et al. (2013) models, which used the MAPPINGS IV code (bottom row, middle panel of the Figure). The Dopita et al. (2013) model includes updated atomic data, elemental abundance measurements, and modeling prescriptions. The input ionizing spectrum uses the STARBURST99 population synthesis model (Leitherer et al. 2014) for a stellar population with a constant SFR, observed at an age of 4 Myr, with a Salpeter IMF with an upper-mass cutoff of 120 M_⊙. The nebular region also assumes spherical geometry, isobaric photoionization, and that the distribution of electron velocities allows for an extended tail to higher energies (a so-called "κ" distribution). These models span the range of R₂₃–O₃₂ observed in the SDSS and CLEAR data, although the models are unable to produce the highest ratios (e.g., the models are limited to R₂₃ ≲ 1 while galaxies with R₂₃ > 1 are evident in the SDSS and CLEAR data).

The Kewley et al. (2013) models in Figure 6 (bottom row, left panel) show the effects of using the PEGASE 2 (Fioc & Rocca-Volmerange 1999) population synthetic models that include harder ionizing spectra (e.g., Kewley et al. 2013 add the spectra of planetary nebular nuclei for stars with high effective temperatures, T_e > 50,000 K to estimate for the effects of the stellar photospheres of massive stars). These models increase the ratios of O₃₂ in response to the increased ionization parameter. Nevertheless, these models also have difficulty achieving the highest R₂₃ values seen in the data. ¹⁷

Figure 6 also shows line ratios for the MAPPINGS V models (Kewley et al. 2019b, 2019a; bottom row, right panel). The MAPPINGS V models include updates with the latest atomic data and relative abundances (see, Nicholls et al. 2017). The input spectrum is based on the STARBURST99 stellar population synthesis models (as above) using models for massive stars (Hillier & Miller 1998; Pauldrach et al. 2001) that are able to better produce the ratios of blue/red supergiants in low-metallicity regions, such as the Magellanic Clouds.

Importantly, the MAPPINGS V models are isobaric, and consider the effects of different values for the ISM pressure, here defined as P/k = n_e T (in units of K cm⁻³, where n_e and T are the nebular electron density and temperature). In the present work, we adopt $\mathrm{log}P/k=6.5$ as this represents the expected conditions in high-redshift galaxies well (see, e.g., Acharyya et al. 2019). For example, Sanders et al. (2016) and Kaasinen et al. (2017) found evidence for higher median electron density, n_e ≃ 250–300 cm⁻³ (where n_e ≈ n_H for ionized gas) for ∼ 1–3 galaxies (see also, Runco et al. 2021). Combined with the expected nebular temperature of ∼10,000–20,000 K (see Andrews & Martini 2013; Sanders et al. 2020, and references therein), this implies a gas pressure of $\mathrm{log}P/k\sim 6.5\,-\,7.0$. Figure 6 shows the MAPPINGS V models assuming $\mathrm{log}P/k=6.5$, but we observe similar results for $6\lt \mathrm{log}P/k\lt 7.5$. Models with $\mathrm{log}P/k$ = 6.5 and 7.0 reproduce the span of the data as illustrated in Figure 6. Models with $\mathrm{log}P/k=6.0$ do not reproduce the data with the highest O₃₂ ratios, while models with higher pressure ($\mathrm{log}P/k$ = 7.5) produce lower R₂₃. We therefore adopt $\mathrm{log}P/k=6.5$ for our analysis while noting that changing this from 6.0–7.5 does not alter our conclusions.

The MAPPINGS V models still have difficulty producing the highest R₂₃ values seen in the data (i.e., those with $\mathrm{log}{R}_{23}\gtrsim 1$ in both SDSS and CLEAR; see Figure 6). This effect has been seen in other studies. To explain this offset could require stellar populations with enhanced α/Fe ratios (e.g., Sanders et al. 2016; Steidel et al. 2016; Strom et al. 2022), or a change in nebular geometry (e.g., "density" bounded nebula; Brinchmann et al. 2008; Nakajima & Ouchi 2014; Kashino & Inoue 2019; or clumpy geometries; Jin et al. 2022). This highlights the need for improvements in photoionization models to fully account for the range of line emission observed in high-redshift galaxies. We plan to investigate this in a future study.

We use the code, "Inferring the gas-phase metallicity (Z) and Ionization parameter" (IZI) developed by Blanc et al. (2015) to estimate the metallicity (Z, which we take to be the nebular oxygen abundance, 12 + log(O/H)) and ionization parameter (q). ¹⁸ IZI is a Bayesian code that computes posterior likelihoods for Z and q by comparing the measured emission-line fluxes for our galaxies to predictions from the photoionization models. IZI is flexible in the sense that it can use any combinations of lines, including the summed fluxes of multiple lines (which is useful in our case where emission lines are blended at the resolution of the WFC3 grisms).

We use MAPPINGS V models (Kewley et al. 2019b) assuming an isobaric pressure in the nebula of $\mathrm{log}P/k=6.5$ (K cm⁻³ for the reasons in Section 4.1). We adapted the output of the MAPPINGS V models into the format required for IZI. We assume a "flat" prior for both $\mathrm{log}Z$ and $\mathrm{log}q$, spanning the ranges $-1.3\leqslant \mathrm{log}Z/{Z}_{\odot }\leqslant 0.3$ and $6.5\lt \mathrm{log}q/(\mathrm{cm}\,{{\rm{s}}}^{-1})\lt 8.5$. We tested alternative pressure values from $\mathrm{log}P/k=6-7.5$ and find no substantial differences in our conclusions. We also tested the use of different priors (including a prior with the shape of the local- and high-redshift MZR; Maiolino et al. 2008; Andrews & Martini 2013; Henry et al. 2021; Sanders et al. 2021). These priors only change slightly the shape of the posteriors of lower-mass galaxies (shifting the inter-68th percentile range to higher values of metallicity by <0.1 dex) while leaving the median values nearly unchanged (at fixed stellar mass). Moreover, adopting a prior affects the MZR, changing the average metallicity of galaxies at $\mathrm{log}{M}_{* }/{M}_{\odot }$ = 9.5 by <0.01 dex (Section 5.1). We therefore adopt the results from the "flat" prior to avoid any bias inflicted by the prior information, but this has a minimal impact on our conclusions.

We focus on results that use the sum of the emission from [O II]λ λ3726, 3729 doublet, the sum of the emission from [O III]λ λ4959, 5008, and the emission from Hβ as these lines are available over our full redshift range of our sample, and they constrain both the ionization state of the gas and trace the majority of nebular oxygen. In all cases, we correct the line emission for dust attenuation as in Section 4.1. In Appendix A, we compare these results to the case where we also include the emission from Hα+[N II]λ λ6548, 6583 and [S II]λ λ6716, 6731 with the emission from [O II], [O III], and Hβ for galaxies at 1.1 < z ≲ 1.6 where all of these lines are also covered by the WFC3 grisms. Adding these lines does not change our interpretation, though it does provide additional confidence in the results derived from only the [O II], [O III], and Hβ lines.

We show examples of the results in Figures 7 and 8 for galaxies in GOODS-N and GOODS-S, respectively. For each galaxy in the figures, we show the image (ACS F775W, WFC3 F105W, and F160W) along with the 1D spectrum from the G102 + G141 grisms. The rightmost panels show the posteriors on gas-phase metallicity (12 + log O/H) and ionization $\mathrm{log}q$ (in units of centimeters per second) derived from IZI for each galaxy. Because the images show the same observed bands, color differences indicate the presence of the strong emission lines, spectral breaks, or spatially variant dust effects as they redshift through the filters. For example, for z ≳ 1.6 the 4000 Å/Balmer break shifts redward of the F105W band. This accounts for the redder appearance of some galaxies (such as GS 28878). In other cases, the [O III] emission line shifts into the F160W filter for z ≳ 1.8, which accounts for the redder appearance of other galaxies (such as GN 32485).

Zoom In Zoom Out Reset image size

In what follows, we report the mode for the metallicity and ionization using the $P(\mathrm{log}Z)$ (where we use $\mathrm{log}Z\,=12+\mathrm{logO}/{\rm{H}}$ here) and $P(\mathrm{log}q)$ distributions, along with the 16th and 84th percentiles to indicate the uncertainties. We derive the latter using the HDI, which is the smallest region that contains 68% of the probability density (see Bailer-Jones et al. 2018; Estrada-Carpenter et al. 2020). These are indicated by the vertical lines in the $P(\mathrm{log}Z)$ and $P(\mathrm{log}q)$ panels for each galaxy in Figures 7 and 8.

Zoom In Zoom Out Reset image size

Figure 9 shows our results for the metallicity (12 + logO/H) derived from the R₂₃ line ratio. The data points show the results for individual CLEAR galaxies with metallicities and ionization derived using IZI from the strong emission lines. The curves in Figure 9 show calibrations from the literature, derived using different methods and galaxy samples. It should be noted that the errors on the data points are highly covariant (as the metallicity values are derived from the R₂₃ ratios). To illustrate this, the right panel of Figure 9 shows the equivalent 1σ and 2σ error ellipses derived from the covariance between these parameters for a fiducial galaxy in the sample with $12+\mathrm{log}({\rm{O}}/{\rm{H}})=8.8$. This accounts for the lack of perceived scatter in the results (the scatter is covariant as indicated by the error ellipse, and is directed along the observed sequence between R₂₃ and metallicity).

Zoom In Zoom Out Reset image size

The R₂₃–metallicity calibrations include those based on nearby star-forming galaxies (e.g., SDSS galaxies at z ∼ 0.1; Maiolino et al. 2008; Curti et al. 2017). Maiolino et al. (2008) used a combination of direct measurements (for low-metallicity galaxies) with metallicities inferred from photoionization models for higher-metallicity galaxies in SDSS. These results show that R₂₃ is famously "double-valued" with a maximum and inflection point around log R₂₃ ≃ 1 at 12 + logO/H ≃ 8.2. On the high-metallicity branch of R₂₃, other calibrations find lower metallicity at fixed R₂₃ using direct T_e metallicity methods (Curti et al. 2017), but those authors caution that the [O III] λ4363 emission exhibits contamination from [Fe ii] emission in higher-metallicity regions. In addition, other studies have argued that some of the offset may result from a contribution to the emission from diffuse interstellar gas (Sanders et al. 2017), but at these redshifts, this effect is expected to be small (Sanders et al. 2021). Yet other studies find that the assumption about the ionization of the gas leads to biases that cause the metallicity derived from R₂₃ to be undervalued (Berg et al. 2021).

The effect of the (isobaric) pressure of the nebular gas is also important. Kewley et al. (2019b) used the suite of predictions from MAPPINGS V with a gas pressure of $\mathrm{log}P/k=5$, valid for $8.53\lt 12+\mathrm{log}({\rm{O}}/{\rm{H}})\lt 9.23$. These calibrations have a dependence on ionization parameter. These models produce the calibration illustrated in Figure 9 (thick gray line). These are consistent with the Curti et al. (2017) calibration. Increasing the pressure to $\mathrm{log}P/k=7$ has a strong impact on logR₂₃ for metallicities $12+\mathrm{log}({\rm{O}}/{\rm{H}})\gtrsim 8.5$ (see Kewley et al. 2019b, their Figure 9), increasing logR₂₃ by nearly ∼1 dex at $12+\mathrm{log}({\rm{O}}/{\rm{H}})\sim 9$. Because we use an isobaric pressure of $\mathrm{log}P/k=6.5$, this explains the offset in our calibration and that of Kewley et al. (2019b) and Curti et al. (2017), illustrated in Figure 9. Our calibration is more consistent with that derived independently by Strom et al. (2018).

We fit a quadratic function to the R₂₃–Z relation derived for the CLEAR galaxies of the form,

$$\begin{eqnarray}&&\mathrm{log}{R}_{23}=A\times {\left(\mathrm{log}Z-\mathrm{log}{Z}_{0}\right)}^{2}+\mathrm{log}{R}_{0},\end{eqnarray} \tag{ 3 }$$

where $\mathrm{log}Z\equiv \,12+\mathrm{log}({\rm{O}}/{\rm{H}})$ and A = −1.07 ± 0.03, $\mathrm{log}{R}_{0}\,=1.041\pm 0.004$, and $\mathrm{log}{Z}_{0}=8.228\pm 0.006$. We include the statistical uncertainties when performing the fit. However, we have not included the effects of the covariance between R₂₃ and Z (as discussed above), and therefore the uncertainties on these parameters are overestimated.

This relation we observe for CLEAR is consistent with other studies of high-redshift galaxies, which also show larger R₂₃ at fixed metallicity compared to calibrations derived for nearby galaxies. ¹⁹ For example, Strom et al. (2018) fit strong emission lines measured in z ∼ 2.3 galaxies using predictions from photoionization models that allow for a range of stellar and nebular metallicity, ionization parameter, and N/O abundance. They then parameterize R₂₃–Z as a quadratic expression, which is also shown in Figure 9. The Strom et al. (2018) relation is similar to the one for CLEAR, with an offset of 0.05–0.1 dex. The results from Strom et al. (2018) are similar to those from Maiolino et al. (2008) (although Strom et al. 2018 argue that the Maiolino et al. (2008) and other calibrations based on direct T_e measurements need to be revised upwards by 0.24 dex).

Figure 9 also shows that for the CLEAR sample, the ionization of the gas increases as the metallicity decreases. The ionization constraints (primarily from [O III]/[O II]) within the IZI fitting allow us to determine a likelihood that galaxies fall on the upper or lower branch of the R₂₃–Z relation. Galaxies on the lower-metallicity branch of R₂₃ have significantly higher ionization than galaxies on the upper-metallicity branch of R₂₃. Physically, the change in ionization causes an increase in the fraction of doubly ionized oxygen (O⁺⁺) at the expense of singly ionized oxygen (O⁺), which therefore leaves the numerator in the definition of R₂₃ mostly unchanged. However, the change in ionization should then be apparent in the ratio of O₃₂, which we discuss in the next subsection.

Figure 10 shows there is a tight correlation between the O₃₂ ratio and the ionization parameter, q, derived by modeling the emission lines of the CLEAR galaxies with the MAPPINGS V photoionization models (Section 4.2). The correlation is expected, as an increase in ionization parameter corresponds to an increase in the ratio of O⁺⁺ to O⁺. Low-redshift star-forming galaxies typically have ionization parameters in the range, $7.3\lt \mathrm{log}q/(\mathrm{cm}\,{{\rm{s}}}^{-1})\lt 7.6$ (Moustakas et al. 2010; Poetrodjojo et al. 2018), while H ii regions and super-star clusters in starburst galaxies (e.g., M82) have ionization parameters as high as $\mathrm{log}q/(\mathrm{cm}\,{{\rm{s}}}^{-1})\sim 8.2-8.7$ (Smith et al. 2006; Pérez-Montero 2014), which is observed in other low-redshift, extreme star-forming galaxies (e.g., Berg et al. 2021; Olivier et al. 2021). For the CLEAR galaxies at 1.1 < z < 2.3, the range of ionization parameter extends over $7.3\lesssim \mathrm{log}q/(\mathrm{cm}\,{{\rm{s}}}^{-1})\lesssim 8.5$, spanning the full range seen in star-forming galaxies in the local universe. This is consistent with other studies of high-redshift galaxies that show evidence for increased ionization (e.g., Sanders et al. 2016; Strom et al. 2018).

Zoom In Zoom Out Reset image size

We fit a linear relation between log O₃₂ and $\mathrm{log}q$ using

$$\begin{eqnarray}&&\mathrm{log}q=A\times {\mathrm{logO}}_{32}+\mathrm{log}{q}_{0}.\end{eqnarray} \tag{ 4 }$$

We fit the relation using a Gaussian mixture model (linmix; Kelly 2007), which yields A = 0.86 ± 0.07 and $\mathrm{log}{q}_{0}=7.53\,\pm 0.02$. This line is shown in Figure 10 (along with 400 random draws from the posterior). This linear fit is very similar to that from Strom et al. (2018), who derived their relation from a sample of z ∼ 2.3 star-forming galaxies with independent photoionization modeling.

Closer inspection of the CLEAR galaxies in Figure 10 shows that while [O III]/[O II] correlates with ionization parameter, there is a secondary effect. The effect is subtle (with galaxies with lower R₂₃ shifting above the relation at lower ionization and below the relation at higher ionization). This effect is likely a result of the fact that many of the galaxies in our sample lie at the peak of the R₂₃–Z relation (near R₂₃ ∼ 1 in Figure 9). R₂₃ correlates with gas-phase metallicity; this translates to a secondary effect in $12+\mathrm{log}({\rm{O}}/{\rm{H}})$. Empirically, this secondary dependence on R₂₃ comes from the relative strength of Hβ compared to [O II] + [O III]. Galaxies with stronger Hβ push toward lower R₂₃. This is predicted by the photoionization modeling (see Figure 7 of Kewley et al. 2019b), and physically is a result of the fact that gas cooling is more efficient in higher-metallicity environments. While the CLEAR galaxies in Figure 10 show an apparent ceiling, with logO₃₂ ≲ 0.75, this seems to be only a property of our sample. For example, some extreme galaxies at z ∼ 0.01 with low metallicity ($12+\mathrm{log}({\rm{O}}/{\rm{H}})\simeq 7.4\mbox{--}7.6$) have even higher line ratios (logO₃₂ ≃ 1.1–1.3) with high-ionization parameters ($\mathrm{log}U\simeq -2.4$ to −1.8, Berg et al. 2021). This is consistent with our relation (Equation (4)), which would predict $\mathrm{log}U\,\simeq -2.0$ to −1.8 for these galaxies. Taken together, this is evidence that while to first order the log O₃₂ ratio traces gas ionization strongly, there is a secondary dependence on metallicity that contributes to the scatter in this relation, and this persists at high redshift (1 ≲ z ≲ 2).

Figure 11 shows the relation between stellar-mass, and gas-phase metallicity (the "MZR") derived for the galaxies in CLEAR at 1.1 < z < 2.3 compared to some relations in the literature at low and high redshift. We also show results for SDSS DR14 using results derived from the same set of photoionization models and emission-line fluxes as for CLEAR (see Section 4.2). Figure 11 shows that the gas-phase metallicities for the SDSS galaxies derived for these models mostly follows those derived by other methods (notably, Tremonti et al. 2004 and Maiolino et al. 2008; the latter is illustrated in Figure 11). Comparing these results, there is evolution from the relation from SDSS (z ∼ 0.2) compared to CLEAR. Galaxies at higher redshift have lower metallicities, and this has been observed by multiple studies (using myriad methods to derive gas-phase abundances; e.g., Maiolino et al. 2008; Maiolino & Mannucci 2019; Henry et al. 2021; Sanders et al. 2021, and references therein).

Zoom In Zoom Out Reset image size

To measure the evolution in the MZR, we parameterize the relation using the prescription of Maiolino et al. (2008),

$$\begin{eqnarray}&&12+\mathrm{logO}/{\rm{H}}=-0.0864{\left(\mathrm{log}{M}_{* }-\mathrm{log}{M}_{0}\right)}^{2}+{K}_{0}.\end{eqnarray} \tag{ 5 }$$

Here, logM₀ is a scale stellar mass (in units of solar masses) when the relation achieves metallicity K₀. Figure 11 shows the relations derived by Maiolino et al. for SDSS and the AMAZE samples at z = 0.07, 0.7, and 2.2 (as labeled), using an independent strong-line calibration (Kewley & Dopita 2002) calibrated against SDSS DR4 observations.

We fit Equation (5) to the results for our SDSS and CLEAR samples. For both SDSS and CLEAR, we have modeled the same set of emission lines ([O III]λ λ4959+5007, [O II] λλ3726+3728, Hβ) with the same photoionization models (see Section 4.2). These are independent from the calibration used by others (e.g., Maiolino et al. 2008). For SDSS, because we have sufficient dynamical range in stellar mass, we fit for both M₀ and K₀. However, because the stellar-mass distribution of CLEAR is strongly focused on $\mathrm{log}{M}_{* }/{M}_{\odot }\sim 9.3\mbox{--}9.9$, we fix K₀ = 9.00 at the value derived from SDSS and equal to that obtained by Maiolino et al. (2008) at z = 2.2. For CLEAR, we also fit for M₀ for subsamples of galaxies split in redshift for 1.1 < z < 1.5 and 1.5 < z < 2.3.

It is worth noting that the MZR relation we derive for SDSS (Figure 11, solid gray line) agrees well with the relation derived by Maiolino et al. (2008; dashed gray line). This comparison is important because we have used a different photoionization model, and different choices of strong emission lines. Maiolino et al. (2008) calibrated their metallicities using photoionization models that assume lower n_e , which result in lower pressure. Nevertheless, the calibrations agree well for R₂₃ versus $12+\mathrm{log}({\rm{O}}/{\rm{H}})$ for high metallicities ($12+\mathrm{log}({\rm{O}}/{\rm{H}})\gt 8.8$). At lower metallicities, our calibration shifts to higher $12+\mathrm{log}({\rm{O}}/{\rm{H}})$ at fixed R₂₃ (see Figure 9). This accounts for the slight increase in the median of the SDSS distribution we observed around masses $\mathrm{log}{M}_{* }/{M}_{\odot }\sim 9$ compared to the fitted relation from Maiolino et al. (2008). The differences emphasize the importance of calibrating the relation between strong emission lines and metallicity in order to study the absolute evolution in the MZR. The comparison we measure here is differential (in that we use the same set of photoionization models for the galaxies at all redshifts), but this ignores possible evolution in the physical conditions in the galaxies (in which case one should use photoionization models whose physical properties—especially the galaxy density/pressure—also evolve with time). We plan to explore these effects in a future study.

Table 2 shows the derived values of M₀ (and K₀) for our fits to the SDSS and CLEAR samples. By fixing K₀, we see a steady increase in logM₀ with increasing redshift. This corresponds to a decrease in the typical gas-phase metallicity with increasing redshift (at fixed mass). Using the results of the analytic fits, at a stellar mass of $\mathrm{log}{M}_{* }/{M}_{\odot }=9.0$, we observe a decrease in 12 + logO/H of 0.3 dex from z ∼ 0.2 to z ∼ 1.35 and an additional decrease of 0.2 dex z ∼ 1.9. This implies that galaxies like a progenitor of the Milky Way (Papovich et al. 2015) had metallicities of $12+\mathrm{log}({\rm{O}}/{\rm{H}})\simeq 8.3$ at z = 1.9 and 8.5 at z = 1.4 (i.e., the Milky Way progenitor had roughly Z_gas ≈ 0.4–0.6 Z_⊙ at 9–10 Gyr in the past). This is consistent with direct measurements of abundances in stars of this age within the Milky Way (e.g., Bergemann et al. 2014), and implies we are building a coherent picture of the chemical enrichment of galaxies like our own.

In Figure 11, the lower two panels show the MZR with the CLEAR galaxies color-coded by sSFR and ionization parameter ($\mathrm{log}q$). There is an apparent dependence on sSFR, in that galaxies with higher sSFR have lower gas-phase metallicity (and because we show this as a function of specific SFR, this relation is at fixed mass by construction). The dependence of the MZR on SFR has been observed previously both at high and low redshift, and several studies have argued that the "fundamental" MZR–SFR relation is independent of redshift (e.g., Ellison et al. 2008; Mannucci et al. 2010; Andrews & Martini 2013; Sanders et al. 2018, 2021; Cresci et al. 2019; Curti et al. 2020; Henry et al. 2021).

Similarly, there is an apparent trend between the galaxy MZR and ionization. Figure 11 shows that galaxies in CLEAR with higher stellar mass are generally only found to have lower ionization ($\mathrm{log}q$): indeed, galaxies with the lowest ionization ($\mathrm{log}q/(\mathrm{cm}\,{{\rm{s}}}^{-1})\lesssim 7.3$) are only found with higher stellar masses ($\mathrm{log}{M}_{* }/{M}_{\odot }\gt 10.2$). Galaxies with the highest ionization ($\mathrm{log}q/(\mathrm{cm}\,{{\rm{s}}}^{-1})\gtrsim 8.1$) are generally only found at lower stellar masses, $\mathrm{log}{M}_{* }/{M}_{\odot }\lt 9.9$. There is also qualitative evidence that at fixed stellar mass, galaxies with higher metallicity have lower-ionization parameters. This is similar to the relation between SFR and the MZR discussed above, and implies that ionization, metallicity, and sSFR are intertwined.

Figure 12 shows the relation between stellar-mass and nebular ionization parameter (i.e., the "mass–ionization-relation," or MQR) derived for the galaxies in CLEAR at 1.1 < z < 2.3. In comparison, we also show the MQR for galaxies in SDSS analyzed using the same set of emission lines and photoionization models used for the analysis of the CLEAR galaxies (see Section 4.2). The MQR clearly evolves from z ∼ 0.2 (from SDSS) to z ∼ 1.3 and to z ∼ 1.8 (from CLEAR). At fixed stellar mass, galaxies have a higher ionization parameter at a higher redshift, where the effect is stronger for galaxies of lower stellar mass. This extends trends seen both at lower redshift and higher stellar masses (see also Kewley et al. 2015; Kaasinen et al. 2018; Strom et al. 2022).

Zoom In Zoom Out Reset image size

We parameterize the MQR using a simple quadratic relation inspired by the MZR (Maiolino et al. 2008),

$$\begin{eqnarray}&&\mathrm{log}q=0.0571\times {\left(\mathrm{log}{M}_{* }-\mathrm{log}{M}_{0}\right)}^{2}+\mathrm{log}{q}_{0}.\end{eqnarray} \tag{ 6 }$$

Here, logM₀ is a scale stellar mass (in units of M_⊙) when the relation achieves ionization logq₀, and q is measured in units of centimeters per second. For our SDSS sample, we fit for M₀ and $\mathrm{log}{q}_{0}$. For CLEAR, we fix logq₀ to the value derived for the SDSS sample ($\mathrm{log}{q}_{0}=7.282$) because the stellar-mass distribution of CLEAR peaks at $\mathrm{log}{M}_{* }/{M}_{\odot }\sim 9.3-9.9$. For CLEAR we also fit M₀ for galaxies in subsamples of redshift, 1.1 < z < 1.5 and 1.5 < z < 2.3.

Table 3 shows the best-fit values and their uncertainties for M₀ (and Q₀) for our fits to the SDSS and CLEAR samples. There is a steady increase in logM₀ with increasing redshift. This corresponds to an increase in the typical nebular ionization with increasing redshift (at fixed mass). Using these fits, at a stellar mass of $\mathrm{log}{M}_{* }/{M}_{\odot }=9.0$, we observe an increase in $\mathrm{log}q$ of 0.29 ± 0.05 dex from z ∼ 0.2 to z ∼ 1.35 and an additional increase of 0.14 ± 0.04 dex from z ∼ 1.3 to z ∼ 1.8.

At a stellar mass of $\mathrm{log}{M}_{* }/{M}_{\odot }=10.0$, the evolution is weaker, as we observe a decrease in ionization parameter of 0.19 ± 0.03 dex from z ∼ 0.2 to z ∼ 1.3. At higher stellar masses $\mathrm{log}{M}_{* }/{M}_{\odot }\gtrsim 10.5$, there is very little evidence that the MQR evolves from z ∼ 0.2 to z ∼ 1.1−2.3, but the sample is limited by smaller numbers of massive galaxies. For example, Kaasinen et al. (2018) found an increase in the ionization parameter of ≃0.3 dex for galaxies with $\mathrm{log}{M}_{* }/{M}_{\odot }\simeq 10.4-10.9$ from z ∼ 0.2 to 1.5, only slightly stronger than the trend we see here. Our results are also similar to the independent measurements derived by Strom et al. (2022) at z ∼ 2.3.

One potential source of concern in the interpretation of the MQR (in Figure 12) is if our sample is biased by sources with strong emission lines. Our samples are selected with m(F105W) < 25 mag, so sources with strong emission lines could be overrepresented (particularly near our magnitude limit). To test this scenario, we used the measurements of the EW for strong lines ([O II], Hβ, and [O III]) for sources in our sample (measured from their G102 and G141 spectra) that have redshifts that place these lines in the F105W passband. We then correct the F105W magnitude (using, e.g., Equation (2) of Papovich et al. 2001), and exclude any object with m(F105W) > 25 mag. This removes only nine sources (a loss of <5% of the sample). Excluding these nine sources, we refit the MQR using Equation (6), and find that logM₀ is reduced by 0.05 and 0.04 dex for CLEAR at z ∼ 1.3 and 1.8, respectively (this is less than the statistical uncertainty in Table 3). Therefore, the evolution in the MQR is not seriously impacted by the presence of strong emission lines in the sample.

The bottom panels of Figure 12 show the dependence of the MQR on sSFR and on the gas-phase metallicity. There is an overall trend of increasing metallicity with increasing stellar mass, which is a consequence of the MZR (see Section 5.1). There is no identifiable trend between metallicity with ionization at fixed mass: for example, galaxies with $\mathrm{log}{M}_{* }/{M}_{\odot }\,=9.4-\,9.8$ span a wide range of $12+\mathrm{log}({\rm{O}}/{\rm{H}})$ and $\mathrm{log}q$. Qualitatively, in this stellar-mass range, galaxies with the highest ionization parameters have lower metallicity (and vice versa), but this statement is limited by the size of our sample. We return to these points in Section 6.

Another interesting question is to what extent the change in ionization parameter is driven by changes in the SFR (or the changes in SFR at fixed stellar mass, which is the sSFR). Kaasinen et al. (2018) considered the correlation between ionization parameter and sSFR for SDSS galaxies and galaxies at z ∼ 1.5. They concluded that higher $\mathrm{log}q$ is directly linked to an increase in SFR (see also Brinchmann et al. 2008).

To investigate the relation between nebular ionization and sSFR, we focus on CLEAR galaxies in the stellar-mass range, $\mathrm{log}{M}_{* }/{M}_{\odot }=9.2-\,10.2$ (where the majority of our sample resides). Figure 13 shows the trend between ionization and sSFR for these galaxies in CLEAR. We fit a linear relation,

$$\begin{eqnarray}&&\mathrm{log}q=A\times (\mathrm{log}[\mathrm{sSFR}]+9.0)\,+\,\mathrm{log}{q}_{0},\end{eqnarray} \tag{ 7 }$$

where q is the ionization parameter (in units of centimeters per second), and the sSFR is estimated from the stellar population fits to the galaxy SEDs (see Section 2.5) in units of yr⁻¹. Table 4 provides the fitted values and their uncertainties for A and logq₀.

Zoom In Zoom Out Reset image size

Table 4. Fitted Parameters for the Linear Relation between Specific SFR and Nebular Ionization for CLEAR Galaxies Using Equation (7)

Sample	Slope, A	logq0	ra	p-value b
Full sample	0.50 ± 0.10	7.67 ± 0.02	0.29	8.2 × 10−5
$\mathrm{log}{M}_{* }/{M}_{\odot }$= [9.2, 10.2]	0.31 ± 0.10	7.70 ± 0.02	0.27	9.2 × 10−4
$\mathrm{log}{M}_{* }/{M}_{\odot }$= [9.2, 10.2]	0.27 ± 0.12	7.71 ± 0.02	0.21	0.013
& logsSFR/yr−1 > −9.5

Notes.

^a Pearson's correlation coefficient. ^b Two-side probability (p-value) of correlation occurring by chance (assuming x and y are drawn from independent Gaussian distributions).

Download table as:  ASCIITypeset image

In all cases, there is a significant correlation between nebular ionization and specific SFR. For fits to the full galaxy sample (not shown in the figure) and for galaxy subsamples, there is a significant correlation (see Table 4). The full galaxy sample gives a slope of A = 0.50 ± 0.10. We also estimate the significance using Pearson's correlation coefficient, which gives r = 0.29 (with a p-value of 8.2 × 10⁻⁵, and 3.9σ for a Gaussian distribution). Limiting the fit to the subsample of galaxies with stellar mass in the range $\mathrm{log}{M}_{* }/{M}_{\odot }$ = [9.2, 10.2] yields a weaker correlation coefficient (r = 0.27), yet still significant (p = 9.2 × 10⁻⁴). Restricting the fit to the subsample with $\mathrm{log}{M}_{* }/{M}_{\odot }$ = [9.2, 10.2] and log sSFR/ yr⁻¹ <−9.5) yields a correlation coefficient of only r = 0.21 (with a p-value of 0.013, about 2.5σ). While still significant, it is slightly weaker than the relation when considering the full sample, indicating that the range of galaxy stellar masses and/or SFR differences likely accounts at least partly for the correlation. Kaasinen et al. (2018) also observed a correlation between ionization parameter and sSFR with a very similar slope for measurements from the stacked spectra for z ∼ 1.5 galaxies at higher stellar masses. Motivated by these results, we explore this relation more below (in Section 6.4).

The evolution in the MZR (Figure 11) has been observed for star-forming galaxies previously, and has been used to constrain the evolution of metals and feedback effects in galaxies as a function of redshift and stellar mass (see, e.g., Maiolino & Mannucci 2019). The most recent measurements find that over the mass range $\mathrm{log}{M}_{* }/{M}_{\odot }=9.5-10$, the evolution of the gas-phase metallicity evolves by ΔlogO/H = 0.2–0.3 dex from z ∼ 2.3 to z ∼ 0.1 (e.g., Henry et al. 2021; Sanders et al. 2021). This is generally consistent with our findings in CLEAR, where we see that galaxies at z = 1.35 (1.90) have metallicity lower by ΔlogO/H = 0.25 (0.35) dex compared to SDSS galaxies at z ∼ 0.2 at stellar masses, $\mathrm{log}{M}_{* }/{M}_{\odot }=9.4-9.8$.

Some comparisons to other studies of the MZR are useful as they highlight systematics in the analyses. Sanders et al. (2021) studied the evolution of the MZR to z > 3 using Keck/MOSFIRE observations of ≃ 450 galaxies. They find that the low-mass slope of the MZR is consistent with no evolution, following ${\rm{\Delta }}(\mathrm{logO}/{\rm{H}})/{\rm{\Delta }}\mathrm{log}{M}_{* }\sim 0.30$ from z ∼ 0 to 3.3. However, the normalization of the MZR evolves strongly, with Δ(logO/H)/Δz ≃ −0.11 (at fixed M_*) with a small uncertainty (≃0.02 dex). They argue this is consistent with the idea that at fixed stellar mass, the galaxy gas fractions and metal removal efficiencies increase at higher redshift.

We find stronger evolution with redshift in the normalization in our CLEAR and SDSS samples, Δ(logO/H)/Δz =− 0.21 ± 0.02 dex from z = 0.2 to z = 1.3 and Δ(logO/H)/Δz = −0.34 ± 0.04 dex from z = 1.35 to z = 1.90. This is higher than that from Sanders et al. (2021) by ∼ 0.2 − 0.3 dex and is likely related to the use of strong-line metallicity calibrators. Sanders et al. (2021) took average metallicities derived from multiple strong-line indicators, while we derive metallicities from the same set of lines fit to the MAPPINGS V photoionization models. Sanders et al. (2021) also used models that allow for increasing α/Fe, which could account for some of the offset in the evolution. We take this difference as an estimate of the systematics in the strong-line calibrators (e.g., see Kewley & Ellison 2008, who show the normalization of different calibrators can vary by as much as 0.7 dex).

The study of Henry et al. (2021) is more similar to the analysis presented here. They used measurements from stacked HST/grism spectra of more than 1000 galaxies at z ∼ 1.3–2.3 (along with higher-quality spectra of ∼50 individual galaxies) to measure the evolution of the MZR and derived gas metallicities using strong lines ([O II], Hβ, [O III], Hα+[N II]) with the calibration from Curti et al. (2017). They derive O/H abundances at z ∼ 1.3–2.3 that are consistent with those from Sanders et al. (2018) at z ∼ 2.3, yielding $12+\mathrm{log}({\rm{O}}/{\rm{H}})=8.28\pm 0.02$ (${8.37}_{-0.02}^{+0.01}$) for $\mathrm{log}{M}_{* }/{M}_{\odot }\,=$ 9–9.5 (9.5–10). Comparing to the z ∼ 0.1 results from Curti et al. (2017), they find that the normalization of the MZR evolves by ΔlogO/H ≃ 0.3 dex at a fixed mass $\mathrm{log}{M}_{* }/{M}_{\odot }\,=$ 9–9.5. This yields an evolution of ${\rm{\Delta }}\mathrm{log}({\rm{O}}/{\rm{H}})/{\rm{\Delta }}z\,\simeq \,$−0.2 dex from z ∼ 0.1 to 1.8, consistent with the evolution we derive from our analysis of the CLEAR and SDSS samples. The fact that Henry et al. (2021) measured the same absolute gas-phase metallicity as Sanders et al. (2018) at z ∼ 2 but different evolution again indicates that systematics are important, primarily in the absolute normalization of the MZR, both at high and low redshifts.

In summary, our CLEAR results add to the evidence that the MZR evolves by Δ(logO/H)/Δz ∼ 0.2 dex from z = 0.2 to z = 1.3 and that this rate of evolution in the MZR may increase at higher z (we observe Δ(logO/H)/Δz = 0.3 from z = 1.35 to z = 1.90). The strength of our result is that we have used the same set of photoionization models to convert the strong emission-line ratios to constraints on the gas-phase metallicities and ionization parameters, both for galaxies in CLEAR (z ∼ 1–2.3) and SDSS (at z ∼ 0.2). However, this strength is also a weakness because using the same set of photoionization models assumes that the physical conditions in the low-redshift galaxies (SDSS) are the same as in the higher-redshift galaxies (CLEAR). For example, we use the MAPPINGS V models with higher pressure, P_e /k = n_e T_e = 10^6.5 cm⁻³ K, for the high-redshift galaxies given the expected temperatures and particle density of the H ii regions (T_e ∼ 10⁴ K and n_e ≃ 300 cm⁻³; see Sanders et al. 2016; Kaasinen et al. 2017; Strom et al. 2017; Runco et al. 2021), where the electron density is an order of magnitude larger than observations at z ∼ 0.1 (e.g., Sanders et al. 2016). As discussed in Kewley et al. (2019b), adopting a lower pressure for the SDSS galaxies would decrease the metallicity of high-Z galaxies at fixed R₂₃. This is evident in Figure 9, which shows the difference between our calibration and that from Kewley et al. (2019b; who used $\mathrm{log}{P}_{e}/k\simeq 5$), and in Appendix A, which shows the differences between our MZR and those from Henry et al. (2019; who used the calibration from Curti et al. 2017, which is similar to that of Kewley et al. 2019b). Adopting lower pressure for the SDSS galaxies would lower the metallicity of high-Z galaxies by ∼0.1 dex (see, also, Sanders et al. 2021). Understanding these potential sources of systematic bias is crucial to obtaining an accurate measurement of the redshift evolution in the MZR.

It is therefore remarkable that even in lieu of the systematic uncertainties, the MZR evolution in Figure 11 shows agreement between measurements from SDSS and CLEAR from different studies Therefore, while there is work needed to understand the impacts of the assumptions about the metallicity calibrations from the strong emission lines, there appears to be some consensus on the absolute evolution of the MZR.

Figure 12 shows evidence for evolution in the MQR for star-forming galaxies. We quantify the evolution in the MQR over the redshift range z ∼ 0.2 to z ∼ 2.3 using our analysis of the galaxies in CLEAR and SDSS. At fixed stellar mass, $\mathrm{log}{M}_{* }/{M}_{\odot }=9.6$, the differential evolution between the SDSS galaxies and CLEAR galaxies corresponds to ${\rm{\Delta }}\mathrm{log}q/{\rm{\Delta }}z\,\simeq 0.2$ dex from z ∼ 0.2 to z ∼ 1.9.

Previous studies have seen evidence for this evolution, primarily in terms of an increase at high redshifts in the strength of emission-line ratios that are sensitive to the ionization parameter. Sanders et al. (2018) measured the evolution of O₃₂ as a function of stellar mass and redshift, finding a 0.5 dex evolution in O₃₂ from z ∼ 2.3 to ∼0 for galaxies at $\mathrm{log}{M}_{* }/{M}_{\odot }\simeq \,9.5$. This corresponds to a change in ionization parameter of ${\rm{\Delta }}(\mathrm{log}q)\approx 0.4$ dex (using Equation (4)). Kaasinen et al. (2018) measured the ionization parameter from stacked spectra of massive galaxies ($\mathrm{log}{M}_{* }/{M}_{\odot }\simeq 10.4\mbox{--}11$) at z ∼ 1.5. They showed these exhibit an increase in the ionization parameter of ≃0.4 dex at fixed stellar mass from z ∼ 0.2 to 1.5. This is consistent with the evolution we measure in CLEAR to z ≃ 1.9 and to z ∼ 0.2 from SDSS.

Figure 12 (bottom-left panel) also shows that at fixed stellar mass, galaxies in CLEAR span a range of specific SFR from log sSFR ∼ − 9.3 to −8.7, which corresponds to about a factor of ∼4. Qualitatively, the figure shows that galaxies with higher sSFRs favor higher ionization parameters. This is reminiscent of the well-studied trend between the MZR and SFR in that galaxies with lower metallicity favor higher SFRs at fixed stellar mass reported in other studies (see Maiolino & Mannucci 2019; Sanders et al. 2021; Henry et al. 2021).

Is the MQR a consequence of the evolution in the MZR, or of the evolution of the SFR–MZR relation? To investigate this, Figure 12 (bottom-right panel) shows that at fixed stellar mass, higher-metallicity galaxies in CLEAR generally have lower-ionization parameters, but the scatter is large. Figure 12 also shows that higher ionization parameters correlate with lower metallicities. However, we argue the situation is more nuanced: in Section 6.4 we show that galaxies in CLEAR, when matched in stellar mass and metallicity, show a wide range of ionization parameter. Therefore, the metallicity is not solely the cause of the evolution in the MQR.

Figure 14 shows the relation between the ionization parameter ($\mathrm{log}q$) and gas-phase metallicity ($12+\mathrm{log}({\rm{O}}/{\rm{H}})$) for the CLEAR galaxies. The figure compares these results to those from MOSDEF (Topping et al. 2020) at z ∼ 2 and to measurements for individual H ii regions in galaxies (Pérez-Montero 2014; Espinosa-Ponce et al. 2022). The metallicity and ionization of the CLEAR galaxies follow the physical parameter space seen in these other samples. Therefore, the general relation between ionization and metallicity seems to describe star formation in galaxies at both low and high redshifts.

Zoom In Zoom Out Reset image size

Inspecting Figure 14 more closely, there is also evidence that the distribution of $12+\mathrm{log}({\rm{O}}/{\rm{H}})$ for the z > 1 galaxies skews to higher ionization parameters at fixed metallicity (at $12+\mathrm{log}({\rm{O}}/{\rm{H}})\gtrsim 8.4$). This is more apparent in the KDE distributions: compared to the nearby H ii regions (compared to both Pérez-Montero 2014 and Espinosa-Ponce et al. 2022). The implication is that galaxies at higher redshift favor higher ionization parameters (see, also, Kewley et al. 2015; Strom et al. 2017; Kaasinen et al. 2018; Topping et al. 2020; Runco et al. 2021; see, also, Sanders et al. 2020, who discuss how a decrease in the effective temperature of stars can lead to an increase in ionization parameter). Here we see this trend exists for high-redshift galaxies even at fixed metallicity and stellar mass.

What causes the increase in ionization parameter? To understand the answer, we divided a sample of galaxies into bins of $\mathrm{log}q$. We selected galaxies at 1.1 < z < 2.3 from our CLEAR sample in a narrower range of stellar mass, $9.3\leqslant \mathrm{log}{M}_{* }/{M}_{\odot }\leqslant 9.7$, and we required that the galaxy gas-phase metallicities be $12+\mathrm{log}({\rm{O}}/{\rm{H}})\gt 8.3$ (to remove low-metallicity galaxies from consideration). We then divided this sample into a high-ionization subsample of galaxies with $\mathrm{log}q\gt 7.8$ (31 galaxies) and a low-ionization subsample of galaxies with $\mathrm{log}q\lt 7.8$ (32 galaxies). The cuts in stellar mass and ionization have the effect of making both the stellar-mass distribution and metallicity distribution approximately the same for the high- and low-ionization subsamples (see below, and Figure 15). So, we are able to correlate trends between ionization and other galaxy properties.

Zoom In Zoom Out Reset image size

We then stacked the WFC3 G102 and G141 dust-corrected spectra for these subsamples (following the methods in Section 3.1). Figure 15 shows these stacked spectra for the "low-ionization" ($\mathrm{log}q\lt 7.8$) and "high-ionization" subsamples ($\mathrm{log}q\gt 7.8$). The differences in the spectra are immediately clear (pun intended). The two samples have very different emission-line intensities, while the stellar continua of the two stacks are nearly identical. Figure 15 also shows the relative differences between the two spectra. The strongest difference is in the [O III] emission, which is nearly twice as strong at the peak of the line for the high-ionization galaxies. The Balmer lines are also stronger by ≈30% stronger at the peak of the lines for the high-ionization galaxies. The [O II] and [Ne III] lines both show an increase in the higher-ionization subsample, while the [S II] lines do not.

The bottom row of panels in Figure 15 shows the distributions of stellar mass, gas-phase metallicity, SED-derived sSFR, and the Hα/Hβ ratios for the high- and low-ionization galaxy subsamples. Both subsamples were selected to have roughly the same stellar-mass and gas-phase metallicities, and the panels show there are no substantive differences. Formally, a Kolmogorov–Smirnov (K-S) test and a Mann–Whitney-U (MWU) test applied to the distributions return p-values of 0.60 and 0.34, respectively, for stellar mass, and 0.45 and 0.71, respectively, for metallicity. However, the specific SFR distributions show evidence they are different. The K-S and MWU tests return p-values 0.07 and 0.08, respectively. This is apparent in Figure 15 as a shift in the sSFR distribution of the high-ionization subsample toward higher sSFR. The final (right) panel of Figure 15 shows the distribution of the Balmer decrement, but this includes only those galaxies with z < 1.63 for which Hα is detected (14 and 23 galaxies in the high- and low-ionization subsamples, respectively). The mode of the high-ionization subsample is consistent with E(B − V) = 0, and the low-ionization subsample shows slightly higher attenuation with a mode of E(B − V) = 0.2, but both have low overall dust attenuation. The K-S and MWU tests yield p-values of 0.45 and 0.69, respectively, providing no evidence they are drawn from different parent distributions (but this is in part because of the smaller sample sizes).

Why then do the galaxies in the high-ionization subsample have such stronger emission lines when all other galaxy properties appear to be similar? The key may be in the fact that we see a correlation between these (broadband-derived) specific SFR and ionization parameters (see Figure 13). Figure 16 shows the region around Hβ and Hγ in the stacked spectra of the high- and low-ionization subsamples. We focus on these two lines as they are the strongest Balmer emission lines that are unblended at the WFC3 G102 and G141 spectral resolution (e.g., Hα is blended with [N II]; H is blended with [Ne iii] 3968). Both of these Balmer lines (Hβ and Hγ) are stronger by ≈20%–30% at the peak of the lines in the high-ionization galaxies.

Zoom In Zoom Out Reset image size

We measured the strength of these emission lines using pPXF applied to the stacked spectra (following Section 3.2). Table 1 reports the emission-line EWs. For all of the Balmer lines, the high-ionization subsample has substantially stronger lines, with the EW of Hβ and Hγ higher by a factor of ≃1.7–1.9 than those of the low-ionization subsample. Similarly, Figure 17 here shows that correlation between both the Balmer-line EW against the ionization parameter ($\mathrm{log}q$) and against the [O III]/[O II] ratio. In Figure 17, we show the Hβ EW as it is measured, and we have increased the Hγ EW according to the theoretical Hβ/Hγ Case-B ratio.

Zoom In Zoom Out Reset image size

The stronger Balmer emission (as evidenced by higher EW) in the high-ionization galaxies indicates that they have a higher production rate of H-ionizing photons (at fixed galaxy stellar mass). This translates to higher specific SFRs. For example, Reddy et al. (2018) showed that log EW(Hβ) = $0.32\times \mathrm{log}$ sSFR for galaxies at z ∼ 2 in MOSDEF. This translates to a factor of two higher sSFR for the high-ionization galaxies in our sample compared to the low-ionization galaxies. We illustrate this in Figure 17, which shows the Hβ EW compared to both the O₃₂ ratio and the ionization parameter (${\rm{\Delta }}\mathrm{log}q$) in the high-ionization galaxies relative to the low-ionization galaxies. The lines in Figure 17 show fits,

$$\begin{eqnarray}&&{\mathrm{logO}}_{32}=1.4\mathrm{logEW}({\rm{H}}\beta )-2.3,\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}&&\mathrm{log}q=1.2\mathrm{logEW}({\rm{H}}\beta )+5.8,\end{eqnarray} \tag{ 9 }$$

for the stacked spectra from the CLEAR subsamples). Using the relation from Reddy et al. (2018), we find that q ∼ sSFR^0.4 for our CLEAR galaxy samples. Figure 17 also shows how these relations compare to the $\mathrm{log}q$–sSFR (derived from broadband fitting, see Figure 13, and converted to EW(Hβ) using this relation), which is similar. The conclusion is that the ionization parameter and specific SFR are correlated for galaxies at fixed stellar-mass and gas-phase metallicity.

Our result for the CLEAR galaxies is consistent with findings from some previous studies. Kewley et al. (2015) showed strong evidence that the [O III]/[O II] ratio is correlated with the Hβ EW in their sample of 0.2 < z < 0.6 galaxies (but they did not differentiate by stellar mass). This is also seen by Kaasinen et al. (2018), who argue that the high-ionization parameters for galaxies in their z < 0.3 sample are driven by the ratio of the number density of hydrogen-ionizing photons to the gas density, and not by (lower) metallicities. In a study of SDSS galaxies, Kashino & Inoue (2019) similarly found that the ionization parameter is predominantly controlled by the specific SFR, with q ∼ sSFR^0.43 (converting to our units), approximately equal to the relation for the CLEAR galaxies.

We are now ready to address the question, what is the origin of the correlation between gas ionization parameter (q) and specific SFR? First, we consider several explanations for the correlation, including selection and physical effects.

6.4.1. Is a Selection Effect in Stellar Mass Responsible?

6.4.2. Is a Selection Effect in Gas-phase Metallicity Responsible?

Figure 14 shows there is a correlation between gas-phase metallicity ($12+\mathrm{log}({\rm{O}}/{\rm{H}})$) and ionization parameter ($\mathrm{log}q$). The distributions of metallicities of the high- and low-ionization are again highly similar (Figure 15, and discussion in Section 6.3). A potential selection effect arises because the $12+\mathrm{log}({\rm{O}}/{\rm{H}})$ values are derived from what surmounts to the O₃₂ and R₂₃ line ratios. R₂₃ is famously double-valued, and the fitting method we employed could skew the metallicities to higher values (as the line ratios are unable to get "over the R₂₃ hump" to lower-metallicity values). Indeed, this accounts for the "ridgeline" in the modes of the metallicity–ionization-parameter distributions seen in Figure 14. (This is in essence asking if our modeling prefers the "high-Z" solution).

To test if this potential bias affects the correlation, we used the [Ne III] λ3868/[O II] ratio an as alternative tracer of gas-phase metallicity. Multiple studies have demonstrated the correlation between [Ne III]/[O II] and $12+\mathrm{log}({\rm{O}}/{\rm{H}})$ (e.g., Maiolino et al. 2008; Bian et al. 2018; Pérez-Montero et al. 2021; Sanders et al. 2021). [Ne III]/[O II] is particularly useful as the lines are close in wavelength (mitigating extinction effects) and are accessible at rest-frame optical wavelengths. While the [Ne III] ratio is often too weak to detect in individual galaxies ($\mathrm{log}[\mathrm{Ne}\,{\rm\small{III}}]/[{\rm{O}}\,{\rm\small{II}}]\sim -1$; see Table 1, although see Backhaus et al. 2022), it is detectable in the stacked spectra.

Figure 15 shows that the [Ne III]/[O II] ratios are similar in the stacked spectra of the low- and high-ionization samples. This is supported by the measured [Ne III]/[O II] ratios measured from the stacked spectra in Table 1. Figure 18 shows the [Ne III]/[O II] against $12+\mathrm{log}({\rm{O}}/{\rm{H}})$ derived from the stacked spectra (compared to the CLEAR stacks for the full sample in bins of stellar mass). Based on the [Ne III]/[O II] ratios, there is only a small difference between the values for the high and low subsamples ($\mathrm{log}[\mathrm{Ne}\,{\rm\small{III}}]/[{\rm{O}}\,{\rm\small{II}}]\,\simeq \,$−1.0 for both), and this is likely attributed to the difference in the ionization (see Kewley et al. 2019b; Maiolino & Mannucci 2019). This is important as the [Ne III] fluxes are not used in the measurements of $12+\mathrm{log}({\rm{O}}/{\rm{H}})$. Indeed, we measured [Ne III]/[O II] in a stacked spectrum of CLEAR galaxies selected to have low metallicity, $12+\mathrm{log}({\rm{O}}/{\rm{H}})$ < 8.2, which shows they have much higher ratios, $\mathrm{log}[\mathrm{Ne}\,{\rm\small{III}}]/[{\rm{O}}\,{\rm\small{II}}]\,\simeq$−0.6 ± 0.1. This is consistent with other studies of low-Z galaxies at these redshifts (e.g., Sanders et al. 2021), and this is substantially higher than that observed in the high- and low-ionization subsamples. Therefore, we disfavor the explanation that metallicity is responsible for the correlation between $\mathrm{log}q$ and sSFR.

Zoom In Zoom Out Reset image size

6.4.3. Physical Connections between the Specific SFR and Ionization Parameter

There are physical reasons to expect a correlation between the specific SFR and the ionization parameter. Such mechanisms need to relate the sSFR (or specifically, the production of ionizing photons per unit mass) to the number density of gas particles. Brinchmann et al. (2008) showed this can account for the emission-line ratios of low-redshift galaxies (from SDSS, at z ∼ 0.1). They further argue the elevated ionization parameters result from higher gas densities possibly combined with higher escape fractions of hydrogen-ionizing photons (f_gas). Similar correlations between specific SFR and ionization parameter (and/or line ratios such as [O III]/[O II]) have been observed previously at low and high redshifts (e.g., Nakajima & Ouchi 2014; Kewley et al. 2015; Sanders et al. 2016; Bian et al. 2016; Kaasinen et al. 2018; Kashino & Inoue 2019). In our CLEAR sample, we see that this correlation persists even when accounting for stellar mass and metallicity (see Figures 13 and 15). Therefore, there is strong evidence for a physical connection between specific SFR and ionization.

A physical connection exists in the case where the H ii regions are radiation bounded, where the ionization parameter can be written as (Charlot & Longhetti 2001; Nakajima & Ouchi 2014; Stasińska et al. 2015),

$$\begin{eqnarray}&&q\approx {\left(\ \displaystyle \frac{3Q\ {n}_{{\rm{H}}}\ {\epsilon }^{2}\ {\alpha }_{B}^{2}}{4\pi }\ \right)}^{1/3},\end{eqnarray} \tag{ 10 }$$

where Q is the rate of hydrogen-ionizing photons, α_B is the case-B hydrogen recombination coefficient, n_H is the number density of hydrogen, and is the volume filling factor of the gas. ²⁰ At present there is little evidence for changes in in any environment (see, Brinchmann et al. 2008), and we do not consider this further here. The ionization parameter therefore depends on Q and n_H, which observations show are higher at higher redshift (for the case of n_H; see, e.g., Sanders et al. 2016 and Strom et al. 2017) or are expected to be higher at higher redshift (for the case of Q, given the observed bluer colors and lower metallicities of galaxies; see, e.g., Shivaei et al. 2018). This is a major driver of the evolution in the MQR (Figure 12). The geometry of H ii regions and star-forming nebula in distant galaxies may be more complicated. This could include nebulae that are either "matter bounded" (also called "density bounded"), where the region of ionization (the Strömgren sphere) extends beyond the nebula (see, e.g., Nakajima & Ouchi 2014) and/or there are denser clouds with a low covering fraction (e.g., Naidu et al. 2022). We consider these effects in the discussion that follows.

Production rate of ionizing photons. Q depends on the age of the stellar population, metallicity, and relative number of massive stars (typically those with ≳ 10 M_⊙). Therefore, to increase Q requires either an increase in the relative number of massive stars and/or evolution in the properties of the stellar populations. One obvious possibility is that there is a change in the IMF: either an increase in the upper-mass cutoff, or high-mass slope. Currently evidence for such a changes are inconclusive (e.g., Finkelstein et al. 2011; Narayanan & Davé 2013), so we do not consider it further. However, it will remain important to consider for future studies.

Several studies modeling the UV stellar continua and optical emission-line strengths have advocated for higher Q values in some star-forming galaxies at both high and low redshifts, particularly when the metallicity of the stellar continuum is low (Z ∼ 0.1 Z_⊙; e.g., Topping et al. 2020; Berg et al. 2021; Olivier et al. 2021). This may be exacerbated by super-solar α/Fe ratios, which lead to harder ionizing spectra (at fixed O/H). There is growing evidence for increased α/Fe in galaxies at z ∼ 2 (e.g., Steidel et al. 2016; Strom et al. 2018; Shapley et al. 2019; Sanders et al. 2020; Topping et al. 2020; Runco et al. 2021). One limitation of our work is that the MAPPINGS models do not currently provide for variations in α/Fe, and it will important to test how this impacts the trends in our data set. Nevertheless, for this to drive our observations, α/Fe would need to vary with sSFR at fixed [O/H], which would itself be an important discovery. Although our current data set is insufficient at present, this may be testable in the future by simultaneously modeling the rest-UV continuum spectra and measurements of the optical emission lines of galaxies (e.g., Topping et al. 2020; Olivier et al. 2021).

Gas Density. There is considerable evidence that the gas densities of galaxies at z ∼ 2 are considerably higher than at low redshift (e.g., Shirazi et al. 2014; Sanders et al. 2016; Acharyya et al. 2019; Runco et al. 2021). This been argued to contribute to the elevated emission-line ratios in high-redshift galaxy samples (Brinchmann et al. 2008; Nakajima & Ouchi 2014; Kewley et al. 2019a). For changes in the gas density to account for our observations, these also must correlate with the specific SFR. One physical explanation for such a correlation is an extension of the Kennicutt–Schmidt relation, where the SFR density scales with gas density, ${\rho }_{\mathrm{SFR}}\sim {\rho }_{\mathrm{gas}}^{N}$ (see, e.g., Bacchini et al. 2019). The exponent scales from N = 1.4 (for the classical Kennicutt–Schmidt relation) to 2 depending on timescales over which star formation occurs (see Madore 1977; Larson 1981; Kennicutt et al. 2007; Krumholz 2014; Elmegreen 2015; Bolatto et al. 2017; Bacchini et al. 2019). There are observations of galaxies at low redshift that indicate higher O₃₂ ratios for galaxies with smaller sizes (implying higher densities; e.g., Z. Ji & M. Giavalisco 2022, in preparation). Regardless, this provides a physical connection between the strength of H-recombination lines (tracing the SFR density) and the ionization parameter (tracing gas density in H ii nebula).

A connection between specific SFR and ionization parameter is therefore expected, and borne out in our CLEAR data set. Testing if the gas density drives this correlation will be feasible using observations of emission lines whose ratios are dependent on gas density. For example, the ratio of the S⁺ emission lines, [S ii] λ6716/[S ii] λ6731, varies by a factor of ∼2 as the gas density changes from n_e ≃ 10 to 1000 cm⁻³ (Ryden & Pogge 2021). Currently, HST/WFC3 grism data have insufficient resolution to deblend [S II] for galaxies at z ≲ 1.5. However, future spectroscopy at higher spectral resolution would be able to test for density variations directly.

Ionizing Radiation Escape Fractions. Several studies of low-redshift galaxies (z ≲ 0.1) show evidence that the escape fraction H-ionizing photons, f_esc, is correlated with the ionization parameter, for example where ${f}_{\mathrm{esc}}\propto {({{\rm{O}}}_{32})}^{2}$ (Chisholm et al. 2018; Izotov et al. 2018; Z. Ji & M. Giavalisco 2022, in preparation). This means that the geometry of the H ii regions is likely an important factor, where nonuniform covering fractions and/or "matter-bounded" nebulae can lead to higher escape fractions, f_esc > 0 (e.g., Nakajima & Ouchi 2014; Naidu et al. 2022). At a fixed ionization parameter, a nonzero f_esc requires a higher number density of ionizing photons relative to the gas density. Using a suite of simulations, Giammanco et al. (2005) showed that increasing the escape fraction of hydrogen-ionizing photons from f_esc = 0 to 0.5 increases $\mathrm{log}q$ by as much as 1 dex. This can be boosted in the case of a "density-" (or "matter-") bounded nebula, where the O⁺ region in the H ii region is truncated, while the O⁺⁺ region, located closer to the ionizing source, is not (e.g., Nakajima & Ouchi 2014). This would lead to a potential correlation between f_esc and the [O III]/[O II] ratio (see also Brinchmann et al. 2008 and Kashino & Inoue 2019). There is recent evidence for this in the findings of Naidu et al. (2022), who argue that high [O III]/[O II] line emission of z ∼ 2 galaxies corresponds to a short-lived phase where massive stars have cleared sightlines in the nebular gas. This produces a geometry with a smaller covering fraction of dense gas, permitting higher f_esc (and higher [O III]/[O II]). As pointed out by Brinchmann et al. (2008, and others above), such an increase in f_esc at higher redshift is an important factor in interpreting the emission-line ratios of high-redshift galaxies. There is an increasing number of observations (either direct measurements or inferences) that the escape fraction of H-ionizing photons is higher at high redshift (Vanzella et al. 2016, 2018; de Barros et al. 2016; Shapley et al. 2016; Bian et al. 2017; Ji et al. 2020; Begley et al. 2022), where inferences for reionization require f_esc ∼ 0.1−0.2 (e.g., Ouchi et al. 2009; Robertson et al. 2013; Bouwens et al. 2015; Ishigaki et al. 2018; Finkelstein et al. 2019).

The fact that we observe a correlation between ionization parameter and specific SFR in our CLEAR galaxies then implies there may be a correlation between specific SFR and f_esc. This prediction stems from the fact that the ionization parameter and escape fraction are expected to correlate (Nakajima & Ouchi 2014), and we observe a correlation between specific SFR and $\mathrm{log}q$ (Figure 13) and between Hβ EW and $\mathrm{log}q$ (Figure 17). Brinchmann et al. (2008) speculated that such a correlation could account for offsets in the emission-line ratios of SDSS galaxies (see, also, Nakajima & Ouchi 2014 and Kashino & Inoue 2019). Currently there are only weak constraints on a correlation between specific SFR or ionization parameter and f_esc at high redshifts, though some recent observations are suggestive (but sample sizes remain small; e.g., Bassett et al. 2019 and Naidu et al. 2022). Future observations of the escape fraction in high-z galaxies like those in our CLEAR sample (either with UV spectroscopy or imaging; see Siana et al. 2010), probing a large range in specific SFR, [O III]/[O II] ratio, would provide the data to test this directly.

We summarize this section by postulating that multiple effects likely contribute to the correlation between specific SFR and ionization parameter. While the current data set is insufficient to differentiate these effects, we favor the explanation that both an increase in the gas density (n_H) and the escape fraction of H-ionizing fractions (f_esc) drive the trend in specific SFR and ionization parameter, as these have the strongest physical bases. This will be testable directly with future data.