Crucial Factors for Ly{\alpha} Transmission in the Reionizing Intergalactic Medium: Infall Motion, HII Bubble Size, and Self-shielded Systems

Using the CoDa II simulation, we study the Ly$\alpha$ transmissivity of the intergalactic medium (IGM) during reionization. At $z>6$, a typical galaxy without an active galactic nucleus fails to form a proximity zone around itself due to the overdensity of the surrounding IGM. The gravitational infall motion in the IGM makes the resonance absorption extend to the red side of Ly$\alpha$, suppressing the transmission up to roughly the circular velocity of the galaxy. In some sight lines, an optically thin blob generated by a supernova in a neighboring galaxy results in a peak feature, which can be mistaken for a blue peak. Redward of the resonance absorption, the damping-wing opacity correlates with the global IGM neutral fraction and the UV magnitude of the source galaxy. Brighter galaxies tend to suffer lower opacity because they tend to reside in larger HII regions, and the surrounding IGM transmits redder photons, which are less susceptible to attenuation, owing to stronger infall velocity. The HII regions are highly nonspherical, causing both sight-line-to-sight-line and galaxy-to-galaxy variation in opacity. Also, self-shielded systems within HII regions strongly attenuate the emission for certain sight lines. All these factors add to the transmissivity variation, requiring a large sample size to constrain the average transmission. The variation is largest for fainter galaxies at higher redshift. The 68\% range of the transmissivity is similar to or greater than the median for galaxies with $M_{\rm UV}\ge-21$ at $z\ge7$, implying that more than a hundred galaxies would be needed to measure the transmission to 10\% accuracy.


Introduction
How was the intergalactic medium (IGM) reionized at z 6? This question is connected to the unknown properties of early galaxies and is, therefore, one of the key questions of modern astronomy. Cosmic microwave background radiation observed by the Planck satellite constrains the midpoint of reionization to be at z ∼ 8 (Planck Collaboration et al. 2020) 1 . Another powerful probe of cosmic history is the Lyα emission from high-z star-forming galaxies. Those galaxies are expected to generate strong emission near Lyα, which the neutral gas in the IGM can easily scatter. Thus, the statistics of the Lyα line strength is closely related to the global neutral fraction of the IGM (e.g., Malhotra & Rhoads 2004;Dijkstra 2014).
In order to constrain the reionization history from observed Lyα emission statistics precisely, one crucially needs to understand the IGM opacity and how it evolves during reionization. Transmission of Lyα photons in the IGM depends on the neutral hydrogen density, velocity, and temperature of the intergalactic gas. Ultraviolet (UV) radiation from starforming galaxies forms expanding HII regions around them, where some Lyα emission can be transmitted owing to low HI density. Thus, capturing the statistics of these HII regions, which can be up to tens of comoving megaparsecs (cMpc) large, is a key to the modeling the IGM transmissivity during the reionization era (e.g., Furlanetto & Oh 2005). Additionally, small-scale velocity and density structures below the megaparsec scale around the source galaxy are also important factors for the IGM transmissivity (Bolton & Haehnelt 2013;Mesinger et al. 2015;Kakiichi et al. 2016).
The IGM opacity of Lyα photons can be classified into two components: resonance and damping-wing opacity. At the line resonance, the cross section is so high that even a tiny trace of neutral fraction in highly ionized gas (well below 1%) can make the IGM opaque to Lyα photons. Thus, photons on the blue side of Lyα typically face high opacity at some point in the rest frame of the IGM and are scattered off the line of sight. On the other hand, photons on the red side of the resonance are much more likely to be transmitted. The IGM opacity for these photons mostly comes from the damping-wing cross section, which extends far away from the resonance. This damping-wing cross section is much smaller than the resonance cross section, but the HI density in neutral regions can be high enough to create a substantial optical depth. Therefore, this damping-wing opacity is expected to depend on the global neutral fraction of the IGM during reionization.
Our goal in this study is to clarify what aspects of the IGM density/velocity field modulate the resonance and damping-wing absorption and quantify how they depend on the progress of reionization. To this end, we shall analyze the Cosmic Dawn II (CoDa II) simulation of Ocvirk et al. (2020), which covers a suitable dynamic range to compute the IGM opacity during reionization and the relevant physics. This work parallels Gronke et al. (2021) in that both works are based on similar IGM transmission calculations from the same dataset. Whereas Gronke et al. (2021) mainly focused on the transmission of the resonance to understand the recently observed blue peaks from some LAEs at z 6.5 (Songaila et al. 2018;Matthee et al. 2018;Meyer et al. 2021), we focus on the red side that is relevant to the evolution of the LAE number density during reionization.
The rest of the paper is as follows. In Section 2, we describe the CoDa II simulation and how we calculate the IGM transmission. In Section 3, we present our main results. In Section 4, we discuss our results in light of the current understanding of LAEs. For the numerical calculations, we adopt the cosmology parameters used by the CoDa II simulation for self-consistency. The parameters are Ω m,0 = 0.307, Ω b,0 = 0.048, and h = 0.678 (Planck Collaboration et al. 2014).

