phase abundance patterns, and enhanced contributions from AGNs and shocks at high redshift (see e.g. Kewley et al. 2013, for a review). Initial results from the MOSFIRE Deep Evolution Field (MOSDEF; Kriek et al. 2015) survey, suggested that the observed offset is primarily caused by an enhanced N/O ratio at fixed oxygen abundance, in addition to higher physical densities in high-redshift galaxies compared to local systems (Masters et al. 2014;Shapley et al. 2015;Sanders et al. 2016b). Results from the Keck Baryonic Structure Survey (KBSS; Steidel et al. 2014) instead suggested that the offset is driven primarily by a harder intrinsic ionizing spectrum at fixed oxygen abundance (Steidel et al. 2016;Strom et al. 2017). Updated results from the MOSDEF survey support the explanation of a harder ionizing spectrum at fixed oxygen abundance (Sanders et al. 2019(Sanders et al. , 2020;;Shapley et al. 2019;Runco et al. 2020).

Furthermore, a harder stellar ionizing spectrum at fixed oxygen abundance (i.e. nebular metallicity, which traces α/H) arises naturally in the case of α enhancement of the most massive stars, which produce the bulk of the ionizing radiation in star-forming galaxies. For a given oxygen abundance, such α-enhanced stars will have lower Fe/H than their local counterparts with a solar abundance pattern. A lower Fe/H corresponds to a harder ionizing spectrum. The αenhancement described above is expected in high-redshift galaxies due to their young ages and rising star-formation histories (e.g. based on SED-fit ages in previous papers, as well as preference for rising SFHs in high-redshift galaxies; Papovich et al. 2011;Reddy et al. 2012) resulting in enrichment primarily from Type II supernovae (Steidel et al. 2016;Matthee & Schaye 2018).

The rest-optical emission lines observed in high-redshift starforming galaxies are strongly affected by the intrinsic ionizing spectrum primarily produced by the most massive stars. Several properties of the massive stars modulate the production of ionizing photons, including stellar metallicity, IMF, and stellar binarity (Topping & Shull 2015;Stanway, Eldridge & Becker 2016;Steidel et al. 2016;Stanway & Eldridge 2018). In addition to controlling the ionizing spectrum, the formation of massive stars is regulated by gas accretion on to galaxies, and in turn regulates the resulting chemical enrichment of the ISM through the deposition of metals by supernova explosions and stellar winds, and the ejection of metals through star-formation-driven galaxy-scale outflows. Strom et al. (2018) investigated the relationship between properties of the ionized ISM and factors contributing to the excitation state (including stellar metallicity) for a sample of z ∼ 2 star-forming galaxies, finding that typically, z ∼ 2 galaxies have lower stellar metallicity compared to their nebular metallicity. This investigation relied purely on strong rest-optical emission lines. However, the emission-line spectrum depends on a host of physical properties in addition to the ionizing spectral shape and is thus a highly indirect tracer of the intrinsic ionizing spectrum. In particular, breaking the degeneracy between the ionization parameter and the ionizing spectral shape is challenging when only the strongest rest-optical emission-line ratios are available. Sanders et al. (2020) used a small sample of star-forming galaxies at z ∼ 1.5-3.5 with direct-method nebular metallicities to connect factors affecting the ionizing spectrum with properties of the ISM. The tight constraints on stellar metallicity facilitate breaking the degeneracy between ionization parameter and hardness of the ionizing spectrum. Another avenue for breaking this degeneracy is directly constraining the ionizing spectrum from massive stars using rest-UV spectra of these high-redshift galaxies.

While directly observing the ionizing spectrum within highredshift galaxies is extremely challenging, information about the massive star population can be determined based on the non-ionizing rest-UV stellar spectrum that is far easier to observe. Specifically, features including the C IV λλ1548, 1550 and He II λ1640 stellar wind lines, and a multitude of stellar photospheric features are sensitive to the properties of massive stars (Leitherer et al. 2001;Rix et al. 2004;Crowther et al. 2006). For example, Halliday et al. (2008) utilized Fe III absorption lines to measure a sub-solar stellar metallicity for a composite spectrum of z ∼ 2 galaxies. Sommariva et al. (2012) tested the ability of additional photospheric absorption line indicators to accurately determine the stellar metallicity of the massive star population of high-redshift galaxies. Recently, instead of using integrated regions of the rest-UV spectrum, Cullen et al. (2019) instead fit stellar population models to the full rest-UV spectrum of multiple composite spectra to investigate galaxy properties across 2.5 < z < 5.0 and a stellar mass range of 8.5 < log(M * /M ) < 10.2.

Crucially, recent studies have concurrently studied the production of the ionizing spectrum with the rest-optical properties of highredshift galaxies by utilizing combined rest-UV and rest-optical spectroscopy (Steidel et al. 2016;Chisholm et al. 2019;Topping et al. 2019). Steidel et al. (2016) constructed a composite spectrum of 30 z ∼ 2.4 star-forming galaxies from KBSS, and found that their rest-UV properties could only be reproduced by stellar population models that include binary stars, have a low stellar metallicity (Z * /Z ∼ 0.1), and a different, higher, nebular metallicity (Z neb /Z ∼ 0.5). By analysing a single composite rest-UV spectrum, Steidel et al. (2016) only probed average properties of their high-redshift galaxy sample. With single rest-UV and rest-optical composite spectra it is not possible to probe the rest-UV spectral properties as a function of the location in the BPT diagram. In contrast, Chisholm et al. (2019) fit linear combinations of stellar population models to 19 individual galaxy rest-UV spectra at z ∼ 2, and determined light-weighted properties. Using these individual spectra, Chisholm et al. (2019) found using mixed age populations that galaxies in their sample had comparable stellar and nebular metallicities. In Topping et al. (2019) we compared the properties of two composite spectra one of which included galaxies lying along the local sequence and the other including galaxies offset towards high [N II] λ6584/H α and [O III] λ5007/H β. This analysis indicated that galaxies offset from the local sequence had younger ages and lower stellar metallicities on average, resulting in a harder ionizing spectrum in addition to a higher ionization parameter and a higher α-enhancement, all of which contributed to the difference in BPT diagram location. Intriguingly, we found that even high-redshift galaxies coincident with local star-forming sequence on the BPT diagram were αenhanced (O/Fe ∼ 3 O/Fe ) relative to their local counterparts, suggesting that coincidence in the BPT diagram does not necessarily imply similar ISM conditions.

In this paper, we improve the methods presented in Topping et al. ( 2019) by expanding the stellar population synthesis models to consider more complex star formation histories (SFHs), and including additional rest-optical emission lines to the photoionization modelling. Furthermore, we test the capability of the models to be fit to individual galaxy spectra that have lower SNR than composite spectra, and analyse a sample of ∼30 individual galaxies with combined rest-UV and rest-optical spectra.

The organization of this paper is as follows: Section 2 describes our observations, data reduction, and methods. Section 3 presents the results of our analysis, and Section 4 provides a summary and discussion of our key results. Throughout this paper, we assume a cosmology with m = 0.3, = 0.7, H 0 = 70 km s -1 Mpc -1 , and adopt solar abundances from Asplund et al. (2009, i.e. Z = 0.014, 12 + log(O/H) = 8.69).

Our analysis utilizes rest-optical spectroscopy of galaxies from the MOSDEF survey (Kriek et al. 2015) at z ∼ 2.3, observed using the Multi-Object Spectrometer for Infra-Red Exploration (MOSFIRE; McLean et al. 2012McLean et al. ) over 48.5 nights during 2012McLean et al. -2016. This restoptical spectroscopic sample is composed of ∼1500 near-infrared spectra at moderate resolution (R ∼ 3500) of H-band selected galaxies targeted to lie within three distinct redshift intervals (1.37 ≤ z ≤ 1.70, 2.09 ≤ z ≤ 2.61, and 2.95 ≤ z ≤ 3.80). Based on the scatter between photometric and spectroscopic redshifts of the MOSDEF targets, the actual redshift ranges slightly differ from the initial target ranges. Therefore, we define the true redshift ranges as 1.40 ≤ z ≤ 1.90, 1.90 ≤ z ≤ 2.65, and 2.95 ≤ z ≤ 3.80. In addition to rest-optical spectra from the MOSDEF survey, MOSDEF targets have extensive ancillary data sets from the CANDELS (Grogin et al. 2011) and 3D-HST (Momcheva et al. 2016)

A full description of the rest-UV data collection and reduction procedures will be described in a future work, but we provide a brief description here. We selected a subset of MOSDEF galaxies for rest-UV spectroscopic follow-up using the Low-Resolution Imaging Spectrograph (LRIS; Oke et al. 1995). Target priorities were determined using the following prescription. Highest priority was given to objects from the MOSDEF survey that had detections in all four BPT emission lines (H β, [O III], H α, [N II] λ6584) with ≥3σ . Then, objects were added to the sample with detections in H β, [O III] λ5007, and H α with ≥3σ , and an upper limit in [N II] λ6584. With decreasing priority, the remaining targets were selected by having a spectroscopic redshift measurement from the MOSDEF survey, objects from the MOSDEF survey without a successfully measured redshift, and objects not targeted as part of the MOSDEF survey, but that were part of the 3D-HST survey catalogue (Momcheva et al. 2016) with properties within the MOSDEF survey photometric redshift and apparent magnitude ranges. In total, these targets comprise a sample of 260 galaxies. Of those targets with spectroscopic redshifts from the MOSDEF survey, 32, 162, and 20 were in the redshift ranges 1.40 ≤ z ≤ 1.90, 1.90 ≤ z ≤ 2.65, and 2.95 ≤ z ≤ 3.80, respectively. The remaining galaxies, with either a spectroscopic redshift not from MOSDEF, or a photometric redshift, made up 9, 31, and 6 galaxies in the redshift intervals 1.40 ≤ z ≤ 1.90, 1.90 ≤ z ≤ 2.65, and 2.95 ≤ z ≤ 3.80, respectively. For this analysis, we excluded the few objects identified as AGN based on their X-ray and rest-frame near-IR properties.

Rest-UV spectra were obtained using Keck/LRIS during ten nights across five observing runs between 2017 January and 2018 June. Our target sample totals 260 distinct galaxies on nine multiobject slit masks with 1. 2 slits in the COSMOS, AEGIS, GOODS-S, and GOODS-N fields. In order to obtain continuous wavelength coverage from the atmospheric cut-off at 3100 Å up to a median red wavelength limit of ∼ 7650 Å, we observed all slit masks using the d500 dichroic, the 400 lines mm -1 grism blazed at 3400 Å on the blue side, and the 600 lines mm -1 grating blazed at 5000 Å on the red side. This setup yielded a resolution of R ∼ 800 on the blue side, and a resolution of R ∼ 1300 on the red side. The exposure times ranged from 6 to 11 h on different masks, with a median exposure time of ∼7.5 h for the full sample. The data were collected with seeing ranging from 0. 6 to 1. 2 with typical values of 0. 8.

We reduced the red-and blue-side data from LRIS using custom IRAF, IDL, and PYTHON scripts. First, we fit polynomials to the edges of each 2D slitlet, and transformed it to be rectilinear. The subsequent steps required slightly different treatment for the red and blue spectral images. We flat fielded each image using twilight sky flats for the blue side, and dome flats for the red side images. Then we cut out each slitlet for each object in every flat-fielded exposure. Following this step, the blue-side slitlets were cleaned of cosmic rays, and background subtracted. These images were registered and median combined to create a stacked 2D spectrum. In order to prevent overestimation of the background due to the target, we measured the trace of each object in the stacked 2D spectrum, and masked it out for a second-pass background subtraction (Shapley et al. 2006). As the red-side images were more significantly affected by cosmic rays, we reduced the red side using slightly different methods. For the red-side slitlets, we constructed a stacked 2D spectrum by first registering and median combining the images using minmax rejection to remove cosmic rays. We used this stacked image to measure the object traces in each slitlet. We then recomputed the background subtraction in the individual images with the object traces masked out, as the stacked image is too contaminated by sky lines to achieve a good background subtraction. After the second pass background subtraction, the individual red-side slitlets were combined to create the final stacked image.

We extracted the 1D spectrum of each object from the red and blue side stacked slitlets. We calculated the wavelength solution by fitting a fourth-order polynomial to the red and blue arc lamp spectra, which resulted in residuals of ∼ 0.035 Å and ∼ 0.3 Å, respectively. We repeated this step on a set of frames that had not had sky lines removed. Using the resulting sky spectra, we measured the centroid of several sky lines and shifted the wavelength solution zero-point until the sky lines aligned with their known wavelengths. This shift typically had a magnitude of ∼ 4 Å throughout the sample. We applied an initial flux calibration based on spectrophotometric standard star observations obtained during each observing run. We checked the flux calibration by comparing spectrophotometric measurements calculated from our objects to measurements in the 3D-HST photometric catalogue, and applied a multiplicative factor to correct our calibration. Throughout our sample, this factor had median magnitude of ∼ 30 per cent. Following the final flux calibration, we ensured that continuum levels on either side of the dichroic at 5000 Å were consistent. A small number of galaxies in our sample required additional wavelength-dependent continuum correction near ( 50 Å) the dichroic cutoff. To apply this correction we defined a target continuum, F targ cont using regions of the spectrum just outside of the area that required correction. We then fit a linear function to the part of the observed spectrum to be corrected, F obs cont . Finally, we multiplied this piece of the spectrum by F targ cont /F obs cont to recover the correct shape of the spectrum near the dichroic cutoff.

While the full sample described above consists of 260 galaxies across three distinct redshift intervals (1.40 ≤ z ≤ 1.90, 1.90 ≤ z ≤ 2.65, and 2.95 ≤ z ≤ 3.80), we focus on a subset of this sample composed of galaxies in the central redshift window that have detections in four primary BPT lines (H β, [O III] λ5007, H α, [N II] λ6584) at ≥3σ from the MOSDEF survey. This 'LRIS-BPT' sample comprises 62 galaxies, each of which has a systemic redshift measured from rest-optical emission lines.

Fig. 1 shows the galaxy properties of objects in the LRIS-BPT sample. Topping et al. (2019) constructed composite spectra from two subsets of the LRIS-BPT sample, each subset being defined based on their location on the BPT diagram (Fig. 1, left-hand panel). These two subsets, referred to as the high and low samples, comprise galaxies that are offset from the local BPT sequence, and those that lie along the sequence, respectively. Fig. 1 (righthand panel) shows the SFR measured from dust-corrected Balmer lines and stellar mass for the galaxies in our sample. The SFR and stellar mass of this sample are characterized by a median and inner 68th percentile log(SFR/(M / yr)) = 1.53 ± 0.44 and log (M * /M ) = 10.02 ± 0.52, respectively. These median values are consistent with those of the sample of galaxies comprising the central redshift range (1.90 ≤ z ≤ 2.65) of the full MOSDEF survey, which are characterized by median and inner 68th percentile of log( SFR/(M /yr)) = 1.36 ± 0.50 and log (M * /M ) = 9.93 ± 0.60. These comparisons of SFR and stellar mass suggest that our LRIS-BPT sample is an unbiased subset of the full z ∼ 2 MOSDEF sample.

For this analysis, we used the version 2. Eldridge 2018). Notably, these stellar population models incorporate the effects of stellar rotation, quasi-homogeneous evolution, stellar winds, and binary stars. These effects can have a substantial effect on the spectrum of a model stellar population, and in particular, the EUV spectrum produced by massive, short-lived stars. The BPASS models are evaluated with multiple initial mass functions (IMFs), including the Chabrier (2003) IMF, and IMFs with high-mass (M ≥ 0.5M ) slopes of α = -2.00, -2.35, and -2.70. In addition, the models using each IMF were formulated with a high-mass cutoff of 100 M and 300 M . For this analysis, we only considered models that assume a Chabrier ( 2003) IMF, and a high-mass cutoff of 100 M . Finally, the BPASS models have been computed with ages between log (Age/yr) = 6.0-11.0 in increments of 0.1 dex, and stellar metallicities of Z * = 10 -5 -0.04. While we considered all available stellar metallicities in our analysis, we restricted the ages to log (Age/yr) = 7.0-9.6. At ages younger than log (Age/yr) = 7.0 we would be probing time-scales shorter than the dynamical time-scale of the galaxies, and therefore could not accurately attribute physical properties to the entire galaxy simultaneously (Reddy et al. 2012;Price et al. 2020). This minimum age restriction is further motivated by the H α EW distribution of our sample, which is characterized by a median value of 150 Å, consistent with galaxies that have ages of 10Myr (Erb et al. 2006). Even the most extreme galaxies in our sample have H α EWs of ∼ 800 Å, which is consistent with our range of stellar population model ages (see Reddy et al. (2018) for the H α EW distribution of the full MOSDEF sample). Additionally, the ionizing spectrum of stellar populations with constant star formation reaches an equilibrium at an age ∼10 7 yr. Finally, at the lowest redshift galaxy in our sample, the age of the Universe was ∼4 Gyr, so including older stellar populations is not necessary.

