Probing Intergalactic Neutral Hydrogen by the Lyman Alpha Red Damping Wing of Gamma-Ray Burst 130606A Afterglow Spectrum at $z=5.913$

The observed spectrum is fitted by models including absorptions by neutral hydrogen along the sightline, imprinted on the baseline power-law with an index $\beta$ ($f_{\nu }\propto {\nu }^{\beta }$). We consider three components of H\emissiontypeI absorption: H\emissiontypeI in the host galaxy, diffuse H\emissiontypeI in IGM, and an intervening DLA that is often found in high-redshift quasar/GRB spectra.

The host H\emissiontypeI component is assumed to have a single Gaussian radial velocity dispersion ${\sigma }_v$ with a column density $N_{H\text{emissiontype}I}^{host}$. In our baseline analysis, we fixed the redshift of host H\emissiontypeI  as $z_{GRB}=5.9131$, from the three low-ionization metal lines found in the same VPH900 spectrum used for the damping wing fit. (We found Si\emissiontypeII $\lambda$1260.42, O\emissiontypeI $\lambda$1302.17, and C\emissiontypeII $\lambda$1334.53 and their redshifts are $z=5.91305\pm 0.00002$, $5.91307\pm 0.00005$, and $5.91270\pm 0.00002$, respectively. The C\emissiontypeII redshift may be affected by a slightly blending nearby line.) The widths of these lines are $\sim$5 Å FWHM that is slightly extended beyond the spectral resolution (3.9 Å), indicating a velocity dispersion of $\sim$100 km/s FWHM (or 43 km/s in $\sigma$). Therefore we perform the fit in the limited range of ${\sigma }_v=$ 0–100 km/s, though the fit results are almost insensitive to ${\sigma }_v$ (see next section).

We consider a component of constant-density diffuse IGM with the damping wing formula of Miralda-Escude (1998) with a neutral fraction $f_{H\text{emissiontype}I}$ extending in a redshift range from $z_{IGM,l}$ to $z_{IGM,u}$. We adopted the GP optical depth ${\tau }_{GP}=3.97\times {10}^5{f}_{H\text{emissiontype}I}\left[\right(1+z)/7{]}^{3/2}$ calculated from the latest estimates of cosmological parameters (Komatsu et al. 2011) and the primordial helium abundance (Peimbert et al. 2007). Though the distribution of neutral hydrogen in IGM should be clumpy and inhomogeneous in reality, we test this simple model to get a rough estimate for mean $f_{H\text{emissiontype}I}$ required to explain the observed damping wing. Dark Ly$\alpha$ troughs are found at $z\gtrsim 5.67$ along the sightline to this GRB (Chornock et al. 2013), and hence we set the lower bound of the H\emissiontypeI distribution as $z_{IGM,l}=5.67$, but our results do not significantly change by choosing different values, because the contribution to the damping wing from H\emissiontypeI  around $z_{IGM,l}$ is small. For the upper bound, we will test two models; one is assuming $z_{IGM,u}=z_{GRB}$, while the other treats $z_{IGM,u}$ as a free parameter.

If we find an evidence for absorption by intervening H\emissiontypeI gas outside the host galaxy, it may be a result of an intervening DLA rather than diffuse IGM. Low-ionization metal absorption lines are generally associated to DLAs, and this would give a constraint on the DLA scenario. For the GRB 130606A spectrum, a metal absorption system at $z=5.806$ is found by Si\emissiontypeII, O\emissiontypeI, and C\emissiontypeII lines (Chornock et al. 2013), which is a good candidate of the DLA redshift. We will thus perform a fit to the observed damping wing with the DLA component at $z_{DLA}=5.806$, with a column density of $N_{H\text{emissiontype}I}^{DLA}$.