Cosmic Dawn II simulation
We compute the IGM transmission of Lyα photons during the Epoch of Reionization based on the simulated galaxies and IGM of the CoDa II simulation. CoDa II is a fully coupled radiation-hydrodynamic simulation of the reionization era in a cubic volume of [64 h −1 cMpc] 3 on a 4096 3 regular grid, which resolves atomically cooling halos down to 10 8 M with 200 dark matter particles. The simulation was initialized at z = 150 and run down to z = 5.8 using the RAMSES-CUDATON code (Ocvirk et al. 2016). The global ionization fraction of the IGM reaches 99% at z ≈ 6.2 in the simulation. The simulated galaxies reproduced the observed UV luminosity function at z 6. The details of this simulation are elaborated on in Ocvirk et al. (2020). Here, we briefly list the features that are relevant to our work.
In CoDa II, photoionization is computed self-consistently by casting rays from stellar particles, which form from gas Figure 1. xy plane of the neutral hydrogen density map centered on the most massive galaxy in the CoDa II simulation at z = 8 (top), 7 (middle), and 6 (lower panels). We zoom in from the entire simulation box of (64 h −1 cMpc) 2 to the (8 h −1 cMpc) 2 and (1 h −1 cMpc) 2 subboxes from left to right. The map is color-coded according to the neutral hydrogen density, but the blue/red color coincides with neutral/ionized regions in most cases. In the middle and right column panels, the partially transparent arrows describe the peculiar motion of gas with respect to the galaxy. The white circles in the left panels describe the maximum integration distance in Equation (1), smax = 24 h −1 Mpc, for the local optical depth τL described in Section 2.2. The yellow circles in the right panels show the virial radius of the galaxies at the center, which is the minimum integration distance, smin = Rvir. cells that satisfy certain star-formation conditions. The CoDa II simulation captures a statistically meaningful number of HII regions with its large volume and resolves small-scale structures around galaxies needed to compute the transmission accurately. Figure 1 illustrates this with the HI density in the IGM around the UV-brightest galaxies at z = 6, 7, and 8 from the full-box scale in the left column down to the virial radius scale in the right column. We also note that CoDa II uses the exact speed of light in its radiative transfer calculation and, thus, is free from possible problems arising from us-ing the reduced speed-of-light approximation Ocvirk et al. 2019).
The UV luminosity of each galaxy is computed from the age of the stellar particles based on stellar synthesis models. We name galaxies according to their ranking in UV luminosity at 1600Å(i.e., M UV ≡ M AB1600 ) within the snapshot to which it belongs. For example, we refer to the 100th UV-brightest galaxies as galaxy #0100. The brightest sample galaxy (galaxy #0001) of the z = 7 snapshot is M UV = −23.1 in the UV magnitude, for example, and there  (vα). The blue, black, and red lines describe the cross-section for three different gas temperatures, T = 10 2 , 10 4 , and 10 6 K, respectively. are hundreds of galaxies up to M UV ∼ −19, which makes CoDa II suitable for understanding LAEs observed at z 7. Note that galaxy #n at z = 7 and galaxy #n at z = 8 are generally not the same because the ranking in UV magnitude changes over time.

IGM Transmission Calculation
We calculate the IGM transmissivity for photons emitted around the Lyα frequency from galaxies in the z = 6, 7, and 8 snapshots. In the simulation, the IGM is almost fully ionized at z = 6, while it is 50% and 13.2% ionized at z = 7 and 8, respectively. Throughout this work, we assume neutral hydrogen is the only opacity source given that the dust and metals from star-formation are unlikely to exist in meaningful amount in the IGM during reionization.
The transmissivity, T ≡ e −τ , for a frequency ν e at the rest frame of a source is given by the optical depth τ (ν e ) = smax smin n HI (r)σ α (ν a (r), T (r))a(z)ds, where s is the comoving distance between the source and the IGM, a(z) is the scale factor at redshift z, r = r s + sγ is the comoving coordinate of the IGM along the line-of-sight directionγ from the source location r s , n HI is the neutral hydrogen number density, σ α (ν, T ) is the ensemble-averaged Lyα cross section of gas with temperature T , and ν a (r) = ν e 1 − v pe (r) + sH(z)a(z) c describes the frequency shift in the rest-frame of the IGM due to the IGM peculiar motion v pe and the cosmic expansion rate H(z). Here, v pe is calculated with respect to the peculiar motion of the galaxy, which we obtain by taking the densityweighted average of the velocity field with the virial radius of the host halo. We define the virial radius as the distance within which the mean density becomes 200 times the cosmic mean (R vir ≡ r 200 ). The Lyα cross-section is given by where T 4 ≡ T /[10 4 K] and φ(x) is the Voigt function as a function of the dimensionless frequency Here, ν α ≈ 2.46 × 10 15 Hz is the resonance frequency of Lyα, and ∆ν D = ν α 2k B T m p c 2 = 4.28 × 10 −5 T 0.5 4 ν α is the thermal broadening frequency (Doppler width), k B is the Boltzmann constant, and m p is the proton mass. The Voigt function is given by where a V ≈ 4.7 × 10 −4 T −0.5 4 is the Voight parameter. For fast computation, it is convenient to use the analytic fitting formula provided by Tasitsiomi (2006).
The shape of σ α (ν, T ) is shown in Figure 2 for T = 10 2 , 10 4 , and 10 6 K. The profile of σ α (ν) has a core at the center with the damping wing extending outward. The FWHM is roughly 2, 20, and 200 km s −1 for T = 10 2 , 10 4 , and 10 6 K, respectively. The gas temperature is approximately T = 10 4 K for the ionized IGM and T 10 2 K for the neutral IGM. Some gas is shock-heated to ∼ 10 6 K by supernova feedback, but such gas rarely exists outside the virial radius, where we shall calculate the IGM transmissivity.
For each galaxy, we calculate τ (ν e ) for 2000 random lines of sight to account for the sight-line variation. For each sight line, we evaluate Equation (1) in two segments to obtain τ = τ L +τ DW . Here, τ L is the local contribution within a distance of 24 h −1 cMpc, which is far enough from the source galaxy to include the surrounding HII bubble. τ DW is the dampingwing absorption by distant IGM beyond that distance.
For the local contribution τ L , we capture the local IGM dynamics and the inhomogeneity of reionization by sampling n HI (r), T (r), and v pe (r) along each sight line by interpolating the mesh data of the simulation. In this work, we regard the IGM as the gas outside the virial radius (R vir ) and integrate to 24 h −1 cMpc from the galaxy by setting s min = R vir and s max = 24 h −1 cMpc in Equation (1). We describe this integration range by white circles in Figure 1. Within this range, we do not consider the time evolution of IGM and calculate the opacity from the still snapshot of the target redshift. The circles in the left panels of Figure 1 describe the integration range for τ L , which is large enough to capture the surroundings of the ionized bubble at z = 7 and 8 created by the central galaxy.
For the large-scale opacity contribution from the neutral gas in the IGM, τ DW , we integrate in Equation (1) from s min = 24 h −1 cMpc to z = 6, where the reionization has ended in the simulation. For this large-scale contribution, we only consider the damping-wing opacity from neutral IGM and ignore the peculiar motion. For the neutral hydrogen density, we use the globally averaged value given bȳ wherex HI is the global average neutral fraction and X = 0.76 is the mass fraction of hydrogen in baryons. In the damping wing, the Voigt profile can be approximated as where accounts for the initial frequency offset of emitted photons from the Lyα in the velocity unit: v α ≡ −c∆ν/ν α . Here, we have assumed H ≈ H 0 Ω m (1 + z) 3 given that z 1 during reionization. We calculatex HI by interpolating in redshift the globally averaged values within each snapshot. In the case of τ DW , the integration starts from s min = 24 h −1 cMpc and ends at s max = 263 (482) h −1 cMpc for galaxies at z = 7 (8). We find that roughly 10% (20%) of the photons are scattered due to this large-scale damping-wing opacity at z = 7 (8). Figure 3 shows the IGM transmissivity T = e −τ as a function of frequency in the velocity unit (v α ) for galaxies #0001, #0014, and #0250 in the z = 6, 7, and 8 snapshots. Note that the positive side of the x-axis (v α > 0) is the red side of Lyα. The transmissions for different sight lines are shown as multiple partially transparent curves to depict the sight-line variation. The blue solid line shows the median value at each frequency, and the cyan lines bracket the 68% range around the median.

