#  JWST Measurements of Neutral Hydrogen Fractions and Ionized Bubble Sizes at
$z=7-12$
Obtained with Ly$\alpha$ Damping Wing Absorptions in 26 Bright Continuum
Galaxies

Hiroya Umeda Institute for Cosmic Ray Research, The University of Tokyo, 5-1-5
Kashiwanoha, Kashiwa, Chiba 277-8582, Japan  Department of Physics, Graduate
School of Science, The University of Tokyo, 7-3-1 Hongo, Bunkyo, Tokyo
113-0033, Japan<EMAIL_ADDRESS>Masami Ouchi National Astronomical
Observatory of Japan, 2-21-1 Osawa, Mitaka, Tokyo 181-8588, Japan  Institute
for Cosmic Ray Research, The University of Tokyo, 5-1-5 Kashiwanoha, Kashiwa,
Chiba 277-8582, Japan  Kavli Institute for the Physics and Mathematics of the
Universe (WPI), University of Tokyo, Kashiwa, Chiba 277-8583, Japan  Kimihiko
Nakajima National Astronomical Observatory of Japan, 2-21-1 Osawa, Mitaka,
Tokyo 181-8588, Japan  Yuichi Harikane Institute for Cosmic Ray Research, The
University of Tokyo, 5-1-5 Kashiwanoha, Kashiwa, Chiba 277-8582, Japan
Department of Physics and Astronomy, University College London, Gower Street,
London WC1E 6BT, UK  Yoshiaki Ono Institute for Cosmic Ray Research, The
University of Tokyo, 5-1-5 Kashiwanoha, Kashiwa, Chiba 277-8582, Japan  Yi Xu
Institute for Cosmic Ray Research, The University of Tokyo, 5-1-5 Kashiwanoha,
Kashiwa, Chiba 277-8582, Japan  Department of Astronomy, Graduate School of
Science, The University of Tokyo, 7-3-1 Hongo, Bunkyo, Tokyo 113-0033, Japan
Yuki Isobe Institute for Cosmic Ray Research, The University of Tokyo, 5-1-5
Kashiwanoha, Kashiwa, Chiba 277-8582, Japan  Department of Physics, Graduate
School of Science, The University of Tokyo, 7-3-1 Hongo, Bunkyo, Tokyo
113-0033, Japan  Yechi Zhang Institute for Cosmic Ray Research, The University
of Tokyo, 5-1-5 Kashiwanoha, Kashiwa, Chiba 277-8582, Japan  Department of
Astronomy, Graduate School of Science, The University of Tokyo, 7-3-1 Hongo,
Bunkyo, Tokyo 113-0033, Japan

###### Abstract

We present volume-averaged neutral hydrogen fractions $x_{\rm H\,\textsc{i}}$
and ionized bubble radii $R_{\rm b}$ measured with Ly$\alpha$ damping wing
absorption of galaxies at the epoch of reionization. We combine JWST/NIRSpec
spectra taken by CEERS, GO-1433, DDT-2750, and JADES programs, and obtain a
sample containing 26 bright UV-continuum ($M_{\rm UV}<-18.5~{}{\rm mag}$)
galaxies at $7<z<12$. We construct 4 composite spectra binned by redshift, and
find the clear evolution of softening break towards high redshift at the rest-
frame $1216$ Å, suggesting the increase of Ly$\alpha$ damping wing absorption.
We estimate Ly$\alpha$ damping wing absorption in the galaxy spectra with
realistic templates including Ly$\alpha$ emission and circum-galactic medium
absorptions. Assuming the standard inside-out reionization picture having an
ionized bubble with radius $R_{b}$ around a galaxy embedded in the
intergalactic medium with $x_{\rm H\,\textsc{i}}$, we obtain $x_{\rm
H\,\textsc{i}}$ ($R_{\rm b}$) values generally increasing (decreasing) from
$x_{\rm H\,\textsc{i}}={0.54}^{+0.13}_{-0.54}$ to ${0.94}^{+0.06}_{-0.41}$
($\log R_{\rm b}={1.89}^{+0.49}_{-1.54}$ to ${-0.72}^{+1.57}_{-0.28}$ comoving
Mpc) at redshift $7.12^{+0.06}_{-0.08}$ to $10.28^{+1.12}_{-1.40}$. The
redshift evolution of $x_{\rm H\,\textsc{i}}$ indicates a moderately late
reionization history consistent with the one previously suggested from the
electron scattering of cosmic microwave background and the evolution of UV
luminosity function with an escape fraction $f_{\rm
esc}{\color[rgb]{0,0,0}\sim 0.2}$. Our ${R_{\rm b}}$ measurements suggest that
bubble sizes could be up to a few dex larger than the cosmic average values
estimated by analytic calculations for a given $x_{\rm H\,\textsc{i}}$, while
our $R_{\rm b}$ measurements are roughly comparable with the values for merged
ionized bubbles around bright galaxies predicted by recent numerical
simulations.

Galaxy evolution (594), Galaxy formation (595), High-redshift galaxies (734),
Reionization(1383)

##  1 Introduction

Cosmic reionization scenarios have been investigated by estimating the neutral
hydrogen fraction $x_{\rm H\,\textsc{i}}$ of the universe at different
redshifts. Most powerful way to probe the neutral hydrogen at the epoch of
reionization is to look at Ly$\alpha$ damping wing absorption. Numerous works
have investigated $x_{\rm H\,\textsc{i}}$ from Ly$\alpha$ damping wing
absorption using the brightest sources such as quasars (QSOs) and gamma-ray
burst (GRBs) at $z\sim 6-8$ (e.g., Mortlock et al., 2011; Bañados et al.,
2018; Schroeder et al., 2013; Davies et al., 2018; Greig et al., 2019; Wang et
al., 2020; Totani et al., 2006, 2014). With the bright rest-frame UV-
continuum, we can detect Ly$\alpha$ damping wing absorption by the
intergalactic medium (IGM). $x_{\rm H\,\textsc{i}}$ value at different
redshift have been estimated by fitting the damping wing absorption profile
for different $x_{\rm H\,\textsc{i}}$ values. However, the number density of
QSOs and GRBs are expected to be very small at the higher redshift, inhibiting
us from accessing the ionization state of IGM at $z\gtrsim 7$ (e.g., Wang et
al., 2019; Matsuoka et al., 2023).

To overcome the problem of lacking the bright objects at the high redshift, we
investigate the fainter but much more numerous galaxies as alternative target
to infer $x_{\rm H\,\textsc{i}}$ from Ly$\alpha$ damping wing absorptions. To
directly observe the Ly$\alpha$ damping wing absorption in the spectrum of
galaxies at $z>7$, we need to detect their rest-frame UV continuum.
Spectroscopically detecting the rest-frame UV continuum at such a high
redshift is very challenging because it requires deep spectroscopic data in
near infrared (NIR). However, recent observations by the James Webb Space
Telescope (JWST; Gardner et al., 2023; Rigby et al., 2023; Rieke et al., 2023;
Jakobsen et al., 2022; Ferruit et al., 2022; Böker et al., 2023) have shown
that it is capable of deep NIR spectroscopic observations. As presented in
recent works, low-resolution (i.e., $R\sim 100$) NIRSpec/PRISM spectroscopic
data is especially suited for detecting faint rest-frame UV continuum (e.g.,
Curtis-Lake et al., 2023; Nakajima et al., 2023; Harikane et al., 2023;
Arrabal Haro et al., 2023a, b; Hsiao et al., 2023; Williams et al., 2023;
Roberts-Borsani et al., 2022).

Curtis-Lake et al. (2023) are the first to demonstrate the power of the PRISM
data to estimate Ly$\alpha$ damping wing absorption in the UV continua of high
redshift galaxies. Curtis-Lake et al. (2023) have demonstrated the potential
of the NIRSpec PRISM data to estimate the value of $x_{\rm H\,\textsc{i}}$.
However, their constraint on $x_{\rm H\,\textsc{i}}$ is not very strict
because the redshift for the object is not constrained due to the lack of
emission line detections. Recently, Hsiao et al. (2023) have incorporated
Ly$\alpha$ damping wing absorption in the spectral fitting of
spectroscopically confirmed $z=10.17$ lensed galaxy and constrained $x_{\rm
H\,\textsc{i}}>0.9$ with a one sigma significance. In our work, we also take
advantage of the multiple PRISM spectroscopy data with the redshift
confirmation obtained by multiple surveys to systematically infer the $x_{\rm
H\,\textsc{i}}$ values at various redshift above $z=7$ to capture the redshift
evolution of the damping wing profile and corresponding $x_{\rm
H\,\textsc{i}}$ values.

This paper is constructed as follows. In Section 2, we describe our galaxy
sample and the reduction of corresponding spectroscopic data. We construct
composite spectra binned by redshift and fit for the Ly$\alpha$ damping wing
absorption seen in the composite spectra in Section 3. We discuss the
implications of the $x_{\rm H\,\textsc{i}}$ estimates on understanding cosmic
reionization history and future prospect in Section 4. We summarize this paper
in Section 5. In this paper, we use the Planck cosmological parameter sets of
the TT + TE + EE+ lowE + BAO + lensing result (Planck Collaboration et al.,
2020a): $h=0.6766$, $\Omega_{m}=0.3103$, $\Omega_{\Lambda}=0.6897$,
$\Omega_{b}h^{2}=0.02234$, and $Y_{p}=0.248$. All magnitudes are in the AB
system (Oke & Gunn, 1983).

##  2 Data and Sample