We constructed stellar population models that assume different star-formation histories (SFH) by combining the BPASS models, which describe a coeval stellar population, using

where t max is the age of the population, t 0 is the time between the age of the first model (i.e. 10 6 yr), and the onset of star formation, t 0 is the SFR at the time of the first model, f (λ) t 0 is the luminosity density of the model spectrum per unit stellar mass at the time of the first model, t i is the star-formation rate of the population at time t i , f (λ) tmax-t i is the luminosity density of the model spectrum per unit stellar mass with age t maxt i (i.e. the model that began t i yr prior to the final age, t max ), and (t it i -1 ) is the time between subsequent model spectra. For the case of a constant SFH, all of the SFR weightings, t i , are set to unity. Each spectrum was constructed using models with the same stellar metallicity (i.e. none of our model spectra represent a combined population of multiple stellar metallicities). In addition to a constant SFH, we considered several models with a 'delayed-τ ' SFH of the form SFR ∝ t × e -t/τ , with log (τ /yr) = 7, 8, 9, 10. With this set of models we covered three schematically different SFHs. These SFHs allowed us to explore models when star formation is rising (t < τ), falling(t > τ), and peaked (t ∼ τ ).

We processed the stellar population model spectra using the photoionization code CLOUDY (Ferland et al. 2017). Using this code, we input an ionizing spectrum from BPASS and, given a set of ISM properties, calculated expected emission line fluxes. We compared the simulated line fluxes to the observed rest-optical emission lines of galaxies in our sample to infer properties of the ISM. We assumed a fixed electron density of n e = 250 cm -3 , which is representative of the galaxies in our sample (Sanders et al. 2016a). In addition, while we vary the nebular oxygen abundance, we assume solar abundance ratios for most elements. However, we adopt the log (N/O) versus Table 1. Summary of model grid parameters. The age and stellar metallicity values correspond to BPASS models we fit to our observed spectra. For each combination of age and stellar metallicity, we computed a set of photoionization models with the listed nebular metallicity and ionization parameter values. On this scale log (Z neb /Z ) = 0 corresponds to 12 + log(O/H) = 8.69. log (Age/yr) 7.0, 7.3, 7.4, 7.5, 7.7, 8.0, 8.5, 9.0, 9.6 Z * 10 -5 , 10 -4 , 0.0001, 0.001, 0.002, 0.003, 0.004, 0.006, 0.008, 0.01, 0.014, 0.02, 0.03, 0.04 log (Z neb /Z ) -1.3, -1.0, -0.8, -0.6, -0.5, -0.4, -0.3, -0.2, -0.1, 0.0, 0.1, 0.2 log (U) -3.6, -3.4, -3.2, -3.0, -2.8, -2.6, -2.4, -2.2, -2.0, -1. For each BPASS model, we ran a grid of CLOUDY models for a range of nebular metallicity (-1.6 ≤ log (Z neb /Z ) ≤ 0.3) and ionization parameter (-3.6 ≤ log (U) ≤ -1.4). We have made several updates to the parameters of the model grids described in Topping et al. (2019), in order to more finely sample the parameter space in regions of interest, including amending the abundances of all elements to the values reported by Asplund et al. (2009). Table 1 summarizes the parameters, and lists each value for which we compute a model. An additional component of the photoionization models is the nebular continuum. The nebular continuum contributes a relatively small amount of flux to the UV spectrum, compared to the stellar component (see Steidel et al. 2016; Fig. 3). We explicitly compute the nebular continuum for BPASS models listed in Table 1, however, it changes smoothly with age, so we interpolate the nebular contribution for the remaining BPASS models.

