The Characteristic Shape of Damping Wings During Reionization

Spectroscopic analysis of Ly$\alpha$ damping wings of bright sources at $z>6$ is a promising way to measure the reionization history of the universe. However, the theoretical interpretation of the damping wings is challenging due to the inhomogeneous nature of the reionization process and the proximity effect of bright sources. In this Letter, we analyze the damping wings arising from the neutral patches in the radiative transfer cosmological simulation suite Cosmic Reionization on Computers (CROC). We find that the damping wing profile remains a tight function of volume-weighted neutral fraction $\left_{\rm v}$, especially when $\left_{\rm v}>0.5$, despite the patchy nature of reionization and the proximity effect. This small scatter indicates that with a well-measured damping wing profile, we could constrain the volume-weighted neutral fraction as precise as $\Delta \left_{\rm v} \lesssim 0.1$ in the first half of reionization.


INTRODUCTION
The epoch of reionization (EoR) brought about a major change to the global properties of the intergalactic medium (IGM) within the first billion years of the universe.Thanks to the numerous data obtained by JWST, we are now in an excellent position to understand this frontier in astrophysics.
One of the fundamental questions surrounding the EoR concerns the timing and duration of reionization, which is not yet wellconstrained.Several methods have been employed to measure reionization, each with its own unique strengths and limitations.One of the earliest constraints came from cosmic microwave background experiments that utilized the Thompson scattering effects of free electrons.For example, Planck Collaboration et al. (2020) has constrained the midpoint of reionization redshift to be  mid = 7.82±0.71(see also Zahn et al. 2012;Dunkley et al. 2011).However, accurately measuring the entire history of reionization from start to end using Thompson scattering effect on CMB alone is challenging due to its integrated nature.Another powerful method to characterize the full cosmic dawn and reionization history is through the 21cm line emission from neutral hydrogen.However, at such long wavelengths, the foreground is orders of magnitude brighter than the signals, making data reduction notoriously difficult (e.g., Kern & Liu 2021;Mertens et al. 2020).
Another alternative of measuring the entire reionization history is to use the Ly absorption in front of bright sources at different redshift during the EoR (e.g., Miralda-Escudé 1998;Totani et al. 2016;Greig et al. 2017;Davies et al. 2018;Greig et al. 2019Greig et al. , 2022;;Gnedin & Prada 2004;Mason et al. 2018;Mason & Gronke 2020).As a strong resonant line, Ly is sensitive to any trace of neutral hydrogen, allowing us to detect the very end of reionization (neutral fraction 10 −4 ) (e.g., Bosman et al. 2018;Zhu et al. 2021;Bosman et al. 2022).Moreover, when there are neutral patches left in the IGM, the Ly absorption line displace a large damping wing, reaching ★ E-mail: hqchen@cita.utoronto.cathousands of km/s in the spectrum where the flux is suppressed.Miralda-Escudé (1998) shows that assuming a uniform reionization model, the damping wings have characteristic shape and can be used in constrain neutral fraction for bright background sources like gamma-ray bursts (GRBs).
However, reionization is a patchy process.The neutral fraction does not drop uniformly everywhere.Rather, some regions become highly ionized first, while other regions, shielded from ionizing sources, remain neutral until much later.In fact, many seminumerical codes based on excursion-set formalism (e.g., Furlanetto et al. 2004;Mesinger & Furlanetto 2007) treat every point in the universe as either neutral or ionized, therefore the term neutral fraction is meaningful only when averaged over certain volume or mass.In the literature of reionization, the term neutral fraction most commonly refers to the volume-weighted neutral fraction over the entire universe  HI v (e.g., Zhu et al. 2021).
Given the patchy nature of reionization, one natural question is whether the variance of the damping wing profile is too large to differentiate universes with different  HI v , or if it is small enough that the characteristic shape still holds.Another complication arises from the fact that bright sources, such as quasars, which provide high-resolution spectra for analysis, emit a large amount of ionizing radiation themselves.This radiation can alter the local morphology of reionization.Many semi-numerical methods can create a map of ionized bubbles created from typical galaxies (e.g., Furlanetto et al. 2004;Mesinger & Furlanetto 2007;Mesinger et al. 2011;Zahn et al. 2011), but unusual bright sources like quasars are not modeled.Does the removal of neutral patches close to bright sources like quasars significantly change the shape of the damping wing?
In this Letter, we use a radiative transfer cosmological simulation suite Cosmic Reionization on Computers (CROC) to address the above questions.Such a study is timely as more and more bright sources at  > 6 are spectroscopically followed-up and available to be used in constraining reionization history.The letter does not intend to describe the full process of extracting neutral fraction from data, but serves to estimate the optimal precision of neutral fraction measurement achievable using Ly damping wing.

