MAMMOTH-MOSFIRE: Environmental Effects on Galaxy Interstellar Medium at z ∼ 2

In HST cycle 28, we conduct the MAMMOTH-Grism slitless spectroscopic survey (GO-16276, PI: X. Wang) to obtain deep WFC3 G141 exposures in these three massive protocluster fields (X. Wang et al. 2022b). The core region of each protocluster is covered by five individual WFC3 pointings, each at three-orbit depth of G141 slitless spectroscopy. The corresponding HST footprints are shown in Figure 1. These observations are designed to capture key emission lines from galaxies at z ∼ 2.2, notably [O III], Hβ, Hγ, and [O II], efficiently confirming protocluster membership. The G141 spectrum of one example galaxy is shown in Figure 2. Using these G141 data, we identify ∼50 emission-line galaxies in each field as protocluster members.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Using images from HST and ground-based telescopes, we estimate the stellar masses of these protocluster member galaxies with BAGPIPES (A. C. Carnall et al. 2018), adopting the initial mass function from G. Chabrier (2003) and the extinction law from D. Calzetti et al. (2000). Our results of mass are shown in Table 3. Further information will be forthcoming in Y. Yang et al. (2025, in preparation).

Table 2. Details of Our Keck MOSFIRE Spectroscopy of the Protocluster Member Galaxies in the BOSS1244, BOSS1441, and BOSS1542 Fields

Field/Slit Mask	Nobja	R.A.	Decl.	Position Angle	Filter	Number of Repeatsb	Observation Date
		(J2000.0)	(J2000.0)	(deg)		(s)
BOSS1244/mask1.3.3	6/22	12h43m3638	35d52m4573	331	K	10	04/15/2022
BOSS1244/mask2.2.3	7/23	12h43m2944	35d56m0233	133	K	11	04/15/2022
BOSS1244/mask5.3	9/26	12h43m3199	35d54m1880	11.5	K	11	05/06/2023
BOSS1244/mask6.6	7/18	12h43m2451	35d55m1484	344	K	8	05/23/2023
BOSS1441/mask2.6	5/22	14h41m3081	40d02m4435	70	K	10	05/06/2023
BOSS1542/mask1.3.1	5/22	15h42m4605	38d58m5116	354	K	13	04/15/2022
BOSS1542/mask2.2	3/22	15h42m5017	38d51m5513	355	K	10	05/23/2023
BOSS1542/mask3.3	1/17	15h42m5086	38d51m5544	309	K	13	05/23/2023

Notes. The observations were taken in 3 full nights from 2022A_U016 and 2023A_U139 (PI M. Malkan). ^aThe number on the right refers to the total number of objects put on each mask, whereas the number on the left corresponds to the number of Hα emitters analyzed in this work. ^bOne repeat of the exposure set comprises one standard ABA’B’ dither pattern, with each exposure having an integration time of 180 s, excluding overheads. This is the recommended exposure strategy for K filter spectroscopy, with the default slit width of 07 and the MCDS16 sampling mode.

Download table as:  ASCIITypeset image

Table 3. Measured Properties of Individual Sources in Galaxies Residing in Extreme Overdensities at z ∼ 2.2

Protocluster Field	ID	zspec					log(ne)	Stellar Mass	12+log(O/H)b	SFR	[S II]
							(cm−3)			(M* yr–1)	λ6716/λ6731
BOSS1244	gal00764	2.22				–a		10.52	–
BOSS1244	gal01335	2.24				–		10.28
BOSS1244	gal01394	2.22						9.65
BOSS1244	gal01998	2.21						10.45
BOSS1244	gal02474	2.19				–	–	9.93
BOSS1244	gal02570	2.25					–	10.4
BOSS1244	gal02611	2.23				–		9.72
BOSS1244	gal02771	2.24	–	–			–	9.79
BOSS1244	gal90044	2.24	–					10.1
BOSS1244	gal01467	2.21						10.15
BOSS1244	gal01855	2.25					–	10.76
BOSS1244	gal01890	2.25		<−1.84		–	–	9.99
BOSS1244	gal02992	2.35					–	10.31
BOSS1244	gal03061	2.23					–	10.21
BOSS1244	gal03331	2.23		<−1.74			–	9.95
BOSS1244	gal00306	2.22	⋯					10.96	–			–
BOSS1244	gal00327	2.23					–	9.86			–
BOSS1244	gal00613	2.23					–	10.1
BOSS1244	gal01417	2.25					–	10.3	–
BOSS1244	gal01878	2.06				–		10.12
BOSS1244	gal02327	2.24				–		10.21
BOSS1244	gal02588	2.25						10.36
BOSS1244	gal02950	2.23				<−1.23	–	10.53			–
BOSS1244	gal03040	2.25	⋯	<−2.71			–	10.16	–		–	–
BOSS1244	gal03340	2.23						10.18
BOSS1244	gal01070	2.23		⋯			⋯	9.96
BOSS1244	gal02510	2.23						10.34
BOSS1441	gal03745	2.25				–	–	10.4
BOSS1441	gal04511	2.21						10.1
BOSS1441	gal04563	2.31	<1.2				–	9.89
BOSS1441	gal04794	2.34					–	10.64			–
BOSS1441	gal04932	2.31				–	–	10.2
BOSS1542	gal04669	2.27	–					10.2	–			–
BOSS1542	gal04816	2.35	–				–	9.12	–			–
BOSS1542	gal05204	2.20				–	–	9.52
BOSS1542	gal90082	2.23	–					10.1
BOSS1542	gal01338	2.25					–	10.2
BOSS1542	gal02830	2.23	–					9.61	–			–
BOSS1542	gal02872	2.24	–			–	–	10.24	–			–
BOSS1542	gal04010	2.25	–			–		10.59
BOSS1542	gal01145	2.25	–	–		–	–	10.06	–			–
BOSS1542	gal01704	2.34	–					10.06
BOSS1542	gal03125	2.31					–	10.1