To compute a composite spectrum, we first interpolated each of the individual galaxy spectra on to a common wavelength grid. We chose the sampling of this common wavelength grid to be 0.8 Å, which corresponds to the rest-frame sampling of our spectra at the median redshift of our sample. Then, we continuum normalized each of the individual galaxy spectra. In the process of continuum fitting, we first extracted regions of the spectra that are not contaminated by absorption lines, in the windows defined by Rix et al. (2004). We then fit a cubic spline to the median flux values within each window to define the continuum level. Because of this approach, we did not need to consider effects that smoothly affect the continuum (e.g. reddening by dust). Following this step, at each wavelength, we median combined all spectra that had coverage at that wavelength. Only 3/62 of the galaxies in our sample do not cover the full wavelength range that we are fitting, therefore a composite spectrum at any given wavelength is well sampled by the individual galaxy spectra it comprises. We defined the error spectrum as the standard deviation of all contributing spectra at each wavelength.

To fit the BPASS stellar population synthesis models to our individual galaxy and composite spectra, we masked out regions of the observed spectra that include components not present in the models (e.g. strong interstellar absorption lines). Then, we continuum normalized both the observed and BPASS model spectra. We then interpolated the BPASS models on to the wavelength grid of the galaxy spectra. We did not smooth either the models or the observed spectra as their resolutions were comparable with values of ∼ 1 Å in the rest-frame. Following this step, we calculated the χ 2 statistic for each BPASS model in the grid, and determined which age and stellar metallicity produced the minimum χ 2 value. We determined the uncertainties in these parameters by perturbing the observed spectrum and calculating which age and stellar metallicity best fit the observed spectrum. In the case of an individual galaxy spectrum, this perturbation simply constitutes adding in noise to each wavelength element pulled from a normal distribution with a standard deviation defined by the magnitude of the error spectrum at that wavelength element. For a simulated composite spectrum, we selected a new sample of galaxies from the initial composite spectrum sample with replacement. Then, each galaxy spectrum was perturbed using the method described above before being combined. After repeating this process 1000 times, we defined the best-fitting value and upper and lower 1σ uncertainties as the median, 16th and 84th percentile of the distribution.