We use spectroscopic datasets obtained in multiple observation programs; the
Cosmic Evolution Early Release Science (CEERS; ERS-1345. PI. S. Finkelstein;
Finkelstein et al., 2023), GO-1433 (PI. D. Coe; Hsiao et al., 2023), DDT-2750
(PI. P. Arrabal Halo; Arrabal Haro et al., 2023b), and the JWST Advanced Deep
Extragalactic Survey (JADES; GTO-1210; PI. N. Lützgendorf; Bunker et al.,
2023). We only used the PRISM data ($R\sim 100$) that covers $0.6-5.3~{}\mu$m.
The total exposure time for the PRISM spectra from CEERS, GO-1433, DDT-2750,
and JADES data are 0.86, 3.7, 5.1, and 9.3-27.7 hours, respectively. As
presented in Nakajima et al. (2023) and Harikane et al. (2023), all the PRISM
data from CEERS, GO-1433, and DDT-2750 we use are reduced with the JWST
pipeline version 1.8.5 with the Calibration Reference Data System (CRDS)
context file of jwst1028.pmap or jwst1027.pmap with some improvements on the
flux calibration, noise estimate, and the composition. For the JADES data, we
use the publicly available reduced data from JADES group
111https://archive.stsci.edu/hlsp/jades. Please see Bunker et al. (2023) for
the detailed reduction processes of JADES data. We then select the galaxies
with systemic redshifts above $z>7$ confirmed with emission lines. For the
$z<9$ galaxies in CEERS, we adopt the spectroscopic redshift reported by
Nakajima et al. (2023) in which the authors determine spectroscopic redshift
via fitting [O iii] and H$\beta$ triplet simultaneously. For the galaxies at
$z>9$ obtained from observational program beside JADES, we adopt the
spectroscopic redshift reported by Harikane et al. (2023). In Harikane et al.
(2023), the authors derive the spectroscopic redshift for $z>11$ object via [O
ii] (and [Ne iii]) for Maisie’s galaxy (CEERS2_588). For MACS0647_JD, Harikane
et al. (2023) use H$\gamma$ to determine spectroscopic redshift. For JADES
objects, we adopt redshift reported in the public catalog released by the
JADES team. The determination of the spectroscopic redshift listed in the
public JADES is described in Bunker et al. (2023). In Bunker et al. (2023),
the authors use emission lines detected with a signal to noise ratio of
$S/N>5$ to determine the spectroscopic redshift. Although the uncertainty of
the spectroscopic redshift is not given in the literature, Jones et al. (2023)
report that the spectroscopic redshift determined from PRISM data has a median
offset of 0.00388 with a standard deviation of 0.00628. In Section 4.5, we
discuss on the effect of the uncertainty in the redshift confirmation with low
resolution PRIMS data to our analysis. Please see Nakajima et al. (2023),
Harikane et al. (2023), and Bunker et al. (2023) for further details on the
data reduction and spectroscopic redshift determination. We further narrow the
sample by selecting galaxies with data available in wavelength corresponding
to rest-frame $1200-1250$ Å to investigate Ly$\alpha$ damping wing features.
We obtain a total of 34 galaxies as our galaxy sample candidates. Our sample
candidates consist of 34 galaxies from $z=$7.045 to 11.40, with the redshift
distribution shown in Figure 1. In Table 1, we summarize properties of the
sample candidate galaxies. We apply additional sample selection by setting an
absolute UV magnitude cut of $M_{\rm UV}<-18.5$ to omit faint ($M_{\rm
UV}>-18.5$) galaxies distributed mainly at $z<8$ in our sample candidates. We
also discard three of the object which Nakajima et al. (2023) only report
5$\sigma$ upper limit value for UV magnitude for Hubble Space Telescope (HST)
observations given that these objects are not detected by HST and lack
JWST/NIRCam data. After this magnitude cut selection, we have a final galaxy
sample consisted by 26 galaxies from $7.045<z<11.40$. Out of all 26 galaxies,
spectroscopic data of 18, 1, 2, and 5 galaxies are obtained in CEERS, GO-1433,
DDT-2750, and JADES programs, respectively. We present the complete PRISM
spectra of galaxies in our final sample in the left panels of Figure 2 in the
order of the highest to lowest redshift galaxies. The black line and grey
region in each panel represents the observed flux and the error of each galaxy
spectra normalized by an arbitrary unit (a.u.), respectively. We also show the
wavelengths of Ly$\alpha$ determined from the systemic redshifts by the red
vertical dashed lines. Next to the spectra of each galaxy, we present the
emission line detected in the PRISM spectra with the highest $S/N$ for each
galaxy in the right hand panels. The systemic redshift values and the name of
the emission line shown for each galaxy are written in between the panels.

The spectra of the galaxy in our final sample show various features around
rest-frame Ly$\alpha$ wavelength. We show two contrasting example in Figure 3.
The 2D and 1D spectra of CEERS2_588 (00021842) with systemic redshift at 11.04
(7.980) is shown in the left (right) panel. The red dotted lines represent the
systemic redshift determined by emission line. The spectra of CEERS2_588 shows
redshifted, softened Ly$\alpha$ break relative to systemic redshift, while we
see blueshifted flux with 4.2$\sigma$ significance beyond the rest-frame 1216
Å for 00021842. In Section 3.4, we physically interpret the diverse spectral
features around the rest-frame 1216 Å via spectral fitting.

Figure 1: Redshift and UV-magnitude distribution of our galaxy sample. In the
bottom panel, circles represent the redshift and UV-magnitude values for the
galaxy observed by CEERS, GO-1433, DDT-2750, and JADES. The circle indicate
the best-fit values while the upper circle with an upper arrow corresponds to
one sigma upper lower-limit for the UV-magnitudes. The filled circle indicates
the galaxy used in our analysis, while the open circle indicates the galaxies
omitted from our sample after a selection based on the UV-magnitudes. The top
(right) panel shows the number counts at each redshift ($M_{\rm UV}$) bins for
our galaxy sample. The fraction of bins represented by black represents the
final sample after the absolution UV-magnitude selection. Figure 2: Full
spectra of galaxy samples. (Left): The x-(y-)axis corresponds to the observed
wavelength (arbitrary normalized flux density). The black lines and grey
regions represent the fluxes of spectra and their uncertainties, respectively.
The red dashed line represents the observed Ly$\alpha$ wavelength for each
object. (Right) The highest $S/N$ emission lines detected from each spectra in
our sample. The type of emission line shown and the name of the corresponding
galaxy is shown between the right and left panels. The x-axis represents the
difference between the rest-frame wavelength $\lambda_{\rm rest.}$ and the
central wavelength $\lambda_{\rm cen.}$ of the optical emission line. Table 1:
Sample Properties

Name | R.A. | Decl. | $z_{\rm spec}$ | $M_{\rm UV}$ | Ref.
---|---|---|---|---|---
| (deg) | (deg) |  | (mag) |
(1) | (2) | (3) | (4) | (5) | (6)
00020961 | 53.13423 | -27.76891 | 7.045 | -18.840 | l,m
CEERS_00542 | 214.831624 | 52.831505 | 7.064 | -19.89 | h
CEERS_00044 | 215.001115 | 53.011269 | 7.104 | -19.36 | h
CEERS_00534 | 214.859117 | 52.85364 | 7.115 | -20.1 | h
CEERS_00829 | 214.861594 | 52.876159 | 7.167 | -19.57 | h
CEERS_80374 | 214.898074 | 52.824895 | 7.174 | -18.09 | h
CEERS_00439 | 214.825364 | 52.863065 | 7.179 | -19.25 | h
CEERS_00498 | 214.813045 | 52.834249 | 7.179 | -20.18 | h
10013905 | 53.11833 | -27.76901 | 7.197 | -18.800 | l,m
CEERS_01038 | 215.039697 | 52.9015971 | 7.201 | -19.2 | h
00008079 | 53.15283 | -27.80194 | 7.260 | -17.950 | l,m
10013682 | 53.16746 | -27.77201 | 7.275 | -17.000 | l,m
CEERS_01163 | 214.990468 | 52.9719902 | 7.455 | $>-20.78$ | h
CEERS_80432 | 214.812056 | 52.746747 | 7.477 | -20.05 | h
CEERS_80372 | 214.927798 | 52.850003 | 7.483 | -19.27 | h
CEERS_80239 | 214.896054 | 52.869853 | 7.487 | -17.92 | h
CEERS_80445 | 214.843115 | 52.747886 | 7.509 | -19.66 | h
CEERS_00689 | 214.999053 | 52.9419767 | 7.552 | $>-20.93$ | h
CEERS_00686 | 215.150862 | 52.9895618 | 7.752 | $>-20.42$ | h
CEERS_01023 | 215.188413 | 53.0336473 | 7.781 | -21.06 | h
CEERS_01027 | 214.882994 | 52.8404159 | 7.834 | -20.73 | h
00021842 | 53.15682 | -27.76716 | 7.980 | -18.590 | l,m
CEERS_00003 | 215.005189 | 52.99658 | 8.005 | -18.64 | h
CEERS_01149 | 215.089714 | 52.9661828 | 8.184 | -20.83 | h
00008013 | 53.16446 | -27.80218 | 8.473 | -17.540 | l,m
EGS_z910_44164 | 215.218762 | 53.0698619 | 8.615 | -21.63 | f,j
|  |  |  |  | h,d
CEERS_90671 | 214.961276 | 52.842364 | 8.637 | -18.7 | b,d
CEERS_1019 | 215.035391 | 52.8906618 | 8.683 | -22.07 | h,j,k
|  |  |  |  | i,g,d
CEERS1_6059 | 215.011706 | 52.988303 | 8.873 | -20.75 | c,d,h
10058975 | 53.11243 | -27.77461 | 9.433 | -20.355 | l,m
00006438 | 53.16735 | -27.80750 | 9.689 | -19.270 | l,m
MACS0647-JD† | 101.982213 | 70.2432667 | 10.17 | -20.3 | d,e
CEERS2_588‡ | 214.906625 | 52.9115278 | 11.04 | -20.4 | d
Maisie’s Galaxy‡ | 214.943167 | 52.9424417 | 11.40 | -20.1 | a,d
(CR2-z12-1) |  |  |  |  |

Note. — (1): Name. (2): Right ascension of the object’s coordinate. (3):
Declination of the object’s coordinate. (4): Spectroscopic redshift. (5):
Absolute UV magnitude. (6): References for spectroscopic redshifts and UV
magnitude:a.) Arrabal Haro et al. (2023a), b.) Arrabal Haro et al. (2023b),
c.) Fujimoto et al. (2023a), d.) Harikane et al. (2023), e.) Hsiao et al.
(2023), f.)Larson et al. (2022), g.) Larson et al. (2023), h.) Nakajima et al.
(2023), i.) Sanders et al. (2023), j.) Tang et al. (2023), k.) Zitrin et al.
(2015), l.) Bunker et al. (2023), m.) Jones et al. (2023)
All spectroscopic data are from CEERS program except for the spectroscopic
data of the object name with † and ‡ are from GO-1433 and DDT-2750 programs,
respectively.

Figure 3: 2D (top) and 1D (bottom) spectra of CEERS2_588 (left) and 00021842
(right). The black line and grey region represent the observed flux and its
uncertainty, respectively. The red dashed line represents the wavelength of
Ly$\alpha$ at the systemic redshift determined by non-Ly$\alpha$ emission
lines.

##  3 Methods

###  3.1 Ly${\alpha}$ Damping Wing

To calculate the optical depth associated with Ly$\alpha$ damping wing
absorption, we assume that galaxy resides at the center of the fully ionized
bubble with radius $R_{b}$ and the surrounding IGM with a constant $x_{\rm
HI}$ until the end of EoR, $z_{n}$. We fix the redshift corresponding to the
end of reionization at $z_{n}=5$ in our calculations. Note that changing
$z_{n}$ by $\pm 1$ have a negligible effect on IGM Ly$\alpha$ damping wing
absorption profile.

To calculate the optical depth of Ly$\alpha$ damping wing at each observed
wavelength $\lambda_{\rm obs}$, $\tau(\lambda_{\rm obs})$, we adopt the
equation 2 of Totani et al. (2006) that is based on the formulation of
Miralda-Escudé (1998):

$\displaystyle\tau(\lambda_{\rm obs})$ $\displaystyle=\frac{x_{\rm
H\,\textsc{i}}\Lambda_{\alpha}\lambda_{\alpha}\tau_{0}(z_{s})}{4{\color[rgb]{0,0,0}\pi^{2}}c}\left(\frac{1+z_{\rm
obs}}{1+z_{s}}\right)^{3/2}$ (1)
$\displaystyle\times\left[I\left(\frac{1+z_{b}}{1+z_{\rm
obs}}\right)-I\left(\frac{1+z_{n}}{1+z_{\rm obs}}\right)\right],$

where $\Lambda_{\alpha}$, $\lambda_{\alpha}$, $c$, and $\tau_{0}$, are the
decay constant for the Ly$\alpha$ resonance, the rest-frame Ly$\alpha$
wavelength, the speed of light, and the Gunn-Peterson optical depth,
respectively. $z_{b}$ is defined as the redshift corresponding to the edge of
the ionized bubble with radius $R_{b}$, whereas $z_{\rm obs}$ is a value
defined as $\lambda_{\rm obs}/\lambda_{\alpha}-1$. $I(x)$ is an integration
defined as follows:

$\displaystyle I(x)$
$\displaystyle=\frac{x^{9/2}}{1-x}+\frac{9}{7}x^{7/2}+\frac{9}{5}x^{5/2}+3x^{3/2}$
(2) $\displaystyle+9x^{1/2}-\frac{9}{2}\ln{\frac{1+x^{1/2}}{1-x^{1/2}}}.$

We calculate $\tau_{0}$ using the following formula also given by Miralda-
Escudé (1998):

$\tau_{0}(z_{s})=\frac{3\lambda_{\alpha}^{3}\Lambda_{\alpha}n_{0}}{8\pi
H(z_{s})},$ (3)

where $n_{0}$ and $H(z)$ are the hydrogen number density and the Hubble
parameter as the function of redshift. We calculate $n_{0}$ using the critical
density of the universe, baryon density parameter, and primordial helium
abundance. We show how the absorption profile of the $z=8$ changes with
different $x_{\rm HI}$ and $R_{b}$ values in Figure 4. As shown in Figure 4,
increasing $x_{\rm HI}$ value extends the absorption feature to a redder side
of 1216 Å while increasing $R_{b}$ value shifts the absorption profile to the
bluer wavelength. Notably, fluxes on the bluer side of the 1216 Å can be
transmitted if the bubble size is sufficiently large (i.e., $R_{b}>10$
comoving Mpc; cMpc).

Figure 4: Example of Ly$\alpha$ damping wing profile for $z=8$. The top and
bottom axes are the rest frame and observed wavelengths, respectively. The
y-axis corresponds to the transmission fraction of the flux at each
wavelength. The solid (dotted) black line corresponds to the absorption
profile with $x_{\rm H\,\textsc{i}}=1.0$ (0.1) and $R_{b}=0$ cMpc. The red
dotted line corresponds to the absorption profile with $x_{\rm
H\,\textsc{i}}=1.0$ and $R_{b}=10$ cMpc. The gray vertical line corresponds to
the rest-frame Ly$\alpha$ wavelength.

###  3.2 Stacking Spectra

We construct composite spectra by stacking individual spectra to maximize the
signal to noise ratio. We divide the galaxy sample into four subsamples by
binning according to the redshifts so that each subsample has at least six
galaxies. We label the subsamples as #1, #2, #3, and #4 bins in the order of
the lower to higher redshift. The redshift range for #1, #2, #3, and #4 bins
are 7.045-7.179, 7.197-7.781, 7.834-8.683, and 8.873-11.40, respectively. We
define representative redshift values $\langle z\rangle$ of each redshift bin
by the average redshift value among the galaxies in each subsample.

For stacking individual spectra of each subsample, we have shifted the galaxy
spectra into rest-frame wavelength based on the systemic redshift and
resampled to a common wavelength pixels with a pixel size of 25 Å, a value
considerably similar to that of the resolution of PRISM at $<$2 $\mu$m ($R\sim
30$), using linear interpolation. We have taken error weighted mean flux
values at wavelength pixel to calculate the stacked spectra. We estimate the
uncertainty of the stacked spectra via 10000 iterations of bootstrap
resampling. During each resampling, we fluctuate the observed flux at each
pixel with Gaussian noise accordingly to the original error spectrum to
consider the measurement error as well. We present the stacked spectra from
all subsamples in Figure 5. Fluxes of the composite spectra are normalized at
the rest frame 1450 Å. As shown in Figure 5, we can clearly see the softening
break feature at the rest-frame 1216 Å towards the high redshift, suggesting
increasing Ly$\alpha$ damping wing absorptions due to increasing $x_{\rm
H\,\textsc{i}}$ with redshifts.

Figure 5: Stacked spectra of our galaxy sample binned into four redshift bins.
The dark blue, pale blue, pink, and red solid lines (shades) represent the
stacked spectra (uncertainties) for the redshift bins at $\langle
z\rangle=$7.12, 7.44, 8.28, and 10.28, respectively. The wavelength is shifted
to the rest frame and the flux is normalized at 1450 Å at which the IGM
absorption have negligible impact on the spectrum. The black vertical dotted
line represent the position of the rest-frame 1216 Å. The horizontal solid
line indicates the position of the flux value at zero. Table 2: Subsamples
Properties

ID | $N$ | $\langle z\rangle$ | $z_{\rm lower}$ | $z_{\rm upper}$
---|---|---|---|---
(1) | (2) | (3) | (4) | (5)
#1 | 7 | 7.12 | 7.045 | 7.179
#2 | 6 | 7.44 | 7.197 | 7.781
#3 | 7 | 8.28 | 7.834 | 8.683
#4 | 6 | 10.28 | 8.873 | 11.40

Note. — (1): ID for a subsample. (2): Size of a subsample. (3): Representative
(i.e., the average) redshift of galaxies in a subsample. (4): Lowest redshift
of galaxies in a subsample. (5): Highest redshift of galaxies in a subsample.

###  3.3 Spectral Fitting

To further confirm the redshift evolution of Ly$\alpha$ damping wing
absorption, we measure the IGM Ly$\alpha$ damping wing absorption. At the same
time, we also consider the effect on the shape of Lyman break from the
intrinsic galaxy’s property such as stellar mass. Because physical properties
of galaxies differ from one to another, we fit an individual observed galaxy
spectrum to modeled galaxy spectra with various physical properties and
Ly$\alpha$ absorptions. For the fitting, we use the observed spectra in the
wavelength range above the rest-frame 1000 Å. We included data at shorter
wavelength value than the rest-frame Ly$\alpha$ (i.e., 1216 Å) to consider the
effect of higher flux transmission due to large surrounding ionized bubble. To
compare the observed spectrum to model spectrum, we first generate the model
spectra incorporating nebular emission, Ly$\alpha$ emission and absorption.
After we implement all the physical ingredients to the model spectra, we
smooth the model spectra with line spread function of PRISM provided by STScI
to account for the instrumental broadening. 222https://jwst-
docs.stsci.edu/jwst-near-infrared-spectrograph/nirspec-
instrumentation/nirspec-dispersers-and-filters

We use a SED fitting code based on prospector (Johnson et al., 2021) to
generate model galaxy spectra with various stellar population and nebular
properties. We use the Binary Population and Spectral Synthesis (BPASS;
Eldridge et al., 2017) model as an isochrone library. We allow different non-
parametric star formation history (SFH), stellar mass ($M_{\star}$),
ionization parameter ($\log U$), and dust attenuation. We adopt the star
formation history in the same manner as Harikane et al. (2023). We also allow
different stellar mass in the range of $6<M_{\star}<11$, which fully covers
the mass range of $z>7$ galaxies estimated in Nakajima et al. (2023) and
Harikane et al. (2023). For the dust attenuation, we apply the attenuation law
by Calzetti et al. (2000) and use an optical depth at 5500 Å as a free
parameter. We also include nebular emission continuum and lines, but we treat
Ly$\alpha$ emission line separately as we describe later in this section. We
fix stellar/nebular metallicity to 10% solar metallicity as suggested for JWST
high-z galaxy by recent work of Nakajima et al. (2023).

For Ly$\alpha$ absorption, we consider Ly$\alpha$ absorption by both IGM and
circum-galactic medium (CGM). We adopt the formulation given by equation 3.1
for the IGM absorption with $x_{\rm H\,\textsc{i}}$ and $\log R_{b}$ as free
parameter. For the H i absorption by CGM, we consider the Voigt profile with
different neutral hydrogen column density ($\log N_{\rm H\,\textsc{i}}$),
velocity dispersion ($\sigma_{v}$), and velocity offset relative to the
systemic redshift ($\Delta_{v,\rm DLA}$). Moreover, we let the equivalent
width of the intrinsic Ly$\alpha$ emission line ($\log EW_{\alpha,0}$) freed
in the range of $-3<\log EW_{\alpha,0}/{\rm\AA}<3$. In this way, we can
consider both extreme and negligible (i.e., low escape fraction) contribution
of Ly$\alpha$ emission to the shape of Lyman break. We discuss about the
validity of Ly$\alpha$ emission line strengths in Section 4.2. We fix the
velocity offset and dispersion of the Ly$\alpha$ emission line to +200 km/s
and 250 km/s, respectively.

We run Markov Chain Monte Carlo (MCMC) method to obtain the posterior
probability distribution functions (PDFs) of free parameters. We use a package
created based on emcee (Foreman-Mackey et al., 2013) called ptemcee (Vousden
et al., 2021) for MCMC. ptemcee incorporates parallel tempering technique that
realizes faster sampling convergence in case of multi-modal posterior
probability distribution functions. We run four temperature parallel tempering
MCMC sampler of 80 walkers for each temperature with 5000 steps. We discard
first 1000 steps of each walker to diminish the dependence of posterior PDFs
to initial parameter distributions.

We infer the best-fit parameters from the posterior PDFs. Because a median
value does not well represent the skewed posterior distribution, we define the
best-fit parameter by the mode (i.e., a peak of the posterior distribution)
and uncertainty by 68th percentile highest posterior density interval (HPDI;
i.e., a narrowest interval containing 68%) of the posterior PDFs. To obtain
mode of 1D marginalized posterior PDF for each parameter, we smooth the
distribution by Kernel Density Estimation smoothing using scipy package
(Virtanen et al., 2020) and choose the mode within an HPDI of 1D posterior
distribution. We also add up 1D posterior probability distribution within each
redshift bin (hereafter, stacked posterior PDF) and obtain the best-fit value
in the same manner for each redshift bin.

###  3.4 Fitting Results

In Figure 6, we show the stacked posterior PDF of the parameter related to the
Ly$\alpha$ damping wing absorption for each redshift bin. The three contours
with different darkness correspond to the 1, 2, and 3$\sigma$ from the darkest
to the faintest order, respectively. The histograms represent the marginalized
probability density function of the $x_{\rm H\,\textsc{i}}$ and $R_{b}$. The
three dashed lines in each histogram represent the mode and lower/upper
boundary of 68th percentile HPDI, respectively. We summarize the best-fit
template and the best-fit parameters in Table 3. For the results of the
stacked spectra in #4 bin, we obtain very high $x_{\rm H\,\textsc{i}}$ and
small $R_{b}$ values (i.e., $x_{\rm H\,\textsc{i}}\sim 1$,
$R_{b}~{}{\color[rgb]{0,0,0}<1}$ cMpc) as similar as the ones reported for a
$z=10.17$ galaxy in Hsiao et al. (2023). For the $z\sim 7$ galaxies, our
$x_{\rm H\,\textsc{i}}$ is not well constrained as one can see from the
marginalized PDF. The reason for tighter constraint on $\log R_{b}$ compared
to that for $x_{\rm H\,\textsc{i}}$ could be attributed to the strong
dependence of $\log R_{b}$ value to the position of observed Lyman break
relative to the systemic redshift.

In Figure 7, we also present examples of our spectrum fitting for individual
objects. From the top to bottom panels, we show the fitting results for
CEERS_1019, CEERS2_588, and 00021842, respectively. For CEERS_1019, strong IGM
absorption with $x_{\rm H\,\textsc{i}}={0.91}^{+0.09}_{-0.45}$ and an ionized
bubble size of $\log R_{b}/{\rm cMpc}={0.82}^{+0.61}_{-0.45}$ explains
moderately softened break feature around the rest-frame 1216 Å of the observed
spectrum. For CEERS2_588, in addition to strong IGM absorption (i.e., $x_{\rm
H\,\textsc{i}}\sim 1$) and small ionized region (i.e., $\log R_{b}/{\rm
cMpc}<0$), Ly$\alpha$ damping wing absorption by CGM with H i column density
of $\log N_{\rm H\,\textsc{i}}/{\rm cm^{-2}}>22$ surrounding the galaxy
explains strongly softened break feature around the rest-frame 1216 Å. In
contrast, the sharp and blue-shifted break feature at the rest-frame 1216 Å
for 00021842 can be explained with Ly$\alpha$ emission together with high
Ly$\alpha$ transmission environment realized with small IGM absorption (i.e.,
$x_{\rm H\,\textsc{i}}\sim 0$) and a large ionized bubble (i.e., $\log
R_{b}/{\rm cMpc}>1$).

Table 3: Best-fit Estimate for $x_{\rm H\,\textsc{i}}$ and $\log R_{b}$

ID | $\langle z\rangle$ | $x_{\rm H\,\textsc{i}}$ | $\log R_{b}$
---|---|---|---
|  |  | [cMpc]
(1) | (2) | (3) | (4)
#1 | $7.12^{+0.06}_{-0.08}$ | ${0.54}^{+0.13}_{-0.54}$ | ${1.89}^{+0.49}_{-1.54}$
#2 | $7.44^{+0.34}_{-0.24}$ | ${0.69}^{+0.30}_{-0.38}$ | ${2.31}^{+0.13}_{-2.11}$
#3 | $8.28^{+0.41}_{-0.44}$ | ${0.92}^{+0.08}_{-0.56}$ | ${0.70}^{+0.72}_{-1.25}$
#4 | $10.28^{+1.12}_{-1.40}$ | ${0.94}^{+0.06}_{-0.41}$ | ${-0.72}^{+1.57}_{-0.28}$

Note. — (1): ID for each subsample. (2): Mean redshift and the lower/upper
boundary of a subsample. (3): Mode value and upper/lower limit of HPDI of
$x_{\rm H\,\textsc{i}}$ for the stacked marginalized PDFs in each redshift
bin. each redshift bin. (4): Mode value and upper/lower limit of HPDI of $\log
R_{b}$ for the stacked marginalized PDFs in each redshift bin.

Figure 6: The stacked posterior probability distribution function (PDF) of
Ly$\alpha$ emission and Ly$\alpha$ absorptions by H i gas in the CGM/IGM for
each subsample. The results for #1 to #4 bin are shown in order of top left to
right bottom. Three shaded regions with different darkness in the bottom left
panel correspond to the 1, 2, and 3$\sigma$ regions of the sampled 2-D PDF
from the darkest to the faintest, respectively. The solid and dotted red lines
represent the mode value and HPDI of the stacked marginalized PDF for each
sampled parameters.

Figure 7: Fitted spectrum for the case of moderate, strong and weak Ly$\alpha$
damping wing absorption feature seen in CEERS_1019, CEERS2_588, and 000012842
from the top to bottom panels, respectively. The black line and shade
represent the observed spectrum and its error. The red (blue) line and shade
represent the model spectrum and its uncertainty after (before) Ly$\alpha$
absorptions applied and smoothed with PRISM resolution. The black dotted line
represent the rest-frame 1216 Å at the systemic redshift. The inset panel
shows zoom-in view of the spectrum around the rest-frame 1216 Å. The best-fit
values and uncertainties for $x_{\rm H\,\textsc{i}}$, $\log R_{b}$, and $\log
N_{\rm H\,\textsc{i}}$ for each objects are also written inside each panel.
The x-axis at the top and bottom of each panel represents the observed and
rest-frame wavelength, respectively. The constraints for $x_{\rm
H\,\textsc{i}}$, $\log R_{b}$, and $\log N_{\rm H\,\textsc{i}}$ are also
written inside each panel.

##  4 Discussion

###  4.1 Ly$\alpha$ Absorption by CGM

In Figure 8, we summarize the best-fit estimation of neutral hydrogen column
density of H i gas in CGM by redshift. In the same figure, we also present the
measured neutral hydrogen column density of $z\sim 0$ metal-poor galaxies
(Hernandez et al., 2020), $z\sim 0$ green pea galaxies (McKinney et al.,
2019), gamma ray bursts (Tanvir et al., 2019), and recent $z>8$ bright
galaxies (Heintz et al., 2023). Measured neutral hydrogen column density spans
up to the range of damped Ly$\alpha$ absorbers (DLA; $N_{\rm
H\,\textsc{i}}>2\times 10^{20}$). Moreover, neutral hydrogen column density
reaches up to $\gtrsim 10^{22}~{}{\rm cm}^{-2}$ for some objects, suggesting a
existence of extreme reservoir of neutral hydrogen gas around some of the
$z>7$ galaxies as suggested by Heintz et al. (2023). We compare the neutral
hydrogen column density measurement of overlapping two galaxy between our
sample and Heintz et al. (2023)’s and find that our measurements are
consistent within errorbars with Heintz et al. (2023)’s. This consistency
further solidifies the validity of our analysis and the existence of high-z
extreme DLA.

Figure 8: Neutral hydrogen column density of Ly$\alpha$ absorption by CGM at
different redshifts. The red circles (bars) represent $\log N_{\rm
H\,\textsc{i}}$ measurements (uncertainties). The blue upward triangles, gray
squares, and yellow downward triangles, and black diamonds represent the $\log
N_{\rm H\,\textsc{i}}$ measurements (uncertainties) of $z\sim 0$ green pea
galaxies (McKinney et al., 2019), $z\sim 0$ metal poor galaxies (Hernandez et
al., 2020), gamma ray bursts (Tanvir et al., 2019), and $z>9$ Lyman break
galaxies (Heintz et al., 2023), respectively.

###  4.2 Escape Fraction

Figure 9: Relation between $E(B-V)$ and Ly$\alpha$ escape fraction. The red
circle and bars (arrows) represent best-fit and errors (upper limits) for our
galaxy sample. The grey circle represent the measurement of $z\sim 0$ galaxies
presented in left panel of Figure 10 from Atek et al. (2014). The black dotted
line and yellow shades represent the best-fit relation and its errors between
$E(B-V)$ and $f_{\rm esc,Ly\alpha}$ for $z\sim 0$ galaxies as presented in
Atek et al. (2014).

As the existence of extremely high neutral hydrogen column density CGM is
suggested around some of the $z>7$ bright galaxies, we consider how much
Ly$\alpha$ emission can escape to the IGM (i.e., Ly$\alpha$ escape fraction
$f_{\rm esc,Ly\alpha}$). We estimate $f_{\rm esc,Ly\alpha}$ by comparing the
Ly$\alpha$ emission flux after the absorption by CGM and the intrinsic
Ly$\alpha$ emission flux. For the intrinsic Ly$\alpha$ emission flux, we use
cloudy prediction calculated by Byler et al. (2017) as Prospector adopts. We
note that we do not apply IGM absorption to calculate $f_{\rm esc,Ly\alpha}$.
We also estimate $E(B-V)$ using sampled posterior PDFs of optical depth at
5500 Å by dust attenuation. In Figure 9, we plot estimated $E(B-V)$ and
$f_{\rm esc,Ly\alpha}$ and compare with local relation discussed in Atek et
al. (2014). As shown in Figure 9, best-fit value for our $f_{\rm
esc,Ly\alpha}$ are below unity, which means that Ly$\alpha$ emission strength
applied in our analysis does not violate prediction from photoionization
models. Moreover, our estimates for $f_{\rm esc,Ly\alpha}$ are roughly
consistent with the local relation between $E(B-V)$ and $f_{\rm
esc,Ly\alpha}$, although some of our $f_{\rm esc,Ly\alpha}$ values are smaller
than the local trend. Because emission line flux is smoothed in low-resolution
PRISM spectra, $f_{\rm esc,Ly\alpha}$ cannot be well constrained for smaller
values. Moreover, galaxy sample used to derive the relation in Atek et al.
(2014) is partially selected based on the Ly$\alpha$ emission detection. In
contrast, we construct our galaxy sample based on the brightness of continuum
(i.e., UV magnitude), so it is natural to have our $f_{\rm esc,Ly\alpha}$ to
lie on or below the local trend of Atek et al. (2014). In any case, we have
confirmed that our treatment of Ly$\alpha$ emission is physical and consistent
with observational trend.

###  4.3 Cosmic Reionization History

Figure 10: Redshift evolution of $x_{\rm H\,\textsc{i}}$. The red circle and
bars represent the $x_{\rm H\,\textsc{i}}$ estimates and the uncertainties
from our work. The pink circle (square) and bars show the $x_{\rm
H\,\textsc{i}}$ estimate and uncertainty from $z=6-7$ Ly$\alpha$ luminosity
function (LAE clustering analysis) with data obtained with Subaru/Hyper
Suprime Cam (Umeda, 2023). The orange dashed, blue solid, and purple dotted
line and the shaded with corresponding colors are the model and its
uncertainty for the redshift evolution of $x_{\rm H\,\textsc{i}}$ predicted by
Naidu et al. (2020), Finkelstein et al. (2019), and Ishigaki et al. (2018),
respectively. The grey circles and white squares represent the $x_{\rm
H\,\textsc{i}}$ estimation from Ly$\alpha$ luminosity function (Ouchi et al.,
2010; Konno et al., 2014; Inoue et al., 2018; Konno et al., 2018; Morales et
al., 2021; Goto et al., 2021) and LAE clustering analysis (Ouchi et al.,
2018), respectively. The white and grey triangles represent the $x_{\rm
H\,\textsc{i}}$ estimation by QSOs Gunn Peterson troughs (Fan et al., 2006)
and Ly$\alpha$+Ly$\beta$ dark fraction (Mesinger et al., 2015), respectively.
The grey and white diamonds represent the $x_{\rm H\,\textsc{i}}$ estimation
from Ly$\alpha$ damping wing absorption of QSOs (Schroeder et al., 2013;
Davies et al., 2018; Greig et al., 2019; Wang et al., 2020) and GRBs (Totani
et al., 2006, 2014). The black diamonds represent the $x_{\rm H\,\textsc{i}}$
estimation from Ly$\alpha$ damping wing absorption profile of individual UV
spectrum of galaxies (Curtis-Lake et al., 2023; Hsiao et al., 2023). The grey
squares represent the $x_{\rm H\,\textsc{i}}$ estimation from Ly$\alpha$ of
Lyman-break galaxies (Mesinger et al., 2015; Hoag et al., 2019; Mason et al.,
2019; Jung et al., 2020; Whitler et al., 2020; Bruton et al., 2023; Morishita
et al., 2023). The white pentagon represents the $x_{\rm H\,\textsc{i}}$
estimation from the cosmic microwave background observations (Planck
Collaboration et al., 2020b).

In Figure 10, we plot our ${x_{\rm H\,\textsc{i}}}$ constraints and the
uncertainties in red circle and bars, respectively. We also show literature
values of $x_{\rm H\,\textsc{i}}$ at different redshift including the
constraint suggested by the electron scattering from the cosmic microwave
background observations (Planck Collaboration et al., 2020b) and Ly$\alpha$
damping wing measurements of QSO spectra (Schroeder et al., 2013; Davies et
al., 2018; Greig et al., 2019; Wang et al., 2020), GRB spectra (Totani et al.,
2006, 2014), and recently obtained JWST galaxy spectra in others works
(Curtis-Lake et al., 2023; Hsiao et al., 2023). Optical depth measurements
using QSOs Gunn Peterson troughs around $z\sim 5-6$ place the constraints on
the end of the epoch of reionization (Fan et al., 2006). We additionally show
other $x_{\rm H\,\textsc{i}}$ constraints from galaxy observations including
the ones inferred from Ly$\alpha$ luminosity function (Ouchi et al., 2010;
Konno et al., 2014; Inoue et al., 2018; Konno et al., 2018; Morales et al.,
2021; Goto et al., 2021), Ly$\alpha$ emitter (LAE) clustering analysis (Ouchi
et al., 2018), Ly$\alpha$+Ly$\beta$ dark fraction Mesinger et al. (2015), and
Ly$\alpha$ emission lines of Lyman-break galaxies (Mesinger et al., 2015; Hoag
et al., 2019; Mason et al., 2019; Jung et al., 2020; Whitler et al., 2020;
Bruton et al., 2023; Morishita et al., 2023). For the $x_{\rm H\,\textsc{i}}$
estimate from Curtis-Lake et al. (2023), we adopt the values shown in Figure 3
of the corresponding paper.

We overplot three distinct evolutionary models by Naidu et al. (2020),
Finkelstein et al. (2019), and Ishigaki et al. (2018) (hereafter, Naidu model,
Finkelstein model, Ishigaki model, respectively) in orange dashed, blue solid,
and purple dotted lines, respectively. We adopt Model II described in Naidu et
al. (2020) as Naidu model. Note that the difference between redshift evolution
of $x_{\rm H\,\textsc{i}}$ in Model I and II in Naidu et al. (2020) is
insignificant compared to relatively large errorbars for current $x_{\rm
H\,\textsc{i}}$ measurements. These evolutionary models differ in the
prescription for the escape fraction $f_{\rm esc}$ and the faint end slope
$\alpha$ of the UV luminosity functions. Naidu model assumes that $f_{\rm
esc}$ correlates with the star formation surface density $\Sigma_{\rm SFR}$
with a shallow faint end slope (i.e, $\alpha>-2$), resulting in a late and
rapid cosmic reionization history driven by a small number of bright galaxies.
Finkelstein model assumes that $f_{\rm esc}$ anti-correlates with halo mass
$M_{\rm halo}$ and a steep faint end slope (i.e, $\alpha<-2$), resulting in an
extended $x_{\rm H\,\textsc{i}}$ evolution driven by many faint galaxies.
Ishigaki model assumes the constant $f_{\rm esc}$ value at $\simeq 0.17$ with
a steep UV faint end slope (i.e, $\alpha<-2$). Ishigaki model results in a
moderately late cosmic reionization consistent with the one suggested by an
electron scattering of the cosmic microwave background. We summarize the
prescriptions for each model in Table 4. As shown in Figure 10, the best
estimate of our $x_{\rm H\,\textsc{i}}$ consistent with all three evolutionary
models and the constraints inferred from the electron scattering of CMB.

