also considered the [O III]λ5008/Hβ versus [S II]λλ6718,6733/Hα diagram (originally introduced in Veilleux & Osterbrock 1987 and referred to hereafter as the '[S II] BPT diagram') when investigating the ionization mechanism of a galaxy. Indeed, the shape of the ionizing spectrum as well as the typical ionization parameter are different in gas excited by an AGN as opposed to by hot stars. Accordingly, AGNs and star-forming galaxies occupy distinct regions within rest-optical emission-line diagrams. These diagnostic diagrams also reveal information about the physical properties of the galaxies. For example, stellar mass (M * ) and metallicity have been found to increase with decreasing [O II]λ5008/Hβ and increasing [N II]λ6585/Hα along the local star-forming sequence in the [N II] BPT diagram (e.g. Masters, Faisst & Capak 2016).

Another diagram commonly used to describe star-forming galaxies is the [O III]λλ4960,5008/[O II]λλ3727,3730 (O 32 ) versus ([O II]λλ4960,5008+[O II]λλ3727,3730)/Hβ (R 23 ) diagram. O 32 and R 23 are rough tracers for ionization parameter and metallicity, respectively. This diagram allows us to probe such physical quantities in star-forming galaxies (e.g. Lilly, Carollo & Stockton 2003;Nakajima et al. 2013). For local galaxies, there is an increase in metallicity from the high-excitation end (high O 32 & R 23 ) to the low-excitation tail (low O 32 and R 23 ) on this diagram (Andrews & Martini 2013;Shapley et al. 2015).

Early studies with Keck/NIRSPEC found that the location of galaxies on the [N II] BPT diagram is redshift-dependent, as galaxies with z > 1 are found to be offset from local Sloan Digital Sky Survey (SDSS; York et al. 2000) galaxies, showing elevated [O III]λ5008/Hβ at fixed [N II]λ6585/Hα (or vice versa; e.g. Shapley et al. 2005;Erb et al. 2006;Liu et al. 2008). Based on new observations with multiobject near-infrared spectrographs on 8-10-m class telescopes, the sample of high-redshift galaxies with measurements of the BPT diagram emission lines now numbers in the hundreds. Of note, two large surveys of the high-redshift BPT diagram include the MOSFIRE Deep Evolution Field (MOSDEF; Kriek et al. 2015) survey and the Keck Baryonic Structure Survey (KBSS: Steidel et al. 2014), which find, based on much more robust statistical evidence, that high-redshift galaxies tend to have elevated [O III]λ5008/Hβ and/or [N II]λ6585/Hα compared to local galaxies (Steidel et al. 2014;Shapley et al. 2015Shapley et al. , 2019;;Sanders et al. 2016;Strom et al. 2017).

Understanding this offset in emission-line ratios is vital because we use strong rest-optical emission lines as empirical tracers for many physical properties (e.g. gas-phase oxygen abundance Pettini & Pagel 2004) of the ISM. Due to the observed offset for z > 1 galaxies, it is unclear if local metallicity calibrations, e.g. Pettini & Pagel (2004), yield accurate metallicities when applied in the high-redshift universe. Therefore, it is essential to gain a complete understanding of why high-redshift galaxies have elevated [O III]λ5008/Hβ and/or [N II]λ6585/Hα compared to the local SDSS sample.

There have been many proposed explanations of this systematic offset. Possible explanations include variations in physical properties of galaxies such as H II region electron densities (or proportionally pressures), density structure, H II region ionization parameter, H II region ionizing spectra at fixed metallicities, gasphase N/O abundance ratio differences, unresolved AGN activity, and shocks (e.g. Brinchmann, Pettini & Charlot 2008;Liu et al. 2008;Wright et al. 2010;Kewley et al. 2013;Yeh et al. 2013;Masters et al. 2014;Coil et al. 2015;Sanders et al. 2016;Steidel et al. 2016;Strom et al. 2017;Freeman et al. 2019;Kashino et al. 2019;Shapley et al. 2019;Topping et al. 2020b). Galaxy selection effects could also be potential factors in this offset (Juneau et al. 2014).

Preliminary results from the MOSDEF survey (Shapley et al. 2015;Sanders et al. 2016) suggested that the offset in emissionline ratios is primarily due to elevated N/O at fixed O/H abundance patterns in z ∼ 2.3 galaxies relative to local ones. Studies using KBSS data (Steidel et al. 2016;Strom et al. 2017Strom et al. , 2018) ) argued for a harder ionizing spectrum at fixed nebular metallicity as the main cause of the offset while other works using the Fiber Multi-Object Spectrograph (FMOS)-COSMOS survey (Kashino et al. 2017(Kashino et al. , 2019) ) or local analogues of high-redshift galaxies (Bian et al. 2020) attribute the observed [N II] BPT offset to a higher ionization parameter. More recent results from the MOSDEF survey using the complete MOSDEF data set now suggest a harder ionizing spectrum drives the BPT offset (Shapley et al. 2019;Sanders et al. 2020b, a;Topping et al. 2020b).

The completed MOSDEF survey provides an ideal data set of highredshift galaxies to help explore the observed offset on the [N II] BPT diagram. This survey provides access to spectra of ∼1500 galaxies in the z ∼ 1.4-3.8 redshift range, with full spectroscopic coverage of all the emission lines needed to complete the [N II] BPT, [S II] BPT, and O 32 versus R 23 diagrams in the z ∼ 1.4-2.6 redshift range.

