# Damping wing absorption associated with a giant Ly$\alpha$ trough at $z<6$:
direct evidence for late-ending reionization

George D. Becker,1 James S. Bolton,2 Yongda Zhu,1,3 and Seyedazim Hashemi1
1Department of Physics & Astronomy, University of California, Riverside, CA,
92521, USA
2School of Physics and Astronomy, University of Nottingham, University Park,
Nottingham, NG7 2RD, UK
3Steward Observatory, University of Arizona, 933 N Cherry Ave, Tucson, AZ
85721, USA E-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

Multiple observations now suggest that the hydrogen reionization may have
ended well below redshift six. While there has previously been no conclusive
proof of extended neutral islands in the $z<6$ intergalactic medium, it is
possible that such islands give rise to the giant Ly$\alpha$ absorption
troughs seen in the spectra of high-redshift quasars. Here we present evidence
that the deepest and longest-known Ly$\alpha$ trough at $z<6$, towards ULAS
J0148+0600 (J0148), is associated with damping wing absorption. The evidence
comes from a window of strong Ly$\alpha$ transmission at the edge of the J0148
proximity zone. We show that the relatively smooth profile of this
transmission window is highly unlikely to arise from resonant absorption
alone, but is consistent with the presence of a damping wing. We further argue
that the damping wing is unlikely to arise from a compact source due to the
lack of associated metal lines, and is more likely to arise from an extended
neutral island associated with the giant Ly$\alpha$ trough. We investigate the
physical conditions that may give rise to the strong transmission window, and
speculate that it may signal an usually deep void, nearby ionizing sources,
and/or the recent passage of an ionization front.

###### keywords:

intergalactic medium – quasars: absorption lines – cosmology: observations –
dark ages, reionization, first stars – large-scale structure of the Universe

††pubyear: 2023††pagerange: Damping wing absorption associated with a giant
Ly$\alpha$ trough at $z<6$: direct evidence for late-ending reionization–C

##  1 Introduction

The reionization of the intergalactic medium (IGM) was a milestone in cosmic
history. As a global baryonic phase change, it represents one of the most
significant interactions between galaxies and their environments. Determining
when and how reionization occurred is therefore essential for developing a
self-consistent picture of galaxy evolution and IGM physics in the early
Universe (for recent reviews, see Dayal & Ferrara, 2018; Wise, 2019;
Robertson, 2022; Choudhury, 2022; Fan et al., 2023).

Over the past several years, a consensus has begun to emerge that reionization
may have ended significantly later than previously believed. Ly$\alpha$
transmission in the spectra of the first known $z\sim 6$ quasars signaled that
reionization was at least largely complete by that redshift (e.g., Fan et al.,
2006). With larger samples and higher-quality spectra, however, it has become
clear that the high degree of scatter in the large-scale Ly$\alpha$ forest
transmission is not consistent with a fully reionized IGM, at least not one
with a uniform ionizing UV background and temperature-density relation (Becker
et al., 2015; D’Aloisio et al., 2015; D’Aloisio et al., 2018; Davies &
Furlanetto, 2016; Bosman et al., 2018; Bosman et al., 2022; Eilers et al.,
2018; Yang et al., 2020b). Models wherein reionization ended closer to $z\sim
5$, however, can naturally reproduce this scatter while remaining consistent
with other constraints on reionization such as the electron optical depth to
the cosmic microwave background and the apparent attenuation of Ly$\alpha$
emission from quasars and galaxies at $z\gtrsim 7$ (Kulkarni et al., 2019;
Keating et al., 2020a, b; Nasir & D’Aloisio, 2020; Choudhury et al., 2021; Qin
et al., 2021). A strong indication that reionization was still underway at
$z\sim 6$ comes from the rapid evolution of the mean free path of ionizing
photons, which has been measured directly from the Lyman continuum
transmission of quasar spectra (Becker et al., 2021; Bosman, 2021; Zhu et al.,
2023) and indirectly from large-scale Ly$\alpha$ opacity fluctuations (Gaikwad
et al., 2023; Davies et al., 2023b). The increase in the mean free path from
$z\sim 6$ to 5 is more rapid than would be expected for a fully reionized IGM
(D’Aloisio et al., 2020), but is broadly consistent with the photoevaporation
of ionizing sinks at the end of reionization (Keating et al., 2020a; Cain et
al., 2021; Garaldi et al., 2022; Lewis et al., 2022).

While a late-ending reionization scenario is supported by a number of
constraints, obtaining direct observational evidence of neutral gas in the
$z<6$ IGM has been challenging due to the fact that resonant Ly$\alpha$
absorption saturates for even very small hydrogen neutral fractions ($x_{\rm
H\,I}\gtrsim 10^{-4}$) . The best potential tracers of neutral regions are the
“dark gaps” in the Ly$\alpha$ and Ly$\beta$ forests (e.g., Djorgovski et al.,
2001; Becker et al., 2001; Songaila & Cowie, 2002; Fan et al., 2006; Zhu et
al., 2021, 2022). These are extended ($\gtrsim$10 $h^{-1}\,$Mpc) regions of
nearly complete absorption, which in multiple cases have been found to trace
low-density regions based on the number density of surrounding galaxies
(Becker et al., 2018; Kashino et al., 2020; Christenson et al., 2021; Ishimoto
et al., 2022). Dark gaps are naturally explained as absorption by the last
remaining neutral islands at the end of reionization (e.g., Keating et al.,
2020a; Nasir & D’Aloisio, 2020). Observationally, however, it is difficult to
distinguish neutral islands from ionized regions of high opacity that could
arise from large fluctuations in the ionizing UV background (Davies &
Furlanetto, 2016; Nasir & D’Aloisio, 2020).

Figure 1: X-Shooter spectrum of ULAS J0148$+$0600\. The black line shows the
observed flux, in arbitrary $F_{\lambda}$ units, binned to 50 km s-1 pixels
for display. The short yellow lines are spline fits to the flux over the red
side of the Ly$\alpha$+N v emission line and the blueshifted C iv line. The
blue line is a composite continuum constructed from the spectra of 200 lower-
redshift quasars chosen to match the spline fits. The continuous red line and
shaded region show a B-spline fit to the composite and its 1$\sigma$
uncertainty. See text for details. The dotted line shows our adopted power-law
continuum for the non-ionizing UV emission, $L_{\nu}\propto\nu^{-0.6}$
($F_{\lambda}\propto\lambda^{-1.4}$). The depression in the flux with respect
to the composite near $10000$ Å is C iv broad absorption. The 110 $h^{-1}$
comoving Mpc Ly$\alpha$ trough appears as the region of zero flux over 7930 to
8362 Å. The transmission feature examined in this work is centered at
$\sim$8385 Å.

In this work, and in a companion paper by Zhu et al. (in prep), we present the
first direct evidence that at least some of the extended Ly$\alpha$+Ly$\beta$
troughs at $z<6$ do indeed arise from neutral islands. Both works demonstrate
that these troughs are associated with damping wing absorption, the hallmark
of neutral gas. Zhu et al. use stacks to show statistical evidence of damping
wing absorption in the Ly$\alpha$ forest adjacent to dark gaps. Here we show
that one specific dark gap produces a damping wing, and argue that the damping
absorption is more likely to be caused by an extended neutral island than by a
compact absorber such as a damped Ly$\alpha$ system (DLA).

The dark gap we examine is part of the giant Ly$\alpha$ trough towards ULAS
J0148$+$0600 (herein J0148). This trough is remarkable partly because it is
the longest ($\Delta\ell\simeq 1110$ $h^{-1}\,$Mpc) and most opaque
($\tau_{\rm eff}\gtrsim 7$)111$\tau_{\rm eff}=-\log{F}$, where $F$ is the
continuum-normalized transmission. absorption feature known at $z<6$ (Becker
et al., 2018). As we describe below, it is also remarkable because the red end
of the trough is adjacent to an unusually strong and smooth transmission
feature at the end of the quasar proximity zone. We will show that the
characteristics of this feature strongly indicate the presence of damping wing
absorption.

This work shares similarities with studies that have identified damping wings
in the spectra of quasars at $z>7$ (Mortlock et al., 2011; Greig et al., 2017,
2019; Greig et al., 2022; Bañados et al., 2018; Davies et al., 2018; Wang et
al., 2020; Yang et al., 2020a), and recent efforts that have extended this
search down to $z=6$ (Ďurovčíková et al., 2024), but there are two key
differences. First, the damping wing signature is at $z<6$, where no direct
detection of neutral gas in the IGM has yet been made. Second, the evidence
for a damping wing comes from Ly$\alpha$ transmission in the proximity zone,
rather than from the otherwise unabsorbed continuum redward of the Ly$\alpha$
forest (although Davies et al. (2018) include both the proximity zone and the
red-side continuum). Malloy & Lidz (2015) explored stacking $z<6$ dark gaps to
statistically detect damping wings in the presence of strong Ly$\alpha$ forest
absorption. This is the approach taken by Zhu et al. (in prep), while here we
focus on the damping wing from a single Ly$\alpha$ trough. Our work is perhaps
most closely related to Mesinger & Haiman (2004, 2007), who argued that the
distribution of Ly$\alpha$ fluxes near the proximity zone boundaries in three
$z>6.2$ quasars favors cosmological H ii regions within a significantly
neutral IGM. Here we go a step further, showing that the Ly$\alpha$
transmission at the end of the J0148 proximity zone is consistent with damping
wing absorption, while the no-damping-wing case is strongly ruled out.

The rest of the paper is organized as follows. In Section 2 we describe the
J0148 transmission feature of interest and introduce the statistic we use to
test for the presence of a damping wing. We generate mock quasar proximity
zones in Section 3, and show that generating features similar to the one
towards J0148 requires the proximity zone transmission to be modified by
damping wing absorption. In Section 4 we discuss the physical origin of the
damping wing, arguing that it is more likely to arise from an extended neutral
island rather than a compact absorber. We also describe a scenario in which
the transmission window signals the recent passage of an intergalactic
ionization front. Finally, we summarize our results in Section 5. Throughout
we assume a $\Lambda$CDM cosmology with $(\Omega_{\rm
m},\Omega_{\Lambda},\Omega_{\rm b},h)=(0.308,0.692,0.0482,0.678)$ (Planck
Collaboration XVI, 2014). Distances are quoted in units of comoving $h^{-1}$
Mpc, unless otherwise noted.

##  2 The ULAS J0148$+$0600 Ly$\alpha$ transmission window

Figure 2: J0148 proximity zone. The top and bottom panels show the regions of
the Ly$\beta$ and Ly$\alpha$ forests, respectively, corresponding to the
proximity zone, where the wavelengths have been scaled to align the two panels
vertically in redshift (top axis of Ly$\beta$ panel). The continuum-normalized
flux is plotted as a histogram, with the 1$\sigma$ uncertainty in red shading.
The smooth line in both panels shows the damping profile for a neutral island
of mean density with a length of 7.5 $h^{-1}\,$Mpc, the maximum extent allowed
by the bordering Ly$\alpha$ and Ly$\beta$ transmission. The cyan shaded region
in the lower panel shows the difference between the damping profile and the
flux over a 1000 km s-1 region at the red edge of the damping wing.