Notes. ^aA dash (–) indicates that the line is too weak to be detected, and thus no corresponding value is available. ^bThe metallicity data are based on the (F. Bian et al. 2018) calibrations.

Download table as:  ASCIITypeset images: 1 2

We follow our previous work to analyze the HST WFC3/G141 slitless spectroscopic data using the GRIZLI¹⁷ software (X. Wang et al. 2019, 2020, 2022a, 2022b). GRIZLI performs forward modeling of paired grism and pre-imaging exposures, fitting the optimally extracted one-dimensional (1D) spectra with linear combinations of spectral templates to determine the best-fit grism redshifts. This procedure is particularly robust for our targets at the protocluster redshifts (z ∼ 2.2–2.3), owing to the presence of strong nebular emission lines ([O III], Hβ, Hγ, and [O II]). We therefore select galaxies with secure grism redshifts that are spectroscopically confirmed members of the protoclusters. Among the confirmed protocluster members, we further prioritize galaxies with significant emission-line detections—specifically signal-to-noise ratio (SNR) ≳ 2 in [O III], Hβ, and [O II]—by assigning higher target weights during MOSFIRE slit-mask design (see Section 2.3).

To acquire K-band spectroscopy (1.91–2.42 μm) of the confirmed protocluster member galaxies, we carry out the MAMMOTH-MOSFIRE program using the MOSFIRE instrument on the Keck I telescope under observing programs 2022A_U016 and 2023A_U139 (PI: M. Malkan). These MOSFIRE K-band data cover multiple key rest-frame optical nebular emission lines (e.g., Hα, [N II], [S II], [O I]) of our protocluster member galaxies at z ∼ 2.2–2.3, highly complementary to those secured by our HST G141 data. In Table 2, we list the detailed information about Keck observations. To trace our target galaxies, we use the MOSFIRE Automatic GUI-based Mask Application (MAGMA) software, which is designed to configure observation masks and execute the corresponding slit configurations at the telescope. The spectroscopic observation regions, covered by Keck MOSFIRE masks, are indicated by red rectangles () in each panel of Figure 1 within the central density peaks of the three massive protocluster fields. The background maps and contours represent the density distribution, while the HAE candidates, identified through ground-based broad and narrow band imaging, are marked by open black circles.

The MOSFIRE data reduction was performed using the PYPEIT package (J. X. Prochaska et al. 2020a). For wavelength calibration, PYPEIT software uses the OH atmospheric lines instead of arc lamp exposures. The 1D spectra from each exposure were then co-added to enhance the SNR. The sensitivity function was built using a standard star to calibrate the 1D spectra of galaxies.

Emission line ratios are commonly used to study the ionization mechanisms of the ISM (J. A. Baldwin et al. 1981). Measuring these ratios requires accurate determination of emission line fluxes, which in turn depends on identifying the observational wavelengths of the emission lines. This identification relies on reliable redshift measurements for each galaxy. While the protoclusters (BOSS1244, BOSS1441, BOSS1542) have a general redshift of z ∼ 2.2–2.3, the exact redshifts of individual galaxies remain undetermined. Therefore, it is essential to first determine the redshifts of the 43 galaxies.

We determine the redshifts of the galaxies using the peak pixel wavelength of the Hα and [N II] emission lines. Each emission line is modeled using a Gaussian function constructed with the Python package LMFIT (M. Newville et al. 2016). To ensure reliability, we manually inspect the Hα peaks of all galaxies to exclude interference from sky lines or other emission lines. The final redshifts of the galaxies are listed in Table 3.

With the redshifts determined, the specific wavelengths of the [O I], [N II] doublet, Hα, and [S II] doublet emission lines are calculated. The fluxes and uncertainties of these lines are measured by fitting a Gaussian function to each emission line. We performed 1000 iterations of Gaussian fitting based on the uncertainties of the spectra, obtaining both the Gaussian model and the associated errors for each emission line. The median SNR of the Hα emission line fluxes is 31.8, while the uncertainties on emission line flux are based on the 68th percentile confidence intervals, representing a typical observational noise level.

We use the Python package SPECUTILS (Astropy-Specutils Development Team 2019) to fit the continuum of the entire spectrum with a linear model. The emission line component is isolated by subtracting the continuum model from the observed spectrum. The uncertainties of the emission line components are identical to those of the observed spectrum.

For emission-line fitting, we employed six three-parameter Gaussian functions, with their centers corresponding to the wavelengths of the pixels with the peak for the six emission lines. During the fitting process, we constrained the Gaussian widths (σ) of the doublets to be consistent. For example, [N II]λ6548 and [N II]λ6584 share the same σ, as do [S II]λ6717 and [S II]λ6731. We performed 1000 iterations for each Gaussian fit based on the spectral errors. The mean value of the integrated flux from these iterations was taken as the emission line flux, and the standard deviation was used as its uncertainty. The line ratios log([O III]/Hβ), log([N II]/Hα), log([S II]/Hα), and log([O I]/Hα) are listed in Table 3.

We measured the FWHM of the Hαemission lines in galaxy spectra and found the distribution value of 188 ± 55 km s⁻¹. When analyzing G141-band spectra with GRIZLI, a key advantage is that GRIZLI accounts for morphological broadening, which fixes the grism FWHM at about 1100 km s⁻¹. In fitting the grism spectra, we adopt an intrinsic 1D line profile.

We selected galaxies based on the following spectroscopic criteria:

1. The galaxy’s Hα emission line flux has an SNR of at least 20.

2. The redshift falls within the range z ∼ 2.2–2.3 to ensure that the primary emission lines, such as Hα, [S II] and [N II], are within the K-band, and the G141 grism can cover emission lines including [O III], Hβ and [O II].

Figure 3 presents example spectra near the Hα wavelength region for 1 confirmed HAE out of a total of 43. The left panel shows the 1D and 2D spectra of galaxy ID gal02588 in the BOSS1244 field covering emission line Hα and [N II], while the right panel shows emission line [S II]. The bottom panels display the 2D spectra, and the top panels show the corresponding 1D spectra, with the continuum subtracted. The blue shaded region indicates the error range of the observed spectrum, and the red lines represent the emission line models fitted using the Gaussian function, as described in Section 2.4. The 2D spectra clearly show Hα emission lines, which serve as one of the criteria for identifying HAEs. Galaxy ID gal02588 exhibits strong Hα emission along with clearly detected [N II] lines in both the 2D and 1D spectra.

Zoom In Zoom Out Reset image size

To investigate differences in the ionization mechanisms of ISM across galaxies with varying stellar masses, we analyzed both individual galaxy spectra and stacked spectra in different stellar mass bins. Stacked spectra enhance the SNR, enabling more precise measurements of weak emission line fluxes. For example, the [O I] emission line is crucial when studying ISM ionization mechanisms, as it can probe the nature of low-ionization regions and processes such as shocks. However, since it is a weak line, the [O I] signal in individual galaxy spectra is often faint. By stacking spectra within different stellar mass intervals, the improved SNR allows for more robust measurements of average [O I] flux. We divide these galaxies into two mass bins: [9.0, 10.0), and [10.0, 11.0]. Based on the stellar masses listed in Table 3, 12 galaxies fall into the low-mass bin, while 31 galaxies are classified into the high-mass bin.

We then adopt the following stacking procedures, following the method described in X. Wang et al. (2022b),

1. 1.  
   Subtract continuum from the observed spectra. The K-band continuum spectrum is modeled by fitting a polynomial using the SPECUTILS (R. J. Sault 1994; Astropy-Specutils Development Team 2019) software, while the G141-band continuum spectrum is constructed with the GRIZLI software.
2. 2.  
   Normalize the spectrum of each galaxy by the measured Hα emission line flux.
3. 3.  
   Mask bright sky lines from the observed spectra to minimize contamination.
4. 4.  
   Apply redshift correction to recover the rest-frame spectral wavelengths.
5. 5.  
   Recreate the stacked spectra through nonparametric bootstrap resampling, drawing 1000 samples to represent the overall data set. We take the median spectrum as the final stacked spectrum and and the 1σ uncertainty range is defined by the 16th and 84th percentiles.

We additionally stack the spectra of all 43 galaxies.

We fit the emission lines in the three stacked spectra to derive the line fluxes, following the method described in Section 2.4. The stacked spectra and the emission line models are shown in Figure 4. From top to bottom, the panels show the stacked spectrum for the low-mass bin, the high-mass bin, and the 43 full stacked galaxies. The blue lines represent the spectra, while the blue shaded regions indicate the 1σ error ranges. The red dashed lines correspond to the emission line models.

Zoom In Zoom Out Reset image size

Notably, the Hα lines are prominent in all three spectra, whereas the [O I], [N II], and [S II] lines are clearly visible in the high-mass bin and the total stacked spectra.

Electron density is a key parameter in determining the intensity of collisionally excited forbidden emission lines. The flux ratios of density-sensitive emission lines, such as the [O II] and [S II] doublets, can be used to estimate the electron density in both local galaxies and those at high redshifts (e.g., N. Bennert et al. 2006; R. L. Sanders et al. 2015). We note that our spectral resolution is sufficient to fully resolve the [S II] doublet components, but not those of the more closely spaced [O II] doublets.

R. L. Sanders et al. (2015) investigated the electron densities of star-forming galaxies at z ∼ 2.3 in the MOSFIRE Deep Evolution Field survey. They adopted effective collision strengths calculated for an electron temperature of 10,000 K, which represents the typical equilibrium temperature of H II regions. Although assuming a fixed electron temperature of 10,000 K may lead to an overestimation of electron density in metal-rich galaxies and an underestimation in metal-poor galaxies, the uncertainty introduced by this assumption is smaller than the typical measurement uncertainties for high-z galaxies. We follow the standard method described in R. L. Sanders et al. (2015) to measure electron density:

where R = [S II]λ6717/[S II]λ6731. The constants are a = 0.4315, b = 2.107, and c = 627.1.

We calculate the electron densities for all 43 spectra and stacked spectra. However, since the [S II] lines are relatively weak, we ultimately determine the electron densities for only a subset of the sample.

The results, including uncertainties, are presented in Table 3, with stacked results summarized separately in Table 4. We find that the electron density in the full stacked spectrum is about 290 cm⁻³, the high-mass stacked sample is 93 cm⁻³, and the low-mass stacked sample is 1070 cm⁻³. Individual galaxies exhibit electron densities ranging from 37 to 2410 cm⁻³.

Table 4. Measured Properties of the Stacked Spectra

Mass Bin	Ngal	log(M*/M⊙)					log(ne)	12+log(O/H)a	[S II]
							(cm−3)		λ6716/λ6731
Low mass	12	[9.0-10.0]
High mass	31	[10.0-11.0]
Full stack	43	…

Note. ^aThe metallicity data are based on the (F. Bian et al. 2018) calibrations.

Download table as:  ASCIITypeset image

The average electron density in our galaxies is consistent with previous studies. R. L. Sanders et al. (2015) reported electron densities of 260 and 291 cm⁻³, similar to the results of R. Shimakawa et al. (2015) for the z = 2.5 protocluster, and 10 times higher than in local galaxies. The low-mass stacked galaxies have high electron density and the high-mass stacked galaxies have low electron density. Additionally, C. C. Steidel et al. (2014) estimated [S II]ratio to measure an electron density of 243 cm⁻³ from their stacking analysis of 113 spectra at z ∼ 2.3 using Equation (7).

There are various methods to measure gas-phase oxygen abundance, including the “direct” (T_e) method and the empirical metallicity calibrations derived from strong emission lines (e.g., X. Liu et al. 2008; S. F. Newman et al. 2013; C. Ly et al. 2014; C. C. Steidel et al. 2014). However, these methods are usually calibrated for local galaxies and may not be directly applicable to high-redshift galaxies. For our z ∼ 2 galaxies, we adopt the “direct” metallicity method using O32 (F. Bian et al. 2018). F. Bian et al. (2018) utilized stacked spectra of star-forming galaxies at z ∼ 2, detecting the weak auroral line [O III]λ4363 to determine the direct electron temperature (T_e). Based on this approach, they established an empirical relation between the direct oxygen abundance and strong-line ratios, such as O32. This calibration, particularly relevant for high-redshift star-forming galaxies, is:

where .

The metallicities of 43 galaxies and those derived from the stacked spectra are given in Tables 3 and 4, respectively. For seven galaxies, the lack of emission lines required to measure the O32 ratio prevented a metallicity determination. The metallicities of the remaining galaxies are relatively similar, spanning a range of 8.15–8.60. The metallicities derived from the stacked spectra are also consistent across different mass bins, with the low-mass bin showing a metallicity of 8.27 and the high-mass bin yielding 8.29.

We show the mass–metallicity relation of our sample in Figure 5. We observe that our sample galaxies are located below the MOSDEF sample (R. L. Sanders et al. 2021), suggestive of more efficient gas dilution in our protocluster environments (Z. Li et al. 2022; X. Wang et al. 2022b). Additionally, we include the protocluster sample from Yang et al. (2025, in preparation), which also shows that protocluster galaxies exhibit lower metallicity than normal field galaxies. By fixing the intercept and fitting our sample against the protocluster field galaxies, we find a difference of 0.03 in the slope compared to the samples from Yang et al. (2025, in preparation). This suggests a minor selection bias in the protocluster sample selection.

Zoom In Zoom Out Reset image size

We find that 80% of the galaxies in the sample have metallicities below 8.5. This preference for a metal-poor ISM aligns with the expectation that star-forming galaxies at z ∼ 2 have lower metallicities compared to their local counterparts.