In this work, we improve upon the previous z > 1 BPT offset studies based on early MOSDEF data by Shapley et al. (2015) and Sanders et al. (2016). We now have the full MOSDEF sample in hand and apply a more careful spectral energy distribution (SED) fitting method that incorporates corrections to broad-band photometric measurements for rest-optical emission-line fluxes. Accordingly, we derive unbiased age and stellar mass measurements, and can conduct a more thorough and complete investigation of the emission-line ratio properties of the MOSDEF sample. As an example of the level of improvement enabled by the full MOSDEF sample, the initial work by Shapley et al. (2015) explored the location of 53 star-forming galaxies at 1.9 ≤ z ≤ 2.7 on the [N II] BPT diagram, and investigated the M * , specific SFR (sSFR), and SFR surface density ( SFR ) of the sample. In this study, the sample increases to 180 star-forming galaxies on the [N II] BPT diagram, and we also investigate additional galaxy parameters such as the galaxy effective radius (R e ), SFR, n e , and stellar population age. Sanders et al. (2016) investigated how a sample of 53 star-forming galaxies at 1.9 ≤ z ≤ 2.7 translated from the [N II] BPT diagram to the [S II] BPT diagram and O 32 versus R 23 diagram. In the current study, that number is more than doubled to 122. Our more comprehensive analysis enables a better understanding of the observed offset between the MOSDEF sample and local galaxies on the [N II] BPT diagram. Through the analysis outlined above, we aim to determine the underlying physical causes of this offset.

The paper is organized as follows. In Section 2, we briefly describe the MOSDEF survey, and review our survey sample selection and data reduction. Section 3 presents the results of this study. Section 4 discusses the results and how they relate to photoionization models. Section 5 presents a summary of key results and looks ahead to future analyses. We adopt the following abbreviations for emissionline ratios used frequently throughout the paper.

O 32 and R 23 have already been defined, but we include them here for completeness. All emission-line wavelengths are vacuum wavelengths. Throughout this paper, we adopt a Lambda cold dark matter cosmology with H 0 = 70 km s -1 Mpc -1 , m = 0.3, and = 0.7. Also, we assume the solar abundance pattern from Asplund et al. (2009).

The MOSDEF survey was a 48. The sample is H-band magnitude limited, and is located in the wellstudied CANDELS and 3D-HST legacy fields (Grogin et al. 2011;Koekemoer et al. 2011;Momcheva et al. 2016): AEGIS, COSMOS, GOODS-N, GOODS-S, and UDS. The ancillary data from these fields enables us to probe other properties of the sample (e.g. galaxy size and stellar mass).

In this study, we restrict the sample to the central redshift bin, spanning the redshift range to 1.9 ≤ z ≤ 2.7 to reflect the scatter between target photometric and measured spectroscopic redshifts (i.e. the fact that not all spectroscopic redshifts, when actually measured, fell precisely within the target 2.09 ≤ z ≤ 2.61 window). For galaxies in this redshift range, [O II]λλ3727,3730, Hβ, [O III]λλ4960,5008, Hα, [N II]λ6585, and [S II]λλ67186733 features are captured within the J, H, and K near-IR filter wavelength ranges. For MOSFIRE, the J, H, and K filters have typical respective wavelength coverages of 1. 142-1.365, 1.450-1.826, and 1.897-2.427 μm and spectral resolutions of R = 3000, 3650, and 3600.

We imposed further restrictions to reach our final sample. Galaxies containing an AGN were identified and removed from the sample based on their IR colours, X-ray luminosity, or if N2 > 0.5 (Coil et al. 2015;Azadi et al. 2017;Leung et al. 2017). We also require an S/N ≥ 3 for each emission-line in our analysis. These criteria yielded a sample of 180 galaxies for initial classification in the [N II] BPT diagram, of which 122 galaxies are additionally detected in [S II]λλ6718,6733 and [O II]λλ3727,3730 and have size measurements from the Hubble Space Telescope (HST). The redshift distribution of the 180 MOSDEF star-forming galaxies at z ∼ 2.3 with [N II] BPT diagram classifications is shown in Fig. 1.

We correct Hα and Hβ line fluxes for stellar Balmer absorption as described in Kriek et al. (2015) and Reddy et al. (2015). We also dust correct emission-line ratios for which the member features differ significantly in wavelength, including O 32 and R 23 . For such corrections, we assumed the Cardelli, Clayton & Mathis (1989) dust attenuation curve and an unreddened Hα/Hβ ratio of 2.86. In the case of emission-line ratios for which the lines are close in wavelength (O3, N2, and S2) dust correction was not applied.

In this study, we investigate additional galaxy photometric, spectroscopic, and structural properties for the MOSDEF sample in order to try to better understand the observed rest-optical emission-line properties. The Hao et al. (2011) calibration for a Chabrier ( 2003) initial mass function (IMF) and solar metallicity is used to estimate SFR(Hα) from stellar-Balmer-absorption-corrected, dust-corrected, and slit-loss-corrected Hα luminosities (Kriek et al. 2015;Reddy et al. 2015;Shivaei et al. 2015). We used the SED-fitting code, FAST (Kriek et al. 2009) to obtain key stellar population parameters including stellar mass and age. For this modelling, as a default, we assumed star formation histories of the delayed-τ form, where SFR ∝t × e -t/τ . Here, t is the time since the onset of star formation (i.e. age), and τ is the characteristic star formation time-scale. Given the range of best-fitting τ values, the meaning of absolute ages is not necessarily clear. In order to obtain a better gauge of the relative maturities of the galaxy stellar populations in our sample, we used both normalized ages (i.e. t/τ ), and also the age (t) obtained from constant star formation (CSF) models. Galaxy sizes, R e , are taken as the F160W galaxy half-light radii from the van der Wel et al. (2014) catalogue, which were estimated using single-component Sérsic profile fits to the two-dimensional light distribution of galaxies in the CANDELS and 3D-HST fields. We combine SFR(Hα) and R e to estimate the SFR surface density, SFR as