In Figure 10, we also overplot the $x_{\rm H\,\textsc{i}}$ values at $z\approx
7$ estimated by Umeda (2023) from Ly$\alpha$ luminosity function and
clustering properties of LAEs as pink circle and square, respectively. Umeda
(2023) derive Ly$\alpha$ luminosity function and LAE spatial angular
correlation function using the $z=6-7$ photometric LAE catalogs constructed by
Ono et al. (2021) from the Public Data Release 2 data of the Hyper Suprime-Cam
Subaru Strategic Programs (Aihara et al., 2019). Umeda (2023) estimate $x_{\rm
H\,\textsc{i}}$ from Ly$\alpha$ luminosity function and angular correlation
function of LAEs in a similar manner as Konno et al. (2018), Ouchi et al.
(2018), and Goto et al. (2021). The $z=6-7$ LAE catalog of Ono et al. (2021)
is the largest in sample size with a total of 3561 galaxies, thus the $x_{\rm
H\,\textsc{i}}$ estimates from Umeda (2023) have advantages in terms of the
reduced statistical uncertainty for estimated Ly$\alpha$ luminosity function
and angular correlation functions. For the $x_{\rm H\,\textsc{i}}$ estimate
from Ly$\alpha$ luminosity function, Umeda (2023) first compared the
luminosity density at $z=5.7$ and $z=6.6$ to derive IGM transmission of
Ly$\alpha$. Then, Umeda (2023) estimate $x_{\rm H\,\textsc{i}}$ from the IGM
transmission using the analytical relation between IGM transmission of
Ly$\alpha$, $x_{\rm H\,\textsc{i}}$, and the ionized bubble size from Dijkstra
(2014) and the semi-numerical relation between $x_{\rm H\,\textsc{i}}$ and
$R_{b}$ from Furlanetto & Oh (2005). For the $x_{\rm H\,\textsc{i}}$ estimate
from LAE clustering analysis, Umeda (2023) compared the derived galaxy bias
evolution from $z=5.7$ to $z=6.6$ to the semi-numerical prediction of galaxy
bias at $z=6.6$ with different $x_{\rm H\,\textsc{i}}$ by Furlanetto et al.
(2006). By comparing the $x_{\rm H\,\textsc{i}}$ estimate at $z>7$ with the
one at $z<7$ estimated by Umeda (2023), we can confirm smooth monotonic
increase of $x_{\rm H\,\textsc{i}}$ towards the higher redshift. Other $x_{\rm
H\,\textsc{i}}$ constraints at $z>7$ from the literature are mostly consistent
with our measurements within error. To summarize, our $x_{\rm H\,\textsc{i}}$
measurements at $z>7$ agree with various constraints, including the ones
previously inferred from the luminosity function and clustering properties of
LAEs at $z<7$, the electron scattering of CMB, and the evolution of UV
luminosity function assuming a constant $f_{\rm esc}\sim 0.2$.

Table 4: $x_{\rm H\,\textsc{i}}$ Evolution Models

Model | Naidu | Ishigaki | Finkelstein
---|---|---|---
$f_{\rm esc}$ | $\propto\Sigma_{\rm SFR}^{0.4}$ | 0.17 | $\propto M_{\rm halo}^{-1}$
$\alpha$ | $>-2$ | $<-2$ | $<-2$

Note. — Naidu, Ishigaki, and Finkelstein models corresponds to $x_{\rm
H\,\textsc{i}}$ evolution model presented in Naidu et al. (2020), Ishigaki et
al. (2018), and Finkelstein et al. (2019), respectively. We adopt Model II of
the Naidu et al. (2020) as Naidu model.

###  4.4 Bubble Size Evolution

Beside $x_{\rm H\,\textsc{i}}$, the size of the ionized bubble is another
tracer of cosmic reionization scenarios. Various studies suggest that in the
case of inhomogeneous cosmic reionization driven by galaxies, the ionized
region forms around the galaxies and their typical size correlates with
$x_{\rm H\,\textsc{i}}$ (e.g., Ciardi & Madau, 2003; Furlanetto & Oh, 2005;
Mesinger & Furlanetto, 2008a, b; Lu et al., 2023). The size of ionized bubble
at the EoR has been inferred from the Ly$\alpha$ emission line profiles and
equivalent widths (e.g., Mason & Gronke, 2020; Endsley & Stark, 2022; Hayes &
Scarlata, 2023; Jung et al., 2023; Saxena et al., 2023). We instead use
Ly$\alpha$ damping wing profiles to investigate $R_{b}$ and its evolution with
the redshift. We use the best-fit parameters presented in Section 3 to see the
relation between $R_{b}$ and $x_{\rm H\,\textsc{i}}$. Figure 11 shows the
estimated $R_{b}$ values for each redshift bin from our analysis in Section 3
together with the relation between $x_{\rm H\,\textsc{i}}$ and $R_{b}$
predicted from analytical model by Furlanetto & Oh (2005) in black dotted
lines. Together with our measurements, we show $x_{\rm H\,\textsc{i}}$ and
$R_{b}$ constrained via damping wing absorption analysis in Hsiao et al.
(2023). We confirm that our constraint for the redshift bin of $\langle
z\rangle=10.28$ is consistent with $x_{\rm H\,\textsc{i}}$ and $R_{b}$
measurements for $z=10.17$ object by Hsiao et al. (2023). In contrast, the
recent bubble size measurement by Fujimoto et al. (2023b) for $z=8.51$ object
suggests a bubble size of $R_{b}=71.9\pm 1.0$ cMpc. This bubble size value is
about a dex larger than our constraint at the redshift bin of $\langle
z\rangle=8.28$. As Fujimoto et al. (2023b) discusses in their paper, this
large offset of bubble size measurement could be physically attributed to
specifically strong ionizing radiation by two possible AGN candidates with
only a separation of 380 kpc within the same ionized bubble. We also plot the
$x_{\rm H\,\textsc{i}}$ and $R_{b}$ relation for the rapid and gradual cosmic
reionization history predicted from a semi-numerical simulations by Lu et al.
(2023) in solid and dashed lines, respectively. Lu et al. (2023)’s predictions
take into account of the overlap of the ionized regions formed around UV
bright ($M_{\rm UV}<-20$) galaxies that are comparable to those for the
majority of galaxies in our sample. Incorporation of the overlaps between
ionized yields systematically larger $R_{b}$ at each $x_{\rm H\,\textsc{i}}$
compared with the ones predicted by Furlanetto & Oh (2005), which do not
consider the overlaps. As shown in Figure 11, our $R_{b}$ estimate at each
$x_{\rm H\,\textsc{i}}$ is up to few dex larger than the one predicted by
theoretical relations of Furlanetto & Oh (2005), shown as the dotted line.
However, our $R_{b}$ estimates are consistent with in errorbars with the ones
predicted by Lu et al. (2023) that consider the overlaps of ionized regions,
confirming the importance of considering the overlaps to characterize the
typical ionized bubble shapes during the EoR. Precise $R_{b}$ measurements in
the future would be useful to distinguish different cosmic reionization
scenarios that yield different $x_{\rm H\,\textsc{i}}$ and $R_{b}$ relations.

Figure 11: $x_{\rm H\,\textsc{i}}$ and $R_{b}$ estimates from different
redshift bins. The red squares and bars represent our estimates and the
errors, respectively, of ($x_{\rm H\,\textsc{i}}$, $\log R_{b}$) at the
redshift bin #1, #2, #3, and #4. Black square and arrows corresponds to
1$\sigma$ upper and lower limits of $x_{\rm H\,\textsc{i}}$ and $R_{b}$,
respectively, on MACS0647_JD ($z=10.17$) by Hsiao et al. (2023). Faint grey
circles and bars represent the best-fit value and its error for individual
galaxies. Solid and dashed lines correspond to the theoretical prediction by
Lu et al. (2023) on the characteristic ionized bubble radii at different
$x_{\rm H\,\textsc{i}}$ assuming rapid and gradual cosmic reionization
considering overlapping ionized regions, respectively. The dotted line
corresponds to the analytical relation between characteristic ionized bubble
radii and $x_{\rm H\,\textsc{i}}$ without considering overlapping ionized
regions (Furlanetto & Oh, 2005).

###  4.5 Possible Systematics

Although we consider the ionized bubbles surrounding the galaxies and H i
absorption by CGM when we fit the stacked spectra for the Ly$\alpha$ damping
wing absorptions, there are still some possible systematics that we could
consider for more accurate estimation. For example, the results from semi-
numerical simulations suggest that our method assuming homogeneous IGM could
overestimate the $x_{\rm H\,\textsc{i}}$ values by $\sim 30\%$ at the middle
stage of cosmic reionization (Mesinger & Furlanetto, 2008a, b). Keating et al.
(2023) also suggest the simple analytical formulation using Miralda-Escudé
(1998) does miss out some of the physics related Ly$\alpha$ absorption that
are seen in the simulation results (i.e., the contributions of residual H i
gas inside the ionized bubble to the Ly$\alpha$ optical depth). On the
otherhand, the results from different numerical simulation with various box
size yield different predictions on the evolution of key cosmic reionization
parameters (i.e., Keating et al., 2023; Lu et al., 2023), indicating better
interpretation of the Ly$\alpha$ damping wing absorption measurements requires
further improvements in both simulation and observations. Another source of
systematics is a biased sampling of galaxies in physical properties such as
UV-luminosity. Although we could use the simulation prediction suitable for
the observed galaxies (e.g., UV-bright galaxies residing at overdensed region;
Lu et al., 2023) to infer cosmic reionization parameters, we ideally need to
know the optical depths of as many unbiased sight lines as possible to
precisely estimate the volume-averaged neutral fraction of the universe with a
minimum dependency on the simulation models. Moreover, we must investigate the
evolution and diversity of physical properties of the intrinsic UV spectrum of
galaxies that are still not well studied above $z>7$ (e.g., UV-slope, ionizing
photon production efficiency, etc.). In terms of the systematics from data
reduction, the damping wing analysis with low-resolution PRISM data introduces
uncertainties in spectroscopic redshift and the degenerating components around
Lyman break (e.g., Ly$\alpha$ emission, Ly$\alpha$ absorption by IGM/CGM).
Additionally, as Jones et al. (2023) shown, NIRSpec/PRISM spectra have
systematics in the wavelength resolution depending on the compactness of the
object. This could also be a potential source of the systematics in the
redshift determination and damping wing absorption measurement using PRISM
spectra. The redshift uncertainty of $\Delta z\sim\pm 0.01$ as suggested by
Jones et al. (2023) give a slight systematics of $\sim 0.1$ for HPDI of
$x_{\rm H\,\textsc{i}}$ for the individual galaxy. This systematics is
insignificant compared to the current large uncertainty of $x_{\rm
H\,\textsc{i}}$ ($\sim 0.5$). However, to further increase the precision of
Ly$\alpha$ damping wing measurement of individual galaxy, the high $S/N$ and
resolution data is admirable to determine systemic redshift precisely and
break down the unresolved spectral features around rest-frame 1216 Å. The
systematics and statistical uncertainties discussed above should be reduced
via the future works including more observations and improvements in data
reduction of $z>7$ galaxies with JWST/NIRSpec.

