{'fraction_non_alphanumeric': 0.03467646180659228, 'fraction_numerical': 0.019726837803870377, 'mean_word_length': 4.414414075717132, 'pattern_counts': {'":': 0, '<': 39, '<?xml version=': 0, '>': 7, 'https://': 1, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 19, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Ly-alpha Radiative Transfer: A Stokes Vector Approach to Ly-alpha Polarization

Ly-alpha emitting galaxies and giant Ly-alpha blobs (LABs) have been extensively observed to study the formation history of galaxies. However, the origin of their extended Ly-alpha emission, especially of LABs, remains controversial. Polarization signals from some LABs have been discovered, and this is commonly interpreted as strong evidence supporting that the extended Ly-alpha emission originates from the resonance scattering. The Monte Carlo Ly-alpha radiative transfer code LaRT is updated to investigate the polarization of Ly-alpha using the Stokes vector formalism. We apply LaRT to a few models to explore the fundamental polarization properties of Ly-alpha. Interestingly, individual Ly-alpha photon packets are found to be almost completely polarized by a sufficient number of scatterings (N_scatt>10^4-10^5 in a static medium) or Doppler shifts induced by gas motion, even starting from unpolarized light. It is also found that the polarization pattern can exhibit a non-monotonically increasing pattern in some cases, besides the commonly-known trend that the polarization monotonically increases with radius. The polarization properties are primarily determined by the degree of polarization of individual photon packets and the anisotropy of the Ly-alpha radiation field, which are eventually controlled by the medium's optical depth and velocity field. If once Ly-alpha photon packets achieve ~100% polarization, the radial profile of polarization appears to correlate with the surface brightness profile. A steep surface brightness profile tends to yield a rapid increase of the linear polarization near the Ly-alpha source location. In contrast, a shallow surface brightness profile gives rise to a slowly increasing polarization pattern.


INTRODUCTION
Lyα is one of the most powerful tracers to investigate star-forming galaxies and circumgalactic/intergalactic media (CGM/IGM) in the universe. After the prediction by Partridge & Peebles (1967) that star-forming galaxies should be strong Lyα emitters (LAEs), many photometric and spectroscopic surveys have been performed to detect Lyα emission from nearby or high-redshift galaxies (Meier & Terlevich 1981;Cowie & Hu 1998;Rhoads et al. 2000;Östlin et al. 2014;Bacon et al. 2015;Ouchi et al. 2018Ouchi et al. , 2020. The observed Lyα spectra and/or surface brightness profiles of LAEs can be well understood by resonance scattering of Lyα in simplified galactic outflow models. The most well-known one is a thin shell model, in which a central Lyα source is surrounded by a constantly expanding spherical thin shell of atomic hydrogen gas. This shell model has surprisingly well reproduced diverse Lyα line profiles (Ahn 2004;Verhamme et al. 2008;Schaerer et al. 2011;Gronke et al. 2015;Yang et al. 2016;Gronke 2017;Karman et al. 2017), even though it is questionable whether the model can reproduce the observed surface brightness profiles as well. Another model designed to better represent real galaxy halos was able to successfully reproduce both the spectra and surface brightness profiles of LAEs (Song et al. 2020). In this halo model, the distributions of Lyα source and gas are described by exponential functions of radius. The outflowing velocity of the halo is characterized by a piecewise linear function of radius.
The extended Lyα nebulae (also known as Lyα blobs or LABs) observed at z = 2 − 6 can provide clues to galaxy formation in the early universe (Keel et al. 1999;Steidel et al. 2000;Matsuda et al. 2004;Yang et al. 2009Yang et al. , 2010Prescott et al. 2012;Bȃdescu et al. 2017). They are often found in association with high-density regions of LAEs; the connection of LABs with the overdensities of LAEs suggests that they are associated with matter density peaks in the universe and thus likely to evolve into the present-day groups and clusters of galaxies. However, the mechanism powering the extended Lyα emission in LABs remains debated. Many possible mechanisms have been proposed to explain LABs, including (1) rapid cooling of the accreting gas that is heated by galactic outflows during powerful starbursts (Taniguchi & Shioya 2000;Mori et al. 2004;Geach et al. 2005) or by the dissipation of gavitational energy as gas falls toward galaxies (Haiman et al. 2000;Fardal et al. 2001;Dijkstra & Loeb 2009), (2) photoionization by luminous active galactic nuclei (AGNs), young stars, and/or the intergalactic ultraviolet background (Haiman & Rees 2001;Jimenez & Haiman 2006;Gould & Weinberg 1996), and (3) resonant scattering of Lyα photons produced by star forming galaxies and/or AGNs hosted within the nebulae Steidel et al. 2011;Beck et al. 2016). In the former two scenarios, Lyα photons are produced in situ inside the nebulae, whereas they are scattered light by neutral hydrogen gas surrounding the central Lyα source(s) in the last scenario.
It is challenging to discern which scenario among the above is better suited using only the surface brightness profile and spectrum. With this regards, polarization signal can provide critical, additional information about the nature of the Lyα nebulae or the diffuse Lyα halos around high-redshift galaxies, which cannot be offered by the photometric and spectroscopic measurements alone; the resonant scattering of Lyα gives rise to polarization, but other mechanisms would not. Therefore, the polarization phenomenon of Lyα has attracted attention in both observational and theoretical studies to reveal the nature of LABs. The first attempt to detect any polarization signal of the Lyα emission from LABd05 was not sensitive enough (Prescott et al. 2011). Later, subsequent observations have detected clear polarization signatures in Lyα nebulae Humphrey et al. 2013;Beck et al. 2016;You et al. 2017;Kim et al. 2020), which qualitatively accord with the theoretical prediction that the polarization level of the scattered Lyα tends to increase with distance from the central source. The detection of polarization signals that follow the predicted trend is interpreted as strong evidence supporting that the LABs are caused primarily by the resonance scattering of Lyα originating from star-forming galaxies or AGNs. Moreover, Li et al. (2021) showed that the observed Lyα/Hβ ratios in LAB2 support the resonant scattering scenario, and the observed Lyα spectra can be well reproduced using an RT model in a clumpy medium. In contrast, Trebitsch et al. (2016) demonstrated that the radial profile of polarization could be well explained by the scenario in which Lyα photons are produced in the cooling gas surrounding galaxies and then self-scattered by the gas. Their results likely suggest that the polarimetric observations alone are unable to disentangle between the two scenarios. Recently, Herenz et al. (2020) found that ionizing radiation and mechanical heating would be dominant near the embedded galaxies in LAB1, based on the observation of He II λ1640 and non-detection of C IV λλ1548, 1550, while Lyα scattering contributes more significantly at larger distances.
To interpret the above observations, modeling of Lyα radiative transfer (RT) process is essential. This is because Lyα is a strong resonance line and thus undergoes a large number of scatterings by neutral hydrogen atoms in the interstellar medium (ISM) within galaxies, in the CGM that immediately surrounds them, and in the IGM between galaxies. Monte Carlo simulation is the most flexible and straightforward approach in calculating the Lyα RT effects in astrophysical systems. The Monte Carlo simulation technique for Lyα RT has now become almost a standard, and thus many Monte Carlo Lyα RT codes have been developed in many different contexts (Avery & House 1968;Ahn et al. 2000;Zheng & Miralda-Escudé 2002;Cantalupo et al. 2005;Tasitsiomi 2006;Dijkstra et al. 2006;Verhamme et al. 2006;Semelin et al. 2007;Laursen & Sommer-Larsen 2007;Laursen et al. 2009;Pierleoni et al. 2009; Barnes & Haehnelt 2010;Forero-Romero et al. 2011;Yajima et al. 2012;Orsi et al. 2012;Behrens & Niemeyer 2013;Gronke & Dijkstra 2014;Smith et al. 2015;Michel-Dansac et al. 2020;Seon & Kim 2020). Most of them have focused on predicting the emergent Lyα spectrum and surface brightness, except Semelin et al. (2007) and Seon & Kim (2020), who investigated the Wouthuysen-Field (WF) effect of Lyα radiation on 21 cm emission.
In addition to predicting the spectra and surface brightness profiles, RT simulations for the Lyα polarization have also been developed (Ahn et al. 2002;Dijkstra & Loeb 2008;Dijkstra & Kramer 2012;Ahn & Lee 2015;Trebitsch et al. 2016;Chang et al. 2017;Eide et al. 2018;Mas-Ribas & Chang 2020). Lee & Ahn (1998) proposed, for the first time, that the polarimetry of Lyα emission can be a useful probe about the geometry and kinematics of the neutral hydrogen gas in nearby starburst galaxies and primeval galaxies. It was also shown that the scattered Lyα radiation would gain a high level of polarization for a wide range of galactic outflows in high-redshift galaxies (Dijkstra & Loeb 2008). Prior to the epoch of cosmic reionization, the Lyα photons emitted by early galaxies would be resonantly scattered by the neutral IGM in their vicinity, leading to a fairly compact Lyα halo . The scattered Lyα light would be highly polarized, and the polarized, intergalactic Lyα halos may then provide a unique tool for probing the neutral IGM before and during the epoch of reionization . These calculations showed that the degree of polarization increases, even up to ∼ 40 − 60%, with increasing radial distance from the central source.
To date, two methods have been developed to implement polarization in the Lyα RT calculations. The first approach deals with individual "monochromatic" photons that are 100% linearly polarized Dijkstra & Loeb 2008;Trebitsch et al. 2016). The second approach employs the density matrix to represent the polarization state of a partially-polarized "quasi-monochromatic" ensemble of photons (Lee et al. 1994;Lee & Ahn 1998;Ahn et al. 2002;Ahn & Lee 2015;Chang et al. 2017). However, neither approaches are complete, at least from a theoretical point of view. First, instead of adopting a scattering phase function that constantly varies as the Lyα frequency changes, they use two distinct phase functions selected according to the frequency regimes. Second, the change in polarization state due to dust scattering was not taken into account properly. Third, no circular polarization, which can naturally arise by dust grains in some circumstances, was included in the formulae.
The present study was mainly motivated in two respects: (1) incomplete treatment of Lyα polarization in previous theoretical studies, and (2) the controversy over the interpretation of the polarization signals observed in LABs. For this purpose, we developed a new, third approach using the Stokes parameters. The method is superb in dealing with quantum-mechanically-derived phase function compared to the methods mentioned above. It is also able to take the dust effect into account appropriately when dealing with the Lyα polarization. In our first paper (Seon & Kim 2020), we introduced the state-of-the-art Monte Carlo Lyα RT code LaRT 6 to study Lyα-related phenomena, including the WF effect. LaRT was utilized to simultaneously analyze the Lyα spectra and surface brightness profiles of the LAEs at z = 3 − 6 (Song et al. 2020). This paper, the third in this series, extends the code to investigate the Lyα polarization. We also present the results from various simple models to demonstrate the ability of LaRT and discuss intriguing properties on the Lyα polarization obtained from the models. This paper is organized as follows. In Section 2, we present algorithms developed for LaRT, which were not explained in the first paper. Section 3 explains two main factors, the polarization of individual photon packets and the degree of isotropy (or anisotropy) in the radiation field, that determine the ensemble average of polarization signals. These two concepts are essential to understanding the results of this paper. Section 4 describes the results obtained from four types of models. We discuss a few topics that are relevant to our results and Lyα RT in Section 5. In Section 6, we summarize the main results. In Appendices, we derive the scattering matrix for the Lyα scattering by hydrogen atoms (A), and provide approximate formulae for the scattering matrix elements for the scattering by dust (B). We also discuss the algorithm to randomly draw scattering angles from the phase function using the inversion method (C). In addition, we obtain generalized formulae for the density matrix method from our formulae (D) and discuss the technique utilizing 100% polarized photons (E).

MONTE CARLO RADIATIVE TRANSFER METHODS
The basic Lyα RT algorithms of LaRT, except those for polarization and the "peeling-off" technique, were described in the first paper of this series (Seon & Kim 2020). In the following, we describe only the basic definitions and the contents relevant to polarization and peeling-off. We use the polarization RT algorithm that is similar to those of Bianchi et al. (1996) and Peest et al. (2017), which were developed for the dust RT. We assume that a newly emitted photon packet is unpolarized, unless otherwise stated. The dust grains are assumed to be spherical and thus no dichroic attenuation due to elongated dust grains is considered (see Seon 2018 for the polarization caused by the dichroic attenuation). The polarization algorithm described in this paper was already exploited to study the polarization of continuum radiation in the dusty halo of edge-on galaxies (Seon 2018); but, no detailed description of the algorithm was published.
The overall sequence of the RT algorithm is as follows. First, we setup the temperature, density, and bulk velocity fields in a three dimensional grid system. Second, an initial position, a propagation direction, and a frequency of the individual photon packet are randomly generated according to appropriate probability distribution functions. Third, an optical depth τ that the photon will propagate before it interacts with a hydrogen atom or a dust grain is randomly drawn from the exponential distribution function P(τ ) = e −τ . The physical distance is then calculated by numerically inverting the optical depth integral given as a function of distance. Fourth, we choose a scattering agent, either a hydrogen atom or dust grain, from the relative ratio between the optical depths due to hydrogen atoms and dust grains. We then find the velocity vector of an atom that scatters the photon if it is scattered by a hydrogen atom. If the photon is assumed to interact with dust, it is either absorbed or scattered by a dust grain; the decision is made by comparing a uniform random number with the dust albedo a. The absorbed photon can be totally removed from the system or its weight is reduced by a factor of 1 − a, depending on the simulation mode. Unless the photon is removed from the system by a dust grain, we proceed further RT processes for the photon packet. Fifth, we draw the scattering angles (θ, φ) of the photon from an appropriate scattering phase function, and update the Stokes vector and propagation direction (Sections 2.5-2.8). A new frequency is also obtained, if the scatterer is a hydrogen atom, from the velocity of the scattering atom and the old and new direction vectors of the photon. If a dust grain is the scatterer, the frequency is not altered. Sixth, to predict observational spectropolarimetric data, the peeling-off procedure is performed for a predefined virtual detector plane at the observer's location (Section 2.9). The processes are repeated until the photon either escapes the system or is removed by dust. The whole procedures are performed for all photon packets and the resulting outputs are converted into physical units.