The ionization parameter (q) is a fundamental physical quantity defined as the ratio of the ionizing photon flux per unit area () to the local hydrogen number density (n):

where represents the ionizing photon flux, and n is the local hydrogen number density. q measures the density of ionizing photons per unit volume and can also be interpreted as the maximum velocity of an ionization front driven by the local radiation field.

Using photoionization models, q can be derived from the ratio log([O III]λ5007/[O II]λλ3727, 3729), which is highly sensitive to the ionization parameter, as it reflects the relative populations of ions in different ionization states (L. J. Kewley & M. A. Dopita 2002). The relationship between [O III]/[O II] and q is parameterized through polynomial fitting, with coefficients that depend on metallicity. The metallicity is first estimated, and an iterative approach is used to simultaneously solve for q and metallicity:

where R is the flux ratio log([O III]/[O II]), and k₀, k₁, and k₂ are coefficients dependent on metallicity. The parameters corresponding to different metallicities are listed in Table 5. We adopt this relation to calculate q.

Since 80% of the galaxies in our sample have metallicities 12 + log(O/H) < 8.5, we compute the ionization parameter assuming three metallicities: 0.1 Z_⊙, 0.2 Z_⊙, and 0.5 Z_⊙. Our analysis reveals the following:

1. 1.  
   At 0.1 Z_⊙, 80% of the galaxies have log(q) > 7.5.
2. 2.  
   At 0.2 Z_⊙, 100% of the galaxies exhibit log(q) > 7.5.
3. 3.  
   At 0.5 Z_⊙, 100% of the galaxies have log(q) > 7.5.

For the stacked spectra, we applied the same metallicity assumptions and calculated the ionization parameters. We find consistently high log(q) values, log(q) > 7.75, for all cases. These findings indicate that protoclusters exhibit elevated ionization parameters. This enhanced ionization parameter is similar to what is found in field galaxy spectra at z = 1.5 (M. Hayashi et al. 2015; A. J. Battisti et al. 2024).

In addition to electron density, metallicity, and ionization parameter, the BPT diagrams and the O32–R23 excitation diagram (O3R2 diagram) are critical related tools for examining the ionization properties of a galaxy’s ISM. The O32 ratio () and R23 ratio () are highly sensitive to both metallicity and ionization parameter (L. J. Kewley & M. A. Dopita 2002; A. E. Shapley et al. 2015).

High-redshift galaxies (z ∼ 2.2) have higher O32 and R23 values compared to local galaxies, occupying the tail region of the local distribution characterized by low metallicity and high ionization parameters (A. E. Shapley et al. 2015). The high ionization parameter is defined based on our model grid and is illustrated in Figure 6, where less than half of the SDSS sample lies above this value. Accordingly, we consider ionization parameters above 7.5 to be high. Since our grid spans an ionization parameter range from 6.5 to 9.0, we adopt 7.5 as the threshold for a high ionization parameter.

Zoom In Zoom Out Reset image size

In this work, we use the O3R2 diagram to examine the metallicities and ionization parameters of the sample galaxies, as shown in Figure 6. The gray points represent local galaxies from the Sloan Digital Sky Survey (SDSS) Data Release 7 (DR7) catalog (K. N. Abazajian et al. 2009), while the cyan points denote the 43 HAEs identified in this work. The red, purple, and blue rhombic symbols correspond to the stacked spectra: the high-mass stacked spectrum, the low-mass stacked spectrum, and the full stacked spectrum, respectively. The color bar in the figure indicates metallicities, which were derived by F. Bian et al. (2018).

As shown in Figure 6, the high-redshift galaxies have elevated O32 and R23 values compared with local galaxies. This indicates that galaxies in protoclusters possess higher ionization parameters and lower metallicities. Further analysis reveals that low-mass galaxies display higher R23 values compared to high-mass galaxies, suggesting that low-mass galaxies have lower metallicities.

Additionally, we present a grid spanning a range of ionization parameters and metallicities. This grid demonstrates that with increasing metallicity, the grid shows an increase in the O32 value, whereas the R23 index changes nonmonotonically. Due to the nonmonotonic behavior of R23, we choose to use O32 instead of R23 for metallicity calibration. The color bar indicating metallicity values follows the calibration from F. Bian et al. (2018).

The three BPT diagrams are widely used diagnostics for the ionization properties of the ISM in galaxies. All use on their Y-axis; each uses a different X-axis: (N2), (S2), and (O1) (J. A. Baldwin et al. 1981; S. Veilleux & D. E. Osterbrock 1987).

Figures 7–9 present our N2, S2, and O1 BPT diagrams, respectively. The black points represent local galaxies from the SDSS sample, consistent with those shown in Figure 6. The cyan points denote the 43 HAEs confirmed, while the red, purple, and blue rhombic symbols correspond to the stacked spectra for high-mass, low-mass, and total stacked samples, respectively, as in Figure 6.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The left and right panels of Figures 7–9 differ in the grid lines representing two sets of models. The pink grid lines in the left panel correspond to photoionization models considering only H II regions, while the wheat and green grid lines in the right panel represent models incorporating both H II regions and shocks with different metallicity. The detailed discussion of these models is provided in Section 4.1.