Resonance Absorption
The cross section of Lyα is large enough to make even the mostly ionized IGM opaque at the resonance. Photons that escape the circumgalactic medium (CGM) at a frequency bluer than the resonance are redshifted to the line center in the rest frame of the IGM at a certain point and encounter this high opacity of the resonance absorption. Notably, all the transmission curves in Figure 3 show a rapid drop at v α = 100 − 300 km s −1 on the red side of the Lyα, implying that the photons are blueshifted to the resonance in the rest frame of IGM.
This red-side resonance absorption occurs for two reasons. First, the peculiar velocity field of the IGM is dominated by the gravitational infall motion toward the source galaxy. Second, the rising IGM density toward the galaxy keeps the IGM optically thick at the resonance frequency despite the rising ionizing intensity from the source galaxy. Without the infall motion and the overdensity of the IGM, the transmission cutoff should appear at or on the blue side of Lyα (e.g., Mason & Gronke 2020).
The velocity map around galaxies in Figure 1 clearly shows the infall motion from all directions. The infall motion extends to at least several comoving Mpc from the galaxy (see middle column panels) and can extend up to tens of Mpc (see Figure 6 of Iliev et al. (2008)). There also are outflow motions generated by star-formation feedback in hot low-density blobs shown as dark red spots. However, the outflowing gas is mostly confined within the virial radius, which we take as the CGM in this work.
In order to describe how the infall motion affects the IGM transmission, we plot gas density, HI density, peculiar velocity, and optical depth toward the −y direction from galaxy #0001 at z = 7 in Figure 4. The plotted quantities are baryon density, ionizing photon density, HI number density, peculiar velocity, the optical depth of IGM to photons emitted at three different frequencies at the galaxy center, and the local density map around the line of sight.
The peculiar velocity shown as the black solid line in panel d is negative (i.e., infalling) up to ∼ 5 h −1 cMpc from the source. The infall velocity becomes as large as ≈ 300 km s −1 when it peaks near s = R vir and gradually decreases toward larger distances. The infall velocity at a distance s is similar to the circular velocity at that distance: where M h is the total galaxy mass (see the green dotted line). This indicates that the gravitational force from the galaxy is the dominant cause of the infall motion. The peculiar velocity profile is also similar to the analytical models of Santos (2004) and Dijkstra et al. (2007). The blueshift in the rest frame of the IGM due to the infall motion results in a high resonance opacity for some photons on the red side of Lyα. The dashed lines in panel d show that the frequency change of the photons emitted at v α,i = −200, 0, and 200 km s −1 due to cosmic expansion along the line of sight. When the dashed line intersects with the solid line (i.e., v α + v pe = 0), the photons face a sudden increase in the optical depth due to resonance opacity (see panel e). The photons that started at v α,i = 0 and 200 km s −1 would never have encountered the resonance cross-section without the infall motion.
Another reason for the red-side absorption is that the infalling gas near the virial radius remains opaque to Lyα photons despite the intense ionizing radiation from the source. From s = 10R vir ≈ 1.8 h −1 cMpc to s = R vir ≈ 0.18 h −1 cMpc, the gas density in in panel a rises by 10 times the cosmic mean while the ionizing radiation intensity (Γ) rises roughly by a factor of a hundred (panel b). The gas temperature remains close to 20,000 K throughout the interval. The 10 times increase in density cancels the hundred times increase in the intensity because the neutral hydrogen density of highly ionized gas (n HI = x HI n H ) goes as n H (n HII /Γ) ≈ n 2 H /Γ in the ionization equilibrium for constant gas temperature. Therefore, the neutral hydrogen density does not fall toward the galaxy as shown in panel c, keeping the IGM optically thick to Lyα.
The transmission cutoff frequency is set by the maximum infall velocity outside the virial radius R vir . In most cases, the infall velocity peaks around s = R vir (panel d), where the infall velocity is close to the circular velocity at the virial radius, V c ≡ GM h /R vir . Here, the galaxy mass, M h , is the total mass that includes both baryonic and dark matter within R vir . For each transmission curve, we measure the maximum infall velocity, v in , by identifying the location of the transmission cutoff, which coincides with the peak of the transmission curve slope. We plot in Figure 5 the probability  The v in distribution has a fairly large scatter but shows a clear galaxy mass dependence. v in is distributed around 300, 200, and 100 km s −1 for galaxies with their masses 1.09, 0.36, and 0.09 × 10 12 M , respectively.
We show the complete statistics of v in from the entire galaxy sample as a function of halo mass in Figure 6 and as a function of UV magnitude in Figure 7. The error bar for each galaxy denotes the 1σ sight-line variation calculated from 2000 sight lines for each galaxy. v in is similar to V c throughout the mass range of our sample, although it is slightly lower at high mass. For example, v in prefers V c at M h = 10 11 M while it prefers 0.85V c at M h = 10 12 M .
We show in Figure 7 that v in also correlates with the UV magnitude. This is due to the strong correlation between M UV and M h , which hardly evolves in time (see Figure 8). We find that the correlation coefficient between M UV and M h is above 95% for the galaxy samples used in this work. The relation between v in and M UV is well described by a fitting function with a 1σ scatter of ∼ 40 km s −1 (see the cyan dashed line in Figure 7). Nonvanishing transmission on the blue side of Lyα often appears at z = 6, while it is mostly zero at z = 7 and 8. At z = 7, only bright galaxies with M UV −23 can transmit a small fraction of blue-side photons. According to Gronke et al. (2021), this is because the typical HI density in ionized IGM is too high ( 10 −9 cm −3 ) to transmit the blue-side photons at z > 6.5. We refer readers to their work for a detailed analysis on this subject.
We note that the residual IGM neutral fraction at the postoverlap phase (z < 6.2) is a factor of ∼ 5 lower than the estimate of Fan et al. (2006). We checked for a subsample of halos that increasing the calculated optical depths by a similar factor did not affect the transmission redward of the cutoff.