Therefore fits are performed with the free model parameters of $\beta$, $N_{H\text{emissiontype}I}^{host}$, ${\sigma }_v$, $f_{H\text{emissiontype}I}$, $z_{IGM,u}$ and $N_{H\text{emissiontype}I}^{DLA}$, while the overall normalization factor of the model is always chosen to fit the data (i.e., marginalized). The theoretical fluxes are calculated for each data point (i.e., pixel) of the VPH900 spectrum, then corrected for the Galactic extinction ($A_V=0.066$ mag, Schlafly & Finkbeiner 2011), and finally convolved with a Gaussian profile of the instrumental wavelength resolution (3.9 Å FWHM).

We do not consider possible extinction by dust in the host, because standard extinction curves have a monotonic wavelength dependence within the short range analyzed here, and adding nonzero extinction is almost perfectly absorbed by an according change of the index $\beta$, giving a practically identical fit (see Fig. 2 where the absorption profiles of various components are presented). Therefore the index $\beta$ here is not the intrinsic value for the afterglow emission, but should be regarded as an effective value in the analyzed wavelength range including extinction at the host. However, it is also useful to estimate the amount of extinction expected from the hydrogen column density in the host, $N_{H\text{emissiontype}I}\sim 7\times {10}^{19}c{m}^{-2}$. Using the mean relations for the Milky Way (MW, Predehl & Schmitt 1995) and the Small Magellanic Cloud (SMC, Gordon et al. 2003), this column density translates into $A_V=0.036$ and 0.0048, respectively. The value from the MW relation is likely an overestimate for this GRB, because the optical spectra of this GRB indicate a metallicity lower than 0.3 or 1/7 times the solar abundance (Chornock et al. 2013; Castro-Tirado et al. 2013b). It should also be noted that many GRB afterglows show smaller extinction than those expected from soft X-ray absorption or the $N_{H\text{emissiontype}I}$-$A_V$ relations (Schady et al. 2010), implying that the above estimate of $A_V$ is likely an overestimate. For reference, $\Delta \beta$ = $-$0.042 and $-$0.094 are required to represent the reddening of $A_V=0.01$ in 8400–8900 Å by the change of $\beta$, for the MW and SMC extinction curves, respectively.

Other independent observational constraints are available for the power-law index $\beta$ from near-infrared photometric observations of the afterglow. In Fig. 3 we plot the reported optical/near-infrared photometric observations (Nagayama et al. 2013; Butler et al. 2013a,b; Morgan et al. 2013; Afonso et al. 2013), showing a constant value of $\beta \sim -1$ at wavelengths longer than the Lyman break from 36 min to 35 hr after the burst. Using observations at the time close to the Subaru observation, we find $\beta =-0.94\pm 0.05$ from the GROND $JHK$ band data (8 hrs) and $-0.78\pm 0.10$ from the PAIRITEL $JH{K}_s$ band data (9.4 hrs). As mentioned above, extinction at the host should be at most $A_V\sim 0.01$, and the expected change of $\beta$ between the near-infrared bands and our wavelength range should be $\Delta \beta \lesssim 0.05$. The best-fit $\beta$ values obtained by our fits presented below (see Table 2) are roughly consistent with these independent measurements, taking into account the systematic uncertainties about host extinction and in estimating $\beta$ from photometric measurements (e.g., band filters, response function, or zero-point calibration).

The spectrum has been reduced in the range of 8377–8902 Å, but we further removed the wavelength shorter than 8426 Å from the fitting analysis, because this region is dominated by the host H\emissiontypeI  absorption rather than the other interesting components (Fig. 2). This is the region where the damping wing rapidly drops with decreasing wavelength, and a fitting becomes sensitive to the velocity distribution of H\emissiontypeI atoms in the host. Although we included a Gaussian distribution with the ${\sigma }_v$ parameter in our fit, the realistic distribution is unlikely a pure Gaussian (Wolfe et al. 2005). The fit in this range is also sensitive to the uncertainties about the mean redshift of H\emissiontypeI in the host (i.e., $z_{GRB}$) and the instrumental spectral resolution. The lower bound of 8246 Å was determined as the shortest wavelength with which the fitting results are mostly insensitive to ${\sigma }_v$, as can be seen in Table 2. Discernible absorption features were also removed from the fitting as shown in Fig 1. Then 323 pixels remain for the fitting analysis.