The primary feature examined in this paper is an extended window of Ly$\alpha$
transmission centered at $\lambda\simeq 8385$ Å, which falls at the blue end
of the J0148 proximity zone. Here we describe the J0148 spectrum used in this
work, introduce our technique for continuum normalization, and describe the
basic characteristics of the transmission window. We also introduce a
statistic, $\langle\Delta F\rangle$, that will be used to test for the
presence of damping wing absorption.

Our analysis is based on a 10-hour VLT X-Shooter spectrum of J0148 first
presented in Becker et al. (2015) and more recently included in the extended
XQR-30 sample of D’Odorico et al. (2023). We refer the reader to those paper
for details of the data reduction, but note that the J0148 spectrum is
extremely deep, with a continuum-to-noise ratio (based on the continuum
described below) of $\sim$130–220 per 10 km s-1 pixel over the transmission
window and the adjacent Ly$\alpha$ trough. The spectrum is shown in Figure 1.
We adopt a redshift for J0148 of $z_{\rm q}=5.9896$, which is based on [C ii]
158 $\mu$m emission-line measurements (Boseman et al., in prep).

We normalize the spectrum using a continuum that is based on a composite of
lower-redshift objects. The method is described in detail in Appendix A.
Briefly, we identify quasars from the Sloan Digital Sky Survey Data Release 16
(DR16; Lyke et al., 2020) that are good matches to J0148 in terms of their
transmission in two wavelength intervals redward of the Ly$\alpha$ forest. The
intervals cover the red side of the Ly$\alpha$+N v emission line and the C iv
emission line, which are highlighted yellow in Figure 1. Ly$\alpha$ emission
is of particular importance for this work because the transmission window of
interest lies on the blue side of the Ly$\alpha$ line. We also include C iv in
the selection because the strength and systemic velocity offset of C iv
emission are known to be correlated with those of Ly$\alpha$ (e.g., Richards
et al., 2011; Greig et al., 2017; Davies et al., 2018). We select the 200
best-matching DR16 spectra based on mean squared residuals with J0148 within
these intervals. Each spectrum is then corrected for foreground Ly$\alpha$ and
Ly$\beta$ absorption using a statistical approach that includes reduced
absorption within the proximity zone (see Appendix A for details). A composite
is computed by averaging the corrected spectra as a function of their rest-
frame wavelengths (blue histogram in Figure 1). Our final continuum model is a
b-spline fit to this composite (red solid line), with an uncertainty (red
shaded region) that is calculated from the scatter in the individual corrected
spectra, smoothed on 5000 km s-1 scales. At the wavelength of the transmission
window we find a 1$\sigma$ continuum uncertainty of 10%. The uncertainty
estimate is conservative in that part of this scatter will be due to
variations in the amount of Ly$\alpha$ forest absorption between lines of
sight rather than intrinsic variations in the DR16 spectra.

The continuum-normalized J0148 spectrum covering the Ly$\alpha$ proximity zone
and the red end of the giant Ly$\alpha$ trough is displayed in Figure 2. We
also show the corresponding section of the Ly$\beta$ forest. The most
conspicuous feature of the Ly$\alpha$ proximity zone is a strong (maximum
transmission $\simeq$ 0.7) and extended transmission window spanning
$\lambda\simeq 8370$–8400 AA ($z\simeq 5.89$–5.91). Ignoring peculiar
velocities, the wavelength range corresponds to a line-of-sight extent of
$\sim$7 $h^{-1}\,$Mpc. In addition to its amplitude and width, two features of
the transmission window stand out. First, it lies very near to the red end of
the Ly$\alpha$ trough, which is also dark in Ly$\beta$ over $z=5.849$–5.878.
Second, the transmission on the blue side increases relatively smoothly from
zero up to its maximum.

Motivated by the shape of this transmission window and the potential
association of dark Ly$\alpha$+Ly$\beta$ troughs with neutral islands, we
calculate the transmission profile for a neutral island that could reside near
the end of the J0148 proximity zone. We assume that the island is at mean
density, and that it spans the maximum extent allowed by the first
Ly$\alpha$+Ly$\beta$ dark gap blueward of the proximity zone. The Ly$\alpha$
and Ly$\beta$ transmission at the red and blue ends of this gap, respectively,
allow a 7.5 $h^{-1}\,$Mpc island with its near edge roughly 34 $h^{-1}\,$Mpc
from the quasar. No further constraints are imposed by higher-order Lyman-
series transitions. For reference, the central redshift of such an island is
5.862, and the total H i column density is $\log{(N_{\rm H\,I}/{\rm
cm^{-2}})}\simeq 20.5$. We note, however, that the transmission profile will
differ from that of a compact absorber of the same column density due to the
extended nature of the island (Miralda-Escudé, 1998). The Ly$\alpha$ and
Ly$\beta$ transmission profiles for the potential island are shown in Figure
2.

The damping wing of the neutral island roughly traces the blue side of the
broad Ly$\alpha$ transmission window. We also note that the entire Ly$\alpha$
(and Ly$\beta$) proximity zone transmission falls under the damping profile.
The J0148 proximity zone is therefore consistent with the presence of a
damping wing created by an extended neutral island. The central question we
will examine in this paper is whether a damping wing is actually required to
produce the transmission profile of the broad Ly$\alpha$ window, or if such a
feature can arise by chance from resonant absorption alone.

To address this question we will compare the J0148 proximity zone to mock
quasar proximity zones generated with and without damping wing absorption. The
spectra will be evaluated based on their similarity to a damping wing profile
at the blue end of the proximity zone. Specifically, we will measure the mean
difference, $\langle\Delta F\rangle$, between the proximity zone transmission
and the damping wing template for a neutral island similar to the one shown in
Figure 2. In order to focus on the part of the spectrum most impacted by the
damping wing, and to simplify the modeling of the proximity zones, we measure
$\langle\Delta F\rangle$ over a 1000 km s-1 region that starts where the red
side of the damping wing template first exceeds 0.01. The difference between
the J0148 transmission and the damping wing over this region is shaded in
Figure 2. Without any adjustment of the profile parameters we find
$\langle\Delta F\rangle=0.15$ for J0148. Below we test whether this value
indicates that a genuine damping wing is likely to be present.

##  3 Simulated proximity zones

###  3.1 Generating mock proximity zones

Figure 3: Scaling of the H i neutral fraction in the proximity zone models
used in this work. Distances are shown along the top axis, with the bottom
axis showing the corresponding observed wavelength for a quasar at $z_{\rm
q}=5.9896$. The solid lines show the factors by which the neutral fraction is
multiplied as a function of distance from the quasar for four different quasar
luminosities. $L_{0}$, corresponds to the nominal ionizing UV luminosity of
J0148. The neutral fraction scaling is calculated for 0.5, 1.0, 2.0, and 4.0
times this luminosity (see text for details). The vertical dashed line
corresponds to the edge of transmitted Ly$\alpha$ flux in the J0148 proximity
zone. Note that these scalings are only applied to the ionized gas within the
proximity zone; trends out to larger distances are shown for illustration. The
scalings include the effect of ionizing opacity within the proximity zone, as
discussed in Section 3.1. For reference, the dotted line shows the zero-
opacity case for $L_{0}$.

Mock proximity zones are generated using one-dimensional lines of sight drawn
from the $z=6.0$ output of the Sherwood hydrodynamical simulation suite
(Bolton et al., 2017). We use 5,000 lines of sight from the 80-2048 run, which
has a box size of 80 $h^{-1}\,$Mpc and 20483 gas particles with a mass of
$7.97\times 10^{5}~{}h^{-1}~{}{\rm M_{\odot}}$. The lines of sight are not
drawn through the massive halos expected to host bright quasars; however, this
is unlikely to affect our analysis, which is focused on features that are
offset by more than $\sim$30 $h^{-1}\,$Mpc from the quasar position. By
comparison, the virial radius for a $10^{12}~{}h^{-1}~{}{M_{\odot}}$ halo will
be $\sim$0.3 $h^{-1}\,$Mpc at $z=6$ (Barkana & Loeb, 2001). We rescale the gas
temperatures using a power-law temperature-density relation of the form
$T(\rho)=T_{0}(\rho/\langle\rho\rangle)^{\gamma-1}$, adopting $T_{0}=12000$ K
and $\gamma=1.2$, broadly consistent with measurements from Gaikwad et al.
(2020) over $5.3<z<5.9$. Baseline Ly$\alpha$ optical depths are set by
globally rescaling the hydrogen neutral fractions to produce a mean Ly$\alpha$
transmission of 0.01, similar to the measurements of Bosman (2021) at $z=5.9$
(i.e., near the end of the J0148 proximity zone). For each of the samples
described below we then generate two proximity zones per line of sight, one in
each direction, for a total of 10,000 proximity zones per sample.