Damping-wing absorption
On the red side of the transmission cutoff, photons do not encounter resonance absorption, leading to a much higher transmissivity. Here, the opacity comes from the dampingwing absorption from highly neutral segments of sight lines. The opacity contribution from each segment depends sensitively on the distance to the segment from the source: dτ DW ∝ s −2 , as shown in Equation (5).
The red side of the transmission curves in Figure 3 show the damping-wing opacity for different sight lines, halos, and redshifts. This red-side transmission decreases toward high redshift due to increasing damping-wing opacity. For example, the median transmission at v α = 400 km s −1 for galaxies #0001, #0016, and #0150 is nearly 100% at z = 6, while it drops to 80, 53, and 62% at z = 7 and to 44, 42, and 37% at z = 8, respectively. We find a substantial sight-line-tosight-line variation in τ DW for individual galaxies at z > 6. The 1σ range of the transmission bracketed by the cyan lines is up to 30%, implying that the ionization field surrounding the source galaxies is highly anisotropic.

HII Bubble Size Estimation
If IGM is fully ionized inside the HII bubble and fully neutral outside, the damping-wing opacity, τ DW , is given by integrating Equation (5) from the bubble size s min = R b to the end of reionization s max = s(z = 6). Then, τ DW would depend only on the bubble size R b and the reionization historȳ x HI , and one can constrainx HI and R b by measuring τ DW from observations. Often, an approximate expression for Equation (5) is used for convenience. In that case, one assumes that the opacity comes entirely from the neighborhood of the source galaxy before the redshift term in the integral changes significantly (Dijkstra 2014). The redshift-dependent terms can then be taken out from the integral, giving τ DW, app (v α ) = 5.7 z g + 1 8 x HI (z g ) where z g is the redshift of the source galaxy. In this case, the opacity does not depend on the reionization history along the line of sight, and one can trivially convert the measured τ DW tox HI , given a relation betweenx HI and R b computed from reionization models (e.g., Furlanetto & Oh 2005).   In this work, we define R b to be the value that matches the damping-wing optical depth τ DW with the computed IGM transmissivity at v α = 400 km s −1 (hereafter T 400 ). v in is smaller than 400 km s −1 for all sight lines considered in this work. Thus, the optical depth at 400 km s −1 is purely from the damping-wing opacity.
In Figure 9, we show R b for galaxies #0001 and #0100 at z = 7 as the line contours overplotted on the HI density maps. We computed R b for 360 equally spaced directions (n) on the xy plane and connected the locations of R bn to create a closed contour. The green solid line contour is from the original expression of Equation (5), while the dashed line contour is from the approximate expression of Equation (9).
We note that the original and approximate expressions generally give similar results, although the approximate value occasionally overshoots in some directions. The overshoot tends to be larger when R b is larger (i.e., τ DW is small). This Figure 9. Effective HII bubble size R b for directions on the xy plane is shown on top of the HI density field around galaxies #0001 (upper panel) and #0100 (lower panel) in the z = 7 snapshot. The green solid line contour describes the R b calculated from the original expression (Eq. 5), while the blue dashed line is from the approximate expression (Eq. 9). The left panels show the maps in a 32 h −1 cMpc box, while the right panels show in a 2 h −1 cMpc box. The black dotted lines with angles in degree mark the sight lines described in Figure 12. is because assuming the [1 + z]x HI (z) term in Equation (5) to be constant overestimates τ DW by ∼ 0.1 at z = 7, due to the decline ofx HI toward low redshift. For large HII bubbles like the one surrounding galaxy #0001, the error in τ DW can bias R b more because τ DW is generally smaller. However, galaxy #0001 is a rare extreme case, and most of the galaxies are in smaller HII bubbles like galaxy #0100, where the error is insignificant.
The R b contour generally follows the HII bubble outline with some deviations that make the contour spikier than the actual bubble shape. The deviation between the apparent and effective bubble sizes largely owes to the bubble geometry being more complex than what is assumed in Equation (5): "fully ionized inside the HII bubble and fully neutral outside." When another large HII bubble is located in the vicinity of the HII bubble surrounding the source galaxy, for example, R b would be larger than what Equation (5) gives toward the direction of the neighboring bubble.
Another important factor ignored in Equation (5) is selfshielded systems, which are small and dense pockets of neutral gas in mini halos within HII regions. When a sight line encounters such a clump, it can pick up a substantial amount of opacity even inside a HII region. The −y direction (φ = 270 • ) from galaxy #0001 (upper panels of Fig-Figure 10. Effective HII bubble size R b for sample galaxies in the z = 7 (left) and 8 (right) snapshots. Vertical bars on the data points describe the sight-line variation of individual galaxies. On the right y-axis, we show the corresponding IGM transmission at vα = 400 km s −1 . The black solid line marks the median of R b each at given MUV's. The grey shade brackets the 1σ range of the combined sight-line-to-sight-line and galaxy-to-galaxy variation. Individual data points are shown only for relatively rare bright galaxies to avoid crowding the figure with too many data points. ure 9) is an example that is affected by a clump. This clump is described in detail in Figure 4, where it is shown as a blue spot in panel f at s ≈ 2.3 h −1 cMpc. The baryon density in this neutral clump reaches 30 times the cosmic mean at the peak (panel a). The optical depth for the photon emitted at v α,i = 600 km s −1 rises to 0.3 at this clump (magenta line in panel e). R b for this sight line is 0.33 h −1 cMpc, although the neutral region actually starts at s ≈ 4 h −1 cMpc.
We present the R b statistics of the entire sample with M UV < −18.7 at z = 7 and 8 in Figure 10. We also give the corresponding transmissivity at v α = 400 km s −1 on the right-hand-side y-axis. We also tabulate the median values of T 400 for M UV = −21, − 20, and, −19 in Table 1. The scattered data points with vertical bars in Figure 10 highlight the large galaxy-to-galaxy and sight-line-to-sight-line variation in the damping-wing opacity; the 1σ range depicted by the gray shape shows is larger than the median value. R b is expected to evolve with redshift, which is considered to be the main cause of reduced LAE visibility at high redshifts. The data points of Figure 10 shows this expected redshift dependence of R b and T 400 clearly: the median values of R b and T 400 are around 5 h −1 Mpc and 73% at z = 7, while 2.5 h −1 Mpc and 37% at z = 8, respectively.
The M UV dependence is not easily visible from the data points, but the median curve reveals a weak dependence. At z = 7, the median R b is 5.5 h −1 Mpc for galaxies M UV = −21, whereas it is about 1 Mpc smaller for the galaxies with M UV = −19. Similarly, at z = 8, R b is about 3 h −1 Mpc for M UV = −21 and 2 h −1 Mpc for M UV = −19. The corresponding T 400 drops from 73% to 70% at z = 7 in the same M UV interval, while from 37% to 30% at z = 8.

