# JWST NIRSpec High-resolution Spectroscopy of at $z=10.167$:
Resolved [OII] Doublet and Electron Density in an Early Galaxy

Abdurro’uf Department of Physics and Astronomy, The Johns Hopkins University,
3400 N Charles St. Baltimore, MD 21218, USA  Space Telescope Science Institute
(STScI), 3700 San Martin Drive, Baltimore, MD 21218, USA  Rebecca L. Larson
School of Physics and Astronomy, Rochester Institute of Technology, 84 Lomb
Memorial Drive, Rochester, NY 14623, USA  Dan Coe Space Telescope Science
Institute (STScI), 3700 San Martin Drive, Baltimore, MD 21218, USA
Association of Universities for Research in Astronomy (AURA), Inc. for the
European Space Agency (ESA)  Department of Physics and Astronomy, The Johns
Hopkins University, 3400 N Charles St. Baltimore, MD 21218, USA  Tiger Yu-Yang
Hsiao Department of Physics and Astronomy, The Johns Hopkins University, 3400
N Charles St. Baltimore, MD 21218, USA  Space Telescope Science Institute
(STScI), 3700 San Martin Drive, Baltimore, MD 21218, USA  Javier Álvarez-
Márquez Centro de Astrobiología (CAB), CSIC-INTA, Ctra. de Ajalvir km 4,
Torrejón de Ardoz, E-28850, Madrid, Spain  Alejandro Crespo Gómez Centro de
Astrobiología (CAB), CSIC-INTA, Ctra. de Ajalvir km 4, Torrejón de Ardoz,
E-28850, Madrid, Spain  Angela Adamo Department of Astronomy, Stockholm
University, Oscar Klein Centre, AlbaNova University Centre, 106 91 Stockholm,
Sweden  Rachana Bhatawdekar European Space Agency (ESA), European Space
Astronomy Centre (ESAC), Camino Bajo del Castillo s/n, 28692 Villanueva de la
Cañada, Madrid, Spain  Arjan Bik Department of Astronomy, Stockholm
University, Oscar Klein Centre, AlbaNova University Centre, 106 91 Stockholm,
Sweden  Larry D. Bradley Space Telescope Science Institute (STScI), 3700 San
Martin Drive, Baltimore, MD 21218, USA  Christopher J. Conselice Jodrell Bank
Centre for Astrophysics, University of Manchester, Oxford Road, Manchester UK
Pratika Dayal Kapteyn Astronomical Institute, University of Groningen, 9700 AV
Groningen, The Netherlands  Jose M. Diego Instituto de Física de Cantabria
(CSIC-UC). Avda. Los Castros s/n. 39005 Santander, Spain  Seiji Fujimoto
Hubble Fellow  The University of Texas at Austin, Department of Astronomy,
Austin, TX, United States  Lukas J. Furtak Physics Department, Ben-Gurion
University of the Negev, P.O. Box 653, Be’er-Sheva 84105, Israel  Taylor A.
Hutchison NASA Postdoctoral Fellow  Observational Cosmology Lab, NASA Goddard
Space Flight Center, Greenbelt, MD 20771, USA  Intae Jung Space Telescope
Science Institute, 3700 San Martin Drive Baltimore, MD 21218, United States
Meghana Killi Instituto de Estudios Astrofísicos, Facultad de Ingeniería y
Ciencias, Universidad Diego Portales, Av. Ejército 441, Santiago 8370191,
Chile  Vasily Kokorev Kapteyn Astronomical Institute, University of Groningen,
P.O. Box 800, 9700 AV Groningen, The Netherlands  Matilde Mingozzi Space
Telescope Science Institute (STScI), 3700 San Martin Drive, Baltimore, MD
21218, USA  Colin Norman Department of Physics and Astronomy, The Johns
Hopkins University, 3400 N Charles St. Baltimore, MD 21218, USA  Space
Telescope Science Institute (STScI), 3700 San Martin Drive, Baltimore, MD
21218, USA  Tom Resseguier Department of Physics and Astronomy, The Johns
Hopkins University, 3400 N Charles St. Baltimore, MD 21218, USA  Space
Telescope Science Institute (STScI), 3700 San Martin Drive, Baltimore, MD
21218, USA  Massimo Ricotti Department of Astronomy, University of Maryland,
College Park, 20742, USA  Jane R. Rigby Observational Cosmology Lab, NASA
Goddard Space Flight Center, Greenbelt, MD 20771, USA  Eros Vanzella INAF –
OAS, Osservatorio di Astrofisica e Scienza dello Spazio di Bologna, via
Gobetti 93/3, I-40129 Bologna, Italy  Brian Welch Observational Cosmology Lab,
NASA Goddard Space Flight Center, Greenbelt, MD 20771, USA  Department of
Astronomy, University of Maryland, College Park, 20742, USA  Rogier A.
Windhorst School of Earth and Space Exploration, Arizona State University,
Tempe, AZ 85287-1404, USA  Xinfeng Xu Department of Physics and Astronomy,
Northwestern University, 2145 Sheridan Road, Evanston, IL, 60208, USA  Adi
Zitrin Physics Department, Ben-Gurion University of the Negev, P.O. Box 653,
Be’er-Sheva 84105, Israel

###### Abstract

We present JWST/NIRSpec high-resolution spectroscopy G395H/F290LP of , a
gravitationally lensed galaxy merger at $z=10.167$. The new spectroscopy,
which is acquired for the two lensed images (JD1 and JD2), detects and
resolves emission lines in the rest-frame ultraviolet (UV) and blue optical,
including the resolved [O ii] $\lambda\lambda$3726,3729 doublet, [Ne iii]
$\lambda$3870, HeI $\lambda$3890, H$\delta$, H$\gamma$, and [O iii]
$\lambda$4363\. This is the first observation of the resolved [O ii]
$\lambda\lambda$3726,3729 doublet for a galaxy at $z>8$. We measure a line
flux ratio [O ii] $\lambda$3729/$\lambda$3726 $=0.9\pm 0.3$, which corresponds
to an estimated electron density of $\log(n_{e}/\rm{cm}^{-3})=2.9\pm 0.5$.
This is significantly higher than the electron densities of local galaxies
reported in the literature. We compile the measurements from the literature
and further analyze the redshift evolution of $n_{e}$. We find that the
redshift evolution follows the power-law form of $n_{e}=A\times(1+z)^{p}$ with
$A=54^{+31}_{-23}$ cm-3 and $p=1.2^{+0.4}_{-0.4}$. This power-law form may be
explained by a combination of metallicity and morphological evolution of
galaxies, which become, on average, more metal-poor and more compact with
increasing redshift.

Galaxies (573), High-redshift galaxies (734), Early universe (435), Strong
gravitational lensing (1643), Galaxy spectroscopy (2171)

††journal: ApJ††facilities: JWST/NIRSpec††software:  grizli (Brammer et al.,
2022) astropy (The Astropy Collaboration et al., 2018), matplotlib (Hunter,
2007), NumPy (Harris et al., 2020), SciPy (Virtanen et al., 2020), STScI JWST
Calibration Pipeline (jwst-pipeline.readthedocs.io (Rigby et al., 2023) IDL
Astronomy Library: idlastro.gsfc.nasa.gov (Landsman, 1993), piXedfit
(Abdurro’uf et al., 2021, 2022), msaexp, pyneb (Luridiana et al., 2015)

##  1 Introduction

Understanding the physical conditions in star-forming galaxies, including the
interstellar medium (ISM), is essential for a complete picture of the
evolution of the stellar and gaseous content of galaxies. Spectroscopy of
galaxies in the rest-frame ultraviolet (UV) and optical provides a wealth of
emission line diagnostics for characterizing the ISM in galaxies (e.g., see
review by Kewley et al., 2019b). The relative strengths of the emission lines
are mainly driven by the ISM properties, including the chemical abundance,
shape of the ionizing radiation field, ionization state, gas density, and
temperature (e.g., Kewley & Dopita, 2002; Dopita et al., 2006a, b).

Before the launch of the JWST, the observations of rest-frame UV and optical
spectra of high-redshift ($z\gtrsim 3$) galaxies were very challenging because
of the lack of highly sensitive (i.e., deep) near-infrared (NIR) spectroscopy.
Only a few spectroscopic surveys have been conducted to observe rest-frame UV
and optical spectroscopy of high-$z$ galaxies out to $z\sim 4$ using 8–10 m
class telescopes (e.g., Kriek et al., 2015; Le Fèvre et al., 2015; Bacon et
al., 2017; Pentericci et al., 2018; Urrutia et al., 2019).

Since the launch of JWST (Rigby et al., 2023; Gardner et al., 2023), with its
NIR and mid-infrared (MIR) spectroscopic capabilities, the rest-frame UV and
optical emission lines of $z>2$ galaxies have become easily accessible,
enabling the study of the ISM properties in galaxies at high redshift (toward
the reionization era). Several studies use JWST/Near Infrared Spectrograph
(NIRSpec) spectroscopy (Jakobsen et al., 2022; Böker et al., 2023) to study
the properties of ISM in galaxies up to $z\sim 9$ (e.g., Fujimoto et al.,
2022; Curti et al., 2023; Fujimoto et al., 2023; Tang et al., 2023; Bunker et
al., 2023; Heintz et al., 2023a, b; Reddy et al., 2023; Nakajima et al., 2023;
Hsiao et al., 2023a; Isobe et al., 2023; Jung et al., 2023; Cameron et al.,
2023; Sanders et al., 2024; Backhaus et al., 2024). Those studies found that
high-$z$ galaxies have significantly lower gas-phase metallicity (e.g.,
Fujimoto et al., 2023; Heintz et al., 2023a; Nakajima et al., 2023; Williams
et al., 2023) and higher electron density ($n_{e}$; e.g., Reddy et al. 2023;
Isobe et al. 2023) than low-$z$ galaxies.

The electron density of the H ii regions is one of the important properties of
the ISM because it is an excellent tracer of the gas density in the star-
forming regions. This quantity can be estimated using density-sensitive
emission line ratios, including [O ii] $\lambda$3729/$\lambda$3726, [S ii]
$\lambda$6716/$\lambda$6731, C iii] $\lambda$1909/$\lambda$1907 (e.g., Kewley
et al., 2019b). In star-forming regions, the excitation energy between the two
lines of those doublets is proportional to the thermal electron density.
Therefore, the relative excitation rates depend only on the collision
strengths of the electrons (Osterbrock, 1989). This makes the flux ratio
between the two lines a perfect diagnostic for the electron density.

