Recent cosmological observations indicate that the reionized universe may have started at around z = 6, where a significant suppression around Ly α has been observed from the neutral intergalactic medium. The associated neutral hydrogen column density is expected to exceed 10 ²¹ cm ^-2 , where it is very important to use the accurate scattering cross section known as the Kramers-Heisenberg formula that is obtained from the fully quantum mechanical time-dependent second order perturbation theory. We present the Kramers-Heisenberg formula and compare it with the formula introduced in a heuristic way by Peebles (1993) considering the hydrogen atom as a two-level atom, from which we find a deviation by a factor of two in the red wing region far from the line center. Adopting a representative set of cosmological parameters, we compute the Gunn-Peterson optical depths and absorption profiles. Our results are quantitatively compared with previous work by Madau & Rees (2000), who adopted the Peebles approximation in their radiative transfer problems. We find deviations up to 5 per cent in the Gunn-Peterson transmission coefficient for an accelerated expanding universe in the red off-resonance wing part with the rest wavelength Δλ ∼ 10 Å Quasar absorption systems have been excellent tools to investigate the intergalactic medium (IGM), from which it has been well known that the IGM in the nearby universe is highly ionized (e.g., Peebles 1993). Since the universe after the recombination era z ∼ 1100 should be mostly neutral, there must be some epoch when the universe began to be re-ionized. Intensive studies have been performed on the emergence of the first objects that ended the dark age of the universe. Numerical calculations adopting the cold dark matter models predict a reionization epoch at z ∼ 6−12 (e.g., Gnedin & Ostriker 1997). Around this epoch a broad absorption trough in the blue part of Ly α is expected and regarded as a strong indicator of the reionization of the universe, which was predicted by Gunn & Peterson (1965) and indepen- dently also by Scheuer (1965). With the advent of the Hubble Space Telescope and 8 meter class tele- scopes there have been extensive searches for the Gunn- Peterson trough in the spectra of high redshift objects. Remarkable contributions have been made by the Sloan Digital Sky Survey, from which a number of high red- shift quasars with z ranging from 4 to 6 have been found. According to recent results based on Keck spec- troscopy of these high redshift quasars (Becker et al. 2001; Fan et al. 2006a), the flux level drop around Ly α is much higher for a quasar with z = 6.28 than those for other quasars with z < 6, which indicates Corresponding author: H.-W. Lee that the reionization epoch may be found at around z ∼ 6. During the last decade, numerous studies of cos- mic reionization have been discussed by using analytic models (Madau, Haardt & Rees 1999; Miralda-Escudé, Haehnelt & Rees 2000; Wyithe & Loeb 2003), and di- rect numerical simulations (Gnedin 2000; Sokasian et al. 2002; Razoumov et al. 2002; Ciardi et al. 2003; Iliev et al. 2006; Ahn et al. 2012). Comprehensive reviews of the cosmic reionization are found in the literature (see Loeb & Barkana 2001; Ciardi & Ferrara 2005; Fan et al. 2006b; Morales & Wyithe 2010; Bromm & Yoshida 2011; Fan 2012). The exact computation of the flux drop around Ly α requires an accurate atomic physical estimation of the scattering cross section. Recent theoretical works on the calculation of the Gunn-Peterson trough were pro- vided by Miralda-Escudé (1998), who adopted the for- mula that was introduced in a heuristic way using the second order time-dependent perturbation theory by Peebles (1993). The formula is derived based on the as- sumption that the hydrogen atom is a two-level atom, in order to describe the behavior of the scattering cross section that is approximated by the Lorentzian near res- onance and yields a ω ⁴ -dependence in the low energy limit. As Peebles noted explicitly, due to the two-level assumption, the formula provides an inaccurate propor- tionality constant in the low energy limit, even though it correctly gives the ω ⁴ -dependence. The accurate cross section should be obtained from the second order time-dependent perturbation theory treating the hydrogen atom as an infinitely many-level atom including the continuum free states, which is known as the Kramers-Heisenberg formula. The dis- crepancy between the two formulae will be significant in the far off-resonance regions where the contribution from the np ( n > 2) states including the continuum states becomes considerable. Therefore, in order to ob- tain an accurate Gunn-Peterson profile it is essential to investigate the exact scattering optical depth of a medium with a high neutral hydrogen column density. In this paper, we present a reliable atomic physics re- lation that governs the scattering around Ly α by intro- ducing the Kramers-Heisenberg dispersion relation and a simple fitting formula around Ly α . With the well- defined cross section schemes, we compute the Gunn- Peterson trough profiles adopting a representative set of cosmological parameters with our choice of the reion- ization epoch and make quantitative comparisons with previous works. The interaction of photons and electrons is described by the second order time-dependent perturbation the- ory, of which the result is summarized as the famous Kramers-Heisenberg dispersion relation (Kramers & Heisenberg 1925). As is well illustrated in a typical quantum mechanics text (see Heitler 1954; Bethe & Salpeter 1957; Sakurai 1967; Berestetskii et al. 1971), it may be written as where are the polarization state vectors as- sociated with the incident photon and outgoing pho- ton, respectively, Γ _I is the radiation damping term as- sociated with the intermediate state I, ω is the inci- dent angular frequency, ω_IA is the angular frequency of the transition between I and the ground state A , and r _o = e ² /( m_ec ² ) = 2.82 × 10 ⁻¹³ cm is the classical elec-tron radius (e.g., Sakurai 1967). Here, the electron is in the ground state A before scattering and de-excites to the same state. The sum- mation (and integration) should be carried over all the intermediate states I including the infinite number of bound states and the free or continuum states. For hy- drogen, the dipole moment matrix elements can be ex- plicitly written using the hypergeometric function (e.g., Bethe & Salpeter 1957; Karzas & Latter 1961; Berestet- skii et al. 1971). In the blue part of Ly α , it is possible that the scat- tering atom may de-excite to the excited 2 s state by re-emitting a photon with much lower frequency than the incident photon. This Raman scattering is negligi- ble near Ly α due to small phase space available for an outgoing photon. However, this process becomes im- portant as the incident photon energy increases. In the range where the current work is concerned, the Raman scattering process is safely neglected. Despite the existence of explicit analytic expressions for each matrix element that is the Kramers-Heisenberg formula for hydrogen, it is still cumbersome to use the formula as it is. The Lorentzian function as an exact so- lution of the first order perturbation theory (Weisskopf & Wigner 1930) provides an excellent approximation at the line center, where Γ _2p = 6.25 × 10 ⁸ s ⁻¹ is the radiation damping constant associated with the Ly α transition. In prin- ciple, it provides a good description of physical pro- cesses, such as the natural line broadening, transition life time, and decay rates of the excited state during the radiation-atom interaction. However, the linear per- turbation theory cannot give a proper solution for the transition probability for off-center scattering. For extremely dense Ly α clouds, the line center is easily saturated and the broad damping part becomes important. In order to describe the broad damping profiles of the Ly α clouds, Peebles (1993) introduced a heuristic formula based on the two-level approximation, In his qualitative analysis, resonance scattering is smoothly connected to the off-center scattering with the ω ⁴ –dependency of classical Rayleigh scattering. As he pointed out, since the contribution of all possible inter- mediate states is neglected, the two-level approximation results in significant underestimation of the transition probabilities (by a factor of 2 in the low energy limit). In previous studies, fitting relations for Rayleigh scat- tering have been introduced. Through their quantal analysis, Gavrila (1967) provided the fitting polynomial with the threshold frequency (λ > 1410 Å), where σ_T = 8π r ₀ ² /3 is the Thomson scattering cross sec- tion. In particular, Ferland (2001) used Gavrila’s fit in his photoionization code Cloudy. Earlier, Baschek & Scholz (1982) has also introduced a similar expression in the wavelength range λ > 2000 Å. However, their fit- ting relations are not smoothly connected to resonance scattering (the Lorentzian), and exhibit about 4% oc- casional errors. Moreover, no simple polynomial fits to the blue wing of Ly α are available. Instead, we investigate an analytic correction func- tion which corrects the classical approximation in prac- tical applications. Adopting a non-linear regression, a reliable correction function and its coefficients have been numerically determined. Between the wavelength interval 1100 Å < λ <8000 Å which covers the whole practical wavelength domain of Ly α , our correction function has the form of where σ_L is the Lorentzian function (Equation (2)), and f ( ω ) is deviation between the Kramers-Heisenberg for- mula and the cross section of the classical Rayleigh scat- tering. The deviation function has been fitted by the analytic form with the coefficients where x = 1 − ( ω/ω ₀ ). Our fitting relation satisfies two asymptotic properties: At the resonance frequency ( ω/ω ₀ → 1), it approaches the Lorentzian function ( f ( ω ) → 0). For the long wavelength limit of the Ly α scattering, it properly corrects the classical Rayleigh scattering (<0.5% error). Moreover, it also extends to the blue wing of Ly α with excellent accuracy (< 1% error). Our formula for the full quantum mechanical cross section has been compared with the classical approx-imations and the fitting relations. In Figure 1 , the Ly α cross section of off-center region (top) and near- resonance (bottom) is presented. The solid line denotes the Kramers-Heisenberg relation, the Lorentzian (dot- ted), the Peebles approximation (dashed), Gavrila’s fit (dot-dashed), and our fitting function (crosses), respec- tively. The Lorentzian provides an excellent approx- imation at the resonance. However, its slow damp- ing rate (∼ ω ⁻² ) leads to overestimation of scatter- ing probabilities at low energies. This asymptotic be- havior implies that the radiation damping phenomenon can not be properly described by linear perturbation theory. In the bottom panel, it is notable that the deviation of the Lorentzian is slightly anti-symmetric with respect to the line center. The Lorentzian un- derestimates transition probabilities in the wavelength region λ ₀ < λ < 1353 Å. The largest deviation in this interval is about 3.5% underestimation at λ = 1274Å, which may cause an asymmetry problem inden-tification of the DLAs line profiles (Lee 2003). Es-pecially, 1% deviation of the Lorentzian from the in- trinsic cross section is presented at a wavelength shift Δλ = λ − λ _α = ±3.3 Å where the corresponding cross section is σ = 3.8 × 10 ⁻²¹ cm ² . This implies that the Voigt fitting has only 99% accuracy for a column den- sity smaller than 3 × 10 ²⁰ cm ⁻² . In the case of the Peebles approximation (Equa- tion (3)), resonance scattering is smoothly connected to the Rayleigh part. At low energies ( ω ≪ ω _α ), it produces ω ⁴ -dependence which corresponds to the clas- sical result. On the other hand, it also becomes the Lorentzian at the line center ( ω → ω ₀ ). However, the two-level assumption leads to underestimation of the transition probability by neglecting all possible inter- mediate states. Basically, oscillator strength is a good measure of contribution of individual excited state. For low energy photons ( ω < ω ₀ ), the contribution of the 2 p state is comparable to that of the whole remaining states. Therefore, the Peebles approximation cannot provide an accurate description of the radiation damp- ing phenomenon. Combining the accurate Ly α cross section with stan- dard cosmology, we compute intergalactic absorption profiles of the reionization epoch. The presence of neu- tral hydrogen during reionization is signalled by the suppressed spectra around Ly α , that is, the Gunn– Peterson troughs (Gunn & Peterson 1965; Scheuer 1965). The ionizing radiation from the first lumi- nous objects was scattered by the intervening neutral IGM. Then, even a small amount of neutral fraction (∼10 ⁻⁴ ) can produce a considerable Ly α opacity (Loeb & Barkana 2001; Fan et al. 2006b). We compute the Gunn-Peterson optical depth de- fined by using the Kramers-Heisenberg formula and make com- parisons with previous works performed by Miralda- Escudé (1998); Loeb & Rybicki (1999); Madau & Rees (2000). Here, λ _obs is the observed wavelength, z_r, z_s are the redshifts of the complete reionization of the universe and the reionizing source, respectively, and n ( z ) = n ₀ (1+ z ) ³ is the homogeneous neutral hydrogen density at redshift z . We choose z_r = 6, z_s = 7 as in Madau & Rees (2000), but do not consider the prox- imity effect of the ionizing source. The Gunn-Peterson optical depth can be written as where Ω _M , Ω _Λ are the density parameters due to matter and the cosmological constant and the characteristic hydrogen column density defined as Following previous work (Madau & Rees 2000), we choose the present Hubble constant and the hydro- gen number density H ₀ = 50 km s ⁻¹ Mpc ⁻¹ , n ₀ = 2.4 × 10 ⁻⁷ cm ⁻³ so that N _HI,0 = 4.3 × 10 ²¹ cm ⁻² . In Figure 2 , the Gunn-Peterson transmission co-efficient defined as T_GP ≡ exp(−τ _GP ) is presented. We assumed two cosmological models: (top panel) the Einstein-de Sitter universe (Ω _M , Ω _Λ )=(1, 0) and (bottom panel) the Λ-cold dark matter universe (Ω _M , Ω _Λ )=(0.35, 0.65). In terms of the characteristic Gunn-Peterson optical depth τ _GP,0 defined as our choice of parameters in the case of the top panel yields the value of τ _GP,0 = 3 × 10 ⁵ at z_s = 7 as was adopted in the work of Madau Λ Rees (2000). The horizontal axis represents the logarithm of the normal- ized wavelength ratio δ defined as thus δ = 0 corresponds to the resonance wavelength of Ly α . The dotted lines represent the Gunn-Peterson transmission coefficient based on Equation (3). The deviation between the two formulae is notable around δ = 10 ⁻² , where the deviation is about 3 per cent in the top panel and 5 per cent in the bottom panel. Peebles’ approximation turns out to be pretty good for contem- porary applications. However, the deviation between the two formulae will increase as n ₀ or N _HI,0 increases, because the discrepancy of the Kramers-Heisenberg for- mula and the Peebles approximation becomes larger as the frequency is further away from the line center. Near resonance, both formulae are excellently approximated by the same Lorentzian. Therefore, no significant de- viation is expected when the neutral medium is of low column density ≲ 10 ²¹ cm ⁻² . It is notable that an ac- curate treatment of atomic physics is more important in an accelerated expanding universe where the universe was more compact than the universe without the cos- mological constant. In practice, the Gunn-Peterson test provides a unique tool which is highly sensitive to the line-of-sight distri- bution of neutral medium. Indeed, the Ly α absorp- tion profiles can be affected by the line-of-sight distri- bution of neutral medium, or by a complex geometry of ionization fronts of H _II bubbles. Adopting recent cosmological concordance values, we examine the effect of density distribution on the Gunn-Peterson profiles ( Figure 3 ). Firstly, assuming a homogeneous distri- bution of the neutral IGM along the redshift interval [ z_r , z_s ]=[6, 8], the Ly α red wings are computed by using two cross section schemes of the Kramers-Heisenberg relation (solid) and the Peebles approximation (dotted). In addition, we also consider inhomogeneity of the line- of-sight distribution of the neutral medium by adjusting the fraction of neutral hydrogen ( f _HI = n _HI / n _H ). The (red) circle represents the transmitted flux based on the Kramers-Heisenberg formula but with the neutral fraction of f _HI = 0.2 in the redshift interval [6.0, 7.6] and with f _HI = 1.0 in [7.6, 8]. This slab-shaped distri- bution where the neutral medium is highly localized is expected to be the last stage of cosmic reionization at z = 6 ∼ 8. In the blue part spectra of the rest frame of the source, transmitted fluxes are swept out by the high optical depth of the near-center scattering. Then, the neutral IGM at lower redshift than source location can affect the Lyα opacities in the far off-center of red damping wing with a small transition probability. This implies that a cumulative effect by the train of neu- tral hydrogen along the line-of-sight can be present in the far off-center region of the red damping wing. It also suggests that the classical approximation can not provide an accurate measurement of the cumulative ef- fect or a proper continuum level of transmission. If the Lyα emitter is embedded in a mostly neutral IGM, the Gunn-Peterson profile will give a substantial informa- tion of the density distribution near the source. In this paper, we have investigated the behavior of the scattering cross section around Lyα in a quantitative way, where the deviation from the Lorentzian becomes significant as the incident frequency gets further away from the line center. Therefore, in an analysis of the Gunn-Peterson trough profile, which is associated with a neutral medium with a high H _I , column density, an inaccurate treatment of the atomic physics of hydrogen may introduce significant errors in estimating impor- tant cosmological parameters including the epochs of the emergence of the first objects and the completion of the reionization of the universe. Voigt profile fitting has been very successfully ap- plied to quasar absorption systems with a broad range of H _I column densities. However, the deviation of the true scattering cross section from the Lorentzian ex-ceeds 1% when the relevant column density becomes N _HI ≥ 3 × 10 ²⁰ cm ⁻² that is the typical column den- sity of a damped Lyα absorber. This is especially im- portant in some damped Ly α systems that may possess N _HI > 10 ²¹ cm ⁻² (e.g., Turnshek & Rao 1998). Recently, numerous theoretical works of reionization processes have been studied by using analytic mod- els (Madau, Haardt & Rees 1999; Miralda-Escudé, Haehnelt & Rees 2000; Wyithe & Loeb 2003), and numerical simulations (Gnedin 2000; Sokasian et al. 2002; Razoumov et al. 2002; Ciardi et al. 2003; Iliev et al. 2006; Ahn et al. 2012). The intergalactic medium underwent a dramatic change in the chemical states and the morphological structures throughout the whole epoch of cosmic reionization (Loeb & Barkana 2001; Fan et al. 2006b). Especially, the neutral column den- sity distribuion from the DLAs is known to be highly inhomogeneous in the redshift range 1 ≤ z ≤ 5 af-ter reionization (Noterdaeme et al. 2012; Zafar et al. 2013). Moreover, recent spectroscopic observations of Lyman break galaxies (LBGs) in the redshift interval 6 ≤ z ≤ 7 indicate a constant decrease of the Lyα emit- ter density (Ota et al. 2010; Pentericci et al. 2011). In their works, this decrease was interpreted as due to a combination effect by the early evolutionary stages of galaxies and by the local distribution of the remaining neutral medium.