Behavior of the Lyα Damping Wings as a Function of Reionization Topology

After the Universe's rapid expansion following the Big Bang, it began to cool down, leading to the recombination of free electrons and protons to form neutral hydrogen and helium. Once recombination ceased, the Universe entered a period known as the "Dark Ages," during which it continued to expand and cool. This period continued until the Universe underwent a phase transition: neutral hydrogen began to reionize, marking the Epoch of Reionization (EoR).

The first stars, galaxies, and black holes, as the strong sources of ultraviolet radiation, initiated the reionization of the surrounding intergalactic medium (IGM), a state that has persisted to the present day. Although the precise timing of the reionization epoch remains uncertain, it is estimated to have begun around a redshift of z ∼ 15 (≈250 million years after the Big Bang) and ended at z ∼ 5.5 (≈1 billion years after the Big Bang) (S. E. Bosman et al. 2022; P. Gaikwad et al. 2023; F. B. Davies et al. 2024). During and before this epoch, the Universe was predominantly composed of neutral hydrogen, making hydrogen emission and absorption features a key source of observational data. The two most informative signals are the Lyα forest and the 21 cm line emission.

However, we currently lack direct detections of the 21 cm signal (e.g., C. M. Trott et al. 2020; HERA Collaboration et al. 2023; A. Acharya et al. 2024), and studies have shown that the Lyα forest provides only approximate information about the end of reionization (S. E. Bosman et al. 2022). Due to the large resonant cross-section of Lyα, nearly all the flux is absorbed even when the neutral hydrogen fraction (x_{H I}) is as low as 10⁻⁴ (J. E. Gunn & B. A. Peterson 1965). However, the Lyα absorption profile also exhibits broad tails extending to both sides of the peak. Hence, at a high neutral fraction, we can observe these tails far from the peak, the so-called "damping wing," which can be observed redward of rest-frame Lyα in the spectra of reionization-epoch objects (J. Miralda-Escudé 1998).

Existing detections of this damping wing signal in the spectra of the most distant quasars (D. J. Mortlock et al. 2011; B. Greig et al. 2017; F. B. Davies et al. 2018; F. Wang et al. 2020; B. Greig et al. 2022) showed that the Universe at z ∼ 7.0–7.5 is about half neutral. These measurements are in agreement with complementary constraints from the evolution of Lyα emission lines from galaxies, which are suppressed by the same damping wing absorption (e.g., C. A. Mason et al. 2018; A. M. Morales et al. 2021).

At present, the best constraints on reionization come from the statistical comparison of quasar damping wings with suites of cosmological simulations boxes with a range of neutral fractions (e.g., B. Greig et al. 2017; F. B. Davies et al. 2018; C. A. Mason et al. 2018); however, except for the discussion in B. Greig et al. (2019), these neutral fraction constraints from damping wings never considered the variation of the model for the ionizing emission from galaxies.

In B. Greig et al. (2019) it was hinted that the choice of the source model could shift the neutral fraction inference. However, it was concluded that it may not be important when considering the damping wings from individual quasars, as the uncertainty of the neutral fraction—largely due to cosmic variance—was large.

In addition to introducing a potential systematic bias in the inferred neutral fraction, the standard method of using a single source model also fails to incorporate the effect of reionization topology on the intrinsic scatter of the damping wing signal. This cosmic variance may become an effective tool to constrain reionization parameters when an ensemble of damping wings is considered (D. Ďurovčíková et al. 2024). Hence, there is a need for a comprehensive analysis of the damping wing signal with a wider choice of parameters encompassing the effects of source and astrophysical parameters on the reionization topology.

Our study addresses this challenge by simulating ensembles of damping wing profiles. The aim is to examine the dependence of the damping wing profile and the sightline-to-sightline scatter on a set of source and astrophysical parameters, demonstrating how variations in reionization topology due to changes in these parameters can be observed in the statistical properties of the damping profiles.

In Section 2 we define our astrophysical parameters, the quasar model, and the simulation box parameters. In Section 3 we illustrate the workings of our models and summarize all our results. In Section 4 we perform various convergence tests for the damping wing signal across different box sizes, with the same set of parameters and grid resolution, to determine the optimum box length for our studies. Finally, in Section 5 we present our understanding of these results and the interesting features we observe. We assume a flat ΛCDM cosmology throughout our work, based on the results from Planck (N. Aghanim et al. 2020), with cosmological parameters h = 0.68, Ω_m = 0.3, Ω_b = 0.045, and σ₈ = 0.8.

In this section, we construct a suite of reionization simulations reflecting a wide range of astrophysical parameters to explore variations of the quasar damping wing signals at a fixed global neutral fraction. We describe the components of our models in detail.

2.1. 21cmFAST