We expanded the model grid used in Topping et al. (2019), which included only stellar population models assuming a constant SFH. We repeated fitting the model grids to the high and low stacks using our updated models that assume different SFHs. For each SFH, we found that the results are consistent with those of Topping et al. (2019). In particular, we fit models that assume a 'delayed-τ ' SFH, with log (τ /year) = 7, 8, 9, and 10. Fig. 2 shows the best-fitting age and stellar metallicity of the high and low stacks for each model Fractional difference between the best-fitting model spectra from models using a delayed-τ SFH compared to the model spectrum using a constant SFH. Best-fitting models fit to the high and low stacks are displayed in the top and bottom panels, respectively. The regions masked out using 'mask1' from Steidel et al. (2016) defined to include contamination from non-stellar sources is shown in grey. On average, the models assuming a delayed-τ SFH differ from those with a constant SFH at the few per cent level.

grid. For each SFH, we find that the high stack has lower stellar metallicity, and a younger stellar age compared to the low stack. While this trend between the properties of the high and low stacks persists for each SFH we considered, the exact values of the stellar metallicity and age differ between the assumed models. In particular, the delayed-τ model with log (τ /yr) = 7 has the youngest best-fitting age and stellar metallicity, and both the age and stellar metallicity increase when assuming an increasing τ .