We estimate the sSFR using SFR(Hα) and M * , both of which are described above. The electron density, n e , is estimated with the [O II] and [S II] emission-line doublets using the method described in Sanders et al. (2016). In that study, the [O II] and [S II] electron densities agree within the uncertainties; both reliable tracers of the density within H II regions. [O II] is the preferred choice because the doublet typically has a higher S/N. We only use [S II] measurements when we do not have adequate [O II] data. This scenario arises for the following reasons:

(i) Low S/N. An S/N ≥ 3 is required for each emission-line in the doublet to estimate n e ; however, a less restrictive requirement of S/N ≥ 3 for the overall doublet is sufficient for a galaxy to be included in the [S II] BPT and O 32 versus R 23 diagrams. We only require S/N ≥ 3 for the individual doublet members when investigating the n e of the sample. Green points indicate 1.9 ≤ z ≤ 2.7 MOSDEF galaxies with S/N ≥ 3 for all four plotted emission-lines. Large magenta stars represent measurements of stacks from composite spectra, binned by stellar mass, for 1.9 ≤ z ≤ 2.7 MOSDEF galaxies. The grayscale 2D histogram indicates local SDSS galaxies. The cyan curve is the fit to the MOSDEF sample from Shapley et al. (2015). The orange curve is a fit to the z ∼ 0 star-forming locus (Kewley et al. 2013). The red curve is the best fit to the z ∼ 2.3 galaxies from KBSS (Steidel et al. 2014). The black curve is the maximum starburst line from Kewley et al. (2001). The blue curve is the empirical AGN-/star-forming galaxy dividing line from Kauffmann et al. (2003).

(ii) Rejection of the [O II] doublet based on visual inspection. Such cases occur when skylines affect the [O II] doublet but not the [S II] doublet, spurious detections, and poor fits.

We set the lower limit for n e to be 1 cm -3 .

Throughout this study, we compare our high-redshift MOSDEF sample to local galaxies. For this comparison, we use archival data from the SDSS Data Release 7 (DR7; Abazajian et al. 2009). We obtain galaxy properties and emission-line measurements from the MPA-JHU DR7 release of spectrum measurements 1 . SDSS galaxies are selected within the 0.04 ≤ z ≤ 0.10 redshift range. We impose similar restrictions on our SDSS sample to those applied to the MOSDEF sample by requiring that each emission-line used in the analysis has an S/N ≥ 3. For the SDSS sample, we remove AGN using equation (1) from Kauffmann et al. (2003). Galaxies are also identified as having an AGN component if N2 > 0.5. These criteria result in a comparison sample of 103 422 SDSS galaxies when considering the [N II] BPT diagram alone, and 74 726 SDSS galaxies when considering galaxies with simultaneous detections across all three emission-line diagrams analyzed in this work ([N II] BPT, [S II] BPT, and O 32 versus R 23 ).

We start by investigating the locations of z ∼ 2.3 star-forming galaxies in the [N II] BPT diagram based on the complete MOSDEF sample (Fig. 2). We include the 180 galaxies with ≥3σ detections for all four emission-lines (Hβ, [O III]λλ4960,5008, Hα, and [N II]λ6585) and the corresponding SDSS sample with the same four emission lines 1 https://wwwmpa.mpa-garching.mpg.de/SDSS/DR7/ detected. Similar to previous MOSDEF studies (e.g. Shapley et al. 2015;Sanders et al. 2016;Shapley et al. 2019;Topping et al. 2020b) and other studies from the literature (e.g. Shapley et al. 2005;Erb et al. 2006;Steidel et al. 2014), there is a systematic offset observed for the high-redshift MOSDEF sample from the local sequence (Fig. 2). The MOSDEF galaxies appear on average to be shifted towards the AGN region of the diagram with elevated N2 and/or O3 values, with some galaxies on the AGN side of the Kauffmann et al. (2003) AGN/SF boundary. There are even a small number of galaxies past the maximum starburst line from Kewley et al. (2001) as well.

To search for any biases in the sample based on the ≥3σ detection requirement in all four [N II] BPT lines, we construct spectral stacks for all MOSDEF galaxies at 1.9 ≤ z ≤ 2.7 with Hα emission detected at S/N ≥ 3. There are four spectral stacks, divided into bins of stellar mass (see Sanders et al. 2018 for a full description of the method for constructing composite spectra). The emission-line ratios measured from stacked spectra follow the distribution of data points measured from MOSDEF galaxies with individual detections in all [N II] BPT lines. Therefore, the sample of individual detections seems to represent the parent z ∼ 2.3 MOSDEF data set with minimal bias. The four stacks, similarly to the individually detected galaxies, are offset from the local SDSS sample with elevated N2 and/or O3 values.

Fig. 2 shows that on the [N II] BPT diagram, part of the MOSDEF sample sits on or near the local SDSS sequence, while the remainder of the galaxies lie off the SDSS sequence and are shifted towards the local AGN region. We now investigate if the location of a high-redshift galaxyeither on or offset from the local SDSS sequenceis connected with other physical properties of the galaxy. For this analysis, we examine how galaxies on the [N II] BPT diagram populate the [S II] BPT and O 32 versus R 23 diagrams, and therefore require the sample to have ≥3σ detections for all emissionlines on these diagrams: [O II]λλ3727,3730, Hβ, [O III]λλ4960,5008, Hα, [N II]λ6585, and [S II]λλ6718,6733. There are 123 MOSDEF galaxies that meet this criterion. One galaxy was removed because it does not have a R e measurement, which brings the final sample to 122 MOSDEF galaxies.