21cmFAST is a seminumeric simulation tool that produces 3D cosmological boxes of various physical fields in the early Universe (A. Mesinger et al. 2011; J. Park et al. 2019; S. G. Murray et al. 2020; Y. Qin et al. 2020; J. B. Muñoz et al. 2022). For our case, we generated boxes with density and ionization fields subject to our range of astrophysical parameters. 21cmFAST creates evolved IGM density fields by first generating Gaussian cosmological initial conditions, which are then perturbed to later times using second-order Lagrangian perturbation theory (R. Scoccimarro 1998).

In 21cmFAST, the ionization field is constructed using a method motivated by excursion set theory. The core idea is to determine regions where the number of ionizing photons exceeds the number of neutral hydrogen atoms (S. R. Furlanetto et al. 2004). This is done using a filtering process that smooths the density field with a spherical top-hat filter (a caveat to this formalism is that the photons are not conserved; see T. R. Choudhury & A. Paranjape 2018). The ionization state of any given cell is then determined by comparing the number of photons, integrated over time and produced by UV rays with the number of neutral atoms within regions of decreasing filters of radius R. The size of these filters starts from the maximum photon horizon, R_mfp, down to the individual pixel resolution of a single cell, R_cell. The cells are flagged as fully ionized if the following condition is satisfied at any filter scale R (B. Greig & A. Mesinger 2018):

$$\begin{eqnarray}&&\zeta {f}_{{\rm{coll}}}(x,z,R,{{\rm{M}}}_{{\rm{\min }}})\geqslant 1.\end{eqnarray} \tag{ 1 }$$

Here, f_coll represents the fraction of collapsed matter residing within halos in the region R with mass greater than ${{\rm{M}}}_{{\rm{\min }}}$ that has collapsed to form bound structures, which is calculated using the Press–Schechter formalism (W. H. Press & P. Schechter 1974; J. R. Bond et al. 1991; C. Lacey & S. Cole 1993; R. K. Sheth & G. Tormen 1999).

The parameter ζ represents the UV ionizing efficiency of dark-matter halos within the filter size. It is the cumulative result of several factors: the halo mass, the number of stars within the halo, the number of ionizing photons they produce, the fraction of these photons that escape the halos, and the proportion of photons that successfully ionize surrounding gas. In the following subsections, we will discuss these terms in greater detail.

The final parameter, ${M}_{\min }$, defines the minimum mass of halos capable of supporting star formation. Below this threshold, star formation becomes inefficient, thereby limiting the production of ionizing photons by such galaxies. This suppression of star formation may result from supernova feedback, photoheating feedback, or inefficient cooling (M. L. Giroux et al. 1994; P. R. Shapiro et al. 1994; L. Hui & N. Y. Gnedin 1997; R. Barkana & A. Loeb 2001; V. Springel & L. Hernquist 2003; A. Mesinger & M. Dijkstra 2008; T. Okamoto et al. 2008; E. Sobacchi & A. Mesinger 2013a, 2013b).

Our 21cmFAST simulations estimate this suppression using a redshift independent duty cycle for a given halo mass M_h (J. Park et al. 2019):

$$\begin{eqnarray}&&{f}_{\,\rm{duty}\,}({M}_{h})=\exp \left(-\displaystyle \frac{{M}_{\min }}{{M}_{h}}\right).\end{eqnarray} \tag{ 2 }$$

Thus, for halos with a mass close to ${M}_{\min }$, only a fraction f_duty of them are forming stars with an efficiency of f_star, the rest do not. This parameter in conjunction with the global neutral fraction controls the ionized bubble size distribution within our simulation boxes.

2.2. Parameter Space

During the epoch of reionization, voxels in the simulation are ionized by UV photons that escape from dark-matter halos and ionize the neutral hydrogen in the IGM. A halo's efficiency in ionizing the surrounding IGM depends on the number of ionizing photons produced within the halo (proportional to the number of stars it contains) and the fraction of those photons that escape into the IGM. These factors are collectively represented by the parameter ζ, as described above.

The standard 21cmFAST approach considers a simplified model for calculating the stellar mass of galaxies and the number of ionizing photons in our simulations. The stellar mass of a galaxy, M_star, can be expressed as a function of its host halo mass, M_h as follows (M. Kuhlen & C.-A. Faucher-Giguère 2012; P. Dayal et al. 2014; P. S. Behroozi & J. Silk 2015; S. Mitra et al. 2015; S. J. Mutch et al. 2016; G. Sun & S. R. Furlanetto 2016; B. Yue et al. 2016; J. Park et al. 2019)

$$\begin{eqnarray}&&{M}_{\mathrm{star}}({M}_{h})={f}_{\mathrm{star}}\left(\displaystyle \frac{{{\rm{\Omega }}}_{b}}{{{\rm{\Omega }}}_{m}}\right){M}_{h}.\end{eqnarray} \tag{ 3 }$$

Here f_star denotes the fraction of galactic baryons present in stars. The f_star also has a power-law dependency on halo mass (P. S. Behroozi & J. Silk 2015; J. Park et al. 2019), which comes from the fact that gas is being exchanged between the IGM and the halo due to feedback processes:

$$\begin{eqnarray}&&{f}_{\mathrm{star}}={f}_{\mathrm{star},10}{\left(\displaystyle \frac{{M}_{h}}{1{0}^{10}{M}_{\odot }}\right)}^{{\alpha }_{\,\rm{star}\,}},\end{eqnarray} \tag{ 4 }$$

where f_star,10 is the normalization factor representing the fraction of galactic gas in stars normalized to a halo mass of 10¹⁰M_⊙, and α_star represents the power-law index of the dependence of f_star on halo mass.

We can express the other important quantity, f_esc, which governs the fraction of ionizing photons escaping from star-forming galaxies into the IGM, in a manner similar to f_star (J. Park et al. 2019),

$$\begin{eqnarray}&&{f}_{\mathrm{esc}}={f}_{\mathrm{esc},10}{\left(\displaystyle \frac{{M}_{h}}{1{0}^{10}{M}_{\odot }}\right)}^{{\alpha }_{\,\rm{esc}\,}},\end{eqnarray} \tag{ 5 }$$

where f_esc,10 is the normalization factor representing the fraction of ionizing UV escape fraction normalized to 10¹⁰M_⊙ halo mass fraction, and α_esc represents the power-law dependency of f_esc on halo mass. Both f_star and f_esc have a physical upper limit of ≤1 by definition. This parameter plays a critical role in understanding the effects of feedback processes. Given the uncertainty surrounding both the nature and the precise value of f_esc, it can be used to calibrate the value of the neutral fraction when varying other source parameters.

Together these quantities tell us how much stellar mass is present in any galactic dark-matter halo and how many ionizing photons they eject into the IGM. We can now calculate the ionizing UV efficiency of a galaxy (B. Greig & A. Mesinger 2017):

$$\begin{eqnarray}&&\zeta ({M}_{h})=30\left(\displaystyle \frac{{f}_{\mathrm{esc}}}{0.12}\right)\left(\displaystyle \frac{{f}_{\mathrm{star}}}{0.05}\right)\left(\displaystyle \frac{{N}_{\gamma /b}}{4000}\right)\left(\displaystyle \frac{1.5}{1+{n}_{\,\rm{rec}\,}}\right),\end{eqnarray} \tag{ 6 }$$

where N_γ/b is the number of ionizing photons per baryon produced in the stars and n_rec is the number of times hydrogen atom recombines. For our case, the rate of recombination is so small in the IGM compare with Hubble expansion that we can neglect n_rec in the last term.

In the following, we will abbreviate the normalization terms f_star,10 and f_esc,10 as f_star and f_esc, respectively, unless stated otherwise.

We aim to compare the effects of the aforementioned parameters on reionization topology at a fixed neutral fraction. In our models, we adjust f_esc to calibrate the desired global neutral fraction for a given combination of parameters. We expect f_star to be degenerate with f_esc, something that we also see in our results, although a subtle distinction remains due to differences in the halo mass corresponding to the upper limit value of 1. All of the parameters described thus far are astrophysical in nature and can be tuned in our simulations, thereby impacting the reionization topology on a global scale. Table 1 lists these parameters, along with their respective ranges and fiducial values used in this study.

Table 1. Astrophysical Parameters Space: This Table Lists All the Parameters Used in Our Study Along with Their Respective Ranges, Fiducial Values, and Parameter Types

Parameter	Range	Fiducial Model	Parameter Type
xH I (Global mean neutral fraction)	(0.25, 0.75)	0.5	IGM parameter
${M}_{\min }$ (Minimum halo mass to support star formation)	(2Mp, 200Mp)	20Mp	Source parameter
αesc	(−1,0)	−0.5	Source parameter
αstar	(0, 1)	0.5	Source parameter
${\mathrm{log}}_{10}{{\rm{f}}}_{{\rm{star}}}$	(−2,−0.25)	−1.125	Source parameter
tq (Quasar lifetime)	(0 Myr, 30 Myr)	1Myr	Quasar parameter
Mqso (Host halo mass)	(≈1 × 109M⊙, ≈4 × 1012M⊙)	≈4 × 1011M⊙	Quasar parameter

Download table as:  ASCIITypeset image

2.3. Quasar Model

After obtaining the reionization topology, we conduct our analysis by selecting random halos of a specific mass and observing a line of sight toward each halo, assuming it hosts a quasar. In contrast to the calculation of the ionization field described above, for this we require instantiating discrete halos within the seminumerical simulation box. We use the method from A. Mesinger & S. Furlanetto (2007) to locate dark-matter halos using excursion set theory on the initial conditions. These halo positions are then perturbed using velocity fields, calculated via linear perturbation theory, to obtain the corrected locations at the desired redshift. Although the ionized field we derive is based solely on the density field, we use this halo catalog to position our quasars. This approach enables us to include lower-mass halos in our damping wing analysis, which would otherwise require much higher particle resolution. Given that quasars are likely to reside in halos with masses around 10¹²M_⊙, we explored a range from 10⁹M_⊙ to 10¹²M_⊙ to fully capture the dependence on host halo mass (e.g., A.-C. Eilers et al. 2024, E. Pizzati et al. 2024) and probe the regime of damping wings that affect lower-mass galaxies as well (L. C. Keating et al. 2024).