##  5 Conclusion

We present our constraints on the $x_{\rm H\,\textsc{i}}$ and $R_{b}$ at
$z=7-12$ via analysing Ly$\alpha$ damping wing absorption of 26 galaxies. We
gathered JWST/NIRSpec’s PRISM data from multiple surveys to create the galaxy
sample. Out sample consists of 26 galaxies and the systemic redshifts of all
galaxies have been confirmed with the emission lines including including [O
iii]5007, H$\beta$, and H$\gamma$. We then construct 4 composite spectra
binned by the redshifts and confirm the softening break feature near the rest-
frame 1216 Å, suggesting the strengthening of Ly$\alpha$ damping wing
absorptions toward the high redshifts. We then fit the galaxy spectra by the
stellar model spectra with Ly$\alpha$ damping wing absorptions that changes
with $x_{\rm H\,\textsc{i}}$ and $R_{b}$.

  1. 1.

Stacked spectra at four redshift bins show the the transition from sharpened
to softened break feature around rest-frame 1216 Å towards the high redshift,
suggesting the increasing $x_{\rm H\,\textsc{i}}$. We confirm that the median
estimates for $x_{\rm H\,\textsc{i}}$ ($R_{b}$) values for the stacked spectra
at $\langle z\rangle=$7.12, 7.44, 8.28, and 10.28 show increasing (decreasing)
trend toward the higher redshift bins.

  2. 2.

The best-fit estimates for $x_{\rm H\,\textsc{i}}$ at $z>7$ by Ly$\alpha$
damping wing absorptions of galaxies’ UV continua is consistent with
moderately late reionization history that is consistent with the $x_{\rm
H\,\textsc{i}}$ constraints previously suggested from CMB observations and the
evolution of UV luminosity function with $f_{\rm esc}\sim 0.2$. We note,
however, that the errors of our $x_{\rm H\,\textsc{i}}$ estimates are still
large. We need to reduce statistical errors as well as possible systematics to
narrow down the range of possible cosmic reionization scenarios.

  3. 3.

We compare estimated $x_{\rm H\,\textsc{i}}$ and $R_{b}$ values with the
theoretical predictions of the typical ionized bubble radii at different
$x_{\rm H\,\textsc{i}}$. Our $R_{b}$ estimates lie up to a few dex above the
cosmic average values estimated from analytic calculations by Furlanetto & Oh
(2005), while our $R_{b}$ estimates are comparable with the values for
overlapping ionized bubbles around bright galaxies predicted by numerical
simulations of (Lu et al., 2023).

## Acknowledgements

We thank Andrei Mesinger, Anne Hutter, Adi Zitrin, Akio Inoue, Hidenobu
Yajima, Hayato Shimabukuro, Daichi Kashino, Satoshi Kikuta, Yoshinobu
Fudamoto, Kazuyuki Omukai, Shintaro Yoshiura, Akinori Matsumoto, and Minami
Nakane for the valuable discussions on this work. This work is based on
observations made with the NASA/ESA/CSA James Webb Space Telescope. The data
were obtained from the Mikulski Archive for Space Telescopes at the Space
Telescope Science Institute, which is operated by the Association of
Universities for Research in Astronomy, Inc., under NASA contract NAS 5-03127
for JWST. These observations are associated with programs ERS-1345 (CEERS),
GO-1433, DDT-2750, and GTO-1210 (JADES). The authors acknowledge the CEERS,
GO-1433, and DDT-2750 teams led by Steven L. Finkelstein, Dan Coe, Pablo
Arrabal Haro, and Nora Lützgendorf, respectively, for developing their
observing programs with a zero-exclusive-access period. We thank the JADES
team for publicly releasing reduced spectra and catalog from the JADES survey.
We thank Fergus Cullen and the team for publicly releasing the composite UV
spectra from VANDELS survey. This publication is based upon work supported by
the World Premier International Research Center Initiative (WPI Initiative),
MEXT, Japan, and KAKENHI (20H00180, 21J20785, 21K13953, 21H04467, 23KJ0646,
23KJ0589) through Japan Society for the Promotion of Science. This work was
supported by JSPS Core-to-Core Program (grant number: JPJSCCA20210003). This
work was supported by the joint research program of the Institute for Cosmic
Ray Research (ICRR), University of Tokyo. This research was supported by FoPM,
WINGS Program, the University of Tokyo.

## References

  * Aihara et al. (2019) Aihara, H., AlSayyad, Y., Ando, M., et al. 2019, PASJ, 71, 114, doi: 10.1093/pasj/psz103
  * Arrabal Haro et al. (2023a) Arrabal Haro, P., Dickinson, M., Finkelstein, S. L., et al. 2023a, arXiv e-prints, arXiv:2304.05378, doi: 10.48550/arXiv.2304.05378
  * Arrabal Haro et al. (2023b) —. 2023b, arXiv e-prints, arXiv:2303.15431, doi: 10.48550/arXiv.2303.15431
  * Atek et al. (2014) Atek, H., Kunth, D., Schaerer, D., et al. 2014, A&A, 561, A89, doi: 10.1051/0004-6361/201321519
  * Bañados et al. (2018) Bañados, E., Venemans, B. P., Mazzucchelli, C., et al. 2018, Nature, 553, 473, doi: 10.1038/nature25180
  * 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
  * Bruton et al. (2023) Bruton, S., Lin, Y.-H., Scarlata, C., & Hayes, M. J. 2023, arXiv e-prints, arXiv:2303.03419, doi: 10.48550/arXiv.2303.03419
  * Bunker et al. (2023) Bunker, A. J., Cameron, A. J., Curtis-Lake, E., et al. 2023, arXiv e-prints, arXiv:2306.02467, doi: 10.48550/arXiv.2306.02467
  * Byler et al. (2017) Byler, N., Dalcanton, J. J., Conroy, C., & Johnson, B. D. 2017, ApJ, 840, 44, doi: 10.3847/1538-4357/aa6c66
  * Calzetti et al. (2000) Calzetti, D., Armus, L., Bohlin, R. C., et al. 2000, ApJ, 533, 682, doi: 10.1086/308692
  * Ciardi & Madau (2003) Ciardi, B., & Madau, P. 2003, ApJ, 596, 1, doi: 10.1086/377634
  * Curtis-Lake et al. (2023) Curtis-Lake, E., Carniani, S., Cameron, A., et al. 2023, Nature Astronomy, doi: 10.1038/s41550-023-01918-w
  * Davies et al. (2018) Davies, F. B., Hennawi, J. F., Bañados, E., et al. 2018, ApJ, 864, 142, doi: 10.3847/1538-4357/aad6dc
  * Dijkstra (2014) Dijkstra, M. 2014, PASA, 31, e040, doi: 10.1017/pasa.2014.33
  * Eldridge et al. (2017) Eldridge, J. J., Stanway, E. R., Xiao, L., et al. 2017, PASA, 34, e058, doi: 10.1017/pasa.2017.51
  * Endsley & Stark (2022) Endsley, R., & Stark, D. P. 2022, MNRAS, 511, 6042, doi: 10.1093/mnras/stac524
  * Fan et al. (2006) Fan, X., Strauss, M. A., Becker, R. H., et al. 2006, AJ, 132, 117, doi: 10.1086/504836
  * 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. (2019) Finkelstein, S. L., D’Aloisio, A., Paardekooper, J.-P., et al. 2019, ApJ, 879, 36, doi: 10.3847/1538-4357/ab1ea8
  * 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. (2023a) Fujimoto, S., Arrabal Haro, P., Dickinson, M., et al. 2023a, arXiv e-prints, arXiv:2301.09482, doi: 10.48550/arXiv.2301.09482
  * Fujimoto et al. (2023b) Fujimoto, S., Wang, B., Weaver, J., et al. 2023b, arXiv e-prints, arXiv:2308.11609, doi: 10.48550/arXiv.2308.11609
  * Furlanetto & Oh (2005) Furlanetto, S. R., & Oh, S. P. 2005, MNRAS, 363, 1031, doi: 10.1111/j.1365-2966.2005.09505.x
  * Furlanetto et al. (2006) Furlanetto, S. R., Zaldarriaga, M., & Hernquist, L. 2006, MNRAS, 365, 1012, doi: 10.1111/j.1365-2966.2005.09785.x
  * Gardner et al. (2023) Gardner, J. P., Mather, J. C., Abbott, R., et al. 2023, PASP, 135, 068001, doi: 10.1088/1538-3873/acd1b5
  * Goto et al. (2021) Goto, H., Shimasaku, K., Yamanaka, S., et al. 2021, ApJ, 923, 229, doi: 10.3847/1538-4357/ac308b
  * Greig et al. (2019) Greig, B., Mesinger, A., & Bañados, E. 2019, MNRAS, 484, 5094, doi: 10.1093/mnras/stz230
  * Harikane et al. (2023) Harikane, Y., Nakajima, K., Ouchi, M., et al. 2023, arXiv e-prints, arXiv:2304.06658, doi: 10.48550/arXiv.2304.06658
  * Hayes & Scarlata (2023) Hayes, M. J., & Scarlata, C. 2023, arXiv e-prints, arXiv:2303.03160, doi: 10.48550/arXiv.2303.03160
  * Heintz et al. (2023) Heintz, K. E., Watson, D., Brammer, G., et al. 2023, arXiv e-prints, arXiv:2306.00647, doi: 10.48550/arXiv.2306.00647
  * Hernandez et al. (2020) Hernandez, S., Aloisi, A., James, B. L., et al. 2020, ApJ, 892, 19, doi: 10.3847/1538-4357/ab77c6
  * Hoag et al. (2019) Hoag, A., Bradač, M., Huang, K., et al. 2019, ApJ, 878, 12, doi: 10.3847/1538-4357/ab1de7
  * Hsiao et al. (2023) Hsiao, T. Y.-Y., Abdurro’uf, Coe, D., et al. 2023, arXiv e-prints, arXiv:2305.03042, doi: 10.48550/arXiv.2305.03042
  * Inoue et al. (2018) Inoue, A. K., Hasegawa, K., Ishiyama, T., et al. 2018, PASJ, 70, 55, doi: 10.1093/pasj/psy048
  * Ishigaki et al. (2018) Ishigaki, M., Kawamata, R., Ouchi, M., et al. 2018, ApJ, 854, 73, doi: 10.3847/1538-4357/aaa544
  * 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
  * Johnson et al. (2021) Johnson, B. D., Leja, J., Conroy, C., & Speagle, J. S. 2021, ApJS, 254, 22, doi: 10.3847/1538-4365/abef67
  * Jones et al. (2023) Jones, G. C., Bunker, A. J., Saxena, A., et al. 2023, arXiv e-prints, arXiv:2306.02471, doi: 10.48550/arXiv.2306.02471
  * Jung et al. (2020) Jung, I., Finkelstein, S. L., Dickinson, M., et al. 2020, ApJ, 904, 144, doi: 10.3847/1538-4357/abbd44
  * 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
  * Keating et al. (2023) Keating, L. C., Bolton, J. S., Cullen, F., et al. 2023, arXiv e-prints, arXiv:2308.05800, doi: 10.48550/arXiv.2308.05800
  * Konno et al. (2014) Konno, A., Ouchi, M., Ono, Y., et al. 2014, ApJ, 797, 16, doi: 10.1088/0004-637X/797/1/16
  * Konno et al. (2018) Konno, A., Ouchi, M., Shibuya, T., et al. 2018, PASJ, 70, S16, doi: 10.1093/pasj/psx131
  * Larson et al. (2022) Larson, R. L., Finkelstein, S. L., Hutchison, T. A., et al. 2022, ApJ, 930, 104, doi: 10.3847/1538-4357/ac5dbd
  * Larson et al. (2023) Larson, R. L., Finkelstein, S. L., Kocevski, D. D., et al. 2023, arXiv e-prints, arXiv:2303.08918, doi: 10.48550/arXiv.2303.08918
  * Lu et al. (2023) Lu, T.-Y., Mason, C., Hutter, A., et al. 2023, arXiv e-prints, arXiv:2304.11192, doi: 10.48550/arXiv.2304.11192
  * Mason & Gronke (2020) Mason, C. A., & Gronke, M. 2020, MNRAS, 499, 1395, doi: 10.1093/mnras/staa2910
  * Mason et al. (2019) Mason, C. A., Fontana, A., Treu, T., et al. 2019, MNRAS, 485, 3947, doi: 10.1093/mnras/stz632
  * Matsuoka et al. (2023) Matsuoka, Y., Onoue, M., Iwasawa, K., et al. 2023, arXiv e-prints, arXiv:2305.11225, doi: 10.48550/arXiv.2305.11225
  * McKinney et al. (2019) McKinney, J. H., Jaskot, A. E., Oey, M. S., et al. 2019, ApJ, 874, 52, doi: 10.3847/1538-4357/ab08eb
  * Mesinger et al. (2015) Mesinger, A., Aykutalp, A., Vanzella, E., et al. 2015, MNRAS, 446, 566, doi: 10.1093/mnras/stu2089
  * Mesinger & Furlanetto (2008a) Mesinger, A., & Furlanetto, S. R. 2008a, MNRAS, 385, 1348, doi: 10.1111/j.1365-2966.2007.12836.x
  * Mesinger & Furlanetto (2008b) —. 2008b, MNRAS, 386, 1990, doi: 10.1111/j.1365-2966.2008.13039.x
  * Miralda-Escudé (1998) Miralda-Escudé, J. 1998, ApJ, 501, 15, doi: 10.1086/305799
  * Morales et al. (2021) Morales, A. M., Mason, C. A., Bruton, S., et al. 2021, ApJ, 919, 120, doi: 10.3847/1538-4357/ac1104
  * Morishita et al. (2023) Morishita, T., Roberts-Borsani, G., Treu, T., et al. 2023, ApJ, 947, L24, doi: 10.3847/2041-8213/acb99e
  * Mortlock et al. (2011) Mortlock, D. J., Warren, S. J., Venemans, B. P., et al. 2011, Nature, 474, 616, doi: 10.1038/nature10159
  * Naidu et al. (2020) Naidu, R. P., Tacchella, S., Mason, C. A., et al. 2020, ApJ, 892, 109, doi: 10.3847/1538-4357/ab7cc9
  * Nakajima et al. (2023) Nakajima, K., Ouchi, M., Isobe, Y., et al. 2023, arXiv e-prints, arXiv:2301.12825, doi: 10.48550/arXiv.2301.12825
  * Oke & Gunn (1983) Oke, J. B., & Gunn, J. E. 1983, ApJ, 266, 713, doi: 10.1086/160817
  * Ono et al. (2021) Ono, Y., Itoh, R., Shibuya, T., et al. 2021, ApJ, 911, 78, doi: 10.3847/1538-4357/abea15
  * Ouchi et al. (2010) Ouchi, M., Shimasaku, K., Furusawa, H., et al. 2010, ApJ, 723, 869, doi: 10.1088/0004-637X/723/1/869
  * Ouchi et al. (2018) Ouchi, M., Harikane, Y., Shibuya, T., et al. 2018, PASJ, 70, S13, doi: 10.1093/pasj/psx074
  * Planck Collaboration et al. (2020a) Planck Collaboration, Aghanim, N., Akrami, Y., et al. 2020a, A&A, 641, A6, doi: 10.1051/0004-6361/201833910
  * Planck Collaboration et al. (2020b) —. 2020b, A&A, 641, A6, doi: 10.1051/0004-6361/201833910
  * Rieke et al. (2023) Rieke, M. J., Kelly, D. M., Misselt, K., et al. 2023, PASP, 135, 028001, doi: 10.1088/1538-3873/acac53
  * Rigby et al. (2023) Rigby, J., Perrin, M., McElwain, M., et al. 2023, PASP, 135, 048001, doi: 10.1088/1538-3873/acb293
  * Roberts-Borsani et al. (2022) Roberts-Borsani, G., Treu, T., Chen, W., et al. 2022, arXiv e-prints, arXiv:2210.15639, doi: 10.48550/arXiv.2210.15639
  * Sanders et al. (2023) Sanders, R. L., Shapley, A. E., Topping, M. W., Reddy, N. A., & Brammer, G. B. 2023, arXiv e-prints, arXiv:2303.08149, doi: 10.48550/arXiv.2303.08149
  * Saxena et al. (2023) Saxena, A., Robertson, B. E., Bunker, A. J., et al. 2023, arXiv e-prints, arXiv:2302.12805, doi: 10.48550/arXiv.2302.12805
  * Schroeder et al. (2013) Schroeder, J., Mesinger, A., & Haiman, Z. 2013, MNRAS, 428, 3058, doi: 10.1093/mnras/sts253
  * Tang et al. (2023) Tang, M., Stark, D. P., Chen, Z., et al. 2023, arXiv e-prints, arXiv:2301.07072, doi: 10.48550/arXiv.2301.07072
  * Tanvir et al. (2019) Tanvir, N. R., Fynbo, J. P. U., de Ugarte Postigo, A., et al. 2019, MNRAS, 483, 5380, doi: 10.1093/mnras/sty3460
  * Totani et al. (2006) Totani, T., Kawai, N., Kosugi, G., et al. 2006, PASJ, 58, 485, doi: 10.1093/pasj/58.3.485
  * Totani et al. (2014) Totani, T., Aoki, K., Hattori, T., et al. 2014, PASJ, 66, 63, doi: 10.1093/pasj/psu032
  * Umeda (2023) Umeda, H. 2023, Master’s thesis, The Univ. of Tokyo, Tokyo, Japan
  * 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
  * Vousden et al. (2021) Vousden, W., Farr, W. M., & Mandel, I. 2021, ptemcee: A parallel-tempered version of emcee, Astrophysics Source Code Library, record ascl:2101.006.  http://ascl.net/2101.006
  * Wang et al. (2019) Wang, F., Yang, J., Fan, X., et al. 2019, ApJ, 884, 30, doi: 10.3847/1538-4357/ab2be5
  * Wang et al. (2020) Wang, F., Davies, F. B., Yang, J., et al. 2020, ApJ, 896, 23, doi: 10.3847/1538-4357/ab8c45
  * Whitler et al. (2020) Whitler, L. R., Mason, C. A., Ren, K., et al. 2020, MNRAS, 495, 3602, doi: 10.1093/mnras/staa1178
  * Williams et al. (2023) Williams, H., Kelly, P. L., Chen, W., et al. 2023, Science, 380, 416, doi: 10.1126/science.adf5307
  * Zitrin et al. (2015) Zitrin, A., Labbé, I., Belli, S., et al. 2015, ApJ, 810, L12, doi: 10.1088/2041-8205/810/1/L12