How to Search for Islands of Neutral Hydrogen in the z\sim 5.5 IGM

# How to Search for Islands of Neutral Hydrogen in the z∼5.5 Igm

Matthew Malloy11affiliation: Department of Physics & Astronomy, University of Pennsylvania, 209 South 33rd Street, Philadelphia, PA 19104, USA & Adam Lidz11affiliation: Department of Physics & Astronomy, University of Pennsylvania, 209 South 33rd Street, Philadelphia, PA 19104, USA
###### Abstract

Observations of the Lyman-alpha (Ly-) forest may allow reionization to complete as late as , provided the ionization state of the intergalactic medium (IGM) is sufficiently inhomogeneous at these redshifts. In this case, significantly neutral islands may remain amongst highly ionized gas with the ionized regions allowing some transmission through the Ly- forest. This possibility has the important virtue that it is eminently testable with existing Ly- forest data. In particular, we describe three observable signatures of significantly neutral gas in the IGM. We use mock quasar spectra produced from numerical simulations of reionization to develop these tests. First, we quantify how the abundance and length of absorbed regions in the forest increase with the volume-averaged neutral fraction in our reionization model. Second, we consider stacking the transmission profile around highly absorbed regions in the forest. If and only if there is significantly neutral gas in the IGM, absorption in the damping wing of the Ly- line will cause the transmission to recover slowly as one moves from absorbed to transmitted portions of the spectrum. Third, the deuterium Ly- line should imprint a small but distinctive absorption feature slightly blueward of absorbed neutral regions in the Ly- forest. We show that these tests can be carried out with existing Keck HIRES spectra at , with the damping wing being observable for and the deuterium feature observable with additional high-resolution spectra for .

cosmology: theory – intergalactic medium – large scale structure of universe

## 1. Introduction

It has been nearly half a century since Gunn & Peterson (1965) pointed out that the lack of prominent absorption troughs, blueward of the Ly- line in quasar spectra, implies that intergalactic hydrogen is highly ionized. Only in the year 2001 were complete “Gunn-Peterson” absorption troughs finally revealed in the Ly- forest of high redshift () quasars discovered using the Sloan Digital Sky Survey (SDSS) (Fan et al., 2001; Becker et al., 2001; Djorgovski et al., 2001). Although these prominent absorption troughs were discovered more than a decade ago, the precise interpretation of the observations, and their implications for the reionization history of the universe, remain unclear. One difficulty here relates to the large optical depth to Ly- absorption: near , the optical depth is in a fully neutral IGM at the cosmic mean density (Gunn & Peterson, 1965). Based on this, it is common to infer that the IGM must be highly ionized below , at which point quasar spectra do show some transmission through the Ly- line. In addition, it is clearly hard to discern whether the gas above – that does show complete absorption in the Ly- line – is mostly neutral or is only neutral at the level of about one part in ten-thousand or so (e.g. Fan et al. 2006); in either case, the Ly- line will be completely absorbed.

However, if reionization is sufficiently inhomogeneous and ends late, there may be some transmission through the Ly- forest before reionization completes (Mesinger, 2009; Lidz et al., 2007). Theoretical models of reionization show that the IGM during reionization resembles a two-phase medium, containing a mixture of highly ionized bubbles along with mostly neutral regions. The ionized bubbles grow and merge, eventually filling essentially the entire volume of the IGM with ionized gas; the redshift at which this process completes is highly uncertain and still awaits definitive empirical constraint. In principle, the ionized bubbles may allow transmission through the Ly- forest even when some of the IGM volume is still in fact filled by neutral regions, i.e., before reionization completes. This calls into question the conventional wisdom described above – that the presence of transmission through the forest necessarily implies reionization completed by (Lidz et al., 2007; Mesinger, 2009) – strictly speaking, this conclusion follows only in the unrealistic case of a homogeneously-ionized IGM.

Indeed, some portions of the Ly- forest are completely absorbed, while other portions of the forest at these redshifts show transmission through the Ly- line. Quantitatively, if one counts only the fraction of pixels with some transmission through the forest as “certain to be ionized”, the volume-averaged neutral hydrogen fraction need only be smaller than at , and smaller than at (McGreer et al., 2011). These constraints are conservative since even mostly-ionized gas will give rise to some completely absorbed regions at these redshifts, but it is nevertheless interesting to ask whether some of the absorbed regions could in fact come from remaining “islands” of mostly neutral hydrogen gas in the IGM. The dark pixel fraction constraints of McGreer et al. (2011) certainly leave plenty of parameter space open for reionization completing at .

In fact, there are hints – albeit indirect ones – that significant amounts of neutral gas may remain in the IGM at these late times and so we believe that investigating this possibility amounts to more than closing a remaining “loophole” in the analysis of the Ly- forest. For example, recent measurements of the rest-frame ultraviolet galaxy luminosity function suggest a relatively low ionizing emissivity at , even for seemingly generous assumptions about the escape fraction of ionizing photons () and allowing significant extrapolations down the faint end of the luminosity function; e.g. the preferred model of Robertson et al. (2013) (that matches these observations) has at . In addition, the fraction of Lyman-break galaxies with detectable Ly- emission lines shows evidence for a rapid drop between which may require a significant neutral fraction at (e.g., Schenker et al. 2012; Pentericci et al. 2014, although see Bolton & Haehnelt 2013; Taylor & Lidz 2013). The inferred neutral fraction here would be easier to accommodate if there is still some neutral gas at . Furthermore, Becker et al. (2014) recently discovered an impressive Mpc/ dark region in the Ly- forest. This may result from an upward opacity fluctuation – driven by a fluctuating ultraviolet radiation field in a mostly ionized IGM – but this striking observation invites contemplating the more radical possibility that diffuse neutral regions remain in the IGM at this late time. Finally, Mesinger & Haiman (2007) and Schroeder et al. (2013) argue that the proximity zones of quasars at show evidence for damping wing absorption and a significant neutral fraction, further motivating the search for neutral gas at slightly later times.

Perhaps more importantly, we can design robust observational tests for the presence of neutral islands in the IGM, and either definitively detect neutral hydrogen at these redshifts, or significantly improve on the existing upper limits from McGreer et al. (2011). Towards this end, we study three possible tests for identifying neutral islands in the IGM, each of which can be applied using existing Ly- forest spectra. The presence of some transmission through the Ly- forest at allows us to consider tests that can not be applied at still higher redshift where the forest is completely absorbed (asides for in the “proximity zones” close to the quasar itself). We develop these tests using mock quasar spectra extracted from the numerical reionization simulations of McQuinn et al. (2007a). The first test we consider has been studied before (e.g. Fan et al. 2006, Mesinger 2010, McGreer et al. 2011), but is the most model dependent: the abundance and size distribution of “dark gaps”, i.e., regions of saturated absorption in the Ly- forest. Here we focus on the plausible impact of inhomogeneous reionization on the dark gap statistics. The second test utilizes the fact that the natural line width of the Ly- line gives rise to extended damping wing absorption, in the case that highly neutral gas is present in the IGM (Miralda-Escude, 1998). As a result, the transmission recovers more slowly around significantly neutral absorbed regions than around absorbed yet ionized regions. We find that this signature can be detected in partly neutral models by examining the stacked profile around extended absorbed regions. Note that, in contrast to previous work, here we propose to search for the damping wing signature in typical regions of the IGM, as opposed to in the proximity zones of quasars (Mesinger & Haiman 2007; Schroeder et al. 2013), or redward of Ly- at the source redshift. Our third test involves the stacked profile of extended absorbed regions in the Ly- forest. If these regions are significantly neutral, there should be a feature from absorption in the deuterium Ly- line just blueward (but not redward) of absorbed regions.

The outline of this paper is as follows. In §2, we briefly discuss which range of (volume-averaged) neutral fractions are physically plausible at . In §3 we describe the simulations used and the process for generating mock spectra. We discuss how the dark gap size distribution may be used to constrain the neutral fraction in §4. In §5, we describe how quasar spectra may be stacked in order to reveal the presence of deuterium and HI damping wing absorption in an idealized scenario and discuss adapting this approach for more realistic spectra in §6. We then apply this approach to mock quasar spectra in §7, discuss the constraining power of the stacking approaches in §8, and conclude in §9. Throughout, we consider a CDM cosmology parametrized by , , , , , and , (all symbols have their usual meanings), broadly consistent with recent Planck constraints from Ade et al. (2013).

## 2. Viability of Transmission Through a Partially Neutral IGM

Ideally, this study would make use of mock Ly- forest spectra extracted from fully self-consistent simulations of reionization, in which the efficiency of the ionizing sources and other relevant parameters are tuned so that reionization completes at . Unfortunately, large-scale reionization simulations that simultaneously resolve the properties of the gas distribution, as well as the sources and sinks of ionizing photons, while capturing large enough volumes to include a representative sample of the ionized regions, are still quite challenging. Here, we instead explore more approximate, yet more flexible, models. As we describe in more detail in the next section, we make use of the reionization simulations of McQuinn et al. (2007a) to describe the size and spatial distribution of the ionized and neutral regions during reionization. Inside of the ionized regions, we rescale the simulated photoionization rates, adjusting the intensity of the UV radiation field to match the observed mean transmitted flux through the Ly- forest. For simplicity, we assume that the intensity of the UV radiation field in the ionized regions is uniform and comment on the possible impact of this approximation where relevant.