For this analysis, we split the final sample of 122 galaxies into two groups using the functional form presented in Shapley et al. (2015), and adjusting the y-intercept so that our sample is divided in half. The equation used to split the MOSDEF sample into two groups of 61 galaxies each is

The new curve is 0.055 dex lower in [O III]λ5008/Hβ compared to equation (1) in Shapley et al. (2015). This splitting of the sample on the [N II] BPT diagram is shown in the upper left-hand panel of Fig. 3, where galaxies above the curve are indicated with blue symbols and those below with red symbols. Hereafter, we refer to the galaxies above the curve on the [N II] BPT diagram as the high sample, and the group below the curve as the low sample, following the nomenclature in Topping et al. (2020b). We plot the two samples on both the [S II] BPT diagram (upper right-hand panel of Fig. 3) and the O 32 versus R 23 diagram (bottom panel of Fig. 3). We also include binned median lines for both the high (blue) and low (red) populations. These binned medians are binned by log 10 (O3N2) for the [N II] BPT diagram, log 10 (O3)log 10 (S2) for the [S II] BPT diagram, and log 10 (O 32 ) + log 10 (R 23 ) for the O 32 versus R 23 diagram. The binning schemes above were adopted because they divide our sample into subgroups segregated roughly along the local star-forming sequence. For both the high and low samples, each of which contains 61 galaxies, there are four equally sized bins (three bins of 15 galaxies and one bin of 16).

For completeness, we checked for potential bias between the smaller sample of 122 galaxies with detections across all three emission-line diagrams and the parent sample of 180 galaxies defined based on detections in the [N II] BPT diagram alone. Of the 58 galaxies that were removed, 32 are above the curve dividing our sample (equation 8) and 26 are below it on the [N II] BPT diagram. Because these two groups are approximately equal in size, we can expect that the smaller subset with [O II] and [S II] also detected (122 galaxies) has approximately the same average BPT offset as the larger [N II] BPT only sample (180 galaxies). Therefore, we can conclude that the cuts to create the sample of 122 galaxies with detections across all three emission-line diagrams do not introduce biases when compared to the larger parent sample of 180 galaxies.

We find separation between the high and low samples on both the [S II] BPT and O 32 versus R 23 diagrams. For the [S II] BPT diagram, the high sample has a systematically higher S2 at fixed O3 (or vice versa). For the O 32 versus R 23 diagram, the high sample has a systematically higher R 23 value at fixed O 32 .

Our new results on the [S II] BPT and O 32 versus R 23 diagrams update the early MOSDEF results from Sanders et al. (2016). In this earlier work, based on a significantly smaller sample of only 53 galaxies, subsamples split in the [N II] BPT diagram were found to be well mixed in the [S II] BPT and O 32 versus R 23 diagrams. In contrast, our results are consistent with those of Strom et al. (2017), based on the KBSS survey. In that study, the sample was split into two groups based on [N II] BPT diagram locationlarge and small offset from the local SDSS sequence. These subsamples were found to remain separated on the [S II] BPT diagram. Strom et al. (2017) did not investigate if this segregation remained on the O 32 versus R 23 diagram. Section 4 contains a more in-depth comparison with the results from Strom et al. (2017) and a discussion of the implications for the underlying causes of the observed shift in the high-redshift [N II] BPT diagram.

The fact that there is a segregation in all three plots in Fig. 3 implies that there are some key differences between galaxies offset from the local SDSS sequence and galaxies that overlap with it. We now investigate several galaxy physical properties to uncover any differences between the galaxies that fall into these two categories. In this study, as described in Section 2.1, we focus on: M * , log(t/τ ) of the stellar population inferred from a delayed-τ star formation model, t of the stellar population assuming a CSF history, SFR(Hα), sSFR, SFR , R e , and n e . We have measurements of these properties for all 122 galaxies in the sample except for n e . As discussed in Section 2.1, a more stringent requirement of an S/N ≥ 3 for each component of the [O II] or [S II] doublets is needed to estimate n e while an S/N ≥ 3 for the combined doublet is needed to be plotted on the [S II] BPT and O 32 versus R 23 diagrams. Accordingly, while our fiducial sample for analysis contains 122 galaxies, only 90 of them (44 high and 46 low) have reliable n e estimates.

We show how the high and low samples divide in the space of each of these parameters using histograms, marking the median value of the parameter with solid lines for each high and low population (Fig. 4). Showing the data in this format clearly highlights how the high and low samples separately distribute in each of the galaxy parameters and, accordingly, which of the galaxy parameters are correlated with the location of galaxies on the [N II] BPT diagram. The median values for the high and low populations are listed with 1σ uncertainties derived from bootstrap resampling in Table 1. We also list the probabilities, based on the Kolmogorov-Smirnov (K-S) test, of the null hypothesis that the high and low samples are drawn from the same parent distribution.

While the distributions of the high and low samples significantly overlap for some galaxy parameters, there are several parameters for which the high and low distributions are measurably offset in the median and spread. These results indicate that high sample of MOSDEF galaxies tend to have a smaller R e and M * , but a larger sSFR and SFR compared to the low sample. We also note that even though R e , sSFR, and SFR are correlated with M * , at fixed M * the high sample is still found to be smaller with larger sSFR and SFR values compared to the low sample. The variation in median SFR between the two samples can be attributed to the difference in median R e . At fixed size, the median SFR(Hα) values for the high and low populations agree within the uncertainties. With smaller significance, the high sample is also younger (smaller median t/τ and CSF t). We find no mutual correlation for either n e or SFR with the location of galaxies on the [N II] BPT diagram.