SIMULATION
We use CROC simulations1 (Gnedin 2014) to study the damping wings arising from patchily ionized IGM.The CROC project uses the Adaptive Refinement Tree (ART) code (Kravtsov 1999;Kravtsov et al. 2002;Rudd et al. 2008) to reach high spatial resolution (base grid length = 39 ℎ −1 ckpc, peak resolution ∼ 100 pc in physical units).CROC simulations include relevant physics such as gas cooling, heating, star formation, stellar feedback and on-the-fly radiative transfer (Gnedin & Abel 2001).The main ionization sources in the simulations are star particles which are formed in dense gas in galaxies.
In this project, we primarily use the uniform-grid data in one of the 40 ℎ −1 cMpc runs (CROC B40F) alongside with Rockstar (Behroozi et al. 2013) halo catalogs to locate dark matter halos.The uniform-grid data contain gas properties of neutral fraction, density, temperature in each base grid cell.They are saved frequently (with increments in expansion factor Δ = 0.001) so that we can sample a large range of  HI v and study the entire reionization process.In Figure 1, we show the neutral fraction map at three different redshifts overlaid with halos of different masses.

RESULTS
To simulate the damping wing profiles, we first draw skewers (sightlines) starting from massive halos.We locate halos from Rockstar halo catalogues in the uniform-grid box.At each redshift, we select the 100 most massive halos and draw 10 skewers of length 200 cMpc/ℎ uniformly distributed in a 3D sphere.In the left panel of Figure 2, the brown line shows the neutral fraction along one example skewer drawn from a snapshot where the volume-weighted neutral fraction is  HI v = 0.5.To study the universe with different neutral fractions  HI v , we use skewers drawn at different redshifts of the same simulation run.When calculating Ly absorption, we keep the neutral fraction and temperature of each cell unchanged while scale the physical length and density to a certain redshift   by  and  −3 , where  is the expansion factor, respectively.The results shown in this paper are calculated for   = 6.54.
An unusually bright source like a quasar could push the I-front farther away.To mimic such an extra ionizing effect, when calculating the damping wing, we first draw a random number from a uniform distribution between [0, 40] cMpc/ℎ and remove all neutral gas within this distance.This procedure aims to examine the maximum variance of the damping wing shape.Then we convolve the Voigt profiles from the rest of the neutral cells ( HI > 0.5) along the skewer.In the left panel of Figure 2, we show this procedure in a skewer drawn from the box with  HI v = 0.5: the faint blue, orange and green vertical lines show three random positions within which we remove all neutral gas, and the solid profiles are the damping wing arising from the remaining neutral gas, integrated until 200 cMpc/ℎ.We find that despite the length of the first neutral patch are different, after convolution with all neutral patches behind it, the shapes of the profiles are very similar.This is more evident in the right panel, where we compare these profiles after aligning them at the starting position (the first point where transmission drops to zero).
In Figure 3, we plot the median of the aligned damping wing profiles in snapshots of different  HI v using solid lines, with each colored band showing the 68% scatter.For  HI v >= 0.5, the scatter of the wing profile is very small despite the patchy nature of reionization, and profiles with Δ  HI v = 0.25 are clearly separated.We also compare the damping wing profiles with the ones that created without randomly cutting the inner region (dash-dotted lines).If the inner neutral regions are not excised, the median damping wing is slightly stronger, but well within the scatter.The scatter of the no-cut case is almost identical to the previous case and thus not shown.We also calculate the damping wings assuming a uniform density and uniform reionization scenario (every cell has the same neutral fraction  HI =  HI v and every cell contributes to the damping wing), and the results are shown as dotted lines.Compared with the patchy ionization scenario with inner region excised, the damping wing is in general stronger, especially for 0.5  HI v 0.75, but the differences are still small compared to the scatter.

Cosmic variance
Due to the small size (40 ℎ −1 cMpc) of the simulation box, one might question whether the small scatter shown in the last section still holds when considering cosmic variance.To investigate this, we repeat the procedure in another box (CROC B40C).Both simulations have the same physics but different initial conditions ("DC modes").As a result, B40F reionized the latest (reionization midpoint  mid = 7.4) while B40C reionized earliest ( mid = 8.2) in all the six 40 ℎ −1 cMpc CROC realizations.Therefore, the density environments and halo distribution in these two box should differ maximally in all six realizations, and comparing damping wings in these two boxes helps us understand the stochasticity due to cosmic variance.In Figure 4, we compare the damping wings in box B40C with B40F of the previous section.We find that the mean and the scatter are almost identical, suggesting that the damping wings indeed have a characteristic shape as a function of  HI v .