Before proceeding further, however, it is worth considering which (volume-averaged) neutral fractions are physically plausible near . In order to get transmission through the Ly- forest, at least some of the hydrogen needs to be highly ionized. This requires the mean free path of ionizing photons to be relatively large, although we should keep in mind that the attenuation length will vary spatially during and after reionization, and so this quantity needs to be large only across some stretches of the IGM. This in turn demands some minimum separation between the neutral islands, because otherwise the neutral islands themselves will limit the mean free path and prevent a sufficiently intense UV radiation field from building up between the islands. Hence, it may be inconsistent to have remaining neutral islands in the IGM, yet still have some transmission through the Ly- forest. Here we briefly quantify this reasoning; we will be content with only a rough estimate, as our focus here is more on designing empirical tests. Further theoretical exploration here might be valuable, however, perhaps along the lines of Xu et al. (2014).

Quantitatively, previous studies infer that a photoionization rate on the order of is required to match the mean transmitted flux in the Ly- forest (e.g., Kuhlen & Faucher-Giguère 2012, Bolton & Haehnelt 2007).111Theses studies assume that reionization is complete at these redshifts. If the universe is in fact partly neutral, then a higher photoionization rate should be required in the ionized regions. In our rough estimate here, we neglect this given the other significant uncertainties involved. If we demand that the photoionization rate between the neutral islands needs to be in this ballpark to allow transmission through the forest, we can translate this into a required minimum average separation between the neutral islands, given an assumed ionizing emissivity. The average ionizing emissivity is likely on the order of photons per atom per Gyr (e.g. Bolton & Haehnelt 2007). This is close to the value required simply to balance recombinations and maintain the ionization of the IGM at the redshifts of interest. This emissivity is also comfortable with that inferred from the above photoionization rate and measurements of the mean free path to ionizing photons (Bolton & Haehnelt 2007, although Becker & Bolton 2013 recently argued for a slightly larger value), as well as the UV emissivity implied by measurements of the galaxy luminosity function (e.g. Robertson et al. 2013).

In this context, it is useful to note that:

 ΓHI =εHIσHI,limλmfpββ+1.5, (1)

where is the average proper ionizing emissivity, is the photoionization cross section at the Lyman limit, is the mean free path of ionizing photons at the Lyman limit, and is the intrinsic, unhardened spectral index of the ionizing radiation. This expression assumes that the mean free path to ionizing photons propagating through a clumpy IGM scales as (Zuo & Phinney, 1993). Inserting typical numbers we find:

 ΓHI =5.0×10−13sec−1[ε%HI3 photons/atom/Gyr] (2) ×[β2][3.51.5+β][1+z6.5]3[λmfp9.1 pMpc]. (3)

In other words, to get transmission through the forest for plausible values of the ionizing emissivity, we require the mean separation between neutral islands to be . This is a minimal requirement in that it assumes the neutral islands set the mean free path, when in fact Lyman limit systems and cumulative absorption in the mostly ionized gas may also play a role. On the other hand, the required minimum separation between the neutral islands would go down if a smaller suffices to allow transmission through the forest, or if the ionizing emissivity is in fact higher. However, the mean free path to ionizing photons has recently been measured at to be (Worseck et al., 2014), only somewhat larger than our assumed here. While there are still uncertainties, and while the measured mean free path scales steeply with redshift (), viable models are unlikely to have neutral islands spaced much more closely than this.

We can then use this requirement on to get some sense of which volume-averaged neutral fractions are plausible at . In the simulation outputs considered here (see §3), the mean separation between neutral islands is pMpc, 5.3 pMpc, and 2.7 pMpc for , 0.22, and 0.35, respectively. The first case certainly satisfies the requirement described above, the second case is just a bit on the small side, while the third case is uncomfortably small. Given the uncertainties in this argument, and the possibility that the neutral islands are a bit larger than in our simulation (which would increase their mean separation at fixed filling factor), we consider all three cases, but refrain from considering still more neutral models. We regard the latter case () as an extreme scenario intended mostly for illustration.

Finally, it is worth keeping in mind that any remaining neutral islands will likely be photoionized on a short timescale. For example, using Eq. 1 in Lidz & Malloy (2014) with , , and , the redshift interval over which the volume average ionized fraction transitions from to is only . However, it is possible that we are catching this – likely brief – phase in Ly- forest spectra and the possibility of testing this remains tantalizing.

## 3. Simulations and Mock Spectra

With the above discussion to frame the range of possibilities, we move to describe the numerical simulations used in this analysis and our approach to constructing mock Ly- forest absorption spectra before reionization completes. We use simulated density and ionization fields generated from a dark matter simulation of McQuinn et al. (2007a) which tracks dark matter particles in a simulation volume with a co-moving sidelength of . We assume that the gas closely follows the dark matter. In this work, we focus on redshift , but consider several possible neutral fractions. In practice, we obtain ionization fields with higher (lower) neutral fractions by using simulation outputs at higher (lower) redshifts. This should be an appropriate approximation since the statistical properties of ionized regions at a given neutral fraction are most sensitive to the neutral fraction and are relatively insensitive to the redshift at which the neutral fraction was attained (see McQuinn et al. 2007b and Furlanetto et al. 2004).

We generate mock quasar spectra according to the usual “fluctuating Gunn-Peterson” approach (e.g., Croft et al. 2002), with a few refinements to capture the main effects of incomplete reionization. First and foremost, we do not assume a fully ionized IGM. The transmission in the Ly- forest is sensitive to the precise ionized fractions in the ionized phase of the IGM. In order to simplify our study, as mentioned in the previous section, we rescale the simulated photoionization rates in the ionized regions to match the observed mean transmitted flux through the Ly- forest. We do this assuming ionization equilibrium, and a constant value of the UV background (with a photoionization rate per atom of ) within ionized regions. Specifically, simulated pixels with are considered highly ionized while less ionized pixels are considered fully neutral.222After effectively thresholding the ionization field in this way, we end up with neutral fractions which are higher than in the original simulation. Throughout the paper, we refer to increased, thresholded neutral fractions. This simplified approach allows us to consider a range of different possibilities for the ionization state of the IGM quickly. We comment on the shortcomings of this approach when appropriate. The optical depth of a given pixel, , in the simulation can then be found by summing over contributions from neighboring pixels (Bolton & Haehnelt 2008):

 τα(i) =cσαδRπ1/2∑jn%HI(j)b(j)H(a,x), (4)