We have divided the z ∼ 2.3 MOSDEF sample according to location in the [N II] BPT diagram. When separated in this manner, our sample also shows segregation in the [S II] BPT and O 32 versus R 23 diagrams. These results update earlier MOSDEF work (Shapley et al. 2015;Sanders et al. 2016), and are in agreement with more recent MOSDEF studies (e.g. Shapley et al. 2019;Topping et al. 2020b). Our findings for the population segregation on the [S II] BPT diagram are also in agreement with those of Strom et al. (2017), in which a similar segregation in the [S II] BPT diagram is found for z ∼ 2.3 star-forming galaxies from the KBSS survey when separated by [N II] BPT location. Strom et al. (2017) did not consider the relative positions of the KBSS equivalent of our high and low galaxies on the O 32 versus R 23 diagram. However, these authors do show how photoionization models with different input parameters (i.e. varying ionizing spectrum at fixed nebular metallicity) vary across the O 32 versus R 23 parameter space. Specifically, it is shown that models with harder ionizing spectra at fixed nebular metallicity are offset towards higher R 23 at fixed O 32 . Strom et al. (2017) additionally show that models with harder ionizing spectra produce higher O3 at both fixed S2 and N2. These separations are at least qualitatively similar to that observed between our high and low samples.

The suggested implication of these segregations on both [S II] BPT and O 32 versus R 23 diagrams is a harder ionizing spectrum at fixed O/H, relative to local galaxies (e.g. Strom et al. 2017;Shapley et al. 2019;Topping et al. 2020b). A harder ionizing spectrum at fixed nebular metallicity in z ∼ 2.3 galaxies may arise due to α-enhancement (i.e. supersolar O/Fe values) in the massive stars exciting the ionized gas in star-forming regions. Such abundance patterns may arise naturally in high-redshift galaxies, given the young median ages of their stellar populations (e.g. Steidel et al. 2016;Sanders et al. 2020a;Topping et al. 2020b).

We have also shown in Fig. 4 that the high sample is associated with more compact (i.e. smaller R e and higher SFR ), and intense (i.e. higher sSFR) star formation than the low sample. These differences between the two populations exist at fixed stellar mass as well. It will be important to explore the links between these global galaxy properties and the abundance patterns of massive stars using realistic galaxy formation simulations (e.g. FIRE-2; Hopkins et al. 2018).

We use photoionization models to explain simultaneously the joint distributions of high and low galaxies in the [N II] and [S II] BPT and is a reasonable assumption given that the high sample has a lower median stellar metallicity (i.e. a harder median ionizing spectrum and greater ionizing photon production efficiency) in addition to more concentrated star formation (i.e. higher SFR ) compared to the low sample. The high population also appears to be more α-enhanced, due to it having a lower Z * but a higher 12+log(O/H) compared to the low population. Based on the variation of best-fitting 12+log(O/H) and log(U) values between the high and low median data points, we conclude that quantifying both of these physical parameters, in addition to stellar metallicity (i.e. the hardness of the ionizing spectrum), is important when explaining the distribution of highredshift galaxies in the [N II] BPT diagram. It is worth noting that we have shown the high population to be less massive (i.e. smaller M * ) but more metal-rich (i.e. larger 12+log(O/H)) compared to the low population. This combination of physical properties is the opposite of what is expected based on results from previous studies investigating the mass-metallicity relationship (e.g. Tremonti et al. 2004;Steidel et al. 2014;Sanders et al. 2018Sanders et al. , 2020a)), which show that M * and 12+log(O/H) have a positive relationship. However, the differences in mass and metallicity between the high and low samples reflect the scatter in the mass-metallicity relationship.

To check the validity of our assumption that the N/O ratio is consistent between the high and low populations, we estimate log(N/O) for the sample median points based on the [N II]λ6585/[O II]λλ3727,3730 tracer. We use the calibration from Strom et al. (2018) to convert from log ([N II] λ6585 / [O II] λλ3727, 3730) to log(N/O), and find that log(N/O) = -1.13 ± 0.05 and -1.08 ± 0.04 for the high and low samples, respectively. The consistency in log (N/O) between the high and low samples validates our Cloudy+BPASS model assumptions and rules out N/O variations as the primary driver of the [N II] BPT offset. This conclusion agrees with results from other current MOSDEF studies (e.g. Shapley et al. 2019;Sanders et al. 2020a, b;Topping et al. 2020a) as well as other z ∼ 2 studies (e.g. KBSS; Steidel et al. 2014;Strom et al. 2017Strom et al. , 2018)).

For the most part, these results are in agreement with Topping et al. (2020b). The best-fitting 12+log(O/H) and log(U) values from both studies are in agreement within 1σ , except for the log(U) of the two high populations, which agree within 2σ . However, this study finds a larger difference in log(U) between the high and low samples, 0.20 dex, compared to Topping et al. (2020b), 0.07 dex. On the other hand, Topping et al. (2020b) find a larger variance in 12+log(O/H) between the two populations, 0.10 dex, compared to this study, 0.06 dex. These differences are due to the slightly different median restoptical line ratios for the high and low populations in Topping et al. (2020b) and in the current work.

To characterize the variation of physical properties within each of our high and low populations, we apply the same methodology of finding the best-fitting 12+log(O/H) and log(U) values, but in this case for the two endpoints of the binned medians for the high and low samples in the [N II] BPT diagram. We use the same binning method for the [N II] BPT diagram (i.e. four equally sized bins based on O3N2 strength; see Section 3.2 above for complete details on binning) when calculating the median O3, N2, S2, and O 32 emission-line ratios for the endpoints of both binned median lines. Uncertainties for the line ratios of the median endpoints are estimated using the same methodology of bootstrap resampling coupled with perturbation of the data points of individual galaxies according to their error bars, as described above. On the [N II] BPT diagram, the binned medians move primarily in the N2 direction, therefore, we refer to the endpoints of both binned median lines as 'large-N2' and 'small-N2'. It is important to note that the large-N2 endpoints of both high and low binned median lines have larger S2, and smaller O 32 compared to the small-N2 endpoints. In other words, while the binning method was not the same on all three diagrams in Fig. 3 (i.e. the four bins do not necessarily contain the same galaxies), the large-N2 bins on the [N II] BPT diagram have significant overlap with (and therefore roughly correspond to) the larger-S2, and smaller-O 32 endpoints of the binned median lines on the [S II] BPT and O 32 versus R 23 diagrams.