Self-shielded Systems
Self-shielded systems can add a large opacity to certain sight lines to create a large discrepancy between the effective and actual bubble size as mentioned in Section 3.2. In this section, we estimate the quantitative impact of these systems on the transmissivity by removing these systems in the transmissivity calculation. For this, we first identify highly neutral segments of slight lines with n HI > 10 −5 cm −3 . We then assume all those segments shorter than 0.3 h −1 cMpc are the compact neutral clumps within the HII region and lower the HI density of those segments to 10 −5 cm −3 when evaluating Equation (1). This way, we exclude the damping-wing opacity of the self-shielded system like the one in Figure 4.
In Figure 11, we compare T 400 before and after removing the self-shielded systems for 2000 sight lines of galaxies #0001, #0021, and #0100. The data points that are above the x = y line are the sight lines that are affected by the neutral clumps. The figure shows that some sight lines have substantially higher transmission when clumps are removed. We find that removing the self-shielded clumps has a small impact on the global statistics of the transmission. Removing the clumps increases the median value of T 400 for galaxies #0001, #0021, and #0100 increase from 80.1% to 80.9%, from 82.9% to 83.5%, and from 62.6% to 63.8%, respectively. Also, only 0.78%, 0.69%, and 2.3% of the sight lines experience more than a 10% increase in T 400 , respectively.
However, if we limit the statistics to low-transmission sight lines with T 400 < 0.3, the impact appears much stronger, finding that T 400 increases by more than 10% for 33%, 100%, and 12% of the sight lines, respectively. This result suggests that unusual nondetections of Lyα in UV bright galaxies can be partially attributed to the self-shielded systems blocking the sight line. Figure 11. Transmission at vα = 400 km s −1 before and after removing the self-shielded systems in the opacity integral of Equation (1) as explained in Section 3.3. The results for 2000 sight lines of galaxies #0001, #0021, and #0100 are shown in the left, middle, and right panels, respectively.
We note that our method of removing self-shielded systems is not perfectly accurate. For example, if the edges of two adjacent HII bubbles are closer than 0.3 h −1 cMpc, our method can erroneously regard the thin HI region between the edges as an isolated neutral clump. However, such complications on the bubble edge would have small impacts on the opacity compared to the impact of dense neutral clumps within the edge of a HII region, which is generally closer to the source. We leave more rigorous and comprehensive analyses of the self-shielded systems for future studies.

Transmissivity of the Lyα emission line
We adopt a model spectral energy distribution for the Lyα emission line profile from the CGM to estimate the quantitative impact of IGM absorption. The adopted model is This model assumes that the profile is given by a Gaussian which the peak location and the FWHM are given by the circular velocity of the galaxy V c = GM h /R vir as motivated by recent observations of LAEs at z < 6 (Yang et al. 2016;Verhamme et al. 2018). 2 Here, we only model on the ratio between the intrinsic and transmitted flux for an arbitrary normalization for the intrinsic flux. Also, we do not consider any emission from the blue side of Lyα as the transmission is nearly zero there at z = 7 and 8. We note that, in reality, the spectral shape of the Lyα emission entering the IGM is expected to vary from sight line to sight line due to the complex dynamics of the CGM (e.g., Verhamme et al. 2012, Song, H. et al. in preparation). We take the product of the transmission curve and the intrinsic spectral shape to obtain the IGM transmitted flux: Then, the transmitted fraction of the line emission is given by

Notable Cases
We first visually inspect the transmission curves and pick several individuals that give useful insights into IGM transmissivity. In Figure 12, we plot the F in , F out , and T for several notable sight lines on the xy planes of galaxies #0001 and #0100 at z = 7. We show these sight lines as black dotted lines in Figure 9.
In Figures 12a and b, We compare two sight lines in almost opposite directions (φ = 5 • versus 190 • ) to galaxy #0100 at z = 7. Although the former case has a larger distance to the neutral region (R b = 11 versus 3.7 h −1 cMpc), it has a lower transmitted fraction X Lyα (35% versus 46%) because it has a higher infall velocity (183 versus 114 km s −1 ). This comparison highlights that the infall velocity is another crucial factor for the IGM transmission and its sight-line variation.
In Figure 12b and c, we compare two sight lines that are only 3 • apart (φ = 190 • versus 193 • ). Despite the small difference in direction, the latter case has a much stronger damping-wing absorption (R b = 3.7 versus 0.26 h −1 cMpc), which results in a more than four times lower transmission (X Lyα = 46% versus 11%). This dramatic difference is due to a self-shielded system that is located in the direction of φ = 193 • from the galaxy (lowerright panel of Figure 9). Figure 12d shows an interesting case in which the transmitted spectrum is double peaked even though the intrinsic spectrum is single peaked. This "pseudo" blue peak occurs due to a leakage feature in the transmission curve, which ap- Figure 12. Spectral shape of the model's intrinsic Lyα emission from the circumgalactic medium (dotted line) and transmitted emission through the intergalactic medium (solid line) along with the transmission curve in the inset panel for several sight lines on the xy planes of galaxies #0001 and #0100 at z = 7. Panels a, b, and c are for the azimuthal angles of φ = 5 • , 190 • , and 193 • from galaxy #0100, respectively, and Panel d is for 24 • from galaxy #0001, respectively (see Figure 9). The maximum infall velocity vin, HII bubble size R b , and Lyα line transmissivity XLyα for each sight line are given in the corresponding panel.
pears as a sharp peak in the inset panel. This sight line is in the φ = 24 • direction from galaxy #0001 (upper-right panel of Figure 9), which goes through an extremely low HI density blob with n HI 10 −10 cm −3 near the virial radius shown as a dark red spot on the HI density map. We find that the temperature of this blob is ∼ 10 6 K, implying that the blob was created by a supernova explosion in a neighboring galaxy. This double-peak feature appears for other adjacent sight lines that go through the optically thin spot. We find similar features for the sight lines between φ = 15 • and 30 • for this galaxy on the xy plane. From visual inspections of several other galaxies, we find a few percent of sight lines exhibit similar features. We a leave more comprehensive analysis of this pseudo-double-peak feature for future studies.

