Hydrodynamic Response of the Intergalactic Medium to Reionization II: Physical Characteristics and Dynamics of Ionizing Photon Sinks

Becker et al. 2021 measured the mean free path of Lyman limit photons in the IGM at $z=6$. The short value suggests that absorptions may have played a prominent role in reionization. Here we study physical properties of ionizing photon sinks in the wake of ionization fronts (I-fronts) using radiative hydrodynamic simulations. We quantify the contributions of gaseous structures to the Lyman limit opacity by tracking the column density distributions in our simulations. Within $\Delta t = 10$ Myr of I-front passage, we find that self-shielding systems ($N_{\rm HI}>10^{17.2}$ cm$^{-2}$) are comprised of two distinct populations: (1) over-density $\Delta \sim 50$ structures in photo-ionization equilibrium with the ionizing background; (2) $\Delta \gtrsim 100$ density peaks with fully neutral cores. The self-shielding systems contribute more than half of the opacity at these times, but the IGM evolves considerably in $\Delta t \sim 100$ Myr as structures are flattened by pressure smoothing and photoevaporation. By $\Delta t = 300$ Myr, they contribute $\lesssim 10 \%$ to the opacity in an average 1 Mpc$^3$ patch of the Universe. The percentage can be a factor of a few larger in over-dense patches, where more self-shielding systems survive. We quantify the characteristic masses and sizes of self-shielding structures. Shortly after I-front passage, we find $M=10^{4} - 10^8$ M$_\odot$ and effective diameters $d_{\rm eff} = 1 - 20$ ckpc$/h$. These scales increase as the gas relaxes. The picture herein presented may be different in dark matter models with suppressed small-scale power.


