A Significantly Neutral Intergalactic Medium Around the Luminous z = 7 Quasar J0252–0503

The earliest luminous quasars, powered by billion solar-mass supermassive black holes (SMBHs), can be used not only to constrain the physics of SMBH accretion and the assembly of the first generation of massive galaxies in the early universe, but also to obtain critical information on the physical conditions of the intergalactic medium (IGM) during the epoch of reionization (EoR). Although more than 200 z > 6 quasars have been found in the past few decades (e.g., Fan et al. 2001; Willott et al. 2010; Wu et al. 2015; Bañados et al. 2016; Jiang et al. 2016; Matsuoka et al. 2016; Wang et al. 2016; Reed et al. 2017), only several tens of them are at z > 6.5 (e.g., Venemans et al. 2015; Mazzucchelli et al. 2017; Wang et al. 2017, 2019; Reed et al. 2019; Yang et al. 2019) and just six are currently known at z > 7 (Mortlock et al. 2011; Bañados et al. 2018; Wang et al. 2018; Matsuoka et al. 2019a, 2019b; Yang et al. 2019). The limited number of known high-redshift quasars is due to the combination of a rapid decline of quasar spatial density toward higher redshifts (e.g., Wang et al. 2019), the lack of deep wide-field near-infrared surveys, and the presence of a large number of contaminants from Galactic cool dwarf populations in the photometric quasar selection process. Near-infrared spectroscopic observations of these known quasars indicate that billion or even ten billion solar-mass SMBHs are already in place in these luminous quasars (e.g., Wu et al. 2015; Shen et al. 2019). The existence of these SMBHs in such a young universe challenges our understanding of the formation and the growth mechanisms of SMBHs (e.g., Volonteri & Rees 2006; Pezzulli et al. 2016; Davies et al. 2019; Wise et al. 2019).

Observations of the Lyman series forests in z ≳ 6 quasars indicate that the IGM is already highly ionized by z ∼ 6 (e.g., Fan et al. 2006; Bosman et al. 2018; Eilers et al. 2018, 2019; Yang et al. 2020), although the final completion of reionization might extend down to z ∼ 5.5 (e.g., Becker et al. 2015; Davies et al. 2018a; Kulkarni et al. 2019; Keating et al. 2020; Nasir & D'Aloisio 2020). However, the Lyman series forests are very sensitive to neutral hydrogen and saturate even at low IGM neutral fraction (i.e., $\langle {x}_{{\rm{H}}{\rm{I}}}\rangle \gtrsim {10}^{-4}$). On the other hand, if the neutral fraction is of order unity, one would expect to see appreciable absorption redward of the wavelength of the Lyα emission line, resulting in a damping wing profile (e.g., Miralda-Escudé 1998) due to significant optical depth on the Lorentzian wing of the Lyα absorption. The first quasar with a damping wing detection is ULAS J1120+0641 (Mortlock et al. 2011) at z = 7.09, although different analyses yielded different constraints on $\langle {x}_{{\rm{H}}{\rm{I}}}\rangle$ (Bolton et al. 2011; Mortlock et al. 2011; Bosman & Becker 2015; Greig et al. 2017; Davies et al. 2018b), ranging from $\langle {x}_{{\rm{H}}{\rm{I}}}\rangle \sim 0$ to $\langle {x}_{{\rm{H}}{\rm{I}}}\rangle \sim 0.5$ at z ∼ 7.1. Recently, the spectrum of quasar ULAS J1342+0928 (Bañados et al. 2018) at z = 7.54 shows a robust detection of the damping wing signal (Bañados et al. 2018; Davies et al. 2018b; Greig et al. 2019; Ďurovčíková et al. 2020), yielding $\langle {x}_{{\rm{H}}{\rm{I}}}\rangle \sim 0.2\,-0.6$ at z ∼ 7.5. Compared to other probes of reionization history, such as cosmic microwave background (CMB) polarization (Planck Collaboration et al. 2018) and Lyα emission line visibility in high-redshift galaxies (e.g., Ouchi et al. 2010; Mason et al. 2018), a main advantage of IGM damping wing measurement is that it can be applied to individual quasar sight lines, thereby constraining not only the average neutral fraction, but also its scatter in different regions of the IGM. However, the damping wing experiment is only feasible at very high redshifts where the IGM is relatively neutral, and current damping wing analyses have been limited to these two sight lines due to the lack of bright quasars at z ≳ 7. Thus, it is crucial to investigate the damping wing experiment along more z > 7 quasar sight lines.

In this paper, we present the detection of strong IGM damping wing absorption along the line of sight to a luminous z = 7 quasar DES J025216.64–050331.8 (hereinafter J0252–0503; Yang et al. 2019), using new high-quality optical/near-infrared spectroscopic observations; we also use the new spectrum to measure the mass and Eddington ratio of the central SMBH. In Section 2, we describe our photometric and spectroscopic observations for J0252–0503. In Section 3, we present the luminosity, BH mass, and Eddington ratio measurements of J0252–0503. In Section 4, we discuss the reconstructions of the unabsorbed spectrum of the quasar and our constraints on the neutral fraction in the IGM at z = 7 by modeling IGM Lyα absorption. Finally, in Section 5 we summarize our results and briefly discuss the implications for the cosmic reionization history and BH growth constraints with larger quasar samples at z ≳ 7 in the future. Throughout this paper, we assume a flat ΛCDM cosmology with h = 0.685 (Betoule et al. 2014), Ω_b = 0.047, Ω_m = 0.3, Ω_Λ = 0.7, and σ₈ = 0.8. All photometry in this paper is in the AB system.

J0252–0503 (Yang et al. 2019) was selected as a quasar candidate using photometry from the Dark Energy Survey (Abbott et al. 2018) and the unblurred coadds of Wide-field Infrared Survey Explorer (Lang 2014) data. It was spectroscopically identified as a quasar at z = 7.02 based on the strong Lyα break using observations from Magellan/LDSS-3. However, the lack of a near-infrared spectrum for this quasar precluded detailed analyses in the discovery paper.

We obtained a high-quality near-infrared spectrum with the Near-Infrared Echellette Spectrometer (NIRES¹⁴ ; Wilson et al. 2004) mounted on the Keck II telescope. NIRES is a prism cross-dispersed near-infrared spectrograph with a fixed configuration that simultaneously covers the Y, J, H, and K bands in five orders from 0.94 to 2.45 μm with a small gap between 1.85 and 1.88 μm. The mean spectral resolving power of NIRES is R ∼ 2700 with a fixed 055 narrow slit. We observed J0252–0503 with NIRES for a total of 4.8 hr of on-source integration on three nights: 1.4 hr on 2018 August 12, 1.0 hr on 2018 September 3, and 2.4 hr on 2018 October 1 (UT). The observations were separated into multiple 300 s or 360 s individual exposures with the standard ABBA dither pattern. We also observed the flux standard star Feige 110. We reduced the NIRES data using a newly developed open-source Python-based spectroscopic data reduction pipeline (PypeIt¹⁵ ; Prochaska et al. 2019; Prochaska et al. 2020). Basic image processing (i.e., flat-fielding) followed standard techniques. Wavelength (in vacuum) solutions for individual frames were derived from the night sky emission lines. Sky subtractions were performed on the 2-D images by including both image differencing and a B-spline fitting procedure. We used the optimal spectrum extraction technique (Horne 1986) to extract 1-D spectra. We flux calibrated the individual 1-D spectra with the sensitivity function derived from the standard star Feige 110. We then stacked the fluxed 1-D spectra from each night and fitted a telluric absorption model directly to the stacked quasar spectra using the telluric model grids produced from the Line-By-Line Radiative Transfer Model (LBLRTM¹⁶ ; Clough et al. 2005). Finally, we combined all spectra obtained on different nights to produce the final processed 1-D spectrum.

Gemini GMOS-S (Hook et al. 2004) observations for J0252–0503 (previously described by Yang et al. 2019) were performed in two wavelength setups both with the R400 grating to cover the small wavelength gaps between detectors, with one setup centered at 860 nm and the other centered at 870 nm. These two setups yields a wavelength coverage of 0.6–1.1 μm and spectral resolution of $R\sim 1300$. Each setup was exposed for an hour. The GMOS data were also reduced with PypeIt. The spectra were flux calibrated with the sensitivity function derived from flux standard star GD71 and telluric absorption was corrected using the same method as the NIRES data reduction. In order to combine the NIRES and GMOS spectra, we scaled the NIRES co-added spectrum to the GMOS flux level using the median in the overlapping wavelength region from 9800 to 10,200 Å. The flux level of the NIRES spectrum is only ∼6% lower than that of GMOS spectrum and the shapes of these two spectra are perfectly matched. Finally, we computed the stacked spectrum in the overlap region after binning the NIRES spectrum to the GMOS wavelength grid.

Since the flux calibration is crucial for the damping wing analyses, we also obtained near-infrared Y, J, H, and K-band photometry with UKIRT/WFCam on 2018 November 27. The on-source times were 8 minutes in each band. The data were processed using the standard VISTA/WFCAM data-flow system by M. Irwin (Irwin et al. 2004). The magnitudes of J0252–0503 were measured to be Y = 20.33 ± 0.07, J = 20.19 ± 0.07, H = 20.02 ± 0.07, and K = 19.92 ± 0.08. We then scaled the combined NIRES and GMOS spectrum by carrying out synthetic photometry on the spectrum using the WFCAM J-band filter response curve to match the J-band photometry for absolute flux calibration. The magnitudes measured from the J-band scaled spectrum in the Y, H, and K bands are 20.36, 20.09, and 19.93 mag, respectively. The consistency of magnitudes derived from the fluxed spectrum and UKIRT observations indicates that the spectrophotometric calibration of the spectrum is accurate to within 10%. Finally, we corrected for Galactic extinction using the dust extinction map derived by Schlegel et al. (1998). The spectrum was then deredshifted with the systemic redshift z = 7.000 ± 0.001, derived from the IRAM NOrthern Extended Millimeter Array (NOEMA) observations of the far-infrared [C ii] emission line.¹⁷ The final spectrum used for the following analyses is shown in Figure 1. Note that in Figure 1, the spectrum is plotted after being rebinned to 200 $\mathrm{km}\,{{\rm{s}}}^{-1}$ pixels.

Zoom In Zoom Out Reset image size

In order to derive the rest-frame ultraviolet (UV) properties of J0252–0503, we fit a pseudocontinuum model which includes a power-law continuum, iron (Fe ii and Fe iii) emission (Vestergaard & Wilkes 2001; Tsuzuki et al. 2006), and Balmer continuum (e.g., De Rosa et al. 2014) to the line-free region of the calibrated and deredshifted spectrum. This pseudocontinuum model is then subtracted from the quasar spectrum, leaving a line-only spectrum. We then fit the Mg ii broad emission line in the continuum-subtracted spectrum with two Gaussian profiles. To estimate the uncertainties of our spectral measurements, we use a Monte Carlo approach (e.g., Shen et al. 2019) to create 100 mock spectra by randomly adding Gaussian noise at each pixel with its scale equal to the spectral error at that pixel. We then apply the exact same fitting algorithm to these mock spectra. The uncertainties of measured spectral properties are then estimated based on the 16% and 84% percentile deviation from the median.

The pseudocontinuum model is shown in Figure 1 and an enlargement of the Mg ii region fitting is shown in the right insert panel of Figure 1. The fitting procedure yields a power-law continuum of ${f}_{\lambda }\propto {\lambda }^{-1.67\pm 0.04}$, from which we measure the rest-frame 3000 Å luminosity to be $\lambda {L}_{3000\mathring{\rm A} }$ = (2.5 ± 0.2) × 10⁴⁶ erg s⁻¹, implying a bolometric luminosity of L_bol = 5.15$\times \,\lambda {L}_{3000\mathring{\rm A} }$ = (1.3 ± 0.1) × 10⁴⁷ erg s⁻¹ (Shen et al. 2011). The rest-frame 1450 Å magnitude is measured to be M₁₄₅₀ = −26.63 ± 0.07. The FWHM and equivalent width (EW) of the Mg ii line are measured to be FWHM${}_{\mathrm{Mg}{\rm{II}}}=3503\pm 205$ km s⁻¹ and EW_{Mg ii} = 18.83 ± 0.92 Å, respectively. The Mg ii emission line is blueshifted by ${{\rm{\Delta }}}_{v,\mathrm{Mg}{\rm{II}}}=(712\pm 50)\,\mathrm{km}\,{{\rm{s}}}^{-1}$ relative to the systemic redshift determined from the [C ii] line, similar to other luminous z ∼ 7 quasars in which Mg ii blueshifts range from a few hundred to ∼1000 km s⁻¹ (e.g., Venemans et al. 2016; Mazzucchelli et al. 2017; Bañados et al. 2018; Decarli et al. 2018).

We adopt the empirical relation obtained by Vestergaard & Osmer (2009) to estimate the black hole mass of J0252–0503, which yields ${M}_{\mathrm{BH}}=(1.39\pm 0.16)\times {10}^{9}\,{M}_{\odot }$. Note that the quoted black hole mass uncertainty does not include the systematic uncertainties of the scaling relation, which could be up to ∼0.5 dex (Shen 2013). By comparing the bolometric luminosity estimated above with the Eddington luminosity, which is ${L}_{\mathrm{Edd}}=1.3\times {10}^{38}\times {M}_{\mathrm{BH}}$, we measure the Eddington ratio of J0252–0503 to be λ_Edd = 0.7 ± 0.1. Note that the uncertainty quoted here does not consider the systematic uncertainties introduced by both single-epoch BH mass estimators and monochromatic bolometric corrections. The Eddington ratio of J0252–0503 is slightly lower than that of the other three luminous z ≥ 7 quasars: ${\lambda }_{\mathrm{Edd}}={1.5}_{-0.4}^{+0.5}$ for J1342+0928 (Bañados et al. 2018) at z = 7.54, ${\lambda }_{\mathrm{Edd}}={1.2}_{-0.5}^{+0.6}$ for J1120+0641 at z = 7.09 (Mortlock et al. 2011), and ${\lambda }_{\mathrm{Edd}}=1.25\pm 0.19$ for J0038–1527 at z = 7.02 (Wang et al. 2018). If J0252–0503 has been accreting at such Eddington ratio since z ∼ 20 with a radiative efficiency of 10%, it would require a seed BH of $\sim {10}^{5}\,{M}_{\odot }$, which significantly exceeds the predicted mass range from stellar remnant BHs and requires more exotic seed formation mechanisms like direct collapse BHs. Even if it was accreting at the Eddington limit, J0252–0503 would still require the seed BH to be more massive than $\sim {10}^{4}\,{M}_{\odot }$. This indicates that J0252–0503 is one of the few quasars that put the most stringent constraints on SMBH formation and growth mechanisms.

Among the six public known z ≥ 7 quasars, two objects already have had damping wing analyses performed (Bolton et al. 2011; Mortlock et al. 2011; Bosman & Becker 2015; Greig et al. 2017, 2019; Bañados et al. 2018; Davies et al. 2018b; Ďurovčíková et al. 2020). Two other quasars are too faint (M₁₄₅₀ ≳ − 25) for damping wing analyses with current facilities (Matsuoka et al. 2019a, 2019b), and another is a broad absorption line (BAL) quasar in which strong absorption precludes determination of the intrinsic quasar spectrum (Wang et al. 2018). Thus, J0252–0503 is the only known bright, non-BAL quasar at z ≥ 7 of which a damping wing analysis has not been performed yet. In order to examine whether the damping wing is present in the spectrum of J0252–0503, we need to know the intrinsic quasar spectrum in the Lyα region (i.e., before IGM attenuation). In the past few years, several methods have been proposed for constructing the quasar intrinsic spectra, including stacking of low-redshift quasar spectra with similar emission line properties (e.g., Mortlock et al. 2011; Simcoe et al. 2012; Bañados et al. 2018), using the principal component analysis (PCA) decomposition approach (Davies et al. 2018b, 2018c), constructing the covariant relationships between parameters of Gaussian fits to Lyα line and those of Gaussian fits to other broad emission lines (Greig et al. 2017, 2019), and using the neural network method (Ďurovčíková et al. 2020). In this paper, we adopt both the empirical composite method and the PCA method to construct the intrinsic spectrum for J0252–0503 as detailed below.

4.1. Empirical Composite Spectra from Analogs

Since there is a lack of spectral evolution of quasars from low redshifts to high redshifts (e.g., Shen et al. 2019), the large sample of SDSS/BOSS quasars at lower redshifts provides a good training set for constructing a high-redshift quasar intrinsic spectrum. First, we use a composite spectrum constructed from a sample of low-redshift quasar analogs to model the intrinsic spectrum. Because the C iv line properties, and especially the line's blueshift, appear to be strongly connected with differences in the quasar spectral energy distribution (e.g., Richards et al. 2011), we select quasar analogs from SDSS/BOSS DR14 quasar catalog (Pâris et al. 2018) by matching the C iv blueshifts to J0252–0503. As most SDSS/BOSS quasars do not have [C ii] redshifts, we measure the relative blueshifts between the C iv and Mg ii lines. This limits us to selecting quasars in the redshift range 2.0 < z < 2.5 in order to get Lyα, C iv, and Mg ii line properties from BOSS spectra. We also excluded quasars marked as BAL and those without Mg ii redshift measurements in the catalog. This yields 85,535 quasars in total.

Before measuring the line properties from these quasars, we first fit a power-law continuum to the quasar spectrum and subtract it from the data. Instead of fitting the C iv and Mg ii lines directly, we use a more robust nonparametric scheme proposed by Coatman et al. (2016) to measure the line centroids of C iv and Mg ii lines from the continuum-subtracted spectra. The relative blueshift between these two lines is then defined as

$$\begin{eqnarray}&&{\rm{\Delta }}v=c\times \left(\displaystyle \frac{1549.06-{\lambda }_{\mathrm{half},{\rm{C}}{\rm{IV}}}}{1549.06}-\displaystyle \frac{2798.75-{\lambda }_{\mathrm{half},\mathrm{Mg}{\rm{II}}}}{2798.75}\right),\end{eqnarray} \tag{ 1 }$$

where c is the speed of light and ${\lambda }_{\mathrm{half},{\rm{C}}{\rm{IV}}}$ (${\lambda }_{\mathrm{half},\mathrm{Mg}{\rm{II}}}$) is the rest-frame wavelength that bisects the cumulative total line flux of C iv (Mg ii). We applied this procedure to the spectra of both J0252–0503 and the 85,535 SDSS/BOSS quasars. The blueshift in J0252–0503 is measured to be 4090 km s⁻¹. We then select quasars with blueshifts between 3000 and 5000 km s⁻¹ and mean spectral signal-to-noise ratios (S/Ns) per pixel in the C iv and Mg ii regions greater than 4 and 2, respectively. These S/N limits were chosen to yield enough sight lines to compute a composite. After this, we visually inspected the continuum normalized spectra and removed quasars that have BAL features, proximate damped Lyα systems (PDLAs) and strong intervening absorbers on top of the emission lines. We also reject objects with Mg ii line measurements that are strongly affected by sky line residuals and remove targets that have strongly different C iv and Mg ii line profiles than J0252–0503 (objects were removed if the line peaks differ by more than three times the spectrum error vector of J0252–0503). In the end, our master quasar analog sample consists of 83 SDSS/BOSS quasars.

Before constructing the composite spectrum, each spectrum was divided by its best-fit power-law continuum. Each spectrum was weighted by the average S/N of that spectrum when computing the composite. Then we multiplied the power-law fit from J0252–0503 with the constructed continuum normalized composite, obtaining the composite spectrum shown in Figure 1. In order to understand the uncertainties of the composite spectrum and minimize the bias introduced by visual checks, we resampled our parent sample with bootstrapping to construct another 100 composites, which are shown as thin orange lines in the insert panel of Figure 1. Overall, the constructed composite matches the J0252–0503 spectrum very well across the whole spectral range, except for the Lyα line region. From the left inset of Figure 1, we can clearly see that these composites have higher fluxes redward of the Lyα emission line (from 1216 to 1250 Å in rest frame) than J0252–0503, indicating strong absorption in the spectrum of J0252–0503.

4.2. PCA

Strong correlations between various broad emission lines of quasars from the rest-frame UV to the optical are known to exist (e.g., Richards et al. 2011). Taking this into account, in principle one can predict the shape of the Lyα line based on the properties of other broad emission lines. Davies et al. (2018c) developed a PCA predictive approach based on a training set of ∼13,000 quasar spectra from SDSS/BOSS quasar catalog (Pâris et al. 2017) to predict the "blue-side" (rest frame 1175–1280 Å) quasar spectrum from the "red-side" (rest frame 1280–2850 Å) spectrum. In brief, we performed a PCA decomposition of the training set truncated at 10 red-side and 6 blue-side basis spectra for each quasar. Then we derived a projection matrix relating the best-fit coefficients in the red-side and a template redshift to the coefficients in the blue-side (Suzuki et al. 2005; Pâris et al. 2011). With this matrix, we can then predict the blue-side coefficients and thus the blue-side spectrum from a fit to the red-side coefficients and template redshift of a given quasar spectrum.

We quantify the uncertainties of this prediction by testing the full predictive procedure on every quasar in the training set and computing their relative continuum error (see Davies et al. 2018c, for more details). We assume a multivariate Gaussian distribution for the relative continuum error, with the covariance matrix determined from the prediction errors measured for similar quasars, i.e., the 1% nearest neighbors, as the uncertainties of the prediction.

The advantage of this PCA method compared to the composite spectrum discussed in Section 4.1 is that the PCA approach takes into account the properties of all broad emission lines in the red-side rather than just the properties of the C iv line. In addition, we can quantify uncertainties in the blue-side spectrum predictions by testing the method on the input training set.

In the upper panel of Figure 2, we show the red-side PCA fit and blue-side prediction for J0252–0503 on top of the GMOS+NIRES quasar spectrum. In the bottom-left panel of Figure 2, we show a zoom-in of the Lyα region overlaid with both the blue-side PCA model and the composite spectrum constructed in Section 4.1. From this zoomed-in plot, we can see that the intrinsic quasar spectrum predicted by the PCA model agrees very well with the composite spectrum. Both models suggest that there is a strong damping wing absorption imprinted on the Lyα emission line of the quasar. Since these two models are consistent with each other, we will only use the PCA continuum model for the following analyses so that we can make use of its well quantified uncertainties.

Zoom In Zoom Out Reset image size

4.3. Modeling the Damping Wing as a Single DLA

The smooth damped absorption profile can be imprinted by either an intervening high-column-density gravitationally-bounded DLA system (${N}_{{\rm{H}}{\rm{I}}}\gt {10}^{20}\,{\mathrm{cm}}^{-2}$) or substantially neutral gas in the IGM. However, DLA systems in the quasar vicinity are very rare at high redshifts. Among more than 250 known z ≳ 5.7 quasars, only a few of them have been identified to be associated with such absorbers close to the quasar redshifts (e.g., D'Odorico et al. 2018; Bañados et al. 2019; Farina et al. 2019; Davies 2020), suggesting that the probability of the strong redward absorption seen in J0252–0503 being caused by a DLA is low. DLA systems are usually associated with a number of metal lines such as Si ii λ1260, λ1304, λ1526, O i λ1302, C ii λ1334, C iv λ1548, λ1550, Mg ii λ2796, λ2803, and a series of Fe ii lines. Thus, one way to distinguish a DLA damping wing from an IGM damping wing is to search for associated metal absorption features.

First, we need to determine the redshift of a potential DLA system. To do so, we fit a Voigt profile to the transmission spectrum, which is normalized by the PCA continuum model. Since the Doppler parameter, b, does not strongly affect the Lyα profile (e.g., Crighton et al. 2015), we fixed the b value to be b = 10 km s⁻¹ and use the Markov Chain Monte Carlo sampler (emcee; Foreman-Mackey et al. 2013) to jointly fit the redshift and H i column density of a DLA model. During the fit, we masked the narrow absorption at v ∼ 0 km s⁻¹. This absorption could be caused by neutral gas inflow since we did not find any associated metal absorption from the quasar spectrum, and it is located at a slightly higher redshift than the quasar if it is caused by neutral hydrogen. The best-fit parameters for the system are determined to be ${N}_{{\rm{H}}{\rm{I}}}\,={10}^{21.04\pm 0.04}\,{\mathrm{cm}}^{-2}$ and z_DLA = 6.939 ± 0.002. In order to qualify the uncertainties of these parameters caused by the continuum model, we then fit DLA models to 100 transmission spectra normalized by the 100 PCA draws shown in the bottom-left panel of Figure 2. The median values and the mean deviation of 16% and 84% percentiles from the median for N_{H i} and z_DLA are measured to be ${N}_{{\rm{H}}{\rm{I}}}={10}^{21.04\pm 0.10}\,{\mathrm{cm}}^{-2}$ and z_DLA = 6.940 ± 0.003. To take both the fitting uncertainty and the PCA continuum uncertainty into account, we take ${N}_{{\rm{H}}{\rm{I}}}={10}^{21.04\pm 0.14}\,{\mathrm{cm}}^{-2}$ and z_DLA = 6.940 ± 0.004 as our fiducial parameters for the DLA model, where the uncertainties are the sum of the uncertainties from the emceefitting on the transmission spectrum and the distribution of the 100 draws. This potential DLA system (if it exists) is ∼2200 km s⁻¹ away from the quasar systemic redshift, which seems unlikely to be associated with the quasar host galaxy. This best-fit DLA model is shown as the yellow dashed line in the bottom-right panel of Figure 2. We caution that the resolution of our spectrum is low in the Lyα region, so the DLA fitting procedure might overestimate the N_{H i} if there are some narrow Lyα transmission spikes in the quasar proximity zone that are unresolved in our spectrum.

We then searched for metal absorption lines at z ∼ 6.94 in the J0252–0503 spectrum. In the end, we did not find any evidence for metal-line absorption at redshifts close to the potential DLA system within Δz ∼ ±0.04, or ∼ ± 1500 km s⁻¹, 10 times wider than the redshift uncertainty of the potential DLA system. We also did not find any metal absorption features at the quasar systemic redshift (i.e., from the quasar host galaxy). We then calculated rest-frame EW 1σ limits for each expected metal absorption line as follows: ${W}_{{\rm{r}},\mathrm{Si}{\rm{II}}1260}\leqslant 0.029\,\mathring{\rm A}$, ${W}_{{\rm{r}},{\rm{O}}{\rm{I}},1302}\,\leqslant 0.024\,\mathring{\rm A}$, ${W}_{{\rm{r}},{\rm{C}}{\rm{II}}1334}\leqslant 0.025\,\mathring{\rm A}$, ${W}_{{\rm{r}},{\rm{C}}{\rm{IV}}1548}\leqslant 0.019\,\mathring{\rm A}$, ${W}_{{\rm{r}},\mathrm{Fe}{\rm{II}}2586}\,\leqslant 0.067\,\mathring{\rm A}$, ${W}_{{\rm{r}},\mathrm{Fe}{\rm{II}}2600}\leqslant 0.049\,\mathring{\rm A}$, ${W}_{{\rm{r}},\mathrm{Mg}{\rm{II}}2796}\leqslant 0.040\,\mathring{\rm A}$. The EW limits were measured by summing over the normalized pixels over an aperture spanning $\pm 2{\sigma }_{\mathrm{inst}}$ from the center of each line, where ${\sigma }_{\mathrm{inst}}=47\,\mathrm{km}\,{{\rm{s}}}^{-1}$ was derived from the NIRES instrumental resolution. In order to derive the column densities for the selected iron, we carried out a curve of growth analysis for four different b-parameters following Simcoe et al. (2012). The curve of growth analysis is shown in Figure 3. Based on the solar abundance (Lodders 2003) and the column densities derived by fixing $b=10\,\mathrm{km}\,{{\rm{s}}}^{-1}$, we find that the metallicity of the potential DLA system is most tightly constrained by C iv. However, whether high-redshift DLAs exhibit C iv is still debated (e.g., D'Odorico et al. 2018; Cooper et al. 2019). Thus we use Mg ii which sets the second most stringent constraint on the DLA abundance with $[\mathrm{Mg}/{\rm{H}}]\lt -4.0$ (3σ). The DLA abundance 3σ limits are estimated to be $[\mathrm{Si}/{\rm{H}}]\lt -3.6$, $[{\rm{O}}/{\rm{H}}]\lt -3.6$, $[{\rm{C}}/{\rm{H}}]\lt -3.7$, and $[\mathrm{Fe}/{\rm{H}}]\lt -3.3$ based on Si ii λ1260, O i λ1302, C ii λ1334, and Fe ii λ 2600, respectively. Since the b-parameter could be as low as $b=8\,\mathrm{km}\,{{\rm{s}}}^{-1}$ at high redshifts (D'Odorico et al. 2018), we also estimate the [Mg/H] based on b = 5 km s⁻¹ and find that $[\mathrm{Mg}/{\rm{H}}]\lt -3.7$ (3σ).

Zoom In Zoom Out Reset image size

To further investigate the properties of a possible DLA, we compute the composite stack of the heavy-element transitions shown in Figure 3 by assuming that there is a metal-poor DLA at z_DLA = 6.94. We stacked the transmitted flux at the expected wavelength using an inverse-variance weighted mean. The composite stack of metal lines is shown in Figure 4, which shows no significant absorption within Δv ∼ 1500 km s⁻¹. Note that the absorption feature at v ∼ 1450 km s⁻¹ in the stack is caused by the Fe ii 2344 transition from a z = 3.5425 absorber (see below). The 1σ limit for the dimensionless EW, W = W_λ / λ (Draine 2011), for the stack is measured to be W ≤ 7.3 × 10⁻⁶ (1σ). This corresponds to a limit of $[{\rm{O}}/{\rm{H}}]\lt -4.1$ (3σ) after scaling it to the cross section and relative abundance of O i. We also compute a set of DLA models by adapting b = 10 km s⁻¹ and solar abundance pattern (Lodders 2003) with varying metallicities. The DLA transmission spectra are computed in the same wavelength grid and same resolution as the spectrum of J0252–0503. The composite stack of these DLA models for different metallicities is also over-plotted in Figure 4. The composite stack with $[{\rm{Z}}/{\rm{H}}]\lt -4.3$ matches the observed stack at 3σ level, consistent with our curve of growth analysis of the observed composite metal transitions within 0.2 dex. From Figure 3, we note that most of the metal transitions are in the linear region of the curve of growth unless $b\lesssim 5\,\mathrm{km}\,{{\rm{s}}}^{-1}$. Thus the metallicity constraint does not change too much by varying the b-parameter. By varying b from 5 to $20\,\mathrm{km}\,{{\rm{s}}}^{-1}$, we can constrain the metallicity of the potential DLA system to be $[{\rm{Z}}/{\rm{H}}]\lt -4.5\sim -4.0$. Our analysis indicates that this potential DLA system would be among the most metal-poor DLA systems known (e.g., Cooke et al. 2011; Bañados et al. 2019). This suggests that the strong damped absorption is very unlikely to be caused by a DLA.

Zoom In Zoom Out Reset image size

In addition, we also searched for absorbers at lower redshifts to make sure that the damped Lyα absorption is not contaminated by lower redshift absorbers. We identify five strong Mg ii absorption systems at z = 4.8793, z = 4.7144, z = 4.2095, z = 4.0338, and z = 3.5425. These systems also exhibit associated Fe ii lines. The z = 4.8793 system also has associated Al ii λ1670 absorption line which falls into the damped absorption region which is masked in the following damping wing analysis. However, this line is very narrow (see the bottom-right panel of Figure 2) and thus would not be responsible for the smoothed absorption profile on much larger scales. These analyses indicate that the damping wing absorption in the J0252–0503 spectrum is more likely to be imprinted by the neutral IGM rather than by a DLA system or other intervening absorbers, especially considering the fact that J1120+0641 and J1342+0928 also have similar (though slightly weaker) absorption profiles that are not associated with metals (e.g., Simcoe et al. 2012; Bañados et al. 2018).

4.4. Constraints on the IGM Neutral Fraction from a Strong Damping Wing at z = 7

In order to quantitatively assess the damping wing strength and constrain the volume-averaged neutral hydrogen fraction at z = 7, we applied the methodology from Davies et al. (2018b) to this quasar sight-line. We refer the reader to Davies et al. (2018b) for a detailed description. In brief, we model the reionization-era quasar transmission spectrum with a multiscale hybrid model. This model combines large-scale seminumerical reionization simulations around massive dark matter halos computed in a (400 Mpc)³ volume with a modified version of 21cmFAST (Mesinger et al. 2011; F. Davies & S. Furlanetto 2020, in preparation), density, velocity, and temperature fields of 1200 hydrodynamical simulation skewers from a separate (100 Mpc/h)³ Nyx hydrodynamical simulation (Almgren et al. 2013; Lukić et al. 2015), and 1D ionizing radiative transfer, which models the ionization and heating of the IGM by the quasar (Davies et al. 2016). We then construct realistic forward modeled representations of quasar transmission spectra after accounting for the covariant intrinsic quasar continuum uncertainty from the PCA training. Finally, we use a Bayesian statistical method to recover the joint posterior probability distribution functions (PDFs) of $\langle {x}_{{\rm{H}}{\rm{I}}}\rangle$ based on these mock transmission spectra.

The damping wing strength not only depends on the $\langle {x}_{{\rm{H}}{\rm{I}}}\rangle$, but also strongly depends on the quasar lifetime, t_Q, due to the ionization of pre-existing neutral hydrogen along the line of sight by the quasar. In order to measure $\langle {x}_{{\rm{H}}{\rm{I}}}\rangle$, we conservatively explore a very broad t_Q range with a flat log-uniform t_Q prior covering ${10}^{3}\mathrm{yr}\lt {t}_{{\rm{Q}}}\lt {10}^{8}\mathrm{yr}$. We then compute the posterior PDF for $\langle {x}_{{\rm{H}}{\rm{I}}}\rangle$ by marginalizing over the entire model grid of t_Q, which is shown in Figure 5. The peak of the PDF leans to the high $\langle {x}_{{\rm{H}}{\rm{I}}}\rangle$ end. This is consistent with what we have seen in Figure 2, where we show a quasar transmission spectrum model within a $\langle {x}_{{\rm{H}}{\rm{I}}}\rangle =1.0$ IGM with a quasar lifetime of t_Q = 10^6.3 yr. The median and the central 68% (95%) confidence interval for $\langle {x}_{{\rm{H}}{\rm{I}}}\rangle$ are estimated to be $\langle {x}_{{\rm{H}}{\rm{I}}}\rangle ={0.70}_{-0.23}^{+0.20}{(}_{-0.48}^{+0.28})$ from the posterior PDF. As a comparison, we also show the PDFs from the other two z > 7 quasar sight lines in Figure 5. Although the redshift of J0252–0503 is lower than the other two quasars, the damping wing in J0252–0503 is the strongest one.

Zoom In Zoom Out Reset image size

In Figure 6, we plot the $\langle {x}_{{\rm{H}}{\rm{I}}}\rangle$ constraints from all three quasar damping wings. In this figure, we also show the $\langle {x}_{{\rm{H}}{\rm{I}}}\rangle$ constraints from the Lyα+Lyβ forest (Fan et al. 2006), as well as Lyα+Lyβ dark gaps (McGreer et al. 2015). All three z ≥ 7 quasars for which a damping wing analysis can be done with current facilities and methodology show evidence of damping wing absorptions, suggesting that the IGM is substantially neutral at z ≥ 7. These constraints are consistent with the integral constraints of $\langle {x}_{{\rm{H}}{\rm{I}}}\rangle$ measured from the electron scattering optical depth of the CMB (Planck Collaboration et al. 2018) shown as the underlying shaded region in Figure 6. They are also in broad agreement with recent calculations (e.g., Robertson et al. 2015; Naidu et al. 2020) and simulations (e.g., Kulkarni et al. 2019) of the cosmic reionization history, as well as constrains from gamma-ray burst damping wings (Totani et al. 2006, 2016; Greiner et al. 2009), the detections of Lyα emissions from high-redshift galaxies (e.g., Ouchi et al. 2010; Mason et al. 2018), and Lyα luminosity functions (e.g., Kashikawa et al. 2006; Konno et al. 2018).

Zoom In Zoom Out Reset image size

In this paper we present high-quality near-infrared spectroscopic observations of a bright z = 7 quasar, J0252–0503, to constrain the cosmic reionization with quasar damping wing modeling and the SMBH growth with BH mass and Eddington ratio measurements.

We measure the mass of the central SMBH to be ${M}_{\mathrm{BH}}=(1.39\pm 0.16)\times {10}^{9}\,{M}_{\odot }$ based on the single-epoch virial method. The Eddington ratio of J0252–0503 is measured to be λ_Edd = 0.7 ± 0.1, slightly lower than that of the other three z ≥ 7 quasars with similar luminosities. If J0252–0503 has been accreting at such Eddington ratio since z ∼ 20 with a radiative efficiency of 10%, it would require a seed BH of $\sim {10}^{5}\,{M}_{\odot }$, which significantly exceeds the predicted mass range from stellar remnant BHs and requires more exotic seed formation mechanisms like direct collapse BHs. J0252–0503, along with the other three luminous z > 7 quasars hosting billion solar-mass SMBHs, places the strongest constraints on early BH assembly mechanisms.

In order to investigate whether a damping wing is present in the spectrum of J0252–0503, we explored two different methods to construct the intrinsic spectrum of J0252–0503. The Lyα region of a composite spectrum computed from a sample of C iv blueshift-matched low-redshift quasar analogs is consistent with the prediction made by a PCA nonparametric predictive approach. Both methods suggest that a strong damping wing absorption is present in the J0252–0503 spectrum. We modeled the damping wing profile produced by either a single-component DLA system or a significantly neutral IGM. However, there is no significant detection of metals at the potential DLA system redshift over a wide range of ±1500 km s⁻¹, suggesting that the strong damping wing in the J0252–0503 spectrum is most likely imprinted by a significantly neutral IGM unless the metallicity of the putative DLA is more than 10,000 times lower than the solar metallicity.

To constrain the IGM neutral hydrogen fraction, $\langle {x}_{{\rm{H}}{\rm{I}}}\rangle$, at z = 7 with the damping wing in J0252–0503, we applied the hybrid model developed by Davies et al. (2018b) to our PCA continuum prediction for J0252–0503. Our analysis shows that the damping wing in J0252–0503 is the strongest one yet seen in z ≥ 7 quasar spectra. By marginalizing over quasar lifetime with a log-uniform prior in the range of ${10}^{3}\lt {t}_{{\rm{Q}}}\lt {10}^{8}\,\mathrm{yr}$, we measure the median and the central 68% (95%) confidence interval for $\langle {x}_{{\rm{H}}{\rm{I}}}\rangle$ to be $\langle {x}_{{\rm{H}}{\rm{I}}}\rangle ={0.70}_{-0.23}^{+0.20}{(}_{-0.48}^{+0.28})$ at z ∼ 7. The recent study by D'Aloisio et al. (2020) suggests that unrelaxed gaseous structures may exist in the postreionization IGM, meaning that the mean free path of ionizing photons is shorter compared with a model that assumes the gas is fully relaxed. The mean free path in the quasar proximity zone, however, should still be quite long due to the strong ionizing radiation of the central luminous quasar (McQuinn et al. 2011; D'Aloisio et al. 2018; Davies 2020). Thus our constraints on $\langle {x}_{{\rm{H}}{\rm{I}}}\rangle$ based on damping wing analysis should not be strongly affected by unrelaxed baryons in the proximity zone.

Despite the limited precision of quasar continuum reconstructions and the degeneracy of $\langle {x}_{{\rm{H}}{\rm{I}}}\rangle$ and quasar lifetime, the damping wing is still highly effective in constraining the reionization history. Although the currently available sample of quasar sight lines at z ≳ 7 is very small, more luminous z ≳ 7 quasars are expected to be found in the next few years through ongoing quasar searches (e.g., Bañados et al. 2018; Wang et al. 2018; Matsuoka et al. 2019b; Reed et al. 2019; Yang et al. 2019). Moreover, the Euclid wide survey will be online soon, and will discover more than 100 quasars at z > 7 (Euclid Collaboration et al. 2019). In addition, the Near-Infrared Spectrograph on the James Webb Space Telescope will provide much higher quality spectroscopic data for more precise quasar damping wing analyses. Thus, we expect that quasar damping wing analyses will have the capability to place increasingly strong constraints on the cosmic reionization history during the next several years.

Support for this work was provided by NASA through the NASA Hubble Fellowship grant #HST-HF2-51448.001-A awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Incorporated, under NASA contract NAS5-26555. J.Y., X.F., and M.Y. acknowledge support from the US NSF grant AST-1515115 and NASA ADAP grant NNX17AF28G. Research by A.J.B. is supported by NSF grant AST-1907290. X.-B.W. and L.J. acknowledge support from the National Key R&D Program of China (2016YFA0400703) and the National Science Foundation of China (11533001 & 11721303). A.C.E. acknowledges support by NASA through the NASA Hubble Fellowship grant #HST-HF2-51434 awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA, under contract NAS5-26555. B.P.V. and F.W. acknowledge funding through the ERC grant "Cosmic Gas."

The data presented in this paper were obtained at the W.M. Keck Observatory, which is operated as a scientific partnership among the California Institute of Technology, the University of California and the National Aeronautics and Space Administration. The Observatory was made possible by the generous financial support of the W.M. Keck Foundation. This work was supported by a NASA Keck PI Data Award, administered by the NASA Exoplanet Science Institute. Data presented herein were obtained at the W. M. Keck Observatory from telescope time allocated to the National Aeronautics and Space Administration through the agency's scientific partnership with the California Institute of Technology and the University of California. This research based on observations obtained at the Gemini Observatory (GS-2018B-FT-202), which is operated by the Association of Universities for Research in Astronomy, Inc., under a cooperative agreement with the NSF on behalf of the Gemini partnership: the National Science Foundation (United States), National Research Council (Canada), CONICYT (Chile), Ministerio de Ciencia, Tecnología e Innovación Productiva (Argentina), Ministério da Ciência, Tecnologia e Inovação (Brazil), and Korea Astronomy and Space Science Institute (Republic of Korea). UKIRT is owned by the University of Hawaii (UH) and operated by the UH Institute for Astronomy; operations are enabled through the cooperation of the East Asian Observatory.

The authors thank Percy Gomez and Greg Doppmann for their expert support and advice during our NIRES observing runs. Some of the data presented herein were obtained using the UC Irvine Remote Observing Facility, made possible by a generous gift from John and Ruth Ann Evans.

The authors wish to recognize and acknowledge the very significant cultural role and reverence that the summit of Maunakea has always had within the indigenous Hawaiian community. We are most fortunate to have the opportunity to conduct observations from this mountain.

Facilities: Gemini(GMOS) - , Keck(NIRES) - , UKIRT(WFCam) - , NOEMA - .

Software: astropy (Astropy Collaboration et al. 2013), emcee (Foreman-Mackey et al. 2013), matplotlib (Hunter 2007), PypeIt (Prochaska et al. 2019).