Global Statistics
We calculate X Lyα for 2000 sight lines for the galaxies with M UV −19. We show the full statistics of the results in Figure 13, where we show X Lyα as a function of M UV with the galaxy-to-galaxy and sight-line-to-sight-line variations. Similarly to in Figure 10, the error bars on the data points show the 68% range for individual galaxies. The solid Figure 13. Lyα transmissivity XLyα for sample galaxies in the z = 6 (left), 7 (middle), and 8 (right panel) snapshots. Each data point with an error describes the sight-line variation of an individual halo. The black solid line marks the median of XLyα at given MUV. The gray shades bracket the 1σ range of the combined sight-line-to-sight-line and galaxy-to-galaxy variation. Individual data points are shown only for relatively rare bright galaxies to avoid crowding the figure with too many data points.
line marks the median value for the corresponding M UV , and the gray shade covers the 68% range with the galaxy-togalaxy and sight-line-to-sight-line variations combined. We list the 1σ range for M UV = -22, -21, -20, and -19 at z = 6, 7, and 8 in Table 2.
First, X Lyα decreases toward high z. For example, the median at M UV = −21 is 59.3, 38.4, and 19.2 % at z = 6, 7, and 8, where the IGM is 99.9%, 50%, and 13.2% ionized, respectively. This is in accordance with the common expectation that X Lyα is a useful probe of the IGM neutral fraction.
It is also notable that X Lyα is well below 100% at z = 6, where the IGM is nearly 100% ionized. The IGM absorption at z = 6 should be entirely due to the resonance scattering by infalling gas given that the IGM is nearly fully ionized. Thus, we take the z = 6 results as the fully ionized limit and calculate the relative change in X Lyα at higher redshifts. That is, we take the ratio of the median at z = 7 or 8, X Lyα z=7 or 8 , to that at z = 6, X Lyα z=6 , in order to separate out the damping-wing opacity from the neutral IGM. For example, the median transmission fraction of galaxies with M UV = −21 is suppressed to 65% and 32% at z = 7 and 8, respectively, compared to that at z = 6. We list this ratio for M UV = −22, −21, −20, and −19 in Table 3. Both Figure 13 and Table 3 show that UV-fainter galaxies show a stronger suppression in X Lyα than the brighter ones do. For example, the median transmission fraction of M UV = −21 galaxies drops to 35% from z = 6 to 8, while that of M UV = −19 galaxies drops to 15%.
The M UV dependence of X Lyα is often attributed to the difference in HII bubble size between bright and faint galaxies. Interestingly, we find that the bubble size only partially explains the dependence. According to Table 1, the IGM transmissivity at a fixed frequency of v α = 400 km s −1 changed from 73.1% to 70.3% between M UV = −21 and −19 at z = 7. This difference is smaller than the difference in X Lyα between the same M UV 's (65% versus 55%). This extra difference in X Lyα is attributed to the difference in the velocity offset of transmitted flux (Figure 7), caused by the M UV dependence of the IGM infall velocity v in and the peak location of our intrinsic line profile. That is, brighter galaxies tend to transmit redder photons due to the stronger infall motion, and these redder photons are subject to smaller damping-wing opacity for a given HI density.
We illustrate our understanding of the M UV dependence of X Lyα in Figure 14, where we show the median (transmitted) emission line profile for bright (M UV < −22) and faint (M UV > −19) galaxy groups at z = 6, 7, and 8. We first construct the median transmission curve for each group using the median of v in and R b for each group. Then, multiply the median transmission by the intrinsic emission profile obtained from Equation (10) using the median circular velocity. The line flux relative to z = 6 drops to 70% (37%) for the bright group while to 58% (12%) for the faint groups at z = 7 (8). The figure shows both bubble size and the wavelength of the transmitted photons create an M UV dependence in the damping-wing opacity as we explained above.
We also highlight the large scatter in the X Lyα − M UV relation. Due to the complicated nature of Lyα radiative processes, galaxies at the same redshift with the same brightness can transmit a significantly different number of Lyα photons from sight line to sight line and from galaxy to galaxy. As X Lyα decreases toward low brightness and high redshift, the variation in X Lyα becomes larger relative to the median. The 1σ range of X Lyα for M UV = −21 at z = 7 is 59.3 +19.8 −22.8 %, while that for M UV = −20 at z = 8 is X Lyα = 7.3 +9.2 −5.2 % (see Tab. 2). The ratio of the 1σ range to the median is 43/59 ≈ 0.7 in the former case, while 14.4/7.3 ≈ 2 in the latter case. This means that one needs to measure Lyα line  −19) galaxy groups, respectively. The left, middle, and right panels show the results for z = 6, 7, and 8 galaxies, respectively. The median transmission curve is constructed using the median infall velocity vin and effective HII bubble size R b for each group, the values of which are shown on the plot. The intrinsic profile is given by Eq. 10, where the circular velocity Vc is given by the median M h of each group. The integrated line flux at z = 7 or 8 with respect to z = 6 is also shown on the plot for each group in percent. strengths of at least 50 (400) galaxies to constrain X Lyα to 10% in the former (latter) case assuming perfectly accurate measurements.

Summary and Discussion
We have computed the IGM transmission for the Lyα emission line from star-forming galaxies during reionization from the Cosmic Dawn II simulation to understand what aspects of the IGM determine the transmission. Thanks to the star-formation physics and the subsequent radiative transfer implemented in a large box with high resolution, we were able to compute the IGM transmissivity for galaxies with M AB1600 −22 while accounting for the small-scale physics around the virial radii of the galaxies. We find that the IGM infall motion, HII bubble size, and self-shielded neutral gas systems are the critical factors. We discuss the implications of our results for the current understandings of LAEs below.