For each of the assumed SFHs, we recovered the same qualitative trend found in Topping et al. (2019), according to which the high stack had a younger age and lower stellar metallicity relative to the low stack. This result assumes that both the high and low stacks are both described by a similar SFH, which is a reasonable assumption based on the SFHs determined for galaxies in these two stacks estimated from SED fitting. We also investigated if the stellar population models yielded any constraint on the form of the SFH for each stack. We tested this question by measuring the minimum χ 2 value for the best-fitting model of each SFH. For the high and low stacks, none of the SFHs were preferred over the others, suggesting that a given UV stellar spectrum is not unique to a particular SFH. A KS test revealed that the high and low stacks have τ distributions with a 45 per cent probability of being drawn from the same parent distribution. Fig. 3 compares the best-fitting spectrum of a constant SFH model to models with a delayed-τ SFH. The best-fitting models for the high stack are nearly identical, and at some wavelengths are different at the few per cent level. The low stack models vary more, in particular for the log (τ /yr) = 7 model, which has some signatures of a young population not seen in the other best-fitting models. In particular, this model has slightly enhanced C IV λλ1548, 1550 and He II λ1640 emission compared to the other models. However, the majority of the models are in agreement, with differences of only up to ∼ 10 per cent in a few wavelength elements. These tests have shown that we cannot determine which SFH best characterizes the rest-UV spectra. However, using models that assume different SFHs does not qualitatively affect our results.

While the composite spectra provide useful constraints due to their high SNR, the SNR of individual galaxy spectra can be much lower. For example, the high and low stacks have SNR/pixel ∼ 18, while the individual spectra have a median SNR/pixel ∼ 4.5, where the sampling of our spectra are ∼ 0.8 Å pixel -1 in the rest frame. These values correspond to SNR per spectral resolution element ∼ 25.5 and SNR per spectral resolution element ∼ 6.3 for the composite and median individual spectra, respectively. We tested how the SNR of a spectrum affects the best-fitting stellar properties by manually introducing noise to one of our composite spectra, refitting the models, and checking if biases arise as the spectrum drops in quality. Fig. 4 displays how the best-fitting stellar metallicity and age change as a function of the amount of added noise. For this composite, the best-fitting stellar metallicity retains an unbiased estimate of the value obtained in the high-SNR limit down to an SNR per spectral resolution element ∼ 5.6 (SNR/pixel∼4). Below this value, the measurement of stellar metallicity is biased high. However, even above this limit, the stellar metallicity uncertainty increases with decreasing SNR. The best-fitting age remains consistent throughout the range of SNR per spectral resolution element, yet as the SNR decreases below SNR per spectral resolution element ∼ 5.6, the uncertainty grows to ≥2 dex, leaving the age unconstrained.

This test showcases how biases in the best-fitting stellar parameters may occur in lower SNR spectra. In order to further quantify this effect, we repeated the process used on the composite spectrum, except on BPASS models for which the 'true' parameters are known. We added noise selected from a normal distribution to the BPASS models at each wavelength element. We repeated this for a combination of ages and stellar metallicities to determine if these biases exist throughout the range of parameter space our SFHs span. Fig. 5 shows how the best-fitting age and stellar metallicity changes as noise is introduced into the model spectra. At all stellar metallicities, a low SNR/resolution element introduces a positive bias in stellar Best-fitting age and stellar metallicity for all galaxies in the LRIS-BPT sample. For completeness, galaxies with SNR per spectral resolution element ≤ 5.6 are displayed as faint grey symbols. The sample comprises galaxies with ages in the range 7.0 ≤ log (Age/yr) ≤ 9.6, with the majority of galaxies having stellar metallicities of 0.0005 ≤ Z * ≤ 0.004. metallicity leading to an overestimate of Z * relative to the true value. At the lowest SNR, this bias can be up to ∼ 150 per cent. The best-fitting stellar age also changes at low SNR/resolution element, however in contrast to the bias of the stellar metallicity, the trend of the bias depends on the 'true' age. At low SNR/resolution element, models with an old input age are biased younger, and models with a young age are biased towards older values. While these biases exist at low SNR/resolution element for both age and stellar metallicity, the input values of Z * and age can be accurately determined for spectra with SNR/resolution element ≥∼5.6. Attempting to model the spectra with lower SNR leads to high uncertainties and systematic biases.

Based on our tests on the BPASS models, we can achieve accurate age (i.e. time since the onset of star formation) and stellar metallicity measurements for individual spectra with SNR per spectral resolution element ≥ 5.6. Fig. 6 displays these values for all galaxies in the LRIS-BPT sample, highlighting those with sufficiently high SNR to yield unbiased Z * and age. The stellar metallicity for the individual galaxies ranges between 0.001 ≤ Z * ≤ 0.006, consistent with the best-fitting metallicities found for our composite spectra. Fig. 7 shows some of our individual galaxy spectra and the best-fitting BPASS model to each spectrum. This figure shows the variety of stellar parameters recovered by our fitting by depicting a low, medium, and high stellar metallicity. Importantly, based on the stellar metallicity and age we are able to constrain the shape of the ionizing spectrum for each of these individual galaxies. However, for some galaxy spectra, the stellar metallicity probability distribution uniformly spanned the parameter space, and we could not assign a stellar metallicity in such cases.

