# Scattering of $Ly\alpha$ Photons through the Reionizing Intergalactic Medium: I. Spectral Energy Distribution

Hyunbae Park, Hyo Jeong Kim, Kyungjin Ahn, Hyunmi Song, Intae Jung, Pierre Ocvirk, Paul R. Shapiro, Taha Dawoodbhoy, Jenny G. Sorce, Ilian T. Iliev

► To cite this version:

Hyunbae Park, Hyo Jeong Kim, Kyungjin Ahn, Hyunmi Song, Intae Jung, et al. Scattering of $Ly\alpha$ Photons through the Reionizing Intergalactic Medium: I. Spectral Energy Distribution. The Astrophysical Journal, 2022, 931, 10.3847/1538-4357/ac69e4. insu-03707834

HAL Id: insu-03707834

https://insu.hal.science/insu-03707834

Submitted on 29 Jun 2022

**HAL** is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire **HAL**, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Scattering of Lyα Photons through the Reionizing Intergalactic Medium: I. Spectral Energy Distribution

Hyunbae Park<sup>1,2,3</sup> <sup>●</sup>, Hyo Jeong Kim<sup>4,5</sup>, Kyungjin Ahn<sup>4</sup> <sup>●</sup>, Hyunmi Song<sup>6</sup> <sup>●</sup>, Intae Jung<sup>7,8,9</sup> <sup>●</sup>, Pierre Ocvirk<sup>10</sup>, Paul R. Shapiro<sup>11</sup> <sup>●</sup>, Taha Dawoodbhoy<sup>11</sup>, Jenny G. Sorce<sup>12,13</sup>, and Ilian T. Iliev<sup>14</sup> <sup>●</sup>

## Abstract

During reionization, a fraction of galactic Lyα emission is scattered in the intergalactic medium (IGM) and appears as diffuse light extending megaparsecs from the source. We investigate how to probe the properties of the early galaxies and their surrounding IGM using this scattered light. We create a Monte Carlo algorithm to track individual photons and reproduce several test cases from previous literature. Then, we run our code on the simulated IGM of the CoDaII simulation. We find that the scattered light can leave an observable imprint on the emergent spectrum if collected over several square arcminutes. Scattering can redden the emission by increasing the path lengths of photons, but it can also make the photons bluer by upscattering them according to the peculiar motion of the scatterer. The photons emitted on the far blue side of the resonance appear more extended in both frequency and space compared to those emitted near the resonance. This provides a discriminating feature for the blueward emission, which cannot be constrained from the unscattered light coming directly from the source. The ionization state of the IGM also affects the scattered light spectrum. When the source is in a small H II region, the emission goes through more scatterings in the surrounding H I region regardless of the initial frequency and ends up more redshifted and spatially extended. This can result in a weakening of the scattered light toward high z during reionization. Our results provide a framework for interpreting the scattered light to be measured by high-z integral-field-unit surveys.

*Unified Astronomy Thesaurus concepts: Reionisation (1383); Early universe (435); Intergalactic medium (813); Observational cosmology (1146)*

## 1. Introduction

A substantial fraction of high-z galaxies during reionization emit strongly in Lyα due to the recombination in the photoionized or collisionally ionized gas in star-forming regions (Dijkstra 2014). Along with the Lyman break feature used to detect Lyman break galaxies (LBGs), the Lyα emission is the main tool for identifying high-redshift galaxies.

UV photons with energies above the Lyman limit (13.6 eV) are a major source of Lyα radiation (Partridge & Peebles 1967). While star-forming clumps in the interstellar medium are considered the main source of Lyα, the ionized intergalactic medium (IGM) can also work as a diffuse Lyα source at large scales (e.g., Fernandez & Komatsu 2006). Strong ionizing radiation from active galactic nuclei or highly star-forming galaxies can also “illuminate” nearby cold non-star-forming

clouds and turn them into Lyα emitters (LAEs; Cantalupo et al. 2012; Rosdahl & Blaizot 2012).

High-z LAEs are being considered as a promising probe of reionization for next-generation surveys. During reionization, the IGM would remain neutral in regions that are far from bright ionizing sources. In these H I regions, the Lyα emission would be suppressed due to the extended damping-wing cross section of a hydrogen atom. Indeed, a steep decline in the LAE number density is found above z=6, in contrast to its more gradual evolution at lower redshifts. This decline appears steeper for fainter galaxies ($M_{UV} \gtrsim -20$) (Malhotra & Rhoads 2004; Fontana et al. 2010; Ouchi et al. 2010; Pentericci et al. 2011; Stark et al. 2011; Curtis-Lake et al. 2012; Ono et al. 2012; Malleroy et al. 2012; Treu et al. 2013; Tilvi et al. 2014; Zheng et al. 2017; Endsley et al. 2021), while the decline appears less dramatic for the brighter ones (Matthee et al. 2015; Castellano et al. 2018; Jung et al. 2019; Harikane et al. 2019; Higuchi et al. 2019; Jung et al. 2020, 2021; Tilvi et al. 2020; Hu et al. 2021), consistent with theoretical expectations for the “inside-out” reionization scenario (e.g., Malhotra & Rhoads 2004; Sadoun et al. 2017; Mason et al. 2018; Katz

<figure><img src="image_3.png" /></figure>

Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

<footnote n="1"><sup>1</sup>Lawrence Berkeley National Laboratory, CA 94720-8139, USA; hyunbae.park@lbl.gov</footnote>

<footnote n="2"><sup>2</sup> Berkeley Center for Cosmological Physics, UC Berkeley, CA 94720, USA</footnote>

<footnote n="3"><sup>3</sup>Kavli IPMU (WPI), UTIAS, The University of Tokyo, Kashiwa, Chiba 277-8583, Japan</footnote>

<footnote n="4"><sup>4</sup>Department of Earth Sciences, Chosun University, Gwangju 61452, Republic of Korea; hyojeong@gist.ac.kr, kjahn@chosun.ac.kr</footnote>

<footnote n="5"><sup>5</sup>Gwangju Institute of Science and Technology, Gwangju 61005, Republic of Korea</footnote>

<footnote n="6"><sup>6</sup>Department of Astronomy and Space Science, Chungnam National University, Daejeon 34134, Republic of Korea</footnote>

<footnote n="7"><sup>7</sup>Department of Physics, The Catholic University of America, Washington, DC 20064, USA</footnote>

<footnote n="8"><sup>8</sup>Astrophysics Science Division, Goddard Space Flight Center, Greenbelt, MD 20771, USA</footnote>

<footnote n="9"><sup>9</sup>Center for Research and Exploration in Space Science and Technology, NASA/GSFC, Greenbelt, MD 20771, USA</footnote>

<footnote n="10"><sup>10</sup>Observatoire Astronomique de Strasbourg, 11 rue de l'Université, F-67000 Strasbourg, France</footnote>

<footnote n="11"><sup>11</sup>Department of Astronomy, University of Texas, Austin, TX 78712-1083, USA</footnote>

<footnote n="12"><sup>12</sup>Université Paris-Saclay, CNRS, Institut d'Astrophysique Spatiale, F-91405, Orsay, France</footnote>

<footnote n="13"><sup>13</sup>Leibniz-Institut für Astrophysik, An der Sternwarte 16, D-14482 Potsdam, Germany</footnote>

<footnote n="14"><sup>14</sup>Astronomy Centre, Department of Physics & Astronomy, Pevensey III Building, University of Sussex, Falmer, Brighton BN1 9QH, UK
Received 2022 February 11; revised 2022 April 9; accepted 2022 April 22; published 2022 June 2</footnote>
et al. 2019; Gronke et al. 2021; Park et al. 2021; Qin et al. 2022; Smith et al. 2022).

Since the IGM has a negligible amount of dust in the high-z universe, most of the scattered light would make it to the present-day universe in the form of extended diffuse light around the source galaxies. We expect nearly 100% of the sky is covered by this kind of diffuse light (Wisotzki et al. 2018), as suggested by recent observations (Leclercq et al. 2017). Lyα blobs (LABs; e.g., Francis et al. 1996; Steidel et al. 2000; Matsuda et al. 2004, 2012; Yang et al. 2014; Kim et al. 2020) are considered to be the scattered light characterized by high Lyα luminosity (~$10^{43}-10^{44}$ erg s⁻¹) and their spatial extent of ~30–200 kpc, although the emission may also be powered by collisional excitation and photoionization due to dynamical interactions between galaxies and the nearby IGM (e.g., Haiman et al. 2000; Dijkstra & Loeb 2009; Faucher-Giguère et al. 2010). Many LABs are in proximity to LAEs (Matsuda et al. 2004), but some are associated with LBGs (Steidel et al. 2000) or active galactic nuclei (e.g., Bunker et al. 2003; Colbert et al. 2011).

The scattered light can extend to megaparsecs with low surface brightness (Zheng et al. 2011). Croft et al. (2018) reported a direct detection of the emission up to 15 comoving Mpc (cMpc hereafter) around a QSO at $z \sim 2-3.5$. Simulation studies find that such large-scale Lyα emission is spatially correlated with the location of LAEs and depends on the IGM ionization state and the intrinsic emission spectrum of the source (e.g., Jeeson-Daniel et al. 2012; Visbal & McQuinn 2018). Recent narrowband surveys are starting to statistically detect the scattered light at $z \sim 6$ as a cross-correlation signal between the LAEs and the surface brightness, although the reionization effect has not been confirmed yet (Kakuma et al. 2021; Kikuchihara et al. 2021).

Therefore, understanding the transfer process of Lyα radiation is essential for constraining reionization from scattered light. The randomness of the scattering process requires a Monte Carlo-type calculation to reproduce the observation. Theoretical studies have evolved from assuming a simplified (e.g., static and symmetric) configuration around a source (Harrington 1973; Neufeld 1990; Loeb & Rybicki 1999; Ahn et al. 2000, 2001, 2002; Zheng & Miralda-Escudé 2002; Tasitsiomi 2006) to more realistic configurations (nonstatic, three-dimensional, or both: Ahn et al. 2002; Zheng & Miralda-Escudé 2002; Cantalupo et al. 2005; Dijkstra et al. 2006; Tasitsiomi 2006; Verhamme et al. 2006; Semelin et al. 2007; Laursen et al. 2009; Yajima et al. 2012). Dust can substantially reduce the Lyα escape fraction $f_α$ (Hansen & Oh 2006; Verhamme et al. 2006, 2008; Schaerer & Verhamme 2008; Yajima et al. 2012; Hutter et al. 2014), but this effect is unimportant in the IGM, which is presumably dust-free. Due to the extreme dynamic range of the interstellar density field, it is still challenging to obtain converged results from modern high-resolution hydrodynamic simulations (Camps et al. 2021). However, the IGM-scale calculation is relatively free from this issue, owing to the much milder dynamic range of the IGM density.

Our primary goal is to understand how the properties of the scattered Lyα light depend on the ionization states of the IGM and the properties of the source galaxies so that we can constrain them from future observations. To this end, we (1) develop our own 3D Monte Carlo Lyα transfer code, and (2) simulate the Lyα photons’ scattering process in the IGM during

reionization. Our Lyα transfer code includes the essence of the already existing calculation schemes by, e.g., Zheng & Miralda-Escudé (2002), Cantalupo et al. (2005), Tasitsiomi (2006), Dijkstra et al. (2006), Verhamme et al. (2006), Semelin et al. (2007), Laursen et al. (2009), and Yajima et al. (2012). On top of this, we shall make extra efforts to implement the cosmological redshift effect and to interpolate the discrete mesh quantities. Our work revisits some of the work by Zheng et al. (2010) with a higher-resolution data set and an inhomogeneous ionization field of the IGM.

The base field for a realistic Lyα radiative transfer (RT) calculation is given by the Cosmic Dawn II simulation (CoDaII; Ocvirk et al. 2020). CoDaII solves fully coupled radiation transfer, hydrodynamics, and gravity to reproduce the density/velocity/ionization/temperature fields during the reionization era in a cosmological volume of [64 h⁻¹ Mpc]³ on a 4096³ mesh. CoDaII reproduces the observed statistical properties of galaxies at $z \gtrsim 6$ (Ocvirk et al. 2016) and fits into the current constraints on reionization (Ocvirk et al. 2020). Its mesh-type output data make it suitable for Monte Carlo Lyα RT calculation.

This paper is organized as follows. In Section 2, we describe the Monte Carlo Lyα transfer code and the relevant equations. In Section 3, we reproduce the known analytic solutions of test problems and validate the accuracy of our code. In Section 4, we present the results of applying our code to the CoDaII data. We summarize and discuss our results in Section 5.

## 2. Method

### 2.1. Basic Equations

The distance that a photon propagates until being scattered is determined by the optical depth of its path. The optical depth to the scattering ($\tau_s$) is drawn from the exponential probability distribution of $P(\tau_s) = e^{-\tau_s}$. For a photon emitted at a frequency $\nu$ in a direction $\hat{n}_i$ from a location $\mathbf{r}$, the optical depth for a propagation distance $s$ can be calculated from the H I number density $n_{\text{HI}}$, the gas temperature $T$, and the bulk gas velocity $V_{\text{pe}}$. Specifically, the thermal velocity of H atoms in the propagation direction, $v_\parallel \equiv v_{\text{th}} \cdot \hat{n}_i$, and the bulk motion in the propagation direction, $V_\parallel \equiv V_{\text{pe}} \cdot \hat{n}_i$, enter the equation:

$$ \tau_{\nu}(s) \equiv \int_{0}^{s} ds' n_{\text{H1}}(\mathbf{r}') \int_{-\infty}^{\infty} dv_{\parallel} P(v_{\parallel}; \mathbf{r}') \sigma(\nu'(\mathbf{r}')). \quad (1) $$

Here, $\mathbf{r}' = \mathbf{r} + s'\hat{n}_i$ is the photon location after the propagation, $P(v_\parallel)$ is the probability distribution of $v_\parallel$, and $\sigma(\nu')$ is the Lyα cross section as a function of the photon frequency in the H atom frame $\nu'$. In the H atom frame, the frequency is shifted from the original value at emission, $\nu$, according to the peculiar motion of the atom and the cosmological redshift during the propagation:

$$ \nu'(\mathbf{r}') = \nu - \nu \frac{v_{\parallel} + V_{\parallel}(\mathbf{r}') + s'H(z)}{c}, \quad (2) $$

where $H(z)$ is the cosmological expansion rate and $c$ is the speed of light. For a gas temperature $T$, the thermal velocity distribution is given by

$$ P(v_{\parallel}; \mathbf{r}') = \frac{1}{\sqrt{\pi} v_{\text{th}}(\mathbf{r}') \exp\left(-\frac{v_{\parallel}^2}{v_{\text{th}}^2(\mathbf{r}')}\right)}, \quad (3) $$
where $v_{\text{th}}(\mathbf{r}') = \sqrt{2k_B T(\mathbf{r}')/m_H}$ is the mean thermal velocity of hydrogen atoms at the location $\mathbf{r}'$, $k_B$ is the Boltzmann constant, and $m_H$ is the mass of a hydrogen atom. The Ly$\alpha$ scattering cross section is given by

$$ \sigma(\nu) = f_{12} \frac{\pi e^2}{m_e c} \frac{\Delta\nu_L/2\pi}{(\nu - \nu_0)^2 + (\Delta\nu_L/2)^2}, \quad (4) $$

where $f_{12} = 0.4167$ is the Ly$\alpha$ oscillator strength, $e$ is the electron charge, $m_e$ is the electron mass, $\Delta\nu_L = 9.936 \times 10^7$ Hz is the natural line width, and $\nu_0 = 2.466 \times 10^{15}$ Hz is the Ly$\alpha$ frequency.

Plugging Equations (3) and (4) into Equation (1) gives

$$ \tau_{\nu} = 34.61 \left( \frac{T}{10^4 \text{ K}} \right)^{-0.5} \times \int_0^s \left( \frac{ds'}{\text{kpc}} \right) \left( \frac{n_{\text{H}}(\mathbf{r}')}{1.899 \times 10^{-7} \text{ cm}^{-3}} \right) \mathcal{H}(a, x). \quad (5) $$

Here, $\mathcal{H}$ is the Voigt function defined as

$$ \mathcal{H}(a, x) = \frac{a}{\pi} \int_{-\infty}^{\infty} \frac{e^{-y^2}}{(x - y)^2 + a^2} dy, \quad (6) $$

where $a \equiv \Delta\nu_L/2\Delta\nu_D = 4.702 \times 10^{-4}(T/10^4 \text{ K})^{-0.5}$ is the ratio of the natural line width to the Doppler line width $\Delta\nu_D = \nu_0(v_{\text{th}}/c)$ and

$$ x \equiv \frac{\nu - \nu_0(1 + V_{\parallel}/c + Hs'/c)}{\Delta\nu_D} \quad (7) $$

is the dimensionless frequency in the gas frame.

In practice, the integral form of the Voigt function (Equation (6)) becomes a nuisance in the numerical calculation, and therefore we instead use a fitting formula given by Equations (7) and (8) of Tasitsiomi (2006), which gives an error of less than 1% for $T \gtrsim 2$ K.

When the optical depth of the photon path reaches $\tau_s$, the photon is scattered by an H atom in a new direction $\hat{n}_f$. In this work, the new direction $\hat{n}_f$ is randomly drawn assuming the scattered photon is isotropically distributed.<fnref n="15" />

During the scattering event, the scattering atom experiences a small recoil $\delta v \sim h\nu_0/(m_H c)$ of the order of a few m s$^{-1}$ depending on the difference between the incoming and outgoing directions of the photon. In the rest frame of the scattering atom, the energy transfer from this recoil is $\sim m_H(\delta v)^2$, which has a negligible impact on the energy of the scattered photon. In the global frame, however, the energy transfer is $\sim m_H v_{\text{atom}} \delta v$, where the atom velocity $v_{\text{atom}} = |v_{\text{atom}}|$ is of the order of a few km s$^{-1}$ and makes a significant change to the photon energy, which is described by

$$ x_f = x_i - \frac{v_{\text{atom}} \cdot \hat{n}_i}{v_{\text{th}}} + \frac{v_{\text{atom}} \cdot \hat{n}_f}{v_{\text{th}}} + g(\hat{n}_i \cdot \hat{n}_f - 1) \quad (8) $$

(see also Section 7.3 of Dijkstra 2017). Here, $g = h\nu_0/(m_H c v_{\text{th}}) \approx 2.6 \times 10^{-4}(T/10^4 \text{ K})^{-0.5}$ is the recoil factor in the atom frame, and the atom velocity $v_{\text{atom}}$ is the sum

of the bulk velocity of the gas and the thermal motion of the atom: $v_{\text{atom}} = V_{\text{pe}} + v_{\text{th}}$. We find that excluding the recoil term does not affect the results as was argued by Adams (1971), although we include it in our calculation.

In order to evaluate Equation (8), we need the thermal motion of the scattering atom $v_{\text{th}}$. The dimensionless thermal velocity parallel to the incident photon direction, $u_{\parallel} \equiv v_{\text{th}} \cdot \hat{n}_i/v_{\text{th}}$, is drawn from the probability distribution function,

$$ f(u_{\parallel}) = \frac{a}{\pi \mathcal{H}(a, x)} \frac{e^{-u_{\parallel}^2}}{(u_{\parallel} - x)^2 + a^2}, \quad (9) $$

which accounts for the simultaneous weighting by the thermal motion and the scattering cross section. We write another velocity component perpendicular to $\hat{n}_i$ as $u_{\perp} \equiv |v_{\text{th}} \times \hat{n}_i|/v_{\text{th}}$ and draw from a Gaussian probability distribution:

$$ f(u_{\perp}) = \frac{1}{\sqrt{\pi}} e^{-u_{\perp}^2}. \quad (10) $$

We generate a random azimuthal angle $\phi$ from the flat distribution between 0 and $2\pi$ for the perpendicular component to obtain the three-dimensional thermal velocity $v'_{\text{th}} = v_{\text{th}}(u_{\perp} \cos \phi, u_{\perp} \sin \phi, u_{\parallel})$. We then apply to this vector a 3D rotation that moves $\hat{n}_i$ to the z direction to obtain the thermal velocity in the global frame $v_{\text{th}}$ to be used in Equation (8).

We note that discreteness in the physical quantities is unavoidable for calculations with numerical simulation outputs. In particular, the discreteness in the velocity field can easily lead to inaccurate results because of the steep dependence of $\sigma(\nu)$ on $V_{\parallel}$. Thus, we perform a 3D linear interpolation based on the eight nearest cell centers that enclose the location to calculate the physical quantities mentioned above. We find that enforcing continuity in the velocity field in this way dramatically reduces numerical artifacts throughout our calculation.

## 2.2. Monte Carlo Simulation of Ly$\alpha$ Scattering

Given the stochasticity of the Ly$\alpha$ resonant scattering process, we adopt the Monte Carlo method for our calculation. We generate individual photons and track their paths as they propagate in space. The major steps of the Ly$\alpha$ scattering simulation using the Monte Carlo method are described below.

1. Read 3D gridded data of gas density ($\rho$), hydrogen ionization fraction ($x_H$), peculiar velocity ($V_{\text{pe}}$), and temperature ($T$) fields from reionization simulation output.

2. Create a photon with the initial frequency $\nu_i$, position $\mathbf{r}_i$, and direction $\hat{n}_i$.

3. Draw a random optical depth $\tau_s$ from the scattering probability distribution, $P(\tau_s) = e^{-\tau_s}$.

4. Accumulate $\tau$ according to the propagation distance $s$. We propagate 0.1% of the grid size at a time to evaluate Equation (1) while keeping track of the H I density, gas velocity, and temperature on the way.

5. When the scattering happens after the optical depth reaches $\tau_s$, we update the position vector to be $\mathbf{r}_f = \mathbf{r}_i + s\hat{n}_i$.

6. Draw thermal velocity components for the atom that scattered the photon, $u_{\parallel}$ and $u_{\perp}$, from Equations (9) and (10).

<footnote n="15">¹⁵ We note that this is a simplification. In the literature, Rayleigh scattering with a probability distribution $P(\mu) \propto 1 + \mu^2$ with $\mu \equiv \hat{n}_i \cdot \hat{n}_f$ is usually considered. Nevertheless, this simplified, isotropic scheme has passed all our test problems in Section 3, indicating that the effect of anisotropic scattering is smeared out after multiple scattering events. Therefore, we adopt this isotropic-scattering scheme for our calculation.</footnote>
7: Draw a new direction vector $\hat{n}_f$ from an isotropic probability distribution.

8: Calculate the new frequency $\nu_f$ and the new direction vector $\hat{k}_f$ of the scattered photon using $V_{pe}$ and $\nu_{th}$.

9: Replace the initial frequency and the initial direction vector by their final quantities from previous steps. Repeat Steps 3–8 until the photon escapes the simulation box.

10: Sample the location ($r_{es}$), direction ($\hat{k}_{es}$), and frequency ($\nu_{es}$) at the final scattering before escape.

11: Repeat Steps 2–10 until accumulating enough photons to draw statistics (typically $10^6$).

We calculate the Ly$\alpha$ scattering in the source comoving frame and do not use any particular code acceleration method such as the core-skipping acceleration scheme (e.g., Ahn et al. 2000; Semelin et al. 2007). Due to the serial nature of the calculation, the code can easily be parallelized to multiple cores with shared memory. The computational cost depends sensitively on the optical depth of the system as well as other parameters such as simulation box size, mesh grid, etc. In the application to the CoDaII data presented in Section 4, the calculation for $10^6$ photons takes several minutes to an hour depending on how neutral the IGM is around the galaxy.

## 2.3. Cosmological Redshift

In our code validation tests presented in Section 3, the physical sizes of the systems are small enough that cosmological redshift is negligible, and all the scattered photons are within the observational field of view. In this case, we can directly use the final frequency at the escape of the system ($\nu_{es}$) to obtain the observed spectrum of the scattered light.

For the application to the CoDaII simulation data in Section 4, however, the scattered light extends to several comoving megaparsecs and beyond, making cosmological redshifting an important factor in the spectrum of the scattered photons. Also, some of the scattered light may not be sampled depending on the transverse distance between the photon and the source, $r_\perp$, and the survey design.

In order to calculate the observed spectrum, one must align the photons on the same sky plane with the source. In order to do so, we apply a path length correction to the final frequency of the last-scattered photons ($\nu_{es}$) to obtain the frequency on the sky plane of the source galaxy ($\nu_{obs}$):

$$ \nu_{\text{obs}} = \nu_{\text{es}} + \nu_{\text{es}} c^{-1} H r_{\text{es}} \cdot \hat{k}_{\text{es}}, \quad (11) $$

where $r_{es}$ is the location of last scattering and $\hat{k}_{es}$ is the final photon direction. Then, the transverse distance to the source on the sky plane is given by

$$ r_{\perp} = r_{es} \sqrt{1 - (\hat{r}_{es} \cdot \hat{k}_{es})^2}. \quad (12) $$

The additional path length term, $r_{es} \cdot \hat{k}_{es}$, in Equation (11) and the projected distance are illustrated in Figure 1. We calculate $\nu_{obs}$ and $r_\perp$ for each sampled photon, assuming the observer is in the final photon direction. We then combine the statistics of all the photons escaping in different directions, effectively averaging the observations of one halo from all sightlines.

## 3. Code Validation Test

We test the Ly$\alpha$ scattering code for several simplified cases studied by previous works. We create mesh-type initial conditions for those models and run our Monte Carlo Ly$\alpha$

<figure><img src="image_4.png" /><figcaption>Figure 1. Schematic description of the path length correction of Equation (11) for calculating the observed frequency $\nu_{obs}$ from the final frequency at the last scattering event.</figcaption></figure>

scattering code. We also reproduce some of the results in Section 4 with another well-tested Ly$\alpha$ scattering code by Seon et al. (2022) and find a good agreement. In particular, the first scattering location (defined in Section 4) is perfectly reproduced by their code because it does not involve any randomness.

### 3.1. Static Homogeneous Slab

In the static homogeneous slab test, we locate a monochromatic source at the center of the slab. The slab is static and homogeneous with neutral hydrogen column density $N_{H\,I}$. We try $N_{H\,I} = 10^{18}$, $10^{19}$, and $10^{20} \text{cm}^{-2}$, which correspond to $\tau_0 = 6 \times 10^4$, $6 \times 10^5$, and $6 \times 10^6$, respectively. We generate photons at the line center (i.e., $x=0$), and the line center optical depth $\tau_0$ is calculated from the slab center to the edge. The temperature of the slab is set to $10^4 \text{ K}$ everywhere. The photon escaping from the slab is collected to obtain the emergent spectrum. Neufeld (1990) solved the radiative transfer equation for this configuration and obtained an angular mean intensity of

$$ J(\tau_0, x) = \frac{\sqrt{6} x^2}{24 \sqrt{\pi} a \tau_0} \frac{1}{\cosh[\sqrt{\pi^3/54} (x^3 - x_i^3)/a\tau_0]} \quad (13) $$

We plot the emergent spectrum in Figure 2. The solid lines are results from the Ly$\alpha$ code, and the dashed lines are the analytic solution from Equation (13). The simulation reproduces the emergent spectrum of the analytic model very well. In the homogeneous slab case, photons diffuse symmetrically in the frequency domain and show double-peak features. Equation (13) was derived for the optically thick case, so the emergent spectrum fits better as $\tau_0$ increases.

The scattering count also has an analytic solution. Harrington (1973) derived that the mean scattering count is

$$ \langle N_{\text{scat}} \rangle = 1.612 \times \tau_0. \quad (14) $$

In Figure 3, we compare the scattering count from our Ly$\alpha$ scattering code to the analytic solution. Our result converges to the analytic solution from Equation (14) in the high-$\tau_0$ limit.
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Intensity of scattered Lyα photons from the static slab model. The agreement between the simulation (solid lines) and the analytic solution (dashed lines) improves as we increase the optical depth of the slab (blue → green → red) because the analytic solution was derived from the optically thick limit.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 3.** Average number of scatterings until the photons escape the system. The results from our Lyα scattering code are shown for $\tau_0 = 6 \times 10^4, 6 \times 10^5$, and $6 \times 10^6$ as blue dots. The analytic solution from Equation (14) is shown as the green dashed line.</figcaption></figure>

### 3.2. Static Homogeneous Sphere

We test the code for a static homogeneous sphere. We generate a cubic volume of $256^3$ cells of homogeneous density, peculiar velocity, ionization fraction, and temperature fields. The H I column density is varied as $N_{H I} = 2 \times 10^{18}, 2 \times 10^{19}, 2 \times 10^{20}$, and $2 \times 10^{21} \text{cm}^{-2}$, which correspond to line center optical depths of $\tau_0 = 1.2 \times 10^5, 1.2 \times 10^6, 1.2 \times 10^7$, and $1.2 \times 10^8$, respectively; the line center optical depth is defined as the optical depth between the center of the sphere and the surface. The density field is generated to meet the intended column density. The sphere size is set to 10 kpc. The small sphere size allows the code to ignore the effect of cosmic expansion. We place a monochromatic Lyα source at the center of the sphere and set the temperature to $T = 10^4$ K everywhere.

Dijkstra et al. (2006) derived the analytic solution for this configuration. In the solution, angle-averaged mean intensity is

<figure><img src="image_7.png" /><figcaption>**Figure 4.** Emergent SED from the homogeneous and isothermal static sphere for $\tau_0 = 1.2 \times 10^5, 1.2 \times 10^6, 1.2 \times 10^7$, and $1.2 \times 10^8$, shown as the solid lines. The analytic solution from Dijkstra et al. (2006) is shown as the black dashed lines.</figcaption></figure>

given by

$$J(\tau_0, x) = \frac{\sqrt{\pi}}{\sqrt{24 a \tau_0} (1 + \cosh[\sqrt{2\pi^3/27}(|x|^3/a\tau_0)]).} \quad (15)$$

The test results are shown in Figure 4. Except for the case of $\tau_0 = 1.2 \times 10^5$, the spectral energy distribution (SED) from the code accurately reproduces the analytic solution with only small deviations from statistical fluctuations. Because the analytic solution is derived for the optically thick limit, $\tau_0 = 1.2 \times 10^5$ shows a greater deviation than other cases with larger optical depths, as was also reported many times by previous studies (Dijkstra et al. 2006; Verhamme et al. 2006; Semelin et al. 2007; Laursen et al. 2009).

### 3.3. Homogeneous Sphere with Hubble-like Flow

On top of the case of a static homogeneous sphere introduced above, we apply a radially outward motion defined by $v_H = V_{max}r/R_{max}$, where $R_{max} = 10$ kpc is the radius of the sphere. We set the column density from the center to the surface as $N_{H I} = 2 \times 10^{18}$ and $2 \times 10^{20} \text{cm}^{-2}$, corresponding to $\tau_0 = 1.2 \times 10^5$ and $1.2 \times 10^7$, respectively. For each density, we set three different maximum velocities at the surface, $V_{max} = -200, 0$, and $200 \text{ km s}^{-1}$, to explore cases of a collapsing, static, and expanding sphere. The results for this configuration are available from previous works (Loeb & Rybicki 1999; Zheng & Miralda-Escudé 2002; Dijkstra et al. 2006; Tasitsiomi 2006; Verhamme et al. 2006; Semelin et al. 2007; Laursen et al. 2009).

The results are shown in Figure 5. In the expanding sphere, the red side of the spectrum is enhanced while the blue part is completely suppressed. This is because all the photons that are emitted on the blue side are scattered in the outskirts, where the expansion is fast enough to shift the frequency to the resonance. The collapsing sphere works oppositely and enhances the blue part of the spectrum. In Figure 6, we fix the H I column density to $N_{H I} = 2 \times 10^{20} \text{cm}^{-2}$ and vary the expansion velocity ($V_{max} = 0, 20, 200$, and $2000 \text{ km s}^{-1}$). As we increase $V_{max}$ from 0 to $2000 \text{ km s}^{-1}$, the blue part is suppressed, and the red peak is extended to longer wavelengths. The red peak shifts toward longer wavelengths as $V_{max}$
<figure><img src="image_8.png" /><figcaption>**Figure 5.** The emergent spectrum of the Lyα source located at the center of static, expanding, and collapsing homogeneous and isothermal spheres. The gas temperature is set to T = 10⁴ K, and the H I density is set so that the optical depth from the center to the surface is τ₀ = 1.2 × 10⁵ and 1.2 × 10⁷ in the left and right panels, respectively. For each density, we explore the cases of a collapsing (green solid), static (red solid), and expanding (blue solid) sphere defined by the maximum radial speed at the surface of Vmax = -200, 0, and 200 km s⁻¹, respectively. The analytic solution of Dijkstra et al. (2006) for the static case is shown as the black dashed line.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 6.** The emergent spectrum of the Lyα scattering in a static, expanding homogeneous and isothermal sphere. N<sub>H I</sub> = 2 × 10<sup>20</sup> cm<sup>-2</sup>, τ<sub>0</sub> = 1.2 × 10<sup>7</sup>, and T = 10<sup>4</sup> K.</figcaption></figure>

increases up to 200 km s⁻¹, but it shifts back toward the line center above a certain threshold value, as can be seen from the 2000 km s⁻¹ case because the steep velocity gradient allows more photons to escape before redshifting further (Laursen et al. 2009). Similar results can be found in Figure 8 of Laursen et al. (2009) and in the right panels of Figures 2 and 3 in Zheng & Miralda-Escudé (2002).

**4. Application to Reionization Simulation Data**

We apply our Lyα RT code to the CoDaII simulation data set. As introduced above, CoDaII is dedicated to reproducing the IGM during the epoch of reionization by simulating the formation of early galaxies and the ionizing radiation from them. The output data include the ionization fraction (χ), density (ρ), peculiar velocity (v<sub>pc</sub>), and temperature (T) fields of gas on a 4096³ mesh of a 64h⁻¹ cMpc box and basic galaxy properties including absolute UV magnitude (M<sub>UV</sub>) and halo mass (M<sub>h</sub>). Our goal is to find how the observed Lyα line shape

would change if we collected the scattered Lyα emission around a UV-bright galaxy along with the unscattered light coming directly from the galaxy in high-z integral-field-unit (IFU) surveys.

We limit the scope of this work to the SED of the scattered light from two galaxies in the z = 7 snapshot, where the IGM is 50% ionized in the simulation. The two galaxies that we analyze are the first and 504th brightest galaxies in the snapshot, which we name as galaxy #0001 and galaxy #0504. The UV magnitudes of these galaxies are M<sub>UV</sub> = -23.1 and -19, and their total masses are M<sub>h</sub> = 1.1 × 10<sup>12</sup> and 6.7 × 10<sup>10</sup> M<sub>☉</sub>, respectively. The former is surrounded by a relatively large H II region (≥5 cMpc), while the latter is surrounded by a smaller one (~2 cMpc). The latter case is meant to represent the early stage of reionization, while the former represents the late stage or the post-reionization regime. At z = 7 in the simulation, the snapshot has a mixture of both cases, making it suitable for exploring both regimes from a single snapshot. We shall provide a more comprehensive analysis with more galaxies for other physical quantities such as the surface brightness profile in our future work.

For each sample galaxy, we trim out the 16h⁻¹ cMpc box with the galaxy at the center. We calculate the gas-density-weighted mean peculiar velocity within r<sub>200</sub> and subtract it from the entire velocity field to work in the source galaxy frame. Then, we generate photons at r<sub>200</sub> from the source galaxy and initialize the photons to propagate radially outward.<sup>16</sup> Given the limited spatial resolution of the simulation below circumgalactic scales, we do not attempt to simulate RT within r<sub>200</sub> and instead focus on the scattering process in the IGM. When a sample photon reaches 8h⁻¹ cMpc from the source, we assume it has escaped the system, and we sample the final frequency (ν<sub>obs</sub>) and the transverse distance to the source r<sub>⊥</sub>, as described in Section 2.3.

Since the gas temperature is not constant in the simulation, the dimensionless frequency x is not convenient for describing the results. Thus, we instead use the wavelength offset from

<footnote n="16"><sup>16</sup> We shall test this case of radial emission against one of nonradial emission in Section 4.3.2 and in the Appendix.</footnote>
Lyα in the velocity unit defined as

$$v_{\alpha} = -c \frac{v_{\text{obs}} - v_0}{v_0} \qquad (16)$$

In this unit, 1 Å roughly corresponds to 250 km s⁻¹ at z = 7.

### 4.1. Scattering Location

As demonstrated in Park et al. (2021, hereafter P21), the residual H I density in the H II regions during the epoch of reionization is generally high enough to keep the IGM opaque at the Lyα resonance ($n_{H I} \gtrsim 10^{-9} \text{ cm}^{-3}$) even in the near-zone of UV-bright galaxies (see also, e.g., Iliev et al. 2008). Thus, a photon emitted on the blue side of the resonance in the IGM frame will eventually redshift to the resonance after the propagation distance $r_s$ given by

$$H_{\delta} + v_{\alpha,i} + v_{\text{pe},r} = 0, \qquad (17)$$

where H is the cosmic expansion rate, $v_{\alpha,i}$ is the initial wavelength of the photon at emission, and $v_{\text{pe},r} = \hat{r} \cdot v_{\text{pe}}$ is the radial peculiar motion of the IGM. Unless there is a highly neutral region on the way, the photon would propagate freely until reaching $r_s$ from the source and be scattered for the first time after its emission.

Due to the gravitational field of the source galaxy, the IGM generally infalls toward the galaxy (i.e., $v_{\text{pe},r} < 0$). As a result, some photons that are emitted on the red side of the resonance in the source frame can be on the blue side in the IGM frame if $v_{\alpha,i} < -v_{\text{pe},r}$. These photons would also redshift to the resonance after propagating a distance $r_s$ from the source.

According to P21, $v_{\text{pe},r}(r)$ is given approximately by $-(GM_h/r)^{0.5}$ with some variations due to the gravitational field of the neighboring density structures. Thus, $r_s$ for a given $v_{\alpha,i}$ forms a near-spherical “first scattering” surface surrounding the source galaxy.

#### 4.1.1. Large H II Bubble Case

In a large H II region, most photons starting blueward of the resonance in the IGM frame ($v_{\alpha} + v_{\text{pe},r} < 0 \text{ km s}^{-1}$) would propagate uninterrupted until they redshift to the resonance. In this case, we expect Equation (17) to accurately give the first scattering location. We test this hypothesis for galaxy #0001, which is surrounded by a relatively large H II region extending beyond $5h^{-1}$ cMpc from the galaxy in most directions, as shown in Figure 7.

We show the first scattering surface on the xy plane for $v_{\alpha,i} = -800, -400, -200, 0$, and $200 \text{ km s}^{-1}$ as line contours in Figure 7. The contour is highly circular because the IGM infall motion is nearly isotropic, and it is larger for smaller $v_{\alpha,i}$ because bluer photons travel a greater distance $r_s$ to reach the resonance. The initially red photons with $v_{\alpha,i} = 200 \text{ km s}^{-1}$ also form the first scattering contour because the gravitational infall motion around galaxy #0001 exceeds $200 \text{ km s}^{-1}$, making them blue-side photons in the IGM frame.<fnref n="17" />

We run our Lyα RT code for the photons that are initially emitted in the xy plane of the source and show their first scattering locations as triangles in Figure 7. The figure shows the scattering locations for 360 photons with their initial azimuthal angle between 0° and 359° and zero latitude angle from the plane. The first scattering location generally coincides

with the $r_s$-contour except for some downward directions (i.e., near the -y direction) for photons with $v_{\alpha,i} = -800 \text{ km s}^{-1}$, which enter the neutral region and are scattered before reaching the contour. Other than this case, all the scattering locations fall slightly inside the contour with a small offset of $\lesssim 0.1h^{-1}$ cMpc. This offset occurs because photons are scattered when the optical depth exceeds one, while the $r_s$-contour marks the peak of the IGM opacity, which happens slightly later. The offset is generally small compared to the value of $r_s$, indicating that it is a good description for the first scattering location in H II regions. The first scattering location, in principle, has a distribution according to the optical depth distribution of the propagation. However, almost all the photons are scattered within a thin surface near the $r_s$-contour because the optical depth remains small until the photon approaches the contour, where the optical depth rises steeply. This behavior of the optical depth is described in detail in Figure 4 of P21.

The last scattering locations shown in the lower panels of Figure 7 also coincide with the $r_s$-contours, indicating that the photons do not travel far between the first and last scatters. The agreement for the redward emission ($v_{\alpha,i} = 200 \text{ km s}^{-1}$) is not as good as that for the bluer emissions, but it is still much better than in the case of a small H II bubble (Figure 8), as we will describe below. We thus conclude that $r_s$ well describes the last scattering location of escaped photons, as well as the first scattering location, in large H II regions.

#### 4.1.2. Small H II Bubble Case

Galaxy #0504 has an H II region of roughly $1h^{-1}$ cMpc (see Figure 8), which is smaller than that surrounding galaxy #0001. In such a small H II region, most Lyα photons would enter the H I region before redshifting to or away from the resonance. Thus, the damping-wing cross section of the H I gas creates a large opacity for these photons, and scatters most of them before they reach $r_s$.

In Figure 8, the first scattering location agrees well with the $r_s$-contour for $v_{\alpha,i} = 0$ and $-200 \text{ km s}^{-1}$, of which the $r_s$-contour lies within the H II region. However, the photons with $v_{\alpha,i} = -400$ and $-800 \text{ km s}^{-1}$ are scattered well inside the $r_s$-contour because the H I gas scatters the photons before they redshift to the resonance with its damping-wing opacity. The first scattering location of the $v_{\alpha,i} = -400 \text{ km s}^{-1}$ case closely follows the boundary of the H II region, indicating that the optical depth is rising steeply there for these photons. On the other hand, the first scattering location of the $v_{\alpha,i} = -800 \text{ km s}^{-1}$ case is more scattered between the H II region boundary and its $r_s$-contour. This is because these photons enter the H I region when their frequencies are relatively far from the line center, where the Lyα cross section is low, and the IGM opacity rises more gradually, resulting in the scattering probability being distributed more extensively in space.

The photons emitted at $v_{\alpha,i} = 200 \text{ km s}^{-1}$ do not encounter the resonance because the peculiar infall velocity of this galaxy is $\sim 150 \text{ km s}^{-1}$, which is not strong enough to put those photons blueward of resonance in the IGM comoving frame. 43% of these photons with $v_{\alpha,i} = 200 \text{ km s}^{-1}$ escape the RT simulation volume unscattered, while the rest are scattered within the H I region due to the damping-wing opacity.

The first and last scattering locations are much less correlated than in the case of a large H II region because a significant fraction of the photons are additionally scattered in

<footnote n="17">¹⁷ See Section 3.1 of P21 for detailed analysis of the infall motion.</footnote>
<figure><img src="image_10.png" /><figcaption>**Figure 7.** H I density map of the *xy* plane containing galaxy #0001 and the first (last) scattering location of the Lyα photons that are initially emitted on that plane in the upper (lower) panels. We generate a photon for every 1 degree from the initial azimuthal angle of $\phi = 0^{\circ}$ to $359^{\circ}$ with zero latitude angle from the plane. The *x-* and *y*-axis ticks are all in units of $h^{-1}$ Mpc. The left panels shows a slice of the entire $16h^{-1}$ cMpc box used for the Lyα RT calculation, and the right panels shows a zoomed-in central region that is $4h^{-1}$ cMpc on a side. The black dotted lines in the left panels mark the boundary of the RT calculation ($8h^{-1}$ Mpc from the source), where we assume the photon has escaped the system and sample the photon information. The red/blue regions on the map generally correspond to ionized/neutral parts of the IGM. The red, black, cyan, blue, and magenta triangle symbols show the scattering location for the photons with initial wavelengths $v_{\alpha,i} = 200, 0, -200, -400, \text{and } -800 \text{ km s}^{-1}$, respectively. The line contours connect $r_s$ from Equation (17) for each direction from the galaxy, where we expected the photons to be scattered for the first time. The last scattering positions are the projected locations on the *xy* plane.</figcaption></figure>

the surrounding H I region. Therefore, the $r_s$-contour is a poor description for the scattered light in a small H II region like this one.

## 4.2. Scattered Light SED: Monochromatic Sources

Before looking into the scattered light for extended emission profiles, we first explore monochromatic cases, where the photons are initially emitted at a fixed wavelength of $v_{\alpha,i} = -800, -400, -200, 0$, or $200 \text{ km s}^{-1}$. These monochromatic cases are not realistic, but they are useful for understanding the relation between input and output spectra, because the output spectrum from an arbitrary input spectrum can be constructed from superpositions of the monochromatic cases.

We initialize the photons at random locations on the sphere of $r_{200}$ with radially outward propagation direction, assuming the source emissivity is isotropic. We obtain the SED by taking the probability distribution of $v_\alpha$ of all the sampled photons, effectively averaging the observations from all possible

viewing angles. We show these results in Figure 9. We show the SED for all the sampled photons at $r_\perp \le 8h^{-1}$ cMpc (or 4.5 arcmin) and for a fraction of photons that are sampled within $r_\perp = 1h^{-1}$ cMpc (or 0.6 arcmin) from the source galaxy. This way, we account for the impact of having a finite light-collecting area on the observed SED. In practical observations, it would be difficult to collect the scattered light beyond ~1 arcmin because of the radiation from other galaxies in the field and the sky noise (Zheng et al. 2011). We shall refer to the former case as the “total” SED and the latter as the “near-zone” SED. As in the previous sections, galaxies #0001 and #0504 represent large and small H II bubbles around the source galaxies.

### 4.2.1. Large H II Bubble Case

The total SED of galaxy #0001 (thick histogram in the left panel of Figure 9) is narrowly peaked at the resonance for $v_{\alpha,i} = 0$ and $200 \text{ km s}^{-1}$ and becomes wider for smaller $v_{\alpha,i}$ (
<figure><img src="image_11.png" /><figcaption>Figure 8. Same as Figure 7 but for galaxy #0504.</figcaption></figure>

i.e., shorter initial wavelengths). For $v_{\alpha,i} = -400$ and $-800 \text{ km s}^{-1}$, the emergent SED shape is similar to a top-hat extending from $v_{\alpha} = v_{\alpha,i}$ to $-v_{\alpha,i}$.

The near-zone SED is similar to the total SED in shape with slightly lower (~30%) intensity for $v_{\alpha,i} = 0$ and $200 \text{ km s}^{-1}$. For bluer emission ($v_{\alpha,i} = -800$ and $-400 \text{ km s}^{-1}$), nearly all the emission around the line center from the total SED is lost, and only the red and blue tips of the distribution are captured in the near-zone SED.

### 4.2.2. Small H II Bubble Case

In the case of a small H II bubble (right panel of Figure 9), the emergent SED is more extended redward than in the case of a large H II bubble. For $v_{\alpha,i} = -800$ and $-400 \text{ km s}^{-1}$, the blue end of the total SED coincides with $v_{\alpha,i}$, as in the case of a large H II bubble, but the SED is much weaker on the blue side. Instead, the SED is substantially more extended redward going beyond $1000 \text{ km s}^{-1}$.

The near-zone SED is also much weaker and extended redward. For $v_{\alpha} = 0$ and $200 \text{ km s}^{-1}$, the blue end of the near-zone SED lies on the total SED, but the redward emission is mostly not included in the near-zone. The near-zone SED is weaker at all wavelengths, and the double-peaked feature seen in the case of a large H II bubble does not appear in this case.

### 4.2.3. Physical Explanation

Here, we provide a physical explanation for the scattered light SEDs of the monochromatic cases shown above. We first describe the SEDs for the case of a large H II region (left panel of Figure 9). Then, the case of a small H II region (right panel of Figure 9) can be understood by considering additional scatterings in the surrounding H I region.

The schematic in the left panel of Figure 10 describes the scattering location and the paths of photons with $v_{\alpha,i} = -800 \text{ km s}^{-1}$, which are emitted on the far blue side of Ly$\alpha$. These photons travel a relatively large distance until being scattered toward the observer. The peculiar infall velocity of the IGM is roughly $(GM_h/r_s)^{0.5} \sim 70 \text{ km s}^{-1}$ at the scattering surface, which is small relative to the initial offset from Ly$\alpha$ ($800 \text{ km s}^{-1}$). Thus, the IGM infall motion is less important for shaping the SED than the cosmological redshift during propagation. We use a color gradation in the arrows to illustrate how the cosmological redshift effect makes photons emitted at the same wavelength end up at different wavelengths. If a photon is initially emitted toward the observer, its path length would be similar to the direct distance to the observer, and the observed wavelength would be similar to the initial wavelength (i.e., $v_{\alpha} \approx v_{\alpha,i}$). When the photon is emitted in the opposite direction and later scattered toward the observer
<figure><img src="image_12.png" /><figcaption>**Figure 9.** Emergent SED of the scattered photons with $v_{\alpha,i} = -800, -400, -200, 0$, and $200$ km s$^{-1}$ shown as magenta, blue, cyan, black, and red histograms, respectively, in the source frame. The results for galaxies #0001 and #0504 are shown in the left and right panels, respectively. The thick lines are for the flux from all the photons sampled at $r_{\perp} < 8h^{-1}$ cMpc, and the thin lines are for those sampled at $r_{\perp} < 1h^{-1}$ cMpc. For galaxies #0001 and #0504, 9% and 43% of the photons with $v_{\alpha,i} = 200$ km s$^{-1}$ escape the system unscattered and are excluded from the histogram, respectively.</figcaption></figure>

at the scattering surface as in the red path, the path length increases by $2r_s$, redshifting the photon to $v_{\alpha} \approx -v_{\alpha,i} = 800$ km s$^{-1}$. Likewise, the black path shows that the photons initially emitted perpendicular to the observer's direction would travel an extra distance $r_s$ and be observed near the line center.

These example paths show that the emergent wavelength of the scattered photon is given by $v_{\alpha} = v_{\alpha,i} \cos \theta$, where $\theta$ is the angle between the initial and final photon directions. This explains why the total SEDs for $v_{\alpha,i} = -400$ and $-800$ km s$^{-1}$ have the shape of a top-hat extending from $v_{\alpha} = v_{\alpha,i}$ to $-v_{\alpha,i}$; $\cos \theta$ is uniformly distributed between 1 and -1 for an isotropic source. The black photon path also demonstrates why the photons around the line center are not captured in the near-zone SED (Figure 9): they are farther away than $1h^{-1}$ cMpc from the source on the sky plane. The photons can be observed in the near-zone only when $\theta$ is close to either 0 or $\pi$ radians, and those are the ones whose emergent wavelength is either $v_{\alpha} \approx v_{\alpha,i}$ or $-v_{\alpha,i}$.

The right panel of Figure 10 describes the $v_{\alpha,i} = 200$ km s$^{-1}$ case, which represents the photons emitted on the red side of Ly$\alpha$ close to the resonant scattering limit ($v_{\alpha,i} = \sqrt{GM_h/r_s}$). In this case, the peculiar infall velocity at the scattering location ($\sim 200$ km s$^{-1}$) is comparable to $v_{\alpha,i}$, and the IGM peculiar motion can significantly affect the wavelength during scattering events. Here, we use a color discontinuity in the arrows to illustrate how the scattering changes the wavelength depending on the initial direction. When a photon initially propagates away from the observer and is scattered toward the observer at the scattering surface (see the blue path in the figure), the infall motion pointing toward the observer blueshifts the photon enough to place it on the blue side in the emergent spectrum. Conversely, a photon initially headed toward the observer would experience a substantial redshift upon scattering (see the red path). The near-zone SED is not much weaker than the total SED in this case because most of the photons are scattered in

the near-zone of the galaxy due to the relatively small scattering surface.

In the case of a small H II bubble (galaxy #0504), the scattering process within the H II region is similar, but a significant fraction of the photons go through additional scatterings in the surrounding H I region. These scattering events increase the photon path and redshift the photons further, suppressing the blue-side SED and instead enhancing the red-side SED. Also, this redward SED is not captured in the near-zone SED because the scattering makes the photons more extended in space. This explains the difference in the SED between cases of large and small H II regions (Figure 9).

### 4.3. Spectrum of Scattered Light: Realistic Source Cases

Next, we consider more realistic cases where the source SEDs have extended profiles. The emergent scattered light SED for an arbitrary emission profile, $F(v_{\alpha})$, can be obtained by superposing the monochromatic source cases:

$$F(v_{\alpha}) = \int_{\min_{v_{\alpha,i}}}^{v_{\alpha,i}^{\max}} F_m(v_{\alpha}|v_{\alpha,i}) W(v_{\alpha,i}) dv_{\alpha,i}, \quad (18)$$

where $F_m(v_{\alpha}|v_{\alpha,i})$ is the emergent SED of the monochromatic source emitting at $v_{\alpha,i}$, and the weight function $W(v_{\alpha,i})$ is given by the intrinsic emission SED of the source galaxy. For each photon, we draw $v_{\alpha,i}$ from a uniform distribution between $v_{\alpha,i}^{\max} = 1000$ km s$^{-1}$ and $v_{\alpha,i}^{\min} = -1000$ km s$^{-1}$ to cover the extent of the typical Ly$\alpha$ emission spectrum of star-forming galaxies. Then, we calculate the scattered light SED for the intrinsic emission model of our choice by weighting each photon by the intrinsic emission profile, $W(v_{\alpha,i})$, in the probability distribution.

In this work, we consider three models for the intrinsic emission. One model has a broad Gaussian peak centered at the
<figure><img src="image_13.png" /><figcaption>**Figure 10.** Schematic of the photon paths for $v_{\alpha,i} = -800$ (left) and 200 km s⁻¹ (right) emitted from galaxy #0001 (lower panels) and the corresponding final wavelengths seen by the observer (upper panels). The scattering location is given by Equation (17) as elaborated in Section 4.1.1. The color gradation in the arrows in the left panel describes the cosmological redshift due to the extra propagation distance from scattering. The color discontinuity in the arrows in the right panels describes the change in frequency during scattering due to the peculiar velocity of the scattering atom. The pair of dashed lines describe the near-zone defined as the region within $r_{\perp} = 1h^{-1}$ Mpc (or 0.6 arcmin) from the source galaxy.</figcaption></figure>

resonance FWHM of $V_1 = 800$ km s⁻¹:

$$W_{\text{CP}}(v_{\alpha,i}) = \exp\left(-\left[\frac{v_{\alpha,i}}{V_1/2.355}\right]^2\right). \quad (19)$$

Here, the factor 2.355 is to make $V_1$ the FWHM of the profile. The other two models are the red- and double-peak models, where we place an off-center Gaussian peak either on the red side only or on both sides. The weight functions are

$$W_{\text{RP}}(v_{\alpha,i}) = \exp\left(-\left[\frac{v_{\alpha,i} - V_2}{V_2/2.355}\right]^2\right) \quad (20)$$

and

$$W_{\text{DP}}(v_{\alpha,i}) = \exp\left(-\left[\frac{v_{\alpha,i} - V_2}{V_2/2.355}\right]^2\right) + \exp\left(-\left[\frac{v_{\alpha,i} + V_2}{V_2/2.355}\right]^2\right), \quad (21)$$

where we assume $V_2 = 300$ km s⁻¹ for the offset and the FWHM of the peaks.

The central peaks model is similar to the results from recent galaxy-scale radiative transfer simulation studies (e.g., Smith et al. 2021). In those simulations, star-forming clouds initially radiate at the $Ly\alpha$ resonance, and the profile is broadened due to the turbulent and rotational motion of the interstellar medium (ISM) within the source galaxy. The red-peak model is supported by observation at $z \lesssim 3$, where the IGM is considered to be transparent to $Ly\alpha$ photons. The outflows in the

circumgalactic medium are known to suppress the blue-side emission from the star-forming ISM (e.g., Yang et al. 2016). The double-peak model is motivated by recent simulation studies suggesting that high-z galaxies may have more porous ISM due to stronger star formation feedback, allowing more blue-side photons to escape the galaxy. The real shape of the intrinsic emission profile is not well constrained today and needs further studies.

We will use the central peak model to demonstrate how the scattered light is processed in the IGM. Then, we compare the results of the red-peak and double-peak models to assess the impact of the blueward emission on the scattered light, which cannot be seen from direct observation of the unscattered light.

### 4.3.1. Central Peak Model

We show the intrinsic emission SED, the unscattered light SED, and the SED of both unscattered and scattered light for galaxies #0001 and #0504 in Figure 11. We show the total SED sampled from the entire volume that we calculated $Ly\alpha$ RT ($r_{\perp} < 8h^{-1}$ Mpc).

Comparing the unscattered light SED to the intrinsic SED shows that the photons with $v_{\alpha,i} \lesssim 250$ (150) km s⁻¹ are completely scattered by the IGM in the case of a large (small) H II region. As detailed in P21, this truncation wavelength is set by the circular velocity of the halo, $V_c = \sqrt{GM_h/r_{200}}$. Above the circular velocity ($v_{\alpha} > V_c$), the unscattered SED converges to the intrinsic SED in the case of a large H II region, but it still remains significantly lower in the small H II region because the
<figure><img src="image_14.png" /><figcaption>**Figure 11.** The intrinsic SED of the central peak emission model (dashed), the emergent SED of the unscattered light directly coming from the source (thin solid), and the emergent SED of both scattered and unscattered light (thick solid). The left panel is for the photons emitted by galaxy #0001, and the right panel is for those emitted by galaxy #0504.</figcaption></figure>

damping-wing opacity of the H I region scatters a fraction of these photons on the red side.

Since we do not consider any absorption by dust in the IGM, the integrated SED of the intrinsic emission is the same as that of the scattered and unscattered light combined. Comparing the two cases shows that the scattered light generally ends up redder than it was at emission. As we described in the monochromatic cases, this is due to the scattering event increasing the path length for the scattered light. In the case of a small H II region, the scattered light adds more to the red side due to the additional scatterings in the H I region.

### 4.3.2. Red- and Double-peak Models

In Figure 12, we compare the scattered light in the red-peak and double-peak source models for cases of large and small H II regions and for total and near-zone SEDs. We note that the unscattered light SED is the same in both models because the blue-side emission is completely scattered by the IGM on its way and does not appear in the unscattered light SED. We also show the intrinsic emission profile of the double-peak model. For the red-peak model, we take the red-side peak of the double-peak model as described by Equations (20) and (21).

In the large H II region (upper panels), the scattered light adds a blueward tail to the combined SED while the unscattered light is truncated near the circular velocity of the halo ($v_\alpha \sim 250 \text{ km s}^{-1}$). In the total SED case (upper left panel), the blue peak makes a significant difference in the SED: the scattered light from the double-peak case adds a thicker and more extended tail (down to $\sim -500 \text{ km s}^{-1}$) than the red-peak case does ($\sim -250 \text{ km s}^{-1}$). However, the blue-peak contribution is much weaker in the near-zone SED (upper right panel), because the blue-side emission forms more extended scattered light in the sky, as we observed in the monochromatic cases. We also repeat this calculation with nonradial emission at $r_{200}$ in the Appendix to confirm that the results do not depend sensitively on the initial photon direction at emission.

In the small H II region (lower panels), the scattered light tends to be distributed over a wider wavelength range due to the additional scatterings in the H I region resulting in more redshifting of the scattered light. The scattered light adds an extended redward tail in the SED well beyond $v_\alpha = 1000 \text{ km s}^{-1}$, as in the monochromatic cases and in the

central peak model. The blueward tail is also present, but it appears much weaker than in the case of a large H II region. Also, the contribution from scattered light in the near-zone SED (lower right panel) is smaller than in the case of a large H II region for both models, because the scatterings in the H I region spread both red- and blue-side photons to a larger patch of the sky.

## 5. Summary and Discussion

We have developed a Monte Carlo Ly$\alpha$ RT simulation code to trace the Ly$\alpha$ photons emitted from high-z galaxies in the intergalactic medium. The code can run on an arbitrary three-dimensional mesh of density, velocity, ionization, and temperature. We do not use any acceleration schemes often adopted in similar works.

We have tested the code for several problems with analytic solutions, including a monochromatic source in a static uniform isothermal slab (Harrington 1973; Neufeld 1990) and sphere (Dijkstra et al. 2006). We also test the code against the results from already published works for simple geometry and kinematics, where a sphere of gas has Hubble-like radially outward (or inward) motion (Zheng & Miralda-Escudé 2002; Tasitsiomi 2006; Dijkstra et al. 2006; Laursen et al. 2009). The emergent spectra and their dependence on the physical parameters agree with the results from the previous works, indicating that our code is reliable.

We ran our code for volumes near two galaxies in the $z=7$ snapshot of the CoDaII simulation to explore the physics of Ly$\alpha$ photon scattering during reionization. Based on the results, we explained how the scattering location is determined by the initial wavelength and cosmological redshift in the H II regions, and how the damping-wing opacity in the H I region can affect the scattering locations. Then, we explained how the geometry of the scattering location is related to the emergent SED of the scattered light. We also present the results for the near-zone of the galaxy (0.56 arcmin from source) in comparison to the total SED within the extent of the simulation (4.5 arcmin from source).

In an H II region, photons propagate freely until they redshift to the Ly$\alpha$ resonance and are scattered by residual neutral hydrogen atoms. The scattering events can increase the frequency by upscattering the photons with the peculiar motion
<figure><img src="image_15.png" /><figcaption>**Figure 12.** The SEDs of scattered light for the red-peak and double-peak models are shown as the red and black solid lines, respectively. The SED of the unscattered light is shown as the thin black solid line, and that of the intrinsic emission is shown as the gray dashed line. The upper panels are for photons emitted from galaxy #0001, and the lower ones are for those from galaxy #0504. The left and right panels show the results for all the sampled photons and the photons within $8h^{-1}$ Mpc and $1h^{-1}$ Mpc or 4.5 and 0.6 arcminutes from the source on the sky, respectively.</figcaption></figure>

of the gas, and also decrease the frequency due to extra cosmological redshift from the increased path length to the observer. Typically, the photons emitted on the blue side of the resonance experience more redshifting because their path length is more dramatically increased by scattering events. Those emitted on the red side, in contrast, face stronger gravitational infall motion of the IGM, giving a larger boost to their frequencies. The blue-side emission makes a significant difference in the emergent spectrum of the scattered light, which potentially allows us to discriminate different intrinsic emission profiles on the blue side, which is impossible from the unscattered light. However, the difference is much smaller if we limit the light collection to the near-zone ($r_{\perp} < 0.6$ arcmin) of the galaxy because the blue-side emission ends up more diffuse and extended in the sky. These findings are broadly consistent with what was reported by Zheng et al. (2010, See Section 4 of their work).

The ionization state of the IGM is another crucial factor. If the surrounding H II region is small ($\lesssim 2h^{-1}$ cMpc), the damping-wing opacity of the nearby H I region becomes significant even for the photons on the red side of the resonance. The photons go through more scattering events in the H I region, resulting in much more spread in both frequency and space. This suggests that the detectability of the scattered

light would steeply drop toward high $z$ as the H II bubbles are expected to be smaller at earlier times.

Our results provide a theoretical framework for interpreting future observations to constrain the properties of the source galaxies. In realistic observations with IFUs, there will be multiple neighboring galaxies within a field of several square arcminutes around the target galaxy, and the light from the neighboring galaxies would mix with that from the target, making the interpretation nontrivial. The collective spatial intensity map of scattered $Ly\alpha$ photons can be studied statistically, e.g., through the intensity power spectrum that could also probe the physical state of the IGM (Visbal & McQuinn 2018). Our study, focused on individual objects, would still be applicable to very bright objects whose scattered $Ly\alpha$ intensity stands out against diffuse background. Further studies may be needed depending on the specifics of interested surveys.

We note that we treat the ISM and CGM as a black box in this work when simulating the photons from the virial radius of the galaxy. Understanding the intrinsic emission exiting the CGM requires dedicated small-scale simulations and is under active investigation by other numerical studies (e.g., Smith et al. 2021). Our results can flexibly accommodate any intrinsic profile from other studies and produce the corresponding scattered light SED using Equation (18).
<figure><img src="image_16.png" /><figcaption>**Figure 13.** The total (left) and near-zone (right) SEDs for nonradial emission. Similarly to Figure 12, we show the results for red-peak and double-peak models by the red and black thick solid lines, and the case of unscattered light only is shown by the thin black line. The dotted lines are from the radial emission cases of Figure 12.</figcaption></figure>

In future work, we plan to extend our analysis to the surface brightness of the scattered light, which is relevant to future intensity mapping surveys such as SPHEREx. We shall also explore the dependence on viewing angle of the scattered light to address the possible variation in observational constraints.

We thank the anonymous referee and A. Smith for helpful comments on this paper. H.P. was supported by the World Premier International Research Center Initiative (WPI), MEXT, Japan and JSPS KAKENHI grant No. 19K23455. The numerical computations of this work were carried out on the gfarm computing cluster of the Kavli Institute for Physics and Mathematics of the Universe and a high-performance computing cluster at the Korea Astronomy and Space Science Institute. H.S. was supported by the Basic Science Research Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Education (2020R111A1A01069228). K.A. was supported by NRF-2016R1D1A1B04935414, 2021R1A2C1095136. and 2016R1A5A1013277. P.R.S. was supported in part by US NSF grant AST-1009799, NASA grant NNX11AE09G, NASA/JPL grant RSA Nos. 1492788 and 1515294, and supercomputer resources from NSF XSEDE grant TG-AST090005 and the Texas Advanced Computing Center (TACC) at the University of Texas at Austin. T.D. was supported by the National Science Foundation Graduate Research Fellowship Program under grant No. DGE-1610403. J.S. acknowledges support from the ANR LOCALI-ZATION project, grant ANR-21-CE31-0019 of the French Agence Nationale de la Recherche. I.T.I. was supported by the Science and Technology Facilities Council (grant Nos. ST/I000976/1 and ST/T000473/1) and the Southeast Physics Network (SEPNet). I.J. acknowledges support from NASA under award number 80GSFC21M0002.

## Appendix
Dependence on Photon Direction at Emission

In this work, we assume all the $Ly\alpha$ photons emitted at $r_{200}$ are in the radial direction, but this is not true in reality, as some photons would be scattered within the CGM and change their directions. Given that $r_{200}$ is small compared to the distance to the first scattering location of most $Ly\alpha$ photons, we do not expect the results to depend sensitively on this assumption. To

test how the results depend on the initial photon direction, we run our $Ly\alpha$ RT for a case of nonradial emission, in which we randomly draw the initial photon direction from the hemisphere pointing in the radial direction. Thus, the actual distribution of the emission direction would lie somewhere between the perfectly radial case and this semi-isotropic (nonradial) case that we test here.

We repeat our calculation with the nonradial emission for the case of a large H II bubble of galaxy #0001 and show the results in Figure 13. For comparison, we show the results from the case with radial emission (solid lines in the upper panels of Figure 12) as a dotted line of the same color and line thickness. Both the total and near-zone SEDs show only a small difference between the two emission cases, confirming that the initial photon direction is a minor factor in the those results.

## ORCID iDs

Hyunbae Park https://orcid.org/0000-0002-7464-7857

Kyungjin Ahn https://orcid.org/0000-0003-3974-1239

Hyunmi Song https://orcid.org/0000-0002-4362-4070

Intae Jung https://orcid.org/0000-0003-1187-4240

Paul R. Shapiro https://orcid.org/0000-0002-0410-3045

Ilian T. Iliev https://orcid.org/0000-0002-5174-1365

## References

* Adams, T. F. 1971, ApJ, 168, 575

* Ahn, S.-H., Lee, H.-W., & Lee, H. M. 2000, JKAS, 33, 29

* Ahn, S.-H., Lee, H.-W., & Lee, H. M. 2001, ApJ, 554, 604

* Ahn, S.-H., Lee, H.-W., & Lee, H. M. 2002, ApJ, 567, 922

* Bunker, A., Smith, J., Spinrad, H., Stern, D., & Warren, S. 2003, Ap&SS, 284, 357

* Camps, P., Behrens, C., Baes, M., Kapoor, A. U., & Grand, R. 2021, ApJ, 916, 39

* Cantalupo, S., Lilly, S. J., & Haehnelt, M. G. 2012, MNRAS, 425, 1992

* Cantalupo, S., Porciani, C., Lilly, S. J., & Miniati, F. 2005, ApJ, 628, 61

* Castellano, M., Pentericci, L., Vanzella, E., et al. 2018, ApJL, 863, L3

* Colbert, J. W., Scarlata, C., Teplitz, H., et al. 2011, ApJ, 728, 59

* Croft, R. A. C., Miralda-Escudé, J., Zheng, Z., Blomqvist, M., & Pieri, M. 2018, MNRAS, 481, 1320

* Curtis-Lake, E., McLure, R. J., Pearce, H. J., et al. 2012, MNRAS, 422, 1425

* Dijkstra, M. 2014, PASA, 31, e040

* Dijkstra, M. 2017, arXiv:1704.03416

* Dijkstra, M., Haiman, Z., & Spaans, M. 2006, ApJ, 649, 14

* Dijkstra, M., & Loeb, A. 2009, MNRAS, 400, 1109

* Endsley, R., Stark, D. P., Charlot, S., et al. 2021, MNRAS, 502, 6044
Faucher-Giguère, C.-A., Kereš, D., Dijkstra, M., Hernquist, L., & Zaldarriaga, M. 2010, ApJ, 725, 633

Fernandez, E. R., & Komatsu, E. 2006, ApJ, 646, 703

Fontana, A., Vanzella, E., Pentericci, L., et al. 2010, ApJL, 725, L205

Francis, P. J., Woodgate, B. E., Warren, S. J., et al. 1996, ApJ, 457, 490

Gronke, M., Ocvirk, P., Mason, C., et al. 2021, MNRAS, 508, 3697

Haiman, Z., Spaans, M., & Quataert, E. 2000, ApJL, 537, L5

Hansen, M., & Oh, S. P. 2006, MNRAS, 367, 979

Harikane, Y., Ouchi, M., Ono, Y., et al. 2019, ApJ, 883, 142

Harrington, J. P. 1973, MNRAS, 162, 43

Higuchi, R., Ouchi, M., Ono, Y., et al. 2019, ApJ, 879, 28

Hu, W., Wang, J., Infante, L., et al. 2021, NatAs, 5, 485

Hutter, A., Dayal, P., Partl, A. M., & Müller, V. 2014, MNRAS, 441, 2861

Iliev, I. T., Shapiro, P. R., McDonald, P., Mellema, G., & Pen, U.-L. 2008, MNRAS, 391, 63

Jeeson-Daniel, A., Ciardi, B., Maio, U., et al. 2012, MNRAS, 424, 2193

Jung, I., Finkelstein, S. L., Dickinson, M., et al. 2019, ApJ, 877, 146

Jung, I., Finkelstein, S. L., Dickinson, M., et al. 2020, ApJ, 904, 144

Jung, I., Papovich, C., Finkelstein, S. L., et al. 2021, arXiv:2111.14863

Kakuma, R., Ouchi, M., Harikane, Y., et al. 2021, ApJ, 916, 22

Katz, H., Galligan, T. P., Kimm, T., et al. 2019, MNRAS, 487, 5902

Kikuchihara, S., Harikane, Y., Ouchi, M., et al. 2021, arXiv:2108.09288

Kim, E., Yang, Y., Zabludoff, A., et al. 2020, ApJ, 894, 33

Laursen, P., Razoumov, A. O., & Sommer-Larsen, J. 2009, ApJ, 696, 853

Leclercq, F., Bacon, R., Wisotzki, L., et al. 2017, A&A, 608, A8

Loeb, A., & Rybicki, G. B. 1999, ApJ, 524, 527

Malhotra, S., & Rhoads, J. E. 2004, ApJL, 617, L5

Mallery, R. P., Mobasher, B., Capak, P., et al. 2012, ApJ, 760, 128

Mason, C. A., Treu, T., de Barros, S., et al. 2018, ApJL, 857, L11

Matsuda, Y., Yamada, T., Hayashino, T., et al. 2004, AJ, 128, 569

Matsuda, Y., Yamada, T., Hayashino, T., et al. 2012, MNRAS, 425, 878

Matthee, J., Sobral, D., Santos, S., et al. 2015, MNRAS, 451, 400

Neufeld, D. A. 1990, ApJ, 350, 216

Ocvirk, P., Aubert, D., Sorce, J. G., et al. 2020, MNRAS, 496, 4087

Ocvirk, P., Gillet, N., Shapiro, P. R., et al. 2016, MNRAS, 463, 1462

Ono, Y., Ouchi, M., Mobasher, B., et al. 2012, ApJ, 744, 83

Ouchi, M., Shimasaku, K., Furusawa, H., et al. 2010, ApJ, 723, 869

Park, H., Jung, I., Song, H., et al. 2021, ApJ, 922, 263

Partridge, R. B., & Peebles, P. J. E. 1967, ApJ, 147, 868

Pentericci, L., Fontana, A., Vanzella, E., et al. 2011, ApJ, 743, 132

Qin, Y., Wyithe, J. S. B., Oesch, P. A., et al. 2022, MNRAS, 510, 3858

Rosdahl, J., & Blaizot, J. 2012, MNRAS, 423, 344

Sadoun, R., Zheng, Z., & Miralda-Escudé, J. 2017, ApJ, 839, 44

Schaerer, D., & Verhamme, A. 2008, A&A, 480, 369

Semelin, B., Combes, F., & Baek, S. 2007, A&A, 474, 365

Seon, K.-I., Song, H., & Chang, S.-J. 2022, ApJS, 259, 3

Smith, A., Kannan, R., Garaldi, E., et al. 2022, MNRAS, 512, 3243

Smith, A., Kannan, R., Tacchella, S., et al. 2021, arXiv:2111.13721

Stark, D. P., Ellis, R. S., & Ouchi, M. 2011, ApJL, 728, L2

Steidel, C. C., Adelberger, K. L., Shapley, A. E., et al. 2000, ApJ, 532, 170

Tasitsiomi, A. 2006, ApJ, 645, 792

Tilvi, V., Malhotra, S., Rhoads, J. E., et al. 2020, ApJL, 891, L10

Tilvi, V., Papovich, C., Finkelstein, S. L., et al. 2014, ApJ, 794, 5

Treu, T., Schmidt, K. B., Trenti, M., Bradley, L. D., & Stiavelli, M. 2013, ApJL, 775, L29

Verhamme, A., Schaerer, D., Atek, H., & Tapken, C. 2008, A&A, 491, 89

Verhamme, A., Schaerer, D., & Maselli, A. 2006, A&A, 460, 397

Visbal, E., & McQuinn, M. 2018, ApJL, 863, L6

Wisotzki, L., Bacon, R., Brinchmann, J., et al. 2018, Natur, 562, 229

Yajima, H., Li, Y., Zhu, Q., & Abel, T. 2012, MNRAS, 424, 884

Yang, H., Malhotra, S., Gronke, M., et al. 2016, ApJ, 820, 130

Yang, Y., Zabludoff, A., Jahnke, K., & Davé, R. 2014, ApJ, 793, 114

Zheng, Z., Cen, R., Trac, H., & Miralda-Escudé, J. 2010, ApJ, 716, 574

Zheng, Z., Cen, R., Weinberg, D., Trac, H., & Miralda-Escudé, J. 2011, ApJ, 739, 62

Zheng, Z., & Miralda-Escudé, J. 2002, ApJ, 578, 33

Zheng, Z.-Y., Wang, J., Rhoads, J., et al. 2017, ApJL, 842, L22