Significance of Infall Motion for the IGM Transmissivity
Previous works have repeatedly reported the impact of the infalling IGM on the Lyα transmission (Santos 2004;Dijkstra et al. 2007;Iliev et al. 2008;Dayal et al. 2011;Laursen et al. 2011;Bolton & Haehnelt 2013;Sadoun et al. 2017;Weinberger et al. 2019;Garel et al. 2021). Although galactic outflow from supernova explosions may extend beyond the virial radius for some occasions, the outflow is confined within the virial radius in most sight lines (e.g., Shen et al. 2013). Thus, the peculiar motion of IGM outside the virial radius is dominated by the gravitational field of the source galaxy. The gravitationally infalling IGM suppresses a sig-nificant portion of the red-peak emission even after reionization has ended, requiring higher intrinsic Lyα luminosity to explain the observed Lyα luminosity function (Dayal et al. 2011). The absorption by infalling IGM extends to v α ∼ 100 km s −1 in Laursen et al. (2011) and Garel et al. (2021), for example, while it goes to v α 200 km s −1 in our work and Weinberger et al. (2019). This difference likely resulted from larger simulation box sizes of the latter two studies, which involve more massive galaxies with stronger infall motion.
The truncation frequency of the IGM transmission is given by the maximum infall velocity along the line of sight. In most sight lines, the infall velocity peaks around the virial radius of the galaxy. If one uses a distance larger than the virial radius to define the IGM, the maximum infall velocity would decrease. In this case, some of the absorption by the infalling IGM in this work would be regarded as the absorption by the CGM. This is illustrated in Figure 6 of Weinberger et al. (2019), where they define gas at R vir < s < 5R vir , 5R vir < s < 10R vir , and s > 10R vir as the "inner CGM", "outer CGM", and "IGM", respectively. All three components are regarded as the IGM in this work.
We note that AGNs, which the CoDa II simulation does not include, can result in a much higher transmission owing to much stronger outflow and ionizing radiation. Jets from AGNs can disrupt the infall motion well beyond the virial radius, and the HI density is low enough to be transparent at the resonance (< 10 −10 cm −3 ) within the proximity zones. Indeed, low-z observations show that luminous LAEs are often associated with AGNs (Konno et al. 2016). Neighboring galaxies within the proximity zone can also receive a boost in the Lyα transmission as was reported by Bosman et al. (2020). Also, the X-ray from AGNs can change the reionization morphology by partially ionizing the surrounding HI region (Kakiichi et al. 2017).
The relation between v in and the galaxy mass (or UV magnitude) does not evolve much with redshift as the infall motion is purely gravitational in nature. However, v in can still affect the IGM transmission by setting the minimum wavelength of transmitted photons. Stronger infall motion results in the transmission of redder photons, which experience less damp-wing opacity. As a result, Lyα emission from brighter galaxies, which are likely more massive, are less subject to increasing neutral IGM fraction toward higher redshifts.
The infall motion varies depending on the viewing angle due to the complexity of density structures around the source galaxy. This variation adds to the statistical uncertainty of the transmitted flux from a single observation. We find the infall velocity typically has a scatter of ∼ 40 km s −1 from sight line to sight line (Figure 7). The sight-line difference in the infall velocity can overwhelm that in the damping-wing opacity in determining the total IGM transmissivity (e.g., Figure 12a versus b).

Observational Signature of the Infall Motion
Truncation of the red peak by the infalling IGM can make the spectral shape of the transmitted flux blue-skewed in frequency even if the intrinsic emission is symmetric (Figure 12). Such blue-skewed spectra are often observed from recent high-z spectroscopic LAE surveys (see, e.g., Figure 4 of Jung et al. 2020).
The infall motion is a potential explanation for the offset between the observed line and the systemic redshift of the source galaxy. The absorption by infalling opaque IGM will force the peak to be on the red side of v α = v in ≈ V c regardless of the intrinsic emission. Given the correlation between the UV magnitude and the infall velocity (Figure 7), the offset is expected to correlate with the UV magnitude, and this is supported by existing observations (e.g., see Figure 10 of Hashimoto et al. 2019).
The pseudo-blue-peak feature described in Figure 12d is an interesting possible outcome of infalling IGM. Given that the IGM is unlikely to transmit photons from the blue side of Lyα at z 6.5 (Figure 3; also see Gronke et al. 2021), the currently existing double-peak cases at z ∼ 6.5 without exact systemic redshifts (Matthee et al. 2018;Songaila et al. 2018;Meyer et al. 2021) may be this case. The line shape of our case is the most similar to NEPLA4 of Songaila et al. (2018, see their Figure 12), although the overall shape is about two times narrower in our case. The peak separation of double-peaked LAEs ranges from 100 to 300 km s −1 , which is narrower than the reported value of ∼ 750 km s −1 from z = 2 − 3 LAEs by Kulas et al. (2012). The pseudo-bluepeak scenario can explain the narrow separation because the two peaks come from incomplete resonance absorption of a single red peak. We also note that more recent observation with improved spectral resolution by Matthee et al. (2021) do find lots of narrowly separated double peaks at z = 2 − 3. If any of the reported z > 6 double-peaked LAE falls into the pseudo-double-peak case, its systemic redshift will be found blueward of both peaks in future observation 3 .

Observational Impacts of Self-shielded Systems
Before reionization, mini halos below 10 8 M could accrete gas from the cold IGM and develop dense gas clumps inside their potential wells (Shapiro et al. 2003(Shapiro et al. , 2004. Upon reionization, these clumps were exposed to the ionizing background radiation and started evaporating. However, the evaporation took up to several hundred megayears to finish for relatively massive ones (Iliev et al. 2005;Park et al. 2016;D'Aloisio et al. 2020) and, thus, some of the clumps are expected to have existed within HII regions and attenuated Lyα emission for incoming sight lines (Bolton & Haehnelt 2013;Mesinger et al. 2015;Kakiichi et al. 2016). Previous works relied on analytical or subgrid models of these clumps to demonstrate their impact on large-scale Lyα transmission.
We have revisited this subject with the CoDa II simulation, which directly resolves these dense clumps at (proper) kiloparsec scales while simultaneously reproducing the HII bubbles as large as several (proper) Mpc. We report that the self-shielded systems do not seem to have a large impact on the average transmission, but they seem to be responsible for a large fraction of sight lines with unusually low transmission. Thus, Lyα nondetections from UV-bright galaxies can be at least partially attributed to these dense clumps blocking the sight line. We note that the spatial resolution of CoDa II may have underestimated the opacity from these systems in this work given that high-resolution simulation studies (e.g., Nasir et al. 2021) suggest there are a lot of self-shielded absorbers smaller than the cell size of CoDa II (16 h −1 ckpc).
The statistics of these clumps depend on multiple parameters including the ionizing background radiation intensity, local overdensity, timing of reionization, X-ray preheating, and baryon dark matter streaming motion (Mesinger et al. 2015;Park et al. 2016;D'Aloisio et al. 2020;Cain et al. 2020;Park et al. 2021). The first three parameters are especially relevant to relatively massive ( 10 6 M ) systems that can survive the ionizing background radiation for more than 10 7 yr. The IGM opacity would also depend on these parameters if the unresolved self-shielded systems add a significant amount of opacity to the IGM. Thus, it is worth exploring the Lyman-limit opacity of these absorbers in future studies dedicated to these small-scale structures.

M UV Dependence of IGM Transmissivity
A number of past observations have reported that Lyα emission from relatively brighter galaxies tends to be less affected by the neutral IGM during reionization than fainter ones do (e.g., Fontana et al. 2010;Stark et al. 2011;Pentericci et al. 2011;Ono et al. 2012;Mallery et al. 2012;Curtis-Lake et al. 2012;Treu et al. 2013;Tilvi et al. 2014;Endsley et al. 2021, Jung, I. et al. 2021. We confirm this trend from our results in Figure 13 and Table 3). We illustrate our results by calculating the "median emission profile" for for bright (M UV < −22) and faint (M UV > −19) galaxies in Figure 14, which is obtained by applying the median v in and R b to the intrinsic emission profile created for the median M h .
Based on Figure 14, we conclude that the M UV dependence of the transmission is driven by two factors. First, bright galaxies tend to reside in HII regions larger than where the fainter galaxies reside (R b = 4.7 versus 3.5 h −1 cMpc at z = 7) because they can ionize larger volumes and are more likely to form in overdense regions, where numerous neighboring galaxies power the HII bubble together Endsley et al. 2021). Second, bright galaxies tend to transmit in longer wavelengths because they have a larger velocity offset in their intrinsic emission profile (Mason et al. 2018;Garel et al. 2021, also compare V c in Figure 14) and higher IGM infall velocity (v in = 248 versus 138 km s −1 at z = 7) suppressing the transmission up to a longer wavelength. Often the former is considered to be the major reason for the M UV dependence, but we find the latter is comparably responsible for the dependence. We also find that the M UV dependence grows stronger toward high z. We find the median transmission relative to z = 6 is 69% versus 57% at z = 7 while it is 37% versus 11% at z = 8. This highlights the importance of dealing with the dependence carefully in high-z LAE surveys to avoid biasing the IGM neutral fraction estimation.
The infall velocity v in and the circular velocity V c stay nearly constant between z = 6 and 8, and thus do not directly contribute to the IGM transmission evolution. However, the strength of the M UV dependence depends on the wavelength of the transmitted flux set by those two parameters. In the Appendix, we demonstrate the flux-weighted transmission and its M UV dependence can be strengthened or weakened by modest amounts if we shrink or expand the intrinsic profiles, respectively. However, our main conclusion that the IGM transmission for brighter galaxies is less suppressed by the neutral IGM remains solid in this analysis. While an apparent decline in Lyα visibility at z > 6 has been interpreted as a sign of an increasing IGM neutral fraction, the complex nature of Lyα radiative processes in the IGM limits the measurement accuracy ofx HI based on the Lyα emission strength. The IGM infall velocity and size of HII regions vary from sight line to sight line and galaxy to galaxy even for galaxies with the same mass or brightness, and self-shielded systems occasionally block unfortunate sight lines. These uncertainties make it difficult to pin down the average IGM transmissivity at a given redshift with a small number of measurements. Moreover, the intrinsic Lyα emission from the CGM appears to also depend on not only the viewing angle but also various characteristics of galaxies (e.g., Matthee et al. 2021;Tang et al. 2021) as theoretically expected (e.g., Verhamme et al. 2012;Smith et al. 2019, Song, H. et al. in preparation).

Variation in IGM Transmissivity and Its Implication for Observations
We show in Table 4 the ratio between the 68% range and the median value of X Lyα based on Table 2. The ratio is larger for fainter galaxies at higher redshifts. For galaxies with M UV ≥ −21 at z ≥ 7, the ratio is 1, implying that we need at least a hundred measurements of Lyα emission strength to constrain the mean IGM transmissivity to 10% accuracy with those galaxies. In practice, the required sample size is likely much larger because most of the measurements would yield only nondetections, and there will be extra uncertainty in the intrinsic emission from the CGM.
Recent state-of-art LAE surveys are starting to constrain the neutral fraction of the IGM at z ∼ 7.5, but their constraints have not converged yet: Hoag et al. (2019), Jung et al. (2020), and Whitler et al. (2020) givex HI = 88 +5% −10% , 49 +19% −19% , and 55 +11% −13% , respectively. According to our findings, these studies are based on less than a hundred galaxies, which is not enough to suppress the uncertainty from the IGM transmission variation.
For LAEs at z 7, ground-based experiments are expected to provide a large-enough sample soon (e.g., SILVER-RUSH; Ono et al. 2021). For higher redshifts, future space experiments like the Nancy Grace Roman Space Telescope (NGRST), the James Webb Space Telescope (JWST), and Euclid will substantially improve the constraints with highquality observations. Because brighter galaxies have a relatively smaller transmission variation, surveying bright LAEs in a wide field may Figure 15. Same as Figure 14 except that we reduce Vc by 33% in Eq. 10 to shrink the intrinsic emission profile. be more effective than surveying faint galaxies in a narrow field. In this regard, NGRST and Euclid may be more efficient than JWST in constraining the IGM transmissivity at high redshifts, although a definitive conclusion requires accounting for the variation in the intrinsic emission from the CGM. Figure 16. Same as Figure 14 except that we increase Vc by 50% in Eq. 10 to inflate the intrinsic emission profile.