Before JWST, several large spectroscopic surveys investigated the evolution of
the electron density based on [O ii] $\lambda\lambda$3726,3729 and [S ii]
$\lambda\lambda$6716,6731 (i.e., from low-ionization regions). They found that
$n_{e}$ declines from $\sim 200$ cm-3 at $z\sim 3$ to $\sim 30$ cm-3 at $z\sim
0$ (e.g., Steidel et al., 2014; Sanders et al., 2016; Kaasinen et al., 2017;
Kashino et al., 2017; Harshan et al., 2020; Davies et al., 2021; Berg et al.,
2022). Due to the lack of highly sensitive NIR spectroscopy, there was a lack
of $n_{e}$ studies at $z\gtrsim 3$. A few exceptions come from studies that
estimated $n_{e}$ based on doubly-ionized regions, such as Killi et al. (2023)
who observed a galaxy at $z=7.133$ with Atacama Large Millimeter Array (ALMA)
and estimated $n_{e}\lesssim 500$ cm-3 based on the [O iii]
$\lambda$52$\mu$m/[O iii] $\lambda$88$\mu$m ratio. However, this diagnostic
probes a very different ISM region (i.e., dust-obscured) and ionization level
compared to the [O ii] $\lambda\lambda$3726,3729 and [S ii]
$\lambda\lambda$6716,6731 doublets.

After the launch of JWST, several studies have already explored $n_{e}$ out to
$z\sim 9$ using various diagnostics. Isobe et al. (2023) used the
spectroscopic data from public JWST surveys to measure $n_{e}$ of 14 galaxies
at $z=4.02-8.68$. Among the sample, only three galaxies (with the highest
$z=7.87$) observed with NIRSpec high-resolution grating by the GLASS survey
(Treu et al., 2022) showed resolved [O ii] $\lambda\lambda$3726,3729 doublet,
while the doublet remains unresolved in the rest of the sample. They estimated
$n_{e}\gtrsim 300$ cm-3. Reddy et al. (2023) investigated the connections
between the ionization parameter ($U$), electron density ($n_{e}$), and star
formation rate (SFR) based on 48 galaxies at $2.7<z<6.3$. They measured
$n_{e}$ based on the [S ii] $\lambda\lambda$6716,6731 doublet resolved in
NIRSpec medium-resolution spectra from the Cosmic Evolution Early Release
Science (CEERS; Finkelstein et al., 2023; Fujimoto et al., 2023). They found
that galaxies with higher O32 = [O iii] $\lambda$4959,$\lambda$5007/[O ii]
$\lambda\lambda$3726,3729 have average $n_{e}\simeq 500$ cm-3 that are
$\gtrsim 5$ times larger than that of lower-O32 galaxies. They also found a
highly significant positive correlation between $U$ and SFR surface density
($\Sigma_{\rm SFR}$), which appears to be independent of redshift at
$1.6\lesssim z\lesssim 6.3$.

While JWST/NIRSpec prism spectroscopy cannot resolve [O ii]
$\lambda\lambda$3726,3729 and [S ii] $\lambda\lambda$6716,6731 lines due to
low spectral resolution ($R\sim 30-300$), some studies have combined [O iii]
$\lambda$5007 observed by this instrument with [O iii] $\lambda$88$\mu$m
observed with ALMA to indirectly estimate $n_{e}$ of galaxies at $z=9.11$
($\sim 400$ cm-3; Stiavelli et al., 2023) and $z=8.496$ ($\sim 220$ cm-3;
Fujimoto et al., 2022) using photoionization modeling. Until now, no studies
have directly measured $n_{e}$ from resolved [O ii] $\lambda\lambda$3726,3729
doublet at $z>8$, while [S ii] $\lambda\lambda$6716,6731 is redshifted beyond
the wavelength coverage of NIRSpec. The detection of [O ii]
$\lambda\lambda$3726,3729 doublet remains challenging because most galaxies at
such high redshift are faint. Moreover, high resolution ($R>2000$)
spectroscopy is needed to be able to resolve the doublets, which is only
achievable using the NIRSpec high-resolution grating. Recently, the [O ii]
$\lambda\lambda$3726,3729 doublet has been observed in the spectrum of GN-z11
by Bunker et al. (2023), one of the brightest galaxies known at $z>10$ (Oesch
et al., 2016; Tacchella et al., 2023). However, the NIRSpec medium-resolution
grating cannot resolve the single lines.

In this paper, we present NIRSpec high-resolution ($R\sim 2700$) spectroscopy
of , a triply-lensed galaxy at $z=10.167$ (Coe et al., 2013; Hsiao et al.,
2023b, a). Previously, JWST/NIRCam imaging and NIRSpec prism spectroscopy have
been obtained for this galaxy by GO 1433 (PI Coe). The NIRCam imaging resolved
the galaxy into two components A and B, suggesting a possible galaxy merger
(Hsiao et al., 2023b). NIRSpec prism spectra of JD1 and JD2 (the two brightest
lensed images of ) were presented in Hsiao et al. (2023a), revealing 7 strong
emission lines and becoming one of the highest redshift galaxies with detected
rest-frame UV emission lines besides GHZ2/GLASS-z12 (Castellano et al., 2024),
GN-z11 (Bunker et al., 2023), and JADES-GS-z11-0 (Hainline et al., 2024).
Here, we observe JD1 and JD2 with the high-resolution G395H/F290LP setting to
further resolve the lines and obtain better measurements of their fluxes and
widths. We resolve [O ii] $\lambda\lambda$3726,3729 doublet for the first time
at $z>8$ and use it to estimate $n_{e}$.

The paper is organized as follows. We describe the NIRSpec data and its
reduction in Section 2. The measurements of emission lines are described in
Section 3. We present and discuss our results in Section 4, and finally
conclude in Section 5. Throughout, we assume the cosmological parameters of
$\Omega_{m}=0.3$, $\Omega_{\Lambda}=0.7$, and $H_{0}=70\text{ km}\text{
s}^{-1}\text{ Mpc}^{-1}$.

##  2 Data

###  2.1 NIRSpec MSA Spectroscopy

The NIRSpec micro-shutter assembly (MSA) Spectroscopy Ferruit et al. (2022)
observation was conducted on January 14, 2024, as part of GO 4246 (PI
Abdurro’uf), which focuses on at $z=10.17$. We use the high-resolution grating
G395H/F290LP that covers a wavelength range of $2.87–5.14$ $\mu$m with
spectral resolution ranges from $R\sim 1900$ to $3500$ over the wavelength
range. The total exposure time was 2.43 hours. Half of that time was observed
with standard 3-slitlet nods. With the other half, we used single slitlets and
performed two dithers. The slitlets on JD1 and JD2 mainly cover the brightest
Component A (see Figure 1). This is similar to the NIRSpec prism observation
previously performed (Hsiao et al., 2023a). The data are all publicly
available on MAST, but not the reduction we analyzed in this paper which will
be described below.111https://mast.stsci.edu/search/ui/#/jwst.

###  2.2 NIRSpec Data Reduction

Figure 1: NIRSpec G395H/F290LP spectra of lensed images JD1 and JD2 (IDs 3593
and 3349, respectively, in the catalog used for preparing the observation).
Left panel: Slitlet positions overlaid on $1.6^{\prime\prime}\times
1.6^{\prime\prime}$ NIRCam color images (blue:F115W; green:F150W; red:F200W).
Right panel: 2D and 1D NIRSpec high-resolution spectra of JD1 (top) and JD2
(bottom). Several emission lines are detected while the continuum is not.
Emission lines are indicated with blue-colored labels and orange vertical
dashed lines. The gray shaded area represents the uncertainties obtained from
the pipeline. The green lines represent the standard deviations of the
spectral fluxes after removing the emission lines. We multiply the flux
uncertainties by a factor of two to match the overall level of the standard
deviations. This is shown with a red line. These rescaled uncertainties are
used for the analysis in this paper.

The data were reduced with the JWST pipeline using CRDS version 11.17.6 and
jwst_1183.pmap with the 2D spectra obtained from the Mikulski Archive for
Space Telescopes (MAST). The 2D frames were visually inspected for artifacts
and the snowballs and residual cosmic rays were manually masked. Additional
sigma clipping with a threshold of 6$\sigma$ was then applied to the 2D
spectra to remove any additional spurious pixels. The source’s one-dimensional
(1D) spectra were obtained via an optimal extraction (Horne, 1986) with a
spatial weight profile from the trace of the source in the masked 2D spectrum
such that the pixels near the peak of the trace are maximally weighted. To
create the extraction profile, the 2D signal-to-noise (SNR) spectrum was
collapsed in the spectral direction, taking the median value at each spectral
pixel and fitting a Gaussian to the positive trace. While the NIRSpec shutter
includes flux from both Component A and Component B for JD1 (ID 3593), the
orientation is such that Component B falls completely out of the shutter for
JD2 (ID 3349). The source is effectively spatially unresolved at all
wavelengths (intrinsic size is less than the JWST point spread function full
width at half maximum, PSF FWHM), and there was only a single detectable trace
in both observations at the expected location of Component A. An attempt to
extract a spectrum for both Component A and Component B from JD1 yielded no
discernible detections for Component B, which was affected by the location of
the upper negative trace in the 2D spectrum.

Figure 1 shows 2D and 1D spectra of JD1 (top right) and JD2 (bottom right)
along with the slitlet positions overlaid on $1.6^{\prime\prime}\times
1.6^{\prime\prime}$ NIRCam color images (left). We detect several emission
lines, including [O ii] $\lambda\lambda$3726,3729 doublet, [Ne iii]
$\lambda$3870, HeI $\lambda$3890, H$\delta$, H$\gamma$, and [O iii]
$\lambda$4363\. However, we do not detect the spectral continuum. The
positions of emission lines are indicated with orange dashed lines and blue
labels. The gray shaded area represents the flux uncertainties produced by the
pipeline. These uncertainties might be underestimated (see e.g., Ferruit et
al., 2022; Christensen et al., 2023). Therefore, we check it by calculating
the standard deviations of the spectral fluxes after removing the emission
lines ($\pm 50$ Å). The green lines in Figure 1 show the standard deviation
level, which is higher by a factor of $\sim 2$ from the estimated flux
uncertainties from the pipeline. The red lines represent the flux
uncertainties multiplied by a factor of two, which we use for our analysis in
this paper. The reprocessed spectra will be publicly available on our Cosmic
Spring website222https://cosmic-spring.github.io.

##  3 Analysis

Most of the emission lines that were previously detected with the NIRSpec
prism data (see Hsiao et al., 2023a) are also detected with the current
NIRSpec high-resolution data, except [C iii] $\lambda$1907 which is beyond the
wavelength coverage of the current data. With the high-resolution grating, we
can resolve [O ii] $\lambda\lambda$3726,3729 doublet for the first time at
$z>8$. At the central wavelength of [O ii] $\lambda\lambda$3726,3729, the
G395H/F290LP grating has $R\sim 2855$ (Jakobsen et al., 2022) which
corresponds to a rest-frame dispersion of $\sim 1.3$ Å per spectral resolution
element. Besides this, HeI $\lambda$3890 is separated from [Ne iii]
$\lambda$3870, which were blended in the NIRSpec prism spectra. We describe
the emission line analysis in this section.

###  3.1 Stacking and Emission Line Measurements

We stack (i.e., sum) the spectra of JD1 and JD2 to get a good overall S/N
ratio. The stacked spectrum has a total lensing magnification ($\mu$) of
$13.3$ ($8.0+5.3$; see Hsiao et al. 2023a, b). We do not de-lense the spectra
of JD1 and JD2 before the stacking because the true line ratios of are
expected to be preserved. We measure the emission line flux and width in the
individual spectra (JD1 and JD2) and the stacked spectrum. We fit the emission
lines individually with a Gaussian function using piXedfit333The line fitting
function will be included in the future release of piXedfit (Abdurro’uf et
al., 2021, 2022, 2023). The code uses the Markov Chain Monte Carlo (MCMC)
method through the emcee package (Foreman-Mackey et al., 2013) for sampling
the posterior probability distributions of free parameters in the fitting.