The results of Topping et al. (2020a) suggest that a positive correlation exists between N2 and Z * (i.e. stellar metallicity tends to be higher at larger N2 values). Therefore, we vary the Z * values that we assume for the large-N2 and small-N2 endpoints of the binned median lines accordingly. For the low population, we assume stellar metallicities of 0.003 and 0.001 (i.e. ±0.001 from the median Z * ) for the large-N2 and small-N2 points, respectively. Similarly, for the high population, we assume stellar metallicity values of 0.002 and 0.001 for the large-N2 and small-N2 points, respectively. We do not lower the stellar metallicity of the small-N2 point of the high population to our next available Z * model value (0.0001, a factor of 10 lower than the next highest Z * value) because the results from Topping et al. (2020a) do not suggest that such a significant variation in stellar metallicity is observed. It is important to note that while we do vary Z * in accordance with the results from Topping et al. (2020b), the best-fitting 12+log(O/H) and log(U) for the large-N2 and small-N2 endpoints of the high and low populations do not change significantly if we use Z * = 0.001 and 0.002 (i.e. the Z * values used for the median high and low population data points). We will discuss this systematic effect of how our choice in Z * influences the best-fitting 12+log(O/H) and log(U) in more detail below in Sections 4.2.3 and 4.2.4.

Median 12+log(O/H) and log(U) values for the large-N2 and small-N2 endpoints of the high and low populations are included in Table 2. The corner plots showing the distributions of 12+log(O/H) and log(U) for the large-N2 and small-N2 endpoints are shown in the top right-and bottom-hand panels, respectively, of Fig. 6. In both samples, we find an anticorrelation between 12+log(O/H) and log(U) along the binned median lines. Moving from large-N2 to small-N2 (therefore also large-S2 to small-S2 and small-O 32 to large-O 32 ) leads to a lower 12+log(O/H) but a higher log(U). This anticorrelation between log(U) and 12+log(O/H) has been observed in local starforming galaxies and H II regions (Pérez-Montero 2014). In addition, this trend is consistent with the idea that the high O 32 tail of the local sequence on the O 32 versus R 23 diagram includes galaxies with low metallicity, and high ionization parameter (Shapley et al. 2015). Also, on the [N II] BPT diagram both stellar and nebular metallicity decrease along the local sequence (i.e. from low O3 and high N2 to high O3 and low N2).

As described above, thus far we have assumed Z * values of 0.001 and 0.002 for the high and low population median points to match the best-fitting values from modelling the rest-frame UV spectra and small-N2 endpoint of the sample binned medians (bottom centre) are included. We assume Z * values of 0.001 and 0.002 (median points; top left-hand panel), 0.002 and 0.003 (large-N2 points; top right-hand panel), and 0.001 and 0.001 (low-N2 points; bottom panel) for the high and low populations, respectively. The dashed lines on the histograms mark the median value and 1σ uncertainties (16th and 84th percentiles) for the distributions of 12+log(O/H) and log(U), and the three contours on the diagonal mark the 1σ , 1.5σ , and 2σ regions in the 2D 12+log(O/H)-log(U) parameter space. It is shown that moving from large-to small-N2 corresponds to decreasing nebular metallicity and increasing ionization parameter. In each of the three distributions, the high population has a higher 12+log(O/H) and log(U).

of z ∼ 2.3 MOSDEF galaxies (Topping et al. 2020b). While this is a reasonable assumption, we also investigate how relaxing our requirement on stellar metallicity affects the best-fitting nebular metallicity and ionization parameter values of the two populations. For this additional analysis, we use the same χ 2 method described above (i.e. fitting the N2, S2, O3, and O 32 emission-line ratios). However, we not only fit for 12+log(O/H) and log(U), but also treat Z * as a free parameter. We investigate two different model grids: one that allows a small range of Z * for the high and low samples, including the best-fitting stellar metallicities from Topping et al. (2020b) for each population and those that are adjacent to the best-fitting values, and a fully unrestricted method that allows all Z * values used by the BPASS models ranging from 0.00001 to 0.03. This analysis shows that when we relax the constraints provided by rest-UV spectra (both with limited freedom and in a fully unconstrained manner), the rest-optical emission lines tend to favour higher stellar metallicities. When we allow for limited freedom in Z * the stellar metallicities for the high and low samples prefer the largest allowed Z *i.e. 0.002 and 0.003, respectively. When all stellar metallicity values are allowed, the high and low populations favour Z * values about four to five times greater than those found by Topping et al. (2020b) based on fitting rest-UV spectra. These shifts demonstrate the importance of imposing external constraints from a more direct probe of massive stars and stellar metallicity, i.e. the rest-UV continuum, as opposed to relying only on a joint fit of rest-optical emission lines. We also note that, while the best-fitting values of Z * increase when stellar metallicity is allowed to vary, the low population still always favours a higher Z * than the high population.

The best-fitting 12+log(O/H) and log(U) values for the single-Z * and limited freedom Z * models are significantly consistent (i.e. within 1σ ) for both the high and low populations. For the high sample, the best-fitting values for the single-Z * and fully unconstrained Z * models are significantly consistent as well. For the low sample, only some of the best-fitting values for the single-Z * and fully unconstrained Z * models are consistent within 1σ ; however, all values are consistent within 2σ . In addition, similar to what we find for the single-Z * models, we find that O/H increases and log(U) decreases with increasing N2 in both free-Z * models.