As we are observing a line of sight toward a quasar, a sufficiently luminous quasar could influence the reionization topology in its vicinity, potentially changing the strength of its observed damping wing. To study this effect, we assume a quasar with a lifetime t_q, residing in a halo of mass M_qso. During its active lifetime, the quasar carves out an ionized bubble around it, R(t_q). The ionizing photon emission rate for luminous quasars at redshift z ≥ 7 is ${\dot{N}}_{ph}\simeq 1{0}^{57}\,{{\rm{s}}}^{-1}$ (D. J. Mortlock et al. 2011). The Strömgren radius of this bubble is expressed as (P. R. Shapiro & M. L. Giroux 1987; R. Cen & Z. Haiman 2000):

$$\begin{eqnarray}&&{R}_{({t}_{{\rm{q}}})}={\left(\displaystyle \frac{3{\dot{N}}_{ph}{t}_{{\rm{q}}}}{4\pi \langle {n}_{{\rm{H}}}\rangle }\right)}^{1/3},\end{eqnarray} \tag{ 7 }$$

where 〈n_H〉 is the average number density of neutral hydrogen within the sphere.

The above equation assumes a homogeneous reionization process, whereas reionization is inherently inhomogeneous. Moreover, quasar emission is likely anisotropic, meaning the equation does not fully capture the overall shape of the ionized bubble. However, as we are primarily concerned with ionizing photons along the line of sight, we can use the relation to estimate the expansion of the ionized bubble by considering the density and neutral fraction distribution along the line of sight, assuming spherical dilution of the photon flux. Specifically, we calculate the total number of neutral hydrogen atoms on a spherical shell of radius ${R}_{{t}_{{\rm{q}}}}$ and a thickness of 1 pixel. To ionize these hydrogen atoms, an equivalent number of photons is required. Therefore, we equate the total number of photons produced by the quasar over its lifetime, t_q, with the number of hydrogen atoms encountered on a spherical shell, incorporating variations in density and neutral fraction along our line of sight:

$$\begin{eqnarray}&&{\dot{N}}_{ph}{t}_{{\rm{q}}}={\int }_{0}^{{R}_{\mathrm{ion}}}{n}_{{\rm{H}}}(r)4\pi {r}^{2}dr.\end{eqnarray} \tag{ 8 }$$

The integration is performed over the voxels along the line of sight corresponding to the radius of the bubble, neglecting any recombination within this region, which should be negligible for the quasar lifetimes we consider. Here, n_H(r) represents the number density of hydrogen atoms in the voxel spanning the distance between r and r + dr. We thus obtain the radius ${R}_{{t}_{{\rm{q}}}}$. Quasar activity modifies the reionization topology locally, up to the extent of ${R}_{{t}_{{\rm{q}}}}$. We assume the ionizing photon emission rate, ${\dot{N}}_{ph}$, remains constant, meaning the quasar lifetime, t_q, is the sole parameter controlling the distance to the ionization front. Notably, the case of t_q = 0 is analogous to a sightline originating from a typical star-forming galaxy, which produces significantly fewer ionizing photons than a luminous quasar, especially for small host halo masses (M_h ≤ 10¹¹M_⊙). The effects of quasar lifetime are incorporated a posteriori into our calculation, specifically during the computation of random skewers through the simulation box. They are evaluated exclusively along the line of sight for computational efficiency.

2.4. Lyα Damping Wings

We calculate the transmitted flux from the total Lyα damping wing optical depth τ_D at an observed wavelength of [λ_obs = λ_α(1 + z)] (where λ_α = 1215.67 Å is the Lyα rest-frame wavelength) along the length of our random skewer (treated as the line of sight), originating from a halo at redshift z_s. This is done by summing the contribution of τ_D from each neutral hydrogen patch encountered along the length of the skewer using the approximation (A. Mesinger & S. R. Furlanetto 2008),

$$\begin{eqnarray}\begin{array}{rcl}{\tau }_{D}(z) & = & \displaystyle \frac{{\tau }_{\mathrm{GP}}{R}_{\alpha }}{\pi }\displaystyle \sum _{i}\left\{{x}_{\mathrm{HI}}(i)(1+\delta (i)){\left(\displaystyle \frac{1+{z}_{{b}_{i}}}{1+z}\right)}^{3/2}\right.\\ & & \left.\times \left[I\left(\displaystyle \frac{1+{z}_{{b}_{i}}}{1+z}\right)-I\left(\displaystyle \frac{1+{z}_{{b}_{i}}}{1+z}\right)\right]\right\},\end{array}\end{eqnarray} \tag{ 9 }$$