Scattering Cross Section
The fine structure of the n = 2 quantum state of hydrogen should be considered to take the Lyα polarization into account properly. We refer the transition 2 S 1/2 ↔ 2 P o 1/2 (corresponding to the lower frequency) to as "H" and the transition 2 S 1/2 ↔ 2 P o 3/2 (higher frequency) to as "K," as in the Ca II λλ3933, 3968 doublet line. The frequencies for the H and K transitions are denoted by ν H and ν K , respectively. The central frequency of Lyα is ν α = 2.466 × 10 15 Hz and the frequency difference between the fine structure levels is ∆ν HK = ν K − ν H = 1.08 × 10 10 Hz, which is equivalent to a Doppler shift of 1.34 km s −1 . The scattering cross-section of a Lyα photon in the rest frame of a hydrogen atom is derived in Appendix A: where χ 0 = f α πe 2 /m e c. Here, f α = 0.4162 is the oscillator strength, and Γ = A α = 6.265 × 10 8 s −1 the damping constant (the Einstein A coefficient) of the Lyα transition. The multiplication factor 2 for the K transition is due to the difference in the statistical weights of 2J + 1 between the two transitions. Integrating over the one-dimensional Maxwellian velocity distribution of the hydrogen gas at temperature T , the cross section in a reference frame comoving with the gas fluid results in where H(x, a) is the Voigt-Hjerting function given by We define x H , x K and x as the relative frequencies of the photon normalized to the thermal Doppler width ∆ν D = ν α (v th /c): Here, v th = (2k B T /m p ) 1/2 and a = Γ/(4π∆ν D ) are the thermal speed of hydrogen gas and the natural width parameter of H(x, a), respectively. Note that ∆ν HK /∆ν D (T /100 K) −1/2 ; thus, the fine structure splitting is negligible if T 10 2 K. The optical depth τ ν (s) of a photon with frequency ν traveling along a path length s is given by where n(v ) is the number density of the neutral hydrogen atom with the velocity component v parallel to the photon's propagation direction. In this paper, the total amount of gas is measured using the column density N HI of neutral hydrogen or the optical depth τ 0 defined by where φ x = H(x, a)/ √ π is the normalized Voigt profile. This definition of τ 0 is consistent with that of the monochromatic optical depth measured at the line center (x = 0) when the fine structure is ignored. The parallel component u = v /v th of the scattering hydrogen atom is obtained from the following composite distribution function: where If a uniform random number ξ (0 ≤ ξ ≤ 1) is choosen to be smaller than P H , the photon is scattered via the H transition, otherwise, via the K transition. A random parallel velocity component of the scattering atom, given a transition, is obtained by the algorithm developed in Seon & Kim (2020). When considering the fine-structure splitting, we need to compute H(x, a) more than twice as much as the case ignoring the splitting. It is, therefore, essential to use a fast algorithm to evaluate H(x, a). We developed a fast yet accurate function, which is inspired by the approximation formula of Harris (1948) and the function "voigt_king" in the VPFIT program (Carswell & Webb 2014). It was found that the routine is accurate over the domain of a 0.4 (T 1 K for the hydrogen gas) and the whole range of x, with a relative error lower than 10 −4 . This routine was already utilized in Seon & Kim (2020). We also implemented two more Fortran subroutines translated from the Matlab routine of Abrarov & Quine (2015), which are very accurate with an average accuracy of 10 −14 over a broad domain of (a, x).

Definition of the Stokes Vector
We describe the polarization state of a photon packet using the Stokes vector where I represents the intensity of the Lyα radiation field, Q and U describe the linear polarization, and V the circular polarization. In the present formulation, we employ a local reference frame of the photon packet, defined by the propagation direction k and two orthogonal unit vectors m and n in the plane perpendicular to k, as shown in Figure  1. The orthogonal triad of unit vectors (m, n, k) constitutes a right-handed coordinate system. The reference frame is attached to the photon packet and thus comoves and corotates with the photon as it propagates and is scattered off. The Stokes vector is defined relative to the two vectors m and n in the photon packet's local reference frame, so m and n are referred to as the polarization basis vectors. The Stokes parameters are defined by two orthogonal complex electric field components along the m and n directions, as follows: where · · · denotes an ensemble average for the quasi-monochromatic wave. The sign of Stokes V should be appropriately chosen according to the definition of the electric field phase. In our convention, the phase is assumed to be e −iωt for the angular frequency ω = 2πν; V should be multiplied by −1 if the phase of e iωt is adopted. The linear polarization angle χ is defined as the angle between the direction of linear polarization and the basis vector m. Using χ, the Stokes parameters Q and U can be represented as Q = Q 2 +U 2 cos 2χ, The degree of linear polarization is given by There have been, confusingly, many different conventions concerning the polarization angle direction of linear polarization and the handedness of circular polarization, i.e., the signs of the Stokes parameters Q, U, and V . This paper uses the "right-handed" system following the International Astronomical Union recommendation (see Hamaker & Bregman 1996). In the IAU standard, a polarization angle is measured counterclockwise from the North to the East in the sky when looking at the photon source. In LaRT, the polarization angle is defined to be an angle starting from m, rotated about the axis k. The vectors m and n do not always coincide with the North and East in the sky. However, when recording the output signals on a detector plane in the peeling-off stage, we rotate m and n to coincide with the North and East directions in the sky, respectively.
Peest et al. (2017) Figure 2. The E 1 parameter for the scattering matrix of Lyα as a function of wavelength (velocity). The wavelengths for the K and H transitions of Lyα are λ K and λ H , respectively. λ min is the wavelength at which E 1 has the minimum value. E 1 has a negative value in the wavelength range of λc < λ < λ H ; a negative E 1 implies that the polarization pattern resulting from a central source is radial rather than concentric. Interestingly, λc and λ min divide the interval λ H − λ K equally into three parts. V adopted by various authors. 7 For example, Chandrasekhar (1960) measures the polarization angle clockwise when looking at the source (left-handed basis; U Chandra = −U IAU ), while Stokes V is defined in the same way as the IAU standard (V Chandra = V IAU ). Ahn et al. (2002) and Eide et al. (2018) use the same convention as Chandrasekhar (1960) in the definition of the density matrix. This paper and Peest et al. (2017) represent the Stokes vector in the continuously moving, local reference frame of the photon, while Ahn et al. (2002) and Eide et al. (2018) use a fixed (or laboratory) reference frame to define the density matrix. In Appendix D, a more general density matrix formula than that given in Ahn et al. (2002) is obtained from the Stokes vector scheme of the present study; the appendix verifies the equivalence of the two approaches.

Transform of the Stokes and Basis Vectors
After a scattering event, the polarization state is updated by multiplying the initial Stokes vector by a scattering matrix (also known as the Müller matrix), suitable for the type of scattering. The scattering matrix elements depend on the geometrical configuration and the material property; they also depend on the photon frequency. The scattering matrix M(θ) is usually defined so that the basis vector m lies in the scattering plane, constructed by the incoming and outgoing propagation directions k and k of the photon, as shown in Figure 1. The scattering matrix is thus a function of only the scattering (polar) angle θ between k and k (cos θ = k · k ).
Suppose a photon is scattered in a direction with a polar angle θ and azimuth angle φ in a coordinate system defined by (m, n, k). In that case, we first need to rotate the polarization basis vectors (m, n) by the azimuth angle φ to make m lie 7 Peest et al. (2017) adopt the right-handed basis as in this paper. However, their approach differs from ours in that they measure the polarization angle starting from n, regarding n and −m as being the North and East direction. In LaRT, a photon packet carries all three basis vectors (m, n, k). On the other hand, Peest et al. (2017) use only two of them (k and n) in the hope of reducing the number of computations. However, we found that the number of calculations is the same if the computation sequence is well organized, even when using all three basis vectors. It is also more convenient to employ the three basis vectors.
in the scattering plane before applying the scattering matrix M. The Stokes vector should then be expressed in the rotated, reference frame by multiplying it by a rotation matrix L(φ): The final Stokes vector after the scattering event is obtained by mutiplying M(θ) to L(φ)S, as follows: After every scattering event, we also need to update the basis vectors (m, n). The new basis vectors are obtained first by rotating them by φ about k, updating (m, n), and then subsequently rotating by θ about the updated n (already rotated by φ). The resulting new basis vectors are expressed in terms of the initial basis vectors and the scattering angles (θ, φ), as follows: m = cos θ (cos φm + sin φn) − sin θk n = − sin φm + cos φn k = sin θ (cos φm + sin φn) + cos θk.
Note that the Stokes vector S in Equation (15) is automatically expressed in the new local reference frame constructed by (m , n , k ).

Initial Stokes and Basis Vectors
We now describe how to select the initial basis vectors for a newly emitted photon packet. When injecting a photon toward a direction of angles (θ, φ) in the galaxy reference frame, the initial basis vectors can be conveniently obtained by rotating m = (1, 0, 0) T , n = (0, 1, 0) T , and k = (0, 0, 1) T by the angles (θ, φ) using Equation (16). Consequently, the initial basis vectors for the angles (θ, φ) are chosen to be: Note that Peest et al. (2017) use only one polarization basis vector referred to as n and choose an initial basis vector to be n Peest 0 = −m 0 in their Equations (38)-(40). The two polarization basis vectors of Ahn & Lee (2015) can be expressed in terms of our definition, such that e 1 = n 0 and e 2 = m 0 .
In this study, we assume, unless otherwise stated, that initial input photons are unpolarized, i.e., I = 1 and Q = U = V = 0. Alternatively, unpolarized light can also be modeled by adopting 100% linearly polarized, monochromatic photons (P L = 1). In this case, however, the polarization angle χ should be chosen uniform-randomly in the range 0 ≤ χ ≤ 2π. This approach of using 100% polarized photons is equivalent to that of  if the simulation is limited to the Rayleigh scattering (see Appendix E for further discussion).

Scattering Matrix for the Scattering by Hydrogen
The scattering matrix for the scattering of Lyα photons by hydrogen atoms is expressed by where as described in Appendix A. Here, E 1 , E 2 , and E 3 are functions of frequency given by The second equation for E 2 is always true for the transition between a singlet ground state and an upper doublet state with arbitrary angular quantum numbers. But, the third equation between E 3 and E 1 is satisfied only for the transitions occurring between the same quantum levels as for Lyα. The matrix elements S 11 , S 12 , S 22 , and S 33 have been described in Stenflo (1980). The matrix element S 44 for the circular polarization is newly derived in Appendix A.
The new Stokes parameters after a scattering event by a hydrogen atom are given by I = S 11 I + S 12 (Q cos 2φ +U sin 2φ) Q = S 12 I + S 22 (Q cos 2φ +U sin 2φ) The last equation for V implies that the scattering of Lyα by hydrogen atoms does not newly produce any circular polarization if there was no initial circular polarization; but, it will remain once it is created. The scattering matrix for the electron Thomson scattering and Rayleigh scattering is obtained by setting E 1 = 1, E 2 = 0, and E 3 = 1.
The parameter E 1 represents the effect of quantum interference between the fine-structure levels in resonant Lyα scattering. It determines the shape of the scattering phase function as a function of frequency. As shown in Figure 2, it is intriguing that E 1 < 0 for a wavelength range of λ c < λ < λ H , where λ c = (2λ K + λ H )/3. However, in most wavelengths, E 1 is positive. From Equation (21), the scattering phase function for initially unpolarized light (Q = U = 0) is obtained to be Degree of linear polarization for E 1 = −1/2 and E 1 = 1. Here, P * L is defined in Equation (23) to allow to have a negative value to explain the direction of polarization vector. In a spherical system where a point light source is located at the center, a negative P * L indicates the tangential polarization, and a positive P * L the radial polarization.
The phase functions, for instance, for E 1 = −1/2 and E 1 = 1, are shown in the top panel of Figure 3. In most cases (E 1 > 0), scattering occurs dominantly in the forward and backward directions. However, for a negative E 1 , scattering occurs more frequently in perpendicular directions to the incident direction. For a positive E 1 , the phase function is a linear superposition of Rayleigh scattering and isotropic scattering with weights of E 1 and 1 − E 1 . An interesting thing to note is that when E 1 < 0, scattering events give rise to a counterintuitive, radial polarization pattern in the case of a central light source; this type of polarization is referred to as the "negative polarization." In most circumstances, the polarization pattern made by a central source is concentric, in other words, perpendicular to radial directions from the source. The counterintuitive, radial polarization pattern for E 1 < 0 can be demonstrated by examining the degree of linear polarization of the scattered light (for initially unpolarized light, Q = U = 0), which is defined to allow to have a negative value: Thus, a negative E 1 yields a positive Q , meaning a radial polarization pattern (|E m | 2 > |E n | 2 in Equation (11)), as shown in the lower panel of Figure 3. We also note that, regardless of the sign of E 1 (or polarization direction), the maximum polarization always happens at the right-angle scattering case.
In most astrophysical cases, a substantial diffusion of Lyα photons in frequency space is expected because of a large number of resonance scatterings. Therefore, most photons will escape the system only after scattering events occurring at the wing regime where the Rayleigh scattering is predominant and E 1 is positive. In the end, such a radial polarization (negative polarization) pattern that happens when E 1 < 0 is unlikely to be observable in actual situations.
2.6. Scattering Matrix for the Scattering by Dust For scattering by spherical dust particles, the scattering matrix has a form of The matrix elements for dust scattering cannot be expressed as an analytical function of θ, unlike the scattering case by hydrogen atoms. Instead, they should be numerically estimated using the Mie theory for a given grain size distribution and material composition. In general, all 16 elements of the scattering matrix for non-spherical dust grains will not be vanishing. Non-spherical dust is beyond the scope of this paper; however, the present work can be relatively easily extended to non-spherical cases. The scattering matrix elements calculated for the Milky-Way (MW), Small Magenellic Cloud (SMC), and Large MC (LMC) dust models of Weingartner & Draine (2001) are described in Appendix B.
In this paper, we adopt the MW dust model. The gas-to-dust ratio is fixed to that of the MW (∼ 100 by mass).
The above equation for V indicates that circular polarization can arise naturally in dust scattering even without initial circular polarization (V = 0). We also note that the presence of circular polarization can modify the Stokes U parameter and thus the degree of "linear" polarization. Therefore, we need to consider the possibility of the circular polarization of Lyα in dusty media, at least for theoretical completeness. The effect of circular polarization is, however, likely to be negligible in most circumstances.
2.7. Random Sampling of Scattering Angles Regardless of whether the scattering partner is a hydrogen atom or dust grain, the scattering phase function, which gives the probability distribution function to be scattered by angles (θ, φ), is given by Here, S 11 is assumed to be normalized, i.e., S 11 sin θdθ = 1; all other scattering matrix elements are also normalized by the same factor as that for S 11 . One wants to generate two random variables θ and φ simultaneously, according to the joint distribution function P(θ, φ). It is, however, more convenient to obtain θ first and then φ, instead of generating (θ, φ) simultaneously. For this purpose, we first marginalize P(θ, φ) over the azimuth angle φ and obtain a random θ from the marginal distribution of θ: For a chosen θ, the conditional probability distribution function of φ is given by A random variate for φ from this distribution function can be obtained using a rejection method. We first choose a random azimuth angle φ = 2πξ from a uniform random number ξ between 0 and 1. We accept it if ξ ≤ P(φ|θ)/P max for a new independent random variate ξ (0 ≤ ξ < 1). If not, then generate a new random φ until it is accepted. Here, the maximum valueP max of P(φ|θ) is given bȳ The procedure is efficient in that more than 50% (in fact, most) of φ is accepted. An alternative sampling strategy is to choose a scattering angle θ according to P(θ) and draw a "biased" azimuth angle φ using a uniform distribution in the range of [0, 2π), i.e., φ = 2πξ. In this case, we need to adjust the weight of the photon packet to w = wP(φ|θ)/2π, where w and w are weights for the incident and scattered photon packet, respectively. The present study uses the approach of using a rejection method for φ.
For the case of scattering by hydrogen atoms, the scattering angle θ between the incident and scattered directions is governed by the following probability distribution function: where µ ≡ cos θ. We use the inversion method to generate random scattering angles from this distribution function. In other words, a random angle θ is given by inverting the integral for a uniform random number ξ. The integral yields a cubic equation of µ for given ξ and E 1 , which can be solved using Cardano's method. The analytical solution of the cubic equation is summarized, for an arbitrary E 1 (−1/2 ≤ E 1 ≤ 1), as follows: where The solution is detailed in Appendix C. For the scattering by dust grains, the scattering phase function for θ is well approximated by a Henyey-Greenstein (H-G) function or a combination of two independent H-G functions (Witt 1977). The integral of the H-G function is analytically invertible, as described in Witt (1977); a scattering angle θ that follows a H-G function can be obtained by If one adopts a sum of two H-G functions, they can use a composition method to generate scattering angles (θ). Instead of using an analytical approximation for the phase function, the numerical inversion technique can also be used to draw the random scattering angle θ. The scattering matrix elements, including the phase function, for the scattering by dust are described in Appendix B. In the appendix, we also provide approximate equations for the matrix elements of Equation (25).
2.8. Update of the Stokes vector After each scattering event, we update the Stokes vector as follows: The Stokes vector is updated to have a unit intensity because the photon packet carries unit intensity in the Monte-Carlo simulation. When calculating the new Stokes vector S after a scattering event by dust, it is more convenient to use the scattering matrix elements divided by S 11 rather than the elements themselves, as described in Appendix B. If we designed for photon packets to have a weight (especially for the dust scattering case), we also need to update the photon weight appropriately and multiply the Stokes vector by the same factor.
2.9. Peeling-off Technique A Monte Carlo RT simulation aims to obtain a spectral (or spectro-polarimetric) image on a detector to compare with observational data. However, in a simple Monte Carlo simulation, there will be hardly a probability to detect a statistically significant number of photon packets on a virtual detector located far from the system. We, therefore, use the "peeling-off" (also known as "next event estimation" or "shadow rays") technique to obtain high signal-to-noise images in a detector plane (Yusef-Zadeh et al. 1984;Laursen & Sommer-Larsen 2007;Yajima et al. 2012). When using this technique, the probability of a photon escaping in the direction of the observer in every scattering (and emission) event is added to the detector plane. The "peeling-off" technique has been implemented in many Lyα RT codes, but, surprisingly, was not well described. Here, we briefly describe a peeling-off strategy implemented in LaRT. . Coordinate systems of the galaxy (or hydrogen+dust cloud, denoted by G) and the observer (O). The coordinates of the observer measured in the galaxy frame are given by r obs = (x obs , y obs , z obs ) = (d cos α sin β, d sin α sin β, d cos β). The angle γ defines the orientation of the detector plane.