The statistical errors are estimated from the noise spectrum produced by the object+sky electron counts and readout noise, which is random without correlation among different pixels. (Note that the FWHM of spectral resolution is larger than one pixel, but the random statistical errors appear in each pixel in addition to the signal convolved by the spectral resolution.) However, pixel-by-pixel statistical fluctuation of the final spectrum is smaller than the noise estimate, because the fluxes at neighboring pixels are mixed to some extent by the resampling for wavelength calibration in the spectrum reduction. At the same time, correlation between the errors of neighbouring pixels is generated, which may affect a ${\chi }^2$ analysis. This problem can be solved by binning the five connecting pixels, and we confirmed that the fit residuals of the best-fit model (the IGM-low$z_u$ model, see below) do not show any correlation between neighbouring bins within the statistical uncertainties. There are 68 bins after the binning, and we can adopt the standard ${\chi }^2$ analysis to calculate the best-fit values and the statistical confidence regions of the model parameters.

We first examine the simplest model including only H\emissiontypeI in the host galaxy (the host-only model hereafter). The best-fit values of the three model parameters ($\beta$, $N_{H\text{emissiontype}I}^{host}$, and ${\sigma }_v$) and the ${\chi }^2$ value are given in Table 2, and the fit residuals $(f_{obs}-f_{model})/{\sigma }_{obs}$ of each pixel are shown in the top panel of Fig. 4. The 1$\sigma$ statistical errors for one model parameter in Table 2 were calculated by the standard procedure (Press et al. 2007), i.e., finding the parameter range corresponding to $\Delta {\chi }^2\left(p\right)\equiv {\chi }^2\left(p\right)-{\chi }_{min}^2=1$, where ${\chi }_{min}^2$ is the global minimum and $\Delta {\chi }^2\left(p\right)$ is the minimum as a function of the parameter of interest, $p$, with the other model parameters kept free (i.e., marginalized). For convenience, we define the four wavelength ranges I–IV (indicated in Figs. 1 and 4).

The most significant excess of the residuals of the host-only model is found in the wavelength range III at a level of about 0.6% of the flux. The trends of increasing/decreasing residuals with increasing wavelength in the ranges II/IV seem to be connected to the excess in the range III, showing the overall systematic trend of concave curvature spanning in the entire analysis wavelength range of 8400–8900 Å. We have examined various sources of systematic uncertainties, but the 0.6% level curvature in this narrow wavelength range cannot be explained (see Section 4).

The ${\chi }^2$ value of this fit is 87.58, and statistically the chance probability of getting this value is 2.7% (i.e., 2.2$\sigma$ rejection) for the 68 data points and four model parameters (i.e., 64 degrees of freedom). Here we examined this for a consistency check, and later we will find that this model is disfavored by about 3$\sigma$ in comparison with other models. This sigma difference is not surprising, because the total ${\chi }^2$ value and ${\chi }^2$ difference for different models are independent tests. For example, if a model parameter is wrong in a model fit but the parameter is affecting only a limited number of ${\chi }^2$ degrees of freedom, the total ${\chi }^2$ value is insensitive and we need to see the ${\chi }^2$ difference when the model parameter is changed.

The concave curvature can be resolved by an absorption component that becomes rapidly stronger at shorter wavelengths, and hence the redward damping wing of another neutral hydrogen component on the sightline is a good candidate to improve the fit. Therefore we test the next simplest model with diffuse IGM H\emissiontypeI extending to the same redshift as the GRB host, i.e., $z_{IGM,u}=z_{GRB}$ (the IGM-$z_{GRB}$ model hereafter). The best-fit of this model with the new free parameter $f_{H\text{emissiontype}I}$ is presented in Table 2 and Fig. 4. The best-fit is obtained at $f_{H\text{emissiontype}I}=0.086$, and the relative statistical significance against the host-only model can be estimated again in the standard manner (Press et al. 2007) as follows. Comparing the best-fit ${\chi }^2$ values of the host-only and IGM-$z_{GRB}$ models in Table 2, there is a reduction of $\Delta {\chi }^2=8.92$ with one new parameter of $f_{H\text{emissiontype}I}$, where all the other parameters are marginalized. Therefore the latter model is preferred at the statistical significance of ${8.92}^{1/2}=3.0\sigma$ from ${\chi }^2$ statistics with one degree of freedom.

However, the curvature of residuals is still evident for this model in Fig. 4. This may be further improved by allowing lower values of $z_{IGM,u}$, because such a model has a weaker wavelength dependence of IGM absorption, and hence the effect can reach to longer wavelength beyond the Ly$\alpha$ break, while keeping the absorption around the break relatively weak. Therefore we performed a fit treating $z_{IGM,u}$ as a further free model parameter (the IGM-low$z_u$ model, hereafter), and the dependence of the minimum ${\chi }^2$, best-fit $\beta$ and $f_{H\text{emissiontype}I}$ on $z_{IGM,u}$ is shown in Fig. 5, where the other parameters of $N_{H\text{emissiontype}I}^{host}$ and ${\sigma }_v$ are marginalized by choosing values that minimize ${\chi }^2$ for each value of $z_{IGM,u}$. We find a modest ($\Delta {\chi }^2\sim 2$, 1.4$\sigma$) improvement of the fit by decreasing $z_{IGM,u}$ down to $\sim 5.83$, with a higher value of $f_{H\text{emissiontype}I}=0.47$.

It is interesting to note that this redshift is coincident with the upper border of the region where the darkest GP troughs of Ly$\alpha$, $\beta$, and $\gamma$ are found to this sightline ($z\sim$ 5.67–5.83, Chornock et al. 2013), i.e., the region where we expect the largest amount of IGM neutral hydrogen. The physical proper distance between $z=5.83$ and the GRB host ($z=5.913$) is 5.3 Mpc. The best-fit model parameters for $z_{IGM,u}=5.83$ are presented in Table 2, and now the systematic curvature has almost disappeared in the fit residuals (Fig. 4). This IGM-low$z_u$ model has two more parameters of $f_{H\text{emissiontype}I}$ and $z_{IGM,u}$ against the host-only model, and $\Delta {\chi }^2=10.96$ is significant by the chance probability of 0.0042 or $2.9\sigma$ from the ${\chi }^2$ statistics of two degrees of freedom. This result again indicates that the host-only model is disfavored at $\sim 3\sigma$ level. On the other hand, the difference between the two IGM models is not statistically significant (1.4$\sigma$), and hence we cannot strongly reject the IGM-$z_{GRB}$ model. The total ${\chi }^2$ value of the model with $z_{IGM,u}=5.83$ is 76.62, corresponding to a modest chance probability of 10.0% for 62 degrees of freedom.

The best-fit host-only and IGM-low$z_u$ models are shown in Fig. 1 compared with the observed flux. To understand these results qualitatively, we also show the flux attenuation factor by various components of absorption in Fig. 2. The key point is that the IGM absorption signature (especially with smaller $z_{IGM,u}$) has a weaker wavelength dependence than that by hydrogen in the host, which is essential to resolve the curvature found in the residual spectrum of the host-only fit.

Another possibility to improve the host-only fit is a DLA rather than diffuse IGM. Therefore we test the intervening DLA model, in which the DLA component at the redshift of the observed metal absorption system ($z=5.806$) is added to the host-only model, but without including diffuse IGM. The results are presented in Figs. 2, 4, and Table 2. It is found that this model gives a similarly good fit to the IGM-low$z_u$ model, if the H\emissiontypeI column density is $N_{H\text{emissiontype}I}^{DLA}=5.1\times {10}^{20}c{m}^{-2}$, which is about 10 times larger than $N_{H\text{emissiontype}I}$ in the GRB host. We will further examine whether this intervening DLA model is consistent with other spectral features in Section 5.