The ionizing spectrum emitted from the stellar population drives the production of the emergent rest-optical emission line ratios. Therefore, constraining the ionizing spectrum within galaxies is crucial in order to use photoionization models to extract physical properties from the observed nebular emission lines. In particular, using CLOUDY, we set the ionizing spectrum based on the appropriate best-fitting BPASS stellar population model, and vary the nebular metallicity and ionization parameter, recording the emergent line fluxes. We then compare the resulting catalogue of nebular emission line flux ratios to those observed from an individual or composite spectrum. The inferred nebular metallicity and ionization parameter is set by which model best reproduces the observed emission lines defined by the minimum χ 2 calculated for all observed line ratios simultaneously. To understand the uncertainty in these quantities, we perturb the observed emission line flux ratios by their corresponding uncertainties, and recompute the best-fitting nebular parameters. Topping et al. (2019) used an approach that compared the locations of the models and observed galaxies on the [N II] λ6584/H α versus [O III] λ5007/H β BPT diagram. In this analysis, we use a slightly different approach that simultaneously fits the [N II] λ6584/H α, and [O III] λ5007/H β emission line ratios. Furthermore, we include additional strong line ratios, [S II] λλ6717, 6731/H α, and [O II] λ3727/H β, in order to better constrain the nebular parameters. Because [O II] λ3727 and H β are farther apart in wavelength compared to the other line ratios, we apply a dust correction to both lines. We investigate the effect that including these additional observables has on the inferred nebular parameters. Fig. 8 shows an example of how the constraint on nebular metallicity and ionization parameter changes for different sets of nebular emission lines when our method is applied to the low stack composite spectrum. In this example, the inferred nebular properties are consistent when considering different sets of lines, however the ionization parameter and nebular metallicity are better constrained when additional emission lines are included.

One assumption made in the photoionization modelling is the form of the N/O versus O/H relation. The median nitrogen abundance of H II regions in the local Universe has been measured to vary by ∼0.5 dex for 8.2 ≤ 12 + log (O/H) ≤ 8.6, with scatter of ∼0.2 dex at fixed O/H (Pilyugin et al. 2012). This assumption strongly affects the output [N II] λ6584/H α flux in our photoionization models.

These nitrogen abundance variations can result in a disparity of the [N II] λ6584/H α ratio ∼0.5 dex, resulting in a biased inference of Z neb and log (U). Fig. 9 shows the effect of removing [N II] λ6584/H α from our fitting procedure, eliminating the uncertainty surrounding the N/O relation. Without [N II] λ6584/H α, the inferred ionization parameters are well matched to those inferred when using [N II] λ6584/H α, with nearly all galaxies falling on the one-to-one relation. In addition, there is a slight difference at low O/H when [N II] λ6584/H α is not considered, however this difference is small compared to the uncertainties. Based on this result, we conclude that the inclusion of [N II] λ6584/H α in our analysis does not significantly bias our inferred nebular parameters. 2020) found a median log (U) = -2.30 ± 0.06 for z ∼ 2 galaxies at 12 + log(O/H) ∼ 8.0, also consistent with the local U-O/H relation. The measurements of log (U) and 12 + log (O/H) for the galaxies in this high-SNR sample in combination with measurements of lower metallicity galaxies from Sanders et al. (2020) suggest that there may be a trend of decreasing ionization parameter with increasing nebular metallicity consistent with the trend seen in Pérez-Montero (2014). The offset of this sample from the z ∼ 0 excitation sequence in the BPT diagram is thus not driven by higher ionization parameters at fixed O/H (Cullen et al. 2016;Kashino et al. 2017;Bian et al. 2020). Furthermore, this result remains largely the same when considering log (U) and 12 + log (O/H) inferred without [N II] λ6584/Hα. However, the different methods used to measure the oxygen abundances between our high-SNR sample, and those from Pérez-Montero (2014), which used the 'direct' method may introduce systematics. In particular, Esteban et al. (2014) demonstrated that metallicities measured using the direct method are ∼0.24 lower than those which use nebular recombination lines on average. However, metallicities measured from recombination lines are in better agreement with those inferred from photoionization modelling. In addition the ionized gas properties presented by Pérez-Montero (2014) were inferred using softer ionizing spectra compared to those used in this paper. The result is a potential overestimate of the ionization parameter measured for the Lines of constant α-enhancement (i.e. O/Fe) are displayed as solid, dashed, and dotted lines, respectively, for O/Fe , 2×O/Fe , and 5×O/Fe . All galaxies in our sample show evidence for α-enhancement. Additionally, some galaxies are in the regime above 5×O/Fe , which has been suggested as the theoretical limit based on supernova yield models. The galaxies for which the stellar metallicity could not be determined are displayed as lower limits.

local systems compared to galaxies in our high-SNR sample. While addressing this potential systematic could move our measurements away from the Pérez-Montero (2014) relation, the difference in methods used to infer oxygen abundances would act in such a way to remedy the disparity. Therefore, despite these systematics, we conclude that there is no strong offset in the U versus O/H relation of galaxies in our sample compared to the local relation

To connect the stellar and nebular properties of the individual galaxies in our high-SNR sample, we combine the nebular O/H abundance inferred from photoionization modelling with the Fe/H measured from the BPASS model fitting to look at the α-enhancement of individual galaxies. Fig. 11 compares the nebular metallicity and the stellar metallicity for each galaxy in our sample. Noticeably, all of the galaxies in our sample show evidence for α-enhancement, falling significantly above the O/Fe line (solid black line). These values range from ∼1.75O/Fe to ≥5O/Fe . Investigating a small sample of z ∼ 2-3.5 galaxies with direct oxygen abundance measurements, Sanders et al. (2020) found similar results, where most of the galaxies analysed show evidence for α-enhancement, some having values > 5 × O/Fe . A number of objects in our sample fall above the expected theoretical limit from pure Type II SNe of ∼ 5 × O/Fe , assuming a typical IMF (Kobayashi et al. 2006;Nomoto et al. 2006). However, this limit depends on the details of the stellar population and expected Type II SNe yields. For example, the theoretical O/Fe limit increases when calculated assuming a top-heavy IMF. Therefore, different assumptions of the IMF or supernova yields could remove the tension between the theoretical limit and some of our observed galaxies. Finally, the tension between our observed O/Fe and the theoretical limit may be relieved due to uncertainties in the stellar modelling of low-metallicity massive stars. If, in fact, the low-metallicity stars produce harder ionizing radiation than what is included in the BPASS models, our galaxies would be best described by stellar metallicities that are higher than what we have found. A higher best-fitting stellar metallicity could bring the α-enhancement into better agreement with the theoretical limits.