First, we determine the systemic redshift by fitting the strongest line [Ne
iii] $\lambda$3870 as it also has the highest S/N ratio. The spectroscopic
redshift of has been previously determined by Hsiao et al. (2023a) using
NIRSpec prism data to be $10.17$. Here, we put a stronger constraint on the
redshift using NIRSpec high-resolution data. We obtain
$z=10.1672_{-0.0003}^{+0.0003}$, $10.1685_{-0.0004}^{+0.0004}$, and
$10.1674^{+0.0002}_{-0.0002}$ for JD1, JD2, and stacked spectrum,
respectively. For each spectrum, we use the systemic redshift to estimate the
observed central wavelength of the other emission lines and set the initial
positions of the MCMC walkers. To better constrain the FWHM of the lines, we
use the FWHM posterior distribution function derived from the fitting of [Ne
iii] $\lambda$3870 as an additional prior. For [O ii]
$\lambda\lambda$3726,3729 doublet, we set the two lines to have the same FWHM,
while the value is allowed to vary during the fitting process.

For each emission line, we fit the spectrum region within $\pm 100$ Å around
the expected peak wavelength to minimize the influence of noise. We estimate
the continuum level by taking the median of fluxes within the fitting region
after excluding the emission line (i.e., masked around $\pm 40$ Å from the
expected peak wavelength). We then subtract the continuum from the spectrum
and fit the emission line. We scale up the flux uncertainties by a factor of
two to match the overall level of the standard deviation (see Figure 1).

We do not use the measured continuum fluxes from the NIRSpec prims spectra
because of the differences in the slitlet position and coverage. In Hsiao et
al. (2023a), the measurements of emission lines were performed to a stacked
spectrum, which is summed of the JD1 and JD2 spectra from two observations
(Obs 21 and Obs 23). In Appendix C, we perform emission line measurements of
the individual NIRSpec prism spectra and the stacked one using the updated
reductions and fitting method.

###  3.2 Measured Emission Line Properties

Figure 2: Top: Fitting results of the emission lines in the stacked spectrum.
The gray shaded area shows the original flux uncertainties, while the red line
represents a factor of two scale-up, which is used in the fitting. This
matches with the overall level of standard deviations as shown by the green
line. The detected emission lines are indicated with orange vertical dashed
lines and blue labels. Middle and bottom: cutouts showing the detected
emission lines and the best-fit Gaussian model drawn from the posteriors
obtained from the MCMC fitting. The best-fit models are shown by the red lines
(median) and red-shaded areas (range between 16th and 84th percentiles). The
best fits of the individual [O ii] $\lambda\lambda$3726,3729 lines are shown
in blue and green, while the total is shown in red. The flux uncertainties in
the cutouts have been scaled up by a factor of two.

We show the fitting results for the stacked spectrum, JD1, and JD2 in Figures
2, 3, and 4, respectively. Measured properties of the emission lines are
summarized in Tables 1, 2, and 3. In each figure, the full spectrum is shown
in the first row, while cutouts of emission lines are shown in the second and
third rows. The red line and shaded area in each cutout show the Gaussian
model drawn from the median posteriors and the range given by the 16th and
84th percentiles, respectively. For the [O ii] $\lambda\lambda$3726,3729
doublet, the best fits of the individual lines are shown by the blue and green
dashed lines, while the red line shows the total. The flux uncertainties shown
in the cutout spectra have been scaled up by a factor of two.

Figure 3: Same as Figure 2 but for JD1 spectrum.

We detect and resolve [O ii] $\lambda\lambda$3726,3729 doublet with S/N of
$\sim 3-4$ on the individual lines. The flux ratio between the two lines,
$R_{[\rm{OII}]}$$=$[O ii] $\lambda$3729/$\lambda$3726 are
$0.71^{+0.32}_{-0.22}$, $1.04^{+0.27}_{-0.30}$, and $0.91^{+0.29}_{-0.27}$ for
the JD1, JD2, and stacked spectrum, respectively. The flux ratios are
consistent within their uncertainties, though the median values are slightly
different. This is likely caused by the noise and the different slit coverage
on JD1 and JD2 (see Figure 1). The stacked [O ii] $\lambda\lambda$3726,3729
has an observed FWHM of $183_{-15}^{+15}$ km s-1, or $150_{-19}^{+18}$ km s-1
intrinsic after correcting for instrumental line broadening of $105$ km s-1.

Figure 4: Same as Figure 2 but for JD2 spectrum. The fitting results of [O
iii] $\lambda$5007 and HeI $\lambda$4473 are omitted because of the
insignificant detection of the lines.

The strongest line [Ne iii] $\lambda$3870 is detected with a high S/N of
$11.2$ and $6$ in JD1 and JD2, while it is $16.6$ in the stacked spectrum. The
observed systemic FWHM from the stacked [Ne iii] $\lambda$3870 is
$176_{-14}^{+16}$ km s-1, or $144_{-17}^{+19}$ km s-1 intrinsic after
correcting for instrumental broadening of $101$ km s-1. The good detection of
[Ne iii] $\lambda$3870 and [O ii] $\lambda\lambda$3726,3729 allows us to
measure Ne3O2$\equiv$[Ne iii] $\lambda$3870/[O ii] $\lambda\lambda$3726,3729
ratio. From the NIRSpec high-resolution data, we obtain Ne3O2 of
$2.1_{-0.4}^{+0.7}$, $0.7_{-0.2}^{+0.2}$, and $1.2_{-0.1}^{+0.2}$ for JD1,
JD2, and stacked spectrum, respectively. We also determine Ne3O2 from NIRSpec
prism data (see Appendix C) and obtain $1.2_{-0.4}^{+0.6}$,
$1.0_{-0.2}^{+0.3}$, $1.1_{-0.1}^{+0.2}$, $1.7_{-0.3}^{+0.4}$ for the JD1 Obs
21, JD1 Obs23, JD2 Obs21, JD2 Obs23, respectively. Again, the discrepancies
can be caused by the difference in the slit coverage and noise (in the case of
the high-resolution data). Recent studies of high redshift ($3\lesssim
z\lesssim 9$) galaxies using JWST/NIRSpec data measured Ne3O2 in the range of
$\sim 0.2-3$ (e.g., Tang et al., 2023; Reddy et al., 2023; Cameron et al.,
2023; Sanders et al., 2024; Backhaus et al., 2024), in agreement with the
values determined for .

Table 1: Measured Emission Line Properties of the stacked Spectrum

Emission Line | Rest wavelength | Observed wavelength | Observed Fluxa | S/N | Observed FWHMb
---|---|---|---|---|---
| ($\rm\AA$) | ($\rm\mu m$) | ($10^{-20}$ erg s-1 cm-2) |  | (km s-1)
[O ii] | 3727.10 | $4.16165_{-0.00039}^{+0.00035}$ | $55_{-14}^{+14}$ | $3.9$ | $183_{-15}^{+15}$
[O ii] | 3729.86 | $4.16473_{-0.00039}^{+0.00036}$ | $53_{-14}^{+14}$ | $3.9$ | $183_{-15}^{+15}$
[Ne iii] | 3869.86 | $4.32164_{-0.00007}^{+0.00007}$ | $157_{-9}^{+10}$ | $16.6$ | $176_{-14}^{+16}$
He i | 3889.75 | $4.34396_{-0.00030}^{+0.00028}$ | $45_{-9}^{+9}$ | $5.0$ | $173_{-15}^{+15}$
[Ne iii] | 3968.59 | $4.43192_{-0.00050}^{+0.00025}$ | $29_{-29}^{+21}$ | $1.4$ | $177_{-17}^{+37}$
H$\delta$ | 4102.89 | $4.58172_{-0.00126}^{+0.00061}$ | $45_{-35}^{+23}$ | $1.8$ | $162_{-11}^{+15}$
H$\gamma$ | 4341.69 | $4.84842_{-0.00027}^{+0.00025}$ | $113_{-20}^{+20}$ | $5.7$ | $162_{-13}^{+13}$
[O iii] | 4364.44 | $4.87352_{-0.00113}^{+0.00085}$ | $58_{-19}^{+18}$ | $3.1$ | $159_{-14}^{+13}$
He i | 4472.74 | $4.99409_{-0.00560}^{+0.00151}$ | $24_{-18}^{+23}$ | $1.3$ | $151_{-15}^{+13}$

aafootnotetext: Measured fluxes from the observed spectrum that have not been
corrected for a lensing magnification. The stacked spectrum has the total
magnification of $\mu=13.3$ ($8.0+5.3$).

bbfootnotetext: Observed FWHM before correcting for the instrumental
broadening.

Table 2: Measured Emission Line Properties of the JD1 Spectrum

Emission Line | Rest wavelength | Observed wavelength | Observed Fluxa | S/N | Observed FWHMb
---|---|---|---|---|---
| ($\rm\AA$) | ($\rm\mu m$) | ($10^{-20}$ erg s-1 cm-2) |  | (km s-1)
[O ii] | 3727.10 | $4.16086_{-0.00034}^{+0.00030}$ | $34_{-10}^{+9}$ | $3.3$ | $175_{-22}^{+24}$
[O ii] | 3729.86 | $4.16517_{-0.00034}^{+0.00030}$ | $23_{-9}^{+10}$ | $2.4$ | $175_{-22}^{+24}$
[Ne iii] | 3869.86 | $4.32154_{-0.00011}^{+0.00010}$ | $121_{-11}^{+11}$ | $11.2$ | $180_{-22}^{+25}$
He i | 3889.75 | $4.34390_{-0.00051}^{+0.00046}$ | $41_{-11}^{+11}$ | $3.8$ | $184_{-22}^{+23}$
[Ne iii] | 3968.59 | $4.43389_{-0.00045}^{+0.00048}$ | $39_{-9}^{+9}$ | $3.4$ | $188_{-23}^{+25}$
H$\delta$ | 4102.89 | $4.58160_{-0.00038}^{+0.00073}$ | $34_{-17}^{+11}$ | $2.2$ | $162_{-15}^{+20}$
H$\gamma$ | 4341.69 | $4.84823_{-0.00018}^{+0.00015}$ | $63_{-9}^{+9}$ | $7.0$ | $152_{-18}^{+18}$
[O iii] | 4364.44 | $4.87314_{-0.00070}^{+0.00095}$ | $37_{-11}^{+11}$ | $3.4$ | $167_{-21}^{+20}$
He i | 4472.74 | $4.99419_{-0.00138}^{+0.00229}$ | $26_{-19}^{+20}$ | $1.5$ | $152_{-20}^{+22}$

aafootnotetext: Measured fluxes from the observed spectrum that have not been
corrected for a lensing magnification, which is $\mu=8$ for JD1.

bbfootnotetext: Observed FWHM before correcting for the instrumental
broadening.