Detector Plane
LaRT can place an observer in an arbitrary location and make a detector plane to have an arbitrary orientation in the sky. As shown in Figure 4, the location of the observer is defined by a luminosity distance (d) and two rotation angles (α, β), or equivalently, by the coordinates x obs , y obs , and z obs of the observer in the galaxy frame; a third angle γ then defines the orientation of the detector plane. The observer frame, given by (x O ,ŷ O ,ẑ O ), is assumed to initially coincide with the galaxy frame defined by (x G ,ŷ G ,ẑ G ). Next, the observer frame is obtained by sequentially rotating the initial observer frame by α aboutẑ O , by β about the rotated axisŷ O , and then by γ about the rotated axisẑ O . 8 The detector is then located at the coordinate (0, 0, d) in the observer frame. Therefore, the components of a direction vector in the observer frame can be obtained by multiplying its components in the galaxy frame by the following rotation matrix: where To perform the "peeling-off" of a photon packet, we imagine a "virtual" scattering event in the direction of the observer and construct a propagation vector (k peel = (r obs − r ph )/ |r obs − r ph |) of the "peeled-off" photon by connecting the photon (r ph ) to the observer (r obs ). The propagation vector is then transformed from the galaxy frame to the observer frame. The celestial coordinates of the "peeled-off" photon in the detector plane are calculated to be θ and k O z are the coordinate components of k peel in the observer frame. The celestial coordinates are finally binned into a two-dimensional array to obtain an image of a desired observable quantity in the detector plane.
2.9.2. The Stokes vector of peeled-off photon To calculate the Stokes vector carried by a peeled-off photon packet, let us consider a photon packet peeled off into a direction (k peel ) pointing toward the observer. We first need to find the scattering angles θ peel and φ peel for this "virtual" scattering event. For this purpose, the basis vectors (n peel and k peel ) of the peeled-off photon are expressed in terms of the initial basis vectors (m, n, k) and the scattering angles (θ peel , φ peel ), as follows: where n peel , θ peel , and φ peel are the unknowns supposed to be found using (m, n, k) and k peel . Here, it is not necessary to derive m peel = n peel × k peel . From Equation (39), the scattering polar angle θ peel is readily calculated via the scalar product of the incoming and "virtual" scattering direction vectors: which uniquely determines the scattering angle in the range of 0 ≤ θ peel ≤ π. To calculate the azimuth angle φ peel unambiguously, both cosine and sine of the angle are required. They can be obtained from Equation (39): cos φ peel = k peel · m/ sin θ peel sin φ peel = k peel · n/ sin θ peel .
The basis vector n peel can then be obtained by Equation (38). The Stokes vector of the peeled-off photon is given by S = M(θ peel )L(φ peel )S. In the above, the detector orientation was not yet taken into account. In other words, m peel and n peel do not coincide with the axesx O andŷ O of the detector plane (see Figures  4 and 5). We, therefore, need to rotate m peel and n peel by an angle φ det between m peel andŷ O (or between n peel and −x O ). Here, recall that the IAU standard measures a polarization angle counter-clockwise from the North in the sky (ŷ O in the observer frame) when looking at the photon source; hence, φ det is defined as an angle between m peel and y O . By this definition, the desired rotation angle φ det is uniquely determined by It should also be noted that n peel in Equation (38) was calculated in the galaxy frame. Therefore, we need to represent the components of n peel in the observer frame by multiplying the rotation matrix R, defined in Equation (37). The final Stokes vector, represented in the detector plane, is then given by Detector Plane Figure 5. A local coordinate system defined by the radial (er) and tangential (et ) vectors at a projected radius rp in the detector plane. A linearly-polarized electric field vector (E) and its polarization angle (χ) in the local coordinate system are also shown.
which follows the IAU standard.

Frequency of peeled-off photon
To perform raytracing in the "peeling-off" procedure of scattered light, the optical depth along the path to the observer should be calculated at the frequency in the comoving fluid frame. On the other hand, the spectral binning on the detector plane should be performed using the frequency measured in the fixed (observer) frame. Therefore, we need to consider two frequencies, the frequency x tau (for optical depth) in the fluid frame and the frequency x spec (for spectral binning) in the observer frame. If the photon is scattered by a hydrogen atom with a velocity u atom , they are given by Here, x i is the initial frequency (before the scattering event) expressed in the fluid frame, u atom, = k i · u atom the velocity component parallel to the initial direction k i , x i − u atom, the frequency in the rest frame of the scattering atom, k peel the direction vector pointing to the observer, and u fluid the fluid velocity. In the case of scattering by a dust grain, x tau and x spec are given by setting u atom = 0 in the above equation because dust grains are heavy and assumed to have no thermal motion. For the peeling-off of direct light, which escapes the system without undergoing any scattering, we consider the two cases (injection from the comoving fluid frame and a non-comoving fixed frame) separately. In most cases, photons are injected in the fluid frame. LaRT also has another option to inject photons in the fixed frame rather than in the comoving frame. The former option is useful to simulate emissions from the gas. The latter is suitable for the emission originating from an object, such as stars/galaxies, that independently moves from the ISM/IGM. For the case injected in the fluid frame, x tau and x spec are given by where x fluid i is the initial photon frequency in the fluid frame. On the other hand, when a photon is injected in the fixed frame, they are given by where x fixed i is the initial frequency in the fixed frame. In this case, the comoving frequency of a photon injected in the fixed frame is x fluid i = x fixed i − k i · u fluid for the initial propagation vector k i .

Peeled-off fraction
The fraction to be peeled-off toward the detector is Here, τ is the optical depth at the photon's current frequency, measured over the path from the scattering (or injection) position to the observer, and d the luminosity distance. The phase function Φ is normalized to one over the solid angle of dΩ = 4π. This paper assumes isotropically-emitting Lyα sources. In this case, the phase function for direct light, which undergoes no interaction with the medium, is The phase function for scattered light is given by Equation (26), as follows: The peeled-off Stokes parameters should also be multiplied by the same factor as for the peeled-off intensity.

TWO FACTORS THAT DETERMINE POLARIZATION PROFILE
Before describing the details, we here first summarize the essentials to understand the simulation results. We explore the intensity vectors polarization vectors Figure 7. Rise of the radiation field anisotropy and the polarization with increasing radius. The left panel illustrates the rays passing through three different locations in the projected plane when viewed normal to the paper. Each blue arrow denotes a ray of photon packets (intensity) propagating in the direction represented by the arrow. The size of the blue arrows represents the intensity or the number of photons passing through a unit area into a unit solid angle of the given ray. The flux passing through the surface normal to the radial direction increases with radius because of the increase in anisotropy. In the right panel, the red arrows at each location denotes the polarization vector carried by the scattered light of the rays passing through that location toward the observer. In the projected plane, the polarization vector of the scattered light of a ray is perpendicular to the direction of the intensity ray, unless the scattering occurs in the core, and its length is proportional to the polarized intensity of the ray. The degree of polarization is given by P where Er and Et are the radially and tangentially oscillating electric vectors, respectively. Therefore, the ensemble-averaged polarization is dominated by the contribution of the outward rays and in general increases with radius. However, note that the core scattering is cable of producing a signifincantly different radial profile from that shown in this figure, as is discussed in Sections 4.2 and 4.4. degree of polarization calculated in two different ways. One is the quantities averaged over an "ensemble" of photon packets, which experienced many "different" scattering histories, measured at a location in the detector plane. In this case, the Stokes parameters are first averaged, as defined in Equation (11), and then the degree of polarization is calculated using the averaged Stokes parameters. Dijkstra & Loeb (2008) refer to these quantities as the "angle-averaged" ones. However, we prefer to refer to them as the "ensemble-averaged" ones. The other is the polarization quantities of "individual" photon packets, which represent the quasi-monochromatic waves undergone the "same" scattering history. When calculating the mean of the polarization amplitudes for individual photons, we take the average ignoring the orientation of the individual polarization vectors. The ensemble-averaged quantities are observable through an experiment, while those of individual photon packets are not. Dealing with individual photon packets is theoretically advantageous to interpreting simulation results. We distinguish the two concepts if necessary to clarify. However, the term "ensemble-averaged" is often omitted when it is evident in context. We also use "photons" interchangeably with "photon packets." The ensemble-averaged degree of polarization is primarily determined by two factors: the polarization amplitudes of individual photon packets and the degree of anisotropy in the Lyα radiation field. . Radial profiles of the surface brightness and polarization for a static sphere model of T = 10 K (the first and second panels) and T = 10 4 K (the third and fourth panels). The central optical depth varies from τ 0 = 10 3 to 10 7 . The column density N HI of each model is given by the equations shown in the first and third panels.
the polarization amplitudes of individual photons are, the higher their ensemble average (denoted by red arrows) will result. Figure 6(b) illustrates the importance of anisotropy of the radiation field in polarization. In the figure, the lengths of arrows are the same in both panels, while the diversity of direction angle indicates the degree of anisotropy (or isotropy) of the radiation field. Suppose the radiation field is highly isotropic (as in the right panel). In that case, two perpendicular components of the polarization vectors of individual photon packets will significantly be canceled out, resulting in a weak polarization. On the other hand, in a highly anisotropic radiation field (left panel), they will be less canceled out, yielding a strong polarization. The polarization levels of individual photons generally tend to increase with the number of scatterings and the strength of the medium's velocity field, as detailed in Sections 4 and 5.1. The degree of the radiation field isotropy rises with an increasing number of scatterings but decreases in fast-moving media. In a spherical symmetric medium with a central source, it is found that, as the optical depth of the medium (and the number of scatterings) increases, the polarization amplitudes of individual photon packets begin to grow first while the radiation field is kept still anisotropic. Then, later, the radiation field appears to become isotropic after the individual photons gain high enough polarization amplitudes.
It is also found that the radial polarization profile increases monotonically with increasing radius, except for some cases in Section 4.4. The increase of polarization degree with increasing radius is attributable to the rise of anisotropy in the radiation field. As illustrated in Figure 7, the radiation field near the center of a sphere will be isotropic by symmetry. On the other hand, at large radii, the outward radiation flux will dominate over the inward one, leading to anisotropy in the radiation field. Eventually, this anisotropy in radiation at large radii results in an increase in polarization with radius.
In the next section, we also show that the steepness (or slope) of the surface brightness profile correlates well with the steepness of the polarization profile. As optical depth increases, the number of scatterings per photon increases, and the radiation field becomes more isotropic due to the spatial diffusion of photons. Therefore, the surface brightness will tend to flatten with the increase of optical depth. In other words, a highly isotropic radiation field eventually leads to a flatter (or shallower) surface brightness profile. On the other hand, an anisotropic radiation field would produce a steep surface brightness profile. Therefore, we can immediately associate the steepness of the surface brightness and polarization profile with the degree of anisotropy in the radiation field; a steep surface brightness profile tends to accompany a rapid increase of polarization.

RESULTS
This section presents the results of four types of models in total. We first present two basic models as the first and second examples. One is a static spherical medium with various optical depths, adopted to examine the "optical depth effect" on Lyα polarization. The other is a Hubble-like sphere expanding at many different velocities to illustrate the "kinematic or velocity effect." Next, we reproduce the results for the thin shell model of Dijkstra & Loeb (2008). This model not only verifies our code but also illustrates a special case where single-scattering is predominant. Lastly, we predict the polarization signal from the galactic halo model of Song et al. (2020), which was utilized to fit the observed spectra and surface brightness profiles of LAEs. This fourth example demonstrates the diversity of the polarization pattern.
In the examples, we deal with a spherically symmetric geometry. Therefore, it is convenient to represent the Stokes parameters in a local coordinate system defined by the radial and tangential vectors (e r and e t , respectively) at the photon's location in the detector plane, as shown in Figure 5. We use this coordinate system to calculate the radial profile of polarization. The polarization angle χ is also measured in this coordinate system. For this purpose, the basis vectorŝ x O andŷ O of the detector plane are rotated to make them coincide with the radial and tangential vectors, respectively. In the local coordinate system, Q = −1 denotes a polarization vector perpendicular to the radial direction, corresponding to a concentric polarization pattern. On the other hand, Q = +1 represents a polarization vector parallel to the radial direction. Here, it should be noted that the Q = +1 shown in this section is caused by photons traveling perpendicular to the radial direction, not by the negative polarization (E 1 < 0) discussed in Section 2.5; in other words, the Q = +1 case occurs due to the randomness of photon directions caused by multiple scattering.
This study presents the results for the media of T = 10 or 10 4 K. In the first three examples, photons are injected from the center of the sphere. In the last example, the spatial distribution of input Lyα photons is assumed to be described by an exponential function of projected radius. The input line profile is assumed to be a Voigt function with a line width corresponding to the medium's temperature, except for the outflowing thin shell model of Section 4.3 (see Seon & Kim 2020 for the random number generation algorithm of the Voigt profile). The width of the input line profile does not significantly affect the present results unless it is as wide as those observed in active galactic nuclei or the output spectra shown in this paper. The dust effect on the Lyα polarization is discussed in the first and third examples (the static, homogeneous sphere and the outflowing thin shell models). The fourth example (the galactic halo model of Song et al. 2020) also contains dust; however, no discussion on the dust effect is given for this example. The gas-to-dust ratio is assumed to be that of the MW. A Cartesian grid is utilized in the present version of LaRT. The number of cells in the simulations ranges from 200 3 to 400 3 .

