other constraints on reionization such as the electron optical depth to the cosmic microwave background and the apparent attenuation of Ly α emission from quasars and galaxies at z 7 (Kulkarni et al. 2019 ;Keating et al. 2020a , b ;Nasir & D'Aloisio 2020 ;Choudhury, Paranjape & Bosman 2021 ;Qin et al. 2021 ). A strong indication that reionization was still underway at z ∼ 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 α opacity fluctuations (Gaikwad et al. 2023 ;Davies et al. 2024 ). The increase in the mean free path from z ∼ 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 photoe v aporation 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 α absorption saturates for even very small hydrogen neutral fractions ( x H I 10 -4 ). The best potential tracers of neutral regions are the 'dark gaps' in the L y α and L y β forests (e.g. Djorgovski et al. 2001 ;Becker et al. 2001 ;Songaila & Cowie 2002 ;Fan et al. 2006 ;Zhu et al. 2021Zhu et al. , 2022 ) ). These are extended ( 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 ; MNRAS 533, 1525-1540 (2024) 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 ). Observ ationally, ho we ver, it is dif ficult 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 ).

In this work and a companion paper (Zhu et al. 2024 ), we present the first direct evidence that at least some of the extended Ly α + Ly β 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 α 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 α system (DLA).

The dark gap we examine is part of the giant Ly α trough towards ULAS J0148 + 0600 (herein J0148). This trough is remarkable partly because it is the longest ( ℓ 1110 h -1 Mpc) and most opaque ( τ eff 7) 1 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 ;Ba ˜ nados et al. 2018 ;Davies et al. 2018 ;Greig, Mesinger & Ba ˜ nados 2019 ;Wang et al. 2020 ;Yang et al. 2020a ;Greig et al. 2022 ), and recent efforts that have extended this search down to z = 6 ( Ďuro v č íko v á et al. 2024 ), but there are tw o k ey 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 α transmission in the proximity zone, rather than from the otherwise unabsorbed continuum redwards of the Ly α 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 α forest absorption. This is the approach taken by Zhu et al. ( 2024 , see also Spina et al. 2024 ), while here we focus on the damping wing from a single Ly α trough. Our work is perhaps most closely related to Mesinger & Haiman ( 2004, 2007 ), who argued that the distribution of Ly α fluxes near the proximity zone boundaries in three z > 6 . 2 quasars fa v ours cosmological H II regions within a significantly neutral IGM. Here, we go a step further, showing that the Ly α 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 1 τ eff = -log F , where F is the continuum-normalized transmission.

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 cold dark matter cosmology with ( m , , 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. J 0 1 4 8 + 0 6 0 0 LY α T R A N S M I S S I O N W I N D OW

The primary feature examined in this paper is an extended window of Ly α transmission centred at λ 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, F , that will be used to test for the presence of damping wing absorption.

Our analysis is based on a 10-h 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 ∼130-220 per 10 km s -1 pix el o v er the transmission window and the adjacent Ly α trough. The spectrum is shown in Fig. 1 . We adopt a redshift for J0148 of z q = 5 . 9896, which is based on [C II ] 158 µm emission-line measurements (Bosman pri v ate communication).

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 Sk y Surv e y Data Release 16 (DR16; Lyke et al. 2020 ) that are good matches to J0148 in terms of their transmission in two wavelength intervals redwards of the Ly α forest. The intervals co v er the red side of the Ly α + N V emission line and the C IV emission line, which are highlighted yellow in Fig. 1 . Ly α emission is of particular importance for this work because the transmission window of interest lies on the blue side of the Ly α 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 α (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 α and Ly β 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 Fig. 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 σ continuum uncertainty of 10 per cent. The uncertainty estimate is conserv ati ve in that part of this scatter will be due to variations in the amount of Ly α forest absorption between lines of sight rather than intrinsic variations in the DR16 spectra.

The continuum-normalized J0148 spectrum co v ering the Ly α proximity zone and the red end of the giant Ly α trough is displayed in Fig. 2 . We also show the corresponding section of the Ly β forest. The most conspicuous feature of the Ly α proximity zone is a strong (maximum transmission 0.7) and extended transmission window spanning λ 8370-8400 Å ( z 5 . 89-5.91). Ignoring peculiar ve-Figure 1. X-Shooter spectrum of ULAS J0148 + 0600. The black line shows the observed flux, in arbitrary F λ units, binned to 50 km s -1 pixels for display. The short yellow lines are spline fits to the flux o v er the red side of the Ly α + 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 σ uncertainty. See the text for details. The dotted line shows our adopted power-law continuum for the non-ionizing UV emission, L ν ∝ ν -0 . 6 ( F λ ∝ λ -1 . 4 ). The depression in the flux with respect to the composite near 10 000 Å is C IV broad absorption. The 110 h -1 comoving Mpc Ly α trough appears as the region of zero flux over 7930-8362 Å. The transmission feature examined in this work is centred at ∼8385 Å.

Figure 2. J0148 proximity zone. The top and bottom panels show the regions of the Ly β and Ly α forests, respectively, corresponding to the proximity zone, where the wa velengths ha ve been scaled to align the two panels vertically in redshift (top axis of Ly β panel). The continuum-normalized flux is plotted as a histogram, with the 1 σ 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 α and Ly β transmission. The cyan shaded region in the lower panel shows the difference between the damping profile and the flux o v er a 1000 km s -1 region at the red edge of the damping wing. locities, the wavelength range corresponds to a line-of-sight extent of ∼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 α trough, which is also dark in Ly β o v er z = 5 . 849-5.878. Secondly, the transmission on the blue side increases relatively smoothly from zero up to its maximum.

Moti v ated by the shape of this transmission window and the potential association of dark Ly α + Ly β troughs with neutral islands, MNRAS 533, 1525-1540 (2024) 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 α + Ly β dark gap bluewards of the proximity zone. The Ly α and Ly β transmission at the red and blue ends of this gap, respecti vely, allo w 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 H I / cm -2 ) 20 . 5. We note, ho we ver, 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 α and Ly β transmission profiles for the potential island are shown in Fig. 2 .

The damping wing of the neutral island roughly traces the blue side of the broad Ly α transmission window. We also note that the entire Ly α (and Ly β) 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 α 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 e v aluated based on their similarity to a damping wing profile at the blue end of the proximity zone. Specifically, we will measure the mean difference, F , between the proximity zone transmission and the damping wing template for a neutral island similar to the one shown in Fig. 2 . In order to focus on the part of the spectrum most impacted by the damping wing, and to simplify the modelling of the proximity zones, we measure F o v er 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 o v er this re gion is shaded in Fig. 2 . Without an y adjustment of the profile parameters we find F = 0 . 15 for J0148. Below, we test whether this value indicates that a genuine damping wing is likely to be present.

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 5000 line of sight from the 80-2048 run, which has a box size of 80 h -1 Mpc and 2048 3 gas particles with a mass of 7 . 97 × 10 5 h -1 M . The lines of sight are not drawn through the massive halos expected to host bright quasars; ho we ver, this is unlikely to affect our analysis, which is focused on features that are offset by more than ∼30 h -1 Mpc from the quasar position. By comparison, the virial radius for a 10 12 h -1 M halo will be ∼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 ( ρ) = T 0 ( ρ/ ρ ) γ -1 , adopting T 0 = 12000 K and γ = 1 . 2, broadly consistent with measurements from Gaikwad et al. ( 2020 ) o v er 5 . 3 < z < 5 . 9. Baseline Ly α optical depths are set by globally rescaling the hydrogen neutral fractions to produce a mean Ly α 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 photoionized 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 H I ∝ ( / bg ) -1 , where the local photoionization rate, = bg + 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 κ 912 = κ bg 912 ( / bg ) ξ , and therefore increases with distance from the quasar following equation ( 8) in Becker et al. ( 2021 ). We adopt a background opacity value of log ( κ bg / ( 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 ξ 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. 2018D'Aloisio et al. , 2020 ; ;Becker et al. 2021 ). For the background photoionization rate we adopt bg = 0 . 17 × 10 -12 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 ˜ nados et al. 2016 ). The UV spectrum is modelled as a double power law of the form L ν ∝ ν -α , with α = 0 . 6 between 1450 Å and 912 Å, and 1.5 shortwards of 912 Å (see Lusso et al. 2015 ;Becker et al. 2021 , and references therein). This agrees well with our composite-based continuum o v er the non-ionizing UV part of the spectrum (Fig. 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 v ariability, v ariations in ionizing opacity along the line of sight not captured by the model abo v e, fluctuations in the ionizing background due to additional sources, and gas temperature. The neutral fraction at the edge of the proximity zone will also depend on the distance from the quasar. In practice, ho we ver, 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 q , but note that these variations may also result from other factors such as those listed abo v e, some of which could act locally where the broad transmission window lies. We examine these factors further in Section 4.2 . The scaling of x H I for optically thin regions within the proximity zones is shown in Fig. 3 . We refer to the scaling given by the nominal parameters abo v e as L q = L 0 , and compute proximity zones for L q = 0 . 5 L 0 , L 0 , 2.0 L 0 , and 4.0 L 0 . At the distance to the edge of the J0148 proximity zone, ∼33 h -1 Mpc, the neutral fraction scaling ranges from 0.60 (i.e. 60 per cent of the value in the absence of the proximity effect) for 0 . 5 L 0 to 0.05 for 4 L 0 . As sho wn belo w, this range is suf ficient to e v aluate 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 q on the ionizing opacity, following the discussion abo v e. The total optical depth within the proximity zone is τ 1 . 66, 0.98, 0.57, and 0.34 for L q = 0 . 5 L 0 , L 0 , 2 L 0 , and 4 L 0 , respectiv ely. F or reference, we plot the zero-opacity case for L 0 as a dotted line in Fig. 3 . showing the corresponding observed wavelength for a quasar at z 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 the text for details). The vertical dashed line corresponds to the edge of transmitted Ly α 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 .

Finally, we apply a random continuum error to each simulated line of sight. The continuum error is modelled 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 o v er the transmission window at ∼8385 Å (Fig. 2 ). We note that this does not include a possible wavelength dependence in the continuum; ho we ver, the analysis belo w 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.

We e v aluate the proximity models belo w based on ho w frequently they produce F 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. F is calculated as the difference between the proximity zone transmission and this template o v er a 1000 km s -1 region starting where the red wing of the damping profile first exceeds 1 per cent. Conserv ati vely, we minimize F 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 F 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. Transmission profiles for mean-density neutral islands of length 2.5 (dot), 5.0 (dash-double-dotted), 7.5 (solid), 10 (dash-dotted), and 15 (dash) h -1 Mpc. The near edge for all islands is 34 h -1 Mpc from a quasar at z q = 5 . 9896. The shaded region shows the 1000 km s -1 interval used to calculate F . We use the 7.5 h -1 Mpc template for our primary analysis, which is moti v ated by the L y α and L y β transmission near the end of the J0148 proximity zone. Our results do not substantially change, ho we ver, when other templates are used.

Secondly, there may be variations in the x H I scaling that are not well modelled by the radial profiles plotted in Fig. 3 . Thirdly, 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 a v oiding these more nuanced aspects of the modelling. In the future, ho we ver, a more robust modelling 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 o v er the entire proximity zone (Fig. 2 ).

It is worth considering whether the details of our F measurement may impact the results below. For example, we have chosen to measure F o v er a 1000 km s -1 re gion partly because this is the extent to which the J0148 transmission window appears to follow the Ly α damping wing in Fig. 2 . Our main result, that we exclude models without damping wings (see Section 3.3 ), actually becomes somewhat stronger if we measure F 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 , F 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 o v er 8401-8417 Å. It is expected that resonant Ly α forest absorption will become increasingly dominant o v er damping wing absorption at large velocities, making F a less sensitive statistic o v er wider windows. F or 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 F . As described abo v e, our choice of a template generated from a 7.5 h -1 Mpc neutral region at mean density is moti v ated by the spacing between L y α and L y β transmission peaks in the J0148 spectrum. For comparison, in Fig. 4 we show damping profiles for mean-density neutral islands with lengths ranging from 2.5 to 15 h -1 Mpc, which are aligned at their edges nearest to the quasar. We note MNRAS 533, 1525-1540 (2024)

Figure 5. Examples of simulated proximity zones ending at Lyman limit systems with column density log ( N H I / cm -2 ) = 18 . 0. The transmitted flux in each panel has been rescaled by a random amount to emulate continuum errors (see the 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 o v er a 1000 km s -1 re gion. F or each line of sight, the damping profile is shifted in wavelength such that it minimizes F 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, ho we ver (see the text for details). There are three lines of sight plotted for each quasar luminosity scaling in Fig. 3 , where the scaling is indicated in the upper left corner of each panel. The lines of sight correspond to the 5 th , 50 th , and 95 th percentiles (top to bottom) in the F distribution for each L q . The value of F is given in parentheses.

that contiguous islands longer than ∼7.5 h -1 Mpc are not allowed towards J0148 by the Ly β transmission peaks in Fig. 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 F 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 o v er the template island length by taking the minimum F from among the profiles shown in Fig. 4 and again found no substantial change. Finally, we note that for the shape of the red wing, variations in the density of the neutral island are largely degenerate with variations in length.

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 H I / 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 q = 0 beyond the LLS. We also produce a complete absorption trough outside of the proximity zone by arbitrarily increasing x 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 Fig. 5 . We plot three lines of sight for each L q scaling, which are chosen to lie MNRAS 533, 1525-1540 (2024) Figure 6. Cumulative probability distributions of F for simulated lines of sight ending in Lyman limit systems. The solid lines show P ( < F ) for the four quasar luminosity scalings, with colours correspond to Fig. 3 . For reference, the left-to-right ordering of these lines at P ( < F ) = 0 . 5 is L q = 2 L 0 , L 0 , 4 L 0 , 0 . 5 L 0 . The vertical dashed line shows F for J0148.

at the 5 th , 50 th , and 95 th percentiles in F for that scaling (bottom to top for each L q ). Similar to Fig. 2 , the template damping profile is plotted in grey and the shaded region shows the 1000 km s -1 region o v er which F is measured. The full F distributions are plotted in Fig. 6 . The smallest F values tend to be produced by L q = L 0 and 2 L 0 . By comparison, 0 . 5 L 0 lines of sight tend to exhibit only weak transmission peaks at the edge of the proximity zone, while the strong transmission for 4 L 0 tends to prohibit the F region of the damping profile from extending over the transmission peaks. In all cases, ho we ver, the LLS models fail to generate F values equal to or less than the J0148 value in more than 0.01 per cent of cases. We therefore conclude that lines of sight without significant damping wings are highly unlikely to produce the kind of smooth Ly α transmission feature seen at the blue end of the J0148 proximity zone.

Next, we e v aluate whether the J0148 feature could be created via damping wing absorption from a compact source such as a damped Ly α 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 H I / 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 F value is then found for each combination of N H I and L q . In this section, we focus on log ( N H I / cm -2 ) = 20 . 3 because it is the canonical DLA value and close the column density where the probability of producing F 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 H I / cm -2 ) = 20 . 3 are shown in Fig. 7 . Similar to Fig. 5 , lines of sight at the 5 th , 50 th , and 95 th percentiles in F are shown for each L q scaling. The full F distributions are plotted in Fig. 8 . We find that a significant percentage of lines of sight with F values equal or less than the J0148 value are produced for the higher L q values (3 per cent for 2 L q , 18 per cent for 4 L q ). As seen in Fig. 7 , in cases where F 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 F 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 o v er the full proximity zone than is seen in J0148 (e.g. comparing the bottom panel of Fig. 2 to the bottom two panels of Fig. 7 o v er 8400 Å λ 8500 Å). As noted abo v e, however, in this work we conserv ati vely consider only the 1000 km s -1 region within which F 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 α feature in the J0148 spectrum.

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 H I = 1 o v er a 7.5 h -1 Mpc interval starting 34 h -1 Mpc from the quasar. The choices of length and distance are moti v ated by the largest mean-density neutral island allowed by the observed L y α and L y β transmission at the end of the J0148 proximity zone (see Section 2 ). We note that in terms of F , the distance will be largely degenerate with the L q scaling. Example simulated proximity zones ending in neutral islands are shown in Fig. 9 . Similar to Fig. 5 , lines of sight at the 5 th , 50 th , and 95 th percentiles in F are shown for each L q scaling. The full F distributions are plotted in Fig. 10 . The results are broadly similar to those for DLAs in that the higher L q values can produce significant percentages of lines of sight with F values at or below the J0148 value (1 per cent for 2 L 0 , 10 per cent for 4 L 0 ). In these cases (e.g. bottom two panels of Fig. 9 ), the template damping profile used to measure F (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 F 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 F 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 F distributions are shown in Fig. 10 , and we note that we are still calculating F 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 F value similar to the J0148 value, the agreement is not strongly dependent on the island length.

The abo v e analysis demonstrates that a smooth transmission feature at the end of proximity zone similar to the J0148 Ly α transmission MNRAS 533, 1525-1540 (2024) 6 , but for simulated proximity zones ending in compact DLAs with column density log ( N H I / cm -foot_1 ) = 20 . 3. Colours correspond to the quasar luminosity scalings in Fig. 3 . The left-to-right ordering of the solid lines is L q = 4 L 0 , 2 L 0 , L 0 , 0 . 5 L 0 . The vertical dashed line shows F for J0148.

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, ho we ver, that while both scenarios are allowed in terms of their F statistics, a neutral island is fa v oured based on other factors.

There are at least three reasons to fa v our 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. ( 2023 ), 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 o v er the redshifts of the Ly α + Ly β trough at the end of the proximity zone ( z 5 . 85-5.88; Fig. 11 ). The deep X-Shooter data excludes systems with log ( N O I / cm -2 ) 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 ∼ 10 19 . 5 cm -2 in order for the probability of generating a F value equal to or less than the J0148 values to be at least 1 per cent, even for the highest L 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 [O / H] -3 . 0, 2 which is comparable to the most metal-poor DLA known (Cooke, Pettini & Steidel 2017 ). We show the expected line profiles for such an absorber in Fig. 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. 2024 ). The non-detection of Si II λ1260 would increase the metallicity constraint to [O / H] -3 . 1 for the same abundance ratios. The presence of such an ultra-lo w-metallicity, lo w-ionization absorber is clearly ruled 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 F o v er the 1000 km s -1 shaded region at the red edge of the template damping wing.

Figure Similar to Fig. 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. Colours correspond to the quasar luminosity scalings in Fig. 3 . The left-to-right ordering of the line groups is L q = 4 L 0 , 2 L 0 , L 0 , 0 . 5 L 0 . The vertical dashed line shows F for J0148.

out at these redshifts. We also show the expected C IV absorption for log ( N C IV / 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, ho we ver, that some DLAs at these redshifts may have metallicities near [O / H] -3 (e.g. D 'Odorico et al. 2018 ;Ba ˜ nados et al. 2019 ), so the lack of metal lines on its own does not fully rule out a compact absorber.

A second reason to fa v our a neutral island is that the evidence of a damping wing is found adjacent to an extremely long and opaque giant Ly α trough. In our simulated proximity zones, we set the Ly α transmission to zero bluewards of the proximity zone by artificially increasing the H I neutral fraction. In more physically moti v ated models, ho we ver, it has been sho wn 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. ( 2024 ) show that long dark troughs o v er 5 z 6 are statistically associated with damping wing absorption, which is consistent with a model in which the troughs arise from neutral islands (see also Spina et al. 2024 ).

Finally, a deep search for [O III ]-emitting galaxies with the JWST as part of the EIGER surv e y rev eals no sources within the potential DLA redshift window of 5 . 85 z 5 . 88 (Eilers et al. 2024 ). Based on the lack of metal lines, the presence of an extreme Ly α trough, and the lack of any directly detected galaxy, we therefore conclude that the damping wing o v er the J0148 Ly α transmission window is more likely to arise from an extended neutral island than a compact absorber.

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

It is clear from Figs 8 and 10 that the probability of producing a line of sight with F equal to or less than the J0148 value increases 9 , but for simulated lines of sight ending in 7.5 h -1 Mpc neutral islands that have F values similar to the J0148 value. Three lines of sight are shown for quasar luminosity scalings L q = L 0 (top), 2 L 0 (middle), and 4 L 0 (bottom). Lines of sight similar to J0148 are extremely rare in these simple models for L q = L 0 , while they are more common for 2 L 0 and 4 L 0 (Fig. 10 ).

towards higher quasar luminosities. The higher ionization rates produce lower neutral fractions (Fig. 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 q (Fig. 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 F f 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 rewards of the F 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 ∼ 6 quasar proximity zones F values within ±0.05 of the J0148 value. Colours correspond to the quasar luminosity scaling values in Fig. 3 , although only L q = L 0 (green), 2 L 0 (orange), and 4 L 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.

only weakly depend on the host halo mass. The fact that our lines of sight do not start at massive haloes is therefore unlikely to account for this difference in transmission (see also Section 3.1 ). We also note that our composite-based quasar continuum (Fig. 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, Hennawi & Eilers 2020 ;Satya v olu et al. 2023 ).

This inconsistency redwards of the F region, particularly with the L q = 2 L 0 and 4 L 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 5 . 92, coincident with the strong Ly α absorption just redwards of the J0148 window. In order for the probability of producing a F value consistent with J0148 to exceed 1 per cent, similar to the 2 L 0 case in Fig. 10 , such sources would need to roughly triple the ionization rate of our baseline L q = L 0 case,foot_2 which is already seven times the background rate at the centre of the transmission window at z = 5 . 897 ( λ obs = 8385 Å; Fig. 3 ). Other ways of boosting the J0148 window transmission include low densities and/or high gas temperatures, which we examine below.

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 F values within ±0.05 of the J0148 value (i.e. 0.10 to 0.20) for L q = L 0 , 2 L 0 , 4 L 0 . We then compute the median density along these lines of sight separately for each L q . The results are shown in Fig. 13 . The median profiles show clear underdensities from ∼24 to ∼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 q , suggesting that the role of density is especially important when the local ionization rate is lower. The underdensities required for the L q = L 0 case are extremely rare, ho we ver, which Figure 14. Cumulative probability distributions of F for simulated lines of sight ending in neutral islands. The solid lines are the same as those plotted in Fig. 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 q = 4 L 0 (red), 2 L 0 (orange), L 0 (green), and 0 . 5 L 0 (blue). In terms of the change in F , increasing the gas temperature by this amount is roughly equi v alent to doubling the quasar luminosity.

is why only ∼0.1 per cent of these lines of sight have F values equal to or less than the J0148 value of 0.15, and only ∼1 per cent have F f 0 . 2.

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, ho we ver, if the gas was recently reionized. We tested the impact of ele v ated temperatures on the F 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 postionization 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 ion = ṅ ion /n H I , where ṅ ion is the incident flux of ionizing photons and n H I is the neutral hydrogen number density, a front driven by a quasar with L q = L 0 would be moving at ∼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 ∼ 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 ∼30 000 K at this velocity for our adopted ionizing spectral index of α = 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 ∼40 000 K if helium is doubly ionized by the quasar (Bolton et al. 2012 ;Davies, Furlanetto & McQuinn 2016 ). For simplicity, we increase the to 30 000 K along the entire line of sight, but again note that we are e v aluating F only within a narrow region at the end of the proximity zone.

The resulting F distributions are plotted in Fig. 14 . Increasing the temperature to 30 000 K substantially increases the number of lines of sight with F values equal to or less than J0148, roughly MNRAS 533, 1525-1540 (2024) comparable to the effect of doubling the quasar luminosity. We note that most of the change in the F distribution comes from the lower neutral fractions associated with suppressed recombination rates rather than from increased thermal broadening. In terms of testing this scenario, ho we ver, the width of the resonant absorption lines in this region, measured with high-resolution spectroscopy, may provide some indication of the gas temperature (e.g. Bolton et al. 2012 ).

In the preceding sections, we have shown that a smooth Ly α 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 ele v ated 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 α + Ly β 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 α 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 α + Ly β 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 ∼30 Myr. Including the radially varying ionizing opacity described in Section 3.1 would increase the required time by roughly a factor of 2, 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 ∼ 6 quasars based on the sizes of their proximity zones ( Ďuro v č íko v á 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 abo v e, Eilers et al. ( 2024 ) find a group of [O III ]-emitting galaxies just redwards of the J0148 transmission window, as well as an o v erdensity 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 windo w. Dri ven 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 α 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 abo v e, the high transmission of the J0148 window may also be produced by an extremely low-density region and/or enhanced ionization from nearby sources.

We have presented evidence that the longest and deepest known Ly α trough at z < 6 is associated with damping wing absorption. We find that the profile of a strong Ly α 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 α + Ly β trough.

The characteristics of the J0148 transmission window suggest that it is associated with a region of intrinsically low Ly α 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 ele v ated gas temperatures. One possibility is that the windo w 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 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. ( 2024 ) find a statistical association between L y α + L y β troughs and damping wing absorption (see also Spina et al. 2024 ). These works are the first direct evidence that extended Ly α + Ly β troughs at z < 6 arise from islands of neutral gas, and hence that hydrogen reionization is still ongoing at these redshifts.

MNRAS 533, 1525-1540 (2024) 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.

MNRAS 533, 1525-1540 (2024) (DR16) quasar catalogue (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 redwards of 3800 Å, a v oiding 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 o v er wav elengths redwards of the Ly α forest. They are then corrected for intervening Ly α and Ly β 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, ho we ver, is that we statistically correct for foreground L y α and L y β absorption. Provided that there are enough objects contributing to the composite to largely o v ercome the scatter in IGM transmission between lines of sight, this approach, which we describe bello w, allo ws us to a v oid the potentially complicating step of fitting continua o v er the forest of the lower redshift spectra (for a re vie w of continuum reconstruction techniques, see Bosman et al. 2022 ).

Our first step is to fit a spline to the J0148 spectrum redwards of the Ly α forest. This provides a basis for comparison that minimizes the impact of metal absorption lines and noise. We note that J0148 e xhibits e xtended C IV broad absorption (Fig. 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, ho we ver, and thus the presence of the BAL does not impact our continuum estimate. We also note that the red edge of the BAL falls ∼20 000 km s -1 bluewards of the systemic redshift. Any corresponding broad absorption in N V or Ly α would therefore fall bluewards of the proximity zone and not impact our F measurement.

The DR16 spectra are treated flexibly in tw o w ays. First, we allow the redshifts to vary. This accounts for redshift errors in the DR16 catalogue (as well as for J0148, in principle, although it has an accurate systemic redshift from [C II ] 158 µm emission). 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 ∼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) catalogue 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 ±1000 km s -1 in increments of 250 km s -1 . Secondly, we adjust the o v erall slope and scaling of 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 o v er the J0148 transmission window near λ 8370-8400 Å, which falls at the blue edge of the Ly α emission line. The shifted and scaled DR16 spectra are thus compared to J0148 o v er the two wavelength regions redwards of the forest that best capture the variation in the strength and shape of Ly α: the red side of the Ly α + N V complex (1216-12151 Å rest-frame) and the C IV emission line (1500-1600 Å rest-frame). These regions are highlighted in Fig. 1 . The inclusion of C IV is moti v ated by the known correlation between C IV and Ly α 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 α + N V we assume that there is no significant damping wing absorption redwards of 1216 Å. For J0148 this is probably a reasonable assumption given the extent of the MNRAS 533, 1525-1540 (2024) proximity zone and the lack of any obvious DLAs within it (Fig. 2 ). We compute the median values of the J0148 spline and the individual DR16 spectra in bins of 7 Å (rest-frame) o v er these regions, where the binning is chosen to minimize the impact of noise and absorption lines in the DR16 spectrum while preserving the o v erall 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 o v er the forest.

Before computing the composite, each individual spectrum is corrected for foreground L y α and L y β absorption. This is done by diving each pixel in the Ly α and Ly β forests by the expected mean transmission at the corresponding redshift, taking into account the quasar proximity effect. The baseline mean Ly α transmission is taken from Becker et al. ( 2013 ). Baseline Ly β transmission values are scaled from the Ly α 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 α transmission at a given redshift, and then measure the mean Ly β transmission from the resulting Ly β 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 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 ∝ ν -α with α = 0 . 6 o v er 912 Å < λ rest < 1450 Å, and α = 1 . 5 at λ rest < 912 Å (see Lusso et al. 2015 ). Background ionization rates are interpolated from values in Becker & Bolton ( 2013 ) o v er 2 . 4 < z < 4 . 0 and Bolton et al. ( 2005 ) at z = 2 . 0. Using this value of R 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 α and Ly β 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 o v erdensities. It will therefore tend to undercorrect the absorption close to the quasar redshift, and hence leave a flux deficit in the composite near the peak of the Ly α emission line. We find, ho we ver, 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 hi z , is drawn from ∼16 000 spectra o v er 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 pro vide co v erage down to Ly γ at observed wavelengths redwards of 3800 Å. The second composite, F low z , is drawn from ∼32 000 spectra o v er 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 ∼5 per cent o v er the Ly α forest, e xcept near the peak of Ly α, where F hi z has a ∼12 per cent deficit with respect to F 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 α peak) where only F low z is used. The composite shown in Fig. 1 uses F hi z at λ < 7909 Å, F low z o v er 8353 Å < λ < 8499 Å (the Ly α 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 σ outliers from the median. They are also rejected if they are more than 3 σ outliers after smoothing the individual spectra by a Gaussian kernel with 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, ho we ver, that in the forest the scatter includes the sightline-to-sightline variation in the Ly α and Ly β 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 α and Ly β transmission) to the composite flux within a sliding window of 5000 km s -1 . The 1 σ uncertainty in the composite is taken to be the standard deviation in this ratio. We note that this approach is conserv ati ve in that it includes both the scatter in the intrinsic quasar continua and the scatter in the Ly α and Ly β transmission between sightlines, which will be non-ne gligible ev en on 5000 km s -1 scales.

We test our continuum estimation technique by applying it to lower-redshift objects. At z 4, the intrinsic continuum o v er the L y α and L y β forests can be roughly inferred from the peaks of the transmission features, providing a check on the composite results. Moti v ated by our focus on J0148, we choose three z ∼ 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 eq = 4 . 064), J0525 -3343 ( z em = 4 . 435), and J1034 + 1102 ( z eq = 4 . 297). The fits to these objects are shown in Fig. A1 . We note that we are plotting spectra from the VIS arm of X-Shooter only in order to a v oid uncertainties in the scaling between the UVB and VIS arms. Nevertheless, the continuum estimates roughly trace the tops of the Ly α forest transmission peaks. For this work, we have not attempted to e v aluate the accuracy of our technique for a broad range of quasar spectral shapes; ho we ver, the results for these objects provide some confidence that composite shown in Fig. 1 provides a reasonable estimate of the unabsorbed J0148 continuum.

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 × 2048 3 particles, giving a gas particle mass of ∼ 8 × 10 5 h -1 M . 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. The results for our nominal quasar luminosity scaling ( L = L 0 ) are shown in Fig. B1 . 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 F values. The higher resolution of the 80-2048 run compared to the 80-1024 run pushes the F distribution to slightly higher values, largely because taller, sharper transmission peaks prevent the F region of the damping wing template from o v erlapping 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 result