Table 3: Measured Emission Line Properties of the JD2 Spectrum

Emission Line | Rest wavelength | Observed wavelength | Observed Fluxa | S/N | Observed FWHMb
---|---|---|---|---|---
| ($\rm\AA$) | ($\rm\mu m$) | ($10^{-20}$ erg s-1 cm-2) |  | (km s-1)
[O ii] | 3727.10 | $4.16192_{-0.00035}^{+0.00032}$ | $22_{-8}^{+8}$ | $2.8$ | $156_{-31}^{+37}$
[O ii] | 3729.86 | $4.16500_{-0.00035}^{+0.00032}$ | $29_{-8}^{+8}$ | $3.4$ | $156_{-31}^{+37}$
[Ne iii] | 3869.86 | $4.32203_{-0.00017}^{+0.00016}$ | $37_{-6}^{+6}$ | $6.0$ | $143_{-25}^{+33}$
He i | 3889.75 | $4.34394_{-0.00327}^{+0.00172}$ | $7_{-5}^{+7}$ | $1.2$ | $130_{-50}^{+47}$
[Ne iii] | 3968.59 | $4.43354_{-0.00246}^{+0.00064}$ | $11_{-6}^{+5}$ | $1.9$ | $149_{-46}^{+44}$
H$\delta$ | 4102.89 | $4.58226_{-0.00091}^{+0.00074}$ | $12_{-9}^{+8}$ | $1.6$ | $125_{-32}^{+42}$
H$\gamma$ | 4341.69 | $4.84911_{-0.00057}^{+0.00082}$ | $30_{-12}^{+12}$ | $2.5$ | $144_{-43}^{+39}$
[O iii] | 4364.44 | $4.87172_{-0.00149}^{+0.00278}$ | $6_{-4}^{+4}$ | $1.4$ | $102_{-39}^{+44}$
He i | 4472.74 | $4.34394_{-0.00327}^{+0.00172}$ | $7_{-5}^{+7}$ | $1.2$ | $130_{-50}^{+47}$

aafootnotetext: Measured fluxes from the observed spectrum that have not been
corrected for a lensing magnification, which is $\mu=5.3$ for JD2.

bbfootnotetext: Observed FWHM before correcting for the instrumental
broadening.

##  4 Result and Discussion

###  4.1 Measurement of Electron Density

Figure 5: The relationship between [O ii] $\lambda$3729/$\lambda$3726 and the
electron density drawn from the photoionization models generated by Kewley et
al. (2019a). The curves for different temperatures are shown with different
colors. Equations A1 and A2 give the formulas to reconstruct these curves. The
black dashed line represents the model from Sanders et al. (2016), which
perfectly matches the Kewley et al. (2019a) model for $T_{e}=10000$ K. The
vertical gray dashed lines mark the range of electron densities for which the
line ratio is a useful density diagnostic according to Kewley et al. (2019a).
Given the [O ii] $\lambda$3729/$\lambda$3726 line ratio and measured
$T_{e}=17000$ K (Hsiao et al., in prep.), the estimated $n_{e}$ of JD1, JD2,
and the stacked spectrum are shown by the red, blue, and black circles,
respectively.

We measure $n_{e}$ from the line ratio $R_{[\rm{OII}]}$ based on the
photoionization models adopted from Kewley et al. (2019a), which were
generated using the MAPPINGS version 5.1 photoionization code (Binette et al.,
1985; Sutherland & Dopita, 1993; Dopita et al., 2015). We perform
interpolations to the model grids from Kewley et al. (2019a) and construct a
functional formula for converting $R_{[\rm{OII}]}$ to $n_{e}$ that takes into
account the electron temperature ($T_{e}$). We describe this in detail in
Appendix A.

Figure 5 shows the model $R_{\rm[OII]}(n_{e})$ curves for 9 values of $T_{e}$
calculated using the constructed formula (Equations A1 and A2). Different
electron temperatures are indicated with different colors. The vertical gray
dashed lines mark the range of electron densities for which the line ratio is
a useful density diagnostic according to Kewley et al. (2019a). The black
dashed line shows the relation from Sanders et al. 2016 (Eq. 6 therein). It
perfectly matches the curve from Kewley et al. (2019a) for $T_{e}=10000$ K,
which was indeed assumed in Sanders et al. (2016).

Figure 6: Resolved [O ii] $\lambda\lambda$3726,3729 doublet of JD1 (left), JD2
(middle), and the stacked spectrum (right). The two lines are resolved in the
three spectra. Green dashed and blue dashed lines are the best-fit Gaussian of
[O ii] $\lambda$3726 and [O ii] $\lambda$3729 respectively, while the solid
red line represents the total. Derived $R_{[\rm{OII}]}$ of JD1, JD2, and
stacked spectrum are $0.71^{+0.32}_{-0.22}$, $1.04^{+0.27}_{-0.30}$, and
$0.91^{+0.29}_{-0.27}$, respectively. Assuming $T_{e}=17000$ K (Hsiao et al.,
in prep.), we estimate $n_{e}$ of $\log(n_{e})=3.3^{+0.6}_{-0.5}$,
$2.7^{+0.5}_{-0.6}$, and $2.9^{+0.5}_{-0.5}$, respectively.

We use the new formula for calculating the $n_{e}$ of system. From our fitting
(see Sections 3.1 and 3.2), we obtain $R_{[\rm{OII}]}$ of
$0.71^{+0.32}_{-0.22}$, $1.04^{+0.27}_{-0.30}$, and $0.91^{+0.29}_{-0.27}$ for
the JD1, JD2, and stacked spectrum, respectively. Assuming $T_{e}=17000$ K for
as calculated in a companion paper (Hsiao et al., in prep.), we obtain
electron densities of $\log(n_{e}/\rm{cm}^{-3})=3.3^{+0.6}_{-0.5}$,
$2.7^{+0.5}_{-0.6}$, and $2.9^{+0.5}_{-0.5}$, respectively. In Figure 5, those
values are indicated with blue, red, and black circles, respectively. With the
current measured $R_{[\rm{OII}]}$, we check that changing the assumptions of
electron temperature from $T_{e}=10000$ to $30000$ K increases $n_{e}$ by a
factor of $\sim 1.7$.

We also estimate $n_{e}$ using pyneb (Luridiana et al., 2015) task
getCrossTemDen and combine measured $R_{[\rm{OII}]}$ from the stacked spectrum
and [O iii] $\lambda$5007/$\lambda$4363 = $40\pm 5$ line ratio measured in a
companion paper (Hsiao et al., in prep.). We obtain
$\log(n_{e})=2.9_{-0.4}^{+0.3}$, which is consistent with the value derived
above.

We highlight the fitting results of [O ii] $\lambda\lambda$3726,3729 in Figure
6. Even though JD1 and JD2 are lensed images of the same galaxy, there is an
opposite trend of $R_{[\rm{OII}]}$ between JD1 and JD2 if we examine the
median values. However, they are consistent within the range given their
uncertainties. With the MCMC method performed here, we propagate the
uncertainty from the spectra to the posterior probability distribution of
$n_{e}$ in a self-consistent way. We show the posterior distributions of
$n_{e}$ (along with the fluxes of the two lines in [O ii]
$\lambda\lambda$3726,3729 doublet, and FWHM) in Appendix B (Figure 9). The
posteriors of these parameters are converged, indicating a good fitting. The
apparent cut (i.e., envelope) in the 2D posterior between the two line fluxes
is due to the range given by the $R_{[\rm{OII}]}$ within which $n_{e}$ is
physically realistic. From the Equations A1 and A2, this range is between
$0.3729$ and $1.4556$ for $T_{e}=17000$ K. We note that this range is not
imposed from the beginning (i.e., in the fitting process).

$R_{[\rm{OII}]}$ is also a good indicator of the ISM pressure, which is
defined as the pressure within the ionized gaseous nebulae surrounding young
hot stars (e.g., Kewley et al., 2019a). Using the photoionization models of
Kewley et al. (2019a, Figure 10 therein) and assuming $12+{\rm
log(O/H)}=7.79\pm 0.09$ (Hsiao et al., in prep.), we estimate ISM pressure of
$\log(P/\rm{k})=7.3_{-0.3}^{+0.4}$.

###  4.2 Redshift evolution of $n_{e}$

Figure 7: Electron densities of galaxies measured from $z=0$ to 10. Our
measurements for JD1, JD2, and the stacked spectrum are shown by the blue,
red, and black squares, respectively. The rest of the $n_{e}$ measurements are
compiled from the literature. For local galaxies at $z\lesssim 0.2$: Berg et
al. (2022), Davies et al. (2021), Kashino et al. (2017), Kaasinen et al.
(2017), Shirazi et al. (2014), and Mingozzi et al. (2022). For $0.2\lesssim
z\lesssim 4$: Davies et al. (2021), Swinbank et al. (2019), Kashino et al.
(2017), Kaasinen et al. (2017), Harshan et al. (2020), Sanders et al. (2016),
Steidel et al. (2014), Shimakawa et al. (2015), Gburek et al. (2019), Onodera
et al. (2016), and Shirazi et al. (2014). For $4\lesssim z\lesssim 9$: Isobe
et al. (2023), Larson et al. (2023), Fujimoto et al. (2022), Stiavelli et al.
(2023), and Fujimoto et al. (2024). Different colors represent different
studies, while different symbols represent various emission-line diagnostics
used to derive $n_{e}$. Overall, there are three employed emission-line
diagnostics to derive $n_{e}$: $R_{[\rm{OII}]}$ (square), $R_{[\rm{SII}]}$
(plus), and an indirect method using [O iii] $\lambda$88$\mu$m/$\lambda$5007
ratio (cross). While [O ii] $\lambda\lambda$3726,3729 doublet is resolved in
the majority of the studies (indicated with closed squares), it is not
resolved in some studies (indicated with open squares) because of the
insufficient spectral resolution. For studies that reported $n_{e}$ of the
individual galaxies, we show those with transparent symbols, in addition to
the median values shown with thick symbols. The black dashed line represents
the best-fit power-law function in the form of $n_{e}=A\times(1+z)^{p}$ with
$A=54^{+31}_{-22}$ cm-3 and $p=1.2\pm 0.4$. The gray shaded area shows
2$\sigma$ uncertainty around the best-fit line.

Here, we put the measured electron density of into the context of galaxy
evolution by comparing it with the measurements at lower redshifts. We compile
measured electron densities at $z=0-9$ reported in the literature and plot
them in Figure 7. Different symbols represent various emission-line
diagnostics used in deriving $n_{e}$, while different colors represent
different studies. Three emission-line diagnostics have been used in the
literature to derive $n_{e}$: $R_{[\rm{OII}]}$ (square), $R_{[\rm{SII}]}$
(plus), and an indirect method using [O iii] $\lambda$88$\mu$m/$\lambda$5007
ratio (cross). Some observations do not resolve [O ii]
$\lambda\lambda$3726,3729 doublet because of the insufficient spectral
resolution. We indicate the results from those studies with open squares. For
studies that reported $n_{e}$ of the individual galaxies, we show those values
with transparent symbols, in addition to the median values shown with thick
symbols. The measured $n_{e}$ of JD1, JD2, and stacked spectrum are shown with
the blue, red, and black squares, respectively. The rest of the data points
are listed in the following, classified based on the redshift range.

  * •

$z\lesssim 0.2$: Berg et al. (2022) measured $n_{e}$ of 45 local star-forming
galaxies at $0.002<z<0.182$ using data from the CLASSY survey (Berg et al.,
2022; James et al., 2022). They derived $n_{e}$ based on $R_{[\rm{SII}]}$
$=$[S ii] $\lambda$6716/$\lambda$6731 and obtained a wide range of $n_{e}$
$\sim 10-900$ cm-3, with a median of $100$ cm-3. Using the same CLASSY sample,
Mingozzi et al. (2022) explored multiple UV diagnostics for measuring the
properties of ISM, including $n_{e}$, in which they applied multiple line
ratios representing low- and high-ionization regions. We take $n_{e}$ measured
based on $R_{[\rm{SII}]}$ whose $T_{e}$ is derived based on [O iii]
$\lambda$4363/$\lambda$5007\. Davies et al. (2021) measured $n_{e}$ of
galaxies in a wide redshift range by stacking the spectra of the galaxies in
four redshift bins: $z\lesssim 0.1$ (471), $0.6\lesssim z\lesssim 1.1$ (39),
$1.1\lesssim z\lesssim 1.9$ (36), and $1.9\lesssim z\lesssim 2.6$ (65). The
data for local galaxies are taken from SAMI survey (Croom et al., 2012). They
estimated $n_{e}=32_{-9}^{+4}$ cm-3 from $R_{[\rm{SII}]}$. Kashino et al.
(2017) and Kaasinen et al. (2017) measured $n_{e}$ of galaxies at $z\sim 1.5$
and $z\sim 0$. For the local sample, they used data from the SDSS Data Release
7 (Abazajian et al., 2009), specifically the MPA-JHU catalogs (Kauffmann et
al., 2003; Brinchmann et al., 2004). For the local sample, Kaasinen et al.
(2017) obtained a median of $n_{e}$ $26.8_{-0.2}^{+0.2}$ cm-3, while Kashino
et al. (2017) obtained $n_{e}$ in the range of $\sim 10-100$ cm-3. Shirazi et
al. (2014) studied the ISM properties in 14 star-forming galaxies at
$2.6<z<3.4$ and local galaxies matched in stellar mass and sSFR. They used
emission lines from the MPA-JHU catalogs for the local sample, and estimated
$n_{e}$ based on $R_{[\rm{SII}]}$, finding a range of $\sim 60-250$ cm-3.

  * •

$0.2\lesssim z\lesssim 4$: Davies et al. (2021) used data from the KMOS3D
survey (Wisnioski et al., 2015) to measure $n_{e}$ of galaxies in 3 high-$z$
bins (mentioned above). Based on the $R_{[\rm{SII}]}$ of the stacked spectra,
they estimated $n_{e}=101^{+59}_{-85}$, $79^{+120}_{40}$, and
$187^{+140}_{-132}$ cm-3 for the samples at $z\sim 0.9$, $\sim 1.5$, and $\sim
2.2$, respectively. Swinbank et al. (2019) used data from the KROSS survey
(Stott et al., 2016) and performed stacking on the spectra of 529 star-forming
galaxies at $0.6<z<1.0$. They estimated $n_{e}$ from $R_{[\rm{SII}]}$ in the
stacked spectrum and decomposed it into broad-line and narrow-line components.
They estimated $n_{e}$ $=30_{-20}^{+40}$ and $75_{-50}^{+55}$ cm-3 for the
narrow-line and broad-line components, respectively. We take the mean of the
two for our analysis. Kashino et al. (2017) used data from the FMOS-COSMOS
survey (Silverman et al., 2015) to study the ISM properties of galaxies at
$1.4\lesssim z\lesssim 1.7$. They measured an average $n_{e}$
$=220_{-130}^{+170}$ cm-3 based on $R_{[\rm{SII}]}$. Kaasinen et al. (2017)
measured $n_{e}$ of 21 star-forming galaxies at $1.4\lesssim z\lesssim 1.7$
using data from the COSMOS-[OII] survey. The estimated $n_{e}$ based on
$R_{[\rm{OII}]}$ ranges from $\sim 9$ to $1000$ cm-3, with a median of
$114^{+28}_{-27}$ cm-3. Harshan et al. (2020) used data from the ZFIRE survey
(Nanayakkara et al., 2016) to calculate $n_{e}$ of $z\sim 1.5$ galaxies, of
which 6 are in the cluster environment and 2 are in the field. They found
$n_{e}$ in the range of $\sim 57-490$ cm-3 for the cluster galaxies and $\sim
49-160$ cm-3 for the field galaxies. Sanders et al. (2016) measured $n_{e}$ of
97 galaxies at $2.0<z<2.6$ using data from the MOSDEF survey (Kriek et al.,
2015). They obtained a median $n_{e}$ $=225^{+119}_{-4}$ cm-3 based on
$R_{[\rm{OII}]}$. Steidel et al. (2014) estimated $n_{e}$ $=220^{+380}_{-160}$
cm-3 from $R_{[\rm{OII}]}$ of a stacked spectrum of 113 galaxies at
$2.0<z<2.6$. Shimakawa et al. (2015) measured $n_{e}$ of 14 H$\alpha$-emitting
galaxies at $z\sim 2.5$ based on $R_{[\rm{OII}]}$. They obtained a median of
$n_{e}$ $\sim 290$ cm-3. Gburek et al. (2019) measured $n_{e}$ $\sim 220$ cm-3
in a lensed dwarf galaxy at $z=2.59$ based on $R_{[\rm{OII}]}$. Onodera et al.
(2016) derived $n_{e}$ of star-forming galaxies at $3\lesssim z\lesssim 3.7$
based on $R_{[\rm{OII}]}$, finding a median $\sim 420$ cm-3. Shirazi et al.
(2014) obtained $n_{e}$ $\sim 110-2700$ cm-3 in star-forming galaxies at
$2.6<z<3.4$ based on $R_{[\rm{SII}]}$.

  * •

$4\lesssim z\lesssim 9$: Isobe et al. (2023) estimated $n_{e}$ of 14 star-
forming galaxies at $4.02<z<8.68$ using data from public JWST surveys that
used the NIRSpec medium- and high-resolution spectroscopy. They observed [O
ii] $\lambda\lambda$3726,3729 doublet resolved in only three galaxies, but
blended in the other galaxies. They carefully modeled the line spread function
(LSF) of the NIRSpec to deconvolve the [O ii] $\lambda\lambda$3726,3729
doublets. They obtained $n_{e}$ range of $\sim 39-3100$ cm-3. Larson et al.
(2023) studied the ISM properties in a galaxy CEERS_1019 at $z=8.679$ using
NIRSpec medium-resolution data. They carefully fit an unresolved [O ii]
$\lambda\lambda$3726,3729 doublet with two Gaussian models and derived $n_{e}$
$1.9\pm 0.2\times 10^{3}$ cm-3. Fujimoto et al. (2022) measured $n_{e}$ of a
galaxy at $z=8.496$ using indirect method based on the [O iii]
88$\mu$m/$\lambda$5007 line ratio and obtained $n_{e}$ $=220_{-130}^{+230}$
cm-3. Using the same method as Fujimoto et al. (2022), Stiavelli et al. (2023)
estimated $n_{e}$ of MACS1149-JD1 at $z=9.11$ to be
$\log(n_{e})=2.60_{-0.27}^{+0.25}$. With the [S ii] $\lambda\lambda$6716,6731
doublets, Fujimoto et al. (2024) evaluated $n_{e}$ of a sub-$L^{\star}$ main-
sequence lensed galaxy at $z=6.07$, the Cosmic Grapes, to be
$n_{e}=260^{+400}_{-230}$ cm-3.

As shown in Figure 7, there is a clear trend of increasing $n_{e}$ with
redshift from $z=0$ to $10$. To get a functional form out of the evolutionary
trend, we fit the data with a power-law function in the form of
$n_{e}=A\times(1+z)^{p}$ using the MCMC method. For the literature that
reports electron densities of multiple individual galaxies in a redshift bin,
we only use the median value for the fitting. For , we only use the electron
density derived from the stacked spectrum. We obtain $A=54^{+31}_{-22}$ cm-3
and $p=1.2^{+0.4}_{-0.4}$. The median function obtained from the fitting
result is shown with the black dashed line, while the gray shaded area shows
2$\sigma$ uncertainty. The fitting result can track the bulk of the data
distribution quite well. Best-fit $A$ suggests an electron density of $\sim
50$ cm-3 in the local universe. The $n_{e}$ of JD1, JD2, and stacked spectrum
falls within the range of the best-fit line and 2$\sigma$ uncertainty.

Previous studies have also hinted at a redshift evolution of electron density
in the form of a power-law function, $n_{e}\propto(1+z)^{p}$ with $p\sim 1-2$
(e.g., Davies et al., 2021; Isobe et al., 2023). Isobe et al. (2023) compiled
electron densities reported from the literature at $z\lesssim 3$ and combined
with their measurements at $4\lesssim z\lesssim 9$. They found that the
majority of the data lie within the range given by the power-law function with
$p=1$ and $2$. They further investigate the dependence of $n_{e}$ on the
$M_{*}$, SFR, and sSFR, and found that $n_{e}$ seems to increase toward higher
redshift at a given $M_{*}$, SFR, and sSFR, suggesting that $n_{e}$ evolution
is not an implication of the evolution of those global properties. This
somewhat contradicts the finding of Kaasinen et al. (2017) that shows a weak
evolution of $n_{e}$ with redshift for a given SFR, suggesting that the
redshift evolution of $n_{e}$ is mainly due to the cosmic evolution of SFR. To
further investigate the driving factors of $n_{e}$ evolution, we need more
data that cover a wide range of redshift and homogeneous methods in the
measurements of the physical properties, which is beyond the scope of our
current study. We do not investigate the dependence of $n_{e}$ on the global
properties of galaxies because of the inhomogeneity of the methods for
deriving those properties used in the literature.

The values $n_{e}\approx 10^{3}$ cm-3 measured in JD1 and JD2 are in good
agreement with the typical density of the ISM in the vicinity of massive
stars, within star-forming regions, found in cosmological simulations of
$z\sim 10$ galaxies (e.g., see Fig. 12 in Sugimura et al., 2024). However,
understanding the dominant driver of the redshift evolution of $n_{e}$ is more
complicated. Although the mean density of the IGM and virialized halos scale
as $(1+z)^{3}$, the phase of the ISM probed by the [O ii]
$\lambda\lambda$3726,3729 emission is ionized gas around newly formed massive
stars, which is several orders of magnitude denser than the mean gas density
in halos. One intriguing interpretation of the redshift evolution of $n_{e}$
is in terms of the mean metallicity evolution of galaxies with redshift. This
suggestion is based on two observations: i) Simulations of high-redshift
galaxies find that the mean density of star-forming clumps roughly scales as
$\overline{n}\propto Z^{-1}$, a result that is a direct consequence of their
higher temperatures in lower-metallicity gas clouds (see Fig. 3 in Garcia et
al., 2023); ii) The mean metallicity of the ISM in star-forming galaxies at
redshifts $z\sim 0-3$, according to simulations (Davé et al., 2011) and
observations (Yuan et al., 2012), decreases with increasing redshift roughly
as $Z\propto(1+z)^{-1}$, hence the relationship
$n_{e}\propto\overline{n}\propto(1+z)$.

Alternatively, the $n_{e}$ evolution in the form of $n_{e}\propto(1+z)^{p}$
may be related to the morphological evolution of galaxies which influences the
density of the ISM. The increasing ISM density toward high redshift can be
caused by the size evolution of galaxies, which are more compact at high
redshift. The galaxy size has been shown to decrease with redshift in the form
of $r\propto(1+z)^{\sim-1}$ (e.g., Shibuya et al., 2015; Ono et al., 2023;
Ormerod et al., 2024). The virial radius of the dark matter halo also evolves
at a similar rate (e.g., Mo & White, 2002). Assuming disk morphology, which is
compatible with star-forming galaxies, the densities of stellar mass and dark
matter halo mass are expected to evolve as $r^{-2}\propto(1+z)^{\sim 2}$. If
$n_{e}$ is expected to evolve following the stellar mass and halo mass
densities, it evolves as $n_{e}\propto(1+z)^{\sim 2}$. Our result suggests a
slightly lower power index of $1.2^{+0.5}_{-0.4}$.

##  5 Summary

We present JWST/NIRSpec high-resolution spectroscopy (G395H/F290LP grating) of
, a triply-lensed galaxy system at $z=10.167$. We obtained spectra for the
MACS0647–JD1 and –JD2. The new spectra detect and resolve emission lines in
the rest-frame ultraviolet (UV) and blue optical, including resolved [O ii]
$\lambda\lambda$3726,3729 doublet, [Ne iii] $\lambda$3870, HeI $\lambda$3890,
H$\delta$, H$\gamma$, and [O iii] $\lambda$4363\. The [O ii]
$\lambda\lambda$3726,3729 is resolved for the first time at $z>8$, revealing
the two lines and providing an unprecedented opportunity to directly measure
electron density ($n_{e}$) in a galaxy at early cosmic time. We stack the
spectra of JD1 and JD2 to get a higher overall signal-to-noise ratio. We
measure $n_{e}$ from the [O ii] $\lambda$3729/$\lambda$3726 and obtained
$\log(n_{e}/\rm{cm}^{-3})=3.3^{+0.6}_{-0.5}$, $2.7^{+0.5}_{-0.6}$, and
$2.9^{+0.5}_{-0.5}$ for JD1, JD2, and stacked spectrum, respectively. We
compare our $n_{e}$ measurements with those reported in the literature for
galaxies in lower redshifts and study the evolutionary trend of $n_{e}$. We
see a clear trend of increasing $n_{e}$ with redshift. The redshift evolution
follows a power-law function in the form of $n_{e}=A\times(1+z)^{p}$ with
$A=54^{+31}_{-22}$ cm-3 and $p=1.2^{+0.4}_{-0.4}$. This power-law form may be
explained as the result of decreasing gas-phase metallicity of galaxies with
increasing redshift, the morphological evolution of galaxies, or the
combination of the two. It has been known that lower-metallicity star-forming
clouds tend to be denser and the morphological evolution of galaxies happened
in such a way that galaxies were more compact at higher redshift.

##  6 Acknowledgments

This work is based on observations made with the NASA/ESA/CSA James Webb Space
Telescope (JWST). The data were obtained from the Mikulski Archive for Space
Telescopes (MAST) at the Space Telescope Science Institute (STScI), which is
operated by the Association of Universities for Research in Astronomy (AURA),
Inc., under NASA contract NAS 5-03127 for JWST. These observations are
associated with programs JWST GO 4246 and 1433.

A is funded by a grant for JWST-GO-01433 and JWST-GO-04246 provided by STScI
under NASA contract NAS5-03127. TYYH appreciates the support from the
Government scholarship to study abroad (Taiwan). We are grateful and indebted
to the 20,000 people who worked to make JWST an incredible discovery machine.
AA acknowledges support by the Swedish research council Vetenskapsrådet
(2021-05559). PD acknowledges support from the NWO grant 016.VIDI.189.162
(“ODIN”) and warmly thanks the European Commission’s and University of
Groningen’s CO-FUND Rosalind Franklin program. M.K. was supported by the ANID
BASAL project FB210003. AZ and LJF acknowledge support by Grant No. 2020750
from the United States-Israel Binational Science Foundation (BSF) and Grant
No. 2109066 from the United States National Science Foundation (NSF); by the
Ministry of Science & Technology, Israel; and by the Israel Science Foundation
Grant No. 864/23.

##  Appendix A General formula for the relationship between $R_{[\rm{OII}]}$
and $n_{e}$

We perform interpolation to the model grids of Kewley et al. (2019a) to
construct a general functional formula for calculating $n_{e}$ from
$R_{[\rm{OII}]}$ given an electron temperature ($T_{e}$). The model grids of
Kewley et al. (2019a) give the relationship between $n_{e}$ and
$R_{[\rm{OII}]}$ for three different electron temperatures ($\log(T_{e}/\rm
K)=3.5,4.0,$ and $4.5$). The model grids were generated using the MAPPINGS
version 5.1 photoionization code (Binette et al., 1985; Sutherland & Dopita,
1993; Dopita et al., 2015). To get a functional formula out of the curves, we
fit each $R_{\rm[OII]}(n_{e})$ curve with the following functional form
adapted from Sanders et al. (2016):

$n_{e}=\frac{\gamma R_{\rm[OII]}-\alpha\beta}{\alpha-R_{\rm[OII]}}.$ (A1)

We then interpolate the $R_{\rm[OII]}(n_{e})$ curves using the cubic spline
method for 9 values of $\log(T_{e})$ between $3.5$ and $4.5$. For each
interpolated curve, we perform fitting with Eq. A1 to get the coefficients
($\alpha$, $\beta$, and $\gamma$) and investigate their dependencies with
$T_{e}$. We show the relations between those coefficients and $T_{e}$ in
Figure 8. The blue crosses represent the coefficients of the original curves
from Kewley et al. (2019a), while the gray circles are from the interpolated
curves. The relations can be well-fitted with a cubic polynomial function of
the form

$\alpha,\beta,\gamma=a+bx+cx^{2}+dx^{3}$ (A2)

with $x=\log(T_{e})$. The best-fit cubic polynomial functions are shown as red
dashed lines in Figure 8 and their coefficients are summarized in Table 4. We
show models $R_{\rm[OII]}(n_{e})$ for 9 values of $T_{e}$ in Figure 5.

Table 4: Coefficients for the polynomial function in Eq. A2 Coefficient | $a$ | $b$ | $c$ | $d$
---|---|---|---|---
$\alpha$ | $-0.00447$ | $0.28813$ | $-0.06841$ | $0.00506$
$\beta$ | $-7616.41042$ | $7075.85077$ | $-2433.11054$ | $321.56781$
$\gamma$ | $-1661.48935$ | $1481.10495$ | $-507.35920$ | $69.67036$
Figure 8: The dependence of the coefficients in Equation A1 to the electron
temperature ($T_{e}$; in logarithmic scale). The blue crosses represent the
coefficients of the original $R_{\rm[OII]}(n_{e})$ curves from Kewley et al.
(2019a) at $\log(T_{e}/\rm{K})=3.5,4.0,$ and $4.5$. The gray circles show the
coefficients of the interpolated curves at 9 values of $T_{e}$. The red dashed
lines show the best-fit cubic polynomial functions.

##  Appendix B Corner plots of the fitting on [O ii] $\lambda\lambda$3726,3729
doublet

Figure 9 shows the corner plots showing the posterior probability
distributions from the MCMC fitting of the [O ii] $\lambda\lambda$3726,3729
doublet of JD1 (top left), JD2 (top right), and the stacked spectrum (bottom).
Only some parameters are shown in the corner plots, including the individual
line fluxes of [O ii] $\lambda\lambda$3726,3729 doublet, FWHM, and electron
density. Among those parameters, only FWHM is a free parameter, while the
others are dependent parameters. Other free parameters, including the peak
wavelength and amplitude, are not shown in the corner plots.

Figure 9: Corner plots showing the posterior probability distributions of some
parameters in the MCMC fitting of the [O ii] $\lambda\lambda$3726,3729
doublet. The top left, top right, and bottom panels show the corner plots for
JD1, JD2, and stacked spectrum, respectively. Other free parameters, including
the central wavelength and amplitude, are not shown here. The fluxes of [O ii]
$\lambda$3726 and [O ii] $\lambda$3729 are dependent parameters.

##  Appendix C Measured emission lines from the NIRSpec prism spectra

The NIRSpec prism spectra of JD1 and JD2 were obtained by JWST program GO 1433
(PI Coe) and presented in Hsiao et al. (2023a). There were two separate
NIRSpec multi-object spectroscopy (MOS) observations (Obs 21 and Obs 23) using
the low-resolution ($R\sim 30-300$) prism spanning a wavelength coverage of
$0.6-5.3$ $\mu$m. Obs 23 was performed with standard 3-slitlet nods, while Obs
21 was performed with single slitlets. Each observation had a total exposure
time of 1.8 hours (3.6 hours). In Hsiao et al. (2023a), the emission lines
were measured from the stacked spectrum of the four spectra (JD1 Obs 21, JD1
Obs 23, JD2 Obs 21, and JD2 Obs 23). Here, we revisit the analysis and measure
the emission lines in the individual spectra. Moreover, we use the latest
reduced spectra from The Dawn JWST Archive (DJA)444https://dawn-
cph.github.io/dja/, which were reduced using the upgraded
msaexp555https://github.com/gbrammer/msaexp version 0.6.7. The NIRSpec data
reduction is described in Heintz et al. (2023b) and Heintz et al. (2024). We
fit the NIRSpec prism spectra and measure the emission line fluxes using
piXedfit. Besides the individual spectra, we also fit the stacked sum of the
four spectra (i.e., JD1 Obs 21 + JD1 Obs 23 + JD2 Obs 21 + JD2 Obs 23).

The fitting results for the individual spectra are shown in Figure 10
alongside the slitlets configurations overlaid on $1.6^{\prime\prime}\times
1.6^{\prime\prime}$ NIRCam color images. We show the fitting result of the
stacked spectrum in Figure 11. We note that most of the observations target
the brighter component A. Only Obs 21 on JD2 covers both components A and B.
We also note that slitlets do not cover the same area in the two observations,
causing the differences in emission line fluxes. The black line and gray
shaded area represent the observed spectra and their 1$\sigma$ uncertainties.
The blue and red lines represent the best-fit continuum and Gaussian models of
the detected emission lines. The spectra roll off gradually from $\sim 1.5$
$\mu$m (1343 Årest-frame) to the Lyman break observed at $1.36$ $\mu$m (1216
Å), instead of a sharp break. This is due to the Lyman-$\alpha$ damping wing
(see Hsiao et al., 2023a; Heintz et al., 2023b). The continuum modeled here
does not include the modeling of Lyman-$\alpha$ damping wing which makes the
discrepancy between the models and observed spectra around the wavelength
range. The best-fit continuum is obtained from SED fitting on the spectrum
after masking the emission lines and cutting at $1.5<\lambda<5.2$ $\mu$m.

Table 5 summarizes the measured emission line properties. The line fluxes are
not corrected for the lensing magnifications ($8.0$ and $5.3$ for JD1 and JD2,
respectively) and slit loss. The slit loss corrections can be estimated by
taking the ratio between the total aperture photometry in NIRCam filters,
obtained from the grizli pipeline (Brammer et al., 2022) (see Hsiao et al.,
2023a), and pseudo photometry obtained by integrating the NIRSpec spectra
through the NIRCam filters. Then the correction factors at the emission lines
can be estimated by interpolating the wavelength-dependent flux ratio. The
correction factors for the listed lines in Table 5 are the following. JD1 Obs
21 = [2.5, 3.7, 3.9, 4.0, 4.2, 4.6, 4.6, 4.6], JD1 Obs 23 = [2.2, 3.3, 3.4,
3.5, 3.7, 4.0, 4.0, 4.0], JD2 Obs 21 = [2.7, 2.4, 2.4, 2.4, 2.3, 2.0, 2.0,
2.0], and JD2 Obs 23 = [2.0, 2.5, 2.6, 2.6, 2.7, 2.8, 2.8, 2.8].

Figure 10: NIRSpec prism spectra of JD1 and JD2 from two separate observations
(Obs 21 and Obs 23). The observations were performed by the JWST program GO
1433 (PI Coe) and the original data were presented in Hsiao et al. (2023a).
Left: Slitlet configurations overlaid on $1.6^{\prime\prime}\times
1.6^{\prime\prime}$ NIRCam color images. Note that most of the observations
target the brighter component A. Only JD2 Obs 21 covers both components A and
B. Right: NIRSpec prism spectra (black line and gray shaded area) and the
best-fit model continuum (blue) and Gaussian for the emission lines (red). The
SED modeling does not include the Lyman-$\alpha$ damping wing which makes the
discrepancy between the best-fit continuum and the observed spectrum around
the Lyman-$\alpha$ wavelength. Figure 11: The stacked sum of the four spectra
shown in Figure 10 (i.e., JD1 Obs 21 + JD1 Obs 23 + JD2 Obs 21 + JD2 Obs 23)
and its fitting result. The line and symbols are the same as those in Figure
10. Table 5: Measured Emission Line Properties from NIRSpec prism spectra

Emission Line |  $\lambda$ rest | JD1 Obs 21 Flux | JD1 Obs 23 Flux | JD2 Obs 21 Flux | JD2 Obs 23 Flux | Stacked Flux
---|---|---|---|---|---|---
(1) | (2) | (3) | (4) | (5) | (6) | (7)
[C iii] | 1906.68,1908.73 | $95_{-14}^{+15}$ | $98_{-21}^{+21}$ | $194_{-17}^{+17}$ | $129_{-20}^{+21}$ | $428_{-35}^{+34}$
[O ii] | 3727.10,3729.86 | $17_{-4}^{+5}$ | $34_{-7}^{+7}$ | $50_{-5}^{+5}$ | $30_{-5}^{+5}$ | $144_{-12}^{+12}$
[Ne iii] | 3869.86 | $19_{-5}^{+6}$ | $35_{-5}^{+6}$ | $54_{-5}^{+5}$ | $51_{-7}^{+8}$ | $177_{-12}^{+12}$
[Ne iii] | 3968.59 |  $\cdots$ | $14_{-6}^{+8}$ | $33_{-4}^{+5}$ | $20_{-5}^{+5}$ | $88_{-10}^{+10}$
H$\delta$ | 4102.89 | $15_{-5}^{+5}$ | $34_{-7}^{+7}$ | $27_{-4}^{+5}$ | $19_{-5}^{+5}$ | $117_{-17}^{+55}$
H$\gamma$ | 4341.69 | $21_{-4}^{+4}$ | $58_{-8}^{+9}$ | $65_{-6}^{+6}$ | $34_{-6}^{+6}$ | $179_{-12}^{+13}$
[O iii] | 4364.44 |  $\cdots$ | $12_{-3}^{+3}$ | $13_{-3}^{+3}$ | $28_{-7}^{+7}$ | $62_{-5}^{+5}$
He i | 4472.735 |  $\cdots$ |  $\cdots$ | $26_{-8}^{+9}$ | $12_{-6}^{+6}$ | $36_{-12}^{+13}$

Note. — (1) Emission line. (2) Rest-frame wavelength. (3) Measured line fluxes
of the JD1 Obs 21 spectrum. (4) Measured line fluxes of the JD1 Obs 23
spectrum. (5) Measured line fluxes of the JD2 Obs 21 spectrum. (6) Measured
line fluxes of the JD2 Obs 23 spectrum. (7) Measured line fluxes from the
stacked sum of the four spectra (JD1 Obs21 + JD1 Obs23 + JD2 Obs21 + JD2
Obs23). All flux densities are in units of $10^{-20}$ erg s${}^{-1}\,$cm-2
Å-1. The line fluxes are not corrected for the lensing magnifications (8.0 and
5.3 for JD1 and JD2, respectively) and slit loss.