In the N2 BPT diagram (Figure 7), we find that both the stacked and individual galaxy samples are offset above the local galaxy sequence. Specifically, the low-mass stacked sample exhibits a higher [O III]λ5007/Hβ ratio, whereas the high-mass stacked sample shows a lower [O III]λ5007/Hβ ratio. Specifically, the low-mass stacked sample exhibits an offset of 0.45 dex above the local galaxy sequence in the [O III]λ5007/Hβ ratio, the high-mass stacked sample exhibits 0.36 dex offset above the local galaxy sequence, and the full-stacked sample exhibits 0.37 dex offset above the local galaxy sequence. All stacked samples lie along the sequences defined by C. C. Steidel et al. (2014)—dashed yellow line—and A. E. Shapley et al. (2015)—solid green line. Thus protocluster galaxies share the same offset from local galaxies as other high-redshift galaxies at z ∼ 2.

In contrast, the S2 BPT diagram (Figure 8) reveals almost no significant offset between our samples and the local galaxies. However, in the O1 BPT diagram (Figure 9), the stacked samples lie well beyond the L. J. Kewley et al. (2001) demarcation line, which typically separates star-forming galaxies from AGNs. Notably, the high-mass and total stacked samples exhibit elevated [O I]λ6300/Hα ratios, with log([O I]λ6300/Hα) ∼ −1.45.

These results suggest that the ionization conditions in protocluster galaxies cannot be fully explained by simple stellar photoionization diagnostics alone. We therefore conduct a more detailed analysis using more complicated photoionization models.

We employ an ionization model to characterize protoclusters, focusing on the physical parameters of H II regions, including the ionization parameter (q), the slope of the ionizing spectrum of stars illuminating the gas, and the electron density (n_e). Based on the parameters calculated in Section 3, protoclusters exhibit high ionization parameters and low metallicities. The n_e for individual galaxies ranges from 34 to 2500 cm⁻³, while the n_e for the stacked samples ranges from 93 to 1067 cm⁻³. Therefore, we set the n_e of the model range from 50 to 6500 cm⁻³. The ionization parameter log(q) for most protocluster galaxies is above 7.5, and for the stacked samples, log(q) exceeds 7.75. Thus, we assume that protocluster environments have high ionization parameters, with a range for log(q) set from 6.75 to 8.50. We set the metallicity to 0.15 Z_⊙ and 0.4 Z_⊙. The plane-parallel geometry is adopted to simulate H II regions, with the parameters selected based on M. A. Dopita et al. (2013) and M. G. Allen et al. (2008), as summarized in Table 6. The detailed modeling results are presented in Table 7.

It is important to note an apparent discrepancy when comparing Figures 6 and 7. Specifically, the photoionization model grid line corresponding to a high ionization parameter (log(q) = 7.5) lies below the main locus of the SDSS galaxy sample in the O32 diagram (Figure 6). In contrast, the same grid line appears to traverse the central region of the data distribution in the N2 BPT diagram (Figure 7). However, this is not a contradiction. The apparent difference arises because the model grids overlaid on the two diagrams are not identical in their metallicity coverage. The H II region grid in Figure 6 spans a wide metallicity range (from Z = 2.0 Z_⊙ to 0.05 Z_⊙), whereas the grid in Figure 7 is shown for only two discrete metallicity values (0.4 and 0.15 Z_⊙). Upon closer inspection of Figure 7, the track corresponding to a metallicity of 0.4 Z_⊙ is indeed positioned consistently with the physical trends and data distribution shown in Figure 6.

Using the Starburst99 (C. Leitherer et al. 1999) spectrum model in MAPPINGS to fit the ionizing flux, we find that the ionizing continuum slope is 21.26.

The emission line flux ratios of [O III]/Hβ, [O I]/Hα, [N II]/Hβ, and [S II]/Hα obtained from the pure photoionization model are shown in the left panels of Figures 7–9. In the left panels, the green grid lines represent the metallicity model with 0.15 Z_⊙, while the wheat grid lines correspond to the model with 0.4 Z_⊙.

Under these model assumptions, the resulting log([O I]/Hα) ratios are all above −1.83. However, the observed values for the high-quality galaxies and the overall stacked sample are approximately −1.44. This discrepancy suggests that additional radiation sources may be enhancing the flux of the [O I] ine. Although we initially attempted to use a pure photoionization model to explain the [O I] emission line offset, we found that this approach did not adequately match the observed data. Therefore, additional radiation sources must be considered to explain the offset in the line ratios.