where ${\tau }_{{\rm{GP}}}\approx 7.16\times 1{0}^{5}{[(1+{z}_{s})/10]}^{3/2}$ is the Gunn–Peterson optical depth of the IGM (J. E. Gunn & B. A. Peterson 1965), and R_α = Λ/(4πν_α), where Λ = 6.25 × 10⁸ s⁻¹ is the decay constant of Lyα at resonance and ν_α = 2.47 × 10¹⁵ Hz is the frequency of the Lyα transition. x_{H I}(i) is the neutral hydrogen fraction in the ith patch, and δ(i) is the matter overdensity in that patch. Lastly, the integration term I is defined as (D. Mortlock 2016; T. Kist et al. 2025),

$$\begin{eqnarray}&&I(x)=\frac{{x}^{1/2}}{1-x}+2\mathrm{log}\left|\frac{1-{x}^{1/2}}{1+{x}^{1/2}}\right|.\end{eqnarray} \tag{ 10 }$$

In our calculation of optical depth, we do not include the absorption inside the proximity zone and are looking at the IGM damping wing alone, as our primary concern is to study the effects of reionization topology.

2.5. Summary of Parameters

2.6. Grid Parameters

In this section, we present our findings from running the grid of models. All simulations were conducted with boxes of 512 Mpc in length at redshift 7.0. Unless otherwise specified, the parameter values used in our case studies are based on the fiducial model 1.

For illustration purposes, the upper panel of Figure 1 displays a sample slice of the ionized box from our fiducial model. On this panel, the red (green and orange) dots represent halos, and the corresponding red (green and orange) lines indicate the trajectory of light from the quasar along the line of sight. In the lower panel of the same Figure, we show the resulting damping wing signal from the respective color-coded sightline.

Zoom In Zoom Out Reset image size

As expected, we observe the strongest damping of the transmitted flux for the orange halo, which is located in a predominantly neutral region and encounters several neutral islands along the sightline. In contrast, the green and red halos, situated in ionized regions, experience less damping. However, while both the red and green halos reside in ionized regions, the sightline from the red halo traverses more neutral islands compare with that of the green halo. As a result, the red sightline experiences relatively greater damping than the green one. This illustrates the dependence of the Lyα signal on both the local environment of the host halo and the global environment along the sightline.

In Figure 2 we present the slices of the ionized box for 0.25, 0.50, and 0.75 mean global neutral fraction to illustrate how the reionization topology changes with the x_{H I}. In Figure 3 we show the mean transmitted flux for the same values of x_{H I}. Unlike the median transmission flux, which represents a typical damping wing spectra, the mean transmission flux averages multiple damping wing profiles that terminate at various distances, effectively stacking them on top of one another. As a result of this stacking and the substantial cosmic variance, a long-tail behavior is observed in the mean signal. This tail is particularly prominent in cases with greater scatter, such as x_{H I} = 0.25.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Some studies have examined this stacking of profiles in detail (D. Ďurovčíková et al. 2024). However, for our purposes, we are more interested in examining the individual damping wing profiles and their variation with reionization topology. Therefore, we focus on the median signal, which better captures the appearance of a typical damping wing, rather than the mean. Furthermore, the median signal is more robust for high-redshift quasars, as the number of sightlines is limited. From this point onward, unless stated otherwise, our analysis will be based primarily on the median signal.

In Figure 4, we present the median damping wing signals for three different values of x_{H I} (0.25, 0.5, 0.75), originating from halos with fiducial M_qso. The shaded region represents the 68th percentile scatter of the damping wings around the median signal. As expected, for larger values of x_{H I}, the sightlines encounter more neutral regions, resulting in significantly more damping. The shaded region illustrates the variability of the damping wings within the 68th percentile around the median signal.

Zoom In Zoom Out Reset image size

We also observe that for higher values of x_{H I}, it is more probable for a sightline to intersect a neutral island earlier in its journey. Conversely, for lower values of x_{H I}, the sightline can travel a greater distance before encountering any neutral islands. This trend is reflected in the scatter width (ΔSW₆₈) of the damping wing profiles, shown in the lower panel of Figure 4. Notably, ΔSW₆₈ increases as the neutral fraction decreases, particularly at larger distances.

The most massive halos tend to reside deep within ionized regions, as these high-mass halos correspond to the peaks in the density field, which host higher concentrations of both high- and low-mass halos. In contrast, low-mass halos are more evenly distributed throughout the simulation box. This distribution is evident in Figure 5. The left panel shows that high-mass halos (≈9 × 10¹¹M_⊙) are biased toward ionized regions, while the right panel illustrates that low-mass halos (≈4 × 10⁹M_⊙) are dispersed throughout the box. This bias is also reflected in the damping wing signals for these halos. In Figure 6, we observe overall weaker damping for sightlines originating in massive halos, consistent with their tendency to be located in larger ionized regions.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 7 illustrates the effect of t_q on the damping wings. The first three plots show the location of the halo and the corresponding sightline, ordered by increasing t_q values [0, 1, 30] Myr. The red-shaded region represents the size of the ionized bubble. In the fourth plot, we compare the damping wing profiles for the quasar-off case (t_q = 0 Myr) with the quasar-on cases (t_q = 1 and 30 Myr). The t_q = 0 case is particularly relevant for damping wing signals originating from galaxies (e.g., L. C. Keating et al. 2024; H. Umeda et al. 2024).

Zoom In Zoom Out Reset image size

In Figure 8, we explore the effects of quasar lifetime for two extreme cases (t_q = 0 and t_q = 30 Myr) on an ensemble of halos with the fiducial M_qso. As expected, quasar lifetime significantly impacts the median damping wings, being the second most dominant factor after the neutral fraction. This raises the issue of degeneracy between quasar lifetime and neutral fraction when examining any particular damping wing profile (e.g., F. B. Davies et al. 2018). We note that if t_q is longer than 30 Myr, the quasar carves out a large enough ionized region that the damping wing signal is almost entirely erased, even in a fully neutral IGM. Therefore, for practical measurements, we restrict t_q ≤ 30 × 10⁶ yr. Despite the significant variation in the median damping wing signal across this range of t_q, the width of the scatter at fixed distance from the ionization front remains largely unchanged, indicating that t_q has a minimal effect on the scatter. We discuss this further in Section 5, where we evaluate the effective values of x_{H I} required to reproduce the median damping wing signals for different t_q. This behavior can be attributed to the fact that quasar activity primarily alters the reionization topology locally, while the scatter arises from the distribution of neutral hydrogen beyond the edge of the Local Bubble.

Zoom In Zoom Out Reset image size

We also explore the impact of varying ${M}_{\min }$, which dictates the minimum halo mass capable of supporting star formation. For higher values of ${M}_{\min }$, only the most massive halos contribute to the reionization topology. As these halos are far fewer in number compare with their lower-mass counterparts, the reionization topology is primarily governed by larger ionized bubbles. In Figure 9, we present slices of the ionized boxes for ${M}_{\min }=1{0}^{8.58}{M}_{\odot }$ and ${M}_{\min }=1{0}^{10.58}{M}_{\odot }$. As expected, the ionized field for ${M}_{\min }=1{0}^{8.58}{M}_{\odot }$ features finer ionized bubbles and extended neutral regions. In contrast, the slice corresponding to ${M}_{\min }=1{0}^{10.58}{M}_{\odot }$ exhibits larger, more coarsely distributed ionized regions. Consequently, for ${M}_{\min }=1{0}^{10.58}{M}_{\odot }$, we anticipate an overall reduction in the damping effect.

Zoom In Zoom Out Reset image size

In Figure 10, we observe that the change in reionization topology due to varying ${{\rm{M}}}_{{\rm{\min }}}$ has a noticeable impact on the damping wing profiles. Specifically, for ${M}_{\min }=1{0}^{10.58}{M}_{\odot }$, there is indeed less damping compare with ${M}_{\min }=1{0}^{8.58}{M}_{\odot }$, as the neutral islands are smaller. For the same reason, and as discussed in the case of x_{H I}, we observe more scatter in the damping wings for ${M}_{\min }=1{0}^{10.58}{M}_{\odot }$. Although the ${M}_{\min }=1{0}^{8.58}{M}_{\odot }$ box appears more patchy, the likelihood of encountering a neutral region along any random sightline is higher in this box compare with the ${M}_{\min }=1{0}^{10.58}{M}_{\odot }$ box. This is because the ${M}_{\min }=1{0}^{8.58}{M}_{\odot }$ box contains larger and more widely distributed neutral islands than the ${M}_{\min }\,=1{0}^{10.58}{M}_{\odot }$ box.

Zoom In Zoom Out Reset image size

Finally, in Figure 11 we present a comprehensive overview showcasing the effects of all the model parameters on the damping wings. We have already explored the impact of x_{H I}, t_q, and ${M}_{\min }$ in the previous sections. For the remaining parameters—α_star, α_esc, and f_star—we do not observe any significant effect on the median damping wings or the scatter. This lack of impact is likely due to these parameters being degenerate with another crucial factor, the escape fraction f_esc, which is used in the calculation of the ionizing efficiency of high-redshift galaxies, as expressed in the Equation (6).

Zoom In Zoom Out Reset image size

The ionizing efficiency ζ, which governs the global ionization state of the IGM, is directly influenced by f_esc. Because α_star, α_esc, and f_star primarily affect galaxy formation and star formation rates, their impact is indirectly reflected in f_esc. Thus, changes in these parameters are absorbed by adjustments in f_esc, which controls the total number of ionizing photons escaping into the IGM, leading to a muted effect on the observed damping wings.

In conclusion, while x_{H I}, t_q, and ${{\rm{M}}}_{{\rm{\min }}}$ play a direct and substantial role in shaping the damping wing profiles, parameters like α_star, α_esc, and f_star being degenerate with f_esc for the calculation of x_{H I} do not show any variations in our plots.

Throughout this work, we have emphasized the crucial role that large-scale structures play in the analysis. It is then natural to consider what size of the box is required for the strength and variance of the damping wing signal to converge. A small box may fail to capture high-mass halos and large structures essential for analyzing sightline-to-sightline scatter. On the other hand, larger boxes are more computationally expensive, inhibiting exploration of the full parameter space. To identify an optimal box size, we perform a convergence test of box size at fixed spatial resolution.

We begin by comparing the mean optical depth of damping wing profiles for halos at the highest and lowest common mass from each box. By analyzing the mean optical depth, we aim to understand the average reionization topology surrounding a given halo. In Figure 12, we observe that for low halo masses (M_qso ≈ 4 × 10⁹M_⊙), the optical depth damping wing profiles converge within a 10% error limit. This indicates that for the range of box sizes considered here, the topology and corresponding damping profiles for low-mass halos remain consistent and stable.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

However, in Figure 16 (presented in the Appendix) the convergence test for high-mass halos (M_qso ≈ 1 × 10¹²M_⊙) shows an unexpected trend. We observe that the mean optical depth damping profiles increase as we move from our fiducial model with box length (512 Mpc)³ up to (768 Mpc)³. Beyond this size, however, the profiles start to decrease.

To ensure that this behavior is not a result of selecting a nonrepresentative halo (i.e., not the most massive), but instead stems from selecting the common most massive halo across boxes, we replicate the analysis using the 100 most massive halos within each box. This procedure mimics that used by past works on quasar damping wings (e.g., F. B. Davies et al. 2018). Additionally, to account for uncertainties arising from the randomness of the initial conditions, we generate six different initial condition boxes for each box size and average the results. This analysis is presented in Figure 17, where we observe a similar trend. The consistent behavior across different initial conditions and multiple halos confirms that this phenomenon is not due to the choice of halos or initial conditions.

We originally expected that, as box size increases, the average damping effect decreases. This is because high-mass halos are typically located in ionized regions, and larger boxes, containing larger halos, would naturally have bigger ionized bubbles surrounding them. Consequently, damping should continue to decrease, and the deviation from the base model should increase monotonically. However, our simulation results consistently indicate the opposite once the box size exceeds (768 Mpc)³.

One possible explanation for this unexpected behavior could be related to the halo mass function for high-mass halos, which decreases exponentially at the highest masses. The increase in the total number of particles (and thus the volume of the box) may not be enough to generate a sufficient number of high-mass halos. As a result, the reionization topology around this exponential tail of massive halos in larger boxes might still be dominated by a mix of smaller and larger ionized bubbles, rather than exclusively by large-scale bubbles surrounding massive halos.

In all cases, we maintained a constant mean global neutral fraction (x_{H I} = 0.5), which should ensure consistency in the overall ionization state of the simulation. To further test this hypothesis, we eliminate smaller bubbles from forming by setting the minimum halo mass ${M}_{\min }=1{0}^{11}{M}_{\odot }$. This will exclude low-mass halos, effectively removing the contribution of smaller ionized bubbles from the topology.

In Figure 18, we observe that when we eliminate the smaller bubbles, the damping wing profiles converge monotonically from (768 Mpc)³ to (896 Mpc)³ as expected. This finding supports the idea that the issue of convergence likely arises from the interplay between small and large ionized bubbles. When both types of bubbles are present, as in the case of the larger (896 Mpc)³ box, this interplay reduces the effective size of the ionized bubbles surrounding high-mass halos, leading to a decrease in the damping signal transmission instead of the anticipated increase.

To further validate our hypothesis, we can take two approaches: (1) artificially increase the size of the ionized bubbles by turning on the quasar (t_q > 0) or (2) calculate the bubble size distribution across different mass bins and box sizes.

In the first case, we repeat the convergence test for the 100 most massive halos, but this time we set t_q = 1 Myr, effectively increasing the ionized region surrounding each halo. By activating the quasar, we artificially boost the size of the bubbles, thereby mitigating the effect of small ionized bubbles. In Figure 19, we show the results of this convergence test. As expected, the damping wing profiles begin to converge within the 10% error limit, which supports the idea that the interplay between small and large bubbles was indeed responsible for the earlier discrepancies in convergence. This indicates that the effective size of the ionized bubbles from high-mass halos had been reduced when small bubbles were present in the simulation, and this effect was mitigated when the quasar was turned on. Thus, we confirm that for larger boxes, the choice of box size becomes crucial. Care must be taken to ensure that the size of the simulation box adequately captures both large-scale structures and high-mass halos. For all of our models, the fiducial box already assumes t_q = 1 Myr, hence the box length of 512 Mpc provides an optimum compromise between large and fast box volumes.