The simulated proximity zones consist of an extended optically thin photo-
ionized region that terminates at an optically thick absorber, with various
absorbers considered in the following sections. For the optically thin
portions we use the proximity zone model described in Becker et al. (2021) and
Zhu et al. (2023). Briefly, hydrogen neutral fractions along the line of sight
are scaled according to $x_{\rm H\,I}\propto(\Gamma/\Gamma_{\rm bg})^{-1}$,
where the local photoionization rate, $\Gamma=\Gamma_{\rm bg}+\Gamma_{\rm
q}(r)$, combines contributions from the diffuse UV background and the ionizing
flux from the quasar. The radial dependence of the quasar flux includes both a
$1/r^{2}$ geometric dilution and attenuation due to the ionizing opacity
within the proximity zone. The opacity scales with the local photoionization
rate as $\kappa_{912}=\kappa_{912}^{\rm bg}(\Gamma/\Gamma_{\rm bg})^{\xi}$,
and therefore increases with distance from the quasar following equation (8)
in Becker et al. (2021). We adopt a background opacity value of
$\log{(\kappa_{\rm bg}/({\rm cm}^{-1})}=-24.5$. This corresponds to a mean
free path of ionizing photons of 1.0 proper Mpc, consistent with recent
measurements from Becker et al. (2021) and Zhu et al. (2023). For $\xi$ we
adopt 0.67 based theoretical models of the ionizing sinks (Furlanetto & Oh,
2005; McQuinn et al., 2011, see discussions in Davies & Furlanetto 2016;
D’Aloisio et al. 2018, 2020; Becker et al. 2021). For the background
photoionization rate we adopt $\Gamma_{\rm bg}=0.17\times 10^{-12}~{}{\rm
s^{-1}}$, consistent with constraints from Gaikwad et al. (2023) at $z=5.9$,
albeit slightly higher than their nominal value. For the simulated quasars we
adopt a nominal absolute UV magnitude at rest-frame 1450 Å of
$M_{1450}=-27.4$, similar to J0148 (Bañados et al., 2016). The UV spectrum is
modeled as a double power law of the form $L_{\nu}\propto\nu^{-\alpha}$, with
$\alpha=0.6$ betwen 1450 Å and 912 Å, and 1.5 shortward of 912 Å (see Lusso et
al., 2015; Becker et al., 2021, and references therein). This agrees well with
our composite-based continuum over the non-ionizing UV part of the spectrum
(Figure 1).

In principle there are a variety of factors that will impact the neutral
fraction within the proximity zone. These include the full quasar spectral
energy distribution, the quasar lifetime and/or variability, variations in
ionizing opacity along the line of sight not captured by the model above,
fluctuations in the ionizing background, and gas temperature. The neutral
fraction at the edge of the proximity zone will also depend on the distance
from the quasar. In practice, however, we will focus only on a narrow region
at the end of the proximity zone, for which the important quantity is the
local neutral fraction. We therefore sample a range of neutral fractions by
simply varying an “effective” quasar luminosity, $L_{\rm q}$, but note that
these variations may also result from other factors, which we examine further
in Section 4.2. The scaling of $x_{\rm H\,I}$ for optically thin regions
within the proximity zones is shown in Figure 3. We refer to the scaling given
by the nominal parameters above as $L_{\rm q}=L_{0}$, and compute proximity
zones for $L_{\rm q}=0.5L_{0}$, $L_{0}$, 2.0$L_{0}$, and 4.0$L_{0}$. At the
distance to the edge of the J0148 proximity zone, $\sim$33 $h^{-1}\,$Mpc, the
neutral fraction scaling ranges from 0.60 (i.e., 60% of the value in the
absence of the proximity effect) for $0.5L_{0}$, to 0.05 for $4L_{0}$. As
shown below, this range is sufficient to evaluate the role of ionization in
the proximity zone models we consider. We note that the scaling of the neutral
fraction with quasar luminosity is amplified by the impact of $L_{\rm q}$ on
the ionizing opacity, following the discussion above. The total optical depth
within the proximity zone is $\tau\simeq 1.66$, 0.98, 0.57, and 0.34 for
$L_{\rm q}=0.5L_{0}$, $L_{0}$, $2L_{0}$, and $4L_{0}$, respectively. For
reference, we plot the zero-opacity case for $L_{0}$ as a dotted line in
Figure 3.

Finally, we apply a random continuum error to each simulated line of sight.
The continuum error is modeled as a constant factor drawn from a normal
distribution with a mean of 1.0 and a standard deviation of 0.1, similar to
the fractional continuum uncertainty of J0148 over the transmission window at
$\sim$8385 Å (Figure 2). We note that this does not include a possible
wavelength dependence in the continuum; however, the analysis below focuses on
a narrow wavelength range, for which a uniform continuum rescaling should be
adequate. We convolve each spectrum with a Gaussian kernel with full width at
half maximum of 23 km s-1, equal to the measured X-Shooter resolution for
J0148 (D’Odorico et al., 2023). We also rebin the spectra to a pixel scale
similar to the J0148 spectrum, although the results should not depend
sensitively on the choice of binning.

###  3.2 $\langle\Delta F\rangle$ measurement method

We evaluate the proximity models below based on how frequently they produce
$\langle\Delta F\rangle$ values equal to or less than the J0148 value. As with
our analysis of J0148, we use a damping profile template calculated from a 7.5
$h^{-1}\,$Mpc neutral island at mean density. $\langle\Delta F\rangle$ is
calculated as the difference between the proximity zone transmission and this
template over a 1000 km s-1 region starting where the red wing of the damping
profile first exceeds 1%. Conservatively, we minimize $\langle\Delta F\rangle$
for each line of sight by placing the island at the minimum distance from the
quasar where the transmission within this region does not exceed the damping
profile.

We note that we allow the proximity zone transmission to exceed the damping
profile at wavelengths redder than the $\langle\Delta F\rangle$ region, for
multiple reasons. First, there may be wavelength-dependent errors in the
continuum that are not captured by our approach of rescaling the transmission
in each simulated proximity zone by a constant. Second, there may be
variations in the $x_{\rm H\,I}$ scaling that are not well modeled by the
radial profiles plotted in Figure 3. Third, the simulated lines of sight are
drawn from random locations with the simulation box and do not capture the
enhanced densities that may be present in the vicinity of a bright quasar,
although this would mainly impact absorption close to the quasar redshift (see
Section 3.1). Restricting our analysis to a narrow region at the edge of the
proximity allows us to simplify our analysis by avoiding these more nuanced
aspects of the modeling. In the future, however, a more robust modeling of the
transmission and continuum could allow the full proximity zone to be used to
constrain the presence of a damping wing (e.g., Davies et al., 2018). We note
here that for our nominal continuum estimate, the J0148 flux remains below the
damping wing profile over the entire proximity zone (Figure 2).

Figure 4: Transmission profiles for mean-density neutral islands of length 2.5
(dot), 5.0 (dash-double-dot), 7.5 (solid), 10 (dash-dot) and 15 (dash)
$h^{-1}\,$Mpc. The near edge for all islands is 34 $h^{-1}\,$Mpc from a quasar
at $z_{\rm q}=5.9896$. The shaded region shows the 1000 km s-1 interval used
to calculate $\langle\Delta F\rangle$. We use the 7.5 $h^{-1}\,$Mpc template
for our primary analysis, which is motivated by the Ly$\alpha$ and Ly$\beta$
transmission near the end of the J0148 proximity zone. Our results do not
substantially change, however, when other templates are used.

It is worth considering whether the details of our $\langle\Delta F\rangle$
measurement may impact the results below. For example, we have chosen to
measure $\langle\Delta F\rangle$ over a 1000 km s-1 region partly because this
is the extent to which the J0148 transmission window appears to follow the
Ly$\alpha$ damping wing in Figure 2. Our main result, that we exclude models
without damping wings (see Section 3.3 below), actually becomes somewhat
stronger if we measure $\langle\Delta F\rangle$ out to 1200 km s-1, and only
mildly weaker at 1500 km s-1, which roughly corresponds to the red edge of the
J0148 transmission window. Predictably, $\langle\Delta F\rangle$ becomes less
able to discriminate between models if we extend the region further to 2000 km
s-1 due to the contribution from strong absorption over 8401–8417 Å. It is
expected that resonant Ly$\alpha$ forest absorption will become increasingly
dominant over damping wing absorption at large velocities, making
$\langle\Delta F\rangle$ a less sensitive statistic over wider windows. For
this line of sight, therefore, we have chosen 1000 km s-1 in order to provide
a reasonable baseline to test for a damping wing signal while limiting the
contamination from unrelated absorption.

We also note that the results below do not depend sensitively on the damping
wing template used to measure $\langle\Delta F\rangle$. As described above,
our choice of a template generated from a 7.5 $h^{-1}\,$Mpc neutral region at
mean density is motivated by the spacing between Ly$\alpha$ and Ly$\beta$
transmission peaks in the J0148 spectrum. For comparison, in Figure 4 we show
damping profiles for mean-density neutral islands with lengths ranging from
2.5 $h^{-1}\,$Mpc to 15 $h^{-1}\,$Mpc, which are aligned at their edges
nearest to the quasar. We note that contiguous islands longer than $\sim$7.5
$h^{-1}\,$Mpc are not allowed towards J0148 by the Ly$\beta$ transmission
peaks in Figure 2; however, there could be multiple islands separated by small
gaps. The extent of the absorption trough depends strongly on the size of the
island, but $\langle\Delta F\rangle$ depends only on the shape of the red
wing, for which the variation is more modest. We repeated the tests below
using these other profiles and found that the main results do not change
substantially. We also tried marginalizing over the template island length by
taking the minimum $\langle\Delta F\rangle$ from among the profiles shown in
Figure 4 and again found no substantial change.

###  3.3 No damping wing absorption

Figure 5: Examples of simulated proximity zones ending at Lyman limit systems
with column density $\log{(N_{\rm H\,I}/{\rm cm^{-2}})}=18.0$. The transmitted
flux in each panel has been rescaled by a random amount to emulate continuum
errors (see text for details). Dashed lines show the LLS absorption profile,
which has been rescaled by the same amount. The smooth solid line shows the
absorption profile for a 7.5 $h^{-1}\,$Mpc neutral island at mean density. The
shaded region shows the difference between this profile and the transmitted
flux over a 1000 km s-1 region. For each line of sight, the damping profile is
shifted in wavelength such that it minimizes $\langle\Delta F\rangle$ while
not allowing the transmission to exceed the damping profile within the widow.
The transmission is allowed to exceed the damping profile at longer
wavelengths, however (see text for details). There are three lines of sight
plotted for each quasar luminosity scaling in Figure 3, where the scaling is
indicated in the upper left corner of each panel. The lines of sight
correspond to the 5th, 50th, and 95th percentiles (top to bottom) in the
$\langle\Delta F\rangle$ distribution for each $L_{\rm q}$. The value of
$\langle\Delta F\rangle$ is given in parentheses. Figure 6: Cumulative
probability distributions of $\langle\Delta F\rangle$ for simulated lines of
sight ending in Lyman limit systems. The solid lines show $P(<\langle\Delta
F\rangle)$ for the four quasar luminosity scalings, with colors correspond to
Figure 3. For reference, the left-to-right ordering of these lines at
$P(<\langle\Delta F\rangle)=0.5$ is $L_{\rm q}=2L_{0},L_{0},4L_{0},0.5L_{0}$.
The vertical dashed line shows $\langle\Delta F\rangle$ for J0148.

We first test whether features similar to the J0148 transmission window can
occur in the absence of damping wing absorption. For this test, the simulated
proximity zones are truncated at a Lyman limit system (LLS) with neutral
hydrogen column density $\log{(N_{\rm H\,I}/{\rm cm^{-2}})}=18.0$. This
provides sufficient optical depth to truncate the quasar’s ionizing flux but
will not produce a significant damping wing. We assume that the LLS should
reside in a density peak that is at least three times the mean density, but
find that either increasing this to ten times the mean or removing the density
requirement entirely would have little impact on our results. For each line of
sight we randomly select a density peak, and then shift the periodic line of
sight such that this peak lies at $z=5.878$ when peculiar velocities are taken
into account, assuming the quasar is at $z=5.9896$. The LLS redshift is chosen
to match the maximum extent of the J0148 proximity zone transmission, and for
zero peculiar velocity corresponds to a distance from the quasar of 33
$h^{-1}\,$Mpc. We set $\Gamma_{\rm q}=0$ beyond the LLS. We also produce a
complete absorption trough outside of the proximity zone by arbitrarily
increasing $x_{\rm H\,I}$. For this test we remain agnostic about the physical
origin of the trough, although one may arise from an ionized IGM if the UV
background outside of the proximity zone is highly suppressed. Here a trough
is introduced simply to maximize the similarity of a simulated line of sight
to J0148.

We note that a few lines of sight randomly intersect high column density
absorbers that are native to the simulation box. These absorbers naturally
produce damping wing absorption, which can impact the appearance of the
proximity zone. Because we wish to determine whether a feature similar to
J0148 can be produced in the absence of any damping wing absorption, we
suppress the Lorentzian wings by using a Gaussian absorption profile for
ionized gas (although not for the LLS). The impact of damping wing absorption
from higher-column density absorbers is examined below.

Example proximity zones ending in LLSs are plotted in Figure 5. We plot three
lines of sight for each $L_{\rm q}$ scaling, which are chosen to lie at the
5th, 50th, and 95th percentiles in $\langle\Delta F\rangle$ for that scaling
(bottom to top for each $L_{\rm q}$). Similar to Figure 2, the template
damping profile is plotted in grey and the shaded region shows the 1000 km s-1
region over which $\langle\Delta F\rangle$ is measured. The full
$\langle\Delta F\rangle$ distributions are plotted in Figure 6. The smallest
$\langle\Delta F\rangle$ values tend to be produced by $L_{\rm q}=L_{0}$ and
$2L_{0}$. By comparison, $0.5L_{0}$ lines of sight tend to exhibit only weak
transmission peaks at the edge of the proximity zone, while the strong
transmission for $4L_{0}$ tends to prohibit the $\langle\Delta F\rangle$
region of the damping profile from extending over the transmission peaks. In
all cases, however, the the LLS models fail to generate $\langle\Delta
F\rangle$ values equal to or less than the J0148 value in more than 0.01% of
cases. We therefore conclude that lines of sight without significant damping
wings are highly unlikely to produce the kind of smooth Ly$\alpha$
transmission feature seen at the blue end of the J0148 proximity zone.

###  3.4 Compact damped Ly$\alpha$ absorbers

Figure 7: Similar to Figure 5, but for simulated lines of sight ending in
compact DLAs with column density $\log{(N_{\rm H\,I}/{\rm cm^{-2}})}=20.3$.
The dashed line in each panel shows the DLA damping profile, which has been
rescaled by the same amount as the transmitted flux. The smooth solid line is
a template profile for a mean-density 7.5 $h^{-1}\,$Mpc neutral island that
has been shifted to produce the minimum allowable $\langle\Delta F\rangle$
over the 1000 km s-1 shaded region at the red edge of the template damping
wing. Figure 8: Similar to Figure 6, but for simulated proximity zones ending
in compact DLAs with column density $\log{(N_{\rm H\,I}/{\rm cm^{-2}})}=20.3$.
Colors correspond to the quasar luminosity scalings in Figure 3. The left-to-
right ordering of the solid lines is $L_{\rm q}=4L_{0},2L_{0},L_{0},0.5L_{0}$.
The vertical dashed line shows $\langle\Delta F\rangle$ for J0148.

Next we evaluate whether the J0148 feature could be created via damping wing
absorption from a compact source such as a damped Ly$\alpha$ system (DLA)
residing in a galaxy halo. Here we reuse the simulated lines of sight
generated in Section 3.3, but add higher-column density absorbers to the end
of the proximity zone ($z=5.878$). The H i column density is varied from
$\log{(N_{\rm H\,I}/{\rm cm^{-2}})}=19.0$ to 20.5 in increments of 0.1 dex,
thus spanning a broad range of DLAs and sub-DLAs. The minimum $\langle\Delta
F\rangle$ value is then found for each combination of $N_{\rm H\,I}$ and
$L_{\rm q}$. In this section we focus on $\log{(N_{\rm H\,I}/{\rm
cm^{-2}})}=20.3$ because it is the canonical DLA value and close the column
density where the probability of producing $\langle\Delta F\rangle$ values at
or below that of J0148 tends to be greatest. The full range of column
densities is examined in Appendix C.

Example lines of sight ending in compact DLAs with column density
$\log{(N_{\rm H\,I}/{\rm cm^{-2}})}=20.3$ are shown in Figure 7. Similar to
Figure 5, lines of sight at the 5th, 50th, and 95th percentiles in
$\langle\Delta F\rangle$ are shown for each $L_{\rm q}$ scaling. The full
$\langle\Delta F\rangle$ distributions are plotted in Figure 8. We find that a
significant percentage of lines of sight with $\langle\Delta F\rangle$ values
equal or less than the J0148 value are produced for the higher $L_{\rm q}$
values (3% for 2$L_{\rm q}$, 18% for 4$L_{\rm q}$). As seen in Figure 7, in
cases where $\langle\Delta F\rangle$ falls near or below the J0148 value of
0.15, the red wing of the neutral island damping wing template (solid grey
line) tends to closely follow the red wing of the actual DLA profile (dashed
line), even though the template profile is considerably more extended on the
blue side. In these cases, the DLA profile modifies a span of the proximity
zone that otherwise has high transmission within the $\langle\Delta F\rangle$
region, as seen by the nearness of the tops of the transmission peaks to the
DLA profile. We note that these lines of sight also exhibit considerably more
transmission over the full proximity zone than is seen in J0148 (e.g.,
comparing the bottom panel of Figure 2 to the bottom two panels of Figure 7
over $8400~{}{\mbox{\AA}}\lesssim\lambda\lesssim 8500~{}{\mbox{\AA}}$). As
noted above, however, in this work we conservatively consider only the 1000 km
s-1 region within which $\langle\Delta F\rangle$ is measured (although see
Section 4.2.1). We thus conclude that damping wings from compact, high-column
density absorbers can plausibly produce features at the edges of proximity
zones similar to the Ly$\alpha$ feature in the J0148 spectrum.

###  3.5 Neutral islands

Figure 9: Simliar to Figure 5, but for simulated lines of sight ending in 7.5
$h^{-1}\,$Mpc neutral islands. The dashed line in each panel shows the actual
damping profile for the neutral island, which has been rescaled by the same
amount as the transmitted flux. The smooth solid line is a template profile
for a mean-density 7.5 $h^{-1}\,$Mpc neutral island that has been shifted to
produce the minimum allowable $\langle\Delta F\rangle$ over the 1000 km s-1
shaded region at the red edge of the template damping wing. Figure 10: Similar
to Figure 6, but for simulated proximity zones ending in neutral islands.
Solid lines are for the fiducial island length of 7.5 $h^{-1}\,$Mpc. Other
line styles are for lengths 2.5 (dotted), 5.0 (dash-double-dotted), 10 (dash-
dotted), and 15 (dashed) $h^{-1}\,$Mpc. Colors correspond to the quasar
luminosity scalings in Figure 3. The left-to-right ordering of the line groups
is $L_{\rm q}=4L_{0},2L_{0},L_{0},0.5L_{0}$. The vertical dashed line shows
$\langle\Delta F\rangle$ for J0148.

Finally, we consider whether the smooth profile at the edge of the J0148
proximity zone could be caused by the damping wing from an extended neutral
island. Here we generate artificial proximity zones as described in Section
3.3, but set $x_{\rm H\,I}=1$ over a 7.5 $h^{-1}\,$Mpc interval starting 34
$h^{-1}\,$Mpc from the quasar. The choices of length and distance are
motivated by the largest mean-density neutral island allowed by the observed
Ly$\alpha$ and Ly$\beta$ transmission at the end of the J0148 proximity zone
(see Section 2). We note that in terms of $\langle\Delta F\rangle$, the
distance will be largely degenerate with the $L_{\rm q}$ scaling.

Example simulated proximity zones ending in neutral islands are shown in
Figure 9. Similar to Figure 5, lines of sight at the 5th, 50th, and 95th
percentiles in $\langle\Delta F\rangle$ are shown for each $L_{\rm q}$
scaling. The full $\langle\Delta F\rangle$ distributions are plotted in Figure
10. The results are broadly similar to those for DLAs in that the higher
$L_{\rm q}$ values can produce significant percentages of lines of sight with
$\langle\Delta F\rangle$ values at or below the J0148 value (1% for $2L_{0}$,
10% for $4L_{0}$). In these cases (e.g, bottom two panels of Figure 9), the
template damping profile used to measure $\langle\Delta F\rangle$ (solid line)
is reasonably well aligned with the actual damping wing produced by the
neutral island along the line of sight (dashed line). The nearness of the tops
of the transmission peaks to the dashed line again suggests that the
transmission would otherwise be high within the $\langle\Delta F\rangle$
region, a point we return to below. Here we simply conclude that the damping
wing profiles created by large neutral islands can plausibly produce features
with $\langle\Delta F\rangle$ values similar to J0148.

We can also test whether the consistency with J0148 depends sensitively on the
length of the neutral island. To check this, we generate mock proximity zones
ending in islands with lengths 2.5, 5.0, 10, and 15 $h^{-1}\,$Mpc. In each
case the red end of the island is 34 $h^{-1}\,$Mpc from the quasar. The
$\langle\Delta F\rangle$ distributions are shown in Figure 10, and we note
that we are still calculating $\langle\Delta F\rangle$ using the 7.5
$h^{-1}\,$Mpc mean-density template. Although there is some sensitivity to the
length of the island, the changes with respect to our fiducial length of 7.5
$h^{-1}\,$Mpc are modest. We therefore conclude that while the damping wing
from a neutral island can produce a $\langle\Delta F\rangle$ value similar to
the J0148 value, the agreement is not strongly dependent on the island length.

##  4 Discussion: origin of the J0148 damping wing

###  4.1 Compact absorber vs. neutral island

Figure 11: Illustration of limits on metal lines at redshifts over the
Ly$\alpha$+Ly$\beta$ trough at the end of the J0148 proximity zone. Each panel
shows the normalized flux over the wavelengths corresponding to the indicated
ion over absorption redshifts $5.844<z_{\rm abs}<5.883$. Unrelated absorption
lines are plotted in light grey. Shaded regions at the edges contain
significant Ly$\alpha$ and/or Ly$\beta$ transmission. The blue line in the top
three panels shows the expected absorption for a low-ionization absorber with
an O i column density of $\log{(N_{\rm O\,I}/{\rm cm^{-2}})}=13.2$ and
relative abundances similar to known $z\simeq 6$ absorbers (see text for
details). Green lines in the next two panels show the C iv transmission for
$\log{(N_{\rm C\,IV}/{\rm cm^{-2}})}=12.8$. Note the vertical scale in the
upper five panels. Red shaded regions in the lower two panels show the
1$\sigma$ flux uncertainty.

The above analysis demonstrates that a smooth transmission feature at the end
of proximity zone similar to the J0148 Ly$\alpha$ transmission window strongly
requires a damping wing. The damping wing can in principle be created either
by a compact absorber, such as a DLA, or an extended neutral island. Here we
argue, however, that while both scenarios are allowed in terms of their
$\langle\Delta F\rangle$ statistics, a neutral island is favored based on
other factors.

There are at least three reasons to favor an extended neutral island. The
first is the lack of any corresponding metal absorption lines. As noted by
Becker et al. (2019) and Davies et al. (2023a), the J0148 spectrum exhibits no
metal absorption systems at $z>5.5$. In particular, no metal absorption is
seen in either low-ionization (e.g., O i, C ii, Si ii, Mg ii) or high-
ionization (C iv, Si iv) lines over the redshifts of the Ly$\alpha$+Ly$\beta$
trough at the end of the proximity zone ($z\simeq 5.85$–5.88; Figure 11). The
deep X-Shooter data excludes systems with $\log{(N_{\rm O\,I}/{\rm
cm^{-2}})}\gtrsim 13.2$ for a Doppler $b$-parameter of 10 km s-1. In Appendix
C, we show that a compact absorber would need an H i column density of at
least $\sim$$10^{19.5}~{}{\rm cm^{-2}}$ in order for the probability of
generating a $\langle\Delta F\rangle$ value equal to or less than the J0148
values to be at least 1%, even for the highest $L_{\rm q}$ scaling we
consider. Assuming that oxygen is in charge exchange equilibrium with
hydrogen, these limits on O i and H i corresponds to a maximum oxygen
metallicity of ${\rm[O/H]}\lesssim-3.0$,222Using photospheric solar abundances
from Asplund et al. (2009), and assuming no dust depletion. which is
comparable to the most metal-poor DLA known (Cooke et al., 2017). We show the
expected line profiles for such an absorber in Figure 11, where the column
densities for Si ii and C ii are scaled with respect to O i based on the
relative abundances measured by Becker et al. (2012) for low-ionization
absorbers near $z=6$ (see also Sodini et al., in prep). The non-detection of
Si ii $\lambda$1260 would increase the metallicity constraint to
${\rm[O/H]}\lesssim-3.1$ for the same abundance ratios. The presence of such
an ultra-low-metallicity, low-ionization absorber is clearly ruled out at
these redshifts. We also show the expected C iv absorption for $\log{(N_{\rm
C\,IV}/{\rm cm^{-2}})}=12.8$ and $b=30$ km s-1, which is also excluded by the
data. If a compact absorber is creating the damping wing observed at the edge
of the J0148 proximity zone then it must be extremely metal poor. We note,
however, that some DLAs at these redshifts may have metallicities near
${\rm[O/H]}\simeq-3$ (e.g., D’Odorico et al., 2018; Bañados et al., 2019), so
the lack of metal lines on its own does not fully rule out a compact absorber.

A second reason to favor a neutral island is that the evidence of a damping
wing is found adjacent to an extremely long and opaque giant Ly$\alpha$
trough. In our simulated proximity zones we set the Ly$\alpha$ transmission to
zero blueward of the proximity zone by artificially increasing the H i neutral
fraction. In more physically motivated models, however, it has been shown that
such long troughs are most likely to occur when a line of sight intersects one
or more neutral islands (e.g., Kulkarni et al., 2019; Keating et al., 2020a;
Nasir & D’Aloisio, 2020). In a companion paper, Zhu et al. (in prep) show that
long dark troughs over $5\lesssim z\lesssim 6$ are statistically associated
with damping wing absorption, which is consistent with a model in which the
troughs arise from neutral islands.

Finally, a deep search for [O iii]-emitting galaxies with the James Webb Space
Telescope as part of the EIGER survey reveals no sources within the potential
DLA redshift window of $5.85\lesssim z\lesssim 5.88$ (Eilers et al., 2024).
Based on the lack of metal lines, the presence of an extreme Ly$\alpha$
trough, and the lack of any directly detected galaxy, we therefore conclude
that the damping wing over the J0148 Ly$\alpha$ transmission window is more
likely to arise from an extended neutral island than a compact absorber.

###  4.2 Physical characteristics of the high-transmission window

Figure 12: Simliar to Figure 9, but for simulated lines of sight ending in 7.5
$h^{-1}\,$Mpc neutral islands that have $\langle\Delta F\rangle$ values
similar to the J0148 value. Three lines of sight are shown for quasar
luminosity scalings $L_{\rm q}=L_{0}$ (top), $2L_{0}$ (middle), and $4L_{0}$
(bottom). Lines of sight similar to J0148 are extremely rare in these simple
models for $L_{\rm q}=L_{0}$, while they are more common for $2L_{0}$ and
$4L_{0}$ (Figure 10).

We now turn to examining the physical conditions that may give rise to
transmission features such the J0148 window. For illustration, in Figure 12 we
plot a sample of simulated lines of sight with $\langle\Delta F\rangle$
similar to the J0148 value for $L_{\rm q}=L_{0}$, $2L_{0}$, $4L_{0}$. In many
cases the tops of the transmission peaks within the $\langle\Delta F\rangle$
region nearly reach the damping wing envelope produced by the neutral island.
(This is somewhat less true for $L_{\rm q}=L_{0}$ than for $2L_{0}$ and
$4L_{0}$, but cases similar to J0148 for $L_{\rm q}=L_{0}$ are extremely rare;
see Figure 10.) This suggests that the regions adjacent to these neutral
islands are highly transmissive in Ly$\alpha$, which can be achieved via high
ionization rates, low densities, and/or elevated gas temperatures. We consider
each of these factors in the following sections

####  4.2.1 Ionization rate

It is clear from Figures 8 and 10 that the probability of producing a line of
sight with $\langle\Delta F\rangle$ equal to or less than the J0148 value
increases towards higher quasar luminosities. The higher ionization rates
produce lower neutral fractions (Figure 3), and hence more transmission
against which the damping wing profile becomes apparent. In our models the
ionization rate within the proximity zone is determined by a combination of
the quasar luminosity and ionizing opacity, and for simplicity we only vary
$L_{\rm q}$ (Figure 3). While this approach may be adequate to bracket the
range of possible ionization rates, invoking a high quasar luminosity to
produce a window of high transmission may be inconsistent with other
characteristics of the J0148 proximity zone. Notably, if we select the
simulated lines of sight with $\langle\Delta F\rangle\leq 0.2$ (cf. the J0148
value of 0.15) and compute their mean transmission in a 2000 km s-1 region
from 8400 Å (typically just reward of the $\langle\Delta F\rangle$ region) to
8456 Å (10 $h^{-1}\,$Mpc from the quasar), all of the lines of sight have
higher transmission than J0148, even when random continuum errors are
included. Keating et al. (2015) find that the properties of $z\sim 6$ quasar
proximity zones only weakly depend on the host halo mass. The fact that our
lines of sight do not start at massive halos is therefore unlikely to account
for this difference in transmission (see also Section 3.1). We also note that
our composite-based quasar continuum (Figure 1) is highly consistent with the
UV spectrum we adopt to compute the ionizing emissivity (Section 3.1). This
may make a factor of two or more departure from the nominal ionizing
luminosity unlikely, although quasar variability could be playing a role
(e.g., Davies et al., 2020; Satyavolu et al., 2023).

This inconsistency redward of the $\langle\Delta F\rangle$ region,
particularly with the $L_{\rm q}=2L_{0}$ and $4L_{0}$ models, which produce
the largest transmission, suggests that there may be other factors driving the
high transmission of the J0148 window. One possibility is that there are local
sources contributing to the ionization rate near the end of the J0148
proximity zone. Indeed, Eilers et al. (2024) identify a group of [O
iii]-emitting galaxies at $z\simeq 5.92$, coincident with the strong
Ly$\alpha$ absorption just redward of the J0148 window. In order for the
probability of producing a $\langle\Delta F\rangle$ value consistent with
J0148 to exceed 1%, similar to the $2L_{0}$ case in Figure 10, such sources
would need to roughly triple the ionization rate of our baseline $L_{\rm
q}=L_{0}$ case,333The difference in the ionization rate between the $L_{0}$
and $2L_{0}$ cases includes an additional factor of 1.5 due to the difference
in ionizing opacity within the proximity zone; see Section 3.1. which is
already seven times the background rate at the center of the transmission
window at $z=5.897$ ($\lambda_{\rm obs}=8385$ Å; Figure 3). Other ways of
boosting the J0148 window transmission include low densities and/or high gas
temperatures, which we examine below.

####  4.2.2 Density

Figure 13: Median line-of-sight overdensities for lines of sight ending in 7.5
$h^{-1}\,$Mpc neutral islands with $\langle\Delta F\rangle$ values within
$\pm$0.05 of the J0148 value. Colors correspond to the quasar luminosity
scaling values in Figure 3, although only $L_{\rm q}=L_{0}$ (green), $2L_{0}$
(orange), and $4L_{0}$ (red) are shown here. These lines of sight tend to
contain an underdensity near the end of the quasar proximity zone, with a
median depth that increases with decreasing quasar luminosity.

The role of density can be examined using the line-of-sight gas densities for
our mock proximity zones. We select lines of sight with neutral islands that
have $\langle\Delta F\rangle$ values within $\pm$0.05 of the J0148 value
(i.e., 0.10 to 0.20) for $L_{\rm q}=L_{0}$, $2L_{0}$, $4L_{0}$. We then
compute the median density along these lines of sight separately for each
$L_{\rm q}$. The results are shown in Figure 13. The median profiles show
clear underdensities from $\sim$24 to $\sim$32 $h^{-1}\,$Mpc from the quasar.
This closely corresponds to the position of the J0148 transmission window, and
is roughly adjacent to the near edge of the simulated neutral islands. The
depth of the median underdensity increases with decreasing $L_{\rm q}$,
suggesting that the role of density is especially important when the local
ionization rate is lower. The underdensities required for the $L_{\rm
q}=L_{0}$ case are extremely rare, however, which is why only $\sim$0.1% of
these lines of sight have $\langle\Delta F\rangle$ values equal to or less
than the J0148 value of 0.15, and only $\sim$1% have $\langle\Delta
F\rangle\leq 0.2$.

####  4.2.3 Temperature

Figure 14: Cumulative probability distributions of $\langle\Delta F\rangle$
for simulated lines of sight ending in neutral islands. The solid lines are
the same as those plotted in Figure 10. The dashed lines are for the same
lines of sight with the gas temperature at mean density increased from 12,000
K to 30,000 K. The left-to-right ordering of both line styles is $L_{\rm
q}=4L_{0}$ (red), $2L_{0}$ (orange), $L_{0}$ (green), $0.5L_{0}$ (blue). In
terms of the change in $\langle\Delta F\rangle$, increasing the gas
temperature by this amount is roughly equivalent to doubling the quasar
luminosity.

Finally, we consider the role of gas temperature in producing transmission
features similar to the J0148 window. In the preceding models we adopted a
temperature-density relation with a temperature at mean density of 12,000 K
(see Section 3.1). The temperature may be higher locally, however, if the gas
was recently reionized. We tested the impact of elevated temperatures on the
$\langle\Delta F\rangle$ distribution by recalculating the optical depths for
simulated lines of sight with neutral islands after increasing the gas
temperature to a uniform value. The choice of temperature is motivated by
radiative transfer calculations from D’Aloisio et al. (2019), which predict
post-ionization temperatures as a function of the ionization front velocity
and the spectral index of the incident ionizing radiation. Calculating the
ionization front velocity as $v_{\rm ion}=\dot{n}_{\rm ion}/n_{\rm H\,I}$,
where $\dot{n}_{\rm ion}$ is the incident flux of ionizing photons and $n_{\rm
H\,I}$ is the neutral hydrogen number density, a front driven by a quasar with
$L_{\rm q}=L_{0}$ would be moving at $\sim$25,000 km s-1 through mean-density
gas once it reached the edge of the J0148 proximity zone. This includes an
attenuation of the quasar flux by $\sim$$1/e$ due to opacity within the
proximity zone predicted by the model described in Section 3.1. D’Aloisio et
al. (2019) find a post-I-front temperature of $\sim$30,000 K at this velocity
for our adopted ionizing spectral index of $\alpha=1.5$. The incident spectrum
may be somewhat harder due to spectral filtering within the proximity zone,
though the D’Aloisio et al. (2019) relations suggest that this would not
strongly impact the temperature. The temperature could be as high as
$\sim$40,000 K if helium is doubly ionized by the quasar (Bolton et al., 2012;
Davies et al., 2016). For simplicity, we increase the temperature to 30,000 K
along the entire line of sight, but again note that we are evaluating
$\langle\Delta F\rangle$ only within a narrow region at the end of the
proximity zone.

The resulting $\langle\Delta F\rangle$ distributions are plotted in Figure 14.
Increasing the temperature to 30,000 K substantially increases the number of
lines of sight with $\langle\Delta F\rangle$ values equal to or less than
J0148, roughly comparable to the effect of doubling the quasar luminosity. We
note that most of the change in the $\langle F\rangle$ distribution comes from
the lower neutral fractions associated with suppressed recombination rates
rather than from increased thermal broadening.

###  4.3 A recently reionized region?

In the preceding sections we have shown that a smooth Ly$\alpha$ transmission
window such as the one towards J0148 is more likely to occur if the absorber
producing the damping wing is adjacent to a region of the IGM that has a high
ionization rate, low density, and/or elevated temperature. Here we develop a
possible explanation of the transmission window motivated by these factors and
the nearness of the window to both an extended Ly$\alpha$+Ly$\beta$ absorption
trough and a bright quasar.

An association between neutral islands and low-density regions may occur
naturally near the end of reionization. Low-density regions are expected to be
among the last regions of the IGM to be ionized due to the scarcity of
ionizing sources. Indeed, long Ly$\alpha$ troughs such as the one towards
J0148 have been shown to trace low-density regions in terms of the number
density of nearby galaxies (Becker et al., 2018; Kashino et al., 2020;
Christenson et al., 2021). The edge of the J0148 proximity zone, which is
adjacent to a deep Ly$\alpha$+Ly$\beta$ trough, may be part of such a low-
density region.

The boundaries of neutral islands will also be the sites of ionization fronts,
and hence recently heated gas. In order for photoionization heating to
increase the contrast in transmission between the edge of the proximity zone
and the interior, the interior would need to have been reionized significantly
earlier, giving it time to cool. This might have occurred if the quasar turned
on after the interior was already ionized (although we note that if it was
ionized by soft spectra then the quasar might still have doubly ionize helium
in the interior, decreasing the temperature contrast with the edge). We can
estimate whether J0148 would be able to drive an ionization front at the edge
of its proximity zone using a simple Strömgren sphere model. Adopting the
nominal ionizing luminosity parameters from Section 3.1, assuming mean density
gas, and ignoring recombinations, producing an ionized region out to 34
$h^{-1}\,$Mpc would require a quasar lifetime of $\sim$30 Myr. Including the
radially varying ionizing opacity described in Section 3.1 would increase the
required time by roughly a factor of two, although somewhat less time would be
required if part of the proximity zone was already ionized when J0148 began
its luminous phase, or if the region near the edge of the proximity zone is
underdense. A lifetime of tens of Myr is consistent with, albeit at the upper
end of recent estimates of the lifetime of typical $z\sim 6$ quasars based on
the sizes of their proximity zones (Ďurovčíková et al., 2024). It is therefore
plausible that J0148 could be driving an ionization front at the end of its
proximity zone. There may also be additional sources contributing to the
ionizing flux. As noted above, Eilers et al. (2024) find a group of [O
iii]-emitting galaxies just redward of the J0148 transmission window, as well
as an overdensity of [O iii] emitters associated with the quasar itself.

The recent passage of an ionization front through a low-density region thus
provides a possible explanation for the J0148 transmission window. Driven by
ionizing photons from the quasar and/or nearby sources, the front would have
passed through the edge of a neutral island, leaving behind hot gas. Such gas
would be highly transparent to Ly$\alpha$ photons, producing a strong
transmission window that is then modified by the damping wing from the
foreground neutral island. This scenario is consistent with the unusual
strength and shape of the transmission window, its nearness to a giant
absorption trough, and the proximity of a bright quasar. As noted above, the
high transmission of the J0148 window may also be produced by an extremely
low-density region and/or enhanced ionization from nearby sources.

##  5 Summary

We have presented evidence that the longest and deepest known Ly$\alpha$
trough at $z<6$ is associated with damping wing absorption. We find that the
profile of a strong Ly$\alpha$ transmission window at the end of the J0148
proximity zone, adjacent to the red end of trough, requires a damping wing and
is highly unlikely to arise from resonant absorption alone. The lack of either
metal lines or directly detected galaxies further suggests that the damping
wing does not arise from a compact absorber such as a DLA. Instead, the we
find that the damping absorption is more likely to arise from an extended
neutral island, particularly given the close physical association between the
proximity zone transmission window and the extended Ly$\alpha$+Ly$\beta$
trough.

The characteristics of the J0148 transmission window suggest that it is
associated with a region of intrinsically low Ly$\alpha$ opacity. It may be a
deep void, a region where the local ionizing background is enhanced by local
sources in addition to the quasar, and/or a region with elevated gas
temperatures. One possibility is that the window traces hot gas that has been
recently ionized by the quasar and/or other nearby sources.

This work has focused on a somewhat unique line of sight, and our analysis has
been purposefully narrow in addressing the question of whether the
characteristics of an individual transmission feature require the presence of
a damping wing. We expect that a more comprehensive and robust analysis of
$z\lesssim 6$ quasar proximity zones, particularly those that end in extended
troughs, may provide further evidence of neutral islands. In a companion
paper, Zhu et al. (in prep) find a statistical association between
Ly$\alpha$+Ly$\beta$ troughs and damping wing absorption. These works are the
first direct evidence that extended Ly$\alpha$+Ly$\beta$ troughs at $z<6$
arise from islands of neutral gas, and hence that hydrogen reionization is
still ongoing at these redshifts.

## Acknowledgements

We thank the referee, Fred Davies, for helpful comments and suggestions. We
also thank Sarah Bosman for providing the J0148 [C ii] redshift ahead of
publication, as well as Anson D’Aloisio and Chris Cain for helpful
discussions. G.D.B., Y.Z., and S.H. were supported by the NSF through grant
AST-1751404. Y.Z. was also supported by the NSF through award SOSPADA-029 from
the NRAO. J.S.B. was supported by STFC consolidated grants ST/T000171/1 and
ST/X000982/1.

This work is based in part on observations made with ESO Telescopes at the La
Silla Paranal Observatory under program ID 084.A-0390.

The hydrodynamical simulations used in this work were performed with
supercomputer time awarded by the Partnership for Advanced Computing in Europe
(PRACE) 8th Call. We acknowledge PRACE for awarding us access to the Curie
supercomputer, based in France at the Tres Grand Centre de Calcul (TGCC). This
work has also used the DiRAC Data Analytic system at the University of
Cambridge, operated by the University of Cambridge High Performance Computing
Service on behalf of the STFC DiRAC HPC Facility (www.dirac.ac.uk). This
equipment was funded by BIS National E-infrastructure capital grant
(ST/K001590/1), STFC capital grants ST/H008861/1 and ST/H00887X/1, and STFC
DiRAC Operations grant ST/K00333X/1. DiRAC is part of the National
E-Infrastructure.

## Data Availability

The raw data used this article are available from the ESO archive at
http://archive.eso.org/cms.html. A reduced version of the J0148 X-Shooter
spectrum is available through a public github repository of the XQR-30 spectra
at https://github.com/XQR-30/Spectra.

## References

  * Asplund et al. (2009) Asplund M., Grevesse N., Sauval A. J., Scott P., 2009, ARAA, 47, 481
  * Bañados et al. (2016) Bañados E., et al., 2016, ApJS, 227, 11
  * Bañados et al. (2018) Bañados E., et al., 2018, Nature, 553, 473
  * Bañados et al. (2019) Bañados E., et al., 2019, ApJ, 885, 59
  * Barkana & Loeb (2001) Barkana R., Loeb A., 2001, Phys. Rep., 349, 125
  * Becker & Bolton (2013) Becker G. D., Bolton J. S., 2013, MNRAS, 436, 1023
  * Becker et al. (2001) Becker R. H., et al., 2001, AJ, 122, 2850
  * Becker et al. (2012) Becker G. D., Sargent W. L. W., Rauch M., Carswell R. F., 2012, ApJ, 744, 91
  * Becker et al. (2013) Becker G. D., Hewett P. C., Worseck G., Prochaska J. X., 2013, MNRAS, 430, 2067
  * Becker et al. (2015) Becker G. D., Bolton J. S., Lidz A., 2015, PASP, 32, e045
  * Becker et al. (2018) Becker G. D., Davies F. B., Furlanetto S. R., Malkan M. A., Boera E., Douglass C., 2018, ApJ, 863, 92
  * Becker et al. (2019) Becker G. D., et al., 2019, ApJ, 883, 163
  * Becker et al. (2021) Becker G. D., D’Aloisio A., Christenson H. M., Zhu Y., Worseck G., Bolton J. S., 2021, MNRAS, 508, 1853
  * Bolton et al. (2005) Bolton J. S., Haehnelt M. G., Viel M., Springel V., 2005, MNRAS, 357, 1178
  * Bolton et al. (2012) Bolton J. S., Becker G. D., Raskutti S., Wyithe J. S. B., Haehnelt M. G., Sargent W. L. W., 2012, MNRAS, 419, 2880
  * Bolton et al. (2017) Bolton J. S., Puchwein E., Haehnelt M. G., Kim T.-S., Meiksin A., Viel M., 2017, MNRAS, 464, 897
  * Bosman (2021) Bosman S. E. I., 2021, arXiv e-prints, arXiv:2108.12446
  * Bosman & Becker (2015) Bosman S. E. I., Becker G. D., 2015, MNRAS, 452, 1105
  * Bosman et al. (2018) Bosman S. E. I., Fan X., Jiang L., Reed S., Matsuoka Y., Becker G., Haehnelt M., 2018, MNRAS, 479, 1055
  * Bosman et al. (2022) Bosman S. E. I., et al., 2022, MNRAS, 514, 55
  * Cain et al. (2021) Cain C., D’Aloisio A., Gangolli N., Becker G. D., 2021, ApJ, 917, L37
  * Choudhury (2022) Choudhury T. R., 2022, GReGr, 54, 102
  * Choudhury et al. (2021) Choudhury T. R., Paranjape A., Bosman S. E. I., 2021, MNRAS, 501, 5782
  * Christenson et al. (2021) Christenson H. M., Becker G. D., Furlanetto S. R., Davies F. B., Malkan M. A., Zhu Y., Boera E., Trapp A., 2021, ApJ, 923, 87
  * Cooke et al. (2017) Cooke R. J., Pettini M., Steidel C. C., 2017, MNRAS, 467, 802
  * D’Aloisio et al. (2015) D’Aloisio A., McQuinn M., Trac H., 2015, ApJl, 813, L38
  * D’Aloisio et al. (2018) D’Aloisio A., McQuinn M., Davies F. B., Furlanetto S. R., 2018, MNRAS, 473, 560
  * D’Aloisio et al. (2019) D’Aloisio A., McQuinn M., Maupin O., Davies F. B., Trac H., Fuller S., Upton Sanderbeck P. R., 2019, ApJ, 874, 154
  * D’Aloisio et al. (2020) D’Aloisio A., McQuinn M., Trac H., Cain C., Mesinger A., 2020, ApJ, 898, 149
  * D’Odorico et al. (2018) D’Odorico V., et al., 2018, ApJ, 863, L29
  * D’Odorico et al. (2023) D’Odorico V., et al., 2023, MNRAS, 523, 1399
  * Davies & Furlanetto (2016) Davies F. B., Furlanetto S. R., 2016, MNRAS, 460, 1328
  * Davies et al. (2016) Davies F. B., Furlanetto S. R., McQuinn M., 2016, MNRAS, 457, 3006
  * Davies et al. (2018) Davies F. B., et al., 2018, ApJ, 864, 143
  * Davies et al. (2020) Davies F. B., Hennawi J. F., Eilers A.-C., 2020, MNRAS, 493, 1330
  * Davies et al. (2023a) Davies R. L., et al., 2023a, MNRAS, 521, 289
  * Davies et al. (2023b) Davies F. B., et al., 2023b, arXiv e-prints, arXiv:2312.08464
  * Dayal & Ferrara (2018) Dayal P., Ferrara A., 2018, Phys. Rep., 780, 1
  * Djorgovski et al. (2001) Djorgovski S. G., Castro S., Stern D., Mahabal A. A., 2001, ApJ, 560, L5
  * Ďurovčíková et al. (2024) Ďurovčíková D., et al., 2024, arXiv e-prints, arXiv:2401.10328
  * Eilers et al. (2018) Eilers A.-C., Davies F. B., Hennawi J. F., 2018, ApJ, 864, 53
  * Eilers et al. (2024) Eilers A.-C., et al., 2024, arXiv e-prints, arXiv:2403.07986
  * Fan et al. (2006) Fan X., et al., 2006, AJ, 131, 1203
  * Fan et al. (2023) Fan X., Bañados E., Simcoe R. A., 2023, ARAA, 61, 373
  * Furlanetto & Oh (2005) Furlanetto S. R., Oh S. P., 2005, MNRAS, 363, 1031
  * Gaikwad et al. (2020) Gaikwad P., et al., 2020, MNRAS, 494, 5091
  * Gaikwad et al. (2023) Gaikwad P., et al., 2023, MNRAS, 525, 4093
  * Garaldi et al. (2022) Garaldi E., Kannan R., Smith A., Springel V., Pakmor R., Vogelsberger M., Hernquist L., 2022, MNRAS
  * Greig et al. (2017) Greig B., Mesinger A., Haiman Z., Simcoe R. A., 2017, MNRAS, 466, 4239
  * Greig et al. (2019) Greig B., Mesinger A., Bañados E., 2019, MNRAS, 484, 5094
  * Greig et al. (2022) Greig B., Mesinger A., Davies F. B., Wang F., Yang J., Hennawi J. F., 2022, MNRAS, 512, 5390
  * Ishimoto et al. (2022) Ishimoto R., et al., 2022, MNRAS, 515, 5914
  * Kashino et al. (2020) Kashino D., Lilly S. J., Shibuya T., Ouchi M., Kashikawa N., 2020, The Astrophysical Journal, 888, 6
  * Keating et al. (2015) Keating L. C., Haehnelt M. G., Cantalupo S., Puchwein E., 2015, MNRAS, 454, 681
  * Keating et al. (2020a) Keating L. C., Weinberger L. H., Kulkarni G., Haehnelt M. G., Chardin J., Aubert D., 2020a, MNRAS, 491, 1736
  * Keating et al. (2020b) Keating L. C., Kulkarni G., Haehnelt M. G., Chardin J., Aubert D., 2020b, MNRAS, 497, 906
  * Kulkarni et al. (2019) Kulkarni G., Keating L. C., Haehnelt M. G., Bosman S. E. I., Puchwein E., Chardin J., Aubert D., 2019, MNRAS, 485, L24
  * Lewis et al. (2022) Lewis J. S. W., et al., 2022, MNRAS, 516, 3389
  * Lopez et al. (2016) Lopez S., et al., 2016, A&A, 594, A91
  * Lusso et al. (2015) Lusso E., Worseck G., Hennawi J. F., Prochaska J. X., Vignali C., Stern J., O’Meara J. M., 2015, MNRAS, 449, 4204
  * Lyke et al. (2020) Lyke B. W., et al., 2020, ApJS, 250, 8
  * Malloy & Lidz (2015) Malloy M., Lidz A., 2015, ApJ, 799, 179
  * McQuinn et al. (2011) McQuinn M., Hernquist L., Lidz A., Zaldarriaga M., 2011, MNRAS, 415, 977
  * Mesinger & Haiman (2004) Mesinger A., Haiman Z., 2004, ApJ, 611, L69
  * Mesinger & Haiman (2007) Mesinger A., Haiman Z., 2007, ApJ, 660, 923
  * Miralda-Escudé (1998) Miralda-Escudé J., 1998, ApJ, 501, 15
  * Mortlock et al. (2011) Mortlock D. J., et al., 2011, Nature, 474, 616
  * Nasir & D’Aloisio (2020) Nasir F., D’Aloisio A., 2020, MNRAS, 494, 3080
  * Planck Collaboration XVI (2014) Planck Collaboration XVI 2014, A&A, 571, A16
  * Qin et al. (2021) Qin Y., Mesinger A., Bosman S. E. I., Viel M., 2021, MNRAS, 506, 2390
  * Richards et al. (2011) Richards G. T., et al., 2011, AJ, 141, 167
  * Robertson (2022) Robertson B. E., 2022, ARAA, 60, 121
  * Satyavolu et al. (2023) Satyavolu S., Kulkarni G., Keating L. C., Haehnelt M. G., 2023, MNRAS, 521, 3108
  * Simcoe et al. (2012) Simcoe R. A., Sullivan P. W., Cooksey K. L., Kao M. M., Matejek M. S., Burgasser A. J., 2012, Nature, 492, 79
  * Songaila & Cowie (2002) Songaila A., Cowie L. L., 2002, AJ, 123, 2183
  * Wang et al. (2020) Wang F., et al., 2020, ApJ, 896, 23
  * Wise (2019) Wise J. H., 2019, ConPh, 60, 145
  * Yang et al. (2020a) Yang J., et al., 2020a, ApJL, 897, L14
  * Yang et al. (2020b) Yang J., et al., 2020b, ApJ, 904, 26
  * Zhu et al. (2021) Zhu Y., et al., 2021, ApJ, 923, 223
  * Zhu et al. (2022) Zhu Y., et al., 2022, ApJ, 932, 76
  * Zhu et al. (2023) Zhu Y., et al., 2023, ApJ, 955, 115

##  Appendix A quasar continuum estimation

Here we provide further details on the continuum estimation technique used to
normalize the J0148 spectrum. The continuum was constructed from a composite
of lower-redshift quasar spectra drawn from the Sloan Digital Sky Survey Data
Release 16 (DR16) quasar catalog (Lyke et al., 2020). We use spectra with
$S/N>5$ per pixel that are not flagged as broad absorption line objects
(BALs). We also use only wavelength redward of 3800 Å, avoiding bluer
wavelengths where the noise in DR16 spectra tends to increase. As we describe
below, the spectra are chosen based on their similarity to J0148 over
wavelengths redward of the Ly$\alpha$ forest. They are then corrected for
intervening Ly$\alpha$ and Ly$\beta$ absorption and averaged to create an
estimate of the unabsorbed continuum.

In many respects our approach is similar to other matched composite techniques
in the literature, especially those where C iv emission is used to select the
matching spectra (e.g., Mortlock et al., 2011; Simcoe et al., 2012; Bosman &
Becker, 2015). A key difference with these and other approaches, however, is
that we statistically correct for foreground Ly$\alpha$ and Ly$\beta$
absorption. Provided that there are enough objects contributing to the
composite to largely overcome the scatter in IGM transmission between lines of
sight, this approach, which we describe bellow, allows us to avoid the
potentially complicating step of fitting continua over the forest of the
lower-redshift spectra (for a review of continuum reconstruction techniques,
see Bosman et al., 2022).

Our first step is to fit a spline to the J0148 spectrum redward of the
Ly$\alpha$ forest. This provides a basis for comparison that minimizes the
impact of metal absorption lines and noise. We note that J0148 exhibits
extended C iv broad absorption (Figure 1) that obscures the continuum between
the Si iv and C iv emission lines. This part of the spectrum is not used to
select DR16 spectra, however, and thus the presence of the BAL does not impact
our continuum estimate. We also note that the red edge of the BAL falls
$\sim$20,000 km s-1 blueward of the systemic redshift. Any corresponding broad
absorption in N v or Ly$\alpha$ would therefore fall blueward of the proximity
zone and not impact our $\langle\Delta F\rangle$ measurement.

The DR16 spectra are treated flexibly in two ways. First, we allow the
redshifts to vary. This accounts for redshift errors in the DR16 catalog (as
well as for J0148, in principle, although it has an accurate systemic redshift
from [C ii] 158 $\mu$m emission; see Bosman et al., in prep). It also provides
a means of accommodating velocity shifts of the emission lines with respect to
the systemic redshift. This is particularly helpful for J0148, which has a C
iv emission line blueshift of $\sim$3000 km s-1. Allowing the redshifts to
vary provides a pool of potential matches that is significantly larger than
what would be available if the DR16 spectra were fixed at their (potentially
inaccurate) catalog redshifts. The results of this approach are tested below.
For each object, we start by adopting a redshift that places the peak of the C
iv emission line at the same rest-frame wavelength as the J0148 C iv peak. We
then allow the DR16 spectrum to shift in velocity by up to $\pm$1000 km s-1 in
increments of 250 km s-1. Second, we adjust the overall slope and scaling of
the the DR16 spectra to match J0148. We do this by multiplying the DR16
spectrum by a power law in wavelength. At each velocity offset, the power law
is fit to the median fluxes over rest-frame wavelengths 1270–1390 Å and
1600–1750 Å.

Our primary goal is to estimate the continuum over the J0148 transmission
window near $\lambda\simeq 8370$–8400 Å, which falls at the blue edge of the
Ly$\alpha$ emission line. The shifted and scaled DR16 spectra are thus
compared to J0148 over the two wavelength regions redward of the forest that
best capture the variation in the strength and shape of Ly$\alpha$: the red
side of the Ly$\alpha$+N v complex (1216–12151 Å rest-frame) and the C iv
emission line (1500–1600 Å rest-frame). These regions are highlighted in
Figure 1. The inclusion of C iv is motivated by the known correlation between
C iv and Ly$\alpha$ emission (e.g., Richards et al., 2011; Greig et al., 2017;
Davies et al., 2018). We note that by including the red side of Ly$\alpha$+N v
we assume that there is no significant damping wing absorption redward of 1216
Å. For J0148 this is probably a reasonable assumption given the extent of the
proximity zone and the lack of any obvious DLAs within it (Figure 2). We
compute the median values of the J0148 spline and the individual DR16 spectra
in bins of 7 Å (rest-frame) over these regions, where the binning is chosen to
minimize the impact of noise and absorption lines in the DR16 spectrum while
preserving the overall spectral shape. The goodness-of-fit is determined from
the mean squared difference between the binned values. We create composites of
the 200 best-fitting spectra, where the number is chosen to balance the
accuracy of fits to the emission lines with the suppression of scatter over
the forest.

Before computing the composite, each individual spectrum is corrected for
foreground Ly$\alpha$ and Ly$\beta$ absorption. This is done by diving each
pixel in the Ly$\alpha$ and Ly$\beta$ forests by the expected mean
transmission at the corresponding redshift, taking into account the quasar
proximity effect. The baseline mean Ly$\alpha$ transmission is taken from
Becker et al. (2013). Baseline Ly$\beta$ transmission values are scaled from
the Ly$\alpha$ values using a hydrodynamical simulation of the IGM.
Specifically, we use outputs from the 40-2048 run of the Sherwood simulation
suite (Bolton et al., 2017). We rescale the ionization rates to match the mean
Ly$\alpha$ transmission at a given redshift, and then measure the mean
Ly$\beta$ transmission from the resulting Ly$\beta$ optical depths. The
proximity effect is taken into account using a model similar to the one
described in Becker et al. (2021). For each object, we compute the distance,
$R_{\rm eq}$, at which the ionization rate from the quasar would be equal to
the background ionization rate in the absence of attenuation or redshifting of
the quasar photons. Following Becker et al. (2021), we compute the ionizing
luminosity of the quasar based on its absolute magnitude at rest-frame 1450 Å,
assuming a broken power-law of the form $L\propto\nu^{-\alpha}$ with
$\alpha=0.6$ over $912~{}{\mbox{\AA}}<\lambda_{\rm rest}<1450~{}{\mbox{\AA}}$,
and $\alpha=1.5$ at $\lambda_{\rm rest}<912$ Å (see Lusso et al., 2015).
Background ionization rates are interpolated from values in Becker & Bolton
(2013) over $2.4<z<4.0$ and Bolton et al. (2005) at $z=2.0$. Using this value
of $R_{\rm eq}$, we calculate the expected ratio of the total (quasar $+$
background) ionization rate to the background rate as a function of distance
from the quasar. The expected mean Ly$\alpha$ and Ly$\beta$ transmission at
each distance is then computed after rescaling the ionization rate in the
simulation by this ratio.

We note that this method does not take into account the fact that quasars will
tend to live in overdensities. It will therefore tend to under-correct the
absorption close to the quasar redshift, and hence leave a flux deficit in the
composite near the peak of the Ly$\alpha$ emission line. We find, however,
that this effect is significantly reduced when using DR16 spectra at lower
redshifts, where associated absorbers stand out more strongly and can be
rejected more easily (see below). We therefore generate composites from two
samples. The first composite, $F_{{\rm hi}z}$, is drawn from $\sim$16,000
spectra over $2.91<z<4.0$. The upper redshift bound is chosen to limit the
correction for foreground absorption, while the lower bound is chosen to
provide coverage down to Ly$\gamma$ at observed wavelengths redward of 3800 Å.
The second composite, $F_{{\rm low}z}$, is drawn from $\sim$32,000 spectra
over $2.2<z<2.5$, where the redshift range is intended to reduce the impact of
associated absorption. The two composites for J0148 are generally similar to
within $\sim$5% over the Ly$\alpha$ forest, except near the peak of
Ly$\alpha$, where $F_{{\rm hi}z}$ has a $\sim$12% deficit with respect to
$F_{{\rm low}z}$. The final composite is an average of the two at wavelengths
where they each have coverage from at least 100 individual spectra, apart from
1195 Å to 1216 Å in the rest-frame (the Ly$\alpha$ peak) where only $F_{{\rm
low}z}$ is used. The composite shown in Figure 1 uses $F_{{\rm hi}z}$ at
$\lambda<7909$ Å, $F_{{\rm low}z}$ over
$8353~{}{\mbox{\AA}}<\lambda<8499~{}{\mbox{\AA}}$ (the Ly$\alpha$ peak), and
an average of the two at other wavelengths.

We apply a modest rejection scheme when calculating the mean composite values.
At a given rest-frame wavelength, individual pixels are rejected if they are
more than 3$\sigma$ outliers from the median. They are also rejected if they
are more than 3$\sigma$ outliers after smoothing the individual spectra by a
Gaussian kernel with ${\rm FWHM}=1000$ km s-1. The first criterion is mainly
meant to reject sharp features such as bad pixels or strong metal lines, while
the second criterion is meant to exclude extended features such as DLAs. We
note, however, that in the forest the scatter includes the sightline-to-
sightline variation in the Ly$\alpha$ and Ly$\beta$ absorption, making
unwanted features more difficult to detect and the rejection relatively mild.

Finally, the uncertainty in the composite continuum is estimated directly from
the scatter among the individual spectra. We compute the ratio of the
individual fluxes (corrected for Ly$\alpha$ and Ly$\beta$ transmission) to the
composite flux within a sliding window of 5000 km s-1. The 1$\sigma$
uncertainty in the composite is taken to be the standard deviation in this
ratio. We note that this approach is conservative in that it includes both the
scatter in the intrinsic quasar continua and the scatter in the Ly$\alpha$ and
Ly$\beta$ transmission between sightlines, which will be non-negligible even
on 5000 km s-1 scales.

Figure 15: Example continuum fits to three quasars at $z\sim 4$ drawn from the
XQ-100 sample. The spectra are from the VIS arm of X-Shooter. Similar to
Figure 1, the two wavelength regions to which the DR16 spectra are matched are
shown with yellow lines, which are largely obscured by the composites. The raw
composites are plotted as blue histograms. Spline fits to the raw composite
are shown with a smooth red line, while the 1$\sigma$ uncertainty in the
composite is given by the shaded region. Note that the composites generally
trace the tops of the transmission peaks over the Ly$\alpha$ forest.

We test our continuum estimation technique by applying it to lower-redshift
objects. At $z\lesssim 4$, the intrinsic continuum over the Ly$\alpha$ and
Ly$\beta$ forests can be roughly inferred from the peaks of the transmission
features, providing a check on the composite results. Motivated by our focus
on J0148, we choose three $z\sim 4$ objects for illustration with strong C iv
blueshifts that have high-quality X-Shooter spectra from the XQ-100 sample
(Lopez et al., 2016) : J0244$-$0134 ($z_{\rm eq}=4.064$), J0525$-$3343
($z_{\rm em}=4.435$), and J1034$+$1102 ($z_{\rm eq}=4.297$). The fits to these
objects are shown in Figure 15. We note that we are plotting spectra from the
VIS arm of X-Shooter only in order to avoid uncertainties in the scaling
between the UVB and VIS arms. Nevertheless, the continuum estimates roughly
trace the tops of the Ly$\alpha$ forest transmission peaks. For this work we
have not attempted to evaluate the accuracy of our technique for a broad range
of quasar spectral shapes; however, the results for these objects provide some
confidence that composite shown in Figure 1 provides a reasonable estimate of
the unabsorbed J0148 continuum.

##  Appendix B Numerical convergence

Our simulated lines of sight are drawn from the 80-2048 run from the Sherwood
simulation suite (Bolton et al., 2017). This uses an 80 $h^{-1}\,$Mpc box with
$2\times 2048^{3}$ particles, giving a gas particle mass of $\sim$$8\times
10^{5}~{}h^{-1}~{}{\rm M}_{\odot}$. To check whether our results depend on
either box size or resolution, we repeat our analysis of simulated proximity
zones ending in LLSs (Section 3.3) using the 80-1024, 40-2048 and 40-1024 runs
from the same suite. The 40-1024 run has the same mass resolution as the
80-2048 run but is a factor of eight smaller in volume. The 80-2048 and
80-1024 runs, meanwhile, have the same volume but differ in mass resolution by
a factor of eight. The 40-2048 and 40-1024 differ in mass resolution by the
same factor.

Figure 16: Results of numerical convergence tests. The four lines show the
$\langle\Delta F\rangle$ distributions for simulated proximity zones ending in
LLSs for a luminosity scaling of $L=L_{0}$. The solid line is for a 80-2048
run and is the same as in Figure 6. The dotted, dashed, and dash-dotted lines
are for the 80-1024, 40-2048 and 40-1024 runs, respectively.

The results for our nominal quasar luminosity scaling ($L=L_{0}$) are shown in
Figure 16. The 80-2048 and 40-1024 runs produce nearly identical results,
suggesting that the 80-2048 run is not missing a substantial number of rare
voids that would increase the number of lines of sight with low $\langle\Delta
F\rangle$ values. The higher resolution of the 80-2048 run compared to the
80-1024 run pushes the $\langle F\rangle$ distribution to slightly higher
values, largely because taller, sharper transmission peaks prevent the
$\langle F\rangle$ region of the damping wing template from overlapping as
much of the proximity zone transmission (see Section 3.3). A similar though
somewhat smaller shift is seen between the 40-1024 and 40-2048 runs. The
results suggests that further increasing the resolution would tend to disfavor
the LLS model even more strongly, although the effect would be modest. We
conclude that numerical effects are unlikely to have a significant impact on
our results.

##  Appendix C High column-density compact absorbers

Figure 17: Similar to Figure 6, but for simulated proximity zones ending in
compact absorbers with a range of column densities. Results shown here are for
luminosity scaling $L=2L_{0}$. As a reference, the solid line is for LLSs with
column density $\log{(N_{\rm H\,I}/{\rm cm^{-2}})}=18.0$, and is the same as
the one plotted in Figure 6. Other lines are for H i column densities of
$\log{(N_{\rm H\,I}/{\rm cm^{-2}})}=19.0$ (dot), 19.5 (dash-dot-dot), 20
(dash-dot), and 20.5 (dash) Figure 18: Probability that $\langle\Delta
F\rangle$ is equal to or less than the J0148 value for simulated proximity
zones ending in high-column density compact absorbers, as a function of H i
column density. The lower, middle, and upper solid lines are for $L_{\rm
q}=L_{0}$, $2L_{0}$, and $4L_{0}$, respectively. The probability tends to peak
between $\log{(N_{\rm H\,I}/{\rm cm^{-2}})}=20.0$ and 20.5. A horizontal line
at $P=0.01$ is plotted for reference.

In Section 3.4 we focused on compact absorbers with the canonical DLA H i
column density of $\log{(N_{\rm H\,I}/{\rm cm^{-2}})}=20.3$. Here we provide
results for a wider range of column densities. Figure 17 shows $\langle\Delta
F\rangle$ distributions for compact absorbers with column densities ranging
from $\log{(N_{\rm H\,I}/{\rm cm^{-2}})}=19.0$ to 20.5, along with the Lyman
limit case (i.e., no damping wing) from Section 3.3. For clarity, we plot the
$L_{\rm q}=2L_{0}$ scaling only. The $\langle\Delta F\rangle$ distribution
shifts towards smaller values with increasing $N_{\rm H\,I}$ up to
$\log{(N_{\rm H\,I}/{\rm cm^{-2}})}=20.0$, but then shifts back towards higher
values at $\log{(N_{\rm H\,I}/{\rm cm^{-2}})}=20.5$. This behavior reflects
the change in the compact absorber damping wing profile with column density.
Figure 18 plots the probability that a line of sight will have $\langle\Delta
F\rangle$ equal to or less than the J0148 value as a function of $N_{\rm
H\,I}$. Here we plot the results for $L=L_{0}$, $2L_{0}$, and $4L_{0}$;
$0.5L_{0}$ is not shown because it does not produce a significant number of
lines of sight with small enough $\langle\Delta F\rangle$ values. Even $L_{\rm
q}$ requires $\log{(N_{\rm H\,I}/{\rm cm^{-2}})}\geq 19.5$ in order for the
probability to exceed 1%. In each case the probability peaks near
$\log{(N_{\rm H\,I}/{\rm cm^{-2}})}\simeq 20.2$. Our focus on the canonical
DLA column density of $\log{(N_{\rm H\,I}/{\rm cm^{-2}})}=20.3$ is thus a
conservative choice in that it roughly maximizes the probability that a
compact absorber will produce a transmission feature that is similar to J0148
in terms of $\langle\Delta F\rangle$.