Shock excitation can arise from various phenomena, including supernovae, stellar winds, galaxy interactions, and AGN-related activities such as jets (L. J. Kewley et al. 2019). Shocks significantly enhance low-ionization emission lines such as [O I], [N II], and [S II]. Since these shocks can originate from different astrophysical processes, quantifying their total contribution poses considerable challenges. Therefore, this paper adopts a weighted summation approach between H II regions and shocks as a simplified example, demonstrating how combining these two components can improve the overall fit.

where ω is the contribution of shocks to the total emission, a fraction between 0 and 1. Higher ω values indicate a greater influence of shock excitation, while lower values reflect the dominance of photoionization in H II regions.

The weight ω represents the contribution of slow shocks. Specifically, the final spectrum is obtained by adding the spectrum of the shocks, scaled by the weight, to the spectrum of the H II regions. This additional free parameter allows us to simultaneously account for both the ionization driven by star formation and the additional excitation caused by shocks, providing a more comprehensive explanation for the observed emission line ratios in protocluster galaxies. The weight in the figure represents the fractional contribution of shocks to the total emission. Higher weights indicate a more significant role of shock excitation in the emission lines, while lower weights suggest that the emission is predominantly dominated by H II regions.

Shock waves can also enhance certain emission lines. Slow shocks (velocity < 500 km s⁻¹) can enhance emission lines such as [N II] and [O I], while fast shocks (velocity > 500 km s⁻¹) can also enhance the [O III] emission line (L. J. Kewley et al. 2019). Given the strong [O I] emission lines observed in Figure 9, shock excitation could be a significant contributor to the BPT diagram offsets in our sample, so in the shock model, we adopt a slow shock velocity of 200 km s⁻¹, as suggested by L. J. Kewley et al. (2019) and A. J. Cameron et al. (2020). The metallicity and magnetic field strength are chosen based on the values recommended by A. J. Cameron et al. (2020). The relevant parameters for the shock model are also summarized in Table 6. M. G. Allen et al. (2008) reported that magnetic fields influence gas compression. Specifically, in the presence of a magnetic field, the compression is proportional to the Alfvén Mach number, which can limit the degree of compression in preshock gas.

We present the mixed model of H II regions and shocks in the right panels of Figures 7–9. The green grid lines represent 0.15 Z_⊙, while the wheat grid lines represent 0.4 Z_⊙. From the figure, we can observe that comparing two figures, the pure H II region model can not fit the protocluster, while the shock plus H II region model could fit the galaxies in protoclusters.

Using the MAPPINGS V software, we find that high-redshift protoclusters are likely affected by shocks, which may account for 0%–40% of the total Hα emission line flux. This contribution helps explain the observed emission line ratio offsets in protoclusters. As shown in Figure 7, when no shock emission is present, the locus of H II regions aligns with the distribution of SDSS galaxies, with no offset. However, as the weight of slow shocks increases, log([N II]/Hα) shifts from the SDSS locus toward the fit lines reported by C. C. Steidel et al. (2014) and A. E. Shapley et al. (2015), indicating that shocks may explain the observed offsets in high-redshift galaxies. Similarly, the S2 BPT diagram shows that combining H II regions and shocks can cover most of the galaxies in our sample. For Figure 9, when there is no shock contribution, the emission line ratios are consistent with local star-forming galaxies. However, as the shock component increases, the ratios gradually move beyond the L. J. Kewley et al. (2001) demarcation line, eventually aligning well with the observed values in protoclusters.

Contrary to Z. Li et al. (2025), who found that z > 5 protocluster galaxies are more metal-rich than field galaxies of equivalent stellar mass, potentially driven by accelerated star formation, our overdense sample exhibits the opposite trend: protocluster galaxies show systematically lower metallicities compared to equally massive field galaxies. This discrepancy suggests that our protoclusters at z ∼ 2 are at different evolutionary stages, where different mechanisms may be at play.

To distinguish between photoionization from star formation and shock excitation, we constructed a kinematic diagnostic diagram plotting gas velocity dispersion (σ) against the [O I]/Hα emission line ratio.

We generated a grid of theoretical models for comparison using the MAPPINGS V code. The models were computed for two metallicities, 0.15Z_⊙ and 0.4Z_⊙, to bracket the metallicity range of our sample galaxies, which spans from to 8.6. A lower metallicity grid is therefore appropriate for our samples. The shock models were generated with parameters consistent with those listed in Table 6, with the only difference being the shock velocity. For a pure photoionization reference, we modeled an H II region with an ionization parameter of . The shock models cover a range of velocities (v_s = 100, 200, 300, 400 km s⁻¹) to match the observed velocity dispersions, which were derived from the FWHM of the Hα emission line. The uncertainties on the velocity dispersion for each galaxy were propagated from the measurement errors of the FWHM.

Figure 10 presents our observational data overlaid on the model grids. The vast majority of our sample galaxies occupy a distinct region of the diagram that cannot be simply explained by pure photoionization. Instead, their positions align remarkably well with the low-metallicity track that represents a combination of an H II region and a shock component.

Zoom In Zoom Out Reset image size