Practical use
Our simulations show that for a mostly neutral universe (  HI v > 0.5), the scatter in damping wing profiles is small enough to distinguish between Δ  HI v ≈ 0.1.However, measuring the entire damping wing profile is complicated in practice.In this subsection we briefly discuss the prospects of using damping wing to constrain  HI v .Miralda-Escudé (1998) originally proposes GRB afterglows as the best candidates for measuring  HI v with damping wings.Compared with galaxies or quasars, GRBs have many advantages.They are thought to be produced in normal galaxies and thus live in less biased environments (Woosley 1993).The number of integrated ionizing photons they contribute is also very small and unlikely to enlarge the local ionized bubble.In addition, they are intrinsically bright to be spectroscopically followed-up.One challenge of using GRB afterglows is how to model the Damped Ly absorbers (DLAs) in the host galaxies.Lidz et al. (2021) shows that using the empirical distribution from current GRB afterglow spectra, one could model the local DLA distribution and marginalize this nuance parameter.Although Lidz et al. (2021) does not consider the scatter of damping wing profiles, the small scatter we find in CROC simulations supports their forecast that with 20 GRB afterglows with spectra resolution  3000 and signal-to-noise ratio (SNR) 20, one could reach a precision close to 15% in the first half of reionization.
Quasars are the kind of sources we can obtain the highest resolution spectra at  > 6.The current highest quality sample of  > 6 quasar spectra have SNR 50 and  10000 (D'Odorico et al. 2023).Thanks to their high luminosity, the residual neutral fraction in their proximity zone is small enough to allow significant flux on the blue side of Ly line.Such flux offers extra information about the shape of the damping wing.The challenge of using quasars is that by  ≈ 7, a quasar may have enlarged the local bubble significantly.Due to the decrease in quasar radiation with distance, the transmitted Ly flux also decreases.This reduction in flux compromises the constraining power on the starting point of the damping wing, which is crucial for anchoring the shape of the damping wing.In the ideal case, we may catch a quasar in its bright phase, where the integrated ionizing photons emitted by the quasar is still small while the instantaneous luminosity is high enough to create a highly transparent proximity zone.This would allow us to observe details of the Ly forest and measure flux close to the starting point of the damping wing, providing greater constraining power for the shape of the entire damping wing.In addition, similar to the GRB afterglow case, we need to develop a better understanding of how to model the intrinsic quasar continua.Since the scatter in damping wing profile is < 10% at wavelength < −1000 km/s from the starting point of the damping wing, it is ideal to have an accuracy in continuum recovery better than 10% across the quasar Ly emission line (from −4000 km/s to where no transmitted flux presents).
With the successful operation of JWST, we now have the capability to measure spectra from Lyman Break Galaxies (LBGs) or Lymanalpha Emitters (LAEs) (e.g., Bolan et al. 2022).While these sources are more numerous than quasars, their low luminosity limits the achievable spectral resolution.As a result, information about the Ly damping wing is mainly contained in the equivalent width (EW) measurements.However, if we can combine the information of both the LBG/LAE positions and their EWs, it would be promising to constrain the neutral fraction by considering both the damping wing strength and the size of ionized bubbles (e.g., Mason & Gronke 2020;Tang et al. 2023).This avenue will be explored in future work.

CONCLUSIONS
In this paper, we analyze the damping wings arisen from the partially ionized IGM in a self-consistent radiative transfer cosmological simulation suite CROC.We find that when the volume-weighted neutral fraction  HI > 0.5, the shape of the damping wing has a characteristic shape with small scatter ( 10%).This scatter remains small even after an unusually bright source (such as a quasar) erodes a significant amount of neutral gas around it.This is because the damping wing arises from the collective, convoluted Voigt profiles along a large distance (hundreds of comoving Mpc).We also calculate the damping wing profiles in a uniform reionization case, and we find that it lies within the 68% scatter.The small scatter in the damping wing profiles indicates that we can expect an accuracy of Δ  HI v ≈ 0.1 if we could measure the damping wing profile precisely.In reality, there are several complications, notably how to model the intrinsic source spectra and the Ly absorption within the ionized bubble.The profiles we find suggest that in order to achieve the best constraints in neutral fraction, we should aim for an accuracy of continuum fitting better than 10% across the Ly emission line of the source (from −4000 km/s to where no transmitted flux presents).For a very bright source such as a quasar, the complication of Ly absorption inside the ionized bubble could potentially be mitigated by properly modeling the large-scale structure, which we plan to explore in the future.

Figure 1 .Figure 2 .
Figure 1.Neutral fraction map in the simulation (CROC B40F), overlaid with halos of different masses.Red dots represent top 1000 massive halos and black dots are all halos with mass larger than 10 9 M .Each panel show a slice of thickness of 2 cMpc.The colorbar shows the neutral fraction of each cell in logarithmic scale.

Figure 3 .Figure 4 .
Figure 3.Damping wing profiles at different global volume-weighted neutral fraction  HI v .Lines show the median while bands show 68% scatter.The solid lines show profiles created with a randomly excised inner region as described in Sec. 3. The dash-dotted lines show profiles which we do not excise a random inner region and the dotted lines show the cases of uniform reionization.