Effects of Small-scale Absorption Systems on Neutral Islands during the Late Epoch of Reionization

The reionization process is expected to be prolonged by the small-scale absorbers (SSAs) of ionizing photons, which have been seen as Lyman-limit systems in quasar absorption line observations. We use a set of semi-numerical simulations to investigate the effects of absorption systems on the reionization process, especially their impacts on the neutral islands during the late epoch of reionization (EoR). Three models are studied, i.e., the extreme case of no-SSA model with a high level of ionizing background, the moderate-SSA model with a relatively high level of ionizing background, and the dense-SSA model with a low level of ionizing background. We find that while the characteristic scale of neutral regions decreases during the early and middle stages of reionization, it stays nearly unchanged at about 10 co-moving Mpc during the late stage for the no-SSA and moderate-SSA models. However, in the case of weak ionizing background in the dense-SSA model, the characteristic island scale shows obvious evolution, as large islands break into many small ones that are slowly ionized. The evolutionary behavior of neutral islands during the late EoR thus provides us with a novel way to constrain the abundance of SSAs. We discuss the 21 cm observation with the upcoming Square Kilometre Array. The different models can be distinguished by either the 21 cm imaging or the 21 cm power spectrum measurements.


Introduction
Cosmic reionization is a significant phase transition in the global history of the universe, during which the hydrogen atoms in the intergalactic medium (IGM) are ionized by photons that have escaped from galaxies. Recent observations of the cosmic microwave background (CMB) and the highredshift Lyα absorptions indicate a relatively late completion of the reionization process (e.g., Becker et al. 2015;Planck Collaboration et al. 2020;Becker et al. 2021;Bosman et al. 2021;Qin et al. 2021a), with potential neutral islands lasting up to below redshift 6 (Kulkarni et al. 2019;Keating et al. 2020). Increasing attention has been driven to the late epoch of reionization (EoR), investigating the statistics of the large-scale neutral regions (islands; Giri et al. 2019;Nasir & D'Aloisio 2020), and how the islands evolves (Xu et al. 2014(Xu et al. , 2017, and a statistical tool has been proposed to determine the completing redshift of reionization by measuring the 21 cm bias . The 21 cm transition of neutral hydrogen (H I) has the potential to trace the whole history of cosmic dawn and reionization (see, e.g., Furlanetto et al. 2006 for a review). The Experiment to Detect the Global Epoch of Reionization Signature has detected an absorption feature in the global radio spectrum that could be the cosmic dawn signature (Bowman et al. 2018). On the other hand, various lowfrequency interferometers have put increasingly lower upper limits on the 21 cm power spectrum from the EoR (e.g., Mertens et al. 2020;Trott et al. 2020;The HERA Collaboration et al. 2022). The upcoming Square Kilometre Array (SKA) will have the ability to do tomographic imaging of the whole reionization process (Koopmans et al. 2015;Xu & Zhang 2020).
During the late EoR, when the ionized regions become overlapped and interconnected with each other, the characteristic scale of ionized regions exceeds the mean free path (MFP) of the ionizing photons. In addition to the large-scale neutral islands, the propagation of the ionizing photons is significantly affected by small-scale absorbers (SSAs; Park et al. 2016). These absorbers have also been observed in high-redshift quasar absorption lines, in particular the numerous Lyman-limit systems (LLSs), which have sufficient column density to remain neutral. We will refer to all of these as SSAs. These played a significant role in regulating the ionizing background and the reionization process (e.g., Miralda-Escudé et al. 2000;Alvarez & Abel 2012;Xu et al. 2017). Therefore, it is essential to properly model the effects of SSAs when predicting and interpreting the upcoming 21 cm data, especially for the signal from the late EoR (e.g., Sobacchi & Mesinger 2014).
However, the SSAs are currently not resolvable in largescale (>100 Mpc) reionization simulations. They are usually implemented with approximations, either by implementing a sub-grid clumping (see, e.g., Iliev et al. (2007) for fullnumerical simulations, and Sobacchi & Mesinger (2014) for semi-numerical simulations), or by introducing an MFP parameter (see, e.g., Iliev et al. 2014;Shukla et al. 2016 for full-numerical simulations, and see Alvarez & Abel 2012;Xu et al. 2017 for semi-numerical simulations). In particular, Alvarez & Abel (2012) used a semi-numerical approach to investigate the effect of the abundance of SSAs on the largescale morphology and the overall progress of reionization, by varying the MFP parameter, and found that the characteristic bubble size is quite sensitive to the MFP of ionizing photons. Sobacchi & Mesinger (2014) and Shukla et al. (2016) used different approaches, but arrived at the same conclusions that the presence of SSAs impeded the late growth of the ionized regions and suppressed the 21 cm power spectrum on large scales substantially.
These previous studies are primarily focused on the ionized regions, and the evolution of neutral islands are less discussed. One exception is Giri et al. (2019), in which a comprehensive study of the statistical properties of neutral islands was presented. They showed different topology and evolutionary features between the ionized bubbles in the early EoR and the neutral islands in the late EoR. Being fully numerical, however, it has been a challenge to explore a large parameter space, and the effect of SSAs on the island statistics is not well studied.
In the present work, we focus on the island stage, i.e., the last stage in the topological evolution of reionization (Chen et al. 2019), when the SSAs are most important in regulating the reionization process. We use semi-numerical simulations to explore a range of parameters that are consistent with a variety of current observations, i.e., the latest Thompson optical depth to the CMB inferred by the Planck data (Planck Collaboration et al. 2020), the ionization state of the IGM inferred from highredshift galaxy and quasar observations (e.g., McGreer et al. 2015;Mesinger et al. 2015;Sobacchi & Mesinger 2015;Davies et al. 2018;Mason et al. 2018Mason et al. , 2019Jung et al. 2020), the observed number density of LLSs at high redshifts (e.g., Songaila & Cowie 2010;Crighton et al. 2019), and measurements of the high-redshift H I photoionization rate (e.g., Calverley et al. 2011;Wyithe & Bolton 2011;D'Aloisio et al. 2018). We investigate the effect of SSAs on the morphology and size evolution of the neutral islands, discuss the degeneracies between the abundance of SSAs and the properties of ionizing sources, and potential ways to break them, and explore the observational strategy to distinguish the different models in upcoming 21 cm tomographic observations. Various techniques have been developed to measure the MFP of ionizing photons, such as using the incidence rate of LLSs (e.g., Meiksin & Madau 1993;Songaila & Cowie 2010), or directly measuring the opacity from the mean spectra of high-redshift quasars (e.g., Prochaska et al. 2009;Worseck et al. 2014), or averaging the free paths along multiple lines of sight against quasars (Romano et al. 2019). With the increasing sample of high-redshift quasars, measurements of the abundance of LLSs and the MFP have effectively constrained the reionization process (e.g., Cain et al. 2021;Davies et al. 2021). However, it is still challenging to measure precisely the MFP in the reionization epoch. Here we propose that the 21 cm observations can be a novel probe to constrain the abundance of SSAs and the MFP of ionizing photons during the EoR. This paper is organized as follows. We briefly describe the semi-numerical simulations used in our study and the models for SSAs in Section 2. Section 3 contains the theoretical results on the effects of SSAs on the neutral islands, in terms of the size distribution and evolution of islands, and the power spectrum of the 21 cm brightness emitted by the islands. We also discuss the long troughs observed in quasar absorption spectra. In Section 4, we compare the 21 cm imaging observation and the power spectrum measurement, and discuss the prospects of distinguishing different reionization models with the low-frequency array of phase one SKA (SKA1-Low). We conclude in Section 5. Throughout this paper, we assume the standard Λ cold dark matter model with the cosmological parameters of Ω Λ = 0.685, Ω m = 0.315, Ω b = 0.0495, H 0 = 67.4 km s −1 Mpc −1 , σ 8 = 0.811, and n s = 0.965.

The Semi-numerical Simulations
In order to investigate the effect of SSAs on the neutral islands, we use the semi-numerical simulation islandFAST (Xu et al. 2017), which is based on the island model (Xu et al. 2014) and was developed for the isolated islands topology during the late EoR. Being semi-numerical, it is more efficient in exploring the parameter space, so that we can study the effects within the parameter space constrained by current observations.
In the framework of the excursion-set theory (Bond et al. 1991;Lacey & Cole 1993;Furlanetto et al. 2004), the islandFAST filters the evolved density fields over increasingly smaller scales, and determines each cell in the simulation box as ionized (or neutral) if the number of ionizing photons available within the filtering scale, corrected for recombinations, exceeds (or drops below) the number of hydrogen atoms in this region. This ionization/neutral criterion is called a barrier in the excursion-set theory. For the early EoR until the percolation of ionized regions, the islandFAST is identical to the original version of 21cmFAST (Mesinger et al. 2011), in which the inhomogeneous recombination (Sobacchi & Mesinger 2014) and mini-halos (Qin et al. 2020(Qin et al. , 2021b are not yet incorporated, so as to be consistent with the homogeneous ionizing background implemented for the late EoR. Assuming that the number of the ionizing photons emitted is proportional to the total collapse fraction of the region, the ionization condition, i.e., the bubble barrier can be written as H r e c 1 is the ionizing efficiency parameter, in which f esc , f å , N γ/H , andn rec are the escape fraction, star formation efficiency, the number of ionizing photons emitted per H atom in stars, and the average number of recombinations per ionized hydrogen atom, respectively.
For the island stage after percolation, when the neutral islands are isolated, the islandFAST adopts an inverse topology, and the condition for a region of mass scale M to keep from being fully ionized, i.e., the island barrier, is written as where δ M is the mean overdensity of the region under consideration, and X H is the mass fraction of baryons in hydrogen. Here N back is the number of background ionizing photons consumed by the region under consideration, which can be written as wheren b is the mean co-moving baryon number density, R i is the initial island scale corresponding to the mass scale M when the ionizing background is just set up, and R f is final scale of the island at the redshift under consideration. The change in the island scale is due to the ionization by the ionizing background; thus, we have where F(z) is the physical photon number flux of the ionizing background, z back is the redshift at which a global ionizing background is set up, and H(z) is the Hubble parameter. The number density of background ionizing photons can be written as (Xu et al. 2017) is the physical distance between the source at ¢ . Note that in islandFAST, the MFP is implemented by its original definition of attenuating ionizing photons, but not a sharp boundary for the filtering scale, or a sharp cut-off for the traveling distance of ionizing photons, as in previous simulations (e.g., Alvarez & Abel 2012;Iliev et al. 2014;Shukla et al. 2016). In this work, ξ is defined in terms of N γ/H , and we assume the same production rate of ionizing photons per He atom in stars as that per H atom in stars, i.e., N γ/He = N γ/H , so Equation (2) implicitly includes the contribution and consumption of ionizing photons by helium atoms, and the factor m H /X H in the second term in Equation (2) cancels out the same factor when computingn H in Equation (5).
For the island stage, the islandFAST takes a two-step filtering approach; it first finds host islands by identifying firstdown-crossings with respect to the island barrier including the contribution from an ionizing background, and then identifies bubbles-in-island by applying a second filtering step within each host island using the bubble barrier. With an input model of SSAs, the islandFAST self-consistently computes the effective MFP limited by both the large-scale neutral islands and the SSAs, and derives the evolution of the ionizing background simultaneously with the ionization field in adaptive redshift steps.
In the original version of islandFAST (Xu et al. 2017), the ionizing background was assumed to be present only in the island stage. However, the percolation process begins as early as when the universe was only ∼30% ionized (Furlanetto & Oh 2016;Chen et al. 2019), and the ionizing background is expected to have been gradually set up since then, though it may have fluctuated significantly during the percolation process. In this work, we assume that a relatively uniform ionizing background has been present since the "neutral fibers" stage (i.e., when the neutral regions thinned into filamentary structures; Chen et al. 2019), and z back is set by the time when the mean neutral fraction is 30%.
We incorporate three different models for the SSAs, as detailed below. As shown in previous works, more abundant SSAs result in a more delayed and prolonged reionization process (e.g., Alvarez & Abel 2012;Sobacchi & Mesinger 2014). In addition to the SSA models, two key parameters determine the process of reionization, i.e., the minimum mass M min , or equivalently the minimum virial temperature T min vir , of halos that contribute ionizing photons, and the ionizing efficiency parameter ξ. The former regulates the time at which the reionization begins in earnest, while the latter determines the speed of the reionization process. By tuning these two parameters, we control the reionization processes with different SSA models to be consistent with various observations.

The Models for SSAs
During the island stage, the MFP of ionizing photons is limited not only by large-scale underdense islands but also by small-scale overdense absorbers. The effective MFP of ionizing photons λ mfp is given by where λ I is the MFP limited by large-scale underdense islands, and λ abs is the MFP due to small-scale overdense absorbers. In the simulation, λ I is computed on-the-fly using the MFP algorithm (Mesinger & Furlanetto 2007), which is taken as the average length of random vectors starting from an ionized pixel reaching an edge of a neutral island. Because of the relatively low column density of Lyα forest systems and the relative rarity of damped Lyman-alpha systems, the LLSs are the dominant contributor to the IGM opacity, reducing the MFP of ionizing photons, and weakening the ionizing background Shukla et al. 2016). Here the parameter λ abs implicitly includes the contributions from all kinds of unresolved absorbers, i.e., the SSAs. We consider the following three models of SSAs corresponding to different levels of the ionizing background.
islandFAST-noSSA: In the extreme case of no SSA, the MFP of ionizing photons is completely determined by the morphology of large-scale neutral islands, i.e., λ mfp = λ I , the ionization background is large in this case. Though unrealistic, we can use this case for comparison. islandFAST-SC: The number density of LLSs was observed up to redshift ∼6 (Songaila & Cowie 2010). In the second model, we adopt the evolution of the MFP limited by SSAs extrapolated from the fitting formula provided by Songaila & Cowie (2010), expressed as: This model corresponds to a relatively high level of the ionizing background, or equivalently a moderate MFP.
islandFAST-RS: Recent observations have favored a rapid evolution of the MFP near the end of reionization (e.g., Becker et al. 2021), implying significant modulation of the MFP by the large-scale ionization field. The island-FAST-RS model adopts the neutral-fraction-dependent λ mfp given by the high-resolution Aurora radiationhydrodynamical simulations for the EoR (Figure 1 in Rahmati & Schaye 2018). We interpolate to calculate the MFP for a derived neutral fraction from our simulation, and iterate to achieve a consistent ionization field and the MFP for each redshift step. As shown below, this model corresponds to a low level of the ionizing background and a short MFP.
Our simulation box has a co-moving size of 1 Gpc a side and 600 3 cells. For each of the SSA models, we run two sets of simulations corresponding to late and high-mass ionizing sources. As a reference, we also run the 21cmFAST for the same sets of parameters. In the following analyses, the fiducial source parameters for 21cmFAST, islandFAST-noSSA, and islandFAST-SC models are , chosen to be consistent with various observations. The evolutions of the volume neutral fractionx H I from the three islandFAST models, and that from 21cmFAST, are shown in Figure 1. Also shown are the observational constraints from various observations as indicated in the legend and the caption. The reionization processes predicted by the 21cmFAST, islandFAST-noSSA, islandFAST-SC, and island-FAST-RS end at redshifts 6.02, 6.00, 5.74, and 5.30, respectively. The corresponding Thompson optical depths are 0.0545, 0.0543, 0.0541, and 0.0557, respectively, all consistent with the latest Planck results (Planck Collaboration et al. 2020). We note that the reionization history of the late EoR based on island topology does not connect smoothly with the history of early EoR simulated assuming bubble-topology. This is an intrinsic defect of the simulation in describing the complex percolation process. However, as will be shown below (in Section 3.1), when compared at the same mean neutral fractions, the islandFAST-noSSA and 21cmFAST predict consistent island statistics, so this discontinuity at the topological transition would not affect the results on the island statistics. We reserve the improvement on the simulation strategy for the mid-EoR to future works. Here we confirm the previous result that the presence of SSAs prolongs the reionization process. Figure 2 shows the evolution of the MFP for ionizing photons (left panel) and the photoionization rate (right panel) predicted by the three SSA models. In the case of no SSA, the MFP is entirely limited by neutral islands, and as expected, the MFP grows rapidly as the islands being ionized (blue solid line) and the ionizing background also increase rapidly. In the presence of SSAs, the propagation of ionizing photons is effectively blocked, resulting in a significant reduction of the  MFP and the ionizing background, and their growth rate also slows down (green solid line). The MFP gradually approaches the value completely limited by the SSAs (cyan dashed line). If we further enhance the absorption of SSAs to ionizing photons, as in the islandFAST-RS model, the MFP and the ionizing background are further reduced, with only a weak evolution (black solid line). In this case, the reionization is significantly delayed by SSAs, and the island stage lasts for a long time. The magenta dashed line is extrapolated from the direct measurement of the MFP using the drop in the continuum flux of highresolution quasar spectra up to z ∼ 5 (Worseck et al. 2014). It is seen that the islandFAST-SC model is in the best agreement with the MFP measurements by Songaila & Cowie (2010), Worseck et al. (2014), and the data point at z ∼ 5 by Becker et al. (2021), while the islandFAST-RS model is more consistent with the photoionization rate measurements and the recent MFP measurement at z ∼ 6 by Becker et al. (2021). Note that the observational data for the MFP and the ionizing background are still very limited at high redshifts, and the measurements could be biased due to the proximity effect (D'Aloisio et al. 2018;Davies 2020;Becker et al. 2021). There is still large uncertainty in the MFP and the ionizing background beyond the reionization. The islandFAST-SC model and the islandFAST-RS model are both of great significance, bracketing possible scenarios consistent with current observations.

Effects of Absorption Systems
Slices of the ionization field from the different seminumerical simulations are shown in Figure 3. The top, middle, and bottom panels are for islandFAST-noSSA, island-FAST-SC, and islandFAST-RS models, respectively. All slices are 1 co-moving Gpc on a side and 1.67 co-moving Mpc thick, and the three columns correspond to mean neutral fractions of 0.16, 0.10, and 0.01, from left to right, respectively. The neutral islands are shown as black patches, and the ionized regions are left white. It is worth noting that most of the neutral islands are "porous," showing the "bubbles-in-island" effect. Generally, there is no significant morphological difference between the islandFAST-noSSA model and the island-FAST-SC model during the late stage of reionization (¯ x 0.1 H I ), and the remaining islands predicted by them are small in number but relatively large in scale. However, when the absorption of ionizing photons is significant enough, as in the islandFAST-RS model, a larger number of small islands are left in the universe. Different from the other two models, the topology of the neutral islands predicted by island-FAST-RS model changes significantly between¯= x 0.16 H I and 0.10; the large islands break up into more small islands. This implies that, with the number density observation of SSAs up to higher redshifts in the future, we may be able to distinguish the reionization models with different morphology, according to the anticorrelation between the number density of SSAs and the typical size of neutral islands.

The Size Distribution of Neutral Islands
In this section, we explore the effects of absorption systems on the size distribution and its evolution of neutral islands. We characterize the scale of islands with the MFP algorithm (Mesinger & Furlanetto 2007). Note here the differences between the MFP algorithm in characterizing the sizes of neural/ionized regions as the typical length of random vectors reaching a phase transition in the ionization field, and the MFP for ionizing photons to travel freely before being absorbed. As a reference model, we first present the size distributions of neutral regions from the 21cmFAST simulation in Figure 4. Note that for the early stages of reionization, all of the islandFAST models are identical to the 21cmFAST, before adopting the different topology for the island stage. For the island stage, the main difference between the original version of 21cmFAST and the islandFAST-noSSA model is the different assumption of the ionization topology and the different filtering algorithm, but both of them assume no additional absorption systems.
The left panel of Figure 4 shows the evolution of the size distribution adopting the parameters of . The size distribution of neutral regions evolves during the early and mid-stages of EoR; however, when the universe enters the neutral fibers stage (¯ x 0.3 H I ; Chen et al. 2019), the evolution becomes extremely slow, and the characteristic scale stalls at about 10 co-moving Mpc. This is consistent with the intuition that when measuring the size of neutral regions with the MFP algorithm, most random vectors would not be parallel with the stretching direction of neutral fibers, and the derived size would be similar to the neutral islands that the neutral fibers fragment into later. However, the size distribution will show a more significant evolution from the neutral fiber stage to the island stage if one characterizes the size using the Friend-of-Friend algorithm, as illustrated in Giri et al. (2019).
The during the late EoR. We find that the characteristic scale of neutral regions does depend on the properties of ionizing sources. More massive ionizing sources with higher ionizing efficiencies will lead to larger islands during the late stage, which is consistent with Giri et al. (2019).
In Figure 5, we present the size distributions of the neutral islands for the three SSA models, during the island stage of reionization. In each panel, the blue and green sets of curves correspond to source parameters of { } x = = T 10, 10 K , respectively. Though with the even faster reionization process, the evolution and source dependence of the island size distribution from the island-FAST-noSSA model (left panel of Figure 5) are similar to the behavior predicted by 21cmFAST, as shown in Figure 4. We note that the source dependence of the island size distribution is similar among the different SSA models, with the more massive ionizing sources resulting in larger typical sizes of islands. However, the different SSA models predict distinct evolutionary behavior of the island size distribution. We find that the characteristic island scale is almost unchanged and stalls at ∼10 co-moving Mpc in the islandFAST-noSSA model and the islandFAST-SC model, which have a strong ionizing background or a long MFP. However, it shows a clear evolution in the islandFAST-RS model with weak ionizing background or a short MFP. This reflects the competitive roles played by the ionizing background and the ionizing sources inside the islands. When the ionizing background dominates the ionizing budget, the islands will be rapidly ionized from the outside, and it would be hard for the small islands to keep their identities. However, if the inside ionizing sources dominate the  ionizing photons, the islands tend to disintegrate from the inside, which results in the more obvious evolution in the island size and a more prolonged reionization process because of the weaker ionizing background. We compare the size distributions of neutral islands between the three models at the same neutral fractions in Figure 6. The blue solid, green dashed, and purple dotted-dashed curves are for the islandFAST-noSSA, islandFAST-SC, and islandFAST-RS models, respectively, with source parameters of , showing the degeneracy between the source properties and the absorber abundances. Although the effect of lower-mass sources with lower ionizing efficiency could mimic the effect of having more abundant SSAs at a certain redshift, this degeneracy can be easily broken by observing multiple redshifts and making use of the evolutionary behavior of the island size distribution.
We notice that there is no significant difference in the size distribution between the islandFAST-noSSA model and the islandFAST-SC model; although, the level of the ionizing background differs significantly as shown in Figure 2, indicating that the characteristic island scale is not sensitive to the abundance of the SSAs, at least below a certain abundance threshold. In order to estimate this threshold, we perform several islandFAST-SC simulations with the fiducial parameters of { } x = = T 10, 10 K min vir 4 , varying the prefactor in Equation (7); λ abs is taken as 3, 0.8, 0.6, 0.4, and 0.2 times of the fiducial islandFAST-SC model, respectively. Then we measure the characteristic island scale, at which the size distribution peaks, at¯= x 0.01 H I for each simulation. Figure 7 shows the relation between the characteristic island scale and λ mfp . It is seen that as the abundance of SSAs decreases, or the MFP increases, the characteristic island size first increases, and then gradually saturates at about 9 comoving Mpc when λ mfp  10 physical Mpc. The SSA abundance in the fiducial islandFAST-SC model roughly represents the abundance threshold in this case. The horizontal dashed line indicates the scale 15% smaller than that predicted by the islandFAST-noSSA model, and the corresponding λ mfp is about 6.5 physical Mpc (vertical dashed line). Note that at the very late stage (¯ x 0.01 H I ), the SSAs dominate the opacity of the IGM, and λ abs can be approximated as λ mfp , then the SSA abundance can be effectively constrained by measuring the typical scale of the remaining neutral islands.

Long Troughs
Becker et al. (2015) observed an extremely long (∼150 comoving Mpc) opaque Gunn-Peterson trough extending down to redshift 5.5 in the spectrum of quasar ULAS J0148+0600. Several works have proposed that there could be large islands that persist until the late EoR, giving rise to a non-negligible possibility of lines of sight passing through these large islands, which results in the long troughs as observed (Giri et al. 2019;Kulkarni et al. 2019;Keating et al. 2020). Here we also calculate the probability that the observed long absorption trough is caused by the neutral islands in the late EoR. Figure 8 shows the probability for the length of a sightline intersected by an island to be larger than 50 co-moving Mpc, as a function of the mean neutral fraction, from the simulated three models. This represents the lower limit to the islands' contribution to the Lyα absorption trough, where the neutral islands are the only contributor to the long troughs. We find that for all three models, the possibility of having an island scale larger than 50 co-moving Mpc is quite low (of the order of ∼10 −3 ) for the island stage, but the probability increases rapidly toward the neutral fiber and earlier stages. This implies that if the long trough is solely caused by a large neutral region, the reionization has to complete quite late. However, our model does not include partial ionizations, and the bubbles-in-islands significantly reduce the length scale of the intersected lines of sight passing through the porous neutral islands. The residual neutral hydrogen in the bubbles-in-islands may connect the neighboring absorption troughs giving rise to longer troughs. In addition, the damping wing absorption of neutral islands can produce troughs longer than the islands themselves. On the other hand, the residual neutral hydrogen in the ionized regions could also contribute to the Lyα opacity, and may significantly enlarge the absorption troughs. The opacity of the ionized IGM is closely related to the residual neutral fraction of gas at low density   . Therefore, a lower temperature or a lower photoionization rate will lead to the increase of the residual neutral fraction. The large-scale overdense regions are likely to be cooler than average, because they were reionized early (D'Aloisio et al. 2015), and large-scale voids are likely to have a weaker ionizing background (Davies & Furlanetto 2016). These regions could both be mostly ionized, but still have a residual neutral fraction high enough to saturate the absorption of Lyα photons (Keating et al. 2020). We plan to investigate these possibilities by incorporating fluctuations in the ionizing background into future works. Here we calculate the probability of a random 50 co-moving Mpc sightline intersecting an island of any size, and the results are shown in Figure 9. These curves put upper limits on the contribution of . The islandFAST-noSSA, islandFAST-SC, and islandFAST-RS models are shown with a blue dot, green pentagon, and black hexagon, respectively. The horizontal dashed line represents a scale 15% smaller than the characteristic island scale predicted by islandFAST-noSSA, and the vertical dashed line denotes the corresponding λ mfp .  , which result in larger and fewer islands. But thereafter, as the large islands break into many small ones, it provides the highest upper limits of the probability.

The 21 cm Power Spectrum
The effect of SSAs on the morphology of neutral islands may also be reflected in the 21 cm power spectrum, hopefully to be constrained in the near future. The differential 21 cm brightness temperature of hydrogen against the CMB at redshift z, corresponding to the observed frequency ν, can be written as (e.g., Mesinger et al. 2011 is the evolved density contrast, dv r /dr is the co-moving gradient of the line-of-sight component of the co-moving velocity, and T γ and T S are the CMB brightness temperature and the spin temperature of the neutral hydrogen gas, respectively. All of these quantities are evaluated at redshift z = ν 0 /ν − 1, with ν 0 = 1420.4 MHz being the rest-frame frequency of the 21 cm line. During the late EoR that we focus on here, the X-rays emitted by the early X-ray binaries have probably heated the IGM temperature to a level much higher than the CMB temperature (Pritchard & Furlanetto 2007;Pritchard & Loeb 2012), and in combination of the efficient Lyα coupling between the spin temperature and the gas kinetic temperature, they ensure that T S ? T γ . We will assume T S ? T γ in the remainder of the paper.
We compute the 21 cm power spectrum with statistical errors for the three SSA models, and compare them at the same neutral fractions as in Figure 10. The statistical errors are computed by , where N is the number of independent modes in each k-bin. On the small-scale end, as the reionization proceeds to the end, the neutral fibers fragment into islands, and the bubbles-in-island effect becomes less and less prominent, and the small-scale fluctuations tend to be suppressed for all three models. However, the models differ significantly on large scales. Comparing the colored curves for the three models with the same source parameters of {ξ = 10, = T 10 K min vir 4 }, it is obviously seen that the large-scale 21 cm power spectrum (k < 0.1 Mpc −1 ) decreases faster in models with more abundant SSAs, which is in consistent with what was found in Shukla et al. (2016). In the very late stage (¯ x 0.10 H I ), the islandFAST-noSSA model predicts the highest fluctuations of 21 cm signal on large scales, because of the larger and fewer islands as shown in the upper panels in Figure 3. The more abundant SSAs, or the lower level of the ionizing background, result in the larger number of small islands, significantly suppressing the large-scale power in the 21 cm fluctuations, which could potentially be tested with the upcoming 21 cm power spectrum measurements. The black dotted-dashed lines in Figure 10 shows the 21 cm power spectrum for the islandFAST-RS model with { } x = =T 30, 5 10 K min vir 4 . As compared with the purple dotted-dashed lines for the same model but different source parameters, they have a higher amplitude on large scales at all redshifts because of the larger islands resulting from the massive ionizing sources with higher ionizing efficiencies. Similar to the situation in the island size distribution, the effect of SSAs is degenerate with the effect of source properties for a single redshift, but the degeneracy can be broken by observing the evolution of the 21 cm power spectrum on large scales.

21 cm Observations
In the above, we have seen that the SSAs have the effects of prolonging the reionization process, and changing the morphology of the neutral islands during the late stages of reionization, which are reflected in both the size distribution of islands and the 21 cm power spectrum. The 21 cm power spectrum could potentially be detected with the current generation of low-frequency interferometers, such as the Low Frequency Array (van Haarlem et al. 2013), the Murchison Widefield Array (Tingay et al. 2013), and the Hydrogen Epoch of Reionization Array (DeBoer et al. 2017), and will be more precisely measured by the upcoming SKA. The neutral regions may also be directly imaged by the SKA (Koopmans et al. 2015). Here we focus on the SKA, and investigate the possibilities of distinguishing the different SSA models with the power spectrum and images.
Specifically, we assume the phase one configuration of the low-frequency array of SKA (SKA1-Low), which includes 512 . antennae stations, each has an effective diameter of 38 m. We consider the core array with 224 stations within a radius of 500 m, and the inner array of 278 stations distributed in an area of radius 1700 m. We will take into account the limited angular and spectral resolution and the expected thermal noise of the array, but assume that the foregrounds and radio frequency interference removal and calibrations are perfect, though these effects would depend on technologies that are still under development (e.g., Barry et al. 2016).
To simulate the slices of the observed δT b field, we first smooth the δT b mock field to match the telescope resolution. The co-moving size of a resolved pixel perpendicular to the line of sight is where D c (z) is the co-moving angular distance to redshift z, and q l D~L 1.22 max is the angular resolution of the array, in which λ and L max are the observing wavelength and the maximum baseline, respectively. The pixel size along the line of sight is related to the observing bandwidth B by where c is the speed of light. In this study, we try B = 0.1 MHz and B = 1 MHz, corresponding to l z ∼ 1.48 Mpc and l z ∼ 14.8 Mpc at z = 6, respectively. When observing with the narrow band of 0.1 MHz, we use a single slice of depth 1.67 Mpc directly from the simulation box, and when observing with the broad band of 1 MHz, we average the adjacent nine slices to get a two-dimensional slice for each redshift z. Then each slice is convolved with a Gaussian window of width s = l 8 ln 2 xy . Next, we add the thermal noise to the smoothed slices. For the core array and the inner array of SKA1-Low considered in this work, we simply assume a uniform uv-coverage; then, the Gaussian thermal noise inside a resolution element corresponding to a scale k ⊥ can be written as (Koopmans et al. 2015): , 11 where Ω FoV is the field of view of the array, T sys is the system temperature, and t int is the integration time. A eff = π(19 m) 2 is the effective collecting area of each station, which is related to the field of view by Ω FoV = λ 2 /A eff . Here A core is the area over which the antennae are distributed, and the total collecting area is A coll . For the core array of SKA1-Low, ( ) p = A 500 m core 2 , A coll = 224 × A eff , and q D~¢ 6.17 at z = 6 (corresponding to 15.2 co-moving Mpc). For the inner array, we have ( ) p = A 1700 m core 2 , A coll = 278 × A eff , and q D~¢ 1.81 at z = 6 (corresponding to 4.45 co-moving Mpc). At the frequency range relevant to the 21 cm signals from the EoR, the system temperature is dominated by the sky brightness temperature, and we assume T sys = 100 + 400 × (ν /150 MHz) −2.55 K (Koopmans et al. 2015). In the following, we adopt deep surveys with t int = 1000 hr. Then the core array has a noise level of σ T ∼ 1.4 mK or 0.43 mK for B = 0.1 MHz or 1 MHz, respectively, and the inner array will have σ T ∼ 12.7 mK (B = 0.1 MHz) or 4.0 mK (B = 1 MHz). It is found that the increased angular resolution comes at a price, the inner array has a high level of thermal noise as compared to the cosmological signal, and it would be difficult to extract the intrinsic properties of the neutral islands. Therefore, we will use the core array in the following analysis. Figure 11 shows the δT b slices from a single slice from the simulation (top-left plot), those averaged from nine adjacent slices (top-right plot), and their corresponding mock images as observed by the SKA1-Low core array using B = 0.1 MHz (bottom-left plot) or B = 1 MHz (bottom-right plot), respectively. In each plot, we show the slices predicted by the islandFAST-noSSA (top row), islandFAST-SC (middle row), and islandFAST-RS models (bottom row) at fixed neutral fractions of¯= x 0.16 H I , 0.10, and 0.01 from left to right, with the corresponding redshifts marked on the slices. The neutral islands are seen in 21 cm emission, while the ionized regions fluctuate around the zero-brightness in the mock images due to the thermal noise.
It is seen from the top-left plot that the intrinsic δT b slices perfectly follow the ionization field ( Figure 3) as expected, while the averaged δT b slices in the top-right plot show clearly the projection effect. Most small bubbles within the islands disappear, and the boundaries of the islands expand due to the line-of-sight smoothing. With the limited angular resolution of the SKA1-Low core array (15.2 co-moving Mpc in the transverse direction at z = 6), δT b near the boundaries of the islands is lowered, while δT b near the edges of ionized regions becomes slightly positive, making the boundaries slurred and enlarging the size of the observed islands, as shown in the bottom plots of Figure 11.
In order to calculate the size distribution of neutral islands that is possibly extractable from two-dimensional δT b slices, even without any instrumental effects, we apply the MFP method to both the single δT b slices directly from simulations and the averaged slices of depth 15 Mpc, and the results are shown in the top and bottom panels of Figure 12, respectively. In each row, the three panels from left to right are for the islandFAST-noSSA, islandFAST-SC, and island-FAST-RS models, respectively. We see that the characteristic island scales extracted from single δT b slices (top panels) can well reflect the typical sizes from the three-dimensional ionization fields as shown in Figure 5, which we refer to as "intrinsic scale" below. However, the characteristic scales measured from the averaged δT b slices (bottom panels) are much larger than the intrinsic scale. Moreover, when projected onto the two-dimensional δT b slices, the evolutionary behavior of the size distribution is somewhat distorted. The island-FAST-noSSA and islandFAST-SC models show an evolutionary trend from¯= x 0.10 We then apply the MFP algorithm to the mock images as observed by the SKA1-Low core array, with a brightness temperature threshold of T c . We have experimented with different T c as compared to the noise level of σ T , and found that a relatively high threshold of brightness temperature is required in order to reconstruct the intrinsic island scales from narrowband imaging, while a relatively low brightness threshold can be used to extract the size distributions of the averaged δT b fields from broadband imaging. Figure 13 shows the size distributions with brightness temperature threshold of 6 σ T for narrowband imaging (top panels) and those with 3 σ T for broadband imaging (bottom panels). It is found that using narrowband imaging with B = 0.1 MHz, we will be able to extract the intrinsic island scale using a brightness temperature threshold of 6 σ T . A lower threshold will overestimate the island scales. Using the broadband imaging with B = 1 MHz, almost all of the noises would be eliminated with a threshold value of 3 σ T . The reconstructed distributions are very similar to the size distributions shown in the bottom panels in Figure 12, though the extracted characteristic scales are a bit Figure 11. The δT b slices (top plots) and their mock images as observed by the SKA1-Low core array (bottom plots). All slices are 1 Gpc on a side in co-moving scale. The slices in the left column are 1.67 Mpc thick, and the observing bandwidth is 0.1 MHz, and the slices in the right column are 15 Mpc thick, corresponding to an observing bandwidth of 1 MHz. In each plot, the top, middle, and bottom panels are for the islandFAST-noSSA, islandFAST-SC, and islandFAST-RS models, respectively, and the three columns correspond to reionization stages with mean neutral fractions of¯= x 0.16 H I , 0.10, and 0.01, from left to right, respectively.
larger. It is clear that only the islandFAST-RS model shows obvious evolution in the extracted island scale throughout the late EoR.
From the above analysis, we see that the line-of-sight and transverse smoothing, and noise jointly affect the observed size distribution and evolution of islands in 21 cm imaging. One needs a balance between the resolution to resolve islands, in   only the islandFAST-RS model shows an obvious evolution in the measured island scale. However, the simple MFP method with a single threshold may not extract all useful information about the distributions, and the result may be affected by the choice of threshold. It may require applying multiple thresholds to extract this information. There may also be other more sophisticated and robust ways to extract the bubble or island characteristic scales from observation, e.g., the "granulometry" technique (Kakiichi et al. 2017). Such methods should be extensively developed and tested with simulation to ensure correct interpretation of the results.
The different models can also be distinguished with power spectrum measurements. In Section 3.3 we noted that the presence of SSAs mainly affects the large-scale power in the 21 cm brightness temperature fluctuations, which is our primary interest. For precise measurement of the large-scale power, the core array of SKA1-Low is suitable. Under the assumption of uniform uv-coverage, the measurement error on the 21 cm power spectrum due to thermal noise is (Koopmans et al. 2015): where l z is the depth of the survey for a bandwidth B given by Equation (10), and N b is the number of independent beams. Here we assume N b = 1, B = 1 MHz, and t int = 1000 hr. The total error is obtained by adding the statistical error to this thermal noise. We compute the 21 cm power spectra of the three SSA models for the late EoR from the δT b -field, smoothed to match the resolution of the core array of SKA1-Low. Figure 14 shows the measured 21 cm power spectra for three models, with the total measurement errors indicated by the shaded regions. We find that the differences between the three SSA models at the island stage (¯ x 0.10 H I ) are significant with respect to the expected measurement errors for the core array, so at least in principle, the upcoming SKA1-Low survey should also be able to discriminate the EoR models with power spectrum measurements. In our present case, the abundance of the small-scale absorbers and the level of the ionizing background can be distinguished or constrained.
However, we note that the 21 cm power spectrum may be biased or distorted by residual foregrounds, though these are not included in the present work. Also, these spectra are quite featureless, so they may also suffer from parameter degeneracies, e.g., the degeneracy between the SSA abundance and the source properties as discussed in Section 3.3. The imaging observations, which can provide morphological properties of neutral islands and ionized regions, are complementary in this aspect.

Conclusion
In this work, we use a set of semi-numerical simulations to study the effects of the small-scale overdense absorbers on the large-scale morphology of the underdense neutral islands during the late EoR. We consider three SSA models, i.e., the islandFAST-noSSA model with no SSA and an extremely high ionizing background, the islandFAST-SC model with a moderate number density of SSAs and a relatively high ionizing background, and the islandFAST-RS model with the most abundant SSAs and a correspondingly low level of ionizing background. In agreement with previous works (e.g., Alvarez & Abel 2012;Shukla et al. 2016), we find that the presence of SSAs prolongs the reionization process, and affects the morphology of the ionization field, which results in the suppression of the large-scale power in the 21 cm power spectrum.
The morphology of the islands reflects the competitive roles played by the ionizing photons generated inside the islands and those coming from outside. In the islandFAST-noSSA and islandFAST-SC models, the ionization is dominated by the photons from outside, i.e., the ionizing background. Small islands could hardly survive; hence, when the universe enters the neutral fibers stage, the evolution of the typical island scale is very slow and stalls at ∼10 co-moving Mpc. In the case of the islandFAST-RS model, where the SSA abundance is high and the ionizing background is weak, sources from inside dominate the reionization, and the neutral islands will more easily break into small ones. The size distribution then evolves, and the typical size of the islands is smaller. , 0.10, and 0.01, from left to right, respectively. The blue solid, green dashed, and black dotted-dashed lines are for the islandFAST-noSSA, islandFAST-SC, and islandFAST-RS simulations, respectively, with the shaded regions showing the measurement errors expected for the core array of SKA1-Low.