JWST observations of galaxy damping wings during reionization interpreted with cosmological simulations

Spectra of the highest redshift galaxies taken with JWST are now allowing us to see into the heart of the reionization epoch. Many of these observed galaxies exhibit strong damping wing absorption redward of their Lyman-$\alpha$ emission. These observations have been used to measure the redshift evolution of the neutral fraction of the intergalactic medium and sizes of ionized bubbles. However, these estimates have been made using a simple analytic model for the intergalactic damping wing. We explore the recent observations with models of inhomogeneous reionization from the Sherwood-Relics simulation suite. We carry out a comparison between the damping wings calculated from the simulations and from the analytic model. We find that although the agreement is good on the red side of the Lyman-$\alpha$ emission, there is a discrepancy on the blue side due to residual neutral hydrogen present in the simulations, which saturates the intergalactic absorption. For this reason, we find that it is difficult to reproduce the claimed observations of large bubble sizes at z ~ 7, which are driven by a detection of transmitted flux blueward of the Lyman-$\alpha$ emission. We suggest instead that the observations can be explained by a model with smaller ionized bubbles and larger intrinsic Lyman-$\alpha$ emission from the host galaxy.


INTRODUCTION
Absorption lines observed in the spectra of luminous objects, blueward of the Lyman- (Ly) wavelength at the redshift of the source, indicate the presence of intergalactic neutral hydrogen.At redshifts approaching the epoch of reionization, this absorption is observed to saturate and regions with no significant transmission are detected on scales of tens of comoving Mpc (Becker et al. 2015;Zhu et al. 2021).However, it is difficult to use these observations to infer the presence of completely neutral gas, as for gas at the mean cosmic density, complete saturation will occur already at a volume-weighted average neutral hydrogen fraction ⟨ HI ⟩ v ∼ 10 −5 (McQuinn 2016).
More compelling evidence for a significantly neutral intergalactic medium (IGM) is an observation of Ly absorption that extends to wavelengths redward of the source redshift (Miralda-Escudé 1998).When the diffuse IGM has a high volume filling factor of neutral gas, its optical depth becomes large enough that the absorption profile is dominated by the natural line broadening described by the Lorentzian component of the Voigt profile (Meiksin 2009).This means that there is a reasonable probability for photons to scatter in the wings of the line, even several thousands of kilometres away from line centre.
★ laura.keating@ed.ac.ukThis results in an absorption profile known as a damping wing.The strength of the damping wing is proportional to the optical depth of the diffuse IGM and hence the volume-averaged neutral fraction (Gunn & Peterson 1965;Chen 2023), allowing for constraints on the reionization history from individual bright objects.
Multiple quasars observed at  = 7 and above display the signature of IGM damping wings (Mortlock et al. 2011;Bañados et al. 2018;Wang et al. 2020;Yang et al. 2020).By comparing the observed damping wings with theoretical models that also take local ionization by the quasar into account, the neutral fraction of the average IGM at the source redshift can be inferred (Bolton et al. 2011;Greig et al. 2017Greig et al. , 2019;;Davies et al. 2018a).A challenge in these measurements is the determination of the intrinsic Ly emission of the quasar, with different reconstructions leading to different constraints on the neutral fraction (e.g., Greig et al. 2022).Despite this, there is consensus from the different analyses of the observed damping wings that reionization is incomplete above  = 7 (Fan et al. 2022).The constraints are further limited by the redshifts of the most distant quasars currently known, although this is expected to improve with the Euclid wide survey, with more than 60  > 8 quasars with magnitude  < 24 predicted to be discovered (Schindler et al. 2023).
The presence of damping wings in the spectra of long-duration gamma-ray burst (GRB) afterglows can also be used to constrain the timing of reionization (Totani et al. 2006;Chornock et al. 2013;Har-toog et al. 2015).This has some advantages over studies with quasars (Barkana & Loeb 2004).The spectra of GRB afterglows close to the Ly wavelength can be described as a power law, sidestepping the issue of reconstructing an intrinsic emission line profile.The hosts of GRBs are further expected to be star-forming galaxies, which are more numerous than quasars and which likely live in less biased environments of the Universe.The main disadvantage is that these GRB afterglows frequently show evidence for damped Ly systems (DLAs) at the redshift of the host galaxy (e.g., Jensen et al. 2001), which will also produce a damping wing that may eclipse the intergalactic signal.There are also relatively few high signal-to-noise GRB afterglow spectra currently available, but detections by future missions are predicted to provide competitive constraints on the evolution of the IGM neutral fraction (Lidz et al. 2021;Tanvir et al. 2021).
The impact of IGM damping wings can also be detected statistically from populations of galaxies.The luminosity function of Ly emitting galaxies is observed to decline above  = 6 more rapidly than expected from comparison with UV luminosity functions (Konno et al. 2014(Konno et al. , 2018)), indicating that absorption from the IGM may be obscuring the Ly emission of the galaxies.This is complemented by the declining fraction of Lyman-break galaxies that show Ly emission above  = 6 (Stark et al. 2010;Pentericci et al. 2011Pentericci et al. , 2014;;Schenker et al. 2012Schenker et al. , 2014)).Through comparison with theoretical models, these observations can be translated into constraints on the progress of reionization (Mason et al. 2018(Mason et al. , 2019;;Hoag et al. 2019;Jung et al. 2020;Bolan et al. 2022;Wold et al. 2022;Jones et al. 2023).Detections of Ly emission in individual galaxies can also be used to estimate the sizes of ionized bubbles surrounding the host galaxies (Mason & Gronke 2020;Whitler et al. 2023;Witstok et al. 2023).
JWST is now opening up a new window in the study of IGM damping wings in galaxies, with damping wings being routinely observed in the spectra of individual galaxies (Curtis-Lake et al. 2023).By comparing the strength of these damping wings with the analytic model of Miralda-Escudé (1998), the spectra can be used to place constraints on the IGM neutral fraction (Curtis-Lake et al. 2023;Hsiao et al. 2023).Umeda et al. (2023) used the same model to fit the strength of damping wings in stacks of galaxy spectra in different redshift bins, placing constraints on the evolution of the volume-weighted average neutral fraction and ionized bubble size out to  = 12.They found evidence for an IGM neutral fraction that is increasing with redshift, and that is consistent with empirical models of galaxy-driven reionization (Ishigaki et al. 2018;Finkelstein et al. 2019;Naidu et al. 2020).Interestingly, the bubble sizes that Umeda et al. (2023) measured were towards the larger end of what is predicted by theoretical models (Lu et al. 2023).However, as with the GRB afterglows, the interpretation of these results may be complicated by DLAs associated with the host galaxies.Indeed, Heintz et al. (2023) recently found evidence for three  > 8.8 galaxies hosting DLAs with column densities  HI > 10 22 cm −2 .
Interpretation of these observations requires a model that captures all of the relevant physics, such as residual neutral gas in ionized regions of the IGM (Mesinger & Haiman 2004;Laursen et al. 2011;Bolton & Haehnelt 2013;Mesinger et al. 2015;Mason & Gronke 2020), the impact of infalling gas onto the host halo (Dĳkstra et al. 2007;Sadoun et al. 2017;Weinberger et al. 2018;Park et al. 2021) and the inhomogeneity of reionization (Mesinger & Furlanetto 2008;McQuinn et al. 2008;Garel et al. 2021;Gronke et al. 2021;Smith et al. 2022;Qin et al. 2022;Chen 2023).In this paper, we explore recent JWST observations of IGM damping wings in the context of these issues, analysing damping wings constructed from lines of sight through simulations of inhomogeneous reionization from the Sherwood-Relics simulation suite (Puchwein et al. 2023).Specifically, we critically evaluate the use of the Miralda-Escudé (1998) model to constrain reionization with recent JWST observations.We describe our models of the IGM damping wing in Section 2. In Section 3, we compare the results of the damping wings constructed from the reionization simulations to the commonly used Miralda-Escudé (1998) model.Section 4 presents a comparison of our results with damping wings observed in JWST galaxies, focusing specifically on the recent results of Umeda et al. (2023) andHeintz et al. (2023).Finally, in Section 5, we summarise our results and present our conclusions.Throughout the paper, we assume the Planck Collaboration et al. ( 2014) cosmological parameters, namely Ω m = 0.308, Ω Λ = 0.692 and ℎ = 0.678.

The Sherwood-Relics simulations
To generate mock IGM damping wings, we analyse simulations from the Sherwood-Relics simulation suite1 (Puchwein et al. 2023).These are a set of cosmological hydrodynamical simulations that build upon the original Sherwood suite (Bolton et al. 2017) and which were designed to model the evolution of the IGM during and after the epoch of reionization.The simulations were performed with a modified version of the smoothed particle hydrodynamical code pgadget-3 (Springel 2005).The fiducial simulation we analyse has a box size of 40 ℎ −1 cMpc and 2 × 2048 3 gas and dark matter particles, resulting in particle masses of  gas = 9.97 × 10 4 ℎ −1  ⊙ and  dm = 5.37 × 10 5 ℎ −1  ⊙ .The gravitational softening used was  soft = 0.78 ℎ −1 kpc.We further analyse a larger, lower resolution volume with box size 160 ℎ −1 cMpc and 2 × 2048 3 gas and dark matter particles.This simulation has  gas = 6.38 × 10 6 ℎ −1  ⊙ ,  dm = 3.44 × 10 7 ℎ −1  ⊙ and  soft = 3.13 ℎ −1 kpc.Star formation in these simulations is treated using a computationally efficient subgrid model, where any gas with temperature  < 10 5 K and density 1000 times the cosmic mean density is immediately converted into star particles.Although this approach is highly simplified, it should not affect the properties of the low-density IGM in which we are interested (Viel et al. 2004).
To account for the patchy nature of reionization, we utilise a subset of the Sherwood-Relics simulations that were performed with an inhomogeneous UV background.Although it is now becoming feasible to perform fully-coupled radiation-hydrodynamic galaxy formation simulations of cosmic reionization (see Gnedin & Madau 2022, for a recent review), a different approach has been taken in the Sherwood-Relics project.We give only a brief overview of the method here, but note that it is described in greater detail in Puchwein et al. (2023).We first perform a cosmological simulation with p-gadget-3 using the uniform UV background of Puchwein et al. (2019), saving snapshots of this simulation every 40 Myr.We then map these snapshots onto 2048 3 Cartesian grids and post-process them with the GPU-powered radiative transfer code aton (Aubert & Teyssier 2008, 2010).Multiple iterations of the radiative transfer simulations are performed, varying the redshift evolution of the volume emissivity assumed in the simulation until the desired reionization history is reached.We use a simple model for the sources, such that their emissivity is proportional to their halo mass (Iliev et al. 2006;Chardin et al. 2015).For each of the 2048 3 cells on the grid, the reionization redshift and  (Mason et al. 2018(Mason et al. , 2019;;Bolan et al. 2022;Bruton et al. 2023;Jones et al. 2023;Morishita et al. 2023), Ly emitter luminosity functions (Inoue et al. 2018;Morales et al. 2021), Ly emitter clustering (Sobacchi & Mesinger 2015;Ouchi et al. 2018), quasar damping wings (Davies et al. 2018a;Wang et al. 2020;Greig et al. 2022) and dark pixel fractions of the Ly forest (McGreer et al. 2015;Jin et al. 2023).The colour of the points highlights the spectacular redshift reach of JWST (black) compared to other sources (grey).
redshift evolution of the amplitude of the UV background are calculated from the outputs of the radiative transfer code.These data are then used as inputs to a new p-gadget-3 simulation, to capture the effect of a spatially fluctuating UV background due to patchy reionization in the hydrodynamic simulation.
The advantages to this hybrid approach are that it is computationally much cheaper than performing a full radiation-hydrodynamic simulation, while still allowing us to track the hydrodynamic response of the gas to reionization.We can also guarantee the reionization history of our final patchy simulation, which can be difficult to calibrate in self-consistent galaxy formation simulations of reionization.There are of course also some disadvantages to this method, such as the finite spatial and temporal resolution of our ionizing background, but we do not expect this to have a significant effect on the diffuse intergalactic gas that is our primary focus here.
Our fiducial reionization history is shown in Figure 1 (blue line) and compared against different probes of the IGM neutral fraction.Our preferred reionization history is a model where reionization ends at  = 5.3 (where we define the "end" of reionization as the redshift where the IGM first becomes 99.9 per cent ionized by volume).This is in agreement with the tight constraints on the end of reionization obtained from the opacity of the Ly forest (Bosman et al. 2022).This model is 50 per cent ionized at  = 7.1, which sits within the scatter of the various probes of the IGM neutral fraction at that redshift.To explore the effect of changing the reionization history on our simulated damping wings, we also analyse an alternative reionization history (shown by the pink line in Figure 1).In this early reionization model, the IGM is 99.9 per cent ionized at  = 6.6 and 50 per cent ionized at  = 7.8.As this early model is in tension with many of the constraints shown in Figure 1, we present it only as an example of how our results would change in a more extreme scenario, and do not advocate it as a preferred reionization history.

Simulated absorption spectra
To construct our mock damping wings, we extract lines of sight from simulation snapshots at  = 7, 8, 9 and 10.As we are interested in the IGM absorption arising in gas in the foreground of galaxies, we take these lines of sight through the 100 most massive haloes at each redshift, extracting the gas neutral hydrogen density, temperature and peculiar velocity.We orient these lines of sight in the ±, ± and ± directions around the halo and continue them for 20 ℎ −1 cMpc from the halo centre.This results in 600 sightlines in total.The width of each pixel in the sightline is 19.53 ℎ −1 ckpc for our fiducial 40 ℎ −1 cMpc 2048 3 -particle simulation volume and 78.13 ℎ −1 ckpc for our larger, lower resolution 160 ℎ −1 cMpc 2048 3 -particle volume.
As neutral gas at large distances can also contribute to the strength of the IGM damping wing, we stitch on five additional random sightlines length 40 ℎ −1 cMpc through the simulation volume.As we would like to account for the redshift evolution of the IGM along the line of sight, we make use of the on-the-fly sightlines that were extracted from the Sherwood-Relics simulation at redshift frequency Δ = 0.1.Using these data, we stitch on subsequent sightlines taken from increasingly lower redshifts, consistent with the redshift interval corresponding to the comoving distance between the halo position and the beginning of the new sightline.This results in final lines of sight that are each 220 ℎ −1 cMpc long.
From these sightlines, we then sum up the contributions of the gas in front of the halo to the Ly optical depth using the analytic approximation to the Voigt profile presented in Tepper-García (2006).Due to the simplified star formation prescription used in our simulations, we do not expect to recover the distribution of neutral gas within the circumgalactic medium of the halo.We further aim to model only the contribution of the IGM to the damping wing.We therefore exclude the contribution of any gas within the virial radius of the halo by default (see also Weinberger et al. 2018).As the redshift of the observed galaxies is known from their emission lines, we renormalise the gas peculiar velocity along the line of sight such that the halo has a velocity of 0 km s −1 .As the damping wing extends redward of the Ly source frame redshift, we also record the optical depth in the mock spectra redwards of the halo location up to a velocity corresponding to 200 ℎ −1 cMpc "behind" the halo, but of course only accounting for absorption by gas that is in front of the halo in position space.Examples of our simulated IGM damping wings are shown in Figure 2.

The Miralda-Escudé model for the IGM damping wing
As the goal of this work is to contrast the IGM damping wings constructed from cosmological simulations with the analytic Miralda-Escudé (1998) model, we recap the details of that model here for completeness.
It for a uniform IGM at redshift , where Λ  is the Ly decay constant,   is the Ly rest frame wavelength and  () is the Hubble parameter.At high redshift, this can be approximated as The background neutral hydrogen density is defined as Here,  crit is the critical density at  = 0, Ω b is the baryon density,  is the helium mass fraction and  H is the mass of a hydrogen atom.
Using Equation 1, we can then define the optical depth of the IGM damping wing, where  D is the optical depth along the line of sight evaluated at redshift ,  GP is the Gunn-Peterson optical depth defined above,  s is the redshift of the source,  b is the redshift of the edge of the ionized bubble and  n is the redshift where reionization is defined to end.  is defined as Λ    /(4), where  is the speed of light.Finally, the function  () is given by Note this function is only well defined for We therefore only compute the optical depth of the damping wing for redshifts  >  b , and at lower redshifts set the IGM optical depth to  GP ().

Comparison with late reionization model
We show the median and scatter of the IGM damping wings computed from Sherwood-Relics with the fiducial late reionization history with the blue lines and shaded regions in Figure 3.Each panel corresponds to a different redshift.We recover the expected result that as we move to higher redshifts and further into the epoch of reionization, the IGM damping wings become stronger due to the increasing volumeweighted average neutral fraction of the IGM. Figure 3 further shows the results of our comparison between the Sherwood-Relics late reionization model and the Miralda-Escudé (1998) analytic model (denoted by the dashed black line).To make a fair comparison between the two models, we compute the Miralda-Escudé (1998) model assuming parameters from our simulation.We measure the volume-weighted average neutral fraction from the simulation output at the redshift of the halo.Note that this will not take any evolution of the IGM along the line of sight into account, which is taken into account in an approximate manner in our simulated sightlines.We further assume the end of reionization to be the redshift where the IGM is 99.9 per cent ionized by volume.For the size of the ionized bubble, we measure the distance between the position of the halo where we begin our sightlines and the point where the neutral fraction of the gas first exceeds  HI = 0.5 along each line of sight.We then use the redshift corresponding to the median of these bubble sizes along all lines of sight as input to Equation 3. The values used are indicated on the bottom right of each panel of Figure 3.
We find that while the Miralda-Escudé (1998) model does an excellent job of recovering the median IGM damping wing predicted by the Sherwood-Relics simulation at velocities a few 100 km s −1 redward of the systemic velocity, it overpredicts the IGM transmission otherwise.The deviation between the analytic model and the simulations begins slightly redward of the systemic velocity and grows towards bluer wavelengths.The discrepancy is small at higher redshift when the ionized bubble sizes are small and the volume-weighted neutral fraction of the IGM is large, but becomes increasingly apparent towards lower redshift.For example, the Miralda-Escudé (1998) model at  = 7 calculated using the parameters from our simulation predicts that the IGM transmission should be 50 per cent at a velocity of −500 km s −1 , while in Sherwood-Relics no transmitted flux is expected at all at that velocity.

Comparison with additional reionization models
We next explore the IGM damping wings generated from different patchy reionization simulations from the Sherwood-Relics simulation suite.The results are shown in Figure 4 for simulation snapshots at redshifts  = 7 and 10.In each case, we overplot the results from the Miralda-Escudé (1998) damping wing model, calculated from the properties of the corresponding simulation.For easy comparison, the first column repeats the results from the fiducial late reionization model as already discussed in Section 3.1 and shown in Figure 3.The second column contrasts this with simulated IGM damping wings computed from the early reionization model.At fixed redshift, this early model has both lower volume-weighted average neutral fractions and larger bubble sizes than the late model, as indicated in the lower right of each panel of Figure 4.This results in weaker IGM damping wings in the early reioniziation model than in the late reionization model at a given redshift.The other notable difference is that at  = 7 in the early reionization model, there is some transmitted flux visible in the Ly forest blueward of the systemic redshift.However, even in this extreme reionization model, this flux is still much lower than the corresponding prediction for the IGM transmission of the Miralda-Escudé (1998) model.
We next explore how our results depend on the volume of our reionization simulation, as it has been shown that volumes of several hundred cMpc are required to fully capture the patchiness of reionization (Iliev et al. 2014), in contrast to the fiducial 40 ℎ −1 cMpc volume we analyse here.The third column of Figure 4 shows the simulated IGM damping wings computed from a 160 ℎ −1 cMpc simulation volume, with a near-identical reionization history to the late reionization 40 ℎ −1 cMpc volume plotted in the first column  over the redshift range probed here.We find that at a given redshift, although the volume-weighted average neutral gas fraction is the same in the two simulations, the bubbles are larger in the larger volume by approximately a factor of two.The result of this is that the sightlines generated from the larger volume show more transmission in the IGM damping wing profile at lower wavelengths, with even a small amount of transmission occurring just blueward of the systemic velocity.However, this effect is somewhat offset by the larger infall velocities associated with the more massive halo masses that are found in the larger volume.At  = 7, the median mass of the 100 most massive haloes in the 40 ℎ −1 cMpc volume is 9.0×10 10 ℎ −1  ⊙ , compared with 4.3 × 10 11 ℎ −1  ⊙ in the 160 ℎ −1 cMpc volume.This results in a larger offset between the systemic velocity and the median of the sharp cutoff in the Ly transmission in the larger volume than in the smaller volume.
To investigate how the two simulations compare for a more similar population of haloes, we take a different set of lines of sight from the 160 ℎ −1 cMpc volume.This time we select a sample of haloes that have a distribution of masses similar to the 100 most massive haloes in the 40 ℎ −1 cMpc volume.The results of this test are shown in the fourth column of Figure 4. We find that the median Ly transmission profile now looks more similar in the 160 ℎ −1 cMpc and 40 ℎ −1 cMpc volumes, as the velocity offset is now much closer between the two models.We also find a median bubble size that is smaller around the lower mass haloes.However, the scatter in bubble sizes is quite different, with the 68 (95) per cent range in bubble size a factor of two (three) larger in the larger volume simulation, most likely because these lower mass haloes are clustered around more massive haloes that can carve out larger ionized bubbles.The result of this is slightly weaker IGM damping wings in the 160 ℎ −1 cMpc compared to the 40 ℎ −1 cMpc volume, even when probing haloes of the same mass at fixed volume-weighted average neutral fraction.
However, in all cases, we still recover the trend that our simulated IGM damping wings are not described well by the Miralda-Escudé (1998) analytic model blueward of the systemic redshift, independent of the reionization history, simulation volume or host halo mass.

Reasons for difference between analytic and simulated damping wing models
We investigate the root of the difference between the IGM damping wings from the Sherwood-Relics simulations and the Miralda-Escudé (1998) model in Figure 5.The first column of this figure shows the IGM damping wings at  = 10 and  = 7.These are exactly the same as displayed in Figure 3, and are calculated as described in Section 2.2.The second column shows the IGM damping wings calculated for the same sightlines, but now we compute the Ly optical depth without taking the peculiar velocity of the gas into account.The effect of this is to produce a sharp cutoff in the IGM transmission at a velocity of 0 km s −1 , whereas previously this cutoff occurred at a range of different velocities redward of the systemic velocity, due to the range of motions of the local gas falling into the host halo.By disregarding this infalling gas and shifting the cutoff in the IGM transmission to bluer wavelengths, we find that we improve the agreement between the Sherwood-Relics damping wings and the Miralda-Escudé (1998) damping wing model such that there is near perfect agreement for all transmission redward of the halo redshift.However, we still find conflicting predictions for the level of transmission expected blueward of the halo redshift.Next, we investigate the impact of residual neutral gas in the ionized bubbles.We recompute the Ly optical depth along our lines of Absorption spectra as shown in Figure 3, calculated with all the relevant physical effects.Second column: Absorption spectra calculated without the effect of gas peculiar velocities.Third column: Absorption spectra calculated without any residual neutral hydrogen inside the ionized bubbles and neglecting the gas peculiar velocities.Fourth column: spectra calculated without accounting for any redshift evolution of the ionization state along the line of sight, without any residual neutral hydrogen inside the ionized bubbles and neglecting the gas peculiar velocities.
sight, but now assume that all of the gas within the ionized bubble surrounding the host halo is completely ionized.We achieve this by finding the first neutral island along our line of sight (which we define as a pixel with  HI > 0.5), and then setting the neutral fraction in all pixels between the halo and that first neutral island to have  HI = 0.The reason that we do not just remove the residual neutral gas in all bubbles, rather than just the host bubble, is that this will produce additional transmission along the line of sight.See, for example, the additional transmission peak separated from the IGM damping wing in Halo 6 of Figure 2.This would confuse our damping wing signal when we stack our absorption spectra, and hence our comparison with the Miralda-Escudé (1998) model.We therefore chose to isolate the effect of the residual neutral gas only in the vicinity of the host halo.
We show the median transmission profiles and their scatter in the third column of Figure 5, where we have now excluded both the residual neutral gas and neglected the peculiar velocities of the gas.We find that this makes a large difference to the predicted transmission blueward of the host halo redshift, and there is now very good agreement between the damping wings calculated from the Sherwood-Relics simulation and the Miralda-Escudé (1998) model.As has previously been noted in other works, residual neutral gas inside the ionized bubbles should play a large role (see, e.g., Mesinger &Haiman 2004 andBolton &Haehnelt 2007 for a discussion in the context of quasar proximity zones, and Mason & Gronke 2020 in the context of Ly emitting galaxies).Although the gas in the ionized bubbles in the Sherwood-Relics simulation is highly ionized, it still has a neutral fraction of order  HI ∼ 10 −3 (see Figure 2).From Equation 1, this will result in a Ly optical depth  Ly ∼ 500, more than enough to completely saturate the IGM absorption.
After removing both the effects of peculiar velocities and residual neutral gas, we find much better agreement between our simulated IGM damping wings and the Miralda-Escudé (1998) model.There still remains a small difference between the analytic model and the median IGM damping wing predicted from our simulations, with the simulations predicting slightly more transmission at velocities blueward of the host halo redshift.As a final test, we therefore investigate the effect of neglecting evolution of the ionization state of the gas along the line of sight.As described in Section 2.2, we account for this in an approximate way by stitching together lines of sight from simulation outputs at different redshifts.We check what effect this has by recomputing the lines of sight, but this time now stitching together lines of sight from a fixed redshift, such that there is no evolution in the average ionization field across the 220 ℎ −1 cMpc line of sight, although individual pixels will still be ionized or neutral depending on their location within the simulation volume.We find however that neglecting the evolution of the IGM along the line of sight has only a small effect on our simulated damping wings.In summary, we find that the most significant reason for the difference between the predictions for the IGM damping wing from the Sherwood-Relics simulations and the Miralda-Escudé (1998) model is the residual neutral gas within the ionized bubbles.While the effects of peculiar motions of infalling gas and evolution of the IGM along the line of sight towards the observer also play a small role, these are subdominant to the resonant absorption caused by the gas within the bubbles.

COMPARISON OF MODELS AND OBSERVATIONS
We next seek to make contact between the IGM damping wings constructed from our simulations of inhomogeneous reionization, and the the first observations of galaxy damping wings performed with JWST.In particular, we investigate the recently inferred 100 cMpc-sized ionized bubbles during reionization (Umeda et al. 2023) and the observations of proximate DLAs in high-redshift galaxies (Heintz et al. 2023).

Large ionized bubbles during reionization
We first turn our attention to the IGM damping wings measured from stacks of galaxies presented in Umeda et al. (2023).The sample consists of spectra of 26 galaxies at 7 <  < 12 obtained from multiple JWST programs (Arrabal Haro et al. 2023;Hsiao et al. 2023;Finkelstein et al. 2023).These galaxies were divided into four redshift bins with median redshifts  = 7.140, 7.452, 7.960 and 9.801.To see the imprint of the IGM damping wing, it was necessary to compare the observed spectrum with a model for the intrinsic spectrum.To properly account for absorption in the circumgalactic medium of the galaxies and isolate the effect of the IGM, Umeda et al. (2023) used stacks of galaxies in bins of stellar mass at 2.5 <  < 5 from the VANDELS survey (Cullen et al. 2019).By choosing the template in the stellar mass bin that best matched the observed galaxy continuum far redward of the Ly source frame wavelength, and finding the best-fitting Miralda-Escudé (1998) damping wing model in the wavelength region close to Ly, they obtained constraints on the redshift evolution of both the volume-weighted average neutral fraction of the IGM and the sizes of the ionized bubbles.
As we wish to compare the IGM damping wings generated from the Sherwood-Relics simulations, we first repeat the process to generate the template intrinsic spectra.We use the stacked spectra from the stellar mass bins that Umeda et al. (2023) found to best fit their stacked spectra.Namely, we use the log( ★ / ⊙ ) = 8.16−8.70stack at  = 7.140, the log( ★ / ⊙ ) = 8.70 − 9.20 stack at  = 7.452 and 7.960 and the log( ★ / ⊙ ) = 9.50 − 9.65 stack at  = 9.801.We renormalise these spectra by a factor (/1800Å)  , where  and  are parameters that can be varied until the 2 between the template stacks and the observed stacks is minimised.Following Jones et al. (2023), to account for the low spectral resolution of the NIR-Spec PRISM mode 2 , we convolve the spectra with a Gaussian of an appropriate width.We further bin the template spectra to pixels of 25 Å in the rest frame, to match the observed stacks.As in Umeda et al. (2023), we perform this fitting the wavelength range between the wavelength where an IGM damping wing generated with the Miralda-Escudé (1998) model for a neutral IGM with no ionized bubble first reaches a transmission of 90 per cent at the blue end, and 2200 Å rest frame wavelength at the red end.The resulting template spectra are shown in orange in the left panel of Figure 6 (the thin lines are before convolution and binning, and the thick lines are afterwards) and can be compared with the stacked observed spectra from Umeda et al. (2023) shown in black.Having generated the template intrinsic spectra, we then convolve these with models for the IGM damping wing.We first use the best-fitting parameters measured in Umeda et al. (2023), which are volumeweighted average neutral fractions ⟨ HI ⟩ v = (0.46, 0.54, 0.63, 0.83) and bubble sizes  b = (149, 96.1, 16.5, 5.04) cMpc at redshift  = (7.140, 7.452, 7.960, 9.801).We show the resulting curves generated from the galaxy template, convolved with these damping wing models, as the dark red curves in Figure 6.
We next repeat this process, but instead now use the simulated IGM damping wings generated from the Sherwood-Relics late reionization simulation.We generate spectra from the snapshots at  = 7, 7.5, 8 and 10, but rescale the densities and velocities to match the median redshift in each bin of the Umeda et al. (2023) sample.The resulting galaxy templates convolved with the simulated damping wings are shown by the blue curves and light (dark) shaded regions in Figure 6, which represent the median and 68 (95) per cent scatter of the distribution.We find that in the two highest redshift bins, at  = 9.801 and 7.960, there is good agreement between the predictions from the best fit IGM damping wing from Umeda et al. (2023) and the Sherwood-Relics simulations.We do find that the Sherwood-Relics simulation seems to underpredict the strength of the damping wing at  = 9.801, but this is also seen in the Umeda et al. ( 2023) IGM damping wing model, and may be due to the influence of proximate absorbers as discussed in Section 4.3.
We find, however, that there is a significant difference in the expected transmitted flux between the Sherwood-Relics damping wings and the best fit Umeda et al. ( 2023) IGM damping wing model at  = 7.140 and 7.452.The intrinsic galaxy template convolved with the Sherwood-Relics model produces too little transmitted flux in the pixels close to the Ly source frame wavelength.At first glance, this suggests that the damping wings in Sherwood-Relics may be too strong when compared with the observations.Indeed, when comparing the best fit ionized bubble sizes from Umeda et al. (2023) with the median bubble sizes around the hundred most massive haloes in Sherwood-Relics, we find that at  = 7.140 (7.452) the bubble sizes are a factor 13 (11) smaller than the measurements from Umeda et al. (2023).The large bubble sizes measured in that work are driven by the significant detections of transmitted flux in pixels blueward of the Ly source frame wavelength.Fitting for this with the Miralda-Escudé (1998) model for the IGM damping wing allows for this detected flux when the ionized bubbles around the host galaxies are large.
However, based on the results of Section 3, we do not expect to observe any transmitted flux blueward of the Ly source frame wavelength, independent of the size of the ionized bubble, as a result of the residual neutral hydrogen in the ionized IGM.Furthermore, if there were indeed large, 100 cMpc-scale regions of the IGM that were transparent to Ly scattering already around  ∼ 7 galaxies, one would expect to see these transmissive regions in the  > 6.5 Ly forest of the highest redshift quasars, but this is not the case (Jin et al. 2023).An exception would be the highly ionized proximity zones observed around high-redshift quasars, but even around the brightest  ∼ 6 quasars, the sizes of these proximity zones are smaller than the sizes of ionized bubbles measured by Umeda et al. (2023) at  ∼ 7 (Eilers et al. 2017).The amplitude of the transmitted Ly forest flux proposed by Umeda et al. (2023) is also unexpected.Using their best-fit parameters in the Miralda-Escudé (1998) model at  = 7.140, we can measure the mean flux in the galaxy's Ly forest in two 50 ℎ −1 cMpc chunks that contain sections of the ionized bubble.We measure a mean flux ⟨⟩ = (0.97, 0.69) in the two chunks, as we move blueward from the galaxy.This can also be expressed as an effective optical depth  eff = (0.03, 0.37).Values as high as this are not observed in the Ly forest of quasars until below  = 3 (Becker et al. 2013).
It is therefore of interest to investigate other scenarios that could produce transmitted flux blueward of the galaxy's Ly emission.One possibility is that the flux and wavelength calibration of JWST PRISM spectra still present challenges.Residual errors from these may result in additional spurious flux blueward of Ly.The choice of  (Cullen et al. 2019), rescaled to match the continua of the higher redshift galaxies.The thick orange lines show the same spectra but after convolution and binning to coarser pixels.The thick red line shows the template intrinsic galaxy spectra convolved with the best-fitting IGM damping wing model of Umeda et al. (2023).The blue lines and shaded region show the intrinsic galaxy spectra convolved with absorption spectra generated from the Sherwood-Relics (late) simulation.Right column: As before, but now using a template spectrum where we have changed the strength of the Ly emission line by hand, such that it is now increased in amplitude by a factor of three.
the intrinsic galaxy spectrum may also play a role, and this possibility is discussed in Section 4.2 below.

Influence of the intrinsic spectrum on inferred bubble sizes
An alternative explanation for these observations may arise from the chosen intrinsic galaxy template spectrum.If the assumed intrinsic spectrum used is not representative of the true intrinsic spectrum, this may lead to a change in the reionization parameters recovered from the observations.In particular, an underestimate of the strength the Ly emission line will effect not just the pixel that line falls in, but also the neighbouring pixels.This is a result of the low resolution of the NIRSpec PRISM mode, which has  = Δ  ∼ 30 − 40 in the redshift range  = 7−10.We investigate the effect of underestimating the intrinsic Ly emission of the galaxy in the right column of Figure 6.We repeat the process of convolving the Sherwood-Relics IGM damping wings with a template galaxy spectrum.We again assume the same stacks of lower redshift spectra from Cullen et al. (2019), but we now adjust the strength of the Ly line such that it is a factor of three stronger.
This factor of three was chosen arbitrarily as a value to demonstrate the effect of changing the intrinsic spectrum.The Ly equivalent width distribution in the sample used to make the low-redshift stacks has a large scatter, with not all galaxies showing Ly emission and some showing absorption instead.In Cullen et al. (2020), which analysed an expanded sample compared to Cullen et al. (2019), the median rest-frame equivalent width of the sample was  Ly = −4 Å.The strongest Ly emitting galaxies had  Ly = 110 Å. Shibuya et al. (2018) presented observations of  ∼ 6.6 galaxies with equivalent widths  Ly > 200 Å.These observations will also contain the effects of attenuation by the IGM, so the intrinsic Ly emission of these galaxies will be stronger.Multiplying the Ly emission in these stacked spectra by a factor of three to produce intrinsic Ly equivalent widths of at most  Ly ∼ 330 Å in individual galax-ies, and lower values in the stacked spectra, is therefore not totally unreasonable.Indeed, a Ly emitting galaxy with  Ly ≈ 400 Å at  = 7.278 has been discovered with JWST (Saxena et al. 2023).We note that we only increase the amplitude of the line, and do not change its shape, but this should not make a difference given the wide pixels used in the stacked spectra.
The resulting template is shown by the thin purple lines in Figure 6 before convolution and binning, and by the thick purple curves afterwards.Comparing the thick orange and thick purple curves in the figures on the left and right columns, we find that, as expected, increasing the intrinsic Ly emission enhances the intrinsic flux not only in the pixel that corresponds to the source frame Ly wavelength, but also in adjacent pixels, due to the broad instrument profile of the NIRSpec PRISM mode (see also Jones et al. 2023).We further find that we now see an excess of flux in the pixel containing the Ly emission line in our template intrinsic spectrum compared to the observations at  = 7.140 and 7.452.This is in contrast to the template used in Umeda et al. (2023), where the expected and observed flux at Ly were nearly identical, due to the large ionized bubbles that were assumed.If this latter option were the case, it would be a different scenario from what is observed in the spectra of  = 7 − 7.5 quasars, where a significant fraction of the Ly emission line is expected to be absorbed (e.g., Davies et al. 2018a;Wang et al. 2020;Greig et al. 2022).
We find that when we convolve the Sherwood-Relics IGM damping wings with this new template for the intrinsic galaxy spectrum, we now find significant transmitted flux blueward of the Ly source frame wavelength, mimicking the effect of invoking the large ionized bubbles in the Miralda-Escudé (1998) model.This effect is largest in the lowest redshift bins, as the IGM damping wing in the simulations is strong enough at  ≳ 8 to absorb the bulk of the Ly emission.Of course, we have here arbitrarily changed the strength of the intrinsic Ly emission by hand to demonstrate the expected effect.In practice, it should be possible to obtain an estimate of the expected intrinsic Ly emission from the Balmer emission lines of the galaxy (e.g., Hayes 2015).
We further note that we found near-indistinguishable results when comparing the Umeda et al. (2023) observations to the Sherwood-Relics late and early models (plotted in Figure 1).Although there are differences in the damping wings we compute from these two simulations, as shown in Figure 4, the broad instrument profile of the NIRSpec PRISM mode and large pixel size of the stacked spectra makes the two models difficult to differentiate when comparing to the JWST data, and the change in reionization history was much less significant than the change in the intrinsic spectrum.This suggests that it may be difficult to constrain reionization with the currently published spectra, and that larger samples of galaxies will be required.

Proximate high column density absorbers
We next analyse our simulated IGM damping wings in the context of the detection of damped Lyman- absorbers in three galaxies at  = 8 − 11 (Heintz et al. 2023).The spectra of these three galaxies were compared with IGM damping wings calculated using the Miralda-Escudé (1998) with volume-weighted average neutral fraction ⟨ HI ⟩ v = 0.1, 0.5 and 0.1 and assuming an ionized bubble radius  b = 0 (i.e., setting  b equal to  s in Equation 3).Under these assumptions, it was shown that IGM damping wing profiles were not strong enough to explain the observed damping wings.Further comparison with simulated IGM damping wings from Laursen et al. (2019) also indicated that the IGM alone could not explain the observations.However, as noted in that work, the neutral fraction in that simulation is somewhat low, with ⟨ HI ⟩ v = 0.13 at  = 8.8.Heintz et al. ( 2023) concluded that including strong proximate DLAs with H i column densities in the range  HI = 10 22.1 − 10 22.4 cm −2 provided a much better fit to the observations.We test this theory in the context of the IGM transmission curves generated from the Sherwood-Relics late reionization model, which as described in Section 3.3 contain several physical effects that are not present in the Miralda-Escudé (1998) model and which has a reionization history in good agreement with estimates of the evolution of the IGM neutral fraction (Figure 1).We note that we do not expect to produce DLAs in our simulations, due to our simplified star formation model and lack of galactic feedback.We further note that, as described in Section 2.2, we are also not including the contribution of any gas inside the virial radius of the galaxy when we compute our IGM damping wings.We focus our comparison on CEERS-43833, a galaxy at redshift  = 8.7622 (Arrabal Haro et al. 2023;Finkelstein et al. 2023), as Heintz et al. (2023) show a version of this spectrum where the transmission has been normalised by the intrinsic transmission, and so can be easily compared with our simulated damping wings.
We compute our absorption spectra starting from haloes extracted from a simulation snapshot at  = 9, but rescale the gas density and velocity assuming the host halo is at the redshift of CEERS-43833.The simulated IGM transmission curves we calculate are shown in blue in the top panel of Figure 7 and can be compared with the normalised spectrum of CEERS-43833 shown in black.As in Section 4.1, we further convolve the spectra with a Gaussian with a width corresponding to the resolution at the Ly wavelength at the redshift of CEERS-43833.The resulting transmission curves are shown in red in Figure 7.We find that even after convolution, the absorption from the IGM alone is not enough to reproduce the damping wing observed in CEERS-43833.This is similar to the simulated IGM damping wings of Laursen et al. (2019), despite the later reionization model we analyse here.This supports the claim of Heintz et al. ( 2023) that there are proximate DLAs associated with these high-redshift star-forming galaxies.Indeed, when we include the contribution of a DLA with column density  HI = 10 22.1 cm −2 at the redshift of the host galaxy, we find a much better fit to the observed spectrum, as shown in the bottom panel of Figure 7.

Implications
It is clear from these JWST observations that IGM damping wings seen in high-redshift galaxies open up many exciting possibilities for studying the first half of the reionization epoch.However, this observable comes with its own set of challenges.Galaxy spectra are expected to be somewhat simpler to model than quasar spectra as, e.g., the quasar spectra can be complicated by strong N v emission close to the Ly emission (Davies et al. 2018b;Greig et al. 2022).However, the strength of the intrinsic Ly emission will still be important to capture, especially in analyses of low resolution spectra.Moving forward, it may be useful to marginalise over grids of galaxy template spectra convolved with a separate grid of IGM damping wing models, in a spirit similar to what has been carried out in analyses of quasar damping wings.However, as pointed out by Umeda et al. (2023), it will be necessary to understand how to incorporate absorption by the circumgalactic medium of the host galaxy in such a setup.
The incidence rate of proximate DLAs in these galaxies will also add a source of uncertainty.Indeed, two of the three galaxies shown to host DLAs in Heintz et al. (2023) (MACS0647-JD at  = 10.170 and CEERS-16943/Maisie's Galaxy at  = 11.409) are also present in the sample analysed by Umeda et al. (2023).In the future, it will be important to either allow for the contribution of local, high-density neutral hydrogen to the observed damping wings when fitting for the IGM component, as has been done for GRBs (Totani et al. 2006).Alternatively, it may be possible to remove the contaminated objects from the sample, by searching for metal lines redward of Ly that are associated with the DLA.However, depending on the metallicity of these objects, this may have to wait for the upcoming spectrographs on 30-metre class telescopes, such as ANDES or HARMONI.It may also be possible to remove the DLA hosts statistically based on the properties of the galaxy, once we have accumulated a large enough sample of spectra from these  > 7 galaxies.
At the same time, it will be important to continue to make advances on the theoretical side to properly interpret these results.For example, it is likely that these simulations would struggle to reproduce observations of the (rare) double-peaked Ly emitting galaxies at  ∼ 6.5 (Hu et al. 2016;Songaila et al. 2018;Meyer et al. 2021), which may call for more efficient production and escape of ionizing photons than we have assumed in the source model for our simulations or the contribution of AGN in locally enhancing the ionization state of the IGM (Bosman et al. 2020).It will further be important to increase the dynamic range of the simulations, such that larger volumes can be modelled while still resolving the small-scale structure of the IGM, to bridge the gap between our models and the bubble sizes predicted in large semi-numerical simulations of reionization (Lu et al. 2023) to fully understand the implications for Ly transmission.

CONCLUSIONS
We have presented here an analysis of mock IGM damping wings generated from the Sherwood-Relics simulation suite.We have compared our simulated damping wings to the analytic IGM damping wing model of Miralda-Escudé (1998).We found excellent agreement between the simulated and analytic damping wings redward of Ly, but found that the agreement was poor on the blue side of Ly, with the analytic model predicting much more IGM transmission than expected from the simulated damping wings.We showed that this was a result of the residual neutral hydrogen inside in the ionized bubbles in our simulations of inhomogeneous reionization, which is enough to saturate the IGM absorption.
We further compared the simulated damping wings against recent observations of damping wings in high-redshift galaxies seen with JWST.We found that the simulated damping wings were unable to reproduce the significant detections of transmitted flux blueward of the galaxy Ly emission, which have led to claims of 100 cMpc-sized ionized bubbles at  ∼ 7 − 7.5.We suggest that as the residual neutral gas in the IGM at  ∼ 7 is expected to saturate the IGM absorption, the flux that has been observed blueward of Ly may instead be a result of stronger intrinsic Ly emission, which is spread into nearby pixels by the instrument profile of the NIRSpec PRISM mode.We also investigated claims of observed proximate DLAs in high-redshift galaxies, and confirmed that the strong damping wings observed in these galaxies cannot be reproduced by the IGM damping wings generated from our simulations.
Of course, we are only witnessing the first attempts at constraining the first half of reionization with galaxy damping wings.The first results are already incredibly impressive, and show huge promise for revealing the timing of reionization, even before the presence of this neutral gas is directly detected through its 21cm transition in upcoming experiments.However, this work has demonstrated that the exquisite quality of the new JWST data demands comparison with models that include all of the relevant physics.Future progress in this area will come from advances on both the observational and theoretical sides in understanding the interactions between galaxies and the intergalactic medium in the reionization epoch.RTI 4002, and by the Max Planck Society through a Max Planck Partner Group.Support by ERC Advanced Grant 320596 'The Emergence of Structure During the Epoch of Reionization' is gratefully acknowledged.We thank Kasper Heintz for answering a question about the spectrum of CEERS-43833.We thank Volker Springel for making p-gadget-3 available.We also thank Dominique Aubert for sharing the aton code.For the purpose of open access, the author has applied a Creative Commons Attribution (CC BY) licence to any Author Accepted Manuscript version arising from this submission.

Figure 2 .
Figure 2.Examples demonstrating the diversity of our simulated IGM damping wings at  = 8.The map on the top shows the neutral fraction of the gas in a 19.53 ℎ −1 ckpc slice through our simulation.The white circles mark the locations of the six haloes plotted below, with the numbers identifying the different haloes, connecting them to the large-scale ionization field of our reionization simulations.The arrow in the bottom right shows the direction along which the spectra were calculated.The six subplots on the bottom two rows show the results from lines of sight through different haloes.The top panel of each subplot shows the neutral gas fraction along the line of sight.The regions intersecting islands of neutral gas (which we define as having a neutral fraction  HI > 0.5) are highlighted in red.The bottom panel of each subplot shows the corresponding Ly absorption spectrum.The solid black lines show the spectra calculated taking the contribution of all gas (both in ionized bubbles and neutral islands) into account.The dashed red line shows the spectra calculated assuming the gas in the bubbles is completely ionized, and only the neutral islands contribute to the Ly optical depth.

Figure 3 .
Figure 3.Comparison of the simulated Sherwood-Relics damping wings for the late reionization history with the Miralda-Escudé (1998) model.The different panels show results at different redshifts.The blue solid line shows the median Ly transmission for the simulated sightlines as a function of velocity.The light and dark shaded regions encompass the 68 and 95 per cent scatter of the transmission, respectively.The dashed black line is the Miralda-Escudé (1998) model computed assuming the volume-weighted H i fraction and the median bubble size in the simulation as indicated on each panel.

Figure 4 .
Figure 4. Comparison of the simulated damping wings from different simulations of inhomogeneous reionization from the Sherwood-Relics suite with the Miralda-Escudé (1998) model.The different rows show results at different redshifts (top:  = 10, bottom:  = 7).First column: Absorption spectra as shown in Figure 3, calculated from the fiducial late reionization simulation in the high resolution 40 ℎ −1 cMpc volume.Second column: Absorption spectra calculated from the early reionization simulation in the high resolution 40 ℎ −1 cMpc volume.Third column: calculated from the fiducial late reionization simulation in the lower resolution 160 ℎ −1 cMpc volume.Fourth column: Absorption spectra calculated from the fiducial late reionization simulation in the lower resolution 160 ℎ −1 cMpc volume for a similar halo mass range as in the 40 ℎ −1 cMpc volume.

Figure 5 .
Figure 5.Comparison of the simulated Sherwood-Relics damping wings with the Miralda-Escudé (1998) model.The different rows show results at different redshifts (top:  = 10, bottom:  = 7).The different columns compute the absorption spectra with and without different physical effects.First column:Absorption spectra as shown in Figure3, calculated with all the relevant physical effects.Second column: Absorption spectra calculated without the effect of gas peculiar velocities.Third column: Absorption spectra calculated without any residual neutral hydrogen inside the ionized bubbles and neglecting the gas peculiar velocities.Fourth column: spectra calculated without accounting for any redshift evolution of the ionization state along the line of sight, without any residual neutral hydrogen inside the ionized bubbles and neglecting the gas peculiar velocities.

Figure 6 .
Figure 6.Left column: Comparison of modelled IGM damping wings and stacks of galaxy spectra from Umeda et al. (2023).stacked observed galaxy spectra are shown in black.Each row corresponds to a different redshift, at  = 7.140, 7.452, 7.960 and  = 9.801.The thin orange lines show stacked spectra at lower redshift from the VANDELS survey(Cullen et al. 2019), rescaled to match the continua of the higher redshift galaxies.The thick orange lines show the same spectra but after convolution and binning to coarser pixels.The thick red line shows the template intrinsic galaxy spectra convolved with the best-fitting IGM damping wing model ofUmeda et al. (2023).The blue lines and shaded region show the intrinsic galaxy spectra convolved with absorption spectra generated from the Sherwood-Relics (late) simulation.Right column: As before, but now using a template spectrum where we have changed the strength of the Ly emission line by hand, such that it is now increased in amplitude by a factor of three.

Figure 7 .
Figure 7.Comparison of the normalised spectrum of CEERS-43833 (a galaxy at redshift  = 8.7622) taken from Heintz et al. (2023) and the IGM damping wings from the Sherwood-Relics late reionization simulation.The black line in both panels shows the spectrum of CEERS-43833 and the grey shaded regions show the corresponding error.Top: The blue line and shaded regions are the median and scatter in IGM transmission curves calculated from the simulation.The red line and shaded regions are the same curves, but after convolution with the NIRSpec PRISM instrument profile.Bottom: As above, but now the simulated damping wings are computed for both the contribution from the IGM and for a proximate DLA with column density  HI = 10 22.1 cm −2 , as estimated by Heintz et al. (2023).
(Curtis-Lake et al. 2023;Hsiao et al. 2023;Umeda et al. 2023)er to a selection of constraints from the literature on the redshift evolution of the volume-averaged neutral fraction of the IGM.The solid blue line shows our fiducial model, where reionization ends at  = 5.3.The solid pink line shows a more extreme early model, where reionization ends at  = 6.6.For comparison, we show reionization constraints from observations of galaxy damping wings(Curtis-Lake et al. 2023;Hsiao et al. 2023;Umeda et al. 2023), Ly equivalent widths