# Damping Wing-Like Features in the Stacked Ly$\alpha$ Forest: Potential
Neutral Hydrogen Islands at $z<6$

Yongda Zhu,1,2 George D. Becker,2 Sarah E. I. Bosman,3,4 Christopher Cain,5,2
Laura C. Keating,6 Fahad Nasir,3 Valentina D’Odorico,7,8,9 Eduardo Bañados,3
Fuyan Bian,10 Manuela Bischetti,11 James S. Bolton,12 Huanqing Chen,13 Anson
D’Aloisio,2 Frederick B. Davies,3 Rebecca L. Davies,14,15 Anna-Christina
Eilers,16 Xiaohui Fan,1 Prakash Gaikwad,3 Bradley Greig,17,15 Martin G.
Haehnelt,18 Girish Kulkarni,19 Samuel Lai,20 Ewald Puchwein,21 Yuxiang
Qin,17,15 Emma V. Ryan-Weber,14,15 Sindhu Satyavolu,19 Benedetta Spina,4
Fabian Walter,3 Feige Wang,1 Molly Wolfson,22 and Jinyi Yang1
1Steward Observatory, University of Arizona, 933 North Cherry Avenue, Tucson,
AZ 85721, USA
2Department of Physics & Astronomy, University of California, Riverside, CA
92521, USA
3Max-Planck-Institut für Astronomie, Königstuhl 17, D-69117 Heidelberg,
Germany
4Institute for Theoretical Physics, Heidelberg University, Philosophenweg 12,
D-69120, Heidelberg, Germany
5School of Earth and Space Exploration, Arizona State University, Tempe, AZ
85281, USA
6Institute for Astronomy, University of Edinburgh, Blackford Hill, Edinburgh,
EH9 3HJ, UK
7INAF-Osservatorio Astronomico di Trieste, Via Tiepolo 11, I-34143 Trieste,
Italy
8Scuola Normale Superiore, Piazza dei Cavalieri 7, I-56126 Pisa, Italy
9IFPU-Institute for Fundamental Physics of the Universe, via Beirut 2, I-34151
Trieste, Italy
10European Southern Observatory, Alonso de Córdova 3107, Casilla 19001,
Vitacura, Santiago 19, Chile
11Dipartimento di Fisica, Sezione di Astronomia, Università di Trieste, via
Tiepolo 11, 34143 Trieste, Italy
12School of Physics and Astronomy, University of Nottingham, University Park,
Nottingham, NG7 2RD, UK
13Canadian Institute for Theoretical Astrophysics, University of Toronto,
Toronto, ON M5R 2M8, Canada
14Centre for Astrophysics and Supercomputing, Swinburne University of
Technology, Hawthorn, Victoria 3122, Australia
15ARC Centre of Excellence for All Sky Astrophysics in 3 Dimensions (ASTRO
3D), Australia
16MIT Kavli Institute for Astrophysics and Space Research, 77 Massachusetts
Avenue, Cambridge, MA 02139, USA
17School of Physics, University of Melbourne, Parkville, VIC 3010, Australia
18Kavli Institute for Cosmology and Institute of Astronomy, Madingley Road,
Cambridge, CB3 0HA, UK
19Tata Institute of Fundamental Research, Homi Bhabha Road, Mumbai 400005,
India
20Commonwealth Scientific and Industrial Research Organisation (CSIRO), Space
& Astronomy, P. O. Box 1130, Bentley, WA 6102, Australia
21Leibniz-Institut für Astrophysik Potsdam, An der Sternwarte 16, 14482
Potsdam, Germany
22Department of Physics, University of California, Santa Barbara, CA 93106,
USA E-mail<EMAIL_ADDRESS>0000-0003-3307-7525 0000-0003-2344-263X
0000-0001-8582-7012 0000-0001-9420-7384 0000-0001-5211-1958
0000-0003-0294-8674 0000-0003-3693-3091 0000-0002-2931-7824
0000-0002-1620-0897 0000-0002-4314-021X 0000-0003-2764-8248
0000-0002-3211-9642 0000-0003-2344-263X 0000-0003-0821-3644
0000-0002-3324-4824 0000-0003-2895-6218 0000-0003-3310-0131
0000-0002-2423-7905 0000-0002-4085-2094 0000-0001-8443-2393
0000-0001-5829-4716 0000-0001-9372-4611 0000-0001-8778-7587
0000-0002-4314-1810 0000-0002-5360-8103 0000-0001-5818-6838
0000-0003-1634-1283 0000-0003-4793-7880 0000-0002-7633-431X
0000-0003-0766-2499 0000-0001-5287-4242