The discrepancy between the host-only model and the observed data is mainly because of the concave curvature of the data compared with the best-fit model curve, extending in the entire wavelength range analyzed (8400-8900 Å) with an amplitude of 0.6% of the continuum flux (Fig. 4). As mentioned above, the sensitivity curve fit to the standard star is accurate by 0.1% level, and the 0.6% level difference of the spectral shape should be detectable.

To further test possible systematic uncertainties in the relative flux calibration, we constructed another sensitivity curve by using another CALSPEC standard star GD 153, which is one of the primary standards of the CALSPEC system and was observed by Subaru one day before the observation of GRB 130606A. The airmass at the time of taking the VPH900 spectrum of GRB 130606A is in between those for the two standard stars (Feige 34 and GD 153). The shapes of the two sensitivity curves are consistent within 0.2% accuracy, indicating that the observed spectral shape is correct at 0.2% level, taking into account the variation of the standard star templates, airmass, and atmospheric/instrumental conditions over the time scale of one day. The primary standards of the CALSPEC calibration system are the model spectra of the three pure hydrogen white dwarfs including GD 153 (Bohlin 2007), and the model spectra have a smooth power-law profile in 8400–8900 Å, showing no complicated absorption structures. Residuals of the HST data from the model spectra indicate that a curvature of 0.6% level in 8400–8900 Å is highly unlikely to be a result of the uncertainty about the CALSPEC template (see Fig. 3 of Bohlin 2007).

In the wavelength range III, a broad enhancement of the spectral noise is seen, which is a feature by the atmospheric O${}_2$ emission lines. If the spectrum is systematically biased by this effect, it would have an important effect on our fitting result, because the deviation of the data from the host-only model is the largest in this range. To check this we examined the fit residuals of the Subaru data for the standard star with respect to the sensitivity curve fit, and there is no systematic excess in this region. The noise level is decreasing to larger wavelength within the range III, but the fit residuals of the host-only model in Fig. 4 do not show any systematic wavelength dependence in this range. The concave curvature of the residuals ranging in the wide range of 8400–8900 Å cannot be resolved only by the O${}_2$ feature in the range III.

We examined the atmospheric absorptions in 8500–8900 Å using the ATRAN database (http://atran.sofia.usra.edu/), and found that the expected curvature produced by weak atmospheric absorption lines is about 0.1% in this wavelength range. It should be noted that the effect of atmospheric absorption should have been corrected if the effect is the same for the standard star and the GRB afterglow. Although a part of this effect may remain uncorrected by a possible systematic change of atmospheric conditions at the observing times for the star and the afterglow, it is impossible to explain the 0.6% level curvature of the host-only fit residuals.

We also examined the effect of slit/aperture loss of afterglow light. The slit loss of light changes the relative flux of the spectrum by the atmospheric dispersion effect, but we confirmed that this is less than 0.05% within the wavelength range analyzed. We found that the image size along the spatial direction changes with wavelength, which is likely due to wavelength dependent seeing and/or chromatic aberration. Therefore we quantitatively examined the systematic change of the spectrum by aperture loss. (The aperture size for the GRB spectra is smaller than that for the standard star.) We found that this changes the relative flux at 8400 and 8900 Å by $\sim$ 2.6%, but the wavelength dependence is monotonic and well described by a power-law, having an effect of changing $\beta$ by 0.45. The analyzed spectrum has been corrected about this before the model fittings, and hence correction for this effect is not necessary for the $\beta$ values presented in this paper. Since this effect is well described by a power-law, it is unlikely to be the origin of the residuals showing a curvature.

We have removed the wavelength regions showing visually discernible absorption features from our analysis, but weaker absorption lines may affect the fit results. Since it is difficult to test the effect of weak lines that cannot be recognized, we instead performed the fittings including the four relatively weak absorption features at 8492, 8505, 8824, and 8889 Å, which were removed in our baseline analysis. Some of these may be statistical fluctuation rather than real lines. The best-fit ${\chi }^2$ values for the host-only and IGM-$z_{GRB}$ models are now 134.69 and 123.11 for 69 data points, and the preference for the IGM model ($\Delta {\chi }^2=11.58$) is not significantly changed. We also tested possible systematic effects induced by large noise in the regions of atmospheric airglow lines. We repeated the fittings after removing the narrow airglow regions centered at 8435, 8455, 8495, 8505, 8542, 8552, 8830, 8840, 8852, 8870, and 8888 Å. (These regions can be seen as bumps in the error curve shown in Fig. 1.) The ${\chi }^2$ values of the two models are 60.17 and 51.73 for 45 data points, and $\Delta {\chi }^2=8.44$ is still significant at 2.9$\sigma$.

We used the central redshift of $z_{GRB}=5.9131$ for the H\emissiontypeI gas in the GRB host, based on the redshifts of the metal absorption lines measured in the same spectrum. The C\emissiontypeII line shows a relatively low redshift of $z=5.9127$, and the uncertainty about the host redshift can be tested by adopting this value for $z_{GRB}$. We then found the best-fit ${\chi }^2$ values of 86.13 and 79.30 for the host-only and IGM-$z_{GRB}$ models, respectively, and $\Delta {\chi }^2=6.83$ is still significant at 2.6$\sigma$.

We considered only a single power-law for the intrinsic afterglow spectrum before absorption effects, though afterglows sometimes show concave breaks in their spectral energy distribution (SED), which may explain the curvature found in the residuals of the host-only fit. However, if a spectral break is close to the observed wavelength range, we expect from the afterglow theory a significant time evolution of SED by the passage of the break in an observing band (Zhang 2007). As shown in Fig. 3, the reported near-infrared colors are fairly constant from 36 min to 35 hr after the burst, and no evidence for spectral breaks other than the Ly$\alpha$ break is seen. If the 0.6% curvature in 8400–8900 Å is represented by a quadratic term in the $\lg f_{\lambda }$ versus $\lg \lambda$ relation, it would result in a change of $\Delta \beta =1$ by just a 7% change of wavelength $\lambda$. Such a strong break is not expected from the afterglow theory (Zhang 2007), and inconsistent with the agreement of $\beta$ for our fits and near-infrared colors.

Since GRB afterglows are time-variable, we also examined the time variability of the reduced spectrum, by reproducing two spectra using the first and second halves of the integrations (1–4 and 5–9 of Table 1). We confirmed that there is no statistically significant change in the spectra, and the concave curvature of the residuals for the host-only model is seen in both the spectra.

Chornock et al. (2013) also reported an analysis of the damping wing shape of the same GRB afterglow taken by another telescope, and they found no statistically significant preference for the diffuse IGM component. This result is apparently different from ours, although their upper limit ($f_{H\text{emissiontype}I}<0.11$ at $2\sigma$) is consistent with our best-fit of $f_{H\text{emissiontype}I}=0.086$ for the IGM-$z_{GRB}$ model. The origin of this difference is difficult to identify, but a clear difference of their analysis from ours is the index of unabsorbed power-law spectrum: $f_{\lambda }\propto {\lambda }^{-0.01}$, i.e., $\beta =-1.99$. Such a small value is significantly different from those measured by near-infrared colors, and it is not supported either from our $\beta$-free fitting analysis, giving an unacceptably large ${\chi }^2$ of $\sim 470$. The fit by Chornock et al. indeed shows a concave curvature of the observed spectrum compared with the fitted model (see their Fig. 2), and this curvature and their small (red) $\beta$ value could be a result of a power-law fit in a wavelength region where the IGM/DLA component is not negligible.

The expected number of such a DLA along the line of sight can be estimated as follows. A modest extrapolation of the DLA statistics up to $z\sim 5$ (Wolfe et al. 2005; Songaila & Cowie 2010) suggests the DLA number distribution per unit redshift $d{N}_{DLA}/dz\sim 0.72$ at $z=6$ for $\lg N_{H\text{emissiontype}I}\ge 20.3$. Assuming the column density distribution as $N_{DLA}(>N_{H\text{emissiontype}I})\propto N_{H\text{emissiontype}I}^{-0.85}$ again from the lower-$z$ sample, the probability of finding such a DLA in the redshift range between 5.806 and 5.913 is estimated to be $\sim 0.03$. Therefore the chance probability of finding such a DLA is rather small, though it is not negligible.

The existence of such a large column DLA at $z=5.806$ can be further constrained by the spectral profiles around Lyman series absorption at this redshift, because absorption by H\emissiontypeI in the DLA would erase the transmission of light around this redshift. The spectrum around Ly$\alpha$, $\beta$, and $\gamma$ are shown in Fig. 6, with the expected absorption factor by the DLA having $\lg N_{H\text{emissiontype}I}^{DLA}=20.704$. Here, the velocity distribution of H\emissiontypeI in the DLA is modeled as follows. The profiles of the three metal absorption lines (Si\emissiontypeII, O\emissiontypeI, C\emissiontypeII) at this redshift are also shown in Fig. 6, and their widths are consistent with the spectral resolution (3.9 Å FWHM), indicating narrow intrinsic widths. The relative velocity difference of these lines are at most $\sim 50$ km/s. The metal absorption line profile of typical DLAs is a mixture of multiple components with a velocity width of ${\sigma }_v\sim 4$–7 km/s for one component, and the relative velocity width $\Delta v$ of different components is widely distributed in $\sim$20–200 km/s (Wolfe et al. 2005; Dessauges-Zavadsky et al. 2003). Therefore, as a reasonable model, we adopt ${\sigma }_v=5.5$ km/s for each single component and assumed that these components are uniformly distributed in the relative velocity range of $\Delta v=-20$ to +20 km/s with respect to $z=5.806$.

The light transmission around the Lyman series in Fig. 6 must be consistent with the predicted attenuation factor by the DLA, and the transmissions around Ly$\alpha$ and $\gamma$ satisfy this requirement, if the spectral resolution is taken into account. However, a large transmission flux centered at 6976.97 Å is found near Ly$\beta$ with an integrated flux of $(1.6\pm 0.1)\times {10}^{-17}ergc{m}^{-2}{s}^{-1}$. This transmission profile has a width comparable with the spectral resolution, indicating a narrow intrinsic width. Furthermore, the light attenuation by the assumed DLA becomes rapidly stronger beyond 6977 Å, which also requires that the transmission is coming from a narrow ($\lesssim$ 1 Å) wavelength window. This means that the transmission flux density must be $f_{\lambda }\ge 1.6\times {10}^{-17}ergc{m}^{-2}{s}^{-1}$ Å${}^{-1}$, which should be compared with the unabsorbed continuum flux $3.6\times {10}^{-17}ergc{m}^{-2}{s}^{-1}$ Å${}^{-1}$ estimated by the continuum flux at 9800 Å and an extrapolation by using the best-fit power-law $\beta =-0.81$ of the intervening DLA model (Table 2). Here, these flux estimates of the transmission and unabsorbed continuum were made from the 300R spectrum that covers Ly$\beta$ and 9800 Å simultaneously, and hence we do not have to consider the time evolution of the flux.

This result means that the IGM transmission fraction must be larger than 1.6/3.6 = 0.44 at 6977 Å, but this is inconsistent with the transmission allowed by Ly$\beta$ absorption of the DLA that we are now considering, 0.18$\pm$0.05. Hence the particular case that the DLA at $z=5.806$ is responsible for the unexplained residuals of the host-only model is excluded. We examined possible systematic uncertainties about the relative flux calibration at 6977 and 9800 Å (by aperture loss, atmospheric dispersion, guide error, standard star, and wavelength-dependent seeing), and it should be less than 10%. A change of $\beta$ by 0.1 results in just a 3% change of the unabsorbed flux at the Ly$\beta$ wavelength. It should also be noted that the modeling about velocity distribution of H\emissiontypeI in the DLA hardly affects this argument, because the 6977 Å transmission is far from the resonant Ly$\beta$ wavelength and the attenuation is simply determined by the damping wing.

Therefore, the H\emissiontypeI column density of the $z=5.806$ system must be lower than the value required to explain the residual of the host-only model. If we assume that the $z=5.806$ system has a similar metallicity to the GRB host galaxy, a reasonable amount of the H\emissiontypeI column density would be $\lg N_{H\text{emissiontype}I}\sim 19.2$, since the metal-line equivalent widths (Si\emissiontypeII, O\emissiontypeI, and C\emissiontypeII) of the $z=5.806$ system are about 3–9 times lower than those found in the GRB host galaxy (Chornock et al. 2013). Flux attenuation factors assuming this H\emissiontypeI column density are also plotted in Fig. 6, and in this case Ly$\beta$ absorption at 6977 Å is negligible, thus allowing the existence of the observed transmission there. Such a small H\emissiontypeI column density does not affect the fitting of the Ly$\alpha$ damping wing at the GRB redshift.

Finally, we note that the simple diffuse IGM models with uniform $f_{H\text{emissiontype}I}\gtrsim 0.1$ is also mathematically excluded by the light transmission at 6977 Å, because of the large predicted Ly$\beta$ GP optical depth of ${\tau }_{\beta }=7.2\times {10}^4{f}_{H\text{emissiontype}I}$ at $z=5.802$ using the absorption oscillator strength ratio of $f_{\beta }/f_{\alpha }=0.190$. However, in reality the neutral fraction cannot be perfectly homogeneous, and many transmission spikes are expected at locally ionized bubbles along the sightlines, as seen in the spectrum of this GRB (Chornock et al. 2013). Therefore, a transmission spike at 6977 Å does not exclude any realistic models of diffuse IGM absorption.

Finally, even if the large column of $\lg N_{H\text{emissiontype}I}^{DLA}=20.704$ is associated with the $z=5.806$ system, the metallicity of this system must then be extremely low. From the reported equivalent widths of the Si\emissiontypeII, O\emissiontypeI, and C\emissiontypeII lines (Chornock et al. 2013) and this H\emissiontypeI column, we estimate the metallicity of [Si/H], [O/H], and [C/H] to be $-3.68\pm 0.10$, $-3.27\pm 0.10$, and $-3.51\pm 0.12$, respectively, in the optically-thin limit. Here we followed the usual practice for DLA studies and did not apply ionization correction, because these low-ionization atoms are the dominant phase of these elements in interstellar medium. These lines may be saturated, and hence we estimated the metal column density considering the saturation effect, assuming that these lines have a single component of ${\sigma }_v=5.5$ km/s. Then the metallicities are increased to $-{3.48}_{-0.16}^{+0.15}$, $-{2.69}_{-0.31}^{+0.44}$, and $-{3.19}_{-0.26}^{+0.30}$, respectively. (Note that the saturation effect would be weaker than this, if these lines are composed of multiple components.)

The saturation effect is not large for the weak line of Si\emissiontypeII, and [Si/H] of this system must be significantly lower than the lowest value ($\sim -2.7$) found for DLAs at lower redshifts up to $z\sim 5$ (Rafelski et al. 2012). In comparison with a theoretical model of Finlator et al. (2013), the saturation-corrected column density of $\lg N_{O\text{emissiontype}I}=14.7$ for a DLA having $\lg N_{H\text{emissiontype}I}^{DLA}=20.704$ at $z\sim 6$ is outside the $N_{O\text{emissiontype}I}$-$N_{H\text{emissiontype}I}$ distribution of 99% sightlines. We have already argued that the Ly$\beta$ profile excludes the particular case of the DLA at $z=5.806$ as an explanation of the residuals found in the host-only fit, and one may consider the possibility of a DLA at a different redshift to save the intervening DLA scenario. However, in this case the DLA metallicity must be further lower than the values estimated above, because there are no detected metal absorption lines.