## References

  * Abazajian et al. (2009) Abazajian, K. N., Adelman-McCarthy, J. K., Agüeros, M. A., et al. 2009, ApJS, 182, 543, doi: 10.1088/0067-0049/182/2/543
  * Abdurro’uf et al. (2021) Abdurro’uf, Lin, Y.-T., Wu, P.-F., & Akiyama, M. 2021, ApJS, 254, 15, doi: 10.3847/1538-4365/abebe2
  * Abdurro’uf et al. (2022) —. 2022, piXedfit: Analyze spatially resolved SEDs of galaxies, Astrophysics Source Code Library, record ascl:2207.033
  * Abdurro’uf et al. (2023) Abdurro’uf, Coe, D., Jung, I., et al. 2023, ApJ, 945, 117, doi: 10.3847/1538-4357/acba06
  * Backhaus et al. (2024) Backhaus, B. E., Trump, J. R., Pirzkal, N., et al. 2024, ApJ, 962, 195, doi: 10.3847/1538-4357/ad1520
  * Bacon et al. (2017) Bacon, R., Conseil, S., Mary, D., et al. 2017, A&A, 608, A1, doi: 10.1051/0004-6361/201730833
  * Berg et al. (2022) Berg, D. A., James, B. L., King, T., et al. 2022, ApJS, 261, 31, doi: 10.3847/1538-4365/ac6c03
  * Binette et al. (1985) Binette, L., Dopita, M. A., & Tuohy, I. R. 1985, ApJ, 297, 476, doi: 10.1086/163544
  * Böker et al. (2023) Böker, T., Beck, T. L., Birkmann, S. M., et al. 2023, PASP, 135, 038001, doi: 10.1088/1538-3873/acb846
  * Brammer et al. (2022) Brammer, G., Strait, V., Matharu, J., & Momcheva, I. 2022, grizli, 1.5.0, Zenodo, doi: 10.5281/zenodo.6672538
  * Brinchmann et al. (2004) Brinchmann, J., Charlot, S., White, S. D. M., et al. 2004, MNRAS, 351, 1151, doi: 10.1111/j.1365-2966.2004.07881.x
  * Bunker et al. (2023) Bunker, A. J., Saxena, A., Cameron, A. J., et al. 2023, A&A, 677, A88, doi: 10.1051/0004-6361/202346159
  * Cameron et al. (2023) Cameron, A. J., Saxena, A., Bunker, A. J., et al. 2023, A&A, 677, A115, doi: 10.1051/0004-6361/202346107
  * Castellano et al. (2024) Castellano, M., Napolitano, L., Fontana, A., et al. 2024, arXiv e-prints, arXiv:2403.10238, doi: 10.48550/arXiv.2403.10238
  * Christensen et al. (2023) Christensen, L., Jakobsen, P., Willott, C., et al. 2023, A&A, 680, A82, doi: 10.1051/0004-6361/202347943
  * Coe et al. (2013) Coe, D., Zitrin, A., Carrasco, M., et al. 2013, ApJ, 762, 32, doi: 10.1088/0004-637X/762/1/32
  * Croom et al. (2012) Croom, S. M., Lawrence, J. S., Bland-Hawthorn, J., et al. 2012, MNRAS, 421, 872, doi: 10.1111/j.1365-2966.2011.20365.x
  * Curti et al. (2023) Curti, M., D’Eugenio, F., Carniani, S., et al. 2023, MNRAS, 518, 425, doi: 10.1093/mnras/stac2737
  * Davies et al. (2021) Davies, R. L., Förster Schreiber, N. M., Genzel, R., et al. 2021, ApJ, 909, 78, doi: 10.3847/1538-4357/abd551
  * Davé et al. (2011) Davé, R., Finlator, K., & Oppenheimer, B. D. 2011, Monthly Notices of the Royal Astronomical Society, 416, 1354, doi: 10.1111/j.1365-2966.2011.19132.x
  * Dopita et al. (2006a) Dopita, M. A., Fischera, J., Sutherland, R. S., et al. 2006a, ApJS, 167, 177, doi: 10.1086/508261
  * Dopita et al. (2006b) —. 2006b, ApJ, 647, 244, doi: 10.1086/505418
  * Dopita et al. (2015) Dopita, M. A., Ho, I. T., Dressel, L. L., et al. 2015, ApJ, 801, 42, doi: 10.1088/0004-637X/801/1/42
  * Ferruit et al. (2022) Ferruit, P., Jakobsen, P., Giardino, G., et al. 2022, A&A, 661, A81, doi: 10.1051/0004-6361/202142673
  * Finkelstein et al. (2023) Finkelstein, S. L., Bagley, M. B., Ferguson, H. C., et al. 2023, ApJ, 946, L13, doi: 10.3847/2041-8213/acade4
  * Foreman-Mackey et al. (2013) Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306, doi: 10.1086/670067
  * Fujimoto et al. (2022) Fujimoto, S., Ouchi, M., Nakajima, K., et al. 2022, arXiv e-prints, arXiv:2212.06863, doi: 10.48550/arXiv.2212.06863
  * Fujimoto et al. (2023) Fujimoto, S., Arrabal Haro, P., Dickinson, M., et al. 2023, ApJ, 949, L25, doi: 10.3847/2041-8213/acd2d9
  * Fujimoto et al. (2024) Fujimoto, S., Ouchi, M., Kohno, K., et al. 2024, arXiv e-prints, arXiv:2402.18543, doi: 10.48550/arXiv.2402.18543
  * Garcia et al. (2023) Garcia, F. A. B., Ricotti, M., Sugimura, K., & Park, J. 2023, MNRAS, 522, 2495, doi: 10.1093/mnras/stad1092
  * Gardner et al. (2023) Gardner, J. P., Mather, J. C., Abbott, R., et al. 2023, PASP, 135, 068001, doi: 10.1088/1538-3873/acd1b5
  * Gburek et al. (2019) Gburek, T., Siana, B., Alavi, A., et al. 2019, ApJ, 887, 168, doi: 10.3847/1538-4357/ab5713
  * Hainline et al. (2024) Hainline, K. N., D’Eugenio, F., Jakobsen, P., et al. 2024, arXiv e-prints, arXiv:2404.04325, doi: 10.48550/arXiv.2404.04325
  * Harris et al. (2020) Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357, doi: 10.1038/s41586-020-2649-2
  * Harshan et al. (2020) Harshan, A., Gupta, A., Tran, K.-V., et al. 2020, ApJ, 892, 77, doi: 10.3847/1538-4357/ab76cf
  * Heintz et al. (2023a) Heintz, K. E., Brammer, G. B., Giménez-Arteaga, C., et al. 2023a, Nature Astronomy, 7, 1517, doi: 10.1038/s41550-023-02078-7
  * Heintz et al. (2023b) Heintz, K. E., Watson, D., Brammer, G., et al. 2023b, arXiv e-prints, arXiv:2306.00647, doi: 10.48550/arXiv.2306.00647
  * Heintz et al. (2024) Heintz, K. E., Brammer, G. B., Watson, D., et al. 2024, arXiv e-prints, arXiv:2404.02211, doi: 10.48550/arXiv.2404.02211
  * Horne (1986) Horne, K. 1986, PASP, 98, 609, doi: 10.1086/131801
  * Hsiao et al. (2023a) Hsiao, T. Y.-Y., Abdurro’uf, Coe, D., et al. 2023a, arXiv e-prints, arXiv:2305.03042, doi: 10.48550/arXiv.2305.03042
  * Hsiao et al. (2023b) Hsiao, T. Y.-Y., Coe, D., Abdurro’uf, et al. 2023b, ApJ, 949, L34, doi: 10.3847/2041-8213/acc94b
  * Hunter (2007) Hunter, J. D. 2007, Computing in Science & Engineering, 9, 90, doi: 10.1109/MCSE.2007.55
  * Isobe et al. (2023) Isobe, Y., Ouchi, M., Nakajima, K., et al. 2023, ApJ, 956, 139, doi: 10.3847/1538-4357/acf376
  * Jakobsen et al. (2022) Jakobsen, P., Ferruit, P., Alves de Oliveira, C., et al. 2022, A&A, 661, A80, doi: 10.1051/0004-6361/202142663
  * James et al. (2022) James, B. L., Berg, D. A., King, T., et al. 2022, ApJS, 262, 37, doi: 10.3847/1538-4365/ac8008
  * Jung et al. (2023) Jung, I., Finkelstein, S. L., Arrabal Haro, P., et al. 2023, arXiv e-prints, arXiv:2304.05385, doi: 10.48550/arXiv.2304.05385
  * Kaasinen et al. (2017) Kaasinen, M., Bian, F., Groves, B., Kewley, L. J., & Gupta, A. 2017, MNRAS, 465, 3220, doi: 10.1093/mnras/stw2827
  * Kashino et al. (2017) Kashino, D., Silverman, J. D., Sanders, D., et al. 2017, ApJ, 835, 88, doi: 10.3847/1538-4357/835/1/88
  * Kauffmann et al. (2003) Kauffmann, G., Heckman, T. M., Tremonti, C., et al. 2003, MNRAS, 346, 1055, doi: 10.1111/j.1365-2966.2003.07154.x
  * Kewley & Dopita (2002) Kewley, L. J., & Dopita, M. A. 2002, ApJS, 142, 35, doi: 10.1086/341326
  * Kewley et al. (2019a) Kewley, L. J., Nicholls, D. C., Sutherland, R., et al. 2019a, ApJ, 880, 16, doi: 10.3847/1538-4357/ab16ed
  * Kewley et al. (2019b) Kewley, L. J., Nicholls, D. C., & Sutherland, R. S. 2019b, ARA&A, 57, 511, doi: 10.1146/annurev-astro-081817-051832
  * Killi et al. (2023) Killi, M., Watson, D., Fujimoto, S., et al. 2023, MNRAS, 521, 2526, doi: 10.1093/mnras/stad687
  * Kriek et al. (2015) Kriek, M., Shapley, A. E., Reddy, N. A., et al. 2015, ApJS, 218, 15, doi: 10.1088/0067-0049/218/2/15
  * Landsman (1993) Landsman, W. B. 1993, in Astronomical Society of the Pacific Conference Series, Vol. 52, Astronomical Data Analysis Software and Systems II, ed. R. J. Hanisch, R. J. V. Brissenden, & J. Barnes, 246
  * Larson et al. (2023) Larson, R. L., Finkelstein, S. L., Kocevski, D. D., et al. 2023, ApJ, 953, L29, doi: 10.3847/2041-8213/ace619
  * Le Fèvre et al. (2015) Le Fèvre, O., Tasca, L. A. M., Cassata, P., et al. 2015, A&A, 576, A79, doi: 10.1051/0004-6361/201423829
  * Luridiana et al. (2015) Luridiana, V., Morisset, C., & Shaw, R. A. 2015, A&A, 573, A42, doi: 10.1051/0004-6361/201323152
  * Mingozzi et al. (2022) Mingozzi, M., James, B. L., Arellano-Córdova, K. Z., et al. 2022, ApJ, 939, 110, doi: 10.3847/1538-4357/ac952c
  * Mo & White (2002) Mo, H. J., & White, S. D. M. 2002, MNRAS, 336, 112, doi: 10.1046/j.1365-8711.2002.05723.x
  * Nakajima et al. (2023) Nakajima, K., Ouchi, M., Isobe, Y., et al. 2023, ApJS, 269, 33, doi: 10.3847/1538-4365/acd556
  * Nanayakkara et al. (2016) Nanayakkara, T., Glazebrook, K., Kacprzak, G. G., et al. 2016, ApJ, 828, 21, doi: 10.3847/0004-637X/828/1/21
  * Oesch et al. (2016) Oesch, P. A., Brammer, G., van Dokkum, P. G., et al. 2016, ApJ, 819, 129, doi: 10.3847/0004-637X/819/2/129
  * Ono et al. (2023) Ono, Y., Harikane, Y., Ouchi, M., et al. 2023, ApJ, 951, 72, doi: 10.3847/1538-4357/acd44a
  * Onodera et al. (2016) Onodera, M., Carollo, C. M., Lilly, S., et al. 2016, ApJ, 822, 42, doi: 10.3847/0004-637X/822/1/42
  * Ormerod et al. (2024) Ormerod, K., Conselice, C. J., Adams, N. J., et al. 2024, MNRAS, 527, 6110, doi: 10.1093/mnras/stad3597
  * Osterbrock (1989) Osterbrock, D. E. 1989, Astrophysics of gaseous nebulae and active galactic nuclei
  * Pentericci et al. (2018) Pentericci, L., McLure, R. J., Garilli, B., et al. 2018, A&A, 616, A174, doi: 10.1051/0004-6361/201833047
  * Reddy et al. (2023) Reddy, N. A., Topping, M. W., Sanders, R. L., Shapley, A. E., & Brammer, G. 2023, ApJ, 952, 167, doi: 10.3847/1538-4357/acd754
  * Rigby et al. (2023) Rigby, J., Perrin, M., McElwain, M., et al. 2023, PASP, 135, 048001, doi: 10.1088/1538-3873/acb293
  * Sanders et al. (2024) Sanders, R. L., Shapley, A. E., Topping, M. W., Reddy, N. A., & Brammer, G. B. 2024, ApJ, 962, 24, doi: 10.3847/1538-4357/ad15fc
  * Sanders et al. (2016) Sanders, R. L., Shapley, A. E., Kriek, M., et al. 2016, ApJ, 816, 23, doi: 10.3847/0004-637X/816/1/23
  * Shibuya et al. (2015) Shibuya, T., Ouchi, M., & Harikane, Y. 2015, ApJS, 219, 15, doi: 10.1088/0067-0049/219/2/15
  * Shimakawa et al. (2015) Shimakawa, R., Kodama, T., Steidel, C. C., et al. 2015, MNRAS, 451, 1284, doi: 10.1093/mnras/stv915
  * Shirazi et al. (2014) Shirazi, M., Brinchmann, J., & Rahmati, A. 2014, ApJ, 787, 120, doi: 10.1088/0004-637X/787/2/120
  * Silverman et al. (2015) Silverman, J. D., Kashino, D., Sanders, D., et al. 2015, ApJS, 220, 12, doi: 10.1088/0067-0049/220/1/12
  * Steidel et al. (2014) Steidel, C. C., Rudie, G. C., Strom, A. L., et al. 2014, ApJ, 795, 165, doi: 10.1088/0004-637X/795/2/165
  * Stiavelli et al. (2023) Stiavelli, M., Morishita, T., Chiaberge, M., et al. 2023, ApJ, 957, L18, doi: 10.3847/2041-8213/ad0159
  * Stott et al. (2016) Stott, J. P., Swinbank, A. M., Johnson, H. L., et al. 2016, MNRAS, 457, 1888, doi: 10.1093/mnras/stw129
  * Sugimura et al. (2024) Sugimura, K., Ricotti, M., Park, J., Garcia, F. A. B., & Yajima, H. 2024, arXiv e-prints, arXiv:2403.04824, doi: 10.48550/arXiv.2403.04824
  * Sutherland & Dopita (1993) Sutherland, R. S., & Dopita, M. A. 1993, ApJS, 88, 253, doi: 10.1086/191823
  * Swinbank et al. (2019) Swinbank, A. M., Harrison, C. M., Tiley, A. L., et al. 2019, MNRAS, 487, 381, doi: 10.1093/mnras/stz1275
  * Tacchella et al. (2023) Tacchella, S., Eisenstein, D. J., Hainline, K., et al. 2023, ApJ, 952, 74, doi: 10.3847/1538-4357/acdbc6
  * Tang et al. (2023) Tang, M., Stark, D. P., Chen, Z., et al. 2023, MNRAS, 526, 1657, doi: 10.1093/mnras/stad2763
  * The Astropy Collaboration et al. (2018) The Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ, 156, 123, doi: 10.3847/1538-3881/aabc4f
  * Treu et al. (2022) Treu, T., Roberts-Borsani, G., Bradac, M., et al. 2022, ApJ, 935, 110, doi: 10.3847/1538-4357/ac8158
  * Urrutia et al. (2019) Urrutia, T., Wisotzki, L., Kerutt, J., et al. 2019, A&A, 624, A141, doi: 10.1051/0004-6361/201834656
  * Virtanen et al. (2020) Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nature Methods, 17, 261, doi: 10.1038/s41592-019-0686-2
  * Williams et al. (2023) Williams, H., Kelly, P. L., Chen, W., et al. 2023, Science, 380, 416, doi: 10.1126/science.adf5307
  * Wisnioski et al. (2015) Wisnioski, E., Förster Schreiber, N. M., Wuyts, S., et al. 2015, ApJ, 799, 209, doi: 10.1088/0004-637X/799/2/209
  * Yuan et al. (2012) Yuan, T.-T., Kewley, L. J., & Richard, J. 2012, The Astrophysical Journal, 763, 9, doi: 10.1088/0004-637X/763/1/9