where is the Doppler parameter, is the temperature of pixel , is the pixel proper width, is the Ly scattering cross section, is the proton mass, is the Hjerting function, and is the number density of hydrogen atoms at pixel , found using the simulated density field. To calculate the Doppler parameter, we assume that the gas obeys a modified temperature-density relationship

 T(δ)={T0(1+δ)γ−1if % ionized1,000Kif neutral, (5)

where is the matter overdensity in units of the cosmic mean and we choose K and as the temperature at mean density and slope of the temperature-density relation, respectively. For simplicity, we assume the ionized gas lies on the aforementioned temperature-density relation, although there should be significant scatter around this relation close to reionization (e.g. Lidz & Malloy 2014). We do not expect this to impact our conclusions significantly. The neutral gas should be colder than the ionized gas, of course, with a temperature set perhaps by low levels of X-ray pre-heating. Here we adopt K for the neutral gas; this choice is likely a bit large (it was chosen partly for ease in computing the Hjerting function below), but we have checked that we get nearly identical results for colder temperature choices.

The Hjerting function is a convolution of a Lorentzian profile, which incorporates the natural line profile of the Lyman-series lines, with a Maxwell-Boltzmann distribution, which accounts for the effects of thermal broadening on the line profile. The Hjerting function is defined by:

 H(a,x) =aπ∫∞−∞e−y2dya2+(x−y)2,

where , is the damping constant, is the Ly wavelength, is the relative velocity of pixel and pixel in units of the Doppler parameter, defined as , where . The peculiar velocity field is generated by applying linear perturbation theory to the underlying density field.333This was done because the full peculiar velocity field was not readily available, but this approximation should not impact our results. In detail, the natural line profile is only approximately described by a Lorentzian (Lee, 2013), with asymmetric corrections becoming important far from line center. In this study, the precise shape of the damping wing far from line center is unimportant: we make use only of the gradual recovery in transmission around saturated neutral regions, rather than the detailed shape of this recovery, which is also strongly influenced by neighboring neutral regions. We hence expect the Lorentzian form to be a good approximation for our present purposes.

In addition to including absorption from the hydrogen damping wing, we also include absorption from primordial deuterium. As a result of big bang nucleosynthesis, primordial hydrogen should be accompanied by traces of deuterium, with a relative abundance by number of (Cooke et al. 2013). Due to its slightly increased reduced mass, Lyman series transitions in deuterium will be shifted blueward by 82kms compared to the same transitions in hydrogen. We account for deuterium by scaling the number density of hydrogen in a given pixel by the relative abundance and shifting the resulting optical depths blueward by 82kms. Additionally, the Doppler parameter is adjusted to to account for the increase in mass.

In this work, we focus mostly on and adopt a mean transmitted flux at this redshift of , consistent with determinations from e.g., Becker et al. (2001). In some cases, we test the sensitivity of our results to the mean transmitted flux by considering as well. In general, the lower the mean transmitted flux, the more challenging it is for us to identify any remaining neutral islands. On the other hand, the likelihood that neutral islands remain increases towards high redshift and decreasing mean transmitted flux. As mentioned previously, we rescale the simulated photoionization rates in the ionized regions to a uniform value, normalized so that an ensemble of mock spectra matches the observed mean transmitted flux. It is important to note that the mean transmitted flux is a very steep function of redshift near , and that the sightline-to-sightline scatter in this quantity is substantial (Fan et al., 2006), and so one may want to carefully test for sensitivity to the precise redshift binning used.

We use the same approach as described above to generate Ly mock spectra, with , , and . However, in generating Ly mock spectra, we must also account for foreground Ly absorption due to gas at lower redshifts, , where is the redshift of the foreground Ly absorber and is the redshift of the Ly absorber. We will assume we are investigating quasar spectra at for this work, such that the corresponding foreground Ly absorption in the Ly spectra occurs at redshift . We adjust for the foreground Ly absorption to match measurements of the mean transmission from Becker et al. (2013) at these redshifts ( at ).444We generate foreground Ly absorption by considering the absorption from regions in the same simulation box, but demand that they are widely-separated from the high redshift regions of interest (). This enforces that the underlying density fields sourcing the Ly absorption and the foreground Ly absorption are uncorrelated, as should be the case for actual spectra. The optical depth of a pixel in a Ly spectrum is then the sum of the contribution from the foreground Ly absorption and the intrinsic Ly absorption .

In Fig. 1, we show an example mock Ly spectrum for a particular line of sight through the simulation. We show only a portion of the line of sight in order to exhibit smaller-scale features. The top figure shows the Ly transmission when the hydrogen damping wing is neglected (black) and when it is included (dashed red), while the bottom panel shows the underlying thresholded ionization field. We have neglected peculiar velocities in creating this figure in order to facilitate a comparison between the spectrum and the underlying ionization field.

From this figure, we see that the damping wing indeed has a significant effect on the transmission, but that its effect is hard to discern without knowing the damping-wing-less transmission. This is the case for two reasons. First, the forest here is very absorbed and the damping wing absorption becomes mixed with resonant absorption from neighboring ionized regions. Second, the damping wing from a particular neutral region may overlap with the damping wing from another neutral region, altering the shape of the resulting absorption. Specifically, we see that, in the example spectra, the region at is sandwiched between HI regions to the left and right, both within 1000kms. Therefore, the optical depths for the corresponding pixels likely have significant contributions from resonant absorption, damping wing absorption from the HI region to the left, and damping wing absorption from the HI region to the right. While detecting individual instances of damping wing absorption in this case seems impossible, we will show that detecting the presence of damping wing absorption on average should be feasible through the stacking of high-redshift quasar spectra.

## 4. Dark Gap Statistics

With the mock spectra of the previous section in hand, we now consider the size distribution of regions of saturated absorption – dark gaps – and its dependence on the underlying neutral fraction. Using such dark gap statistics has been widely discussed as a potential probe of the IGM neutral fraction (see, e.g., Fan et al. 2006, Mesinger 2010, McGreer et al. 2011). In a fully ionized IGM, the size of dark gaps in quasar spectra should grow with increasing redshift, simply owing to the increasing mean density of the universe and as a result of any decline in the intensity of the UV radiation background. However, once quasar spectra start to probe the tail end of reionization, the increase in dark gap size should accelerate due to the presence of islands of neutral hydrogen.

In Fig. 2, we have plotted the size distribution of dark gaps, , in blue for the mock spectra, assuming a mean transmission of . Additionally, the dashed curves display the two underlying populations of dark gaps: those sourced by ionized gas (magenta) and those sourced by neutral gas (cyan). For clarity, we only show dark gaps larger than (), since smaller saturated regions will be predominantly ionized. Additionally, we have neglected peculiar velocities when generating spectra here. Two important points become apparent from this figure. First, at (), dark gaps transition from being primarily sourced by ionized gas to being primarily sourced by neutral gas. This reinforces our intuition that, in a partially neutral IGM, the largest dark gaps should correspond to the remaining neutral islands. Second, the dark gaps being composed of two different populations gives rise to a bimodality in the size distribution. This suggests that the behavior of the large- tail of the size distribution may offer additional information about the neutral fraction, with a steep decline suggesting a highly ionized IGM and a more gradual decline, or the emergence of a second peak, suggesting a significantly neutral IGM. Such a “knee” in the dark gap size distribution is also mentioned in Mesinger (2010).

Additionally, we can consider the large- tail of the size distribution and its dependence on neutral fraction at a fixed mean transmission. In Fig. 3, we plot an expected histogram of dark gap sizes for 10 spectra for (magenta), 0.05 (cyan), 0.22 (blue), and 0.35 (black), again assuming that . Three trends become obvious from this plot. First, as the neutral fraction is increased (at fixed ), the number of large saturated regions increases and, second, as the neutral fraction is increased, the size of the largest dark gaps also increases. For example, in the model, the largest dark gaps are roughly five times bigger than in the fully ionized model. Additionally, we again see hints of the underlying dark gap size distribution being bimodal as the neutral fraction is increased, supporting the idea that the shape of the dark gap size distribution may be a diagnostic for the underlying neutral fraction.

Given these trends, it should be possible to compare dark gap distributions from observed spectra against models at various neutral fractions and use this to constrain the mean neutral fraction of the IGM. This approach is appealing in that it does not require especially high-resolution or high signal-to-noise spectra. However, it does require comparison with simulated models of the dark gap size distribution and so the conclusions reached will be somewhat model dependent. Additionally, the distributions are dependent on the assumed mean transmission, which is itself uncertain. In particular, estimates of the mean transmission at high redshift may be impacted by continuum fitting errors, given the inherent difficulty in estimating the unabsorbed continuum level in highly-absorbed spectra.

In order to investigate the impact of possible continuum fitting errors, we generate mock spectra in the fully ionized model with but then rescale the flux in each simulated pixel by a multiplicative factor – to mimic the effect of continuum misplacement – such that the measured mean transmitted flux appears to be . This case is shown as the magenta dashed line in Fig. 3. Here the dark gap size distribution is shifted towards sizes than one would expect in an ionized model at . However, the shape of the size distribution is still quite different than in the partly neutral models. Importantly, the dark gap distribution in the ionized model still lacks the distinctive bump at large sizes that is the hallmark of a partly neutral IGM in our models.

## 5. Stacking Toy Spectra

In this section we describe our basic approach of stacking Ly and Ly spectra in order to detect the presence of HI damping wings and absorption due to deuterium, respectively. While the forest is too absorbed at these redshifts to easily detect damping wings or deuterium absorption due to individual neutral regions, here we demonstrate that the presence of such features can be revealed on average by stacking regions of transmission over many spectra.

This section serves as a proof of principle by applying a simplified stacking approach to mock spectra generated using an idealized IGM model. Specifically, we consider an ensemble of sightlines through our simulation box and assume that the IGM is entirely ionized with the exception of a single HI island with mean density and varying length, , inserted randomly along each line of sight. We then generate mock spectra assuming these density and ionization fields. The stacking in this section is always done starting at the HI/HII boundaries of a given HI region moving outward.

### 5.1. HI Damping Wing

Our stacking approach can be clearly demonstrated by considering the damping wing from neutral hydrogen. Due to the natural width of the Ly line, a neutral hydrogen gas parcel should cause Ly absorption over a range of frequencies. Far from line center, this absorption will have an optical depth roughly following (Miralda-Escude & Rees 1997):

 τDWLyα(Δv)≈τGPRαcπ[1Δv−1Δv+v%ext] (6)

where is the Gunn-Paterson optical depth, , is the Ly decay constant, is the separation from the HI/HII boundary in velocity space, is the extent of the hydrogen region in velocity space, and is the speed of light. For a large neutral region, this equation implies that at . This excess absorption is referred to as the hydrogen “damping wing”. While both neutral gas and highly ionized gas can cause absorption in quasar spectra, only a significantly neutral hydrogen patch will result in damping wing absorption, owing to the greatly reduced optical depth in the wing compared to line center. As such, detecting damping wing absorption would be a smoking gun for the presence of significantly neutral hydrogen islands. Note that the transmission profile will differ from the simple form of Eq. 6, owing mostly to neighboring neutral regions, however the gradual recovery to transmission around saturated neutral regions should be a distinctive indicator that highly neutral regions remain in the IGM.

In Fig. 4 we show the results of stacking transmission outside of neutral regions in the toy mock spectra described earlier in this section, neglecting deuterium for the time being. Namely, we show the stacked transmission outside neutral islands of length () in black, () in blue, () in cyan, and stacked transmission neglecting the damping wing in red. Additionally, we have plotted the analytic curves corresponding to Eq. 6 for the various values, shown with dashed curves. We have applied a single multiplicative factor to these curves to account for average resonant absorption from ionized gas. Together, this figure implies that damping wing absorption from isolated neutral regions has a significant impact on quasar spectra, extending past the HI/HII boundaries, which may be observable through stacking as expected from Eq. 6.

In providing a toy example of how the hydrogen damping wing can affect spectra, we have neglected many important challenges that such a measurement would face. For example, we assumed perfect knowledge of the underlying ionization state of the IGM in order to determine where to stack and we assumed that we could discriminate between neutral and highly ionized absorption systems. However, the presence of such a large and potentially-observable feature provides motivation for us to apply the stacking approach in a more realistic manner. In §6 and in Appendix B, we describe several such challenges and subtleties along with potential resolutions.

### 5.2. Deuterium

With the stacking approach of the previous section in mind, we now consider absorption due to deuterium. As noted in §3, primordial hydrogen should be accompanied by traces of deuterium, with a relative abundance of (Cooke et al. 2013). Due to its slightly increased reduced mass, atomic transitions in deuterium are shifted blueward by 82kms compared to the same transitions in hydrogen. This implies that absorption due to neutral hydrogen in the IGM should be accompanied by additional absorption from deuterium, shifted blueward by 82kms. We can estimate the optical depth for Ly absorption in deuterium at cosmic mean density by simply scaling the hydrogen Ly optical depth by the deuterium abundance:

 τD,α =[DH]×τGP≈8.25xHI(1+δ)[1+z6.5]3/2. (7)

Thus, we see that while the relative abundance of deuterium is extremely small, the Gunn-Peterson optical depth is so large that the resulting deuterium optical depth is still of order 10 in Ly.

An appealing aspect of searching for damping wing absorption is that the optical depth in the wing is large enough to cause significant absorption in the presence of neutral islands, but small enough to be negligible for ionized absorption systems. We see this again in the case of deuterium absorption, suggesting that it may be useful as an additional “smoking gun” indicator for underlying neutral hydrogen. However, an obvious problem with detecting deuterium in Ly spectra is that the feature should be narrow and well within the broad range of velocities where the hydrogen damping wing is significant. Specifically, according to Eq. 6, at , the damping wing optical depth for an extended neutral region should be . Therefore, the feature should be completely wiped out in Ly spectra by the hydrogen damping wing.

However, the damping wing optical depth in the Ly line is much smaller. Specifically, according to Eq. 6, the damping wing optical depth scales as

 τDWLyβτDWLyα =τGP,βτGP,α×RβRα=fβλβfαλα×(Γβ+ΓHα)λβΓαλα =f2βf2α⎛⎝1+fHαfβλ2βλ2Hα⎞⎠=.0410, (8)

and should therefore be significantly narrower in Ly than in Ly. In the above expression, , , and are the oscillator strengths of the Ly, Ly, and Balmer- transitions, respectively, with denoting the corresponding wavelengths and denoting the corresponding decay constants. By modifying Eq. 6 for Ly, we see that . Therefore we find that the hydrogen damping wing should not wipe out deuterium absorption features in Ly. Furthermore, while the hydrogen damping wing optical depth is reduced by a factor of roughly when considering Ly, the total optical depth in the deuterium line is only reduced relative to deuterium Ly by , such that the optical depth should still be of order 1 for deuterium Ly. Therefore, not only should a deuterium absorption feature survive the hydrogen damping wing, but it should still have a strong enough optical depth to cause significant absorption if neutral islands in fact remain.

Naturally, it should be very difficult to detect individual deuterium absorption features from the diffuse IGM, as the Ly- spectra will be very absorbed when the universe is neutral enough to produce the features in the first place. However, the feature may nonetheless be observable on average through the stacking of high-resolution quasar spectra. In order to demonstrate the strength of the deuterium absorption feature in stacked spectra, we incorporate deuterium into the same toy sightlines from §5.1 to produce mock Ly spectra, neglecting foreground Ly absorption for the time being. We are then able to stack transmission outside of neutral regions in the spectra, starting at the HI/HII boundaries and moving outward. However, since deuterium absorption will only occur on the blue side of neutral regions, we need only stack those regions of transmission. In fact, this offers a clean test for detecting deuterium. Namely, we can separately stack transmission redward and blueward of neutral regions and compare. Excess absorption on the blue side of neutral regions, on average, could signal the presence of deuterium absorption. This is especially appealing since there should be no sources of contamination that would cause a similar, and significant, red/blue asymmetry.555One source of asymmetry we do find, which can be seen in Fig. 9, results from the fact that, when dealing with realistic spectra, we force there to be transmission in Ly at the locations where stacking begins. This results in a small selection effect, where selected neutral absorption systems have a reduced probability of having nearby neutral regions and have correspondingly-smaller nearby optical depths. For deuterium, this smaller optical depth is shifted blueward, causing less absorption on the blue side of the line for . However, this asymmetry is minor and opposes the asymmetry from deuterium absorption.

In Fig. 5, we show the results of stacking transmission in Ly redward (red) and blueward (black) of the toy neutral regions across the full ensemble of mock quasar spectra. As in Fig. 4, all stacking begins at HI/HII boundaries. We can see the blueward transmission clearly exhibits excess absorption due to deuterium extending roughly 80kms from the HI/HII boundary. Thus, in this idealized scenario, the presence of deuterium in islands of neutral hydrogen leaves a very clear signature in the stacked Ly transmission.

Before proceeding further, we should point out one important caveat here. In our simulated models, the transition between fully neutral and highly ionized regions is, by construction, perfectly sharp. If this transition is more gradual in reality, then the narrow deuterium feature could be overwhelmed by absorption from mostly ionized hydrogen in this transition region. A minimal scale for this transition region is set roughly by the mean free path to ionizing photons through the neutral IGM, which is only . This minimal scale is two orders of magnitude smaller than the scale of the deuterium feature and hence does not present a worry. However, if the edges of the ionized regions tend to experience a reduced ionizing background, this might obscure the deuterium feature, even in the case of a partly neutral IGM. We believe the possibility of detecting this deuterium feature is enticing enough to warrant further investigation.

As was the case in §5.1, we have made several simplifying assumptions and have additionally neglected foreground Ly absorption from the lower-redshift IGM. However, the clear presence of deuterium absorption revealed through the simplified stacking approach provides motivation to also consider applying it to more realistic spectra, as will be discussed in §6.

## 6. Steps of Approach

In §5, we demonstrated the utility of stacking idealized quasar spectra in order to reveal the presence of the HI damping wing and deuterium absorption. The success of this approach in the toy case provides motivation for us to apply it to realistic mock spectra. In doing so, we must confront the simplifying assumptions made in §5.

The most obviously unrealistic assumption made in §5 is that we can precisely identify the HI/HII boundaries underlying our spectra. In practice, we will only have access to the level of transmission at each point along the spectra. However, based on Fig. 5, the recovery from saturated absorption to transmission occurs within in Ly from the edge of the neutral zone, and should therefore provide a relatively good indicator of the HI/HII boundary. Therefore, we choose to identify stacking locations based on where transmission recovers in Ly. To be clear, for the case of the hydrogen damping wing, we are stacking transmission in the Ly forest, but we are choosing where to start the stacking based on features in the Ly forest. A drawback of this approach, when searching for the hydrogen damping wing, is that we are only able to stack regions of the Ly forest with corresponding regions in the Ly forest that are not contaminated by Ly absorption. This effectively reduces the amount of usable spectra, since, for a quasar at , the pure Ly forest will extend , but Ly absorption will contaminate the Ly forest at . If presented with a limited number of spectra, it may be worth searching for the damping wing by using only the Ly regions of the spectra.

By stacking at the precise locations of HI/HII boundaries in §5, we were also ensuring that our sample of absorption systems was all neutral. However, when we modify our approach to begin stacking at locations where transmission recovers from saturated absorption, we may start stacking transmission outside of ionized absorption systems together with transmission outside of neutral absorption systems, diluting our signal. Since the signal we are aiming to find is small to begin with, it is important that we minimize this contamination from ionized regions. To do this, we take advantage of the main argument of §4, namely that regions of saturated absorption sourced by neutral gas should be significantly larger, on average, than those sourced by ionized gas. Therefore, we choose to stack only transmission outside of large regions of saturated absorption. Furthermore, since true neutral regions should cause saturated absorption in Ly, we choose to stack only outside of large saturated regions which are fully absorbed in Ly, where we define “large” to be () in Ly. Note that this choice is tuned for the case of : a different choice may be better for other values of the mean transmitted flux. At any rate, in applying these tests to real data, one would likely vary this size scale across a range of possible values.

Additionally, an appealing feature of the search for deuterium absorption is that it offers a very clean test for its detection, namely a red/blue asymmetry in the transmission outside of plausibly neutral regions. The disparity in the size distribution of saturated regions sourced by neutral and ionized gas suggests a similar test may be possible for the detection of the HI damping wing. Namely, while large regions of saturated absorption are likely to be sourced by neutral gas, small regions of saturated absorption are likely to be sourced by ionized gas. Therefore, to find evidence of excess absorption outside of neutral regions due to the HI damping wing, we compare the stacked transmission outside of large absorption systems, plausibly sourced by neutral gas, to that outside of small absorption systems, likely sourced by ionized gas. A significant amount of excess absorption outside of the former compared to the latter, extending further than any possible density correlations, would suggest the presence of damping wing absorption.

Furthermore, in §5.2, we discussed how the damping wing is greatly weakened in Ly compared to in Ly. Therefore, an additional test for the presence of damping wing absorption could be to take the ratio of the stacked Ly transmission to the stacked Ly transmission, where stacking occurs in the same physical regions in both cases. In the event that there is significant damping wing absorption, this ratio should also slowly recover to some constant value at large . We further discuss and develop this approach in Appendix B.

When dealing with realistic spectra, we must adjust our approach to accommodate the presence of noise (and finite spectral resolution). While noise should average out in stacked regions, the presence of noise will also obfuscate the precise boundaries between saturated absorption and transmission. We choose to handle this by smoothing our noisy spectra over a scale of 100 km/s () and defining any pixel, , with transmission to be consistent with saturated absorption, where denotes the standard deviation of the smoothed noise. We then define regions in the smoothed spectra where the flux goes from to as the transitions from saturated absorption to transmission, and therefore as potential points to start stacking. When stacking transmission, however, we stack the transmission in the unsmoothed spectra.

Another concern is that damping wing absorption sourced by DLAs may erroneously be attributed to a significantly neutral IGM. However, in Appendix A, we estimate the expected rate of DLAs occurring in quasar spectra and find it is small enough to be ignored. Additionally, DLAs may be discriminated from diffuse neutral islands based on the presence of metal lines and the relative sizes of their absorption in Ly and Ly.

Finally, as mentioned previously, we approximate the ionizing background in the ionized regions as uniform and ignore scatter in the temperature density relation. Accounting for these fluctuations might lead to a more gradual recovery in the transmission around absorbed regions – in the case of a fully ionized universe – than in our models. Further investigation of this issue would certainly be required if a gradual recovery is indeed found in real spectra. In Appendix B, we discuss a possible empirical test that may help in this regard.

## 7. Results

Having considered the subtleties of the previous section, we are now ready to apply the three-pronged approach to more realistic mock spectra. In each section, we first consider the ideal case where no noise has been applied to give an idea of the potential constraining power of the different methods. Subsequently, we add realistic levels of noise and consider realistic spectra resolution to give an idea of the constraining power of the approaches applied to Keck HIRES spectra for the deuterium feature and damping wing, and spectra with slightly higher resolution than SDSS for the dark gap size distribution.

### 7.1. Detecting the Damping Wing

We first consider the ability to uncover the presence of the hydrogen damping wing by strategically stacking regions of transmission in the the Ly forest of noiseless mock quasar spectra. As discussed in §6, our aim is to compare the average transmission outside of plausibly neutral absorption systems to the transmission outside of likely ionized systems.

We identify the plausibly neutral absorption systems by requiring the regions be completely absorbed in Ly, and also that the regions of saturated absorption are at least () in Ly. We begin stacking at the point in the Ly spectrum which corresponds to the recovery from absorption to transmission in Ly. We identify the likely ionized absorption systems by requiring that they are below a maximum length in Ly.

In Fig. 6, we show the results of applying this approach to realistic mock spectra generated assuming various ionization states of the IGM. In the top panel, we show the stacked transmission outside of plausibly neutral absorption systems (solid) and likely ionized absorption systems (dashed), using a volume-averaged neutral fraction of (black), 0.22 (blue), 0.05 (cyan), and (magenta). The curves agree with our expectations, namely that transmission outside of neutral regions should recover more slowly and exhibit a rough damping wing shape with a large extent in velocity space. We see that the excess absorption extends farther than the expected from an isolated damping wing. However, as discussed in Appendix C, we find that the spatial clustering of neutral regions is responsible for this effect.

While, for several reasons discussed earlier, the shape of the absorption is distorted compared to Fig. 4, it can be seen for all significantly neutral ionization states. An important check is to apply the stacking procedure to a fully ionized IGM and ensure that we do not make a false detection. The results of this check are shown by the magenta curves in Fig. 6. As we can see, the resulting stacked transmission outside of plausibly neutral regions lacks an overall damping wing shape and stays roughly fixed near the mean transmission.

We can also see that the transmission outside of small absorption systems is very sensitive to the underlying neutral fraction. We expect this, however, since this stacked transmission depends strongly on the average transmission in regions which are not in saturated absorption, denoted . Since the dark pixel covering fraction in our mock spectra increases with , mock spectra with larger neutral fractions must have larger values for to maintain . As such, Fig. 6 shows that the stacked transmission outside of small absorption systems increases monotonically with .

We estimated the stacked transmission from a large ensemble of simulated spectra to produce a smooth estimate of the average transmission around saturated regions in each model. The transmission curves outside of individual absorption systems are, however, quite noisy on their own such that, from saturated region to saturated region, there is significant scatter about the mean-value curves shown in the top panel. In order to estimate how confidently we can distinguish the solid and dashed curves with a reasonable number of quasar spectra, we scale the number of identified absorption systems to what we would expect using spectra. Specifically, we take the difference between the dashed and solid curves and divide by the scatter in each bin. The scatter of each bin is simply the scatter in stacked transmission outside of large absorption systems, scaled by , added in quadrature with the scatter in the stacked transmission outside of small absorption systems, scaled by . Here we scale to estimate the plausible scatter around the mean after estimating the transmission around saturated regions using quasar absorption spectra.

The results of this are shown in the bottom panel of Fig. 6 for the same ionization states. The results appear to be very encouraging, indicating that, assuming noiseless spectra, the solid and dashed curves are statistically-significantly different (even for !). In addition, we see that the difference roughly follows a damping wing shape and remains significant for . We should emphasize that, while the deuterium absorption feature will necessarily be a feature and require high resolution spectra to be seen, the damping wing feature extends an order of magnitude farther in velocity space and should be accessible to lower-resolution spectra.

In Fig. 7, we show the same results as in Fig. 6, but assume a lower mean transmission of , consistent with spectra at (Becker et al. 2001). From the figure, we see that these results are very similar to those for , but with the significance curves peaking at a lower value and with the stacked transmission recovering to a lower mean. Overall, this provides encouragement for applying the approach to higher- spectra, suggesting that a range of physically interesting neutral fractions could be probed.

It is also interesting to consider these results when spectra are generated according to the specifications of existing data. In Fig. 8 we show the same results as in the bottom panel of Fig. 6 except we have adjusted the spectra to mimic HIRES spectra. Namely, we have assumed a spectral resolution with and bins with size (e.g. Viel et al. 2013). Additionally, we have assumed a signal to noise of at the continuum per pixel and that we have 20 such spectra. While we are currently only aware of 10 such spectra, this case is still interesting since spectra with significantly worse spectral resolution should also be adequate for this test.

From this figure, we can see that, despite the degradation of the spectra, the damping wing is still visible with the significance curve peaking at () significance for the () ionization state. However, this figure suggests that, in the case, it is less-clear whether the damping wing is detectable.

An important effect of adding noise to the mock spectra is that it obscures the precise location where spectra should be stacked and also increases the fraction of selected saturated regions which are, in fact, ionized. We find that for the spectra in this section , 40%, and 75% of identified plausibly neutral regions are in in fact ionized for , 0.22, and 0.05, respectively. This is compared to , 10%, and 20% contamination when noise is neglected.

Statistical significances in this section are only estimates. In reality, the statistical significance with which the damping wing can be detected will depend on how extended the significance curves are, along with how correlated the errors in neighboring bins are. We discuss this in §8.

### 7.2. Deuterium Feature Results

We now turn to consider the prospects for identifying deuterium absorption in realistic Ly mock spectra. As discussed in §6, our aim is to identify plausibly neutral absorption systems in the Ly spectra and compare the stacked transmission moving blueward and redward away from the absorption. We identify the plausibly neutral regions in the same manner as for the damping wing.

In Fig. 9, we show the results of applying the stacking approach to the same mock spectra as in the previous section. The top panel shows the mean transmission blueward (solid) and redward (dashed) moving away from plausibly neutral absorption systems for the same ionization states as in the previous section. We can very clearly see excess absorption in the partially neutral spectra, extending , consistent with our expectations from Fig. 5. Additionally, we also find that, in the fully ionized case, the blueward and redward stacked transmission match up very well.

As in the previous section, we can construct a rough significance curve for the difference between the blueward and redward transmission. Specifically, in the bottom panel of Fig. 9 we show the excess blueward absorption in units of the standard deviation of the stacked blueward transmission assuming 20 quasar spectra. We can see that the significance of the red/blue asymmetry extends () and is for all of the partially neutral models considered, with increasing significance for models with higher neutral fractions. Additionally, we see that the curve corresponding to the fully ionized model shows no statistically-significant deviation from red/blue symmetry. Thus, this is indeed a very clean test for the presence of deuterium. However, the signal itself is an order of magnitude smaller in velocity-space extent and is found with significantly less statistical significance than the damping wing signal. Therefore, we expect that high-resolution, high-signal-to-noise spectra will be necessary to search for it.

As in §7.1, we can reproduce Fig. 9 assuming spectra with specifications mimicking Keck HIRES. Unfortunately, we find that, with a signal to noise per pixel in the continuum of 10, the deuterium feature is hard to observe. Because of this, we consider using 20 HIRES-style spectra with a signal to noise per pixel of 30 in the continuum. While this signal-to-noise value is higher than those for existing spectra we found in the literature, it is not unreasonable to assume such spectra may become available in the future. Furthermore, this may provide additional motivation to obtain such spectra. Regardless, after applying the stacking approach with twenty HIRES spectra, we obtain the results shown in Fig. 10. This figure shows that the feature should be observable with modest statistical significance. Specifically, for (0.22) the significance curve peaks at (). Additionally, when these curves are generated assuming MIKE-style spectra, with spectral resolution and velocity bin size , we obtain similar curves as in Fig. 10 but with the signal being statistically significant over a smaller range of velocities.

Again, important effects of adding noise to the mock spectra are that it obscures the precise location where spectra should be stacked and increases the fraction of selected plausibly neutral regions which are, in fact, ionized. We find that for the spectra in this section , 20%, and 40% of identified plausibly neutral regions are in in fact ionized for , 0.22, and 0.05, respectively. This is compared to , 10%, and 20% contamination when noise is neglected. As expected, we find a smaller level of contamination than in the previous section, owing to the increased signal to noise of the spectra used. However, for the case of deuterium, the effect of noise on the stacking location is more apparent. Fig. 9 demonstrates that, without noise, deuterium absorption imprints a feature on stacked noiseless spectra extending , but only extending in stacked noisy spectra, as shown in Fig. 10.

The above results suggest that stacking Ly transmission in high- spectra can indeed be used to detect the presence of primordial deuterium, and hence that of hydrogen, but that high-resolution and high signal-to-noise spectra will be required. Nevertheless, it would certainly be interesting to apply this approach to existing HIRES/MIKE spectra as it provides an additional test, independent of the damping wing search, for the presence of underlying neutral hydrogen in the IGM. As such, a detection with modest levels of statistical significance could lend credence to a claimed detection of the HI damping wing.

### 7.3. Dark Gap Statistics

We now shift our focus away from stacking and toward the distribution in lengths of regions of saturated absorption (dark gaps). As discussed in §4, the dark gap size distribution in quasar spectra should contain information about the underlying ionization state of the IGM. Specifically, in a more neutral IGM, the typical sizes of dark gaps should be larger and the shape of the dark gap size distribution should have a more gradual decline, and possibly show a hint of bimodality, toward large .

We continue this discussion in this section by considering plausible dark gap size distributions that could be observed with moderate-resolution, moderate-signal-to-noise spectra. Specifically, we consider spectra with spectral resolution , bin size , and a signal-to-noise ratio of 10 at the continuum. These spectra are of only slightly better quality than SDSS spectra. Additionally, since we are not concerned with Ly transmission, we are able to use the entire Ly forest for each spectra.

In Fig. 11, we show the resulting dark gap size histogram expected for 20 such spectra for (black), 0.22 (blue), 0.05 (cyan), and 0 (magenta). In generating this figure, we use the same ensemble of mock spectra as in §7.1 and §7.2, except with their spectral resolution and bin size modified as mentioned. We maintain the requirement that .

This figure qualitatively agrees with Fig. 3, where noiseless spectra with finer spectral resolution were used, but shows a shift toward larger due to smoothing the spectra. Additionally, the ionization states are not as easily distinguishable as in Fig. 3, with the distribution looking practically identical to the fully ionized scenario. However, for the other neutral fractions considered, the situation looks very encouraging. The distributions for and 0.35 show a more gradual decline toward large than the fully ionized case and also reveal the clear emergence of a bimodal distribution. Additionally, the largest dark gaps in these ionization states are roughly twice as large as in the fully ionized case.

## 8. Forecasts

Having discussed the results of the proposed stacking approaches applied to realistic mock spectra, we now consider the ability of these methods to constrain the ionization state of the IGM. Specifically, in this section we focus on the ability of each method to rule out the hypothesis of a fully ionized IGM.

In both the case of the HI damping wing and deuterium absorption feature, we would like to compare models representing different ionization states and estimate the between models and the fully ionized model, assuming a reasonable number of spectra. Let denote the mean behavior for a model with given neutral fraction, , as a function of stacked velocity separation and let denote the mean behavior of the ionized model, also as a function of stacked velocity separation. The precise definitions of what is meant by behavior will be discussed later in this section. In this case, the value between two models can be calculated by

 Δχ2⟨xHI⟩ =ΔF⟨xHI⟩⋅C−1⋅ΔFT⟨xHI⟩ (9)

where is the covariance matrix of the model, representing the correlation between stacked pixels, and , with denoting its transpose. For simplicity, rather than estimating the full covariance matrix and its inverse, we approximate pixels at sufficiently wide separations as independent. We then coarsely sample the pixels – on the scale where they can be well approximated as independent – and assume a diagonal covariance matrix for the coarsely sampled pixels. Specifically, we estimate by simply adding up the squared statistical significance of each coarsely-sampled bin, , where is the standard deviation of .

### 8.1. Deuterium

Perhaps it is best to consider the case of the deuterium absorption feature first. In the case of a fully ionized IGM, the transmission looking blueward and redward from absorption systems should be symmetric on average, with excess blueward absorption only occurring when the IGM is significantly neutral. Therefore, we may calculate the between stacked transmission looking redward () and blueward () from plausibly neutral absorption systems for each ionization state to estimate our ability to rule out the hypothesis of a fully ionized IGM in each case.

Thus, in the context of Eq. 9, we have

 ΔF⟨xHI⟩(Δv) ≡F⟨xHI⟩,red(Δv)−F⟨xHI⟩,blue(Δv) (10) C−1ij =δij/σ2⟨xHI⟩,blue(vi) (11)

where is the standard deviation of the stacked transmission blueward of plausibly neutral absorption systems, assuming a given number of spectra, and we have assumed that we have already resampled at sufficient velocity separations such that neighboring bins can be approximated as independent. At this point, the only missing ingredient is the minimum separation between two stacked pixels for them to be considered independent. We calculate the correlation function between stacked pixels in Ly, and find that the correlation has a width of and, as such, we do not expect to get more than one independent bin within the scale of the deuterium feature. Therefore, a rough estimate of the value obtainable in each ionization state can be estimated simply by the peak value in the “significance curves” in Fig. 10. Thus, if the underlying neutral fraction of the IGM is (0.35), then we expect to be able to rule out a fully ionized IGM with () confidence, assuming 20 HIRES/MIKE spectra with signal to noise of 30 per pixel at the continuum. Unfortunately, we do not expect to be able to rule out the hypothesis of a fully ionized IGM if the underlying neutral fraction is .

### 8.2. HI Damping Wing

Assessing the statistical significance with which we can observe the HI damping wing is slightly more complicated than the deuterium case since the test for its detection is not as simple. When faced with actual spectra, we would look for the presence of significant and extended absorption outside of large absorption systems compared to that outside of small absorption systems.

Therefore, the behavior we would like to compare in each case is the stacked transmission outside of plausibly neutral absorption systems () and the transmission outside of likely ionized absorption systems (). Let us denote

 F(Δv) ≡fsmall(Δv)−flarge(Δv) (12)

as the difference in these stacked transmissions where and represent this behavior for the ionization state with neutral fraction and the fully ionized state, respectively. Thus, in the context of Eq. 9, we have

 ΔF⟨xHI⟩(Δv) =Fion(Δv)−F⟨xHI⟩(Δv) (13) C−1ij =δij/σ2F⟨xHI⟩(vi) (14)

where denotes the standard deviation of at . The resulting value indicates the expected significance with which we could rule out a fully ionized IGM if the neutral fraction were, in fact, . Again, for Eq. 14, we have assumed that has been resampled at velocity separations such that the pixels can be treated as independent. We find that the correlation function between pixels of stacked transmission in the Ly forest within the scale of the HI damping wing has . While this scale is large, the excess absorption due to the presence of damping wing absorption leaves a feature extending , leaving us with independent bins within the scale of the feature.

In this manner, we are able to calculate a rough estimate for the values for the ionization states considered thus far. Assuming the same type of spectra as in Fig. 6, namely 20 HIRES spectra with signal to noise in the continuum of 10 per pixel, we find that if the IGM is, in fact, , , or neutral, then we should be able to rule out a fully ionized IGM at , , or , respectively. In the case of , this reduces to , , and , respectively.666This is for the case where we do not attempt to further optimize the analysis for the decrease in transmission. It is possible that further gains could be made, with Fig. 7 representing a best-possible-case scenario. While we are only aware of such spectra that exist at the moment, we still regard this estimate as somewhat conservative. We found that excess stacked absorption due to the damping wing extends thousands of , and therefore it is not necessary to have the state-of-the-art in spectral resolution to measure it. Especially with such extended correlation among neighboring pixels, it is unclear how much is gained by resolution improvements beyond .

## 9. Conclusion

In this work, we developed empirical tests of the possibility that the Epoch of Reionization is not yet complete by . Specifically, we proposed three measurements that can be made with existing, and future, high-redshift quasar spectra in order to investigate this region of reionization history parameter space.

First, we discussed using the dark gap size distribution in quasar spectra as a means of constraining the neutral fraction. We find that not only do the typical sizes of dark gaps increase with neutral fraction but that the shape of the size distribution is also sensitive to the neutral fraction. Specifically, the presence of dark gaps sourced by significantly neutral hydrogen islands introduces a bimodality in the dark gap size distribution. We find that this bimodality should be observable at , provided that , and should not be affected by continuum fitting errors.

Next, we proposed a method for searching for hydrogen damping wing absorption by strategically stacking regions of transmission in the Ly forest. Specifically, we searched for excess absorption in stacked transmission outside of plausibly neutral regions compared to that outside of likely ionized regions. We found that the presence of the hydrogen damping wing will result in excess absorption extending past ionization boundaries of neutral regions. Furthermore, this excess absorption should be detectable with statistical significance for , using 20 HIRES-style spectra with a signal-to-noise value per pixel of 10 at the continuum.

Lastly, we proposed a similar stacking measurement utilizing the Ly forest in order to search for deuterium absorption associated with significantly neutral hydrogen islands at . We proposed searching for this feature by looking for excess absorption in stacked Ly transmission blueward of plausibly neutral regions compared to the corresponding redward transmission. We find that this feature should be observable in principle but will likely require additional high-resolution spectra in order to be detected. Specifically, we find that the feature should be observable at with () confidence using 20 HIRES-style spectra with a signal-to-noise value per pixel of 30 at the continuum if (0.35). While we are not aware of this many available spectra with such specifications, this provides motivation for acquiring such spectra in the future, possibly through the follow-up observation of SDSS quasars.

While we have taken many steps to ensure that the analysis of mock spectra presented in this work is realistic, there are still additional complexities that will be faced when one is presented with actual spectra. For example, we treat all portions of our spectra as being at when, in reality, the redshift will evolve along the lines of sight. In addition, we ignored spatial fluctuations in the UV radiation field and in the temperature density relation. Additional work will certainly be required to definitively interpret future measurements along the lines we suggest here. However, we believe the signatures explored here are well-worth further investigation and should ultimately improve our understanding of the reionization history of the IGM.

## Acknowledgements

We would like to thank Matt McQuinn for providing us with the simulations used in this analysis. We are also grateful to Andrei Mesinger and Jordi Miralda-Escudé for helpful questions and discussions during the ICTP “Cosmology with Baryons at High Redshift” workshop. The authors acknowledge support from the NSF through grant AST-1109156 and NASA through grant NNX12AC97G.

## References

• Becker & Bolton (2013) Becker, G. D., & Bolton, J. S. 2013
• Becker et al. (2014) Becker, G. D., Bolton, J. S., Madau, P., et al. 2014
• Becker et al. (2013) Becker, G. D., Hewett, P. C., Worseck, G., & Prochaska, J. X. 2013, MNRAS, 430, 2067
• Becker et al. (2001) Becker, R. H., et al. 2001, Astron.J., 122, 2850
• Bolton & Haehnelt (2007) Bolton, J. S., & Haehnelt, M. G. 2007, MNRAS, 382, 325
• Bolton & Haehnelt (2013) —. 2013, MNRAS, 429, 1695
• Cooke et al. (2013) Cooke, R., Pettini, M., Jorgenson, R. A., Murphy, M. T., & Steidel, C. C. 2013
• Cooray & Sheth (2002) Cooray, A., & Sheth, R. K. 2002, Phys.Rept., 372, 1
• Croft et al. (2002) Croft, R. A., Weinberg, D. H., Bolte, M., et al. 2002, Astrophys.J., 581, 20
• Djorgovski et al. (2001) Djorgovski, S., Castro, S., Stern, D., & Mahabal, A. 2001
• Fan et al. (2001) Fan, X., et al. 2001, Astron.J., 122, 2833
• Fan et al. (2006) Fan, X.-H., Strauss, M. A., Becker, R. H., et al. 2006, Astron.J., 132, 117
• Furlanetto et al. (2004) Furlanetto, S., Zaldarriaga, M., & Hernquist, L. 2004, Astrophys.J., 613, 1
• Gunn & Peterson (1965) Gunn, J. E., & Peterson, B. A. 1965, ApJ, 142, 1633
• Kuhlen & Faucher-Giguère (2012) Kuhlen, M., & Faucher-Giguère, C.-A. 2012, MNRAS, 423, 862
• Lee (2013) Lee, H.-W. 2013, Astrophys.J., 772, 123
• Lidz & Malloy (2014) Lidz, A., & Malloy, M. 2014, Astrophys.J., 788, 175
• Lidz et al. (2007) Lidz, A., McQuinn, M., & Zaldarriaga, M. 2007, Astrophys.J., 670, 39
• McGreer et al. (2011) McGreer, I. D., Mesinger, A., & Fan, X. 2011
• McQuinn et al. (2007a) McQuinn, M., Hernquist, L., Zaldarriaga, M., & Dutta, S. 2007a, Mon.Not.Roy.Astron.Soc., 381, 75
• McQuinn et al. (2007b) McQuinn, M., Lidz, A., Zahn, O., et al. 2007b, Mon.Not.Roy.Astron.Soc., 377, 1043
• Mesinger (2009) Mesinger, A. 2009
• Mesinger (2010) Mesinger, A. 2010, MNRAS, 407, 1328
• Mesinger & Haiman (2007) Mesinger, A., & Haiman, Z. 2007, Astrophys.J., 660, 923
• Miralda-Escude (1998) Miralda-Escude, J. 1998, Astrophys.J., 501, 15
• Miralda-Escude & Rees (1997) Miralda-Escude, J., & Rees, M. J. 1997
• Pentericci et al. (2014) Pentericci, L., Vanzella, E., Fontana, A., et al. 2014
• Prochaska et al. (2005) Prochaska, J. X., Herbert-Fort, S., & Wolfe, A. M. 2005, Astrophys.J., 635, 123
• Rahmati & Schaye (2013) Rahmati, A., & Schaye, J. 2013
• Robertson et al. (2013) Robertson, B. E., Furlanetto, S. R., Schneider, E., et al. 2013, Astrophys.J., 768, 71
• Schenker et al. (2012) Schenker, M. A., Stark, D. P., Ellis, R. S., et al. 2012, ApJ, 744, 179
• Schroeder et al. (2013) Schroeder, J., Mesinger, A., & Haiman, Z. 2013, MNRAS, 428, 3058
• Taylor & Lidz (2013) Taylor, J., & Lidz, A. 2013
• Viel et al. (2013) Viel, M., Becker, G. D., Bolton, J. S., & Haehnelt, M. G. 2013, Phys.Rev., D88, 043502
• Wolfe et al. (2005) Wolfe, A. M., Gawiser, E., & Prochaska, J. X. 2005, Ann.Rev.Astron.Astrophys., 43, 861
• Worseck et al. (2014) Worseck, G., Prochaska, J. X., O’Meara, J. M., et al. 2014
• Xu et al. (2014) Xu, Y., Yue, B., Su, M., Fan, Z., & Chen, X. 2014, Astrophys.J., 781, 97
• Zuo & Phinney (1993) Zuo, L., & Phinney, E. S. 1993, ApJ, 418, 28

## Appendix A: Contamination from DLAs?

A potential concern is that damping wings from super Lyman-limit systems and damped Ly- absorbers (DLAs) might produce “false positives” and contaminate our search for diffuse neutral islands. Since DLAs and super Lyman-limit systems are mostly associated with galaxies and the circumgalactic medium (see Wolfe et al. 2005 for a review), we would like to distinguish these absorbers from the more diffuse and spatially extended islands of neutral hydrogen that are the subject of our search. In addition, note that it is difficult to fully resolve and model high column density absorbers in cosmological simulations (e.g. Rahmati & Schaye 2013 and references therein) – especially given our present aim of capturing the large-scale features of reionization – and so the impact of these systems is not captured in our present modelling.

Fortunately, we don’t expect these dense absorbers to be a big contaminant, provided we make use of the Ly- forest – which helps owing to the lower cross section in the wing of the line (compared to Ly-) – and confine our search to fairly extended neutral islands. The Ly- line profile for a high column density absorber can be approximated by a Lorentzian, so that the optical depth at velocity offset is:

 τβ,DLA(Δv)=NHIσβ,0πRβ(Δv/c)2+R2β. (1)

For comparison, a fully neutral and isolated absorber of co-moving length produces saturated Ly- absorption over a velocity extent slightly larger than . We can then determine the column density required for a DLA to produce as long a saturated region in the Ly- forest as produced by a neutral island of length . We consider a DLA to produce saturated absorption at velocity separations where . Provided the extent of the absorber is large enough that (which is a good approximation for the extended neutral islands of interest), this critical column density, , is given by:

 NHI,crit=7.2×1021cm2[τβ,DLA3][1+z6.5][Lneut3.8Mpc/h]2.

The fiducial value of in the above equation corresponds to km/s – this is the minimum saturated stretch included in our stacks when we search for neutral regions (see §6). The column density required for a DLA to produce this much saturated absorption is quite large, and the abundance of DLAs with column densities larger than is very small.

Quantitatively, taking the Gamma function fit to the column density distribution of DLAs from Prochaska et al. (2005) 777Specifically, we use their highest redshift bin fit, which includes DLAs between redshifts . (which accounts for the sharp cutoff in the observed abundance of DLAs at high column densities), we find that the number of DLAs with is only . For reference, the redshift extent of the forest between the Ly- and Ly- emission line at these redshifts is roughly , and so such high column density DLAs should be exceedingly rare. Since is only a little smaller than the exponential cut-off in the column density distribution function, the results are rather sensitive to the precise choice of . Given that smaller column-density DLAs might still leak into our stack if they happen to be next to saturated ionized regions, it is worth checking this dependence. However, even choosing cm yields only , which is still smaller than the abundance of neutral islands we seek to detect. From these estimates, we expect very minimal contamination from DLAs leaking into our stacked sample of possible neutral regions. Note also that deuterium Ly- absorption from these high column density absorbers will be in the saturated part of the HI Ly- line, and so DLAs should not contaminate our search for the deuterium signature of neutral islands either.

A separate possible worry is that DLAs could instead contaminate our sample of small saturated regions that we use for comparison purposes (as described in §6). Our small saturated sample is meant to reflect absorption around saturated yet ionized regions, and so should not contain significant damping wing absorption. In principle, wings from any DLAs in this sample could influence the transmission around the small saturated regions. It seems unlikely that this is a significant worry, since the saturated yet ionized regions are likely vastly more abundant than even the small column density DLAs and super Lyman-limit systems. In addition, we can simply inspect the profile of the small saturated sample to see whether it shows any hint of damping wing absorption that might arise from either small isolated neutral regions or DLAs.

Although contamination from DLAs does not appear to be a big worry for our tests, a more detailed examination would certainly be warranted if possible neutral islands are discovered in real data. We may also be able to remove DLA-contaminated regions by flagging spectral regions in the Ly- and Ly- forest that have the same redshift as strong metal absorbers, which generally accompany DLAs (see e.g. Wolfe et al. 2005).

## Appendix B: Further Utilizing the Lyβ Forest

In order to infer the presence of the HI damping wing, we would like to compare the stacked transmission outside of plausibly neutral absorption systems to what that transmission would have been in the absence of the damping wing. Up to this point, we have been using the stacked transmission outside of small absorption systems as a proxy for the latter quantity. From there, we argued that any extended excess absorption outside of large, plausibly neutral absorption systems compared to small, likely ionized absorption systems is indicative of the presence of damping wing absorption.

However, we do have another handle on what transmission would be like in the absence of the HI damping wing and that is the transmission in the Ly forest. As discussed in §5.2, the HI damping wing should be significantly reduced in the Ly forest compared to the Ly forest. Specifically, we saw that the damping wing optical depth in Ly falls to less than one at velocity separations from neutral gas. Therefore, at separations greater than this, stacked transmission in the Ly forest should provide information on what the shape of the Ly transmission would have been in the absence of the damping wing, with foreground Ly absorption only altering the stacked Ly transmission by an overall constant . Using information from stacked Ly transmission has the appeal that it does not require using physically different regions of space in order to estimate the damping-wing-less transmission outside of selected plausibly neutral absorption systems. This provides protection from problems arising from unanticipated differences between the small, likely ionized absorption systems and the large, plausibly neutral absorption systems.

Thus, we would like to find a way to estimate the Ly transmission in the absence of the damping wing by using only the Ly transmission. In principle, this can be done by using simulations to model the relationship between the two and generating a (ionization-state-dependent) mapping that takes a measurement of stacked Ly transmission outside of large absorption systems and maps it to an estimate of the damping-wing-less Ly transmission in the same regions. From there, the ratio of the stacked Ly transmission to this estimate of the damping-wing-less Ly transmission would leave us with an estimate of . In the left-hand panel of Fig. 12, we show the recovered curve after applying this approach to each of the ionization states considered thus far, and then normalizing each curve to peak at 1. Specifically, a mapping between stacked Ly transmission and stacked damping-wing-less Ly transmission for each ionization state was generated using a large ensemble of mock spectra and then applied to groups of 20 spectra. The error bars in the figure show the scatter in the estimated damping wing absorption between realizations of 20 spectra. For ease of viewing, we show only the error bars for and 0.35. This figure demonstrates that the approach works well and recovers a damping wing shape for an IGM with (in the absence of noise).

A few things are worth pointing out about this process. First, the recovered damping wing profiles are only useful to the extent that they provide confidence that we are, in fact, observing neutral hydrogen in the IGM. The stacked profile of the HI damping wing is a complicated entity and, as such, we do not expect to be able to, for example, fit the recovered curves to Eq. 6 and estimate . Secondly, it is comforting to note that not only is no damping shape recovered in the case of , but even if a mapping corresponding to a significantly neutral IGM is applied to a measurement of a fully ionized IGM, we do not recover a damping wing shape. Therefore, we do not expect this approach to yield false positives. Lastly, this process relies on simulations in order to map the stacked Ly transmission to damping-wing-less Ly transmission and is therefore somewhat model-dependent. However, we do not expect the specifics of reionization to significantly impact this mapping and we are also not interested in the fine details of the results here. We are primarily interested in whether damping wing absorption can be measured at all in the case of a somewhat neutral IGM and, as such, we are comfortable with this level of model dependence.

Finally, we find that this mapping is relatively simple. Namely, for velocity separations , the stacked Ly transmission and the stacked damping-wing-less Ly transmission differ by roughly a constant multiplicative factor. Thus, in aiming to recover the shape of the stacked damping wing absorption, it appears to be a very good approximation to simply divide the stacked Ly transmission by the stacked Ly transmission. In the right-hand panel of Fig. 12 we simply take groups of 20 spectra and divide their stacked Ly transmission by their stacked Ly transmission (outside of large absorption systems) and give the result unity amplitude. Qualitatively, the results look very similar to those obtained from the mapping procedure (shown on the left-hand side) but are without any model-dependence. Additionally, we again find that, in the case of a fully ionized IGM, we do not recover a damping wing shape. Thus, this provides another check which may be performed with actual spectra in order to bolster confidence that a damping wing is in fact being observed.

A potential concern for our Ly stacking approach in general could be that, while we make the approximation that ionized regions are exposed to a uniform ionizing background, the ionizing background will in fact be fluctuating spatially. It is then possible that, in scenarios where the ionizing background is weaker closer to the stacking location and stronger farther from the stacking location, an extended recovery could be imprinted on the stacked transmission despite the IGM being fully ionized. If one were not careful, and if these spatial fluctuations occurred on scales comparable to the damping wing feature, a false detection could be possible. One tool we have to protect against this is the fact that the scale of the damping wing is significantly smaller in Ly than in Ly. Therefore, any extended recovery in stacked transmission which occurs over similar scales in Ly and Ly is unlikely to be caused by the damping wing.

## Appendix C: Extended damping wing Absorption from Correlated HI Islands

As mentioned in §7.1, and seen in Fig. 6, 7, & 8, stacked Ly transmission outside of large saturated regions in a significantly neutral IGM displays excess absorption extending significantly past the scale of an isolated damping wing. To explain this, we seek to model the expected transmission outside of a neutral region, incorporating the correlation between the neutral region at the origin and neighboring neutral regions. In order to simplify the calculation, while capturing the main effect, we ignore correlations between the neighbors themselves – i.e., we only include the correlation between the neutral region at the origin and the neighbors and ignore inter-neighbor correlations.

Ignoring correlations between the neighboring neutral regions, we can approximate them as following a Poisson distribution and consider the total absorption contributed by these neutral islands or “clouds” following Zuo & Phinney (1993). Suppose that on average clouds contribute to the absorption at a given region of the spectrum. Let denote the transmission when clouds reside along the line of sight and impact the given spectral region, with being the optical depth of the th cloud. If we allow the clouds to be placed independently and if they have equal optical depths, then:

 ⟨Fk⟩ ≡⟨e−τ1e−τ2⋯e−τk⟩=⟨e−τ1⟩⟨e−τ2⟩⋯⟨e−τk⟩=(e−τ)k. (3)

Using this expression, we can calculate the ensemble-averaged transmission by averaging over all the possible numbers of intervening clouds:

 ⟨F⟩ (4) =e−meme−τ=e−m(1−e−τ). (5)

Let us define the quantity as

 e−m(1−e−τ) ≡e−τeff (6) τeff ≡m(1−e−τ). (7)

We now have the absorption from neighboring clouds, characterized by the parameter .

We just need to adapt this slightly to the problem at hand. Suppose that – with certainty – there is a neutral region located at and let us consider the excess absorption, above random, contributed by neighboring neutral regions. First, the optical depth at that is contributed by a neighboring cloud located at will depend on : clouds closer to will have larger optical depths at that corresponding frequency. We can account for this by substituting , according to Eq. 6. Next, the expected number of HI islands in a region with velocity extent nearby our stacking location can be approximated as

 m ≈dv′⟨nHI⟩(1+ξHI,HI(v′)) (8)

where is the average number of HI islands per interval and is the correlation function between the centers of neutral regions separated by . From here, we can model the effective optical depth at a given velocity separation due to neighboring HI islands, , as

 τeff(Δv) =∫dv′⟨nHI⟩(1+ξHI,HI(v′))[1−e−τ(|Δv−v′|)]. (9)

In other words, the excess effective optical depth from neighboring systems involves the convolution of the absorption profile around each region with the correlation function of the neutral regions. This is analogous to the “two-halo” term in the halo model (e.g. Cooray & Sheth 2002).

Thus, the model for the overall stacked transmission outside of neutral regions, also incorporating the damping wing from the central neutral region, becomes:

 ⟨F⟩(Δv) =e−τDW(Δv)e−τeff(Δv). (10)

This model requires two inputs. First, it requires the correlation function between the centers of neutral regions, , which can be calculated from a model of the underlying ionization field. Second, the optical depth profile for neutral regions, , is a function of the size of the neutral regions, per Eq. 6. While the neutral regions in the mock spectra take on a range of sizes, we make the simplifying approximation here – but not in the body of this paper – that the neutral regions at a given neutral fraction have one typical size and denote this with a corresponding extent in velocity space . This effectively results in the optical depth profile of an individual neutral island in the right hand side of Eq. 9 being described by a piecewise function

 τDW(Δv)=⎧⎪ ⎪⎨⎪ ⎪⎩∞if |Δv|

Next, we would like to compare this model against results using mock spectra. We do this by first generating mock spectra which only include absorption from neutral islands, since this is the only type of absorption incorporated in our model, and stack these mock spectra at the HI/HII boundaries. To be clear, while the model curve described above adopted a fixed for the purpose of calculating a , the mock spectra here are generated using the same simulated ionization fields as throughout the rest of the paper, with a wide range of sizes for the underlying neutral regions.

To obtain a model for the stacked transmission, we first calculate the correlation function between the centers of neutral regions using the mock underlying ionization fields and also choose a value for to be used in the optical depth profile. Additionally, in Eq. 9, the term effectively breaks our integral into two pieces: the first representing the mean absorption from neutral regions and the second representing the excess or reduced absorption at due to the clustering of neutral islands. For our case, we are only concerned with the excess absorption, so we make the replacement .

Therefore, by providing a value for and measuring , we can get a estimate for the mean transmission outside of neutral regions which incorporates absorption from spatially-correlated neighboring regions. In the left panel of Fig. 13, we plot an example of this for . We show the modelled damping wing absorption from the central neutral region in blue, the modelled absorption from neighboring neutral islands and their damping wings in cyan, and the product of these in black. For comparison, we show the stacked transmission in the mock spectra in magenta, shifted by to account for stacking occurring at HI/HII boundaries instead of at the center of neutral regions. All curves have been divided by the mean transmission. After taking , we find good agreement between the above model and the stacked transmission. The precise agreement should be taken with a grain of salt, since the model makes several simplifying assumptions, especially that the neutral regions have a fixed size. However, the model does demonstrate that clustered neutral islands should lead to extended excess absorption, significantly beyond the scale of an individual damping wing.

In the right panel of Fig. 13, we show the comparison between the stacked transmission (solid) and modelled transmission (dashed) for (black), 0.22 (blue), and 0.05 (cyan), where each curve has been multiplied by the mean transmission for clarity. In generating these plots, we have assumed , , and for , 0.22, and 0.05, respectively. We again find a very nice agreement between the modelled and stacked transmission, further confirming that spatially-correlated regions are indeed responsible for the significantly-extended excess absorption.

You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters