# Seurat: SPH scheme extended with ultraviolet line radiative transfer

Makito Abe, Hiroyuki Suzuki, Kenji Hasegawa, Benoît Semelin, Hidenobu

Yajima, Masayuki Umemura

► To cite this version:

Makito Abe, Hiroyuki Suzuki, Kenji Hasegawa, Benoît Semelin, Hidenobu Yajima, et al. Seurat: SPH scheme extended with ultraviolet line radiative transfer. Monthly Notices of the Royal Astronomical Society, 2018, 476 (2), pp.2664-2673. 10.1093/mnras/sty233. hal-02345604

HAL Id: hal-02345604

https://hal.science/hal-02345604

Submitted on 15 May 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished 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
# SEURAT: SPH scheme extended with ultraviolet line radiative transfer

Makito Abe,¹* Hiroyuki Suzuki,¹ Kenji Hasegawa,² Benoit Semelin,³ Hidenobu Yajima⁴,⁵ and Masayuki Umemura¹

¹Center for Computational Sciences, University of Tsukuba, Ten-nodai, 1-1-1 Tsukuba, Ibaraki 305-8577, Japan

²Graduate School of Science, Nagoya University, Furo-cho, Chikusa-ku, Nagoya, Aichi 464-8602, Japan

³LERMA, Observatoire de Paris, Sorbonne Université, PSL research university, CNRS, F-75014 Paris, France

⁴Frontier Research Institute for Interdisciplinary Sciences, Tohoku University, Sendai 980-8578, Japan

⁵Astronomical Institute, Tohoku University, Sendai 980-8578, Japan

Accepted 2018 January 25. Received 2018 January 24; in original form 2017 May 9

## ABSTRACT

We present a novel Lyman alpha (Ly$\alpha$) radiative transfer code, SEURAT (SPH scheme Extended with Ultraviolet line RAdiative Transfer), where line scatterings are solved adaptively with the resolution of the smoothed particle hydrodynamics (SPH). The radiative transfer method implemented in SEURAT is based on a Monte Carlo algorithm in which the scattering and absorption by dust are also incorporated. We perform standard test calculations to verify the validity of the code; (i) emergent spectra from a static uniform sphere, (ii) emergent spectra from an expanding uniform sphere, and (iii) escape fraction from a dusty slab. Thereby, we demonstrate that our code solves the Ly$\alpha$ radiative transfer with sufficient accuracy. We emphasize that SEURAT can treat the transfer of Ly$\alpha$ photons even in highly complex systems that have significantly inhomogeneous density fields. The high adaptivity of SEURAT is desirable to solve the propagation of Ly$\alpha$ photons in the interstellar medium of young star-forming galaxies like Ly$\alpha$ emitters (LAEs). Thus, SEURAT provides a powerful tool to model the emergent spectra of Ly$\alpha$ emission, which can be compared to the observations of LAEs.

**Key words:** hydrodynamics – line: profile – radiative transfer – methods: numerical – galaxies: high-redshift.

## 1 INTRODUCTION

The hydrogen Lyman alpha (Ly$\alpha$) emission is a significant probe to explore high-z young galaxies, as pioneeringly suggested by Partridge & Peebles (1967). So far, a number of distant galaxies have been detected via intensive observations of the Ly$\alpha$ line, the so-called Ly$\alpha$ emitters (LAEs: Iye et al. 2006; Ouchi et al. 2008, 2010; Vanzella et al. 2011; Ono et al. 2012; Shibuya et al. 2012; Finkelstein et al. 2013; Konno et al. 2014; Zitrin et al. 2015). Recently, the most distant LAE has been detected at $z = 8.68$ by Zitrin et al. (2015). It is difficult to reveal the detailed structure and physical properties of distant galaxies even with state-of-the-art observational facilities. However, the Ly$\alpha$ line can, in principle, provide copious information about the internal structure of those galaxies, because the emergent spectra of Ly$\alpha$ line depend on velocity fields and ionization structure in the galaxies (Dijkstra, Haiman & Spaans 2006a,b; Verhamme, Schaerer & Maselli 2006; Yajima et al. 2012a, 2015) and the photon escape fraction is sensitive to the metallicity of the gas (Atek et al. 2008; Verhamme et al. 2008; Yajima et al. 2014).

In addition, the neutral hydrogen in the intergalactic medium (IGM) erodes the Ly$\alpha$ emission from high-z galaxies. Consequently, the cosmic reionization history is constrained by investigating the redshift evolution of the luminosity function (LF) of LAEs (e.g. Kashikawa et al. 2006; Ouchi et al. 2010). Moreover, LAEs can be the main ionizing sources responsible for the cosmic reionization (Yajima et al. 2009, 2014). It also should be noted that the shape of Ly$\alpha$ line emission can constrain the neutral fraction of the IGM by considering the distortion due to the IGM transmission (Santos 2004; Dijkstra, Lidz & Wyithe 2007). Therefore, understanding the Ly$\alpha$ emission from high-redshift galaxies is crucial for revealing not only galaxy formation but also cosmic reionization.

Ly$\alpha$ photons are emitted by the 2P–1S transitions as a result of the excitation of H$_1$ atoms or the recombination in H$_2$ regions (Faucher-Giguère et al. 2010; Yajima et al. 2012b). Due to the large cross-section of Ly$\alpha$ radiation to neutral hydrogen, the interstellar medium (ISM) of LAEs can be readily optically thick to Ly$\alpha$ photons. In general, the transfer in a spectral line arising from the spontaneous bound–bound transition can be treated as resonant scattering, if the transition time-scale is much shorter than the other physical time-scales (e.g. collisional deexcitation) and also the decaying to some other state is negligible. Since the Ly$\alpha$ transition with

<footnote>* E-mail: mabe@ccs.tsukuba.ac.jp</footnote>
a large Einstein's A coefficient satisfies these conditions in the ISM of LAEs, we can treat the absorption and subsequent re-emission of a Lyα photon as a resonant scattering process. When scatterings are included, the radiative transfer equation takes the form of an integro-differential equation that should be solved through an iterative procedure. Moreover, we should pay special attention to the partial frequency redistribution during line scatterings, which is essential in Lyα radiative transfer. Therefore, in order to obtain an exact solution of the Lyα radiation transfer equation, we have to solve iteratively the frequency-dependent radiative transfer equation. Needless to say, it is generally time consuming to integrate such an equation directly. Under the diffusion approximation, analytical solutions for the emergent spectrum were derived for systems with simple geometry composed of pure hydrogen, such as a static and uniform slab or a spherical cloud (Harrington 1973; Neufeld 1990; Dijkstra et al. 2006a). On the other hand, for inhomogeneous or moving gas clouds, Lyα properties have to be calculated by numerical simulations.

A Monte Carlo approach, in which various phenomena are stochastically treated, is often employed to manage the complexity of the partial frequency redistribution in Lyα line transfer (e.g. Zheng & Miralda-Escudé 2002; Verhamme et al. 2006; Semelin, Combes & Baek 2007; Baek et al. 2009; Laursen, Razoumov & Sommer-Larsen 2009; Yajima et al. 2012a; Smith et al. 2015; Yajima & Li 2014; Yajima et al. 2014). Some previous works investigated the Lyα radiation properties of galaxies by using simple models (e.g. Dijkstra, Haiman & Spaans 2006a; Dijkstra, Mesinger & Wyithe 2011; Verhamme, Schaerer & Maselli 2006; Verhamme et al. 2008; Gronke, Bull & Dijkstra 2015). For example, Verhamme et al. (2008) solved the Lyα radiation transfer in spherical expanding shells and reproduced the Lyα properties of observed galaxies with tuned H I column densities and expansion velocities (see also Dijkstra et al. 2011). However, the ISM can have a complicated structure due to stellar feedback, instability of the galactic disc, and interaction with other galaxies. Recent observations with high angular resolution have revealed the inhomogeneous and clumpy ISM of high-z galaxies (e.g. Genzel et al. 2011). Therefore, a solver that allows us to treat Lyα radiation transfer in more complex ISM structures is required.

Recent cosmological simulations have been able to model galaxies and resolve their internal ISM structure (e.g. Wise et al. 2012; Hasegawa & Semelin 2013; Hopkins et al. 2014; Vogelsberger et al. 2014; Schaye et al. 2015). Most of simulations have been conducted using the smoothed particle hydrodynamics (SPH), which is a Lagrangian numerical scheme (e.g. see a review by Springel 2010). An important advantage of SPH is that the spatial resolution automatically augments with increasing local density. Therefore, the SPH method can resolve adaptively the ISM and star-forming regions, which produce Lyα photons. So far, Monte Carlo simulations of Lyα radiative transfer have been combined with SPH by mesh-based schemes, where the SPH densities are assigned on the meshes before solving the Lyα radiative transfer (e.g. Semelin et al. 2007; Yajima et al. 2012a). Laursen et al. (2009), combining cosmological SPH simulations with post-processing transfer on the meshes, investigated the Lyα properties of high-z galaxies and showed that the Lyα flux changes by a factor of 3–6 depending on the viewing angles due to the complex ISM structure (see also, Yajima et al. 2012b). Yajima et al. (2015) calculated the Lyα properties of high-z progenitors of a local Milky Way-like galaxy, and showed that the filamentary accreting gas produces Lyα photons efficiently via the excitation cooling process at $z > 6$. However, the mesh-based schemes of Lyα radiative transfer inevitably smooth out highly

resolved structures in SPH simulations, and might lead to a lack of accuracy in the Lyα photon transfer. In this paper, we develop a novel Monte Carlo scheme for Lyα radiative transfer adapted for SPH simulations, SEURAT (SPH scheme Extended with Ultraviolet line RAdiative Transfer). Some basic parts of SEURAT are the extension of the mesh-based Lyα radiation transfer code LICORICE (Semelin et al. 2007). In SEURAT, SPH particles themselves are directly used to solve the radiation transfer unlike in the previous mesh-based codes. Hence, the transfer of Lyα photons can be pursued without reducing the resolution of SPH simulations.

This paper is organized as follows. In Section 2, we describe the basic physics of Lyα radiation transfer and the algorithm of SEURAT. The results of some standard test calculations are presented in Section 3. In Section 4, we demonstrate the adaptivity of SEURAT for highly inhomogeneous media. Also, we compare the Lyα transfer with SEURAT to that with LICORICE for a model galaxy obtained by a cosmological radiation SPH simulation. Section 5 is devoted to the conclusions.

## 2 THE CODE

We first describe the flow of our Monte Carlo radiative transfer calculations. Once a photon packet (a monochromatic group of photons) is sent from a radiation source, the packet flies in a straight line along a stochastically determined direction until an interaction by an atom occurs (scattering or absorption). The probability distribution function of the path length that a packet can travel without an interaction is $exp(-\tau)$, where $\tau$ is the optical depth of the path. In order to treat radiative interaction events stochastically, we randomly chose an optical depth $\tau_{\text{target}}$, where an interaction takes place, as $\tau_{\text{target}} = -\ln\xi$, with $\xi$ being a uniform random number between 0 and 1. We then integrate the optical depth along the light ray and let the photon packet propagate until $\tau = \tau_{\text{target}}$ as described later in Section 2.2. When the packet reaches $\tau = \tau_{\text{target}}$, we stochastically determine what happens there. Only scatterings occur in the pure hydrogen gas, whereas both absorption and scattering take place in dusty gas. If a scattering occurs at $\tau = \tau_{\text{target}}$, we settle the resultant direction and frequency shift of the scattered photon packet following the method described in Section 2.1, and chose a new target optical depth. On the other hand, if the packet is absorbed, we completely eliminate the packet, or reduce the photon flux by a factor of attenuation. We repeat this flow until the photon packet escapes or vanishes from the system (Sections 2.2 and 2.3).

### 2.1 Treatment of frequency shifts during the Ly α scattering process

The broadening effect due to the quantum uncertainty results in the Lorentzian line profile $\phi(\nu)$ described as

$$ \phi(\nu) = \frac{\Delta v_L / 2\pi}{(\nu - v_0)^2 + (\Delta v_L / 2)^2}, \quad (1) $$

where $\Delta v_L = 9.936 \times 10^7$ Hz and $v_0 = 2.466 \times 10^{15}$ Hz are the natural broadening width and the central frequency of the Ly α line, respectively. The cross-section $\sigma_ν$ of the Lyα scattering in the rest frame of hydrogen atom is

$$ \sigma_{\nu} = f_{12} \frac{\pi e^2}{m_e c} \phi(\nu) = f_{12} \frac{\pi e^2}{m_e c} \frac{\Delta v_L / 2\pi}{(\nu - v_0)^2 + (\Delta v_L / 2)^2}, \quad (2) $$

where $f_{12} = 0.4162$ is the Lyα oscillator strength, $m_e$ is the mass of electron, and $c$ is the speed of light. When we consider the motion of an atom, the line centre frequency shifts to $v_0(1 + v_{\parallel}/c)$ due
to the Doppler shift in the laboratory frame, where $v_\parallel$ denotes the velocity component of the atom parallel to the incident direction of a photon. Therefore, if the atoms have a thermal velocity distribution function (i.e. Maxwellian), the line profile in the laboratory frame is the well-known Voigt profile and can be obtained by summing up the Lorentzian profiles with various velocities;

$$ \phi_H(v) = \frac{\Delta v_L}{2\pi} \int_{-\infty}^{\infty} dv_\parallel \frac{(m_H/2\pi k_B T)^{1/2} \exp(-m_H v_\parallel^2/2k_B T)}{(v - v_0 - v_0 v_\parallel/c)^2 + (\Delta v_L/2)^2}, \quad (3) $$

where T is the temperature, $k_B$ is the Boltzmann constant, and $m_H$ is the neutral hydrogen mass. A normalized frequency, $x = (v - v_0)/\Delta v_D$, is often used for the frequency distribution function, where $\Delta v_D = (v_{th}/c)v_0$ denotes the thermal Doppler broadening width, and $v_{th} = \sqrt{2k_B T/m_H}$ corresponds to the thermal velocity dispersion. With these expressions, the Ly$\alpha$ scattering cross-section is described as

$$ \sigma_{H,x} = f_{12} \frac{\sqrt{\pi} e^2}{m_e c \Delta v_D} H(a, x), \quad (4) $$

where $a = \Delta v_L/(2\Delta v_D)$ is the relative line width. $H(a, x)$ is the Voigt function defined as

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

The Voigt profile is composed of a Gaussian core and power-law wings. Although the Voigt function cannot be integrated analytically, Tasitsiomi (2006) has provided a useful analytical fitting formula for the Voigt function in the form

$$ H(a, x) \sim q\sqrt{\pi} + e^{-x^2}, \quad (6) $$

where

$$ \begin{aligned} q &= \begin{cases} 0 & \text{for } z \le 0 \\ (1 + \frac{21}{x^2}) \frac{a}{\pi(x^2+1)} & \text{for } z > 0 \end{cases} \\ P(z) &= 5.674z^4 - 9.207z^3 + 4.421z^2 + 0.1117z, \\ z &= (x^2 - 0.855)/(x^2 + 3.42). \end{aligned} $$

We adopt this formula to evaluate the Ly$\alpha$ scattering cross-section in our code. It is worth mentioning that the cross-section in the wing is typically ~10⁵ times smaller than that at the line centre. Hence, once a frequency is shifted to the wing part, the photon can easily escape from a system even if the system is quite optically thick at the line centre frequency.

If a photon is scattered by an atom in the direction parallel to the incident one, the frequency of the photon does not vary, because the scattering is coherent in the atom's rest frame. However, if the scattering direction is not parallel, the frequency is shifted owing to the motion of the atom. This results in the partial frequency redistribution. We denote the incident frequency of the photon in the laboratory frame by $v_{in}$, its incoming direction by $d$, the atom's microscopic velocity by $v$, and the bulk (macroscopic) velocity of the fluid by $V$. Using these expressions, the incident frequency in the atom rest frame $v_{atom}$ is written as

$$ v_{\text{atom}} = v_{\text{in}} \left( 1 - \frac{\mathbf{v} \cdot \mathbf{d}}{c} - \frac{\mathbf{V} \cdot \mathbf{d}}{c} \right). \quad (7) $$

After the scattering, the photon flies in a different direction $d'$. Since the scattering is coherent in the rest frame of the atom, the frequency of the scattered photon in the laboratory frame $v_{out}$ turns out to be

$$ v_{\text{out}} = v_{\text{atom}} / \left( 1 - \frac{\mathbf{v} \cdot \mathbf{d}'}{c} - \frac{\mathbf{V} \cdot \mathbf{d}'}{c} \right). \quad (8) $$

It should be noted that the scattering direction is actually dependent on the phase function, which is determined by the excitation state and scattering frequency (e.g. see Tasitsiomi 2006). However, for an optically thick medium, the anisotropy of the phase function is not significant for the emergent spectrum and/or the escape fraction (Semelin et al. 2007; Laursen et al. 2009; Yajima et al. 2012a). Therefore, we assume a simple isotropic scattering in our code.

Since the frequency shift is determined by the velocity component in the scattering direction, its probability depends on the shape of the velocity distribution function. To quantify the frequency shift in our code, we split the velocity components of the atoms into $v_\parallel$ and $v_\perp$ that are, respectively, the perpendicular and parallel components with respect to $d$. For the perpendicular component $v_\perp$, the velocity distribution function is a Maxwellian. Therefore, the distribution function of the perpendicular component normalized by the thermal velocity dispersion ($u_\perp \equiv v_\perp/v_{th}$) is given by

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

On the other hand, for a given incoming relative frequency $x$, the scattering cross-section is the superposition of the Lorentz profiles shifted in terms of the parallel velocities of different atoms, which is described by the Voigt function (5). Hence, the probability that a photon is scattered by an atom with a certain parallel velocity $u_\parallel(\equiv v_\parallel/v_{th})$ is given by

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

In the case of $|x| \ll 1$, photons are selectively scattered by atoms with $u_\parallel \approx x$. As a result, the probability function $f(u_\parallel)$ exhibits a sharp peak at $u_\parallel \approx x$. On the other hand, for large $|x|$, the number of atoms with $u_\parallel \approx x$ exponentially decreases and therefore photons are scattered predominantly in the Lorentz wing. In this case, the distribution function behaves as a Gaussian distribution (e.g. see fig. 2 of Laursen et al. 2009).

In the present code, we randomly generate the perpendicular components $u_\perp$ by the Box-Muller method,

$$ \begin{aligned} u_{\perp,1} &= \sqrt{-\ln R_1} \cos(2\pi R_2) \\ u_{\perp,2} &= \sqrt{-\ln R_1} \sin(2\pi R_2), \end{aligned} \quad (11) $$

where $R_1$ and $R_2$ are two univariates. The parallel component $u_\parallel$ is randomly generated to obey the distribution function (10) in the same way as in Zheng & Miralda-Escudé (2002).

## 2.2 Ray-tracing in SEURAT

The method for integrating optical depths in our code is fairly distinct from those in the previous Monte Carlo schemes of line transfer, since our ray-tracing algorithm is optimized for using SPH particles directly. Ray-tracing schemes of continuum radiative transfer optimized for SPH simulations have hitherto succeeded in handling a wide dynamic range by the Lagrangian description of SPH (e.g. Kessel-Deynet & Burkert 2000; Susa 2006; Altay, Croft & Pelupessy 2008; Pawlik & Schaye 2008; Hasegawa & Umemura 2010).

The ray-tracing algorithm in SEURAT is similar to that in Susa (2006) except that rays are not always cast towards SPH particles. Fig. 1 shows the principle of the ray-tracing algorithm in SEURAT. After determining the direction of the light ray from a radiation source, we search the particle closest to the light ray (particle 1 in Fig. 1) from the list of neighbour particles of the source. Then, we
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Schematic illustration of the ray-tracing scheme in SEURAT. The filled circles and numbered open circles represent the SPH particles. In the ray-tracing, we search the downstream particle which is the closest to the ray. Grid points are determined by the projection of the downstream particles, as shown by red dots, and the optical depth is evaluated at intermediate points between two grid points. We do not use particle 1' as a downstream particle since particle 1 is closer to the ray, although particle 1' has the smaller angle. The dashed sectors represent the region where candidates for the next downstream particle are searched. Note that the radii of the sectors are not necessarily equal to the SPH smoothing length (see Section 4).</figcaption></figure>

create a virtual grid point at the position where the perpendicular line intersects with the ray, and directly assign the physical quantities of the particle such as the density, the temperature, and the velocity component ($d \cdot v$) to the grid point (grid point $s_1$ in Fig. 1). We define the next downstream grid point (grid point $s_2$ in Fig. 1) on the light ray in the same way except that we search from the list of neighbour particles of particle 1. The optical depth is then integrated with the following formula,

$$ \tau_{i+1/2} = \tau_{i-1/2} + \Delta\tau_i, \qquad (12) $$

$$ \Delta\tau_i = \sigma_{x_i} n(r_i) (s_{i+1/2} - s_{i-1/2}), \qquad (13) $$

where $\sigma_{x_i}$ is the scattering cross-section of a Ly$\alpha$ photon at a frequency $x_i$ and $n(r_i)$ is the number density at a particle position $r_i$ that is obtained by the superposition of neighbouring SPH particles as

$$ n(r_i) = \sum_j m_j W(|r_i - r_j|, h_i)/m_H, \qquad (14) $$

where $m_j, h_i$, and $W$, respectively, denote the mass of the jth neighbour particle, the smoothing length of the ith particle, and the kernel function. The accuracy of our ray-tracing scheme will be discussed in Section 4.4. The $s_{i+1/2}$ is the position defined as the intermediate point between the grid points $s_i$ and $s_{i+1}$ on the light ray. Here, we assume that the physical quantities are constant between the back and front of the ith particle. This formulation is employed to avoid the mixing of the physical quantities between adjacent SPH particles. Especially for velocities, if adjacent SPH particles move with the similar speed in opposite directions, the mixing results in canceling out the velocities and underestimating the $d \cdot v$ term artificially. Repeating calculations described by equations (12) and (13) along the light ray, we integrate the optical depth while $\tau_{i+1/2} < \tau_{\text{target}}$. If the condition of $\tau_{i+1/2} > \tau_{\text{target}}$ is satisfied, we determine the position $s_{\text{target}}$, where a scattering occurs by using the physical quantities of ith particle (in Fig. 1, the ith particle corresponds to the particle 3), i.e.

$$ \tau_{\text{target}} = \tau_{i-1/2} + \sigma_{x_i} n(r_i)(s_{\text{target}} - s_{i-1/2}). \qquad (15) $$

On scattering at $s_{\text{target}}$, we stochastically choose a new direction $d'$ to which the scattered photon packet travels, and let the packet propagate until a subsequent scattering occurs. If there is no neighbour

downstream particle during the ray-tracing, the position is regarded as a boundary of the system, and we assume that the photon packet escapes from the system.

It is worth mentioning that the scattering condition $\tau_{i+1/2} > \tau_{\text{target}}$ is often satisfied before the photon packet renews the downstream particle, since the local optical depth is very large. In this case, to ensure continuity of the physical quantities (e.g. velocities) before and after the scattering, the same particle $i$ is used as the origin of the subsequent ray-tracing rather than the particle closest to the packet. Consequently, multiple scattering events frequently occur around a single particle and the photon packet tends to remain within its smoothing kernel length. We will discuss how this behaviour affects to the SPH-based ray-tracing in Section 4.

## 2.3 Dust absorption and scattering

The absorption and scattering by dust grains make a significant effect on the Ly$\alpha$ radiative transfer in metal-enriched media. The optical depth for the dust scattering and absorption can be written as

$$ d\tau_{\text{d},x} = (Q_{s,x} + Q_{a,x})n_d\pi r_d^2 ds \equiv d\tau_{\text{d},s,x} + d\tau_{\text{d},a,x}, \qquad (16) $$

with

$$ n_d = n_H \left( \frac{m_H}{m_d} \right) f_d, \qquad (17) $$

where $r_d$ is the dust grain radius, $Q_{s,x}$ and $Q_{a,x}$, respectively, denote the scattering and absorption efficiencies (so-called Q-value) at a frequency $x$, and $n_d$ is the dust density that is determined by the dust grain mass $m_d$ and the dust-to-gas mass ratio $f_d$. In the range of ultraviolet frequencies, both $Q_{s,x}$ and $Q_{a,x}$ are ~1 (Verhamme et al. 2006), i.e. they are independent of the photon frequency. Then, the total optical depth $d\tau_{\text{tot},x}$ including contributions by hydrogen and dust is given by

$$ d\tau_{\text{tot},x} = (n_{H_i}\sigma_{\text{H},x} + n_d\sigma_d) ds \equiv d\tau_{\text{H},x} + d\tau_{\text{d},s} + d\tau_{\text{d},a}, \qquad (18) $$

where $d\tau_{\text{H},x}$, $d\tau_{\text{d},s}$, and $d\tau_{\text{d},a}$ correspond to the optical depths for the hydrogen scattering, the dust scattering, and the dust absorption, respectively.

The probability that Ly$\alpha$ photons are scattered by hydrogen is given by

$$ p_{\text{H},x} = \frac{n_{H_i}\sigma_{\text{H},x}}{n_{H_i}\sigma_{\text{H},x} + n_d\sigma_d}. \qquad (19) $$

We generate a uniform random number $\xi'$ between 0 and 1, and let the photon packet interact with hydrogen if $\xi' < p_{\text{H},x}$ and otherwise interact with dust.

We invoke two different methods to assess the absorption by dust: one is the 'stochastic elimination method' and the other is the 'flux attenuation method'. When a photon packet interacts with dust, the fraction of scattering is given by the scattering albedo, $\sigma = Q_s/(Q_s + Q_a)$. In the stochastic elimination method, we generate another uniform random number $\xi''$, and eliminate the photon packet if $\xi'' < 1 - \sigma$. On the other hand, in the flux attenuation method, we reduce the number of photons in a photon packet in accordance with an attenuation factor $\exp(-\tau_{d,a})$ during the journey. To avoid double counting in the latter case, we neglect the third term $d\tau_{d,a}$ in equation (18) when we integrate $d\tau_{\text{tot},x}$ up to $\tau_{\text{target}}$.

## 3 TESTS OF THE CODE

In this section, we present several test calculations to demonstrate the validity of the code. Throughout the tests, we model the systems
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Emergent spectra from a static homogeneous spherical cloud for three different optical depths. The line centre optical depth of the cloud is shown by the attached number. The temperature of the clouds is assumed to be 10 K, and the photon source is located at the centre of the cloud. Solid lines represent the results of the numerical calculations, whereas dotted lines denote the analytical solutions derived by Dijkstra et al. (2006a).</figcaption></figure>

with 64³ SPH particles, and pursue the propagation of 10⁵ photon packets emitted at the Lyα line centre frequency. We note that it is unsuitable to define a geometrically sharp boundary of a system in the SPH formalism, since the gas density fields associated to the SPH particles are smoothed with the kernel function, even though a simple geometry such as a spherical cloud or a slab is assumed. Thus, for the test calculations, we regard photons beyond an assumed boundary as escaped ones. Note that previous Monte Carlo codes often adopted the core-skipping method, which artificially avoids the significant number of the scattering events that happen in the core of the line and reduces the computational time (Ahn, Lee & Lee 2002). However, we do not use the core-skipping scheme in this paper, since accelerating the computation is not our aim for now. We focus here on proving the validity of the mesh-free Monte Carlo Lyα radiation transfer calculation.

## 3.1 Static homogeneous spherical cloud

For a static homogeneous spherical cloud, the analytical formula of the emergent spectrum was derived by Dijkstra et al. (2006a) as follows:

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

where $\tau_0$ denotes the line centre optical depth from the centre to the boundary of the cloud. We test three cases of optical depth, $\tau_0 = 10^4, 10^5$, and $10^6$, for a spherical cloud with the gas temperature of 10 K. SPH particles are randomly distributed in the cloud. The source of Lyα photons is located at the centre of the cloud. Fig. 2 shows the comparison between the analytical solutions of Dijkstra et al. (2006a) and our numerical calculations. We can see in Fig. 2 that the numerically calculated emergent spectra (solid lines) are in good agreement with the analytical solutions (dotted lines). Strictly speaking, the agreement between the numerical and analytical emergent spectra tends to be better at large optical depth, since the analytic formula is valid for extremely optically thick media, say, $a\tau_0 \gtrsim 10^3$ (Neufeld 1990). In this test, $T = 10$ K corresponds to

<figure><img src="image_4.png" /><figcaption>**Figure 3.** Emergent spectra from an expanding spherical cloud. The temperature of the cloud is assumed to be 20 000 K, and three different column densities of hydrogen are assumed. The radial velocity is proportional to radius $r$ with the velocity of 200 km s⁻¹ at the edge.</figcaption></figure>

$a \sim 1.5 \times 10^{-2}$. Thus, the agreement for the corresponding $\tau_0 = 10^4$ spectrum is slightly worse when compared to more optically thick cases. The resultant spectra reproduce the double-peaked shapes, and the positions of the peak frequencies move outward with increasing the optical depth of the system. This is because only Lyα photons with frequencies far from the line centre can escape from the clouds as $\tau_0$ increases.

## 3.2 Expanding homogeneous spherical cloud

Next we calculate the emergent Lyα spectrum from an expanding spherical cloud. For such a moving medium, there is no analytical solution. However, for a homogeneous expanding/infalling cloud, it is comparatively easy to formulate a physical interpretation of the emergent spectrum. Hence, this test is useful to demonstrate the validity of the code, and has been used in previous works (e.g. Zheng & Miralda-Escudé 2002; Dijkstra et al. 2006a; Verhamme et al. 2006; Semelin et al. 2007; Laursen et al. 2009; Yajima et al. 2012a). We consider a uniform spherical cloud with a Hubble-like velocity field

$$v(r) = v_{\max} \left( \frac{r}{r_{\max}} \right), \quad (21)$$

where $v_{\max}$ is the radial velocity at the edge of the cloud, and $r_{\max}$ is the radius of the cloud. We assume a gas temperature of 20 000 K and $v_{\max} = 200$ km s⁻¹ (e.g. Zheng & Miralda-Escudé 2002). The neutral hydrogen column densities are set to be $N_H = 2 \times 10^{18}$, $2 \times 10^{19}$, and $2 \times 10^{20}$ cm⁻². They correspond to $\tau_0 = 8.3 \times 10^4$, $8.3 \times 10^5$, and $8.3 \times 10^6$, respectively. The source is located at the centre of cloud.

The emergent spectra are presented in Fig. 3.

The bulk expanding motion of the medium shifts the resonant frequency blueward in the laboratory frame (see Section 2.1). Hence, the blueward photons efficiently interact with hydrogen as they travel towards outer regions of the cloud. As a result, the emergent spectra acquire asymmetric shapes with an emission peak in the redward frequencies. The resultant emergent spectra are in good agreement with those from previous Monte Carlo codes (e.g. see fig. 2 of Semelin et al. 2007). Here, we do not present the infalling
<figure><img src="image_5.png" /><figcaption>Figure 4. The escape fraction from a dusty slab as a function of $(b\tau_0)^{1/3}\tau_{d,a}$. Red circles and blue crosses are, respectively, the results in the stochastic elimination of photon packets' method and in the 'attenuation of flux' method (see Section 2.3 for the detail). The dotted line represents the analytical solution obtained by Neufeld (1990).</figcaption></figure>

case, because the emergent spectrum becomes basically in a mirrored shape with respect to $x = 0$.

### 3.3 Static dusty slab

The escape fraction of Lyα photons from a dusty slab is a standard test, since it also can be solved analytically. Neufeld (1990) derived the analytical formula of the escape fraction $f_{esc}$ as

$$f_{esc} = \frac{1}{\cosh[\zeta'(a\tau_0)^{1/3}\tau_{d,a}]}, \qquad (22)$$

where $\zeta' = \sqrt{3}/(\zeta\pi^{5/12})$ with $\zeta = 0.525$, a fitting parameter. This formula is valid in the case of an optically thick medium with $(a\tau_0)^{1/3} \gg \tau_{d,a}$. For the test calculation, we assume spherical silicate dust grains with a radius $r_d$ of $2.4 \times 10^{-2}$ µm and solid material density of 3 g cm$^{-3}$ (Draine & Lee 1984) so that the dust opacity is equal to that in Verhamme et al. (2006). Note that the dust opacity under this assumption of single-sized grains is equivalent to that for the dust distributions by Mathis, Rumpl & Nordsieck (1977) given the geometrical cross-section in the range from $5.8 \times 10^{-4}$ to 1 µm (Yajima et al. 2017).

Fig. 4 shows the resultant escape fraction. We see in the figure that the numerically obtained escape fraction is concordant with the analytical solution, regardless of the adopted method for the dust absorption.

## 4 ADAPTIVITY OF SEURAT

In the mesh-free radiative transfer adaptive for SPH, we should often treat the propagation of Lyα photons in a local region having a very high optical depth or a large density gradient. Here, we describe the contrivance incorporated in SEURAT, which ensures the high adaptivity for media with high density contrast.

### 4.1 High optical depth regions

As explained in Section 2.2, photon packets are judged to escape from the system if the packets do not find a downstream particle. Therefore, an appropriate neighbour search is required. Otherwise, we wrongly judge the escape of a photon packet and it may impact

<figure><img src="image_6.png" /><figcaption>Figure 5. Conceptual figure for the constitution of the neighbour list. Red and open circles represent the ith particle and its neighbour particles, respectively. Blue circles are the neighbour particles of the jth particle, but not of the ith particle. Magenta depicts the particles belonging to the neighbour lists of both ith and jth particles. Dashed line denotes the smoothing kernel length of each SPH particle. We suppose these particles are in a high-density region and another particle (yellow circle) is located in a low-density region. The smoothing kernel length of the yellow particle is plotted by a solid line. A thick solid line represents the trajectory of the photon packet.</figcaption></figure>

the resultant emergent spectra and/or escape fraction. This point is critical for a mesh-free Monte Carlo based Lyα radiative transfer scheme, but has never been a concern for mesh-based schemes, since downstream grid cells are always recognized.

Fig. 5 shows a schematic view of the latent problems in the neighbour search procedure. In the mesh-free Monte Carlo technique, the photon packets fly searching a neighbour particle. As a result, a ray-tracing algorithm sometimes fails if one uses the 'gather' neighbour lists that are composed of the particles within the smoothing kernel length. This problem frequently arises, for instance, where the local optical depth is very large, since many scatterings occur around a single particle and the photon packet moves around the particle. When a scattering occurs at the boundary of the optically thick region, no downstream neighbour particle is found and then the packet is incorrectly labelled as an escaped packet (see the particle i in Fig. 5). To circumvent this problem, we assign a larger number of neighbour particles by making the searching radius larger than the smoothing length $h_i$. In addition, we extend the downstream particle search up to the 'neighbours of each neighbour particle', if no downstream particle can be found from the neighbour list. We checked that this recipe successfully avoids false escapes in a simple uniform optically thick system. As a result, we correctly solved the Lyα radiative transfer in the standard test calculations as already presented in the previous section (Figs 2, 3 and 4). Note that, from the viewpoint of the computational cost, we avoid the naive double loop processing to track the neighbour list of the neighbour particles, since the two particles usually share some particles in their neighbour lists (depicted by magenta circles in Fig. 5). Instead, we construct the neighbour lists that contains only the 'neighbours of the neighbours' (blue circles in Fig. 5).
## 4.2 Large density gradient regions

Moreover, incorrect escape events are highly probable when the local density gradient is quite large. Let us consider the particle *k* in Fig. 5 at the surface of a dense region surrounded by a low-density ambient gas and a photon packet propagating outward from the region. Since the ‘gather neighbour list’ of the particle *k* is mainly composed of particles in the dense region, we may not find any downstream particle of the particle *k* even if we trace the ‘neighbours of neighbours’. As a recipe coping with this situation, we construct the ‘scatter list’ for the particle *k*, i.e. the list of particles which includes the particle *k* in their gather neighbour list (Hernquist & Katz 1989), and use it instead of the gather neighbour list of the particle *k*. The gather lists of particles in the low density regions (particle *n* in Fig. 5) often contains the particle *k*. Hence, we successfully find downstream particles for the particle *k* from the scatter list. Therefore, we use the scatter list as the second fail-safe against false escapers.

## 4.3 Highly inhomogeneous media

To demonstrate the validity of the neighbour search scheme, we solve the Lyα radiative transfer in a highly inhomogeneous spherical cloud. In this test, we randomly distribute 16 clumps in a uniform spherical cloud. The total SPH number is 32³, and each clump consists of 32² particles. We assume three models for the density of the clumps; $n_{clump}/n_f = 10^2$, $10^3$, and $10^4$, where $n_{clump}$ and $n_f$ denote the gas density for a clump and for the field, respectively. The clump size is comparable to the average separation between field particles in the model of $n_{clump}/n_f = 10^3$. The source is located at the centre of the cloud and the number of photon packets is set to be $10^5$. Fig. 6 shows the positions at which each photon packet is judged to escape. As clearly seen in this figure, if we take only the ‘neighbour of neighbour’ (green dashed line), photon packets escape before they reach the edge of the cloud. This is because, owing to the high density contrast, the neighbour lists of the SPH particles in the clumps hardly involve any SPH particles residing in the low-density field. As a result, photon packets in the high-density regions cannot find a downstream particle and therefore mistakenly label their positions as the boundary of the system. Thus, the ‘neighbours of neighbours’ method is obviously insufficient to avoid incorrect escape of the photon packets. On the other hand, if we take the scatter

list into account in the neighbour search procedure, we correctly solve the travel of photon packets (red solid line). As shown in the panels, all photon packets correctly recognize the edge of the system, regardless of the level of the density contrast. We emphasize that the scheme is valid even when the clump size is smaller than the typical separation of the field particles. Hence, we conclude that the neighbour search scheme in conjunction with the scatter list is indispensable to solve the mesh-free Monte Carlo based radiative transfer.

## 4.4 Comparison between SEURAT and LICORICE

In this section, we compare the transfer in the Lyα lines with SEURAT to that with LICORICE, which is a mesh-based RT code (Semelin et al. 2007).

### 4.4.1 Model galaxy

For this comparison, we use a model galaxy obtained by a cosmological radiation SPH simulation (see Hasegawa & Semelin 2013, for the details of the simulation). The model galaxy has a halo mass of ~5.7 × 10¹⁰ M☉ at z ~ 6.1 and is composed of ~4.0 × 10⁴ SPH particles. The minimum smoothing length is ~0.15 kpc in the simulation. In LICORICE, we generate a uniform grid with 512³ cells of size ~0.20 kpc to resolve the minimum smoothing length of the SPH simulation. The radiative transfer of ionizing photons is solved coupled with hydrodynamics, assuming the Case B recombination. The resultant distributions of neutral hydrogen and temperature are shown in Fig. 7. Here, we do not incorporate dust extinction. As for the Lyα emissivity, we consider the recombination and the collisional excitation of hydrogen atoms. Then, the local Lyα emissivity is, respectively, given by

$$ \epsilon_{\alpha}^{\text{rec}} = f_{\alpha} \alpha_{\text{B}} h v_{\alpha} n_{\text{HII}} n_{\text{e}}, \qquad (23) $$

$$ \epsilon_{\alpha}^{\text{col}} = C_{\text{Ly}\alpha} n_{\text{H}} n_{\text{e}}, \qquad (24) $$

where $f_α$ is the average number of the Lyα photons generated via the Case B recombination, ($f_α = 0.68$, Osterbrock & Ferland 2006), $α_B$ is the Case B recombination coefficient (Hui & Gnedin 1997), $C_{Lyα}$ is the collisional excitation

<figure><img src="image_7.png" /><figcaption>Figure 6. Distribution of the positions at which photon packets escape. From the left to right-hand panel, the panels show the results for the density contrast of $n_{clump}/n_f = 10^2$, $10^3$, or $10^4$, respectively. The vertical axis denotes the number count of the escape positions, and the horizontal axis is the distance from the centre of the cloud normalized by the cloud radius $r_{cloud}$. In each panel, a red solid line represents the results employing the scatter list in addition to the ‘neighbour of the neighbour’ list, whereas a green dashed line is the results dispensing with the scatter list.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>**Figure 7.** Distributions of neutral hydrogen number density (left-hand panel) and temperature (right-hand panel) for a model galaxy.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 8.** Emergent spectra from a model galaxy. The horizontal axis is the wavelength and the vertical axis denotes the number count of photon packets in each wavelength bin. A solid curve represents the result with *SEURAT*, whereas a dashed curve does the result with *LICORICE*. The spectra are normalized such that the integration of the spectrum over frequencies should be unity. In upper panel, the velocity structure is disregarded, whereas in the lower panel the velocity structure is included.</figcaption></figure>

coefficient, $C_{Ly\alpha} = 3.7 \times 10^{-17} \exp(-hv\alpha/k_B T) T^{-1/2} \text{ erg s}^{-1} \text{ cm}^3$ (Osterbrock & Ferland 2006). We cast $10^5$ photon packets.

### 4.4.2 Emergent spectrum

We calculate the emergent $Ly\alpha$ spectra for the model galaxy with *SEURAT* and *LICORICE*. In Fig. 8, the spectra calculated with the two methods are compared. In the upper panel, we set all velocities to

<figure><img src="image_10.png" /><figcaption>**Figure 9.** Velocity profile in a certain direction. Upper panel shows the profile from the centre to the edge of the galaxy. Lower panel represents the profile from a source located at ~36 kpc from the centre. A solid line indicates the velocities in *SEURAT*, which are determined according to the neighbour search procedure, whereas a dashed line is the velocity in *LICORICE*.</figcaption></figure>

zero to evaluate the difference between the two methods in a static problem. As seen in this figure, the emergent spectra agree well with each other. This demonstrates that there is no significant difference between *SEURAT* and *LICORICE* for a static system. The lower panel in Fig. 8 represents the result including the velocity structure. In both cases, we obtain asymmetric spectra shifted to shorter wavelength, which is typical of infalling gas. However, the red peak is slightly weaker in the result with *SEURAT*, and the blue peak is stronger, compared to that with *LICORICE*. Since the agreement is excellent in a static problem, we believe this discrepancy comes from the difference in tracing the velocity structure between the two methods. To verify this, we check the velocity profiles in the two methods. In Fig. 9, we show the profiles measured from different points. The upper panel represents the velocity profile from the centre to the edge of the galaxy. We find that the two profiles behave in a fairly similar way. The lower panel represents the profile from a certain source located at ~36 kpc from the centre. In this figure, we see a significant discrepancy at $0 \lesssim r \lesssim 5$ kpc, where the velocity is negative in *SEURAT*, but is positive in *LICORICE*. This may come from the difference in the way the velocity is computed. In *SEURAT*, a local SPH particle with $v < 0$ is possibly selected when assigning the velocity along the line of sight, whereas it can be smoothed out by the interpolation on the grid in *LICORICE*. Such a discrepancy is anticipated to become apparent especially in low-density regions. Actually, we have confirmed that if we dismiss the sources in
<figure><img src="image_11.png" /><figcaption>**Figure 10.** Cumulative optical depth at the position X (upper panels) and the ratio between *SEURAT* and the *KB* scheme (lower panels). The horizontal axis is the position X normalized by $r_{gal}$, where $r_{gal}$ is the radius of the galaxy (twice the size of the virial radius of the halo). In the upper panels, a solid line represents the result with *SEURAT*, whereas a dashed line is that with the *KB* scheme. The ray-tracing is performed along with the X-axis. Left-hand panels represent the result in the negative direction, whereas right-hand panels display that in the positive direction.</figcaption></figure>

low-density regions of $r > 4.5$ kpc, then the emergent spectra fall in a good agreement. Thus, we speculate that the difference of the velocity assignment in low density regions causes the discrepancy in the spectra.

### 4.4.3 Validity of the ray-tracing scheme in SEURAT

As described in Section 2.2, *SEURAT* directly uses the density at *ith* SPH-particle, following the method adopted in Radiation-SPH (RSPH) (Susa 2006). However, the densities can be evaluated at grid-points along a ray by using the kernels of SPH-particles neighbouring the grid points, as originally implemented by Kessel-Deynet & Burkert (2000) (hereafter *KB* scheme). To see the difference between the two methods for the density assignment, we apply the two methods for the model galaxy and compare the optical depth along a ray. To clarify the difference of the density assignment, we set field velocities to be zero, and neglect the scattering process in this test. We calculate the optical depth at the line centre frequency from the galactic centre to the virial radius along the X axis. Fig. 10 shows the resultant cumulative optical depth $\tau(X)$ as a function of position along the ray. Due to the high neutral hydrogen density, the optical depth steeply rises up to $\tau \sim 10^{10}$ in the compact inner regions of the system ($|X|/r_{gal} \lesssim 10^{-2}$). It is notable that the two optical depths are in a good agreement in the inner regions. However, the agreement of $\tau(X)$ breaks down at $|X|/r_{gal} \gtrsim 10^{-2}$. Actually, the discrepancy appears at a position where the density gradient is large. In *SEURAT*, the neighbour list of particles in sparse regions does not necessarily consist of particles in sparse regions. In other words, some particles can have neighbour particles both in the dense and sparse regions (see the particle *n* in Fig. 5). On the other hand, it is possible that particles only in the sparse regions contribute to the density assignment in the *KB* scheme (the particle *n* and *n'* in Fig. 5). Thus, we see in Fig. 10 that the cumulative optical depth calculated with *SEURAT* tends to be larger than that with the *KB* scheme. However, including the information of dense particles would not be always a shortcoming. Increasing the number of photon packets, the possibility that a

photon packet chooses the trajectory passing through a dense region (thick dashed line in Fig. 5) is enhanced. The ray-tracing scheme in *SEURAT* may effectively chose the trajectories through dense regions even for the relatively small number of photon packets. Also, it is worth mentioning that although the *KB* scheme seems to fit the concept of SPH, it would not be always accurate, since a sufficient number of the superposing particle is not guaranteed in the sparse regions. Anyway, the difference is between $1 \times 10^{10}$ and $1.5 \times 10^{10}$. A few 10 per cent difference in such a large optical depth would not lead to a significant dissimilarity in the line profile.

## 5 CONCLUSIONS

We have developed a new numerical code, *SEURAT*, to solve the Ly$\alpha$ radiative transfer adaptively based on the SPH particle distribution, using a Monte Carlo method. Although the previous Monte Carlo codes have been developed as mesh-based schemes, *SEURAT* is a mesh-free scheme, which directly uses the SPH particles themselves for the integration of optical depth. Hence, *SEURAT* can solve the Ly$\alpha$ radiative transfer without reducing the resolution of the SPH simulations. We have performed the standard test calculations, which include the emergent spectra from a static homogeneous spherical cloud, the emergent spectra from an expanding homogeneous spherical cloud, and the escape fraction from a dusty slab. We have shown that the results reproduce the analytic solutions or the results obtained by previous studies. To solve Ly$\alpha$ radiative transfer in highly inhomogeneous media, special attention has been paid for the neighbour list construction. The neighbour list of the neighbour particle allows us to properly treat local optically thick regions. Furthermore, a scatter list is also required to perform the ray-tracing in systems having large density gradients. We have demonstrated that *SEURAT* reliably searches the neighbour particles, and correctly performs the ray-tracing in significantly clumpy media. As a result, we have confirmed that *SEURAT* is successfully designed to manage extremely optically thick and highly inhomogeneous media. Finally, we have compared the *SEURAT* to the mesh-based code *LICORICE*. We have applied the codes to a model galaxy and calculated the emergent spectra. We have shown that for a static system the two methods produce very similar spectra. When including velocity gradients, the velocity assignment in low-density regions may lead to some moderate differences in the emergent spectra.

So far, the SPH method has been widely employed to study galaxy formation. One of the potential applications of *SEURAT* is to couple it with simulations of galaxy formation to model the high-$z$ LAEs. In modelling LAEs, both resonant scattering with the hydrogen atoms and dust scattering/absorption are essential physical processes. A primary advantage of *SEURAT* is to treat such processes adaptively with the resolution of SPH. We plan to carry out simulations of Ly$\alpha$ radiative transfer coupled with numerical simulations of galaxy formation. The results will be presented in a forthcoming paper.

## ACKNOWLEDGEMENTS

We are grateful to A. Inoue for fruitful discussions. The numerical simulations have been performed with COMA provided by Interdisciplinary Computational Science Program in Center for Computational Sciences, University of Tsukuba, with the K computer provided by the RIKEN Advanced Institute for Computational Science and with Cray XC30 at Center for Computational Astrophysics, NAOJ. This research was supported in part by Grant-in-Aid for Scientific Research (B) No.15H03638 (MU), Grant-in-Aid for Scientific Research (A) No.17H01110 (KH), and Grant-in-Aid for
Young Scientists (A) No.17H04827 (HY) through Japan Society for the Promotion of Science.

## REFERENCES

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

Altay G., Croft R. A. C., Pelupessy I., 2008, MNRAS, 386, 1931

Atek H., Kunth D., Hayes M., Östlin G., Mas-Hesse J. M., 2008, A&A, 488, 491

Baek S., Di Matteo P., Semelin B., Combes F., Revaz Y., 2009, A&A, 495, 389

Dijkstra M., Haiman Z., Spaans M., 2006a, ApJ, 649, 14

Dijkstra M., Haiman Z., Spaans M., 2006b, ApJ, 649, 37

Dijkstra M., Lidz A., Wyithe J. S. B., 2007, MNRAS, 377, 1175

Dijkstra M., Mesinger A., Wyithe J. S. B., 2011, MNRAS, 414, 2139

Draine B. T., Lee H. M., 1984, ApJ, 285, 89

Faucher-Giguère C.-A., Kereš D., Dijkstra M., Hernquist L., Zaldarriaga M., 2010, ApJ, 725, 633

Finkelstein S. L. et al., 2013, Nature, 502, 524

Genzel R. et al., 2011, ApJ, 733, 101

Gronke M., Bull P., Dijkstra M., 2015, ApJ, 812, 123

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

Hasegawa K., Semelin B., 2013, MNRAS, 428, 154

Hasegawa K., Umemura M., 2010, MNRAS, 407, 2632

Hernquist L., Katz N., 1989, ApJS, 70, 419

Hopkins P. F., Kereš D., Oñorbe J., Faucher-Giguère C.-A., Quataert E., Murray N., Bullock J. S., 2014, MNRAS, 445, 581

Hui L., Gnedin N. Y., 1997, MNRAS, 292, 27

Iye M. et al., 2006, Nature, 443, 186

Kashikawa N. et al., 2006, ApJ, 648, 7

Kessel-Deynet O., Burkert A., 2000, MNRAS, 315, 713 (KB)

Konno A. et al., 2014, ApJ, 797, 16

Laursen P., Razoumov A. O., Sommer-Larsen J., 2009, ApJ, 702, 824

Mathis J. S., Rumpl W., Nordsieck K. H., 1977, ApJ, 217, 425

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

Ono Y. et al., 2012, ApJ, 744, 83

Osterbrock D. E., Ferland G. J., eds, 2006, Astrophysics of Gaseous Nebulae and Active Galactic Nuclei. University Science Books, Mill Valley, CA

Ouchi M. et al., 2008, ApJS, 176, 301

Ouchi M. et al., 2010, ApJ, 723, 869

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

Pawlik A. H., Schaye J., 2008, MNRAS, 389, 651

Santos M. R., 2004, MNRAS, 349, 1137

Schaye J. et al., 2015, MNRAS, 446, 521

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

Shibuya T., Kashikawa N., Ota K., Iye M., Ouchi M., Furusawa H., Shimasaku K., Hattori T., 2012, ApJ, 752, 114

Smith A., Safranek-Shrader C., Bromm V., Milosavljević M., 2015, MNRAS, 449, 4336

Springel V., 2010, ARA&A, 48, 391

Susa H., 2006, PASJ, 58, 445

Tasitsiomi A., 2006, ApJ, 645, 792

Vanzella E. et al., 2011, ApJ, 730, L35

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

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

Vogelsberger M. et al., 2014, Nature, 509, 177

Wise J. H., Turk M. J., Norman M. L., Abel T., 2012, ApJ, 745, 50

Yajima H., Li Y., 2014, MNRAS, 445, 3674

Yajima H., Umemura M., Mori M., Nakamoto T., 2009, MNRAS, 398, 715

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

Yajima H., Li Y., Zhu Q., Abel T., Gronwall C., Ciardullo R., 2012b, ApJ, 754, 118

Yajima H., Li Y., Zhu Q., Abel T., Gronwall C., Ciardullo R., 2014, MNRAS, 440, 776

Yajima H., Li Y., Zhu Q., Abel T., 2015, ApJ, 801, 52

Yajima H., Ricotti M., Park K., Sugimura K., 2017, ApJ, 846, 3

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

Zitrin A. et al., 2015, ApJ, 810, L12

This paper has been typeset from a TeX/LA TeX file prepared by the author.