INTRODUCTION
Hubble Space Telescope (HST) measurements of the rest-frame ultraviolet (UV) luminosity function of z > 6 galaxies have provided a first census of reionization sources (e.g. Finkelstein et al. 2015;Bouwens et al. 2015). These observations give us a broad view of plausible reionization histories that are consistent with the observed star formation history of the Universe (e.g. Robertson et al. 2015;Finkelstein et al. 2019;Bouwens et al. 2021). Of equal importance to our understanding of reionization, however, are the sinks of ionizing photons. The sinks shaped how ionization fronts (I-fronts) progressed (Iliev et al. 2005a) and they played an important role in setting the ionizing photon budget required to complete and maintain reionization (Park et al. 2016;D'Aloisio et al. 2020).
The sinks can be characterized by their distribution of H I column densities. At z = 2 − 5 (post-reionization), the column density distribution has been constrained by numerous quasar absorption spectrum studies (Storrie-Lombardi et al. 1994;Songaila & Cowie 2010;Prochaska et al. 2010;Rudie et al. 2013;Kim et al. 2013;Crighton et al. 2019). At these redshifts, much of the opacity is contributed by optically thick absorbers with columns 10 17.2 < N HI < 10 19 cm −2 , the so-called Lyman-limit systems (LLSs). For example, Prochaska et al. (2010) found that ≈ 55% of the Lyman limit opacity at z = 3.7 is produced by N HI ≥ 10 17.5 cm −2 LLSs. Although the exact nature of LLSs is debated, cosmological radiative transfer (RT) simulations successfully reproduce their observed abundance, and suggest that the LLSs correspond to ∆ g ∼ 100 gas at the outskirts of halos, where ∆ g is the gas density in units of the cosmic mean (Mc-Quinn et al. 2011;Altay et al. 2011Altay et al. , 2013. In contrast to this relatively well-studied picture of the post-reionization IGM, the properties of the sinks during reionization are poorly understood. This owes to a lack of observational constraints at z > 6 as well as the computational challenges in simulating the sinks. On the theoretical side, the simulations of Park et al. (2016) and D'Aloisio et al. (2020, Paper I) suggest that the LyC opacity during reionization was considerably more complicated than at lower redshifts. Before a patch of the IGM was reionized, the gas clumped on a hierarchy of scales down to its Jeans mass, which could have been as low as M J ∼ 10 4 M for gas at temperature T ∼ 10 K. After an I-front swept through a region, the photoionized gas expanded in response to the sudden pressure increase, a process that we term relaxation. In addition, I-fronts became stuck within gaseous halos until they were photoevaporated (Shapiro et al. 2004;Iliev et al. 2005b). These processes were limited by the sound speed in the ionized gas (c s ∼ 20 km s −1 ), so they took place over several hundred Myr. During this period the LyC opacity evolved significantly. Moreover, the non-trivial interplay between the hydrodynamic response and selfshielding renders the evolution dependent on the local intensity of the ionizing background, gas density, and redshift of reionization.
On the observational side, Becker et al. (2021) recently extended to z = 6 direct measurements of the mean free path (MFP) from stacked quasar absorption spectra. At z = 5.1 they found λ 912 mfp = 37.71 +5.31 −6.64 cMpc/h, consistent with the previous measurement of Worseck et al. (2014). At z = 6, the short value of λ 912 mfp = 3.57 +3.09 −2.14 cMpc/h measured by Becker et al. (2021) implies a rapid evolution of the IGM opacity between z = 5 − 6. Cain et al. (2021) used radiative transfer (RT) simulations, applied with a new sub-grid model for the sinks, to argue that the rapid evolution favors a late and rapid reionization process driven by faint galaxies. Davies et al. (2021) assessed that a cumulative source output of 6.1 +11 −2.4 ionizing photons per baryon are required for consistency with both the short λ 912 mfp (z = 6) and upper limits on the IGM neutral fraction from McGreer et al. (2015). Both papers concluded that the sinks played a principal role in shaping reionization if λ 912 mfp (z = 6) is as low as the Becker et al. (2021) measurement. 1 In this paper, we use the suite of radiative hydrodynamics simulations published in Paper I to take a more in-depth look at the absorption systems responsible for the LyC opacity during reionization. We will address four main questions: (1) What are the physical properties of the absorption systems that set the opacity?; (2) How do these properties depend on the local environment?; (3) What are the length and mass scales 1 At face value, the ionizing photon budgets in the models of Cain et al. (2021), ≈ 3 photons per H atom to complete reionization, appear discrepant with the larger budget of 6.1 +11 −2.4 found by Davies et al. (2021). However, the models of Cain et al. (2021) allow λ 912 mfp (z = 6) to be 1σ larger than the central value measured by Becker et al. (2021). And they have global neutral fractions ≈ 20% at z = 6, in ∼ 2σ tension with the dark pixel constraints of McGreer et al. (2015). If the calculations of Davies et al. (2021) are adjusted for these allowances, the agreement is quite good. Under these assumptions, they find a cumulative budget of 2.3 photons per baryon at z = 6, compared to 2.2 in the model of Cain et al. (2021). that characterize these systems?; (4) How do the systems evolve during the relaxation process? A detailed understanding of the sinks' physical nature will provide insight into what absorption systems set the ionizing photon budget for reionization. It will also help us understand the early evolution of the LyC opacity and how it fits with the standard model of the post-reionization IGM. Lastly, the sinks are expected to play an important role in setting the spatial structure of reionization (Miralda-Escudé et al. 2000;Furlanetto & Oh 2005;Mc-Quinn et al. 2007;Mao et al. 2019). Since almost all reionization observables are sensitive at some level to its morphology, an accurate model for the sinks is likely critical for confronting simulations with forthcoming observations. One aim of this work is to inform the further development of such models (e.g. Cain et al. 2021).
The structure of this paper is as follows. In §2 we briefly summarize the simulations of Paper I. In §3, we discuss the relationship between photo-ionization rate and neutral hydrogen density. Section 4 lays out contributions to the LyC opacity of the IGM. In §5, we explore demographics of optically thick absorbers. Section 6 quantities physical properties of the structures responsible for optically thick absorbers. We conclude in §7. Throughout we adopt a ΛCDM cosmology with Ω m = 0.31, Ω b = 0.048, H 0 = 100h km s −1 Mpc −1 , with h = 0.68, σ 8 = 0.82, n s = 0.9667 and a hydrogen mass fraction of X Hy = 0.7547, consistent with the latest measurements (Planck Collaboration et al. 2018). Proper distances are denoted with a "p" prefix (e.g. pMpc), while all other distances are reported in comoving units.

RADIATIVE HYDRODYNAMICS SIMULATIONS
We briefly summarize the main features of the fully coupled RT hydrodynamic simulations from Paper I. The simulations were run using a modified version of the RadHydro code, which combines ray tracing RT with Eulerian hydrodynamics (Trac & Pen 2004;Trac & Cen 2007;Trac et al. 2008). All of the runs have N = 1024 3 dark matter particles, gas cells and RT cells, with box sizes of L = 1.024 h −1 Mpc . This provides a gas/RT cell width of 1 h −1 kpc , small enough to capture reasonably the pre-reionization Jeans scale of the gas. In Appendix A, we present a series of numerical convergence tests (see also Appendix A of Paper I). We will discuss these tests in the relevant contexts below.
The simulations do not explicitly model the sources of reionization. Instead, the gas is ionized by external sources, a process that we model by sending plane parallel I-fronts through the simulation volume. The box is divided into N dom = 32 3 cubic RT domains. Source cells are placed on two adjacent sides of each domain and they are turned on at a specified redshift of reionization, z re (i.e. rays are sent from two directions within each domain). As discussed in Paper I, this setup affords us a clean interpretation of the dynamics; the gas is ionized at nearly the same time and at the same impinging radiation intensity. Appendix A of Paper I demonstrates that the domain setup does not significantly alter the IGM clumping and ionization structure of the reionized gas. In addition, Appendix A of the current paper tests the effect of the domain structure on the mass and size distributions of optically thick absorbers. We will discuss these quantities in detail in §6.1. We adopt a powerlaw spectrum with specific intensity J ν ∝ ν −1.5 between 1 and 4 Ry, in 5 frequency bins (ν is frequency). This is motivated by stellar population synthesis models of metal poor populations (e.g. D'Aloisio et al. 2019).
Our small-scale simulations span a range of z re and LyC intensities to sample the patchiness of the global reionization process (see Table 1 of Paper I). The impinging LyC intensity is parameterized by Γ −12 , the H I photoionization rate in the source cells, expressed in units of 10 −12 s −1 . Our runs assume Γ −12 = 0.3 and 3.0. The former is consistent with Lyα forest measurements (e.g. D'Aloisio et al. 2018) just after reionization (z = 5 − 6), or perhaps during its tail end. The latter is intended to model the intensity near an overdensity of sources. The simulations also include a set of "DC mode" runs that model fluctuations away from the mean density on the box scale (Gnedin et al. 2011). The box-scale density is parameterized by δ/σ, the linearly extrapolated density contrast smoothed on the box scale, in units of its standard deviation. 2 In addition to δ/σ = 0, i.e. cosmic mean density runs, we consider δ/σ = ± √ 3, which correspond to present-day linearly extrapolated over-densities of δ 0 = ±5.013. 3

PHOTOIONIZATION RATE VS DENSITY
We begin with the relationship between the photoionization rate and local hydrogen number density (n H ), which characterizes self-shielding in our simulations. Earlier works have studied this relationship in full cosmological simulations of reionization (Rahmati et al. 2013;Chardin et al. 2018). In this section, we examine how the hydrodynamic response of the IGM to photoheating drives evolution in the local relationship between Γ HI and n H during patchy reionization.
In Figure 1 we illustrate several key dependencies of the Γ HI vs. n H relationship. The curves in the left column show the median Γ HI vs. n H at ∆t = 10, 60, and 300 Myr from z re . The top and bottom panels correspond to Γ −12 = 0.3 and 3.0, respectively, both with z re = 8 and δ/σ = 0. For simplicity, we focus on these two simulations but the trends described here hold more generally. Consider the case with Γ −12 = 0.3 (top-left). For all ∆t, Γ HI reaches 0 at approximately the same log(n H /cm −3 ) ≈ −1.9. However, Γ HI falls more steeply for times closer to z re . The bottom-left panel shows the same general trend. A more intense ionizing background simply moves the self-shielding cutoff to higher densities. 4 The right column of Figure 1 shows how the relationship depends on z re for fixed Γ −12 = 0.3. The different curves correspond to z re = 12, 8, and 6 at ∆t = 10 (top) and 300 (bottom) Myr after z re . Shortly after Ifront passage, the decline in Γ −12 is steeper at lower z re . The bottom panel shows, however, that the Γ HI vs. n H curves approach a nearly identical form by ∆t = 300 Myr. As we will now discuss, all of these trends may be understood qualitatively in terms of the density profiles of the sinks, and how they are shaped by the competing effects of relaxation and structure formation.
A toy model can provide insight into the results of Figure 1. Consider a spherically symmetric absorber with a density profile, wheren H (z) is the proper cosmic mean hydrogen number density, r 0 characterizes the size of the absorber, and α sets the steepness of the density profile. We have normalized the profile such that n H = 50n H at r = r 0 . We also take our fiducial r 0 to be 10h −1 kpc. Both the over-density and r 0 are broadly motivated by the characteristic over-densities at which self-shielding occurs and the sizes of absorbers in our simulations (see §4.3 of Paper I, as well as §6 of the current paper). However, the ensuing qualitative discussion does not depend on the particular choice of normalization. We expose this absorber to an external isotropic ionizing background and solve the one-dimensional RT equation to obtain Γ HI vs. . The black/solid, red/dashed, and blue/dot-dashed curves correspond to z = 7.9, 7.5, and 6.0, respectively, ∆t = 10, 60 and 300 Myr after zre . Bottom Left: Same as top-left but with Γ−12 = 3.0. In the text we provide a simple model for how photoevaporation and relaxation change the ΓHI vs. nH relationship, in terms of density profiles and sizes of absorption systems.
Top-Right: The ΓHI vs. nH relationship at a fixed time interval, ∆t = 10 Myr, from zre = 12 (black/solid), 8 (red/dashed), and 6 (blue/dot-dashed). Bottom-right: Same but for ∆t = 300 Myr. The top-right panel shows differences that arise due to structure formation (see text). The bottom-right panel shows that these differences are largely erased by the relaxation/photoevaporative processes.
n H under the assumption that the gas is in photoionization equilibrium. 5 Figure 2 shows results from this model at z = 6 for several combinations of α and r 0 , where the latter is given in comoving units. We have fixed the intensity of the impinging ionizing background such that the hydrogen photoionization rate far away from the absorber is Γ bk = 3 × 10 −13 s −1 . The black curve with 5 In detail, the photo-ionization rate towards the center of the absorber is Γ HI (r) = Γ bk exp(−τ (r)) where the opacity along the radial direction towards the center is given by If the gas everywhere inside the absorber is in photo-ionizational equilibrium with the radiation, then we have Γ HI (r)n HI (r) = α B (T )(1 + χ)(n H (r) − n HI (r)) 2 These equations can be solved iteratively to obtain the profile of Γ HI for an assumed density profile. For simplicity, we assume a uniform temperature of T = 10 4 K.
(α, r 0 ) = (2, 10 kpc/h) corresponds to an isothermal density profile with size representative of the simulated systems that we will discuss in §6. The other curves vary α and r 0 to mimic the effects of relaxation, photoevaporation, and structure formation. Increasing r 0 results in self-shielding setting in at smaller n H , i.e. Γ HI falls off at smaller n H . This occurs because a column density of N HI = 1/σ 912 is reached at smaller n H for larger values of r 0 . (Here, σ 912 is the photoionization cross section of hydrogen at 1 Ry, or ≈ 912Å.) If the black/solid curve is a density peak at high redshift, the green/dotted curve represents that same peak at lower redshift after structure formation has increased its mass. This provides a qualitative understanding of the trends seen in the top-right panel of Fig. 1. For lower z re , structure formation has had more time to grow the sizes of the density peaks that serve as absorbers. We thus see Γ HI falling off at smaller n H .
At fixed r 0 , flattening the density profile (decreasing α) also causes Γ HI to fall off at smaller n H . This effect, combined with possible evolution in r 0 during relax- A toy model that we use to understand how photoevaporation and relaxation drive changes in the ΓHI − nH relation. This model assumes a spherically-symmetric power-law density profile illuminated from the outside by an isotropic ionizing background. The gas is assumed to be in photoionization equilibrium with the background. We show results for z = 6. Illustrative values for the power-law slope and absorber size are displayed. If photoevaporation makes the absorber density profiles shallower, and characteristically larger, for example, self-shielding will set in at lower nH.
ation/photoevaporation, provides a simple picture for the evolution with ∆t seen in the left panels of Fig. 1. The relaxation and photoevaparation processes alter the density profiles of the absorbers, changing α and r 0 and driving the statistical evolution in Fig.1. We note that these changes are complicated; they likely depend on the details of the initial configuration (Shapiro et al. 2004). Finally, while this single-absorber model provides insight into the trends seen in our simulations, we caution that the situation is more complicated in reality with a diverse population of absorbers. For example, shortly after I-front passage much of the absorption may occur in diffuse gas within filaments (c.f. Figs. 3 and 4 in Paper I), in which case spherical symmetry is a poor approximation. The model nonetheless provides some intuition for how both structure formation and the hydrodynamic response to photoheating drive evolution in the Γ HI vs. n H relation.

CONTRIBUTIONS TO THE LYC OPACITY
Paper I explored the evolution of λ 912 mfp in ionized gas during the relaxation process. In this section, we aim to better understand what sets λ 912 mfp by decomposing the Lyman limit opacity into its constitute absorption systems. In Paper I, the MFP was calculated by tracing skewers of length L dom = 32 h −1 kpc along each coordinate axis through the RT domains, and evaluating the outgoing flux f out = e −τ912 along each skewer, where τ 912 = σ 912 N HI . Then the MFP was obtained using λ 912 = −L dom / ln( f out ), where the average is taken over all segments. It is straightforward to show that, in the limit 1 − f out << 1, the effective absorption coefficient in our simulation volume, κ 912 ≡ λ −1 912 , can be expressed as where the integral runs over all column densities and is the number of segments per unit column, per unit length, i.e. the column density distribution. This conveniently re-casts the MFP calculation as an integral over column densities, allowing us to write the absorption coefficient per logarithmic interval in N HI as (Rahmati & Schaye 2018) where we take ∆ log(N HI ) = 0.1. Through eq. (7), we quantify the contribution to the opacity from different column densities, which will give us insight into the properties of the gaseous structures responsible for setting the MFP. Equation (5) neglects spatial correlations between the absorption systems; it is identical to the expression for κ 912 in the widely used model of Poisson distributed absorbers (Paresce et al. 1980). Recognizing that this expression can be derived in the appropriate limit from our definition of λ 912 mfp in Paper I, we find that the former is quite accurate in practice.
In Figure 3 we show, from left to right, the column density distribution (Eq. 6), the absorption coefficient per ∆ log(N HI ) (Eq. 7), and the cumulative fractional contribution to the absorption coefficient (Eq. 5), as functions of N HI , for the 32 kpc/h domain segments from our simulations. The simulation parameters are indicated in the left-most panel in each row. The black, red, and blue curves correspond to ∆t = 10, 60, and The black, red, and blue curves correspond to ∆t = 10, 60, and 300 Myr after zre , respectively. Center: Contribution to the absorption coefficient, κ912, per logarithmic interval in NHI. Right: Cumulative κ912(< NHI) expressed as a fraction of the total κ912. The fractional contribution to the opacity from columns above or below a given value can be read off of these panels. 300 Myr after z re . (For the z re = 6 case, the blue corresponds to ∆t = 244 Myr, since we did not run the simulation past z = 5.) In each panel, the top axis shows τ 912 and the vertical dashed lines denote τ 912 = 1.
Some key qualitative features are evident for all the simulations shown in Figure 3. There is a distinct peak around N HI = 10 19 cm −2 in the column density distribution of recently ionized gas (black curves) in all but the Γ −12 = 3.0 case. This peak disappears by ∆t = 300 Myr except in the over-dense run (δ/σ = √ 3), although it is suppressed significantly there as well. The steady decrease in the abundance of segments with columns 10 15 cm −2 reflects the hydrodynamic response of the gas after ionization. At early times, there is abundant small-scale structure with τ 912 ∼ 1 but it is erased over ∆t ∼ 100 Myr by Jeans smoothing and photoevaporation.
The second column shows that dκ 912 /d log N HI exhibits two prominent peaks shortly after I-front passage (∆t = 10 Myr) -the first at N HI ∼ 10 17 cm −2 and the second at N HI ∼ 10 19 cm −2 . As we will discuss below, these peaks represent two physically distinct types of optically thick absorbers. The relative heights of the peaks -and therefore their relative contributions -are set by the interplay between all three simulation parameters. Generally, lower Γ −12 , lower z re , and larger environmental density (δ/σ) all increase the contribution of the N HI ∼ 10 19 cm −2 peak. However, as the gas relaxes, this peak all but disappears except in the over-dense run. Most N HI 10 17 cm −2 systems do not survive the relaxation/evaporation except where the local over-density is large. 6 From the right-most panels, we can glean the fractional contribution of optically thick systems to the κ 912 of our boxes. The contributions vary significantly with the environmental parameters. For example, in the second panel with (z re , Γ −12 , δ/σ) = (8, 0.3, 0), the τ 912 > 1 systems account for 47% of the opacity at a time ∆t = 10 Myr after I-front passage. At the same ∆t, they account for only 27% when Γ −12 is a factor of 10 higher (bottom panel). In the over-dense run (fourth panel), they account for 68%. Notably, τ 912 > 1 systems still account for 35% of the opacity in this run when ∆t = 300 Myr. This is in contrast to all other runs, for which the contribution from τ 912 > 1 systems becomes ≤ 15 % in the relaxed limit.
By the time a few hundred Myr have passed since z re , the local MFP in an average ∼ 1 h −1 Mpc patch of the IGM is set by low column-density gas. For example, in our simulation with (z re , Γ −12 , δ/σ) = (8, 0.3, 0), approximately 80% of the opacity comes from sight lines with N HI < 10 16 cm −2 when ∆t = 300 Myr. This owes to almost all of the high-column systems being wiped out by the relaxation process. Although the local MFP in this relaxed 1h −1 Mpc patch is set mainly by optically thin gas, the global MFP becomes regulated by optically thick absorbers in over-dense patches, recently reionized patches, or by neutral gas that has yet to be reionized. Indeed, the fourth row of Fig. 3 illustrates that optically thick absorbers have greater longevity in over-dense regions. We conclude by noting that cosmic expansion and an increasing ionizing background intensity (Γ −12 ) both act to increase the contribution of LLSs to κ 912 well after reionization. Future work should explore how the picture presented here transitions to the well-studied opacity structure of the z = 2 − 4 IGM.

DEMOGRAPHICS OF OPTICALLY THICK
ABSORBERS DURING REIONIZATION Figure 3 shows that the distributions of column densities in our simulations are bi-modal throughout much of the relaxation process. Here we illustrate that the bi-modality arises from two distinct types of absorption systems in recently reionized gas.
The sight lines with the highest column densities, which comprise the peaks located around log N HI /[cm −2 ] ≈ 19, almost always contain absorbers with neutral fractions close to unity. Note that these high-N HI peaks disappear over a time scale of ∼ 100 Myr because the absorbers are photo-evaporated by the ionizing background impinging on their outskirts. To illustrate this, the left panels of Figure 4 show an example sight line intersecting an absorption system undergoing photo-evaporation. The skewer is drawn from one of the RT domains in the simulation with (z re , Γ −12 , δ/σ) = (8, 0.3, 0). From top to bottom, we show the gas density (∆ g ; in units of the cosmic mean), the H I fraction (x HI ), the gas temperature (T ), the photo-ionization rate, and the ratio of the photo-ionization rate to the recombination rate, P/R ≡ Γ HI /(R/n HI ), where R = α B (T )n e n HII , and α B is the case B recombination coefficient of hydrogen. The black, red, and blue curves correspond to snapshots in time at ∆t = 10, 60, and 300 Myr after z re , respectively. The total column density of the sight line evolves from log N HI /[cm −2 ] = 18.8 at ∆t = 10 Myr, to 16.3 and 14.7 at ∆t = 60 and 300 Myr, respectively. At ∆t = 10 Myr, the absorber located near the center of the sight line is peaked in density. (Though, the sight line does not necessarily pass through the center of the peak in three dimensions.) The neutral frac- . From top to bottom, the rows show the gas overdensity, H I fraction, temperature, H I photoionization rate, and the ratio of photoionization rate to recombination rate. The last quantity is unity for gas in photoionization equilibrium. Each panel shows the time evolution of the sight line with the black, red, and blue lines corresponding to ∆t = 10, 60,and 300 Myr after zre , respectively. The sight line on the left has column densities of log NHI/[cm −2 ] = 18.8, 16.3, and 14.7 at ∆t = 10, 60, and 300 Myr, respectively. On the right, we have log NHI/[cm −2 ] = 17.3, 16.3, and 15.2, respectively. At ∆t = 10 Myr, we classify the absorber on the left as "evaporating," characterized by its neutral core and the presence of gas out of photoionization equilibrium. We call the absorber on the right "relaxing." This sight line is highly ionized, yet optically thick at ∆t = 10 Myr. Unlike the case on the left, the I-front likely never transitioned to D-type as it swept past this sight line. tion reaches unity near the center, where Γ −12 dips to zero -the result of self-shielding. And the temperature of the peak is ∼ 1, 000 K, much lower than the photoionized gas outside of it. A key feature of this system is that the gas at the boundaries of the neutral core is out of photo-ionization equilibrium, as evidenced by the spikes in P/R shown in the bottom panel. By ∆t = 60 Myr, the density peak has been smoothed out considerably. The gas is highly ionized and has reached photoionization equilibrium, with the equilibration time scale being ∼ 1/Γ HI ∼ 10 13 s = 3 × 10 5 years at the center of the peak (see 3rd row). The gas is evacuated from the system by ∆t = 300 Myr and the absorber is gone. 7 In what follows, we we will refer to systems like these, which exhibit large deviations in P/R from unity, as "evaporating systems." The situation is considerably different in the right panels. The column densities of this sight line are log N HI /[cm −2 ] = 17.3, 16.3, and 15.2 at ∆t = 10, 60 and 300 Myr, respectively. Despite the sight line having an opacity of τ 912 ≈ 1.2 at ∆t = 10 Myr, there is never any fully neutral core associated with the central density peak. In fact, the peak is highly ionized and in photoionization equilibrium at all snapshots shown. However,  Figure 5. The contribution of segments containing evaporating (red/dashed) and relaxing (blue/dashed) absorbers to the differential opacity, dκ912/logNHI , shown as solid black curve. The distributions are from the (Γ−12, zre, δ/σ) = (0.3, 8, 0) simulation. The optically thick absorbers consist of two distinct populations of sinks: (1) density peaks with fully neutral cores (evaporating); (2) highly ionized absorbers in photo-ionization equilibrium with the ionizing background (relaxing). this peak, too, is erased by Jeans pressure smoothing within ∆ = 300 Myr. We will refer to this type of system as a "relaxing" absorber. We note that the key distinction between these two classifications is the timescale over which I-fronts are able to ionize all of the gas in the system. Evaporating systems are those in which the I-fronts get "stuck" as they climb the steep density gradient, transitioning from R-type to D-type (see e.g. Shapiro et al. 2004). In contrast, relaxing systems may have been quickly ionized as I-fronts swept supersonically through the region, perhaps never transitioning to D-type.
To support the picture that the absorbers consist of two distinct populations, Figure 5 shows what happens when we select sight lines based on whether the gas is in photo-ionization equilibrium. The black/solid curve shows dκ 912 /d log N HI at ∆t = 10 Myr from the full sample of sight lines in our simulation with (Γ −12 , z re , δ/σ) = (0.3, 8, 0). As a proxy for evaporating systems, the red/dashed curve corresponds to sight lines with |P/R| ≥ 1.5 anywhere along the sight line. 8 The blue/dashed curve corresponds to those with |P/R| < 1.5. That the two peaks separate cleanly suggests that the high-column peaks in Fig. 3 are comprised almost entirely of evaporating systems with fully neutral gas. Although we show an illustrative example here, we find that this separation occurs for all the simulations in our suite. In the next section, we will explore in more detail the physical properties of these systems.

PROPERTIES OF ABSORPTION SYSTEMS
Here we explore the properties of the 3-dimensional structures responsible for the 1-dimensional absorbers studied above. Accomplishing this requires a prescription for identifying self-shielding structures. Following Rahmati & Schaye (2018), we define a self-shielding structure (or "absorber") to be a simply-connected group of cells with neutral fraction above a threshold x thresh HI . We use a variable threshold (which depends on Γ −12 , z re , δ/σ) chosen to optimally link the population of τ 912 ≥ 1 absorbers to the 3D systems that create them. 9 Appendix B describes our procedure for calculating x thresh HI . Essentially, we optimize x thresh HI to maximize completeness of our sample, while also minimizing the number of optically thin 1-dimensional sight lines through our absorbers. We then form simply connected groups out of the cells above x thresh HI . ≈ 0.001 for Γ −12 = 3 or δ/σ = − √ 3. To illustrate our identification procedure, Figure 6 shows example slices through the density field from our run with (Γ −12 , z re , δ/σ) = (0.3, 8, 0). From left to right, the panels exhibit the effects of relaxation on the density structure. Self-shielding structures are marked by red circles or blue triangles, with the sizes of the symbols proportional to the effective diameters, d eff ≡ 2(3V /4π) 1/3 , where V is volume. Red circles correspond to evaporating systems (with |P/R| ≥ 1.5 somewhere in the system) and blue triangles correspond to relaxing systems (|P/R| < 1.5). Visual inspection of Fig 6 reveals that the self-shielding structures trace out the filaments and tend to be clustered around the nodes. By ∆t = 300 Myr (right panel) the few systems that remain are at the highest density peaks and are classified as evaporating because they possess fully neutral cores. Figure 7 shows the time-evolution of absorber mass (left column) and size (right column) in three of our simulations, where the former includes both dark matter and gas mass. To obtain dark matter masses, we smooth particles onto the hydro grid using cloud-in-cell interpolation. We characterize the sizes of the systems z=7.9 z=7.5 z=6.0 0.5 0.0 0.5 1.0 1.5 2.0 log g Figure 6. Visualization of absorbers identified in our simulation with (zre, Γ−12, δ/σ) = (8, 0.3, 0). See main text for details on our identification procedure. Gas over-density is shown in color scale and each slice is 32h −1 kpc thick. The red circles and blue triangles show evaporating and relaxing systems, respectively. The size of each symbol is proportional to the effective diameter of the corresponding absorber. The same slice is shown at z = 7.9 (left), 7.5 (center) and 6.0 (right), corresponding to ∆t = 10, 60, and 300 Myr. The surviving systems (right panel) are located at the highest density peaks and contain fully neutral cores, hence their classification as evaporating.

Masses & Sizes
by their effective diameters, d eff . Note, however, that many are not spherically symmetric, especially early in the relaxation process.
The typical mass and size of an absorption system shifts from smaller to larger values during the relaxation process. In our runs with Γ −12 = 0.3 (top two rows), the masses and sizes range from M ∼ 10 4 − 10 8 M and d eff ∼ a few − 20h −1 kpc at ∆t = 10 Myr. By ∆t = 300 Myr, most systems have M ∼ 10 7 − 10 9 M and d eff ∼ 5 − 30h −1 kpc. The shift to larger systems with time is straightforward to understand. Just after I-front passage, the gas contains a large population of small self-shielding systems, but they evaporate/relax quickly. The larger systems survive for longer time scales. The middle row shows results from one of our over-dense simulations with δ/σ = √ 3. At ∆t = 10 Myr (z = 7.9), the distribution of absorber masses appears similar to the case with δ/σ = 0. Note, however, that the local density enhancement results in a couple of very massive M ∼ 10 10 − 10 11 h −1 M systems. Indeed, although not visible in the middle-right panel, these rare systems have d eff ∼ 150h −1 kpc and are in the extreme tail of the distribution. By ∆t = 300 Myr, the absorbers in over-dense box span a broader range of masses, with the high-mass tail of the distribution extending up to > 10 10 h −1 M . The density enhancement also results in significantly more surviving absorbers. For example, at ∆t = 10 Myr, we find 4, 533(4, 069) absorbers in the run with δ/σ = 0( √ 3). By ∆t = 300 Myr, the over-dense run has 98 surviving absorbers, a factor of ≈ 2 more than the 52 found in the mean-density box. Turning to the bottom row, higher Γ −12 generally eliminates the systems with lower mass. This occurs because the smaller, milder over-densities are not able to selfshield against the intense background and are erased on a quicker timescale. By ∆t = 300 Myr, the surviving systems are somewhat smaller than the case with Γ −12 = 0.3 because self-shielding occurs at higher densities. We find 1, 519(11) absorbers in the Γ −12 = 3.0 run at ∆t = 10(300) Myr.
We quantify the ratio of gas mass to total mass for the absorption systems shown in Figure 7. From top to bottom, the mean fractions are 0.10, 0.12 and 0.092, respectively, at ∆t = 10 Myr. Note that the simulation with δ/σ = √ 3 has a slightly higher mean fraction compared to the δ/σ = 0 case. On the other hand, the run with enhanced ionizing background, Γ −12 = 3.0, exhibits a reduced fraction. Most of the systems are below the cosmic baryon fraction, Ω b Ωm−Ω b = 0.187. At ∆t = 60 Myr, the fractions drop to 0.075, 0.091 and 0.072. By ∆t = 300 Myr, there are fewer systems and a large scatter in gas fraction. The mean values actually increase to 0.083, 0.145 and 0.187, respectively, owing the a small number of systems with large (≈ 0.2−0.5) gas fractions.
One concern about the results of Fig. 7 is that the distribution of absorber sizes might be affected by the and effective diameters (right) in three of our simulations. Each row corresponds to a different simulation with parameters provided in the legends of the right panels. The black/solid, red/dashed, and blue/dotdashed curves correspond to ∆t = 10 ,60 and 300 Myr, respectively. The effective diameter is d eff = 2(3V /4π) 1/3 , where V is absorber volume. The absorption systems responsible for setting the LyC opacity during reionization can be as small as the Jeans scale of the cold, pre-reionization gas. In a cold-dark-matter dominated universe, the smallest absorbers may be beyond the resolution limits of our simulations. domain structure of our RT. For example, one could imagine that a source plane along the boundary of a domain might intersect an absorber, causing its inner regions to be artificially ionized and splitting the absorber into two. We tested explicitly for such effects. In Appendix A we show that we obtain similar absorber size distributions even when using larger domains. Of greater concern, however, is the fact that our simulations do not capture the smallest absorption systems that exist in an adiabatically cooling IGM before I-front passage. Based on the numerical convergence tests also presented in Appendix A, we caution that we may be over-estimating the typical mass of absorption systems by a factor of ≈ 2, and the typical sizes by ∼ 25%. In a cold-dark-matter dominated universe, the smallest absorbers may even be beyond the resolution limits of our simulations 6.2. Over-densities Figure 8 shows scatter plots of the n HI -weighted mean gas over-density of the absorbers vs. their column densities. To calculate both quantities, we integrate over randomly oriented skewers traced through the peaks of n HI in the absorbers. We average over the skewers to obtain a single value for each absorber. As in Fig. 6, the red circles and blue triangles correspond to evaporating and relaxing absorbers, respectively. For brevity, we show results for our fiducial simulation with (z re , Γ −12 , δ/σ) = (8, 0.3, 0), but the results are qualitatively similar in our other simulations.
The progress of time from top to bottom illustrates how relaxation depletes the abundance of absorbers that exist in the unrelaxed IGM following I-front passage. At all snapshots, ∆ g increases with N HI except for a flattening that occurs around log(N HI /[cm −2 ]) ∼ 17−19. This is qualitatively similar to the self-shielding "plateau" between log(N HI /[cm −2 ]) ∼ 18 − 20 seen in the simulations of McQuinn et al. (2011). At ∆t = 10 Myr (top panel), self-shielding sets in at the characteristic over-density of ∆ g ∼ 25. However, by ∆t = 300 Myr (bottom), almost all of the absorbers have ∆ g 100 and log(N HI /[cm −2 ]) 19.5. Most also possess neutral cores as indicated by their classification as evaporating systems. Although not shown in Fig. 8, larger Γ −12 values shift the scatter plots upwards in ∆ g and further deplete the absorber population. Larger box-scale densities, δ/σ, have the opposite effect.
The unrelaxed absorbers at ∆t = 10 Myr exhibit a higher degree of scatter compared to later times. Notably, there are absorbers that deviate from the main ∆ g vs N HI relation. In fact, there even appear to be under-dense absorbers. These systems result from the complete or partial shadowing of the LyC radiation by nearby over-densities. For example, a region in the shadow of a large density peak may remain neutral until the peak is sufficiently smoothed for the I-front to break through. Although these occurrences are already rare in our simulations, we caution that our 2-direction RT may exaggerate their prevalence. More realistically, radiation could come from many directions as nearby sources turn on during reionization.

Ionization and Thermal states
An informative way to characterize the state of absorbers in our simulations is through the photoionization equilibrium condition in highly ionized gas, where χ = n He /n H ≈ 1.08 accounts for singly ionized He, α B is the case B recombination rate of H, and n H is the cosmic mean H number density. Note that, because of our RT domain setup, Γ HI is approximately uniform throughout highly ionized gas in our simulations. The utility of characterizing the absorbers with equation (8) is that, for gas in equilibrium, the right-hand side should depend only on T , to a good approximation. To the extent that recently ionized gas has a relatively narrow range of temperatures, we therefore expect x HI /∆ g to take on a narrow range of values for relaxing systems. This should not be true for evaporating systems, however.
In Fig. 9 we show x HI /∆ g versus Γ HI , T , and ∆ g , for our identified absorbers. All quantities are n HIaveraged. 10 For brevity the figure contains results only from the simulation with (Γ −12 , z re , δ/σ) = (0.3, 8, 0), but the same basic trends hold in all of our simulations. The top and bottom panels correspond to ∆t = 10 and 300 Myr after z re , respectively. The evaporating systems (red circles) show a large range of values, whereas the relaxing systems (blue triangles) reside in a smaller region of the parameter space, as we had anticipated above.
In the left panels, the relaxing systems mostly have mean Γ HI near the background value of Γ HI = 3 × 10 −13 s −1 . In contrast, the evaporating systems span a range of Γ HI ; at z = 6 they extend down to Γ HI = 10 −18 s −1 . This is consistent with the evaporating systems possessing neutral cores. The trend of x HI /∆ g increasing with Γ HI for evaporating systems reflects that lower densities are less self-shielded against the background. In the top Photoevaporting Relaxing Figure 8. Relationship between nHI-weighted mean overdensity (∆g) and column density (NHI) for self-shielding absorption systems. Absorbers are identified as simply connected regions above a neutral fraction threshold, as described in the main text. The red circle and blue triangles correspond to evaporating and relaxing absorbers, respectively. We show results from our simulation with (zre, Γ−12, δ/σ) = (8, 0.3, 0). From top to bottom, the panels correspond to ∆t = 10, 60, and 300 Myr after zre. The hydrodynamic response of the gas to reionization depletes the total number of absorbers and raises the characteristic density above which self-shielding occurs.
panel, the scatter towards x HI /∆ g ∼ 1 owes to mildly over-dense gas shadowed from the background by nearby density peaks.
In the middle column, the evaporating systems span a range of temperatures, from ∼ 4, 000 − 32, 000 K. The lower temperatures are similar to virial temperatures of halos with M ∼ 10 6 M . The higher temperatures may result from very recently ionized, low-density gas  or from the collapse of dense systems. In contrast, the relaxing systems are all clustered around temperatures characteristic of photoionized gas in equilibrium with the background, T ∼ 14, 000 K. In all panels at z = 7.9 (∆t = 10 Myr), there are transitional systems bridging the evaporating and relaxing populations. The populations become more separated at late times.
In the right columns, we see that the relaxing systems, with characteristic over-densities ∆ g ∼ 100, reside at the lower end of the ∆ g -distribution. The few evaporating systems with ∆ g significantly less than 100 correspond to diffuse neutral gas in the shadows of nearby peaks. A large fraction of evaporating systems exhibit a tight correlation with ∆ g . This owes to the fact that x HI /∆ g ∼ 1/∆ g for largely neutral systems. By z = 6.0, the surviving optically thick absorbers are highly overdense; almost all are fully neutral systems centered on halos. The analysis presented in the three preceding sub-sections broadly supports our earlier discussion of two distinct populations of optically thick absorbers.

CONCLUSIONS
We have studied the properties of ionizing photon sinks during reionization. Our study was based on the suite of high-resolution radiative hydrodynamics simulations of Paper I, which model the self-shielding and hydrodynamic response of the IGM in the wake of I-fronts. Our main findings can be summarized as follows: • The density and ionization structures of the IGM evolve considerably after z re as the gas relaxes and density peaks are photo-evaporated. We have quantified these effects with the Γ HI vs. n H relationship, which is often used to model selfshielding in cosmological simulations without RT (Rahmati et al. 2013;Chardin et al. 2018). The evolution can be understood qualitatively with a simple model that considers the competing effects of gravitational growth and density-peak smoothing by the ionizing background.
• We measured the column-density distributions in our simulations to characterize the absorption systems that contribute most to the LyC opacity. The top and bottom rows correspond to two snapshots in time, z = 7.9 and 6, respectively. The evaporating and systems are shown as red circles and relaxing systems are shown as blue triangles. From left to right, we show the nHI-weighted mean H I photoionization rate, temperature, and gas over-density of absorbers, vs. xHI/∆g. These quantities are calculated by tracing randomly oriented skewers passing through the peak nHI in each absorber. Relaxing systems are relegated to a relatively small region of the parameter space, consistent with highly ionized gas in equilibrium with the ionizing background.
• For ∼ 50 Myr after z re , the τ 912 ≥ 1 absorbers are comprised of two distinct populations; (1) ∆ g 100 systems with fully neutral cores; (2) Highly ionized systems at milder over-densities, which are in photo-ionization equilibrium with the ionizing background. The latter are just dense enough to be self-shielding but were not dense enough at z re to halt the passing I-front. The relative abundances of these two populations depends on z re , Γ −12 , and δ/σ. After a couple hundred Myr, relaxation and photo-evaporation erase the vast majority of these structures. The surviving structures have ∆ 500, fully neutral cores, and are clustered preferentially in large-scale overdensities.
• The characteristic masses and sizes of τ 912 ≥ 1 absorbers evolve from smaller to larger during the relaxation process. Within ∼ 50 Myr after z re , total masses and sizes are typically in the range M = 10 4 − 10 8 M and d eff = a few − 20h −1 kpc. However, by ∆t ∼ 300 Myr, they are M = 10 7 − 10 9 M and d eff = 5 − 30h −1 kpc, respectively. Again, these properties depend on the environmental parameters z re , Γ −12 , and δ/σ.
Finally, we note some implications for the small characteristic masses and sizes of sinks in a cold-dark-matter universe. First, these scales highlight the extreme computational challenge of modelling the sinks in cosmological simulations of reionization, especially during the first ∼ 50 Myr of the relaxation process. This consideration is particularly relevant for modeling the high-z IGM if reionization ended around z = 5, with a significant fraction of the gas still relaxing at those times (e.g. Kulkarni et al. 2019;Keating et al. 2020;Nasir & D'Aloisio 2020). Second, it is worth noting that the broad conclusions presented here could be considerably different in other dark matter cosmologies with a damping scale in the primordial power spectrum, e.g. warm dark matter, fuzzy dark matter, and some varieties of self-interacting dark matter. In this sense, measurements of the IGM opacity near reionization can, in principle, test dark matter models at scales and densities that are not probed by other methods.

A. NUMERICAL CONVERGENCE
In this section we extend the discussion on numerical convergence in Appendix A of Paper I.

A.1. Resolution
To test convergence with respect to grid size, we have run a series of test simulations in a smaller box with L box = 256h −1 kpc. The grid sizes span N = 64 3 to N = 1024 3 in factors of 8. We use N dom = 8 3 to match the (32h −1 kpc) 3 RT domain sizes in the productions runs. In what follows, we adopt (z re ,Γ −12 , δ/σ ) = (8, 0.3, 0). All quantities plotted in Figure 10 are at z = 7.9, corresponding to ∆t = 10 Myr. We choose this early time in the relaxation process because the gas still clumps on small scales soon after I-front passage. The numerical convergence will generally be better at later times after pressure smoothing has time to act.
The top-left panel of Figure 10 compares the cumulative opacity, κ 912 (< N HI ), in our test runs. The curves corresponding to N = 512 3 and N = 1024 3 agree to better than 10%, indicating that the latter may be reasonably converged. The N = 256 3 case (red curve) has the same spatial resolution as our production runs. While the shape of this curve at columns larger than log N HI /[cm −2 ] ≈ 17 is similar to the N = 1024 3 case, the normalization is lower by 20%. We note that this is consistent with the convergence study of λ 912 mfp in Appendix A of Paper I. Importantly, the relative contribution of τ 912 > 1 sight lines to the total κ 912 appears to be well converged at our production resolution. In the N = 256 3 case, τ 912 > 1 sight lines constitute 43 % of the total opacity, compared to 41 % in the N = 1024 3 case. Thus our conclusions in the main text about the relative contributions of gas at different columns appear to be robust to grid-size convergence.
The top-right panel of Figure 10 considers the median relationship between Γ −12 and n H (see §3). The neutral hydrogen density corresponding to half background intensity (i.e. Γ HI /Γ bk = 0.5) appears to be reasonably converged, though we do underestimate the n HI where Γ HI drops to zero. The bottom panels compare the distribution of masses and sizes for self-shielding structures (or "absorbers"). In this section we define the absorbers to be connected groups of cells above a neutral fraction threshold of x th HI = 0.005. We refer the reader to §6 and Appendix B for the motivation of this choice. The bottom-left panel of Fig. 10 shows that lower hydro/RT spatial resolution leads to the absorbers being larger in mass, on average. The mean absorber masses for the N = 64 3 , 128 3 , 256 3 , 512 3 , and 1024 3 runs are 10 7.64 , 10 7.02 , 10 6.64 , 10 6.45 , and 10 6.29 M , respectively. The corresponding mean d eff are 11.41, 5.15, 4.04, 3.52, and 3.08 h −1 kpc. Comparing the distributions indicates that our production simulations are missing smaller self-shielding systems; we overestimate the typical masses of absorbers by a factor of ≈ 2, and the typical size by ≈ 25%. We also find that the number of identified absorbers increases with spatial resolution. For example, in the N = 64 3 run, we find just 7 absorbers, while in the N = 1024 3 run we find 429. This suggests that poorer resolution tends to blend together the fine grain structure that should in reality be separated into separate systems.

A.2. Effect of RT domains on sizes of absorption systems
Our RT domain structure was designed to avoid the difficulties in interpretation that would arise if gas parcels were ionized at different times and intensities within our boxes. However, absorbers may span multiple domains if their typical size is larger than the domain length, L dom = 32h −1 kpc. This could result in over-dense gas within the absorbers being unnaturally ionized by the intersecting source planes. In Paper I we tested the use of different domain sizes and demonstrated that the gas clumping factor does not change significantly when compared at the same global ionized fraction. Here we test specifically for any effects on the mass and size distributions of the absorbers.
For these tests we adopt L box = 256h −1 kpc and N = 256 3 , with (z re ,Γ −12 , δ/σ ) = (8, 0.3, 0). We ran two additional simulations with larger domain sizes, N dom = 4 3 and 2 3 . Figure 11 compares the mass and size distributions of selfshielding structures in these runs against the case with N dom = 8 3 (which has the same L dom = 32h −1 kpc as in our production runs). As in the last section, we adopt a constant neutral fraction threshold of 0.005 for identifying the absorbers. Overall, we find that the distributions are broadly consistent with each other. This occurs because source cells that intersect the dense core of an absorber will tend to stay optically thick. In this case the cells will be grouped as a single system despite the absorber spanning two or more domains. Note, however, that a tail appears at the low-mass end for the N dom = 4 3 and 2 3 cases. This may arise because of an increased prevalence of shadowing when  Figure 10. Numerical convergence with respect to gas/RT grid size for key quantities examined in the main text. Test simulations were run in a smaller box with L box = 256h −1 kpc and grid sizes spanned N = 64 3 to N = 1024 3 in factors of 8. We used N dom = 8 3 to match the (32h −1 kpc) 3 RT domain sizes in the productions runs. All runs used (zre ,Γ−12, δ/σ ) = (8, 0.3, 0). Moving clockwise from top-left, we show the cumulative Lyman limit absorption coefficient as a function of H I column, the median photoionization rate (in units of the background rate) vs. neutral hydrogen number density, the distribution of their effective diameters and the distribution of total masses for self-shielding systems the domains are larger. This would create some low-mass self shielding systems. We conclude from this test that the absorber mass and size distributions presented in §6.1 are likely robust to the RT domains and the particular choice of their sizes.

B. THRESHOLD FOR SELECTING SELF-SHIELDING SYSTEMS
Here we describe our procedure for choosing an H I fraction threshold, x thresh HI , to define the self-shielding systems discussed in §6. There is no ideal procedure for associating three-dimensional structures with optically thick onedimensional sight lines. For example, suppose we form groups out of cells that are above some threshold x thresh HI . Any value that we choose will exclude some sight lines that contain optically thick absorbers, and include some that do not. The essence of our procedure is to pick a threshold that minimizes such errors.
We wish to minimize the combined likelihood that our absorber definition will exclude optically thick segments (Type I errors) and include optically thin segments (Type II errors). For a given threshold HI fraction x thresh HI , we count the fraction of optically thick segments that do not intersect a system, P (Type I), and the fraction of optically thin segments that do, P (Type II)). Then we adjust x thresh HI until the sum P (Type I)+P (Type II) is minimized. Table 1 shows our values of x thresh HI obtained using this method, and the conditional probability of each type of error. The low error probabilities suggest that our method picks out the systems responsible for generating optically thick sight-lines effectively. The optimal threshold is typically in the range 0.001 -0.005, and is most sensitive to the value of the photo-ionization rate, with mild dependence on reionization redshift and time.  Table 1. From left to right, the simulation parameters, the H I threshold that minimizes the sum of type I and II errors, and the minimized error probabilities. The thresholds are calculated at 10, 60 and 300 Myr after zre for each model.