A Static, Homogeneous Sphere
In this example, we investigate the optical depth effect on the polarization property using a static, spherical medium. The optical depth at the line center, measured from the center to the outer edge of the sphere, varies from τ 0 = 10 3 to 10 7 , and the gas temperature is T = 10 K or 10 4 K. Figure 8 shows typical images of the Stokes parameters I, Q, and U, obtained for a model of τ 0 = 10 6 and T = 10 4 K. The figure illustrates the concentric polarization pattern that is generally expected in a configuration with a central source. Figure 9 shows the surface brightness and the degree of polarization as a function of radius. In the figure, we notice the following. First, the surface brightness profile is relatively steep in small optical depths (τ 0 10 4 for T = 10 K and τ 0 10 5 for T = 10 4 K) but flattens quickly as τ 0 increases. Second, as commonly expected, the radial profiles of polarization follow the general trend of increasing with radius. Third, as τ 0 increases to a certain critical value (τ 0 ≈ 10 4 for T = 10 K and 10 5 for T = 10 4 K), the overall polarization level increases rapidly but decreases as τ 0 rises further above that value. Fourth, the radial profile of polarization jumps at r p ≈ 0 suddenly from zero to a non-zero value if τ 0 is relatively small (τ 0 10 4 for T = 10 K and τ 0 10 5 for T = 10 4 K); but, it gradually increases with radius if τ 0 is large. In other words, the polarization profile is steep near r p ≈ 0 for a small τ 0 , while it is shallow for a large τ 0 . Lastly, both the surface brightness profile and the polarization profile appear to begin to change in shape at the same critical optical depth (τ c ≈ 10 4 for T = 10 K and 10 5 for T = 10 4 K), indicating that they are closely associated with each other. In the cases where the optical depth is higher than this value, a steeper surface brightness profile tends to be accompanied by a higher overall level of polarization and a slightly steeper polarization profile.
The flattening of surface brightness with increasing τ 0 (the first result in Figure 9) is caused by the trend that, as τ 0 increases, the number of scatterings increases, and the radiation field becomes more isotropic (see Section 3). The increase of polarization degree with radius (the second result) is due to the enhancement of anisotropy in the radiation field, as illustrated in Figure 7. The third and fourth results in Figure  9 can be understood with the help of polarization properties of individual photon packets. Figure 10 shows the probability density functions of the polarization degree P L , the Stokes parameters Q and U, and the polarization angle χ * = χ − π/2 of individual photon packets, which were calculated at various radius bins, for the model of T = 10 K 9 . To display the distribution of polarization angle conveniently, we use χ * , measured counterclockwise from the tangential vector defined in Figure 5; thus, χ * = 0 denotes the tangential direction. The results for the model of T = 10 4 K are shown in Figure 11.
In this present study, the initial input photons are injected unpolarized. Hence, the variation in the distribution of P L with τ 0 implies an interesting fact that individual photon packets are more polarized as they are scattered more, which are further discussed in Section 5.1. Figures 10 and 11 show that, for relatively low optical depths (τ 0 10 4 − 10 5 ), the rise of polarization levels of individual photon packets with increasing τ 0 is responsible for the increase of ensemble-averaged polarization. However, as τ 0 increases further, the photon packets are fully polarized, and no further increase of the individual polarization levels occurs. Instead, the radiation field becomes more isotropic with increasing τ 0 , leading to the decrease of the ensemble-averaged polarization due to the cancellation of individual polarization vectors. This polarization trend explains the third result in Figure 9 that, as τ 0 increases, the overall polarization level increases at first and then decreases later. In the models with relatively low optical depths, a sudden breakdown of the isotropy in the radiation field at the location right adjacent to the center causes a sudden increase of the ensemble-averaged polarization at that location. In contrast, in the models with relatively high optical depths, the strong isotropy in the radiation field near the center largely cancels the individual polarization vectors, even though most individual photon packets bring ∼100% 9 Figures 10, 11 , 15, 16, 19, 20, 24, and 25 describe the statistical polarization properties of individual photon packets escaping the system. In the figures, the Stokes parameters were calculated in the local reference frame, defined by the last direction vector, at the last scattering position. The peeling-off technique was not utilized in the figures.  . Probability density functions of the degree of polarization P L , polarization angle χ * , and the Stokes Q and U of individual photon packets for a static sphere of T = 10 K. The model optical depth varies from τ 0 = 10 4 (the first row) to 10 7 (the fifth row). The figure shows the distribution functions for several radial bins (0 < rp < 0.1, 0.1 < rp < 0.3, 0.3 < rp < 0.5, 0.5 < rp < 0.7, 0.7 < rp < 0.9, and 0.9 < rp < 1), as denoted in the upper left panel. We use χ * = χ − π/2, measured counterclockwise from the tangential direction, for polarization angle, instead of χ, to make its histogram peak at χ * = 0. In this definition, χ * = 0 represents a tangential or concentric polarization pattern, and χ * = ±π/2 a radial polarization pattern.  . Probability density functions of P L , χ * , and the Stokes Q and U of individual photon packets for a static sphere of T = 10 4 K. The model optical depth varies from τ 0 = 10 4 (the first row) to 10 7 (the fifth row).
polarization. Therefore, the fourth result in Figure 9 is also understood.
More detailed descriptions of the polarization properties of individual photons are as follows. In the models with small optical depths (τ 0 < 10 4 for T = 10 K and τ 0 < 10 5 for T = 10 4 K), most photon packets appear to have P L < 0.5 (the first column in the figures). However, as the optical depth (and the number of scatterings) increases, most photon packets tend to be fully polarized (P L = 1). We also find a weak but noticeable bump shape with a peak at P L ≈ 0.4 in the distribution function of P L for the model of T = 10 4 K; this bump shape is associated with the pot lid shape in the distribution function of Q found in the range |Q| 0.5 (also seen in U). It diminishes as τ 0 increases or T decreases; it is not noticeable in the models of T = 10 K, except for τ 0 = 10 3 , as in Figure 10. In the second column of the figures, the distribution function of polarization angle shows a peak at χ * = 0 when τ 0 is small, while it flattens as τ 0 increases. In  Figure 12. Polarization degree for a static homogeneous spherical medium. The first and second rows show the results for T = 10 K and T = 10 4 K, respectively. The optical depth ranges from τ 0 = 10 4 to 10 7 . In the small ancillary panels, the black dotted, red solid curves are the initial and emergent spectra, respectively. The green curves are polarization degrees as a function of frequency. The main panels show the linear polarization degree in a two-dimensional space of photon frequency and projected radius. Here, the polarization degree is the "ensemble-averaged" one for photons detected at radius rp with frequency x in the detector plane. The color bar denotes the polarization levels.  Figure 13. Dust effect on the radial profiles of the surface brightness and polarization for a static spherical medium of the H I optical depth τ 0 = 10 7 (N HI = 5.5 × 10 18 and 1.7 × 10 20 cm −2 for T = 10 and 10 4 K, respectively). The optical depths due to dust are τ dust = 8.8 × 10 −3 and 2.7 × 10 −1 for T = 10 and 10 4 K, respectively. The blue and red colors denote the model with a gas temperature of T = 10 4 K and 10 K, respectively. The solid lines and dots represent the model without and with dust, respectively. particular, for τ 0 = 10 7 , the distribution function of χ * at r ≈ 0 is almost completely uniform. The results imply that at small radii or as τ 0 increases, the polarization orientations tend to be completely randomized, and the radiation field becomes isotropic. Therefore, even though photon packets were almost fully polarized for a sufficiently large τ 0 , the polarization vectors are largely canceled out due to the isotropy of the Lyα radiation field. The third column shows that the distribution function of Q is dominated by Q = 0 when τ 0 is small, but its peak at Q = 0 diminishes and consequently disappear, as τ 0 increases. We also note that as τ 0 increases, the peak at Q = −1 (concentric polarization) grows first, and then the peak at Q = 1 (radial polarization) follows; eventually, the distribution function of Q becomes symmetric and has two sharp peaks at Q = ±1. This trend is because, as τ 0 increases, the polarization degree of individual photon packets increases first, and then the radiation field becomes isotropic later. The peak at Q = −1 rises as the polarization degree of individual photons increases (but while the radiation field remains strongly anisotropic). As τ 0 increases further, the peak at Q = 1 also begins to grow as a result of the enhanced isotropy of the radiation field. The distribution of U is also dominated by U = 0 for small τ 0 but by |U| = 1 for large τ 0 . Unlike for Q, the distribution function of U is always symmetric in all models discussed in this paper due to the spherical symmetry of the models and by definition in Equation (11). In summary, in the cases with a relatively low optical depth (τ 0 10 4 − 10 5 ), the polarization amplitudes of individual photons are the main factor that governs the radial profile of the ensemble-averaged polarization. On the other hand, if the optical depth is high (τ 0 10 4 − 10 5 ), the degree of isotropy The column density of the medium is N HI = 2 × 10 19 cm −2 (the first and second panels) or 2 × 10 20 cm −2 (the third and fourth panels).
(or anisotropy) in the radiation field is the prime factor that determines the polarization. It is also evident that an isotropic radiation field will result in a shallow or flat surface brightness profile, while an anisotropic radiation field would produce a steep surface brightness profile. This association of the surface brightness profile with the anisotropy of the radiation field explains the last finding in Figure 9. Figure 12 presents the dependence of polarization on the photon frequency and the projected radius, which may be useful to compare with the spatially resolved spectropolarimetric observations. The small ancillary panels show the frequency dependence of the polarization degree averaged over the whole spatial area. The polarimetric spectra are symmetric between the red and blue wavelengths because this model is static. The strongest polarization is found in the line wings of both red and blue sides. This is because in the wings, not only can the photons easily escape, but also they tend to preserve the polarization level of 100% during scattering events, as discussed in Section 5.1. On ther other hand, the core photons tend to be depolarized. In more detail, the radial profile of polarization is frequency-dependent; the farther the frequency is from the line center, the steeper the radial polarization profile. This is because, at frequencies farther from the line center, the Lyα photons tend to escape more easily, making the radiation field more anisotropic. Dijkstra & Loeb (2008) argued that the presence of dust will boost the polarization level as a result of preferential extinction of multiply-scattered photons. A photon scattered multiple times passes a longer length in the medium, leading to a higher possibility of absorption by dust. In other words, the presence of dust grains causes less spatial diffusion of Lyα photons by destroying them, particularly the multiple-scattered ones. We also note that the central parts of the Lyα line are mainly absorbed by dust grains, as shown in Figures 5 and D1 of Seon & Kim (2020). As a result, the presence of dust leads to a slightly steeper surface brightness profile and a more anisotropic Lyα radiation field. This effect is expected to strengthen the polarization level eventually. The dust effect will be more critical with higher τ 0 or N HI . We therefore examined the dust effect on the polarization and surface brightness profiles by assuming the MW dust, as shown in Figure 13. The polarization level for the model with a gas temperature of T = 10 4 K (τ 0 = 10 7 ) is found to increase slightly by the presence of dust. On the other hand, no noticeable change is found in the case of T = 10 K. The figure also shows that the reduction of surface brightness level by dust is more profound in the model with a high temperature (T = 10 4 K). This is because, as noted in Seon & Kim (2020), the ratio of dust extinction to the H I scattering cross-section is proportional to T 1/2 at the Lyα line center. Therefore, the effect is stronger when the gas temperature is higher.