The results from this section provide guidance on which wavelength regimes best constrain different galaxy properties. Supported by other studies (e.g. Topping et al. 2020b), rest-UV spectra are essential for accurately estimating Z * . Rest-optical spectra, even with the combination of multiple emission-line ratios, cannot constrain the stellar metallicity. At high redshift, we have shown that modelling rest-optical spectra alone will lead to overestimates of stellar metallicity relative to more accurate methods that incorporate rest-UV spectral information. However, rest-optical spectra can be used to constrain 12+log(O/H) and log(U). Because the best-fitting values of nebular metallicity and ionization parameter vary based on which rest-optical emission-line diagram is used (see Fig. 5), the combination of line ratios from multiple diagrams will increase the accuracy of the models. These parameters can even be estimated with reasonable precision without detailed knowledge of the ionizing spectrum (i.e. the metallicity of the stellar population).

In summary, using rest-optical emission-line ratios alone without the constraints from fitting rest-UV spectra, we infer systematically higher Z * values for both high and low samples. At the same time, when multiple rest-optical emission lines ratios are measured (N2, S2, O3, and O 32 ), the inferred median nebular parameters (12+log(O/H) and log(U)) and their variation across the BPT diagram, do not depend strongly on the allowed range of Z * . However, the inclusion of constraints on Z * from rest-UV fitting is essential for our understanding of the abundance patterns (i.e. α-enhancement) of z ∼ 2.3 galaxies and we adopt the constrained values of Topping et al. (2020b) for our fiducial modelling procedure (Table 2). Because we see variation in Z * , 12+log(O/H), and log(U) at the median and large-and small-N2 endpoints for the high and low samples, we conclude that constraints on all three of these physical properties are required for fully understanding the observed distribution of z ∼ 2.3 star-forming galaxies in the [N II] BPT diagram. The importance of α-enhancement has also been highlighted by Steidel et al. (2016) and Cullen et al. (2019).

It is also important to note that while the high sample is more α-enhanced with a harder ionizing spectrum compared to the low sample, both populations are α-enhanced and both are characterized by harder ionizing spectra when compared with their equivalents (i.e. galaxies with similar rest-optical line ratios) at low redshift. Notably, such physical differences apply to the low population, even though it overlaps with the local [N II] BPT sequence. Specifically, an overlap in emission-line ratios does not correspond to the same properties for the ionizing radiation field and ionized gas. Therefore, we recommend determining the location of galaxies in multiple restoptical emission-line diagrams, as well as having some knowledge of the ionizing spectrum, when trying to identify local systems that can be used as true z ∼ 2 analogues.

In addition, the differences between the high and low populations found in both this section and in Section 3.2.2 suggest that z ∼ 2.3 galaxies in the high sample that generally have a lower Z * with higher 12+log(O/H) and log(U) values are also typically smaller with higher sSFR and SFR values. Looking ahead, we must understand why there exists a connection between the local properties of ionized star-forming regions and the overall compactness/intensity of star formation. However, exploring this phenomenon is outside the scope of the current work.

It is worth noting that in all three diagrams in Fig. 5, there are galaxies that do not lie within the Cloudy+BPASS grid space. On the [N II] BPT and [S II] BPT diagrams, there are galaxies with O3 values that exceed the extremes of the model grids. Similarly, on the O 32 versus R 23 diagram, we observe MOSDEF galaxies with higher R 23 values than what the Cloudy+BPASS models predict. These offsets are primarily found among the high population. Such extreme galaxies are not specific to this study, as Strom et al. (2017) find similar discrepancies between data and models in the corresponding emission-line diagnostic diagrams.

For the Cloudy+BPASS model grids in this study, two possible solutions could be to either elevate n e or lower Z * . Increasing n e has been shown to cause elevated O3 and N2 values (Kewley et al. 2013). Since the high population (defined by its elevated O3 and/or N2 values) also has elevated R 23 values on average, increasing n e in our Cloudy+PASS models theoretically could push the curves to higher R 23 . The models in this study, which were constructed by Topping et al. (2020b), set n e = 250 cm -3 , which is typical of z ∼ 2.3 galaxies in both the KBSS sample (Strom et al. 2017) and early MOSDEF work (Sanders et al. 2016). It is also comparable with the median electron densities reported in this work. The models employed in Strom et al. (2017) use n e = 300 cm -3 , and similarly cannot reach the highest observed O3 and R 23 values in the KBSS survey.

Using Z * = 0.001, we explore raising n e to 1000 cm -3 in the Cloudy+BPASS models (right-hand column of Fig. 7). Comparing the elevated electron density models with the Z * = 0.001, n e = 250 cm -3 scenario (left-hand column of Fig. 5) shows the differences in the predicted emission-line ratios caused by isolating variations in n e on the [N II] BPT, [S II] BPT, and O 32 versus R 23 diagrams. We find that elevating n e to such an extreme value has only a small effect on the emission-line ratios, and the models still do not encompass the galaxies observed at the highest N2, S2, O3, and R 23 values.

A more promising approach is to lower the stellar metallicity. The left-hand column of Fig. 7 gives the model grids for Z * = 0.0001, which reach the majority of the extreme data points. Therefore, a very low Z * (i.e. a very hard ionizing spectrum) is the more likely than variations in n e to be cause for these outliers.