We check how the α-enhancement changes as a function of distance away from the local BPT sequence for our sample of highredshift galaxies. We calculate the distance away from the local sequence by first finding the line that is perpendicular to the local sequence from Kewley et al. (2013)

where N2 = log([N II]λ6584/Hα) and O3 = log([O III]λ5007/Hβ).

We calculated uncertainties by perturbing each rest-UV spectrum and fitting the grid of BPASS stellar population models to measure the stellar metallicity. The best-fitting stellar population parameters define an ionizing spectrum which we then use as an input to CLOUDY, from which, in turn, we infer the nebular oxygen abundance. We define the 1σ uncertainty as the 16th and 84th percentile of the distribution of α-enhancements resulting from repeating this process 1000 times. Fig. 12 displays how the α-enhancement depends on our measured distance from the local BPT sequence. There appears to be a weak positive correlation between the distance away from the local BPT sequence and α-enhancement. While there is a large amount of scatter in this relation, we do see that the galaxies with the highest O/Fe typically lie farther from the local sequence. We also compare the measurements of individual galaxies in our high-SNR sample with the high and low composite spectra from Topping et al. (2019). These two composite spectra follow the trend of higher O/Fe at a greater BPT distance, and both composite spectra are consistent with the results of our high-SNR sample. While the data are not constraining enough to measure a functional form of a relation between α-enhancement and BPT distance, they are consistent with the result that the offset of high-redshift galaxies relative to local galaxies on the BPT sequence is significantly driven by a harder ionizing spectrum due to α-enhancement.

Since the first evidence that suggested that high-redshift galaxies are offset on the BPT diagram, several hypotheses have been proposed to explain the underlying cause. Among the proposed sources for this offset between local and high-redshift galaxies are harder ionizing spectra at fixed nebular metallicity, higher electron densities, contributions from AGNs and shocks at high redshift, and variations in gasphase abundance patterns. Recently, two prevailing theories suggest that the offset is primarily driven by higher ionization parameters at fixed gas-phase metallicity (Kewley et al. 2015;Kashino et al. 2017;Bian, Kewley & Dopita 2018;Cullen et al. 2018), or that highredshift galaxies exhibit a harder intrinsic ionizing spectrum at fixed nebular metallicity driven by α-enhancement at high redshift (Steidel et al. 2016;Sanders et al. 2020).

To answer this question of the origin of the BPT offset, Sanders et al. (2020) used the 'direct' method to estimate oxygen abundances for a sample of 18 high-redshift galaxies at low nebular metallicities. These authors then used photoionization models very similar to our own to fit for log (U) and Z * after fixing Z neb to match the direct metallicity, finding that the average of the sample lies along log U versus 12 + log(O/H) of local H II regions (Pérez-Montero 2014). This result suggests that the high ionization parameter measured in their sample is due to their low nebular metallicity, and that their sample has consistent ionization parameter with local H II regions that share the same O/H. Furthermore, Shapley et al. (2019) demonstrated that high-redshift galaxies are also offset toward higher [S II] λλ6717, 6731/H α and [O III] λ5007/H β when using the appropriate comparison to local galaxies with low contribution from diffuse ionized gas (DIG) within the ISM. Using photoionization models from Sanders et al. (2016a), Shapley et al. (2019) concluded that both the offset on the [N II] and [S II] BPT diagrams is best explained by a harder ionizing spectrum at fixed nebular metallicity.

The results described in this paper suggest that z ∼ 2 galaxies do not have an elevated ionization parameter compared to local H II regions that share the same 12 + log (O/H). Our analysis illustrates the importance of an independently constrained ionizing spectrum when trying to extract ionization parameter information from strong rest-optical line ratios. Without such a constraint, the degeneracy between ionization parameter and the intrinsic ionizing spectrum can bias inferences of the ionization parameter. It is important to note that the method used to infer oxygen abundances of our sample is different from the method of our local H II region comparison sample (Pérez-Montero 2014), which could introduce systematics in the comparison. The offset between these two methods results in ∼0.24 dex lower oxygen abundance when using the direct method relative to the nebular recombination lines (see e.g. Sanders et al. 2020 for a more detailed discussion). However, the modelling from Pérez-Montero (2014) assumed a softer ionizing spectrum at fixed oxygen abundance compared relative to the BPASS models used for our analysis. Using a softer ionizing spectrum results in an ionization parameter that is systematically higher for a fixed set of nebular emission line ratios (Sanders et al. 2016a). Therefore, correcting for these biases affecting the inferred ionization parameter and oxygen abundance would shift our inferred values along the U versus O/H relation described by Pérez-Montero (2014), and would not lead to a significant difference between the nebular parameters of local H II regions and the z ∼ 2 galaxies in our sample.

Finally, we address the possibility of differing abundance scales between our stellar and nebular fitting procedures using a few different methods. We remeasured the best-fitting stellar metallicities using Starburst99 stellar population synthesis models (Leitherer et al. 2014) and compared the results to those found using BPASS models. In the range of stellar metallicities spanned by the majority of galaxies in our sample, (0.001 Z * 0.004) the Starburst99-and BPASSbased stellar metallicities agree within the uncertainties. However, for some galaxies the stellar metallicity found by BPASS is <0.001, the minimum metallicity in the Starburst99 grid. In addition, both BPASS and Starburst99 models have only been computed using solar abundance patterns. For this analysis, altering the abundance patterns within the models is unlikely to significantly affect our results, as the differences between models of varied stellar metallicities is primarily driven by the Fe abundance (Steidel et al. 2016), while the other elements have subdominant effects. To address the gasphase abundance scale, the ideal case would be to utilize directmethod metallicities. However, making the necessary measurements to calculate direct metallicities is challenging (Sanders et al. 2020), and requires a data set that does not exist for a large sample and in the metallicity range spanned by our galaxies. Instead, we compared the inferred nebular properties presented above using photoionization modelling to properties inferred using an empirical method described by Sanders et al. (2020) using the calibrations from Bian et al. (2018). We find that the empirical measurements of 12 + log(O/H) are on average 0.12 dex lower than those inferred from photoionization modelling. Despite this difference, when using oxygen abundances inferred from empirical calibrations we find that all of our galaxies show evidence for α-enhancement. One other consideration is the possibility of rapid oxygen enrichment of the H II regions resulting in a different nebular abundance relative to that of the massive stars. However, analysis of the Orion star-forming region finds the gas-phase oxygen abundance is in reasonable agreement with the oxygen abundance massive stars (Simón-Díaz 2010; Simón-Díaz & Stasińska 2011). Based on our stellar and nebular results, we find that the offset on the BPT diagram is primarily due to a harder ionizing spectrum resulting from subsolar stellar Fe abundances relative to local galaxies.