Notably, most of the data points are around the parts of the model track corresponding to low-velocity shocks, in the range of 100–200 km s⁻¹. This provides evidence that slow shocks, rather than high-velocity shocks, are a prevalent mechanism contributing to the gas excitation in this protocluster environment.

The confirmed presence of shocks has significant implications for the determination of gas-phase abundances. Emission from shock-heated gas can contaminate key diagnostic lines, such as [N II], leading to an overestimation of the intrinsic nitrogen abundance. Following the framework established by S. R. Flury et al. (2024), we can use our findings to correct for this effect.

Given that most of our galaxies are consistent with shock velocities between 100 and 200 km s⁻¹, we adopted a 200 km s⁻¹ shock model as a representative case to estimate the shock contribution to the [N II] flux. We then corrected the observed [N II]/Hα ratios that would arise purely from photoionization. Assuming a shock contribution could be as high as 40%, for our low-mass, high-mass, and total stacked samples, this would lead to a significant decrease in the [N II]/Hα ratio; this suppression can be as large as 0.48 dex. This correction brings the stacked line ratios of our sample into agreement with those of local star-forming galaxies from the SDSS. As a result, the inferred nitrogen abundances of our high-redshift protocluster galaxies become consistent with those of their low-redshift counterparts, and are also in agreement with the findings of S. R. Flury et al. (2024). It is also a possible resolution for the paper (C. C. Steidel et al. 2014; A. E. Shapley et al. 2015).

This finding suggests a potential solution to the long-standing puzzle of apparent nitrogen overabundance in high-redshift galaxies (C. C. Steidel et al. 2014; A. E. Shapley et al. 2015; S. R. Flury et al. 2024). Our analysis suggests that this observed excess may not entirely reflect intrinsic differences in chemical evolution, but may be significantly influenced by an observational artifact from unaccounted-for shock excitation. This aligns with previous studies by S. F. Newman et al. (2013) and S. R. Flury et al. (2024), which also pointed to the potential role of shock excitation mechanisms in biasing abundance measurements.

The environment plays a key role in galaxy evolution. As pioneered by A. Dressler (1980), galaxies in high-density regions have undergone a different evolutionary path than their counterparts in the field.

An important question arises: which environmental mechanisms are responsible for driving the widespread shocks we have detected? In dense protocluster environments, two primary candidates are gravitational interactions, most notably galaxy mergers, and hydrodynamic processes like ram pressure stripping (RPS).

Galaxy mergers are a well-established mechanism for generating shocks (W. J. Pearson et al. 2024). Our measured shock contribution to the Hα flux (0%–40%) and the average gas velocity dispersion (FWHM ∼ 188 km s⁻¹) are both consistent with the properties observed in local merger-induced shock systems (J. A. Rich et al. 2011), making mergers a plausible candidate.

However, to test if mergers are the dominant driver, we must consider their prevalence. Following the morphological classification catalog from S. Liu et al. (2023) for galaxies in these same protocluster fields, we cross-match our sample with catalog from S. Liu et al. (2023). Out of our 43 galaxies, we find only 7 that exhibit clear merger signatures. While mergers certainly contribute to the overall shock budget, their low observed fraction makes it unlikely that they are the primary mechanism responsible for the widespread shock signatures detected across the majority of our sample. This suggests that a more pervasive mechanism must be at play.

RPS presents another plausible origin for the shocks. The link between ram pressure and shock formation has been established by work such as H. Kang & D. Ryu (2011). Given that RPS is considered a key process in the dense environments of high-redshift protoclusters (A. Abramson & J. D. P. Kenney 2014; V. Strazzullo et al. 2019), it is a strong candidate for driving the shocks in our sample.

The detection of gas tails is a primary method for identifying the influence of RPS. Historically, RPS was first invoked to explain the morphology of head-tail radio galaxies in nearby clusters (G. K. Miley et al. 1972; G. Miley 1980). Observationally, these tails can be traced using different gas phases. Tails of neutral atomic gas (H I) are a classic indicator of RPS. Hot gas tails, visible in X-ray observations, have also been found in massive early-type galaxies. While direct evidence for RPS in molecular gas is still emerging, ionized gas has proven to be an excellent tracer. Numerous one-sided Hα tails, indicative of RPS, have been observed in nearby clusters like Virgo, Coma, and A1367, and similar evidence has been found in some protoclusters (V. Strazzullo et al. 2019). A future study with high-resolution imaging is therefore essential to search for these features and provide a conclusive verdict on the dominant environmental mechanism at play.

In conclusion, our analysis provides evidence that the widespread shocks in these protocluster galaxies are primarily driven by RPS as the most likely dominant mechanism. This conclusion is supported by two findings: (1) the low observed fraction of major mergers is insufficient to explain the prevalence of shocks, and (2) the galaxies are located in the dense cores of protoclusters, where environmental conditions favor RPS as the primary driver of the observed shocks. A definitive confirmation, however, requires high-resolution imaging to directly detect morphological signatures of these processes, such as stripped gas tails or faint tidal features.