In addition to a low Z * , other effects may also lead to a harder ionizing spectrum (e.g. variations in the IMF slope and high-mass cutoff or AGN partially contributing to the integrated emission lines). Systematic uncertainties such as the BPASS models underpredicting the hardness of the ionizing spectrum at a given Z * (particularly at lower values of Z * ), different star formation histories, or uncertainties in the dust corrections applied to O 32 and R 23 could play a role in the observed discrepancies between models and a minority of the sample as well. A comparison of these models with the corresponding panels in Fig. 5 shows how raising the electron density and lowering stellar metallicity affects the predicted emission-line ratios. The correspondence between colours and symbol size and, respectively, 12+log(O/H) and log(U), matches that in Fig. 5. Once again, we adopt an age of 10 8.6 yr because it is the approximate median age of the MOSDEF sample (given a CSF history) and include the local SDSS sample (grey 2D histogram), the high (blue) sample, and low (red) sample with binned medians as shown in Fig. 3. The green and yellow squares, with associated uncertainties, are the median values of the high and low populations as shown in Fig. 5. Aside from at roughly solar and supersolar nebular metallicities, a higher electron density has a minimal affect on the Cloudy+BPASS model grids and cannot reach the elevated O3 and O 32 galaxies. A lower stellar metallicity is able to reach these values, suggesting that the most offset galaxies have the hardest ionizing spectrum.

We present results on the emission-line properties of a sample of 122 star-forming galaxies at 1.9 ≤ z ≤ 2.7 from the complete MOSDEF survey with ≥3σ detections for the [O II]λλ3727,3730, Hβ, [O III]λλ4960,5008, Hα, [N II]λ6585, and [S II]λλ6718,6733 emission lines. To investigate the observed systematic offset of z > 1 star-forming galaxies on the [N II] BPT diagram relative to local systems, the MOSDEF sample is split into the high (offset with elevated O3 and/or N2 values) and low (overlapping with the local SDSS sequence) samples. We compare the location of galaxies in both populations on the [S II] BPT and O 32 versus R 23 diagrams. Additionally, we compare physical properties -SFR, M * , R e , sSFR, SFR , n e , and t (i.e. stellar population age)of the galaxies in both groups using additional CANDELS and 3D-HST ancillary data. Finally, we use Cloudy+BPASS photoionization models to investigate potential differences in the stellar metallicity, nebular metallicity, and ionization parameter of the high and low groups.

The main results are as follows:

(i) On the [S II] BPT and O 32 versus R 23 diagrams, the high sample is offset on average towards higher S2 at fixed O3 (or vice versa) and higher R 23 value at fixed O 32 , respectively, relative to the low sample. These results update earlier work from MOSDEF (Sanders et al. 2016) and are consistent with results from Strom et al. (2017), based on the KBSS survey.

(ii) The high sample has a smaller median R e and M * , but higher median sSFR and SFR compared to the low sample. The observed differences in R e , sSFR, and SFR are present at fixed stellar mass. The high sample is also slightly younger in median age (t); however, this age difference is not significant within the uncertainties. There is not a significant variation in SFR or n e for the two populations. These results imply that the high population is associated with more concentrated and intense star formation than the low population.

(iii) Using photoionization models, we find that the high population has larger 12+log(O/H) and log(U) values compared to the low population. These conclusions hold both when we use the results from Topping et al. (2020b) and assume Z * values of 0.001 and 0.002 for the high and low samples, respectively, as well as when we treat Z * as a free parameter. Also, while the lower Z * /Z neb ratio of the high population implies that this sample is more α-enhanced than the low population, both samples are significantly α-enhanced relative to local star-forming galaxies with similar emission-line ratios assuming typical local galaxies have roughly solar α/Fe. These differences in the properties of the ionizing radiation field are critical to include when describing high-redshift galaxies -even those that overlap the local emission-line sequences.

(iv) Combining the results from the median galaxy properties of the high and low samples with the Cloudy+BPASS model grids leads to the conclusion that galaxies with a harder ionizing spectrum are associated with smaller sizes and higher sSFR and SFR .

(v) A harder ionizing spectrum at fixed nebular metallicity is favoured as one of the key drivers of the [N II] BPT offset. In this study, we find that variation in multiple physical parameters drives the variation in emission-lines in the MOSDEF sample, including Z * , 12+log(O/H), and log(U). However, even in regions where the z ∼ 2.3 galaxies overlap local ones in the space of emission-line ratios, the inferred underlying physical parameters for z ∼ 2.3 galaxies are distinct from those of local galaxies. The full set of these parameters must be understood through modelling both rest-UV and rest-optical spectra, in order to understand the translation between empirical emission-line ratios and key physical quantities such as the gas-phase oxygen abundance.

Understanding the global properties of high-redshift galaxies using large samples at z ∼ 2.3 is essential when comparing them with local galaxies. The use of large high-redshift samples has enabled us to discover more subtle trends in the MOSDEF data that were not found with more limited data sets (e.g. Shapley et al. 2015;Sanders et al. 2016). The work from this study has shown that galaxies with a harder ionizing spectrum also tend to be smaller in size with more intense and compact star formation. Why these factors are linked is not yet clear, and combining large z ∼ 2.3 samples of star-forming galaxies with realistic galaxy formation simulations (e.g. FIRE-2; Hopkins et al. 2018) will be important in finding this connection.

It has also been shown in both this work and in Strom et al. (2017) that models are not yet able to reach the highest observed O3 and R 23 values for z ∼ 2.3 star-forming galaxies on the [N II] BPT, [S II] BPT, and O 32 versus R 23 diagrams. The local sequence does not reach such elevated O3 and R 23 values, and unlike their highredshift counterparts, do not greatly exceed the data space covered by the model grids. This discrepancy indicates a need to develop both stellar population and emission-line models that better match the now large data sets of z ∼ 2.3 star-forming galaxies.