A Hubble-like Flow
We now investigate the velocity effect on the polarization properties of Lyα using a Hubble-like expanding medium. The Hubble-like flow model has been used as a benchmark test to predict the emergent spectra in a moving medium (Laursen et al. 2009;Yajima et al. 2012;. This type of model is also proven to be useful in interpreting the Lyα observations. For instance, to explain the Lyα halos and low-ionization metal absorption lines in Lyman break galaxies, Steidel et al. (2010) suggest a simple model in which Lyα photons are scattered in a spherical symmetric clumpy medium of outflowing at velocity increasing with distance from the galaxy.
In this model, we consider an isothermal and homogeneous sphere, which is isotropically expanding or contracting. The velocity of a fluid element is assumed to be proportional to the radial distance: V (r) = V exp r/r max . Here, V exp is the maximum expansion velocity, representing the velocity gradient in the radial direction. The neutral hydrogen gas has a temperature of T = 10 4 K and a column density of N HI = 2 × 10 19 cm −2 or 2 × 10 20 cm −2 (measured from r = 0 to r = r max ). The sphere expands radially with a maximum velocity ranging from V max = 20 to 3000 km s −1 at the outer edge of sphere (r = r max ). The emergent spectra for a few cases with N HI = 2 × 10 20 cm −2 were presented in Seon & Kim (2020). Figure 14 shows the radial profiles of the surface brightness and polarization obtained from this model. The main results are summarized as follows. First, the polarization rises monotonically with radius. Second, the surface brightness profile is nearly flat in the case of the lowest velocity V max = 20 km s −1 , as in the static model, while it becomes steeper as V max increases. This trend is because the number of scatterings, which photons undergo before escape the system, decreases as the medium's expansion velocity increases. Third, the overall level of polarization is found to be in general higher for the model with a higher expanding velocity. Fourth, the degree of polarization in fast-moving models (V max 1000 km s −1 ) is found to reach almost up to 100% at the boundary. Fifth, in the fast-moving models, the polarization increases abruptly at r ≈ 0, while in rather slowly moving models, it rises gradually. Lastly, the steepening of surface brightness profile and the polarization jumping in the highest velocity models appear to be more prominently when the column density is lower. Similar results were obtained in contracting media with a negative V exp but not presented here.
In a medium moving away from the source, most photons will be recognized to have longer wavelengths in the rest frame of the gas (e.g., Figure 15 in Seon & Kim 2020). Then, the photons are predominantly scattered by the wing (Rayleigh) scattering, increasing the probability of being 100% polarized, as discussed in Section 5.1. Moreover, photons in a fast-moving medium would experience a relatively small number of scatterings. Hence, the resulting radiation field will be highly anisotropic, and the surface brightness profile would be relatively steep (the second result in Figure 14). If once a photon packet acquires 100% polarization through the Doppler shift after only a small number of scatterings, the polarization level of the photon packet would be kept high. Also, the ensemble-averaged polarization level will be very high because of strong anisotropy in the radiation field (the third and fourth results). In contrast, in a slow-moving medium, photon packets will have yet a high probability of being scattered in the core regime; thus, their polarization will repeat the process of decreasing and increasing during multiple scatterings until eventually scattered into the wing regime. This process gives rise to a relatively lower probability of achieving 100% polarization than in fast-moving media. More importantly, the radiation field will be much more isotropic than in fastly moving media, resulting in a flatter surface brightness and a lower ensemble-averaged polarization.
In the end, the model with a faster expansion velocity (e.g., V epx = 200 km s −1 ) will achieve higher polarization than the model with a lower velocity (V exp = 20 km s −1 ). However, if the expansion speed is even faster (V exp 500 km s −1 ), the chance of being scattered is lowered, and thus the fraction of photons that escape the medium after being fully polarized rather decreases, compared to that for V exp = 200 km s −1 . The same effect arises also in the emergent spectrum; the peak in the spectrum for the model of V exp 500 km s −1 is less redshifted (Figure 3 of Seon & Kim 2020 and the small ancillary panels of Figure 17), than the model of V exp = 200 km s −1 . Nevertheless, owing to a strongly anisotropic radiation field, the models of V exp 500 km s −1 give rise to higher levels of ensemble-averaged polarization than the model of V exp = 200 km s −1 . Given that, in the lower-velocity models (V exp 200 km s −1 ), more photon packets are 100% polarized, but the ensemble-averaged polarization is lower, it is evident that the anisotropy in the radiation field is of primary importance to yield a polarization level as high as those seen in the models of V exp 500 km s −1 .
The abrupt jump of polarization at r ≈ 0 in the fastest models (V max 500 km s −1 for N HI = 2 × 10 19 cm −2 and V max 1000 km s −1 for N HI = 2 × 10 20 cm −2 ) results from a sudden breakdown of isotropy in the Lyα radiation field, as in the static models with low optical depths (the fifth result). However, the jump in polarization occurs to a much higher level in an expanding medium than a static medium. Figures 15 and 16 show the polarization properties of individual photons for the Hubble-like expanding media with N HI = 2 × 10 19 cm −2 and 2 × 10 20 cm −2 , respectively. The difference in the jump levels between expanding models and static models can be explained by examining the figures. Figure 15 illustrates that the vast majority of individual photon packets in the expanding media with N HI = 2 × 10 19 cm −2 gains ∼100% polarization and the peak in the distribution function of P L at P L ≈ 0 is very weak or absent, except for the inner most region (0 ≤ r p ≤ 0.1) in fast-moving models. This is because the photons are rapidly Doppler shifted to the wings by the outflowing gas motion. When the column density is high (N HI = 2 × 10 20 cm −2 ; Figure 16), similar results are obtained; but, the probability of having low P L is much lower.  Figure 17. Polarization degree for a Hubble-like expanding medium. The first and second rows show the results for N HI = 2 × 10 19 cm −2 and 2 × 10 20 cm −2 , respectively. In the small ancillary panels, the black dotted and red solid curves are the initial and emergent spectra, respectively. The green curves are polarization degrees as a function of frequency. The main panels show the linear polarization degree in a two-dimensional space of photon frequency and projected radius. The color bar denotes the polarization levels. The thick black lines denote the peak frequency where the spectrum has maximum in each radial bin. The thick dotted lines indicate the frequency at which the spectrum has a half maximum value.
We now examine the difference due to the column density. In general, individual photon packets achieve less polarization in a lower column density model. However, for the same expansion velocity, a lower column density model yields a higher ensemble-averaged polarization because of stronger anisotropy in the radiation field, as shown in Figure 14. The anisotropy of the radiation field plays a more critical role than the polarization level of individual photons; the leading role of the anisotropy explains the last finding in Figure 14. The fact that fewer photons are fully polarized in lower column density models, especially in fast-expanding cases, is illustrated in the distribution function of Q of Figures 15 and 16. In Figure 15 (N HI = 2 × 10 19 cm −2 ), the peak at P L = 1 in the distribution function of P L drops rapidly with increasing V exp , together with the increase in the probability of being 0.5 < P L < 1; on the other hand, in the model of N HI = 2 × 10 20 cm −2 , the peak at P L = 1 drops only mildly. However, the peaks in the distribution functions of χ * and U in Figure 15 (the lower column density model) are much sharper than those in Figure 16 (the higher column density model), indicating that the radiation field in the lower column density model is more anisotropic.
More detailed descriptions of the polarization properties of individual photons are as follows. Figures 15 and 16 show that the overall shapes of the distribution functions of P L , Q, and U change gradually as the expansion velocity increases.
For instance, the shape of the distribution function of P L for the model of V max = 2000 km s −1 varies significantly with the projected radius r p , compared to that of V max = 20 km s −1 , which shows no significant variation. The models expanding faster than V max = 500 km s −1 show a sharp peak at χ * = 0 in the distribution function of χ * . Moreover, their distribution function of U appears to significantly differ in shape from that of Q, unlike the slowly-expanding or static models. The peak at U = 0 in the distribution function of U tends to grow as the expansion speed of the medium increases. A similar trend is also found in the case of static media, but in this case, it occurs as the optical depth decreases. However, the peak in the expanding models is fairly broad, compared to that found in the static models. Also, the peak in U is not accompanied by the peak in Q, except for the case of 0 ≤ r p ≤ 0.1 in the models of N HI = 2 × 10 19 cm −2 , unlike the static models in which the peaks in Q and U seem to coexist. These differences are mainly because the peak of P(U) at U = 0 in the static models with a low optical depth is caused by the dominance of P L = 0 (weak polarization of individual photons), while that in the expanding models is due to the distribution of χ, i.e., anisotropy of the radiation field. In Figure 15, the peak of the distribution function of Q at Q = 0 in the radial bins of r p 0.1 becomes bigger with increasing V exp , because of the lowering of the effective optical depth. This trend in the small radii can also be found by noticing that the probability of being P L 0.2 rises with V exp . On the other hand, for larger radii (r p 0.6), the probability of being P L 0.2 decreases gradually as V exp increases. In Figure 16, we find no apparent peak of the distribution function of Q at Q = 0, unlike the case of Figure 15. This difference is attributed to that photons in higher column density models gain higher polarization amplitudes. Figure 17 shows the change of the polarization degree when the expansion velocity varies, as a two-dimensional function of radius and frequency. The figure also shows the frequencies where the spectrum has maximum (thick solid lines) and half maximum (thick dotted lines) values in each radius bin. The model with the lowest expansion velocity shows no significant difference in polarization from the static media, except for a weak tendency of being more polarized in the red wavelengths. For the intermediate velocity (V max ≈ 200 km s −1 ), the polarization degree shows a relatively weak dependence on frequency in a frequency range of x −10 at a fixed radius (r 0.4); however, it strongly depends on the radius. The degree of polarization is in general weak at r ≈ 0. As the expansion velocity increases, the polarization begins to develop near the central frequencies and r ≈ 0. For instance, for the highest velocity (V max = 2000 km s −1 ) model, the frequency at which the maximum P L is found to change linearly with the radius r; the peak frequency shifts to red as r increases. This is because photons are more easily Doppler-shifted at smaller radii when the velocity gradient is larger. It is also noticeable that the peak polarization, denoted by the reddest color in the figure, occurs at a slightly lower frequency than the peak frequency, indicated by the thick black lines, especially in the highest velocity models. The variation in the polarization properties with the expansion velocity illustrates that the spatially-resolved spectropolarimetric observations would be of great importance to understand the kinematics of hydrogen gas in the galactic halos.

An Outflowing Thin Shell
The galactic superwind model of Dijkstra & Loeb (2008) was developed to explain the Lyα spectra redshifted relative to the systemic velocity, measured using non-resonant lines, of Lyman break galaxies (LBGs) and LAEs. (Ahn et al. 2003;Verhamme et al. 2006Verhamme et al. , 2008. In the model, the redshifted spectra are attributed to the Doppler boost that Lyα photons undergo when they are scattered back toward the observer by an outflowing shell on the far side of the galaxy. The outflow is modeled by a radially expanding thin shell. The shell expands at a constant speed of V exp = 200 km s −1 and has a column density of N HI = 10 19 or 10 20 cm −2 . The initial Lyα spectrum originating from a central source is assumed to be a Gaussian with a velocity width σ = 200 km s −1 , which corresponds to a typical circular velocity of a dark matter halo with a total mass of M tot = 3 × 10 11 M . This example was adopted to verify our code and investigate a case in which single scattering is predominant. Figure 18 shows the resulting emergent spectrum, surface brightness, and degree of linear polarization, consistent with those of Dijkstra & Loeb (2008). 10 The left panel of the figure shows that the spectrum becomes more redshifted and broader 10 We found that a variation of the shell thickness gives rise to small yet recognizable changes. As the thickness decreases, the emergent spectrum, particularly for N HI = 10 20 cm −2 , was found to be less redshifted, the surface brightnesses became flatter, and the maximum degree of polarization as N HI increases. The green line in the same panel denotes the initial input spectrum. The middle panel shows the profile of surface brightness as a function of projected radial distance r p , which is relatively flat up to r p ∼ 0.7r max . We note that a considerable fraction of Lyα flux (49% for N HI = 10 19 cm −2 and 8% for N HI = 10 20 cm −2 ) escapes out of the system without experiencing any scattering. As a result, a sharp central peak is apparent in the surface brightness profile, but it is, unfortunately, not consistent with the shapes observed in LAEs. Note that, in this model, the ISM can be regarded to be concentrated in a small central region together with the Lyα emission source, and the galactic halo is assumed to be geometrically thin and far from the ISM. The presence of a peak in the surface brightness indicates that the galactic wind should be dealt with a spatially extended medium rather than a geometrically thin shell so that the ISM and galactic halo are seamlessly connected without an extreme gap. The right panel shows that the degree of linear polarization tends to increase with radius, as expected. The degree of polarization at the center is zero due to symmetry. In the outermost part, the polarization reaches as high as ∼ 50% for N HI = 10 19 cm −2 and ∼ 20% for N HI = 10 20 cm −2 . Dijkstra & Loeb (2008) found that the probability distribution of the number of scattering events (N scatt ) peaks at N scatt = 1 and rapidly decreases with increasing N scatt . This is because, in this model, most photons reside on the wing of the line profile when they reach the expanding shell; thus, they escape the system without any scattering or after only one scattering. The left panel of Figure 19 shows the average ( N scatt ) and most frequent value (N mode scatt ) of the number of scatterings as a function of projected radius in the detector plane. In the middle panel, we also show the probability distribution function of the number of scatterings, P(N scatt ). The distribution function peaks at N scatt = 1 for N HI = 10 19 cm −2 and N scatt = 1 − 3 for N HI = 10 20 cm −2 , and decreases sharply with increasing N scatt . The median number of scatterings is found to be N median scatt = 3 for N HI = 10 19 cm −2 and N median scatt = 9 for N HI = 10 20 cm −2 , considering only the photons that have experienced scatterings at least once (but not shown in the figure). Even though most photons are scattered only a few times at most, the distribution function of N scatt shows a very long tail. Thus, N scatt is found to be very high. Interestingly, N scatt remains nearly constant as r p increases, except for a slight rise near r p = r max due to the rim effect. The highest optical thickness arises toward the line of sight of r p ≈ r max . The dip at r p ≈ 0.9r max in the model of N HI = 10 19 cm −2 is also a geometrical effect due to a finite shell thickness of 0.1r max , adopted in our calculation. In the Rayleigh regime, scatterings toward the observer at the right angle at r p = 0.9r max are least probable, causing the dip. The dip disappears for a higher N HI because of an increase in the number of scatterings. Dijkstra & Loeb (2008) also noticed that the "ensemble-averaged" polarization level of photons that experienced more than 10 scatterings (N scatt > 10) nearly vanishes. They also found that a group of photons scattered only once gives rise to the highest (ensemble-averaged) polarization, while photons scattered multiple times yield much lower polarization. To confirm the results, we computed decreased slightly. A shell thickness of ∼ 0.1rmax was found to yield surface brightness profiles that best match Dijkstra & Loeb (2008); hence, a thickness of 0.1rmax was adopted in this paper.  Dijkstra & Loeb (2008). The emergent spectrum (Jx), surface brightness (S), and degree of linear polarization (P L ) are shown in the left, middle, and right panels, respectively. In the middle panel, the intensity ratio of the direct (I direc ) to the scattered (Iscatt) intensity at the center (r = 0) is also denoted. The blue and red represent the models with column densities of N HI = 10 19 and 10 20 cm −2 , respectively. Models without and with dust are denoted by lines and circles, respectively. Here, rmax and L are the radial size and Lyα luminosity of the system, respectively. The models containing dust grains have the dust optical depths of τ dust = 1.6 × 10 −2 (N HI = 10 19 cm −2 ) and 1.6 × 10 −1 (N HI = 10 20 cm −2 ).  Figure 19. Results of an expanding thin shell model of Dijkstra & Loeb (2008). The left panel shows the mean ( Nscatt ) and mode (N mode scatt ) of the number of scatterings, which a photon experiences before escaping the system, as a function of projected radius (or impact parameter). The middle panel represents the probability distribution function of the number of scatterings. The right panel shows the degree of linear polarization as a function of the number of scatterings that photons experience. Here, P indiv L denotes the mean value of the "individual" linear polarization amplitudes of photons that are scattered Nscatt times, while P avg L is the "ensemble-averaged" polarization level. The blue and red represent the models with column densities of N HI = 10 19 and 10 20 cm −2 , respectively. the "ensemble-averaged" degree of polarization for photon packets escaping after a specific number of scatterings, as shown in the right panel of Figure 19. For comparison, we also present the polarization degree of "individual" photons. In the figure, P indiv L is the mean polarization amplitude of individual photon packets, measured regardless of the polarization direction, and P avg L the ensemble-averaged polarization, i.e., the polarization calculated after averaging the Stokes Q and U of photon packets that were scattered N scatt times. As in Dijkstra & Loeb (2008), P avg L peaks at N scatt = 1 and decreases as N scatt increases. On the other hand, P indiv L increases with N scatt in a range of N scatt 10, and then decreases; but, as N scatt further increases (N scatt 10 2 ), it increases again until it finally saturates. The saturated level of the polarization of individual photons reaches on average as high as ∼ 85% for N HI = 10 19 cm −2 and ∼ 100% for N HI = 10 20 cm −2 . Figure 20 shows the distribution functions of the polarization angle χ * , the Stokes Q and U, and the linear polarization degree P L of individual photons at different radii in the detector plane. In this figure, only photon packets that experienced scatterings at least once or more were considered. The polarization angle χ * is uniformly distributed at the central part (r p ∼ 0), except for a very sharp peak at χ * = 0.
The peaks at χ * = 0 found in all radius bins are mainly due to photons escaping after only a single or couple of scatterings. The peaks of the distribution functions of χ * , Q, and U found in this model are much stronger and sharper than those shown in other models. This is because the thin shell model is predominated by single-scattering, while the other models are by multiple scattering. The polarization degree of individual photons tends to have P L < 0.5 for r p 0.7r max but increases with increasing r p . At r p ∼ r max , most photons have a polarization of P L ∼ 1. It also appears that, as N HI increases, more photons tend to have the polarization vectors aligned to the tangential direction (χ * = 0) and be 100% polarized.
The distribution functions of Q are similar to those in the static sphere models with τ 0 10 5 , except for two aspects: the presence of a rectangle shape bump and the weakness (or absence) of the peak near Q = 0 in radii r p 0.7r max . The rectangle bump, which is caused by photons escaping after a single scattering (N scatt = 1), is also found in the distribution function of P L . When we ignored the single-scattered photons, the bumps disappeared, leading to smoothly varying distribution functions. The bump tends to move toward higher P L and |Q|, as r p increases. This trend can be understood with the help of the Rayleigh scattering phase function in Equation (22). Suppose photons that escape . Probability density functions of P L , χ * , and the Stokes Q and U of "individual" photons for the expanding thin shell model of Dijkstra & Loeb (2008). They are all measured in a local coordinate system defined by radial and tangential vectors in the detector plane. Note that the ordinates were plotted in logarithmic scale to highlight weak changes in their values. Only the photon packets that have experienced scatterings at least once (Nscatt ≥ 1) are included in calculating the probability distribution functions.  Dijkstra & Loeb (2008). The first and second columns show the results for the models without dust, while the third column shows the results for the model with dust. In the small ancillary panels, the black dotted and red solid curves are the initial and emergent spectra, respectively. The green curves represent polarization degrees as a function of frequency. The main panels show the linear polarization degree in a two-dimensional space of photon frequency and projected radius. The color bar denotes the polarization levels.
after a single Rayleigh scattering. As the projected radius at which a photon is scattered increases from r p = 0 to r max , the scattering angle toward the observer also increases from θ = 0 to π/2. The increase in θ leads to higher P L and |Q| because the right angle scattering gives rise to the highest polarization level. The extent of Q in the rectangular bump corresponds to the range of radii from which the photons escape. In Equation (21), we also note that Q = S 12 I (for Q = U = 0) is negative. This is why the rectangular shapes occur only at negative Q values. The shapes of the distribution functions of Q at large radii resemble those for the Hubble-like expanding media that show no peak at Q = 0, while the sharp peaks at small radii are more like those seen in the static models.
We also note that the polarization degree of individual photon packets shows a rather complex variation of increasing and decreasing, as shown in the right panel of Figure 19.
However, eventually, photons that experienced a sufficiently large number of scatterings end up having a very high polarization level of up to ∼100%. As noted in Section 4.1, if once photon packets were diffused out into the wings, further scatterings lead to an increase in the degree of polarization of individual photon packets. However, photons that have traversed through many different routes will have randomly oriented polarization directions. The randomness of the polarization direction will largely cancel out the Stokes Q and U parameters. In the end, the final "ensemble-averaged" polarization level will be much lowered than those of individual photons. The model with N HI = 10 20 cm −2 produces higher polarizations for individual photon packets than the model with N HI = 10 19 cm −2 (Figures 18  and 20) but gives a lower "ensemble-averaged" polarization because of higher isotropy (P avg L in Figure 18). . Spectrum and polarization degree as a function of wavelength (first column), surface brightness profile (second column), polarization vs. projected radius (third column), and polarization degree as a function of wavelength and projected radius (fourth column) of four LAEs (MUSE 1185, 0082, 6905, and1343), which are obtained using a momentum-driven galactic wind model of Song et al. (2020). In the first and second columns, the filled circles with error bars are the observational data. The blue/red lines in the first, second, and third panels are the best-fit models before/after convolution with the instrumental line or point spead function. The green lines in the first column represent the polarization degree vs. wavelength, convolved with the instrumental line spread function.
The green dashed lines in the second column is the point spread function of MUSE. The spectra are shown as a function of the observational wavelength as well as the rest frame wavelength. The third column show the degree of linear polarization, predicted using the best-fit model. The theoretical spectra and surface brightness/polarization profiles were obtained after applying the aperture sizes and wavelength ranges, respectively, described in Leclercq et al. (2017) and Song et al. (2020). In the third column, the polarization profiles constructed using all photons are also shown in dotted lines; the blue and red lines denote the profiles before and after the convolution with the PSF, respectively. The best-fit model parameters are also shown in the panels of the third column. Here, r UV and r HI are the scale parameters of the distribution functions of Lyα source and gas density. DGR is the dust-to-gas ratio relative to that of MW, τ 0 the optical depth, and z the redshift of galaxy. r peak , V peak , and ∆V define the velocity profile of the galactic wind, as in Equation (50); V peak is the peak velocity, r peak is the radius at which the fluid velocity peaks, and ∆V is the difference between V peak and the velocity at r = rmax. All size variables (r UV ,r HI , and r peak ) are normalized to the size of halo (rmax). The velocity is in units of km s −1 . The rightmost panels are the color bars to represent the polarization level. We additionally computed the thin shell models that include dust grains; the filled circles in Figure 18 show the results. The dust-to-gas ratio was set to that of the MW. In the figure, especially for the higher N HI , the presence of dust appears to make the surface brightness profile slightly steeper and raise the polarization level; no significant enhancement in the polarization is found for N HI = 10 19 cm −2 , as in the static model of τ 0 = 10 6 because of relative weakness of the dust extinction effect. We also found that ignoring single-scattered photons significantly reduced the ensemble-averaged degree of polarization (but not shown in figures); this is because the radiation field's anisotropy in this model is mainly attributed to single scattered photons, and ignoring them decreases the anisotropy of the radiation field.
In Figure 21, we show the dependence of polarization degree on radius and frequency. As in the Hubble-like outflow model, the polarization level is found to increase toward the red wavelength. This is because the redder photons are farther from resonance in the frame of the expanding gas; therefore, they are scattered less, yield a relatively steeper radial profile compared to the bluer photons, and eventually achieve a higher polarization. However, the result appears more complicated if we examine the variation of the frequency dependence of polarization with radius. At small radii, red photons tend to give a higher polarization. In contrast, near r p = r max , higher polarization is found at blue wavelengths. This is because blue photons are scattered more (10 times more than red photons in the radius range r p > 0.9r max ). The radiation field at r ≈ r max is highly anisotropic, independently of the wavelength, because of the geometrical effect. In the boundary region, the "ensemble-averaged" polarization level appears to be primarily determined by the polarization of "individual" photons rather than by the anisotropy of the radiation field. However, despite this opposite trend at large radii, the spatially averaged polarization is higher in red wavelengths because fluxes at smaller radii dominate the total flux. In the right panel, we note that the model with dust reveals a lower maximum polarization than the case without dust (middle panel). This is because of the preferential extinction of blue photons scattered more. However, at small radii, the model with dust gives a higher polarization than the model without dust; therefore, dust grains yield, in the end, a higher polarization in the spatially averaged spectrum than the dustless case. Song et al. (2020) recently simultaneously modeled the emergent spectra and surface brightness profiles of eight LAEs at z = 3 − 6 observed with the Multi-Unit Spectroscopic Explorer (Bacon et al. 2015;Leclercq et al. 2020). In the study, we considered a galaxy as an outflowing halo in which the distributions of Lyα source and gas are described by exponential functions. The outflow of the medium is described by a simplified momentum-driven galactic wind; rather than adopting the functional form in Dijkstra & Kramer (2012) 11 , obtained by solving a momentum equation in a gravitational potential well, we instead considered a linearly increasing and then decreasing velocity profile, as follows:

A Momentum-Driven Galactic Wind Model
V (r) = V peak r/r peak if r ≤ r peak V peak + ∆V (r − r peak )/(r max − r peak ) otherwise, (50) where V peak is the peak velocity, r peak the radius where the velocity peaks and ∆V the difference between velocities at r peak and r max . In the model, we assumed ∆V ≤ 0 to consider only the case in which the gas accelerates and then decelerates.
In the model calculation, Song et al. (2020) adopted uniform distribution functions for the initial input spectrum and spatial distribution of the Lyα source. Then, to reduce the number of simulations, they obtained the final results by multiplying the intermediate results by the desired input profiles (i.e., Voigt and Bessel profiles for the spectrum and UV surface brightness profile, respectively). The resulting spectrum and surface brightness profile were convolved by the spectral and spatial kernels appropriate to the instrument, and were compared with the observational data to find the best-fit parameters. For this paper, however, we directly generated the desired input spectrum and spatial distribution of injected photons using the best-fit parameters instead of adopting a uniform distribution and then adjusting the inputs by multiplying appropriate distribution functions. The present study also took into account the polarization state of Lyα in 11 Adopting the analytical solution of Dijkstra & Kramer (2012) yielded unrealistic spectra with a strong peak at the Lyα line center (x = 0) due to an extremely large velocity gradient at r = 0. This was why we adopted a linearly increasing and then decreasing velocity profile instead of their solution. The first and second columns in Figures 22 and 23 compare the observational spectra and surface brightness profiles with those obtained by simulations for the LAEs. The output spectra and surface brightness profiles calculated in this study are confirmed to reproduce those of Song et al. (2020) well. We also present the polarization profiles predicted from the best-fit models for each LAE in the third and fourth columns of the figures; the third column shows the radial profile of polarization degree, and the fourth column the degree of polarization as a function of wavelength and projected radius. To compare with the observational data, the calculated spectra and surface brightness/polarization profiles were obtained after applying the aperture sizes and wavelength ranges described in Leclercq et al. (2017) and Song et al. (2020). The radial profiles were also convolved with the MUSE point spread function (PSF). In the third column, the polarization profiles obtained using all photons, including those having wavelengths outside the wavelength boundaries, are also shown in dotted lines.
In the figures, it is noteworthy that the polarization profile is fairly diverse in shape. In previous models, the polarization level monotonically increased with radial distance, even though the rate of increase (or slope) was dependent on model parameters. However, in the models of MUSE 6905, 0547, and 0364, the radial profile of polarization is found not to be a simple monotonic function. In particular, the model of MUSE 6905 shows an interesting feature, even after the convolution with the instrumental PSF. The polarization profiles for MUSE 0082, 0171, 0547, and 0364 show a rapid jump near r ≈ 0, although the convolution by the instrumental PSF smoothes out the profiles, erasing the jump shapes. The third column also demonstrates that the polarization profile can be significantly altered, at least in some cases (MUSE 6905 and 0364), depending on the wavelength range adopted to measure the polarization signal.
The polarization property of individual photons is too complex to interpret, unlike the other models. Detailed analysis on the variation of polarization pattern depending on the model parameters is beyond of the scope of this paper. However, it is clear that the polarization profile depends on the kinematic properties of the galactic wind, as shown in Figure  24. The distribution functions of P L and χ * for MUSE 1185 are somewhat similar to those obtained for the Hubble-like flow model. MUSE 0082 shows several spiky peaks, whose cause is not clear, in the distribution function of P L . The distribution functions of χ * for MUSE 1185 and 0082 indicate that the radiation field becomes more anisotropic as the radius increases, yieding the rise of polarization.
The distribution function of P L for MUSE 6905 reveals a more or less similar trend to those for the static models of T = 10 4 K and τ 0 10 5 (Figure 11), except that the static models barely show the dependence on radius. However, unlike the static model, the degree of isotropy of the radiation field shows a non-monotonic trend, decreasing and then later increasing with increasing radius, as shown in the distribution of χ * . This trend in the model of MUSE 6905 is primarily due to the velocity profile of outflow; the outflow monotonically decelerates, starting with the maximum velocity of V exp = 300 km s −1 at the center, as the radius increases. However, except for MUSE 0171, the other outflow models accelerate, starting at rest, until a radius (r peak > 0) and then decelerate. The difference between the models for MUSE 0171 and 6905 is that the outflow of MUSE 0171 expands at a constant speed, while that of MUSE 6905 monotonically decelerates with increasing radius. The peculiarity in the polarization pattern of MUSE 6905 is attributable to this monotonically decelerating velocity profile. Lyα photons are recognized as wing photons at small radii because of the high velocity there. However, the wing photons that gained high polarization at small radii will be identified to be in the core regime when they reach large radii. Then, their polarization levels will be reduced by the core scattering at large radii, as can be seen in the distribution function of P L . Furthermore, the radiation field becomes highly isotropic due to many scatterings, as shown in the distribution function of χ * . Eventually, the ensemble-averaged polarization level decreases at large radii. Note that the model of Song et al. (2020) is a simplified version to avoid an exceedingly abrupt jump in outflow speed. Therefore, the model of MUSE 6905 may be applicable to the systems where a rapid outflow is launched near the galactic center and then decelerate monotonically.

DISCUSSION
The energy levels of hydrogen are splitted not only by the coupling between the electon wavefunctions but also by the interference between the electron and proton wavefunctions. In the present paper, we considered only the fine structure of hydrogen atom. The hyperfine structure was taken into account in Chamberlain (1990), Brasken &Kyrola (1998), andHirata (2006). It is relatively straightforward to extend the present RT algorithm to include the hyperfine structure.
However, the present scope is sufficient for the polarization RT in most cases interested in astrophysics.
In our study, we used the quasi-monochromatic photon packets that represent superpositions of mutually independent, incoherent waves. We found that Lyα photon packets tend to be fully polarized, unlike the dust scattering case, when repeatedly scattered. This result indicates that the quasi-monochromatic Lyα light that has undergone the same scattering history will be fully polarized as if they are purely monochromatic waves. In the following, we further discuss the tendency of photon packets to get 100% polarized, together with other topics relevant to Lyα but not addressed or fully discussed in the previous sections. The topics include the distribution function of the Stokes parameter Q, an implication of the correlation between the surface brightness profile and the polarization profile to the observations of the polarization of LABs, the negative polarization, and a sudden jump of polarization. We also discuss possibilities to break the degeneracies between different models of LAEs using the polarization data, and then a misleading concept relevant to the peeling-off technique.

Arise of 100% Polarization of Lyα by Scatterings
In the above sections, we noticed that "individual" Lyα photon packets gradually approach an almost 100% polarized state when they are repeatedly scattered. Here, we discuss in more detail how this process can occur. We first need to note that the degree of polarization is preserved under a rotation of the axes about the photon's direction vector. Therefore, it is convenient to define the Stokes vector in a reference frame where the polarization basis vector m lies in the scattering plane (i.e., φ = 0). In that coordinate system, the Stokes parameters of scattered light are given by: Now, we can readily show that Here, I 2 − Q 2 − U 2 denotes the intensity of unpolarized light. The above equation indicates that isotropic scattering (E 1 = 0) always completely depolarizes photons (Q = U = 0) regardless of their initial polarization states. 12 It can also be proven that if the incident light was 100% polarized, the scattered light is 100% polarized if and only if E 1 = 1 (Rayleigh scattering). We note that Lyα photons diffuse out in frequency space while being scattered, approaching the wing part of the line profile. If once the photon reaches the wing, where E 1 = 1, the above equation is reduced to the following: Probability distribution function of Q/I for the static, homogeneous sphere model. The gas temperature is T = 10 K in the left panel and T = 10 4 K in the right panel. The optical depth for both panels is τ 0 = 10 7 . The red and blue curves are the distribution functions obtained in the radial bins 0.1 < r < 0.1 and 0.9 < r < 1.0, respectively. The black curve denotes the theoretical distribution function of Equation (55).
Note that the factor (9/4 cos 2 θ) in the right-hand side results in 0.9 when multiplied by the scattering phase function for Rayleigh scattering and averaged over the scattering angle. This indicates that the intensity of unpolarized light gradually decreases; in other words, the degree of polarization will gradually rise while photon packets are being Rayleigh scattered. Moreover, we can immediately recognize that if a photon packet in the wing regime is scattered at the "right" angle (cos θ = 0), the photon will be 100% polarized regardless of its initial polarization state. After then, the photon will keep the polarization state of ∼100%, independently of the scattering angle, unless it happens to be scattered back into the core regime. In this stage, the photon is highly likely to escape the system in a single long-excursion while trapped in the wing (Adams 1972); consequently, the individual photon packets escaping in the wings will have a very high polarization degree of ∼100%. In a static medium, photons need a lot of scatterings to diffuse into the wings. On the other hand, in a fast-moving medium, photons will be scattered into the wings relatively easily with only a small number of scatterings. Therefore, Lyα photon packets will approach a 100%-polarized state in an outflowing or infalling medium quickly, unless the gas motion is too fast.
It is also expected that individual photon packets of broad Lyα lines, as observed in AGNs, could readily be fully polarized because most of the line photons will be wing-scattered even in static media. At the same time, the depolarization effect due to the core scattering may play essential roles in the broad lines if the velocity gradient is considerable in the outflow or inflowing medium.

Distribution Function of the Stokes Parameters Q
It is evident that, at the center (r = 0) of a sphere, the polarization angle of individual photons will be uniformly distributed because of the isotropy of the radiation field. Our results also indicate that the distribution function of the Stokes parameter Q, measured at the central region, is symmetric. As the optical depth and number of scatterings increase, the polarization of individual photon packets approaches 100%.
At the same time, the radiation field tends to be isotropic, particularly even up to relatively large radii in static media. In this case, it can be found that a simple analytic function describes the distribution function of Q.
To derive the probability distribution function for Q, we assume an isotropic radiation field and that the polarization angle χ is uniformly distributed in a range of 0 ≤ χ ≤ 2π, as given by: From Equation (12), we obtain the distribution function of the normalized Stokes parameter Q/I, as follows: Here, the coefficient "4" is adopted to take into account the fact that Q/I has the same value four times in the range 0 ≤ χ ≤ 2π. The obtained distribution function of Q has two peaks at Q = −1 (concentric polarization) and Q = +1 (radial polarization). It is evidently far from a Gaussian or Poisson function that is generally adopted to analyze the intensity I. We also note that the Stokes parameter U/I follows the same distribution function. As we already noted, the individual photon packets tend to achieve ∼ 100% polarization and thus the distribution function of Q (and U) is well approximated by the distribution function of Q/I (and U/I). In other words, P(Q) ≈ P(Q/I). Figure 25 shows that the theoretical distribution function derived above indeed well reproduces those obtained in the static sphere model, especially in the central region. We note that the distribution function of Q appears to be reasonably well described, even near the boundary (r ≈ 1), by the same function. This probability density function can be a starting point to understand the statistical properties of the Stokes parameters Q and U in actual observations. We note that similar shapes with two peaks at Q = ±1 are also commonly found in the distribution functions of Q of Figures 15 and 16. The relative strengths of the two peaks are closely associated with the degree of the isotropy of the radiation field. The peak at Q = 1 weakens while that of Q = −1 becomes stronger as the gas motion boosts up the anisotropy of the radiation field. Related to this point, it is noteworthy that the polarization profiles observed in LABs rise somewhat slowly with radial distance from the center You et al. 2017;Kim et al. 2020). This trend suggests that the Lyα radiation field in the LABs might be more or less isotropic, and the distribution function of Q may have a relatively strong peak at Q = 1. Such a peak in P(Q) at Q = 1 would make the observed polarization signals noisy.

An Implication to the Polarization Pattern in LAB1
In our results, we noticed that there exists a correlation between the surface brightness profile and the polarization profile. In a static spherical model, as the optical depth increases, the overall shape of the surface brightness profile flattens out and, at the same time, the polarization profile becomes shallower near r = 0. This result is associated with the tendency that as the optical depth increases, the Lyα radiation field becomes more isotropic and the polarization vectors of near-100%-polarized individual photons largely cancels out. A similar correlation was also found in a Hubble-like expanding modium. In this case, however, the polarization levels of individual photons are primarily boosted up by the Doppler shift and a faster medium gives rise to a steeply declining surface brightness profile and a rapidly rising polarization pattern.
In the well-studied Lyα blob LAB1, the degree of polarization is observed to increase relatively slowly with radius. Trebitsch et al. (2016) utilized a hydrodynamic simulation in the context of the gravitational cooling and presented radial polarization patterns predicted by assuming central sources and diffuse ionized gas surrounding the sources. The polarization pattern for the central sources was found to be relatively steep, which is more or less similar to those obtained in fastly expanding media discussed in this paper. On the other hand, the polarization obtained using the in-situ diffuse emission gas appeared to show a slowly varying pattern. Based on these results, they showed that the observed polarization in LAB1 is consistent with the pattern predicted by the circumgalactic diffuse gas or by the combination of two source types.
However, it is worth noting that the observed surface brightness profile in Hayes et al. (2011) is shallower than that predicted by Trebitsch et al. (2016). Shallow surface brightness profiles are expected for fast-moving gases or media with small optical depths. In that sense, their hydrodynamic simulation seems to predict a somewhat too fast gas motion or a too-small optical depth. Recalling that a slow rise in polarization usually accompanies a shallower surface brightness profile, the observed polarization and surface brightness patterns accord well, at least qualitatively, with the central engine scenario in which Lyα photons originating from central sources are scattered by the hydrogen gas surrounding the sources yielding the Lyα blob. Quantitative comparisons of the observational data with more extensive theoretical calculations, taking into account inhomogeneous and clumpy media, would be required to understand the origin of these polarization signals.

The Negative Polarization and Polarization Jump
LaRT is superb, compared to the preexisting codes, in that it uses a smoothly and seamlessly varying phase function as frequency changes. Our approach provides a theoretically complete framework in dealing with the polarization signals of resonance doublet lines with the same fine structure as Lyα. The other approaches that use a discrete phase function or 100%-polarized photons cannot reproduce the "negative" polarization signals perpendicular to the radial direction. Regarding the Lyα RT, however, our approach appears to yield no significant difference from the preexisting approaches; this is because the number of scatterings of Lyα is enormous, and the frequency gap between the fine structure is very narrow so that the chance of being negatively polarized is rare, as described in Section 2.5. The present method would be advantageous in studying the resonance lines that have the same fine structure as Lyα but relatively low optical depths, such as C IV λλ1548, 1550, Mg II λλ2795, 2803, and Ca II λλ3934, 3968 doublets. The negative polarization was indeed observed in the Ca II doublet line in the Sun (e.g., Stenflo 1980).
In Sections 4.1 and 4.2, we discovered an abrupt jump in the polarization level at the center of the spherical models. However, it should be noted that the polarization jump may not be readily observable in actual observations. No galaxies are perfectly symmetric, and currently available telescopes do not seem to have sufficient spatial resolutions to resolve the abrupt change in the polarization pattern. The effect of the telescope PSF is demonstrated for the models of LAEs in Section 4.4. Figures 22 and 23 show a rapid jump near r ≈ 0 in the polarization profiles of MUSE 0082, 0171, 0547, and 0364. However, these jump features are found to disappear after the convolution with the telescope PSF mostly. Asymmetries in the Lyα source and medium may also blur the central polarization profile and thus significantly reduce the jump level. Clumpiness in the medium might also tend to smooth out the abrupt change in the polarization pattern. Song et al. (2020) demonstrated the existence of degeneracies between different models of LAEs when analyzing only either the spectrum or the surface brightness profile, and at least some of them could be broken by utilizing both the spectrum and surface brightness profile simultaneously. For instance, analyzing only the surface brightness profile was found to show degeneracies of ∆V − rs HI and τ 0 − V peak for most galaxies. The output spectrum was revealed to be relatively independent of the scale length of the gas density to a large extent. Analysis of the polarization signal, together with the spectrum and surface brightness profile, may provide additional clues, which cannot be revealed by analyzing only spectra and surface brightness profiles, and help us gain a deeper understanding of the Lyα galaxies and blobs (Eide et al. 2018). As discussed in Section 4.2, the variation in the polarization properties with the outflow speed demonstrates the importance of the spatially-resolved spectropolarimetric observations in understanding the kinematics of the galactic halos.

Breaking the Degeneracies between Models of LAEs
The galactic halo model in Section 4.4 illustrates the diversity of the polarization patterns. Therefore, a more detailed understanding of the polarization patterns may help us better constrain the galactic halos' kinematic properties. For this purpose, it would be worth investigating how the polarization pattern changes by varying each model parameter while keeping other parameters fixed, as done for the spectrum and surface brightness profile in Song et al. (2020). However, we would like to defer such a detailed analysis to a future study.

A Note on the Peeling-Off Technique
We here would like to clarify a misconception with regards to the peeling-off technique. Murthy (2016) and Dijkstra (2019) claim that one has to subtract the peeled-off fraction from the photon's weight whenever a peeling-off procedure is performed. Dijkstra (2019) also states that the portion to be subtracted is too small to affect the simulation result, and thus one can, in practice, safely ignore the subtraction. However, in principle, the peeled-off fraction must not be subtracted from the photon's weight.
The basic idea behind this technique is to calculate the number (or probability) of experiments (realizations) in which the photon is scattered towards the observer if a sufficiently large number of identical experiments have been performed. Suppose that we trace the trajectory of a single photon and perform a countless number of Monte Carlo experiments. We construct an ensemble made up of an almost infinite number of Monte Carlo experiments performed under identical conditions and record the whole history of every (emission and scattering) event in each experiment. Next, we choose a subset of the ensemble, which is composed of experiments where the photon's history is exactly the same until a specific time t 0 . Performing a scattering of the photon toward a randomly chosen direction in the next simulation step (t 0 + δt) is equivalent to picking a particular experiment (a "realization") among the subset by drawing a random number. Similarly, performing a peeling-off procedure is simply to count the number of experiments in which the photon is scattered toward the detector. Therefore, the peeling-off procedure is nothing to do with performing the next scattering of the photon by drawing a random scattering direction for a specific simulation realization. Thus the photon weight should not be altered after the peeling-off procedure.
In this sense, the term "peeling-off" may be misleading. The "next event estimation" would be a better term, although it is unfamiliar in the astrophysics community. We also note that if one performs the peeling-off technique toward a detector plane, they should not count the photons that happen to be actually detected in the detector plane by chance. However, fortunately, even if one adopts the wrong method, it is unlikely to affect the result in most cases significantly because the peeled-off fraction would be tiny.

SUMMARY
This paper describes an RT method using the Monte Carlo technique to investigate the Lyα polarization. Our results demonstrate that the spectropolarimetry data can be used to supplement the spectrum and surface brightness profile of the LAEs and LABs to tightly constrain the physical origin of Lyα halos. Furthermore, spatially-resolved spectropolarimetry will provide a deep understanding of the detailed nature of Lyα emitting sources. For this purpose, we need to not only take observations for many different objects but also develop theoretical models for various possible situations. This paper provides mainly a primary tool and fundamental results from the theoretical point of view. We expect more detailed models will be developed in future works.
The principal conclusions of this paper are summarized as follows: • The degrees of polarization of individual photon packets increase while they are scattered repeatedly, approaching the wing regime. If once the photon packets diffuse out to red or blue wavelengths in the Rayleigh scattering regime, they can easily achieve 100% polarization and their polarization level is persistent until they escape, unless they are scattered into the core of the line profile.
• In outflowing or contracting media, the polarization degree of individual photon packets rises fairly quickly by the Doppler shift effect.
• The "ensemble-averaged" degree of polarization, the observable quantity, is primarily determined by (1) the polarization levels of individual photon packets and (2) the degree of the isotropy (or anisotropy) of the Lyα radiation field. They are both eventually controlled by the number of scatterings and the Doppler boosting due to the gas motion. In most cases, the polarization degree of individual photons approaches 100%, unless the medium's optical depth is too low.
• There exists a correlation between the surface brightness profile and the polarization profile. This is because a shallow surface brightness profile is associated with a relatively isotropic radiation field, while a steeply declining profile is due to an anisotropic radiation field. At the center of spherical media, the polarization vectors are canceled out due to the symmetry. If the surface brightness profile is steep, the polarization immediately survives as soon as the isotropy of the radiation field is broken at a location in the vicinity of the source, yielding a rapidly rising polarization profile. On the other hand, for a shallow surface brightness profile, the radiation field is somewhat isotropic even at a location far from the source, so the polarization vector is significantly canceled there, leading to a slowly varying polarization pattern.
• The correlation between the radial profiles of the surface brightness and polarization can help interpret the observed data of LABs. The relatively shallow profiles of the surface brightness and polarization observed in the Lyα blob LAB1 seem to support the central engine scenario of the Lyα halo pretty well. However, further quantitative studies would be required to understand better the origin of the observed polarization patterns of LABs.
• The polarization angle of individual photon packets can be regarded as a proxy to measure the degree of isotropy of the radiation field. Assuming a uniform distribution for the polarization angle, we can derive a distribution function for the Stokes parameter Q, which may be useful in understanding the statistics of observational data.
• Dust grains slightly steepen the surface brightness and increase the degree of polarization by selectively destroying the core photons that hold relatively low amplitudes of polarization. The effect appears more prominent in a medium of a higher temperature.
• We also found that the polarization profile can exhibit a non-monotonically increasing pattern with radius in some models developed for LAEs in Song et al. (2020). The unusual profiles are primarily due to the depolarization by the core scattering in a decelerating medium. Hence, the spatially-resolved polarization observations may play a critical role in understanding the galactic outflow or inflow by providing complementary information on the kinematics of the galactic halo gas, in addition to the spectrum and surface brightness profile. APPENDIX A. SCATTERING MATRIX FOR LYα LINE SCATTERING Stenflo (1980) and Stenflo (1994) provide the formulae for the scattering matrix for the Lyα scattering by hydrogen atoms, but without the term for circular polarization. We, therefore, rederived the formulae, including the circular polarization term. The time reversal term of the Quantum Figure A1. Coordinate system and polarization basis vector to calculate the scattering matrix of Lyα. The spherical unit vectors (ē 0 ,ē ± ) are defined in a fixed reference frame (determined by the scattering plane and the initial propagation vector) and the polarization basis vectors (e 1 , e 2 ) in the photon's local frame. The propagation vectors of the incident and scattered photon are denoted by k and k, respectively. The primed variable refers to the incident photon but the non-primed the scattered photon. See Figure 3.3 of Stenflo (1994) for a more general definition. For simplicity, we choose the angles defined in Stenflo (1994) to be θ S94 = θ S94 = 90 • , φ S94 = 0, and φ S94 = θ.
Field theory, which is negligible, will be ignored in the Kramers-Heisenberg dispersion formula. Stenflo (1994) adopts the convention of E 1 = E n and E 2 = −E m and defines the Stokes parameters as in Equation (11), but replacing (E 1 , E 2 ) by (E m , E n ). They define the scattering plane as the y − z plane (in the photon's local frame), as shown in Figure A1, while we define it as the x − z plane. They define the handedness of circular polarization in the opposite direction. Therefore, their Stokes parameters are expressed as (I S94 , Q S94 , U S94 , V S94 ) = (I IAU , −Q IAU , −U IAU , −V IAU ). The scattering matrix elements are then the same as those obtained for the IAU standard except for only one; M IAU 12 = −M S94 12 for (i, j) = (1, 2) and i j otherwise. We follow the definition of Stenflo (1994) to avoid confusion when deriving the scattering matrix elements and then later transform the obtained results to the IAU standard.
The 4-dimensional coherency vector is defined by where the superscript "S94" indicates the definition of Stenflo (1994). The Stokes vector is related to the coherency vector by S = TD (4) , where T is the matrix to transform the 4-dimensional coherency vector to the Stokes vector. The scattering matrix for the 4-dimensional coherency vector D (4) is given by Equation (2.39) and (9.23) in Stenflo (1994), as follows: where other matrix elements vanish because µ i = µ f and µ i = µ f ± 1 cannot be simultaneously satisfied. The scattering amplitude is given by Equation (8.115) of Stenflo (1994) (see also Stenflo (1998)): Here, the 2 × 3 matrices are the Wigner 3-j symbols, J and µ indicate the total angular momentum quantum number and magnetic quantum number, i and f denote the initial and final levels (both are lower levels), and m represents the intermediate excited levels (upper levels). f JiJm is the oscillator strength between the lower and upper (intermediate) levels. The exponents r im and r f m determine the sign of each term in the summation and are given in Stenflo (1994Stenflo ( , 1997. For Lyα scattering, J i = J f and (−1) rim+r f m = 1. The components of the polarization basis vectors (e 1 and e 2 ) represented in terms of the "complex" spherical unit vectors, e 0 = e z andē ±1 = ∓ ē x ± iē y / √ 2, are defined by ε α q = e α ·ē q , ε α q = e α ·ē q (A5) for α = 1, 2 and q = 0, ±1. Here, the primed quantities, as in ε α q , refer to the incident photon and the unprimed quantities the scattered photon. To derive explicit expressions for w αβ , we adopt a coordinate system and basis vectors similar to those of Stenflo (1994), which are shown in Figure A1. Assuming that the photon is incident in the x-direction and that the scattering takes place in the x−y plane into a scattering angle of θ, the spherical vector components of the polarization vectors of the incident (ε ) and scattered (ε) photons are given by The line profile, in terms of frequency ν, of Lyα is given by Here, we ignored a proportional constant. For Lyα, we use J i = J f = 1/2, J m = 1/2, 3/2, µ i = ±1/2, and µ m = ±1/2 (for J m = 1/2) or ±1/2, ±3/2 (for J m = 3/2). Then, after a straightforward mathematical manipulation of the 3-j symbols, we obtain the following scattering matrix elements for D (4) : The scattering matrix for the Stokes vector S S94 is obtained to be where Best-fit parameters of the empirical formulae, Equation (B1), to the scattering matrix elements numerically-obtained for the dust scattering of Lyα. The albedo a and the asymmetry factor g for the case of using only a single H-G phase function are also shown.   Figure B1. Scattering matrix elements for the MW, LMC, and SMC dust models of Weingartner & Draine (2001), Draine (2003a), and Draine (2003b). Black curves denote numerically obtained results as functions of cos θ. Red curves show the best-fit curves given by Equation (B1). Green curves in the leftmost panels represent the best-fit results for the case of using only one H-G phase function. The best-fit parameters are shown in Table B1. all elements of the scattering matrix for the MW, SMC, and LMC dust models of Weingartner & Draine (2001), Draine (2003a), and Draine (2003b), which are composed of carbonaceous and silicate grains, using a Mie scattering code. 13 In the calculation, we also considered the contribution of free electrons in the optical constants of graphite, as described in Draine & Lee (1984 White (1979) proposed empirical equations to represent the elements of the scattering matrix for the MW dust model of Mathis et al. (1977). We thus attempted to reproduce our numerical results of the scattering matrix elements using their 13 https://www.astro.princeton.edu/~draine/scattering.html equations. However, we found that their equations are not adequate for the dust models of Weingartner & Draine (2001). Instead, we slightly modified the equations and found that the following formulae, obtained through the least-square fitting, reasonably well reproduce the numerical results: S 11 = s 1 φ HG (cos θ; g 1 ) + (1 − s 1 )φ HG (cos θ; g 2 ) S 12 S 11 = −p 1 − cos 2 θ 1 + s 2 2 cos 2 (θ − θ 0 ) S 33 S 11 = s 3 2 cos θ 1 + cos 2 θ + (1 − s 3 ) cos θ + s 4 (cos 2 θ − 1) S 34 S 11 = p c 1 − cos 2 θ 1 + s 5 exp −s 6 θ/π 1 + cos 2 θ 1 + s 5 exp −s 6 θ/π .
Here, φ HG (cos θ; g) = 1 2 1 − g 2 1 + g 2 − 2g cos θ 3/2 (B2) is the H-G phase funtion, p the maximum linear polarization level, and p c the maximum circular polarization level. The best-fit parameters of the empirical equations are shown in Table B1. The table also shows the albedo, and the asymmetry factor g = cos θ for the case of using only one H-G function. Figure B1 compares the numerical results and the best-fit curves of the scattering matrix elements. The deviation between the approximate formulae and the numerical values of S 11 at cos θ ≈ −1 causes no significant difference in RT results because the backward scattering occurs rarely. Note that all we need in the Monte Carlo RT simulation are not the matrix elements themselves but S 11 , S 12 /S 11 , S 33 /S 11 , and S 34 /S 11 . Interestingly, the functional shapes of S 12 /S 11 and S 33 /S 11 are similar to those for the scattering case by hydrogen atoms. However, S 11 for the scattering by dust is totally different from that due to hydrogen atoms; the scattering by dust is predominantly forward-directed, whereas that by hydrogen is relatively more or less isotropic. In the wavelength range dealt with in Lyα RT, the dust extinction and scattering properties are not altered as the photon wavelength changes.

C. SCATTERING ANGLE FOR LYα LINE SCATTERING
The phase function for Lyα scattering by a hydrogen atom depends on the frequency of incident Lyα photon (Stenflo 1980(Stenflo , 1994. This appendix briefly describes an efficient algorithm, implemented in LaRT following the approach described in Seon (2006), to randomly sample scattering angles from the frequency-dependent phase function. The phase function for the angle θ between the incident and scattered direction is: where µ ≡ cos θ. For E 1 = 0, this gives the isotropic scattering, and the scattering angle can be sampled by µ = 2ξ − 1 for a uniform random number ξ between 0 and 1. In general, the scattering angle can be obtained by inverting the cumulative integral µ −1 P(µ )dµ = ξ.
This equation results in the following cubic equation: where p ≡ 4 − E 1 3E 1 and q ≡ 8ξ − 4 E 1 .
By setting µ = 2 |p| 1/2 y, we obtain the standard form to solve a cubic equation  Figure D1. Geometry for the density matrix formae in a fixed coordinate system. The photon with the initial propagation direction k (OP) is scattered in the direction k (OP ). Note that the orientation of polarization basis vectors is defined as left-handed in the density matrix formalism; thus, the third basis vector is given by e 1 = k × e 2 , but not shown in the figure for clarity. The angle (φ or φ ) between the meridian plane (OPQ or OP Q) and the scattering plane (OPP ) is measured in counterclockwise. Those with the prime symbol are for the scattered photon, as oppose to the definition of Chandrasekhar (1960).

D. DENSITY MATRIX FORMAE IN A FIXED FRAME
Instead of the Stokes parameters I, Q, U, and V , Ahn et al. (2002), Ahn & Lee (2015), Chang et al. (2017), Eide et al. (2018), and Chang & Lee (2020) use the density matrix. In their approach, no circular polarization is assumed, i.e., V = 0 and E 1 E * 2 = E * 1 E 2 , where E 1 and E 2 are the electric field components along two directions e 1 and e 2 at right angles to each other. The 2 × 2 Hermitian density matrix D (or coherency matrix) is related to the Stokes parameters by D ≡ ρ 11 ρ 12 where I 1 = |E 1 | 2 and I 2 = |E 2 | 2 are the intensities in the directions e 1 and e 2 , respectively. Therefore, the density matrix is in fact equivalent to the Stokes vector. In this Appendix, we derive the density matrix formulae, given by Ahn et al. (2002) and others, starting from ours.
The most critical difference between the density matrix approach and our method is that they represent the polarization state in a "fixed," left-handed coordinate system; however, this paper expresses it in the photon's local frame. They represent the direction vector of the scattered photon by the polar angles (ϑ and ϕ) measured in a fixed reference frame (i.e., galaxy system). However, in the present approach, we use the angles (θ and φ) measured in the incident photon's local frame. From this point of view, we may refer to the aforementioned "density matrix" approch as the "fixed frame" method and ours as the "local frame" method.
Coordinate systems for the fixed-and local frame methods are illustrated in Figure D1 (see also Figure 8 in Chandrasekhar 1960). In the figure, the polar angles of the incident and scattered photons are (ϑ, ϕ) and (ϑ , ϕ ), respectively. The axes of the photon's local frame are given as follows: e 1 = − sin ϕx + cos ϕy e 2 = cos ϑ cos ϕx + cos ϑ sin ϕy − sin ϑz k = sin ϑ cos ϕx + sin ϑ sin ϕy + cos ϑz, where x, y, and z are axes of the fixed frame. Here, we note that (e 2 , e 1 , k) corresponds to (m, n, k).
To derive the formulae for the density matrix in a fixed frame, we define a column vector: Here, I 2 and I 1 are intensities in two perpendicular directions (e 2 and e 1 ) in a left-handed coordinate system, which correspond to the intensities I m = |E m | 2 and I n = |E n | 2 , respectively, in our convention. The scattering matrix for I can readily be obtained, using the scattering matrix (Equation The transformation matrix of I for a rotation of axes by an angle φ is found, from Equation (14), to be L(φ) =     cos 2 φ sin 2 φ 1 2 sin 2φ 0 sin 2 φ cos 2 φ − 1 2 sin 2φ 0 − sin 2φ sin 2φ cos 2φ 0 (see Equation (190) in Chandrasekhar 1960). The four intensity vector I after a scattering through the angles (θ, φ) are expressed by, in a fixed frame: where the phase matrix can be written as follows: Here, the rotation angle φ denotes the angle between the meridian plane OPQ through and the plane of scattering OPP , and φ the angle between the planes OP Q and OPP in Figure D1. The phase matrix can be decomposed into three terms: P ≡P 1 +P 2 +P 3 P 1 = 2 3 E 1 ×     (mm) 2 (nm) 2 (mm)(nm) 0 (mn) 2 (nn) 2 (mn)(nn) 0 2(mm)(mn) 2(nn)(nm) (mm)(nn) + (nm)(mn) 0 Here, we introduced the abbreviations, similar to those of Chandrasekhar (1960) but for a right-handed system: (mm) = cos θ cos φ cos φ + sin φ sin φ (nn) = cos θ sin φ sin φ + cos φ cos φ (nm) = cos θ cos φ sin φ − sin φ cos φ (mn) = cos θ sin φ cos φ − cos φ sin φ.
We need to multiply our scattering matrix elements by a factor of 8 to obtain Equation (11) of Ahn & Lee (2015) and by a factor of 2/3 for Equation (12). We also note that there was a typographical error in the last term for ρ 22 in Equation (5) of Ahn et al. (2002), Equation (11) of Ahn & Lee (2015), and Equation (21) of Eide et al. (2018). The last term for ρ 22 should read 6 sin 2∆ϕ cos 2 ϑ cos ϑ + 2 sin ∆ϕ cos ϑ sin ϑ sin ϑ ρ 12 .
The numeric factor "6" should come out in front of the parenthesis. The typographical error was fixed in Chang & Lee (2020). The above equations for the density matrix allow us to perform Monte-Carlo simulations for an arbitrary E 1 using an appropriate rejection method.
E. TECHNIQUE USING 100%-POLARIZED LIGHT  describe a Monte-Carlo simulation technique, initially developed by Angel (1969) for Thomson scattering in X-ray, to study the polarization state of Lyα. In their approach, every photon packet is assumed to be 100% linearly polarized along a unit vector perpendicular to its propagation direction. Light is 100% polarized if it is strictly monochromatic. However, it is not generally so because, in actual situations, light is a superposition of mutually-incoherent monochromatic light beams. Partial polarization arises due to the incoherent superposition of mutually uncorrelated photons with different polarization states in a statistical ensemble of photons. The photons are mutually uncorrelated since the atomic processes by which they are created are stochastically independent of each other. Therefore, to obtain unpolarized, quasi-monochromatic light by employing the 100%-polarized photons, they let the photon's polarization direction be randomly oriented. The algorithm was devised for Rayleigh scattering (and additionally, the resonance transition S 1/2 − P 1/2 ; E 1 = 1). Then, it was later extended to the case of the transition S 1/2 − P 3/2 (E 1 = 1/2), which can be regarded as a superposition of Rayleigh and isotropic scattering (Dijkstra & Loeb 2008;Trebitsch et al. 2016).
To obtain a random, new propagation direction after a scattering event,  use the angle Ψ between the polarization and scattered directions, instead of the scattering angle θ defined by the incident and scattered directions. Suppose that β is an angle between the incident photon's polarization vector and the scattering plane, measured counterclockwise when viewed toward the photon. The intensities parallel and perpendicular to the scattering plane are, then, given by I m = I cos 2 β and I n = I sin 2 β. The phase function for 100%-polarized photons can be expressed as a function of Ψ, as follows: P(θ, φ) = I I = 3 4 E 1 cos 2 θ + 1 + E 2 + 3 4 E 1 cos 2 θ − 1 Q I = 3 2 E 1 cos 2 θ cos 2 β + sin 2 β + E 2 = 3 2 E 1 sin 2 Ψ + E 2 , Here, we used the relations Q = I m − I n and cos Ψ = sin θ cos β to obtain the last equation. This indicates that the scattering event can be regarded as a superposition of Rayleigh and isotropic scattering only when E 1 ≥ 0 and E 2 ≥ 0. Otherwise, the phase function cannot be interpreted as a superposition of two probability distribution functions. The frequency range in which the superposition interpretation is invalid (E 1 < 0) is where δν KH ≡ ν K − ν H . This frequency range corresponds to the wavelength range of λ c < λ < λ H in Figure 2. We now provide a method to generate random angles Ψ using the inversion method, as described for the scattering angle θ in Appendix C. Substituting E 2 = 1 − E 1 into Equation (E1), we obtain the distribution function for Ψ: This function results in a cubic equation, which should be solved for cos Ψ, as follows: As in Appendix C, we obtain the following solutions: 2 |p| 1/2 cos cos −1 Q + 4π /3 for E 1 > 0 where In Equation (E5), we, however, note that the solution for E 1 < 0 should not be used in a Monte Carlo simulation using 100%-polarized photons.
One may want to generate random angles Ψ using the composition method for a function composed of the Rayleigh function with a weight E 1 and a uniform distribution function with a weight 1 − E 1 . In this case, random angles following Rayleigh function can be obtained by setting E 1 = 1 in Equation (E5).
A final note we want to mention is that the above technique using 100%-polarized photons is valid only for E 1 ≥ 0. In the case of E 1 < 0, we may assume that E 1 = 0 because the frequency range in which E 1 < 0 is very narrow (especially when the gas temperature is as high as T 10 3 K). If an accurate method is desirable for E 1 < 0, the approach of using partially polarized light, as adopted in LaRT, is required. Another but more elaborated approach is to use the electric field E and the "amplitude" matrix, instead of using the Stokes parameters and the scattering matrix.