We used the combination of rest-UV and rest-optical spectra for a sample of 62 galaxies to investigate the physical conditions within galaxies at z ∼ 2.3. We expanded upon the results of Topping et al. (2019) in which we constructed composite spectra based on location in the [N II] λ6584 BPT diagram, and found that galaxies offset from the local sequence typically had younger ages, lower stellar metallicities, higher ionization parameters, and were more α-enhanced. We expanded the fitting analysis to include additional SFHs and rest-optical emission line fluxes. In addition, we quantitatively determined the rest-UV continuum SNR limit above which we can constrain the stellar metallicity and age of individual galaxies in an unbiased manner. We summarize our main results and conclusions below.

(i) We constructed additional BPASS stellar population models for a variety of SFHs. We repeated the fitting analysis of the two stacked spectra defined by Topping et al. (2019) and found that for each SFH, the stack composed of galaxies offset from the local sequence on the BPT diagram had a younger age and lower stellar metallicity. Additionally, when fitting across all SFHs for a single stack, we do not find any preference for one SFH over another. Therefore, we cannot determine which SFH best characterizes the rest-UV spectra.

(ii) We tested the SNR down to which individual galaxy spectra are suitable to be fit using this type of analysis. Based on the test of perturbing model spectra with known stellar parameters with increasing amounts of noise, we find that the best-fitting stellar metallicity is biased high when the spectrum reaches a SNR/resolution element < 5.6, and the best-fitting age is biased towards the middle of the grid (log(Age/yr) ∼ 8.5), with an uncertainty that fills the parameter space. The best-fitting age and stellar metallicity remain consistent with the true value for spectra with SNR per spectral resolution element > 5.6. Therefore, rest-UV spectral fitting should not be attempted on spectra for which the SNR/resolution element is less than 5.6 in order to avoid biased results.

(iii) Based on the SNR requirements described above, we found that 30 galaxies in our sample satisfied the criteria to be fit on an individual basis. We find that galaxies in this high-SNR sample have ages that have a wide ranges of ages spanning log(Age/yr) ∼ 7.0 -9.6, and that most galaxies have stellar metallicities in the range ∼0.001 < Z * < 0.004.

(iv) We examined how different rest-optical emission lines affect the inferred ionization parameter and nebular metallicity. Previously, we inferred nebular parameters by comparing observed [N II] λ6584/H α and [O III] λ5007/H β with a suite of photoionization models. In this analysis, we tested how adding [S II] λλ6717, 6731/H α and [O II] λ3727/H β to the fitting procedure affects the resultant parameters. In general, adding the additional lines yields results with smaller uncertainties. Additionally, because one assumption we made is the form of the N/O versus O/H relation as an input to our models, we tested fitting the rest-optical lines that are not affected by this assumption, namely [O III] λ5007/H β, [O II] λ3727/H β, and [S II] λλ6717, 6731/H α. We find that when [N II]λ6584/Hα is excluded from fitting, the nebular metallicity and ionization parameter remain consistent with the values inferred when [N II]λ6584/Hα is included.

(v) With the constrained ionizing spectrum for each individual galaxy, we used photoionization models to infer ionization parameters and nebular metallicities for each galaxy in our sample, and find that the inferred ionization parameters (log (U) med = -2.98 ± 0.25) are consistent with those measured in local H II regions that share the same oxygen abundance (12 + log (O/H) med = 8.48 ± 0.11). This result suggests that the offset of high-redshift galaxies on the BPT diagram relative to local galaxies is not due to elevated ionization parameters at fixed O/H.

(vi) Combining the best-fitting stellar metallicities from fitting BPASS model spectra to the nebular metallicities inferred from photoionization modelling, we find that all of our individual galaxies are α-enhanced compared to local galaxies. Furthermore, we find that high-redshift galaxies that show the highest levels of α-enhancement are typically more offset from the local BPT star-forming sequence. The O/Fe values span the range from ∼1.75 ×O/Fe to above the theoretical limit (∼5 ×O/Fe ; Nomoto et al. 2006). This limit could be affected by details of the IMF, including the high-mass slope and upper mass cutoff. Furthermore, the stellar metallicities may change as stellar modelling of the most metal-poor massive stars are better understood. In particular, if the lowest Z * stars actually produce harder ionizing spectra compared to current models, we would infer a higher stellar metallicity for our rest-UV spectra. In combination with (v), this result suggests that the offset of high-redshift galaxies on the BPT diagram relative to the local sequence is likely due to a harder ionizing spectrum at fixed nebular oxygen abundance resulting from elevated O/Fe in high-redshift galaxies.

A combined understanding of the intrinsic ionizing spectrum that excites H II regions and the physical properties of the ISM is a crucial step towards a complete model of high-redshift galaxy evolution. In order to fully understand high-redshift galaxies we must explore how their properties differ from local galaxies, but also how the population of high-redshift galaxies varies within itself. Ultimately, detailed modelling of large numbers of individual galaxies will be required to expand our understanding of galaxies beyond the level of a sample average, but our work here begins to illustrate the range and interplay of nebular and stellar properties observed during the epoch of peak star formation in the Universe. Constraining the ionizing spectrum is crucial in order to use photoionization models to infer nebular properties, as without this constraint solutions are highly degenerate. Rest-UV spectroscopy is the ideal tool to gain insight into the massive star population, and therefore the ionizing spectrum, for individual galaxies at high redshift. This type of analysis is key in order to compare the internal properties of high-redshift galaxies to those of local H II regions and galaxies.