In the latter approach, we calculate the bubble size distribution. As before, we use six different realizations of the initial conditions for each box size, referring to each as a batch for the respective box size. Within each batch, we select 100 halos from each mass bin. For each halo, we generate 10 sightlines in random directions, then measure the neutral fraction x_{H I} along the length of each sightline. If the number of halos in a mass bin is less than 100, we compensate by generating additional random sightlines to ensure that the total number of halo–sightline combinations remains consistent across all mass bins. For each sightline, we clip the length as soon as we encounter a neutral voxel (x_{H I} = 1.0), representing the boundary of the ionized region. We then average over all these sightlines for each halo to calculate the spherical average radius of the ionized bubble around the halo. By repeating this process for all halos in the mass bin, we obtain a distribution of bubble sizes for halos of a particular mass. We repeat this procedure for all mass bins across all boxes in each batch, which results in the curves shown in Figure 13. As we have seen, the average bubble size increases with the halo mass while the scatter decreases Figure 6.

We then select the mass bins corresponding to the 100 most massive halos from each simulation box in a batch, concatenate them, and repeat the above calculation to obtain the distribution of ionized region sizes. In Figure 14, we plot the aforementioned distribution of ionized radii/sizes. From the inset panel of the same Figure, we observe that the mean of the distribution for the (768 Mpc)³ box is higher than that of the (896 Mpc)³ box, thus confirming our hypothesis.

Finally, to conclude, for our fiducial model, we assume t_q = 1 Myr. Hence, for all our analyses, the 512 Mpc box converges within the 10% error limit and is sufficient for our study.

The results above have established an understanding of how each parameter affects the median damping wing profile and the sightline-to-sightline scatter of the damping wings around the median. We observed that parameters influencing the global ionization topology, such as x_{H I} and ${M}_{\min }$, significantly impact the scatter of the damping wings. In contrast, t_q, which alters the local ionization topology, does not substantially affect the overall scatter. In this section, we aim to further discuss this behavior and explore the possibility of handling the degeneracy between x_{H I} and t_q.

As noted in Section 4, the damping wing signals vary rapidly with changes in both x_{H I} and t_q, raising the issue of degeneracy in the median signal. In Figure 4, we observe that variations in the global mean neutral fraction (x_{H I}) have a pronounced effect on both the median and the scatter. This is expected, as any change in x_{H I} necessitates an alteration in the overall reionization topology of the box. However, if the quasar is switched on for a time t_q, it only modifies the neutral fraction within the region described by Equation (7), which is a local phenomenon. Therefore, this local modification contributes minimally to the scatter of longer damping wings, which extend well beyond the radius of the ionized bubble.

We can potentially use this difference in impact to explore the degeneracy between x_{H I} and t_q. As illustrated in Figure 15, when we match the effective value of x_{H I} to obtain the same median damping wing for different values of t_q, we are unable to recover the scatter. Furthermore, we observe that the scatter induced by x_{H I} evolves much more rapidly than that caused by t_q. Thus, for an ensemble of damping wings, studying the scatter in conjunction with the median signal may provide a more robust constraint on both of these parameters.

We also observe that ${M}_{\min }$ has a noticeable effect on the median damping wing signals and a stronger effect on the scatter. As shown in Figure 9, the strong variation in scatter arises from the fact that changes in ${M}_{\min }$ affect the bubble size distribution around the halos. With larger and coarser bubbles, the probability of encountering a neutral island along a random sightline varies significantly with direction. A skewer can travel much longer distances before encountering any neutral patch, thus leading to a strong correlation with the scatter. This suggests that the statistics of damping wing absorption features may provide a novel probe of the ionizing output of the lowest mass galaxies, which can be regulated both by internal (supernovae, winds) and external (photoionization) feedback processes, as well as the escape fraction of ionizing photons.

Finally, in Section 4, we establish that the damping wing profiles computed using a simulation box length of 512 Mpc converges within the 10% error limit with our fiducial parameters, and hence is sufficient for our studies.

With our study, we have shown that the Lyα damping wings are a very powerful tool to not just constrain the neutral fraction (x_{H I}), but also the properties of the source model. We have established that other than x_{H I}, the quasar lifetime (t_q), the mass of the host halo (M_qso), and the minimum mass of the halo that can support star formation (${M}_{\min }$) can have significant effect on the damping wing profile. We also looked at the sightline-to-sightline scatter of damping wings and showed how this scatter can be used to somewhat break the degeneracy between t_q and x_{H I}. This is because quasar activity primarily alters the reionization topology locally, while the scatter arises from the distribution of neutral hydrogen beyond the edge of the Local Bubble. In the near future we are expected to see more than just a few z > 7 quasars. The Euclid mission is poised to uncover hundreds of luminous high-redshift quasars (R. Barnett et al. 2019; J.-T. Schindler et al. 2023). Thus, making this systematic uncertainty not only relevant but possible to measure.

We thank Joseph Lewis for useful discussions on the interpretation of ${M}_{{\rm{\min }}}$. S.E.I.B. is supported by the Deutsche Forschungsgemeinschaft (DFG) through Emmy Noether grant No. BO 5771/1-1.

Figures 16–19 show convergence test results for the 100 most massive halos from each box.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size