(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

Recent quasar absorption line observations suggest that reionization may end
as late as $z\approx 5.3$. As a means to search for large neutral hydrogen
islands at $z<6$, we revisit long dark gaps in the Ly$\beta$ forest in
VLT/X-Shooter and Keck/ESI quasar spectra. We stack the Ly$\alpha$ forest
corresponding to the edges of these Ly$\beta$ dark gaps and identify a damping
wing-like extended absorption profile. The average redshift of the stacked
forest is $z=5.8$. By comparing these observations with reionization
simulations, we infer that such a damping wing-like feature can be naturally
explained if these gaps are at least partially created by neutral islands.
Conversely, simulated dark gaps lacking neutral hydrogen struggle to replicate
the observed damping wing features. Furthermore, this damping wing-like
profile implies that the volume-averaged neutral hydrogen fraction must be
$\langle x_{\rm HI}\rangle\geq 6.1\pm 3.9\%$ at $z=5.8$. Our results offer
robust evidence that reionization extends below $z=6$.

###### keywords:

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

††pubyear: 2024††pagerange: Damping Wing-Like Features in the Stacked
Ly$\alpha$ Forest: Potential Neutral Hydrogen Islands at $z<6$–Damping Wing-
Like Features in the Stacked Ly$\alpha$ Forest: Potential Neutral Hydrogen
Islands at $z<6$

##  1 Introduction

Hydrogen reionization carries key implications for the formation and evolution
of the first stars, galaxies, and supermassive black holes (see, e.g., Fan et
al., 2023 for a review). Cosmic microwave background (CMB) observations
indicate a midpoint of reionization at $z\sim 7-8$ (Planck Collaboration et
al., 2020, see also de Belsunce et al., 2021). Meanwhile, transmitted flux
observed in the Ly$\alpha$ forest at $z\sim 6$ towards high-redshift quasars
(e.g., Fan et al., 2006) has long been interpreted as an indicator of the end
of reionization. However, if reionization ended by $z\sim 6$, galaxies would
have had to produce an extremely large amount of ionizing photons to complete
reionization within a short timeframe (see also Muñoz et al., 2024). This may
pose a significant challenge to our understanding of star formation, the
escape fraction of ionizing photons, and the spectral energy distribution
(SED) of early galaxies (e.g., Robertson et al., 2015; Bouwens et al., 2015;
Finkelstein et al., 2019; Stark, 2016).

A later end to reionization, which would help resolve such tension, is gaining
support from recent observations. Supporting evidence includes large-scale
fluctuations in the Ly$\alpha$ effective optical depth111Defined as $\tau_{\rm
eff}=-\ln{\langle F\rangle}$, where $F$ is the continuum-normalized
transmission flux. measured in quasar spectra (e.g., Fan et al., 2006; Becker
et al., 2015; Eilers et al., 2018; Bosman et al., 2018, 2022; Yang et al.,
2020b); long troughs extending down to or below $z\simeq 5.5$ in the
Ly$\alpha$ and Ly$\beta$ forests (e.g., Becker et al., 2015; Zhu et al., 2021,
2022), potentially indicating the existence of large neutral intergalactic
medium (IGM) islands (e.g., Kulkarni et al., 2019; Keating et al., 2020; Nasir
& D’Aloisio, 2020; Qin et al., 2021); observed underdensities around long dark
gaps traced by Ly$\alpha$ emitting galaxies (LAEs, Becker et al., 2018;
Kashino et al., 2020; Christenson et al., 2021, 2023); the evolution of metal-
enriched absorbers at $z\sim 6$ (e.g., Becker et al., 2019; Cooper et al.,
2019; Davies et al., 2023a, b; Sebastian et al., 2024); and the dramatic
evolution in the measured mean free path of ionizing photons over $5<z<6$
(Becker et al., 2021; Zhu et al., 2023, see also Bosman, 2021; Gaikwad et al.,
2023; Satyavolu et al., 2023; Roth et al., 2024; Davies et al., 2024), which
provides the most unambiguous evidence to date of ongoing reionization at
$z<6$. A late-end reionization scenario is also consistent with numerical
models that reproduce a variety of observations (e.g., Weinberger et al.,
2019; Choudhury et al., 2021; Qin et al., 2021; Gaikwad et al., 2023; Asthana
et al., 2024, Y. Qin et al. in prep.).

A pressing question is whether we can detect neutral hydrogen islands at $z<6$
if they exist. In terms of resonant Ly$\alpha$ and Ly$\beta$ absorption, the
observational signatures of neutral islands may be difficult to distinguish
from those of large-scale fluctuations in the ionizing ultraviolet background;
either one may produce extended regions of nearly zero transmission. At higher
redshifts, where the IGM is more neutral, damping wing features over the
Ly$\alpha$ emission of quasars and galaxies have been observed and used as a
powerful probe of the volume-averaged neutral hydrogen fraction ($x_{\rm HI}$)
of the Universe (e.g., Bañados et al., 2018; Davies et al., 2018; Wang et al.,
2020; Yang et al., 2020a; Greig et al., 2022, 2024; Umeda et al., 2023;
Ďurovčíková et al., 2024). Any signature of damping wing absorption at lower
redshifts might similarly indicate the presence of neutral islands. A
detection of damping wings would also enable constraints of $x_{\rm HI}$ that
would complement the upper limits from the fraction of dark gaps and dark
pixels in the forest (McGreer et al., 2013; Zhu et al., 2022; Jin et al.,
2023).

Malloy & Lidz (2015) proposed to test for neutral islands at $z\sim 5.5$ by
searching for damping wing absorption over the stacked Ly$\alpha$ forest
adjacent to highly absorbed regions. This method has not been implemented in
observations because it requires a large sample of high-quality spectra of
high-$z$ quasars and a relatively complete catalog of intervening metal
absorbers to exclude contamination from damped Ly$\alpha$ systems (DLAs)
arising from galaxies. Thanks to a large sample of quasar spectra provided in
the (Extended) XQR-30 large program (D’Odorico et al., 2023) and the metal
absorber catalog of Davies et al. (2023a), we are now able to carry out the
experiment proposed in Malloy & Lidz (2015). Long dark gaps have been detected
in the Ly$\beta$ forest (Zhu et al., 2022), and they indicate regions of high
IGM opacity that may potentially host neutral islands (e.g., Kulkarni et al.,
2019; Keating et al., 2020; Nasir & D’Aloisio, 2020; Qin et al., 2021).
Therefore, stacking the Ly$\alpha$ forest at the redshifts of these long dark
gaps is a powerful way to implement the Malloy & Lidz (2015) test and search
for IGM damping wing features at $z<6$.

This letter is organized as follows. In Section 2 we describe the data and
Ly$\beta$ dark gaps used in this work. Section 3 presents the stacked
Ly$\alpha$ profile. Sections 4 compares our results to model predictions, and
Section 5 presents our constraints on $x_{\rm HI}$. Finally, we conclude our
findings in Section 6. Throughout this paper, we quote distances in comoving
units unless otherwise noted and assume a $\Lambda$CDM cosmology with
$\Omega_{\rm m}=0.308$, $\Omega_{\Lambda}=0.692$, and $h=0.678$ (Planck
Collaboration et al., 2014). In addition to the stacking method used in this
work, our companion paper, Becker et al. (2024) present the discovery of an
IGM damping wing at $z<6$ towards the individual quasar sightline, J0148+0600.

##  2 Data

Table 1: Ly$\beta$ Dark Gaps Used in this Work

No. | Quasar | $z_{\rm q}^{\rm ref}$ | $z_{\rm blue}$ | $z_{\rm red}$ | $L_{\rm gap}$
---|---|---|---|---|---
1 | ULASJ1319+0950 | $6.1330^{\rm vi}$ | 5.876 | 5.903 | 8
2 | PSOJ060+24 | $6.1793^{\rm ii}$ | $5.833^{\rm a}$ | 5.856 | 7
3 | PSOJ108+08 | $5.9647^{\rm ii}$ | 5.674 | 5.751 | 24
4 | SDSSJ0842+1218 | $6.0763^{\rm iii}$ | 5.784 | 5.830 | 14
5 | SDSSJ2315-0023 | $6.124^{\rm i}$ | 5.790 | 5.883 | 28
6 | SDSSJ2315-0023 | $6.124^{\rm i}$ | 5.897 | $5.937^{\rm b}$ | $\geq 12$
7 | CFHQSJ1509-1749 | $6.1225^{\rm iii}$ | 5.800 | 5.863 | 19
8 | CFHQSJ1509-1749 | $6.1225^{\rm iii}$ | 5.870 | 5.910 | 12
9 | SDSSJ2054-0005 | $6.0391^{\rm vi}$ | 5.751 | 5.790 | 12
10 | SDSSJ0840+5624 | $5.8441^{\rm vi}$ | 5.585 | 5.607 | 7
11 | PSOJ340-18 | $5.999^{\rm i}$ | 5.774 | 5.810 | 11
12 | PSOJ065-26 | $6.1877^{\rm iii}$ | 5.954 | 5.988 | 10
13 | PSOJ007+04 | $6.0008^{\rm iii}$ | 5.741 | 5.780 | 12
14 | ULASJ0148+0600 | $5.9896^{\rm ii}$ | 5.654 | 5.735 | 25
15 | ULASJ0148+0600 | $5.9896^{\rm ii}$ | 5.741 | 5.803 | 19
16 | PSOJ217-16 | $6.1498^{\rm iii}$ | 5.807 | 5.873 | 20
17 | PSOJ217-16 | $6.1498^{\rm iii}$ | 5.920 | $5.961^{\rm b}$ | $\geq 12$
18 | J0408-5632 | $6.0264^{\rm ii}$ | 5.715 | 5.741 | 8
19 | PSOJ359-06 | $6.1718^{\rm iv}$ | 5.886 | 5.917 | 9
20 | PSOJ025-11 | $5.8414^{\rm ii}$ | 5.526 | 5.613 | 28
21 | PSOJ025-11 | $5.8414^{\rm ii}$ | 5.632 | 5.661 | 9
22 | PSOJ158-14 | $6.0681^{\rm iv}$ | 5.764 | 5.836 | 22
23 | PSOJ158-14 | $6.0681^{\rm iv}$ | 5.843 | 5.880 | 11
24 | SDSSJ1250+3130 | $6.137^{\rm v}$ | 5.836 | 5.870 | 10

Notes. Columns: (1) Index of dark gaps; (2) quasars used in Zhu et al. (2021);
(3) quasar redshift with reference; (4) redshift at the blue end of the gap;
(5) redshift at the red end of the gap; (6) dark gap length in units of
$h^{-1}\,\mathrm{Mpc}$. a Dark gap starting at the blue edge of the Ly$\beta$
forest. b Dark gap ending at the red edge of the Ly$\beta$ forest as defined
in Zhu et al. (2022). Quasar redshift references: i. Becker et al. (2019), ii.
S. Bosman et al. (in prep), iii. Decarli et al. (2018), iv. Eilers et al.
(2020), v. Shen et al. (2019), vi. Wang et al. (2010, 2013).

To create a stacked spectrum of the Ly$\alpha$ forest over potentially neutral
regions at $z<6$, we revisit dark gaps detected in the Ly$\beta$ forest in Zhu
et al. (2022). We refer the readers to Zhu et al. (2021, 2022) for details of
the data and dark gap detection. Briefly, Ly$\beta$ dark gaps are detected in
42 spectra of quasars at $5.77\lesssim z_{\rm quasar}\lesssim 6.31$. The
spectra are taken with the Echellette Spectrograph and Imager (ESI) on Keck
(Sheinis et al., 2002) and the X-Shooter spectrograph on the Very Large
Telescope (VLT; Vernet et al., 2011). Among these, 19 X-Shooter spectra are
from the XQR-30 large program (D’Odorico et al., 2023). A dark gap in the
Ly$\beta$ forest is defined as a continuous spectral region in which all
pixels binned to $1{h^{-1}{\rm Mpc}}$ have an observed normalized flux
$F=F_{\rm obs}/F_{\rm c}<0.02$, where $F_{\rm obs}$ is the observed flux and
$F_{\rm c}$ is the continuum flux predicted from Principal Component Analysis
(PCA; see Bosman et al., 2021, 2022; Zhu et al., 2021). The Ly$\alpha$ forest
at the corresponding redshift of a Ly$\beta$ dark gap also needs to be opaque
($F_{\rm Ly\alpha}=F_{\rm obs}/F_{\rm c}<0.05$).

In this work we use only the Zhu et al. (2022) dark gaps at $z<6$. Dark gap
statistics based on mock spectra generated from simulations suggest that
longer dark gaps are more likely to have a higher covering fraction of neutral
hydrogen (Nasir & D’Aloisio, 2020; Zhu et al., 2022). The stack based on
longer dark gaps might also create a stronger damping wing profile in the
corresponding Ly$\alpha$ forest (Malloy & Lidz, 2015). To balance the sample
size and dark gap length, here we use Ly$\beta$ dark gaps with $L\geq
7{h^{-1}{\rm Mpc}}$. Considering the low transmission in the Ly$\alpha$ forest
at these redshifts, we only include quasar spectra with a continuum-normalized
flux error $<0.05$ per pixel in the Ly$\alpha$ and Ly$\beta$ forest.
Therefore, dark gaps toward quasars PSO J308-21 and VIK J2318-3029 are
excluded due to the relatively low S/N of the spectra.

We have carefully checked these dark gaps to ensure they do not contain
intervening metal absorbers within 3000 km s-1 near the edge of dark gaps,
based on the metal line catalog presented in Davies et al. (2023a) and through
visual inspection. This requirement rejects the dark gap toward quasar SDSS
J0818+1722 spanning $z=5.761-5.794$ with $L=10{h^{-1}{\rm Mpc}}$. This dark
gap contains absorption lines of O i, C ii, Si ii, C iv, etc., near the red
edge at $z\simeq 5.79$. We note that, however, including this dark gap
produces no substantial difference in the results. A summary of Ly$\beta$ dark
gaps used in this work is provided in Table 1. Our final sample includes 24
dark gaps, for which the average redshift at the red and blue edges is
$z=5.8$.

##  3 Stacked Ly$\alpha$ Forest

Following Malloy & Lidz (2015), we create a stacked Ly$\alpha$ transmission
profile by aligning the spectra based on the redshift at the edges of
Ly$\beta$ dark gaps.222When aligning the spectra at the redshift corresponding
to the blue edge of Ly$\beta$ dark gaps, we flip the sign of the velocity
offset ($\Delta v$) such that the velocity offset increases when moving away
from the dark gap edges. We firstly bin the continuum-normalized Ly$\alpha$
forest to 80 km s-1 intervals, corresponding to $\sim 0.5{h^{-1}{\rm Mpc}}$ at
$z=5.8$, which is half of the bin size used by Zhu et al. (2022) when
searching for dark gaps. To remove the redshift evolution in the mean
Ly$\alpha$ transmission, we re-normalize the sightlines with a factor equal to
the mean transmission at $z=5.8$ divided by the mean transmission at the
redshift of each pixel according to the measurements in Bosman et al. (2022)
with a linear interpolation. The normalization factor is small and close to 1,
ensuring minimal distortion of the data. The stacked spectrum is then created
by taking a mean value of the flux in each bin. We exclude the quasar
proximity zone effect by conservatively disregarding the Ly$\alpha$ forest
within 11 proper-Mpc towards the quasar following Zhu et al. (2022). We also
mask pixels that are affected by sky-line subtraction residuals or telluric
absorption correction as indicated by peaks in the flux error array.

Figure 1 shows the stacking result along with the Ly$\alpha$ forest
corresponding to each individual Ly$\beta$ dark gap. We see an extended
damping wing-like absorption feature redward of the edge of the dark gaps,
with the transmission gradually recovering to the mean value at a velocity
offset of $\Delta v\sim 1000$ km s-1. Notably, such a damping wing-like
feature is not obvious in any of the individual spectra.333Although strong
transmission spikes appear adjacent to the blue edge of the Ly$\beta$ dark gap
toward ULAS J1319+0950, we have tested that excluding this sightline will not
change the damping wing-like profile significantly.

Following Malloy & Lidz (2015, see also , ), the absorption due to the
presence of a neutral hydrogen island will be extended, with optical depth far
from the line center described by a damping wing. The optical depth at a
velocity offset $\Delta v$ from the edge of a neutral island will be

$\tau_{\rm Ly\alpha}^{\rm DW}(\Delta v)\approx\frac{\tau_{\rm
GP}R_{\alpha}c}{\pi}\left[\frac{1}{\Delta v}-\frac{1}{\Delta v+v_{\rm
ext}}\right]\,,$ (1)

where $\tau_{\rm GP}$ is the Gunn-Peterson optical depth, for which we take a
nominal value of $2.5\times 10^{5}$ at $z\sim 5.8$,
$R_{\alpha}\equiv\Gamma_{\alpha}\lambda_{\alpha}/4\pi c$, and
$\Gamma_{\alpha}=6.265\times 10^{8}$ s-1 is the Ly$\alpha$ decay constant. We
adopt $v_{\rm ext}$ as an empirical parameter to describe the extent of the
transmission profile.

We fit the damping wing profile in Equation 1 to the upper envelope of the
stacked spectrum in the interval between $\Delta v=0$ km s-1 and $1000$ km
s-1. Specifically, we perform a least-squares fitting over the stacked flux,
and require that the fitted curve cannot be significantly lower than the
stacked flux in any pixel (e.g., Prochaska et al., 2003). Here, we set the
threshold to be the smaller of $80\%$ of the stacked flux or $1-\sigma$ white
noise below the flux. The white noise is given by the standard deviation
measured over the transmission-free region in the stacked Ly$\alpha$ forest
over $\Delta v=-1000$ km s-1 to $0$ km s-1. We have tested that using a
different threshold does not change our conclusions. As Figure 1 shows, the
fitted curve well describes the observed extended absorption feature and
yields $v_{\rm ext}=357$ km s-1. We note that, as Figure 1 shows, the stacked
Ly$\alpha$ transmission profile is comparable with that expected for neutral
islands with $L=5.34{h^{-1}{\rm Mpc}}$ in Malloy & Lidz (2015). Similar
damping wing-like signals at $z\lesssim 6$ are also presented in B. Spina et
al. (in prep), although a different method is used in their paper.

We emphasize that our approach does not represent a physical fit to a damping
wing per se, but rather a fit to a mean transmission profile using a damping
wing-like function. The observed mean profile potentially includes damping
wing absorption arising from genuine neutral islands, but will also include
resonant Ly$\alpha$ absorption, which may be modified by on the scales of
interest by UVB fluctuations. As described in Section 4, we assess the
evidence for damping wing absorption by fitting the same function to mock
stacks drawn from simulations.

We note that there are large fluctuations in the stacked transmission profile
at $\Delta v\gtrsim 1000~{}{\rm km\,s^{-1}}$. These are likely statistical
fluctuation resulting from the relatively small number of spectra included in
the stack. Alternatively, these absorption dips may come from individual
strong absorbers in the ionized regions redward of the gaps in the stack. A
clustering of neutral islands on the scales of $\sim 10{h^{-1}{\rm Mpc}}$,
however, may also create absorption features redward the damping wing-like
profile in the stack, although this hypothesis would require more data to
test.

Figure 1: Stacked Ly$\alpha$ forest corresponding to the redshift of edges of
Ly$\beta$ dark gaps (displayed in red). The grey curve represents the fitted
damping wing profile, as described by Equation 1, applied to the upper
envelope of the stacked flux between $\Delta v=0$ km s-1 and $1000$ km s-1.
The individual spectra used in this stack are shown in blue (sorted by their
$\Delta v$ coverage), with the flux scaled by a factor of 0.005 for display
purposes. We have masked regions that may suffer from sky line subtraction
residuals. The letters “r” and “b” following the quasar name indicate that we
are stacking the Ly$\alpha$ forest at the redshift corresponding to the red
and blue edges of the Ly$\beta$ gap, respectively. The stacked Ly$\alpha$
transmission profile is comparable with that expected for neutral islands with
$L=5.34h^{-1}\rm Mpc$ (dashed purple line) in Malloy & Lidz (2015).

##  4 Comparison with Simulations

Figure 2:  (a) Comparison of the stacked Ly$\alpha$ forest at the edges of
Ly$\beta$ dark gaps in mock spectra to the observations. The shaded regions
show the 68% intervals of model predictions from bootstrap trials. (b)
Distributions of the damping wing parameter, $v_{\rm ext}$, in the models. The
light and dark shaded horizontal bars in each subplot show the 68% and 95%
intervals, respectively. The vertical black line denotes $v_{\rm ext}$
acquired from the observation. Based on these comparisons, we can see that
simulations wherein dark gaps are substantially neutral are consistent with
the observed damping wing-like profile.  Figure 3:  (a) Similar to Figure 2,
but comparing the observed stacked Ly$\alpha$ forest with mock data that are
drawn from simulations with different overall neutral fraction from the
Sherwood-Relics patchy reionization simulation (Puchwein et al., 2023). Here,
we include _all_ dark gaps regardless of their individual neutral hydrogen
fraction. We display the shape of the mock stacked profiles from $z=5.4$ and
6.2 snapshots just for illustration purpose because their mean forest
transmission and dark gap length distribution are significantly different from
those at $z=5.8$, and different flux thresholds for dark gap detection are
used at these redshifts. (b) $v_{\rm ext}$ distribution from the patchy
reionization simulation at $z=5.8$. The observed damping wing-like profile is
consistent with the model prediction that has $\langle x_{\rm
HI}\rangle=7.4\%$.

A critical question is whether the observed damping wing-like profile in the
stacked Ly$\alpha$ forest indicates the presence of genuine neutral islands at
$z<6$. To explore the conditions under which such a profile might arise, we
generate mock spectra from three sets of simulations, create stacks based on
the Ly$\beta$ gaps in these spectra, and compare the properties of the mock
stacks to the observations.

In our first model, the sightlines are drawn from a radiative transfer (RT)
simulation conducted in a $(200h^{-1}\rm Mpc)^{3}$ volume box with $N=200^{3}$
RT cells, using the adaptive ray-tracing code described in Cain et al. (2021,
2023a, 2023b). Detailed methodology will be presented in C. Cain et al. (in
prep.). In the simulation, halos and galactic sources are populated based on
abundance matching according to the luminosity function from Finkelstein et
al. (2019), and the hydrogen ionizing emissivity of the sources is assumed to
be proportional to their UV luminosity. This simulation is calibrated to match
the observed mean transmission in the Ly$\alpha$ forest as reported in Bosman
et al. (2022). We account for the attenuation from the foreground Ly$\alpha$
forest when computing the Ly$\beta$ forest flux. Next, we adjust the modeled
spectra to have the same redshift coverage as the observed sample and add
noise according to the flux error array from the observations.

Dark gaps in the simulated sightlines are identified following the same
criteria as in Zhu et al. (2022); specifically, we select gaps that are dark
in both the Ly$\beta$ forest and the Ly$\alpha$ forest. We then stack the
Ly$\alpha$ forest at the edges of Ly$\beta$ dark gaps with lengths $L\geq
7{h^{-1}{\rm Mpc}}$. For our first model, we create stacked Ly$\alpha$ forest
by including dark gaps with a volume-weighted neutral hydrogen fraction of
$\geq 0.7$ measured in the pencil beam cells (denoted as “$x_{\rm
HI,\,gap}\geq 0.7$”), and dark gaps that contains only ionized gas (denoted as
“Ionized dark gap”, respectively. In practice, we exclude dark gaps with any
pencil beam cells that have volume-weighted neutral fraction greater than 0.5.
The fraction of “Ionized dark gap” is 29% out of all $L\geq 7{h^{-1}{\rm
Mpc}}$ Ly$\beta$ gaps, and the mean $x_{\rm HI,\,gap}$ is $2\times 10^{-3}$.

The second simulation we use here employs a fluctuating UVB model, wherein the
UVB fluctuations are boosted by a short mean free path and neutral islands are
not explicitly included (denoted as “ND20-early-shortmfp”; refer to Nasir &
D’Aloisio, 2020; Zhu et al., 2021 for details). In contrast to the radiative
transfer simulations of our first model, Nasir & D’Aloisio (2020) employ
approximate “semi-numeric” methods to model reionization’s effects on the
forest. For this model, we include all Ly$\beta$ dark gaps with $L\geq
7{h^{-1}{\rm Mpc}}$, except for those with all pencil beam cells that have
volume-averaged neutral fraction greater than 0.5. As a result, 99.88% of all
$L\geq 7{h^{-1}{\rm Mpc}}$ Ly$\beta$ gaps are selected, and the mean $x_{\rm
HI,\,gap}$ is $<0.03$. Here, this model serves as an independent reference for
dark gaps that are not dominated by neutral hydrogen islands.

In addition to the scenarios mentioned above, we include simulated sightlines
from a simulation of inhomogeneous reionization from the Sherwood-Relics suite
(described in detail in Puchwein et al., 2023) at different redshifts to
investigate the effects of a varying overall neutral hydrogen fraction. The
simulation we draw the sightlines from has a volume of 160 $h^{-1}$Mpc and was
performed with 20483 particles. An inhomogeneous UV background is included in
the simulation, taken from a pre-existing radiative transfer simulation
performed on snapshots of a cosmological simulation that used the same initial
conditions. This allows for spatial fluctuations in the ionization state of
the gas due to patchy reionization and the associated hydrodynamic response of
the gas to the photoheating. Reionization in the simulation completes at
$z=5.3$. We draw mock spectra from $z=5.4$, 5.8, and 6.2, which have overall
neutral fraction of $\langle x_{\rm HI}\rangle=0.4\%$, 7.4%, and 21.3%,
respectively.

We generate mock datasets from the simulations to test the likelihood that a
profile similar to the data will emerge in small samples. To mitigate the
effect of variations in dark gap length distributions across different cases,
in each bootstrap realization, we randomly draw dark gaps from the
simulations, following exactly the same distribution of dark gap lengths as in
the observed data to create the stack. Nevertheless, the shape of the stacked
Ly$\alpha$ forest does not change significantly even without the correction on
the dark gaps length distribution.

In the Figure 2, we compare the stacked Ly$\alpha$ forest between observations
and the mock stacks from a fixed redshift, i.e., $z=5.8$. As expected, the
stack of dark gaps containing a significant fraction of neutral gas exhibits
an extended and damping wing-like profile in the Ly$\alpha$ forest. The
observed curve broadly aligns with the model profile and falls within the 68%
interval. On the other hand, the “ionized dark gap” case shows a rapid
increase in the stacked Ly$\alpha$ transmission, inconsistent with the
observed damping wing-like profile. The “ND20-early-shortmfp” case also
displays a steep increase in the stacked Ly$\alpha$ transmission, albeit with
a slightly lower maximum flux compared to the “ionized dark gap” cases.

For a quantitative comparison, we fit the upper envelope of the model-
predicted stacked Ly$\alpha$ transmission in each bootstrap realization to
Equation 1 and plot the distribution of the damping wing parameter $v_{\rm
ext}$ in Figure 2(b). Only the scenario where dark gaps contain a large
fraction of neutral gas predicts a broad distribution of $v_{\rm ext}$,
spanning the observed value within the 68% interval. Conversely, scenarios
that fail to reproduce an extended absorption profile predict a very narrow
$v_{\rm ext}$ distribution, confined to $v_{\rm ext}\lesssim 100$ km s-1, with
the observed $v_{\rm ext}$ lying beyond their 95% confidence limits. These
results suggest that the observed damping wing-like profile in the stacked
Ly$\alpha$ transmission likely originates from genuine neutral IGM absorption
at $z<6$, as opposed to fully ionized gas that is highly opaque to due to a
fluctuating UVB or high densities.

As an additional test, we show the redshift evolution of this signal with mock
stacks of spectra created from snapshots of one of the Sherwood-Relics
simulations as described above. For this test we include all dark gaps present
in the simulated lines of sight, without any pre-selection based on their
individual neutral fractions. The results are shown in Figure 3. As the
volume-averaged ionization fraction of the IGM evolves, we observe changes in
both the mean flux of the Ly$\alpha$ forest and the shape of the Ly$\alpha$
profile. For the snapshot at $z=5.4$, when the volume-averaged HI fraction is
0.4%, the Ly$\alpha$ profile is consistent with the “ionized dark gap”
scenario, where the Ly$\alpha$ transmission is boosted on the edge of
Ly$\beta$ dark gaps. With increasing redshift and hence an increasing $\langle
x_{\rm HI}\rangle$, the flux boost becomes smoother and finally approaches a
damping wing-like shape. For the snapshot at $z=5.8$, which has a Ly$\alpha$
forest mean flux consistent with the observational sample, the distribution of
$v_{e}xt$ is consistent with the observation within the 95% limits, possibly
because the neutral fraction is just above the lower limit constraint (see
Section 5). This shows that simultaneously fitting for the mean flux of the
Ly$\alpha$ forest as well as the shape of the Ly$\alpha$ profile at the edge
of Ly$\beta$ dark gaps can be a powerful probe of the end stages of
reionization.

Figure 4: Stacked (mean) Ly$\alpha$ transmission profiles at the redshift of
the edges of Ly$\beta$ dark gaps detected in ideal sightlines from our first
simulation. Panels show (a) the Ly$\alpha$ transmission, (b) the ionization
rate, (c) the IGM temperature, (d) the overdensity, and (e) the volume filling
factor of neutral hydrogen. The solid blue line represents the “$x_{\rm
HI,\,gap}\geq 0.7$” case, where dark gaps are at least $70\%$ neutral. The
solid orange line depicts the case of “ionized dark gaps”. Dotted lines
indicate the stacked profiles when Ly$\beta$ dark gaps are selected without
contamination from foreground Ly$\alpha$ absorption.

We further investigate the physical properties the simulated Ly$\beta$ dark
gaps by dividing the sample based on their neutral fraction. For our first
model, Figure 4 displays the stacked properties of the Ly$\beta$ dark gaps at
the redshift of their edges. As the figure shows, dark gaps with substantial
neutral hydrogen produce an extended absorption profile in the adjacent
Ly$\alpha$ forest. These dark gaps are typically found in regions
characterized by a low ionizing UVB, low temperature, and low densities. Their
properties align with findings in studies such as Gnedin (2022).

Notably, we find that if Ly$\beta$ dark gaps are selected in the absence of
contamination from foreground Ly$\alpha$ absorption, then even the “ionized
dark gap” scenario exhibits extended Ly$\alpha$ absorption. Physically, this
is because the Ly$\alpha$ transmission should recover gradually when moving
away from ionized regions of high Ly$\alpha$ opacity. Contamination from
Ly$\beta$ gaps that are largely created by foreground Ly$\alpha$ absorption
obscures this effect. We find that up to $41\%$ of the total Ly$\beta$ dark
gap length in our first model can arise from the foreground contamination. As
for the “$x_{\rm HI,\,gap}\geq 0.7$” scenario, the selection criteria ensure
that the dark gaps contain genuine neutral islands that are not highly
contaminated. As discussed above, the observed extended absorption profile may
indicate that the long dark gaps in Zhu et al. (2022) are similar to simulated
dark gaps that are created by neutral islands.

Figure 5: Lower limit constraint on the volume-averaged neutral hydrogen
fraction from the damping wing-like profile in the stacked Ly$\alpha$ forest
based on Ly$\beta$ dark gaps over $5.6<z<6.0$. We also plot $\langle x_{\rm
HI}\rangle$ constraints from recent observations for reference: upper limit
constraints from Ly$\alpha$ \+ Ly$\beta$ dark gap statistics (Zhu et al.,
2022), upper limit constraints from the fraction of Ly$\alpha$ \+ Ly$\beta$
dark pixels (Jin et al., 2023), measurements from high-$z$ quasar damping
wings (Greig et al., 2017, 2019, 2022; Bañados et al., 2018; Davies et al.,
2018; Wang et al., 2020; Yang et al., 2020a; Ďurovčíková et al., 2024; Greig
et al., 2024), constraints based on Ly$\alpha$ opacities(Gaikwad et al.,
2023), and measurements from the Ly$\alpha$ damping wing of high-$z$ galaxies
(Umeda et al., 2023). We also include the measurements based on the stacked
Ly$\alpha$ forest using a different method in B. Spina et al. (in prep.). The
orange shaded region plot the posterior from inference based on multiple
reionization observables (Y. Qin et al. in prep.). Data points in this figure
are slightly displaced along $z$ for display purpose.

##  5 Lower Limit Constraint on Neutral Hydrogen Fraction

We can use the observed damping wing-like profile to constrain the hydrogen
neutral fraction in unbiased locations within the IGM. To do this, we generate
dark gaps following the observed gap length distribution for different ranges
of $x_{\rm HI,gap}$, specifically $0.1\leq x_{\rm HI,gap}<0.2$, $0.2\leq
x_{\rm HI,gap}<0.3$, …, up to $0.9\leq x_{\rm HI,gap}<1.0$, and create mock
Ly$\alpha$ forest stacks for each range using the RT simulation used for our
first model. We find that only $x_{\rm HI,gap}$ values greater than 0.4 can
predict $v_{\rm ext}$ consistent with the observed value within the 95% limit.
After convolving $x_{\rm HI,gap}\geq 0.4$ with the fraction of the Ly$\alpha$
forest subtended by dark gaps longer than the mean length of dark gaps used in
this work, $f_{L\geq 15{h^{-1}{\rm Mpc}}}=0.152\pm 0.098$, (Zhu et al., 2022)
we derive a lower limit on the volume-averaged neutral hydrogen fraction of
$\langle x_{\rm HI}\rangle\geq 0.061\pm 0.039$ at $z=5.8$. This constraint is
conservative in that it excludes potential neutral hydrogen in dark gaps with
$L<7{h^{-1}{\rm Mpc}}$. We also emphasize that the simulations we are using
include significant UVB fluctuations, which helps us to account for the
possibility that a fraction of the dark gaps will be produced by ionized gas.

Figure 5 displays our constraint alongside other recent $\langle x_{\rm
HI}\rangle$ measurements. Our constraint on $\langle x_{\rm HI}\rangle$
suggests that the IGM must be at least 6.1% neutral even down to $z=5.8$.
Alongside upper limit constraints from dark gaps and dark pixels (Zhu et al.,
2022; Jin et al., 2023), these observations accommodate a wide range of late-
ending reionization histories, with $\langle x_{\rm HI}\rangle$ constraints
ranging from approximately 6% to 20% near $z=5.8$. We also note that the
inference based on multiple recent reionization observables (Y. Qin et al. in
prep.), which includes the mean transmission in the Ly$\alpha$ forest, galaxy
UV luminosity functions, and CMB observations, aligns well with these $\langle
x_{\rm HI}\rangle$ constraints.

##  6 Summary

We have identified a damping wing-like profile in the stacked Ly$\alpha$
forest associated with Ly$\beta$ dark gaps at $z<6$, using data from the
XQR-30 program and archival spectra. A comparison with mock spectra generated
from simulations suggests that this damping wing-like profile most likely
arises from dark gaps containing neutral hydrogen islands, as opposed to fully
ionized regions whose high opacities are the result of high densities and/or a
locally low UVB. We further ascertain that such a profile necessitates a
neutral hydrogen fraction of at least $6.1\pm 3.9\%$ at $z=5.8$. Combined with
the IGM damping wing discovered toward an individual quasar sightline at
$z=5.878$ (see the companion paper, Becker et al., 2024), our results provide
robust evidence that reionization extends to $z<6$.

## Acknowledgements

Y.Z. and G.D.B. were supported by the National Science Foundation through
grant AST-1751404. Y.Z. was also supported by the NSF through award
SOSPADA-029 from the NRAO. G.K. and M.G.H. have been partially supported by
the DAE-STFC project “Building Indo-UK collaborations towards the Square
Kilometre Array” (STFC grant reference ST/Y004191/1). Parts of this research
were supported by the Australian Research Council Centre of Excellence for All
Sky Astrophysics in 3 Dimensions (ASTRO 3D), through project no. CE170100013.
B.S. and S.E.I.B. are supported by the Deutsche Forschungsge-meinschaft (DFG)
under Emmy Noether grant number BO 5771/1-1. F.W. acknowledges support from
NSF Grant AST-2308258. For the purpose of open access, the authors have
applied a Creative Commons Attribution (CC BY) licence to any Author Accepted
Manuscript version arising from this submission.

This work is based on observations collected at the European Southern
Observatory under ESO programme 1103.A-0817. Some of the data presented herein
were obtained at Keck Observatory, which is a private 501(c)3 non-profit
organization operated as a scientific partnership among the California
Institute of Technology, the University of California, and the National
Aeronautics and Space Administration. The Observatory was made possible by the
generous financial support of the W. M. Keck Foundation. The authors wish to
recognize and acknowledge the very significant cultural role and reverence
that the summit of Maunakea has always had within the Native Hawaiian
community. We are most fortunate to have the opportunity to conduct
observations from this mountain. This research has made use of the Keck
Observatory Archive (KOA), which is operated by the W. M. Keck Observatory and
the NASA Exoplanet Science Institute (NExScI), under contract with the
National Aeronautics and Space Administration.

We respectfully acknowledge the University of Arizona is on the land and
territories of Indigenous peoples. Today, Arizona is home to 22 federally
recognized tribes, with Tucson being home to the O’odham and the Yaqui.
Committed to diversity and inclusion, the University strives to build
sustainable relationships with sovereign Native Nations and Indigenous
communities through education offerings, partnerships, and community service.

The simulations used in this work were performed using the Joliot Curie
supercomputer at the Très Grand Centre de Calcul (TGCC) and the Cambridge
Service for Data Driven Discovery (CSD3), part of which is operated by the
University of Cambridge Research Computing on behalf of the STFC DiRAC HPC
Facility (www.dirac.ac.uk). We acknowledge the Partnership for Advanced
Computing in Europe (PRACE) for awarding us time on Joliot Curie in the 16th
call. The DiRAC component of CSD3 was funded by BEIS capital funding via STFC
capital grants ST/P002307/1 and ST/R002452/1 and STFC operations grant
ST/R00689X/1. This work also used the DiRAC@Durham facility managed by the
Institute for Computational Cosmology on behalf of the STFC DiRAC HPC
Facility. The equipment was funded by BEIS capital funding via STFC capital
grants ST/P002293/1 and ST/R002371/1, Durham University and STFC operations
grant ST/R000832/1. DiRAC is part of the National e-Infrastructure.

Software: Astropy (Astropy Collaboration et al., 2013), Matplotlib (Hunter,
2007), NumPy (van der Walt et al., 2011), SpectRes (Carnall, 2017).

## Data Availability

The raw data used in this work are available from the ESO archive at
http://archive.eso.org and KOA at https://koa.ipac.caltech.edu. The reduced
X-Shooter spectra are available through a public GitHub repository of the
XQR-30 spectra at https://github.com/XQR-30/Spectra. The reduced Keck spectra
can be obtained from the corresponding author upon reasonable request.

## References

  * Asthana et al. (2024) Asthana S., Haehnelt M. G., Kulkarni G., Aubert D., Bolton J. S., Keating L. C., 2024, Late-End Reionization with ATON-HE: Towards Constraints from Lyman-$\alpha$ Emitters Observed with JWST, doi:10.48550/arXiv.2404.06548
  * Astropy Collaboration et al. (2013) Astropy Collaboration et al., 2013, A&A, 558, A33
  * Bañados et al. (2018) Bañados E., et al., 2018, Nature, 553, 473
  * Becker et al. (2015) Becker G. D., Bolton J. S., Madau P., Pettini M., Ryan-Weber E. V., Venemans B. P., 2015, MNRAS, 447, 3402
  * 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
  * Becker et al. (2024) Becker G. D., Bolton J. S., Zhu Y., Hashemi S., 2024, Damping Wing Absorption Associated with a Giant Ly$\alpha$ Trough at $z $<$ 6$: Direct Evidence for Late-Ending Reionization, doi:10.48550/arXiv.2405.08885
  * Bosman (2021) Bosman S. E. I., 2021, Constraints on the Mean Free Path of Ionising Photons at $z\sim6$ Using Limits on Individual Free Paths, doi:10.48550/arXiv.2108.12446
  * 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. (2021) Bosman S. E. I., Ďurovčíková D., Davies F. B., Eilers A. C., 2021, MNRAS, 503, 2077
  * Bosman et al. (2022) Bosman S. E. I., et al., 2022, MNRAS, 514, 55
  * Bouwens et al. (2015) Bouwens R. J., Illingworth G. D., Oesch P. A., Caruana J., Holwerda B., Smit R., Wilkins S., 2015, ApJ, 811, 140
  * Cain et al. (2021) Cain C., D’Aloisio A., Gangolli N., Becker G. D., 2021, ApJ, 917, L37
  * Cain et al. (2023a) Cain C., D’Aloisio A., Lopez G., Gangolli N., Roth J. T., 2023a, On the Rise and Fall of Galactic Ionizing Output at the End of Reionization, doi:10.48550/arXiv.2311.13638
  * Cain et al. (2023b) Cain C., D’Aloisio A., Gangolli N., McQuinn M., 2023b, MNRAS, 522, 2047
  * Carnall (2017) Carnall A. C., 2017, arXiv e-prints, p. arXiv:1705.05165
  * 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
  * Christenson et al. (2023) Christenson H. M., et al., 2023, ApJ, 955, 138
  * Cooper et al. (2019) Cooper T. J., Simcoe R. A., Cooksey K. L., Bordoloi R., Miller D. R., Furesz G., Turner M. L., Bañados E., 2019, ApJ, 882, 77
  * D’Odorico et al. (2023) D’Odorico V., et al., 2023, MNRAS, 523, 1399
  * Davies et al. (2018) Davies F. B., et al., 2018, ApJ, 864, 142
  * Davies et al. (2023a) Davies R. L., et al., 2023a, MNRAS, 521, 289
  * Davies et al. (2023b) Davies R. L., et al., 2023b, MNRAS, 521, 314
  * Davies et al. (2024) Davies F. B., et al., 2024, ApJ, 965, 134
  * Decarli et al. (2018) Decarli R., et al., 2018, ApJ, 854, 97
  * Ďurovčíková et al. (2024) Ďurovčíková D., et al., 2024, Chronicling the Reionization History at $6\lesssim z \lesssim 7$ with Emergent Quasar Damping Wings, doi:10.48550/arXiv.2401.10328
  * Eilers et al. (2018) Eilers A.-C., Davies F. B., Hennawi J. F., 2018, ApJ, 864, 53
  * Eilers et al. (2020) Eilers A.-C., et al., 2020, ApJ, 900, 37
  * Fan et al. (2006) Fan X., et al., 2006, AJ, 132, 117
  * Fan et al. (2023) Fan X., Bañados E., Simcoe R. A., 2023, ARA&A, 61, 373
  * Finkelstein et al. (2019) Finkelstein S. L., et al., 2019, ApJ, 879, 36
  * Gaikwad et al. (2023) Gaikwad P., et al., 2023, MNRAS, 525, 4093
  * Gnedin (2022) Gnedin N. Y., 2022, ApJ, 937, 17
  * 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
  * Greig et al. (2024) Greig B., et al., 2024, MNRAS, 530, 3208
  * Hunter (2007) Hunter J. D., 2007, CSE, 9, 90
  * Jin et al. (2023) Jin X., et al., 2023, ApJ, 942, 59
  * Kashino et al. (2020) Kashino D., Lilly S. J., Shibuya T., Ouchi M., Kashikawa N., 2020, ApJ, 888, 6
  * Keating et al. (2020) Keating L. C., Weinberger L. H., Kulkarni G., Haehnelt M. G., Chardin J., Aubert D., 2020, MNRAS, 491, 1736
  * 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
  * Malloy & Lidz (2015) Malloy M., Lidz A., 2015, ApJ, 799, 179
  * McGreer et al. (2013) McGreer I. D., et al., 2013, ApJ, 768, 105
  * Miralda-Escudé (1998) Miralda-Escudé J., 1998, ApJ, 501, 15
  * Miralda-Escudé & Rees (1998) Miralda-Escudé J., Rees M. J., 1998, ApJ, 497, 21
  * Muñoz et al. (2024) Muñoz J. B., Mirocha J., Chisholm J., Furlanetto S. R., Mason C., 2024, Reionization after JWST: A Photon Budget Crisis?, doi:10.48550/arXiv.2404.07250
  * Nasir & D’Aloisio (2020) Nasir F., D’Aloisio A., 2020, MNRAS, 494, 3080
  * Planck Collaboration et al. (2014) Planck Collaboration et al., 2014, A&A, 571, A16
  * Planck Collaboration et al. (2020) Planck Collaboration et al., 2020, A&A, 641, A6
  * Prochaska et al. (2003) Prochaska J. X., Gawiser E., Wolfe A. M., Cooke J., Gelino D., 2003, ApJS, 147, 227
  * Puchwein et al. (2023) Puchwein E., et al., 2023, MNRAS, 519, 6162
  * Qin et al. (2021) Qin Y., Mesinger A., Bosman S. E. I., Viel M., 2021, MNRAS, 506, 2390
  * Robertson et al. (2015) Robertson B. E., Ellis R. S., Furlanetto S. R., Dunlop J. S., 2015, ApJ, 802, L19
  * Roth et al. (2024) Roth J. T., D’Aloisio A., Cain C., Wilson B., Zhu Y., Becker G. D., 2024, MNRAS, 530, 5209
  * Satyavolu et al. (2023) Satyavolu S., Kulkarni G., Keating L. C., Haehnelt M. G., 2023, Robustness of Direct Measurements of the Mean Free Path of Ionizing Photons in the Epoch of Reionization (arxiv:2311.06344)
  * Sebastian et al. (2024) Sebastian A. M., et al., 2024, MNRAS, 530, 1829
  * Sheinis et al. (2002) Sheinis A. I., Bolte M., Epps H. W., Kibrick R. I., Miller J. S., Radovan M. V., Bigelow B. C., Sutin B. M., 2002, PASP, 114, 851
  * Shen et al. (2019) Shen Y., et al., 2019, ApJ, 873, 35
  * Stark (2016) Stark D. P., 2016, ARA&A, 54, 761
  * Umeda et al. (2023) Umeda H., Ouchi M., Nakajima K., Harikane Y., Ono Y., Xu Y., Isobe Y., Zhang Y., 2023, JWST Measurements of Neutral Hydrogen Fractions and Ionized Bubble Sizes at $z=7-12$ Obtained with Ly$\alpha$ Damping Wing Absorptions in 26 Bright Continuum Galaxies (arxiv:2306.00487)
  * Vernet et al. (2011) Vernet J., et al., 2011, A&A, 536, A105
  * Wang et al. (2010) Wang R., et al., 2010, ApJ, 714, 699
  * Wang et al. (2013) Wang R., et al., 2013, ApJ, 773, 44
  * Wang et al. (2020) Wang F., et al., 2020, ApJ, 896, 23
  * Weinberger et al. (2019) Weinberger L. H., Haehnelt M. G., Kulkarni G., 2019, MNRAS, 485, 1350
  * Yang et al. (2020a) Yang J., et al., 2020a, ApJ, 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
  * de Belsunce et al. (2021) de Belsunce R., Gratton S., Coulton W., Efstathiou G., 2021, MNRAS, 507, 1072
  * van der Walt et al. (2011) van der Walt S., Colbert S. C., Varoquaux G., 2011, CSE, 13, 22


64 | 142 / 1424.39 | 156 / 3197.97 | 172 / 6781.49 | 189 / 14399.77
128 | da fare … | da fare … | $\ast$ | $\ast$

Table 9: Rotated quarter domain. Performance of $\widehat{\mathbf{B}}$,
$\widehat{\mathbf{B}}^{\boldsymbol{G}}$ and $\widehat{\mathbf{Q}}$.

##  11 Conclusions

In this work we proposed a review of several preconditioners suited for space-
time Galerkin isogeometric discretizations of the heat equation. Our
preconditioners are represented by a suitable sum of Kronecker products of
matrices, that makes the computational cost of their construction (setup) and
application, as well as the storage cost, very appealing. In particular,
inspired by the FD technique, the application of the preconditioner
$\widehat{\mathbf{A}}$ exploits an ad-hoc factorization of the time matrices.
The preconditioners $\widehat{\mathbf{P}}$ and $\widehat{\mathbf{Q}}$
factorize in time direction as sum of diagonal blocks and rank-1 block
perturbations, while in space have a diagonal structure. Lastly, the
preconditioner $\widehat{\mathbf{B}}$ admits a full diagonal factorization.

The application cost of all the preconditioners is almost equal to
$O(N_{dof})$ and does not depend on the polynomial degree.

At the same time, the storage cost is roughly the same that we would have by
discretizing separately in space and in time, if we assume $n_{t}\leq
Cp^{d}N_{s}$. Indeed, in this case the memory used for the whole iterative
solver is $O(p^{d}N_{s}+N_{dof}).$

In this review, we have compared the performance of the four preconditioners
$\widehat{\mathbf{A}},\widehat{\mathbf{B}},\widehat{\mathbf{P}}$ and
$\widehat{\mathbf{Q}}$, first on the parametric 3D cube spacial domain, and
then on a rotated quarter of ring. In the first computational domain
$\widehat{\mathbf{A}}$ and $\widehat{\mathbf{P}}$ are direct solvers and
$\widehat{\mathbf{Q}}$ seem to converge, for finer meshes and higher degrees,
to a direct solver. The preconditioner $\widehat{\mathbf{B}}$ is although
stable under mesh refinement and degree elevation.

On the rotated quarter of ring geometry, the performances of
$\widehat{\mathbf{A}}^{\boldsymbol{G}}$ and
$\widehat{\mathbf{B}}^{\boldsymbol{G}}$, that are the versions of
$\widehat{\mathbf{A}}$ and $\widehat{\mathbf{B}}$ with a partial inclusion of
the geometry’s information, are outstanding when compared to their parametric
versions. Thus the inclusion of information for the preconditioners
$\widehat{\mathbf{P}}$ and $\widehat{\mathbf{Q}}$ seems to be a natural way
forward as a future development of this work.

As a final comment, we mention that our methods has a strong potential for
parallelization, and this will be an interesting future direction of study.

## Acknowledgments

## References

  * [1] T. J. R. Hughes, J. A. Cottrell, Y. Bazilevs, Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement, Computer Methods in Applied Mechanics and Engineering 194 (39) (2005) 4135–4195.
  * [2] J. A. Cottrell, T. J. R. Hughes, Y. Bazilevs, Isogeometric analysis: toward integration of CAD and FEA, John Wiley & Sons, 2009.
  * [3] J. A. Evans, Y. Bazilevs, I. Babuška, T. J. R. Hughes, $n$-widths, sup-infs, and optimality ratios for the $k$-version of the isogeometic finite element method, Computer Methods in Applied Mechanics and Engineering 198 (2009) 1726–1741.
  * [4] A. Bressan, E. Sande, Approximation in FEM, DG and IGA: a theoretical comparison, Numerische Mathematik (2019).
  * [5] G. Sangalli, M. Tani, Matrix-free weighted quadrature for a computationally efficient isogeometric $k$-method, Computer Methods in Applied Mechanics and Engineering 338 (2018) 117 – 133.
  * [6] I. Fried, Finite-element analysis of time-dependent phenomena., AIAA Journal 7 (6) (1969) 1170–1173.
  * [7] J. C. Bruch Jr., G. Zyvoloski, Transient two-dimensional heat conduction problems solved by the finite element method, International Journal for Numerical Methods in Engineering 8 (3) (1974) 481–494.
  * [8] J. T. Oden, A general theory of finite elements. II. Applications, International Journal for Numerical Methods in Engineering 1 (3) (1969) 247–259.
  * [9] F. Shakib, T. J. R. Hughes, A new finite element formulation for computational fluid dynamics: IX. Fourier analysis of space-time Galerkin/least-squares algorithms, Computer Methods in Applied Mechanics and Engineering 87 (1) (1991) 35–58.
  * [10] U. Langer, S. E. Moore, M. Neumüller, Space–time isogeometric analysis of parabolic evolution problems, Computer Methods in Applied Mechanics and Engineering 306 (2016) 342 – 363.
  * [11] U. Langer, M. Neumüller, I. Toulopoulos, Multipatch space-time isogeometric analysis of parabolic diffusion problems, in: International Conference on Large-Scale Scientific Computing, Springer, 2017, pp. 21–32.
  * [12] K. Takizawa, T. E. Tezduyar, Space–time computation techniques with continuous representation in time (st-c), Computational Mechanics 53 (1) (2014) 91–99.
  * [13] Y. Ueda, N. Saito, Stability and error estimates for the successive-projection technique with B-splines in time, Journal of Computational and Applied Mathematics 358 (2019) 266 – 278.
  * [14] M. J. Gander, M. Neumüller, Analysis of a new space-time parallel multigrid algorithm for parabolic problems, SIAM Journal on Scientific Computing 38 (4) (2016) A2173–A2208.
  * [15] C. Hofer, U. Langer, M. Neumüller, R. Schneckenleitner, Parallel and robust preconditioning for space-time isogeometric analysis of parabolic evolution problems, SIAM Journal on Scientific Computing 41 (3) (2019) A1793–A1821.
  * [16] A. Mantzaflaris, F. Scholz, I. Toulopoulos, Low-rank space-time decoupled isogeometric analysis for parabolic problems with varying coefficients, Computational Methods in Applied Mathematics 19 (1) (2019) 123–136.
  * [17] J. Bonilla, S. Badia, Maximum-principle preserving space–time isogeometric analysis, Computer Methods in Applied Mechanics and Engineering 354 (2019) 422–440.
  * [18] K. Takizawa, T. E. Tezduyar, Y. Otoguro, T. Terahara, T. Kuraishi, H. Hattori, Turbocharger flow computations with the space–time isogeometric analysis (ST-IGA), Computers & Fluids 142 (2017) 15–20.
  * [19] K. Takizawa, T. E. Tezduyar, T. Terahara, Ram-air parachute structural and fluid mechanics computations with the space–time isogeometric analysis (ST-IGA), Computers & Fluids 141 (2016) 191–200.
  * [20] K. Takizawa, T. E. Tezduyar, T. Terahara, T. Sasaki, Heart valve flow computation with the space–time slip interface topology change (ST-SI-TC) method and isogeometric analysis (IGA), in: Biomedical Technology, Springer, 2018, pp. 77–99.
  * [21] M. Montardini, M. Negri, G. Sangalli, M. Tani, Space-time least-squares isogeometric method and efficient solver for parabolic problems, Mathematics of Computation (accepted for publication) (2019).
  * [22] G. Loli, M. Montardini, G. Sangalli, M. Tani, An efficient solver for space–time isogeometric Galerkin methods for parabolic problems, Computers & Mathematics with Applications 80 (11) (2020) 2586–2603.
  * [23] O. Steinbach, Space-time finite element methods for parabolic problems, Computational Methods in Applied Mathematics 15 (4) (2015) 551–566.
  * [24] R. Stevenson, J. Westerdiep, Stability of Galerkin discretizations of a mixed space-time variational formulation of parabolic evolution equations, arXiv:1902.06279 (2019).
  * [25] R. E. Lynch, J. R. Rice, D. H. Thomas, Direct solution of partial difference equations by tensor product methods, Numerische Mathematik 6 (1) (1964) 185–199.
  * [26] G. Sangalli, M. Tani, Isogeometric preconditioners based on fast solvers for the Sylvester equation, SIAM Journal on Scientific Computing 38 (6) (2016) A3644–A3671.
  * [27] C. A. Dorao, H. A. Jakobsen, A parallel time–space least-squares spectral element solver for incompressible flow problems, Applied Mathematics and Computation 185 (1) (2007) 45–58.
  * [28] M. J. Gander, 50 years of time parallel time integration, in: Multiple Shooting and Time Domain Decomposition Methods, Springer, 2015, pp. 69–113.
  * [29] A. M. Kvarving, E. M. Rønquist, A fast tensor-product solver for incompressible fluid flow in partially deformed three-dimensional domains: Parallel implementation, Computers & Fluids 52 (2011) 22–32.
  * [30] C. De Boor, A practical guide to splines (revised edition), Applied Mathematical Sciences, Springer, Berlin, 2001.
  * [31] T. G. Kolda, B. W. Bader, Tensor decompositions and applications, SIAM review 51 (3) (2009) 455–500.
  * [32] L. C. Evans, Partial Differential equations, American Mathematical Society, Berlin, 2010.
  * [33] L. Beirão da Veiga, D. Cho, G. Sangalli, Anisotropic NURBS approximation in isogeometric analysis, Computer Methods in Applied Mechanics and Engineering 209 (2012) 1–11.
  * [34] M. O. Deville, P. F. Fischer, E. H. Mund, High-order methods for incompressible fluid flow, Cambridge University Press, 2002.
  * [35] K. P. S. Gahalaut, S. K. Tomar, C. Douglas, Condition number estimates for matrices arising in NURBS based isogeometric discretizations of elliptic partial differential equations, arXiv preprint arXiv:1406.6808 (2014).
  * [36] M. Montardini, G. Sangalli, M. Tani, Robust isogeometric preconditioners for the Stokes system based on the Fast Diagonalization method, Computer Methods in Applied Mechanics and Engineering 338 (2018) 162 – 185.
  * [37] A. Tkachuk, M. Bischoff, Direct and sparse construction of consistent inverse mass matrices: general variational formulation and application to selective mass scaling, International Journal for Numerical Methods in Engineering 101 (6) (2015) 435–469.
  * [38] L. Wunderlich, A. Seitz, M. D. Alaydın, B. Wohlmuth, A. Popp, Biorthogonal splines for optimal weak patch-coupling in isogeometric analysis with applications to finite deformation elasticity, Computer Methods in Applied Mechanics and Engineering 346 (2019) 197–215.
  * [39] L. Gao, V. M. Calo, Fast isogeometric solvers for explicit dynamics, Computer Methods in Applied Mechanics and Engineering 274 (2014) 19–41.
  * [40] J. Chan, J. A. Evans, Multi-patch discontinuous galerkin isogeometric analysis for wave propagation: Explicit time-stepping and efficient mass matrix inversion, Computer Methods in Applied Mechanics and Engineering 333 (2018) 22–54.
  * [41] G. Loli, G. Sangalli, M. Tani, Easy and efficient preconditioning of the isogeometric mass matrix, Computers & Mathematics with Applications 116 (2022) 245–264.
  * [42] A. Mantzaflaris, B. Jüttler, B. N. Khoromskij, U. Langer, Matrix generation in isogeometric analysis by low rank tensor approximation, in: Curves and Surfaces: 8th International Conference, Paris, France, June 12-18, 2014, Revised Selected Papers 8, Springer, 2015, pp. 321–340.
  * [43] A. Mantzaflaris, B. Jüttler, B. N. Khoromskij, U. Langer, Low rank tensor methods in galerkin-based isogeometric analysis, Computer Methods in Applied Mechanics and Engineering 316 (2017) 1062–1085.
  * [44] C. Hofreither, A black-box low-rank approximation algorithm for fast matrix assembly in isogeometric analysis, Computer Methods in Applied Mechanics and Engineering 333 (2018) 311–330.
  * [45] C. Schwab, R. Stevenson, Space-time adaptive wavelet methods for parabolic evolution problems, Mathematics of Computation 78 (267) (2009) 1293–1318.
  * [46] J. Henning, D. Palitta, V. Simoncini, K. Urban, An ultraweak space-time variational formulation for the wave equation: Analysis and efficient numerical solution, ESAIM: Mathematical Modelling and Numerical Analysis 56 (4) (2022) 1173–1198.
  * [47] R. Vázquez, A new design for the implementation of isogeometric analysis in Octave and Matlab: GeoPDEs 3.0, Computers & Mathematics with Applications 72 (3) (2016) 523–554.
  * [48] L. Sorber, M. Van Barel, L. De Lathauwer, Tensorlab v2. 0, Available online, URL: www.tensorlab.net (2014).