Implications of $z \sim 6$ Quasar Proximity Zones for the Epoch of Reionization and Quasar Lifetimes

TL;DR

This study analyzes quasar proximity zones at high redshift, finding a shallower evolution than previously thought and good agreement with simulations, while identifying outliers possibly due to short quasar lifetimes or local absorption.

Contribution

It provides a homogeneous analysis of 34 high-redshift quasars, updates systemic redshift measurements, and compares observations with simulations revealing new insights into reionization and quasar lifetimes.

Findings

Proximity zone size evolution is shallower than previously reported.

Good agreement between observed and simulated proximity zones regardless of IGM ionization.

Identification of outliers with very small proximity zones potentially indicating short quasar lifetimes.

Abstract

We study quasar proximity zones in the redshift range $5.77 \leq z \leq 6.54$ by homogeneously analyzing $34$ medium resolution spectra, encompassing both archival and newly obtained data, and exploiting recently updated systemic redshift and magnitude measurements. Whereas previous studies found strong evolution of proximity zone sizes with redshift, and argued that this provides evidence for a rapidly evolving intergalactic medium (IGM) neutral fraction during reionization, we measure a much shallower trend $\propto(1+z)^{-1.44}$. We compare our measured proximity zone sizes to predictions from hydrodynamical simulations post-processed with one-dimensional radiative transfer, and find good agreement between observations and theory irrespective of the ionization state of the ambient IGM. This insensitivity to IGM ionization state has been previously noted, and results from the fact that the definition of proximity zone size as the first drop of the smoothed quasar spectrum below the $10\%$ flux transmission level probes locations where the ionizing radiation from the quasar is an order of magnitude larger than the expected ultraviolet ionizing background that sets the neutral fraction of the IGM. Our analysis also uncovered three objects with exceptionally small proximity zones (two have $R_p < 1$proper Mpc), which constitute outliers from the observed distribution and are challenging to explain with our radiative transfer simulations. We consider various explanations for their origin, such as strong absorption line systems associated with the quasar or patchy reionization, but find that the most compelling scenario is that these quasars have been shining for $\lesssim 10^5$yr.

Tables

Table 1. Table 1 Overview of the observations of the 34 34 34 quasars in our data sample.

SDSS J0002+2550	${00}^h{02}^{m}39\stackrel{s}{.}39$	$+25°50\prime 34\stackrel{\prime \prime }{.}96$	Kakazu	Nov. 2004	$5400$ s
			Cowie	Aug. 2005	$16300$ s
SDSS J0005-0006	${00}^h{05}^{m}52\stackrel{s}{.}34$	$-00°06\prime 55\stackrel{\prime \prime }{.}80$	Becker	Dec. 2002	$1200$ s
			Sargent	Oct. 2010	$15000$ s
CFHQS J0050+3445	${00}^h{55}^{m}02\stackrel{s}{.}91$	$+34°45\prime 21\stackrel{\prime \prime }{.}65$	Willott	Sep. 2008	$6250$ s
			Sargent	Oct. 2010	$6000$ s
			Sargent	Oct. 2012	$9000$ s
SDSS J0100+2802	${01}^h{00}^{m}13\stackrel{s}{.}02$	$+28°02\prime 25\stackrel{\prime \prime }{.}92$	White	Jan. 2016	$3600$ s
ULAS J0148+0600	${01}^h{48}^{m}37\stackrel{s}{.}64$	$+06°00\prime 20\stackrel{\prime \prime }{.}06$	Sargent	Oct. 2010	$11200$ s
ULAS J0203+0012	${02}^h{03}^{m}32\stackrel{s}{.}38$	$+00°12\prime 29\stackrel{\prime \prime }{.}27$	Sargent	Oct. 2010	$6600$ s
CFHQS J0210-0456	${02}^h{10}^{m}13\stackrel{s}{.}19$	$-04°56\prime 20\stackrel{\prime \prime }{.}90$	Sargent	Oct. 2010	$6000$ s
PSO J0226+0302	${02}^h{26}^{m}01\stackrel{s}{.}87$	$+03°02\prime 59\stackrel{\prime \prime }{.}42$	White	Jan. 2016	$11700$ s
CFHQS J0227-0605	${02}^h{27}^{m}43\stackrel{s}{.}29$	$-06°05\prime 30\stackrel{\prime \prime }{.}20$	Willott	Sep. 2008	$10540$ s
SDSS J0303-0019	${03}^h{03}^{m}31\stackrel{s}{.}40$	$-00°19\prime 12\stackrel{\prime \prime }{.}90$	Sargent	Oct. 2010	$6000$ s
SDSS J0353+0104	${03}^h{53}^{m}49\stackrel{s}{.}73$	$+01°04\prime 04\stackrel{\prime \prime }{.}66$	Sargent	Oct. 2010	$13200$ s
			Becker	Jan. 2006	$3600$ s
PSO J0402+2452	${04}^h{02}^{m}12\stackrel{s}{.}69$	$+24°51\prime 24\stackrel{\prime \prime }{.}43$	White	Jan. 2016	$10800$ s
SDSS J0818+1723	${08}^h{18}^{m}27\stackrel{s}{.}40$	$+17°22\prime 52\stackrel{\prime \prime }{.}01$	Becker	Apr. 2005	$2400$ s
SDSS J0836+0054	${08}^h{36}^{m}43\stackrel{s}{.}86$	$+00°54\prime 53\stackrel{\prime \prime }{.}26$	Becker	Mar. 2001	$1731$ s
			Madau	Feb. 2002	$18900$ s
			Cowie	Feb. 2002	$13540$ s
			Djorgovski	Mar. 2002	$8923$ s
			Kulkarni	Jan. 2003	$5400$ s
			Cowie	Jan. 2003	$3600$ s
			Cowie	Jan. 2004	$10800$ s
			Kakazu	Nov. 2004	$10200$ s
			Djorgovski	Dec. 2004	$7200$ s
SDSS J0840+5624	${08}^h{40}^{m}35\stackrel{s}{.}30$	$+56°24\prime 20\stackrel{\prime \prime }{.}22$	Djorgovski	Dec. 2004	$10800$ s
			Becker	Mar. 2006	$1200$ s
SDSS J0842+1218	${08}^h{42}^{m}29\stackrel{s}{.}43$	$+12°18\prime 50\stackrel{\prime \prime }{.}58$	Becker	Mar. 2006	$2400$ s
SDSS J0927+2001	${09}^h{27}^{m}21\stackrel{s}{.}82$	$+20°01\prime 23\stackrel{\prime \prime }{.}64$	Becker	Mar. 2006	$1200$ s
SDSS J1030+0524	${10}^h{30}^{m}27\stackrel{s}{.}11$	$+05°24\prime 55\stackrel{\prime \prime }{.}06$	Becker	May 2001	$1800$ s
			Becker	Jan. 2002	$12000$ s
			Madau	Feb. 2002	$7339$ s
			Cowie	Feb. 2002	$16200$ s
			Djorgovski	Mar. 2002	$16200$ s
			Cowie	Jan. 2003	$5400$ s
SDSS J1048+4637	${10}^h{48}^{m}45\stackrel{s}{.}07$	$+46°37\prime 18\stackrel{\prime \prime }{.}55$	Becker	Dec. 2002	$3600$ s
			Cowie	Jan. 2003	$18000$ s
			Cowie	Mar. 2003	$5400$ s
			Sanchez	Jun. 2003	$5400$ s
			Djorgovski	Dec. 2004	$6600$ s
SDSS J1137+3549	${11}^h{37}^{m}17\stackrel{s}{.}73$	$+35°49\prime 56\stackrel{\prime \prime }{.}85$	White	Jan. 2016	$7800$ s
			Becker	Jan. 2005	$2400$ s
SDSS J1148+5251	${11}^h{48}^{m}16\stackrel{s}{.}65$	$+52°51\prime 50\stackrel{\prime \prime }{.}39$	Becker	May 2002	$8100$ s
			Becker	Dec. 2002	$22800$ s
			Djorgovski	Dec. 2002	$11400$ s
			Kulkarni	Jan. 2003	$3600$ s
			Cowie	Jan. 2003	$16200$ s
			Becker	Feb. 2003	$16800$ s
			Cowie	Jan. 2004	$11700$ s
SDSS J1250+3130	${12}^h{50}^{m}51\stackrel{s}{.}93$	$+31°30\prime 21\stackrel{\prime \prime }{.}90$	Becker	Jan. 2005	$3600$ s
SDSS J1306+0359	${13}^h{06}^{m}08\stackrel{s}{.}27$	$+03°59\prime 26\stackrel{\prime \prime }{.}36$	Becker	May 2001	$900$ s
			Cowie	Feb. 2002	$15300$ s
			Djorgovski	Mar. 2002	$12600$ s
			Cowie	Mar. 2003	$9000$ s
			Cowie	Mar. 2004	$12600$ s
ULAS J1319+0950	${13}^h{19}^{m}11\stackrel{s}{.}30$	$+09°50\prime 51\stackrel{\prime \prime }{.}52$	Steidel	Mar. 2008	$2400$ s
SDSS J1335+3533	${13}^h{35}^{m}50\stackrel{s}{.}81$	$+35°33\prime 15\stackrel{\prime \prime }{.}82$	Becker	Mar. 2006	$1200$ s
SDSS J1411+1217	${14}^h{11}^{m}11\stackrel{s}{.}29$	$+12°17\prime 37\stackrel{\prime \prime }{.}28$	Djorgovski	Apr. 2004	$10200$ s
			Becker	Jan. 2005	$3600$ s
			Cowie	May 2005	$10800$ s
			Cowie	Mar. 2008	$12600$ s
SDSS J1602+4228	${16}^h{02}^{m}53\stackrel{s}{.}98$	$+42°28\prime 24\stackrel{\prime \prime }{.}94$	Cowie	May. 2005	$10620$ s
			Cowie	Aug. 2005	$3600$ s
SDSS J1623+3112	${16}^h{23}^{m}31\stackrel{s}{.}81$	$+31°12\prime 00\stackrel{\prime \prime }{.}53$	Becker	Jan. 2004	$3600$ s
SDSS J1630+4012	${16}^h{30}^{m}33\stackrel{s}{.}90$	$+40°12\prime 09\stackrel{\prime \prime }{.}69$	Sanchez	Jun. 2003	$17700$ s
			Sargent	Oct. 2010	$2000$ s
CFHQS J1641+3755	${16}^h{41}^{m}21\stackrel{s}{.}73$	$+37°55\prime 20\stackrel{\prime \prime }{.}15$	Willott	Oct. 2007	$2400$ s
SDSS J2054-0005	${20}^h{54}^{m}06\stackrel{s}{.}49$	$-00°05\prime 14\stackrel{\prime \prime }{.}80$	Sargent	Oct. 2010	$12000$ s
CFHQS J2229+1457	${22}^h{29}^{m}01\stackrel{s}{.}65$	$+14°57\prime 09\stackrel{\prime \prime }{.}00$	Willott	Sep. 2008	$3600$ s
			Prochaska**spectrum taken with LRIS	Sep. 2016	$3600$ s
SDSS J2315-0023	${23}^h{15}^{m}46\stackrel{s}{.}57$	$-00°23\prime 58\stackrel{\prime \prime }{.}10$	Becker	Jan. 2006	$1200$ s
			Sargent	Oct. 2010	$27000$ s
CFHQS J2329-0301	${23}^h{29}^{m}08\stackrel{s}{.}28$	$-03°01\prime 58\stackrel{\prime \prime }{.}80$	Willott	Oct. 2007	$10800$ s
			Willott	Sep. 2008	$14400$ s

Table 2. Table 2 Overview of our data sample and the measurements of the proximity zone sizes.

PSO J0226+0302	$6.5412$	$100$	[C II]	13	$-27.33$	$10$	$3.64\pm 0.13$	$3.20\pm 0.11$
CFHQS J0210-0456	$6.4323$	$100$	[C II]	12	$-24.53$	$2$	$1.32\pm 0.13$	$3.47\pm 0.34$
SDSS J1148+5251	$6.4189$	$100$	[C II]	15	$-27.62$	$34$	$4.58\pm 0.13$	$3.59\pm 0.10$
CFHQS J2329-0301	$6.417$	$270$	Mg II	9	$-25.25$	$2$	$2.45\pm 0.35$	$4.86\pm 0.70$
SDSS J0100+2802	$6.3258$	$100$	[C II]	17	$-29.14$	$41$	$7.12\pm 0.13$	$3.09\pm 0.06$
SDSS J1030+0524	$6.309$	$270$	Mg II	3	$-26.99$	$26$	$5.93\pm 0.36$	$5.95\pm 0.36$
SDSS J1623+3112	$6.2572$	$100$	[C II]	18	$-26.55$	$9$	$5.05\pm 0.14$	$6.03\pm 0.16$
CFHQS J0050+3445	$6.253$	$270$	Mg II	9	$-26.70$	$18$	$4.09\pm 0.37$	$4.60\pm 0.41$
CFHQS J0227-0605	$6.20$	$1000$	Ly$\alpha$	6	$-25.28$	$3$	$1.60\pm 1.37$	$3.15\pm 2.69$
PSO J0402+2452	$6.18$	$1000$	Ly$\alpha$	16	$-26.95$	$13$	$4.17\pm 1.38$	$4.26\pm 1.40$
CFHQS J2229+1457	$6.1517$	$100$	[C II]	15	$-24.78$	$2^c$	$0.45\pm 0.14$	$1.07\pm 0.33$
SDSS J1250+3130	$6.15$	$1000$	Ly$\alpha$-O I-Si IV	2	$-26.53$	$8$	$6.59\pm 1.38$	$7.93\pm 1.66$
ULAS J1319+0950	$6.1330$	$100$	[C II]	11	$-27.05$	$10$	$3.84\pm 0.14$	$3.77\pm 0.14$
SDSS J2315-0023	$6.117$	$1000$	Ly$\alpha$	5	$-25.66$	$15$	$3.70\pm 1.39$	$6.26\pm 2.36$
SDSS J1602+4228	$6.09$	$1000$	Ly$\alpha$-O I	1	$-26.94$	$21$	$7.11\pm 1.40$	$7.28\pm 1.43$
SDSS J0303-0019	$6.078$	$270$	Mg II	7	$-25.56$	$2$	$2.21\pm 0.38$	$3.88\pm 0.67$
SDSS J0842+1218	$6.069$	$270$	Mg II	10	$-26.91$	$10$	$6.47\pm 0.38$	$6.71\pm 0.39$
SDSS J1630+4012	$6.065$	$270$	Mg II	7	$-26.19$	$15$	$4.80\pm 0.38$	$6.59\pm 0.52$
CFHQS J1641+3755	$6.047$	$270$	Mg II	9	$-25.67$	$4$	$3.98\pm 0.38$	$6.71\pm 0.64$
SDSS J2054-0005	$6.0391$	$100$	[C II]	11	$-26.21$	$17$	$3.17\pm 0.14$	$4.32\pm 0.19$
SDSS J1137+3549	$6.03$	$1000$	Ly$\alpha$-O I-Si IV	2	$-27.36$	$24$	$6.98\pm 1.42$	$6.06\pm 1.23$
SDSS J0818+1723	$6.02$	$1000$	Ly$\alpha$	2	$-27.52$	$8$	$5.89\pm 1.42$	$4.80\pm 1.16$
SDSS J1306+0359	$6.016$	$270$	Mg II	4	$-26.81$	$41$	$5.39\pm 0.38$	$5.80\pm 0.41$
ULAS J0148+0600	$5.98$	$270$	Mg II	14	$-27.39$	$30$	$6.03\pm 0.39$	$5.18\pm 0.33$
SDSS J1411+1217	$5.904$	$270$	Mg II	4	$-26.69$	$32$	$4.60\pm 0.39$	$5.19\pm 0.44$
SDSS J1335+3533	$5.9012$	$100$	CO	8	$-26.67$	$7$	$0.78\pm 0.15$	$0.89\pm 0.17$
SDSS J0840+5624d	$5.8441$	$100$	CO	8	$-27.24$	$28$	$0.88\pm 0.15$	$0.80\pm 0.13$
SDSS J0005-0006	$5.844$	$270$	Mg II	10	$-25.73$	$15$	$2.87\pm 0.40$	$4.73\pm 0.66$
SDSS J0002+2550	$5.82$	$1000$	Ly$\alpha$-O I	1	$-27.31$	$62$	$5.43\pm 1.49$	$4.81\pm 1.31$
SDSS J0836+0054	$5.810$	$270$	Mg II	4	$-27.75$	$120$	$5.06\pm 0.40$	$3.77\pm 0.30$
SDSS J0927+2001	$5.7722$	$100$	CO	8	$-26.76$	$7$	$4.68\pm 0.15$	$5.14\pm 0.16$

Table 3. Table 3 Details of the continuum modeling of each quasar.

object	normalization to unity at rest-frame	set of PCS	number of PCS
SDSS J0002+2550	$1280$ Å	Suzuki (2006)	$5$
SDSS J0005-0006	$1280$ Å	Pâris et al. (2011)	$5$
CFHQS J0050+3445	$1280$ Å	Pâris et al. (2011)	$5$
SDSS J0100+2802	$1280$ Å	Pâris et al. (2011)	$5$
ULAS J0148+0600	$1280$ Å	Pâris et al. (2011)	$5$
CFHQS J0210-0456	$1245$ Å	Pâris et al. (2011)	$5$
PSO J0226+0302	$1280$ Å	Pâris et al. (2011)	$5$
CFHQS J0227-0605	$1265$ Å	Pâris et al. (2011)	$5$
SDSS J0303-0019	$1280$ Å	Pâris et al. (2011)	$7$
PSO J0402+2452	$1280$ Å	Suzuki (2006)	$5$
SDSS J0818+1723	$1280$ Å	Pâris et al. (2011)	$5$
SDSS J0836+0054	$1280$ Å	Suzuki (2006)	$5$
SDSS J0840+5624	$1280$ Å	Pâris et al. (2011)	$5$
SDSS J0842+1218	$1280$ Å	Pâris et al. (2011)	$5$
SDSS J0927+2001	$1280$ Å	Pâris et al. (2011)	$5$
SDSS J1030+0524	$1280$ Å	Pâris et al. (2011)	$5$
SDSS J1137+3549	$1280$ Å	Suzuki (2006)	$5$
SDSS J1148+5251	$1280$ Å	Pâris et al. (2011)	$5$
SDSS J1250+3130	$1280$ Å	Pâris et al. (2011)	$5$
SDSS J1306+0359	$1285$ Å	Pâris et al. (2011)	$5$
ULAS J1319+0950	$1280$ Å	Suzuki (2006)	$5$
SDSS J1335+3533	$1280$ Å	Pâris et al. (2011)	$5$
SDSS J1411+1217	$1280$ Å	Pâris et al. (2011)	$3$
SDSS J1602+4228	$1280$ Å	Suzuki (2006)	$5$
SDSS J1623+3112	$1280$ Å	Pâris et al. (2011)	$5$
SDSS J1630+4012	$1280$ Å	Suzuki (2006)	$5$
CFHQS J1641+3755	$1280$ Å	Pâris et al. (2011)	$3$
SDSS J2054-0005	$1280$ Å	Pâris et al. (2011)	$5$
CFHQS J2229+1457	$1250$ Å	Pâris et al. (2011)	$7$
SDSS J2315-0023	$1280$ Å	Pâris et al. (2011)	$5$
CFHQS J2329-0301	$1230$ Å	Pâris et al. (2011)	$5$

Equations

$$R_{\mathrm{ion}}\approx {\left(\frac{3{\dot{N}}_{\gamma }{t}_{\mathrm{Q}}}{4\pi n_{\mathrm{H}}{f}_{\mathrm{HI}}}\right)}^{1/3},$$

$$|q_{\lambda }⟩\approx |{\mu }_{\lambda }⟩+\sum _{i=1}^m{\alpha }_i|{\xi }_{i,\lambda }⟩,$$

$$\mathbf{\alpha }=\mathbf{P}\cdot {\mathbf{\alpha }}_{\mathrm{red}}.$$

$$R_p\approx 5.57\mathrm{pMpc}\times 10^{-0.4M_{1450}/2.35}.$$

$$R_p\approx 5.03\mathrm{pMpc}\times 10^{-0.4M_{1450}/2.45}.$$

$$R_p\approx 4.71\mathrm{pMpc}\times 10^{-0.4\times (M_{1450}+27)/3.42},$$

$$R_{p,\mathrm{corr}}\approx R_p\cdot 10^{-0.4(-27-M_{1450})/2.35}.$$

$$R_{p,\mathrm{corr}}\approx 4.87\mathrm{pMpc}\times {\left(\frac{1+z}{7}\right)}^{-1.44}.$$

$$R_{p,\mathrm{corr}}\approx 5.66\mathrm{pMpc}\times {\left(\frac{1+z}{7}\right)}^{-2.45},$$

$$R_{p,\mathrm{corr}}\approx 5.06\mathrm{pMpc}\times {\left(\frac{1+z}{7}\right)}^{-1.62}.$$

$$\frac{\mathrm{d}N}{\mathrm{d}z}\approx 2.84{\left(\frac{1+z}{4.5}\right)}^{2.04},$$

$$x_{\mathrm{HI}}(t_{\mathrm{Q}})=x_{\mathrm{HI},0}+(x_{\mathrm{HI},0}-x_{\mathrm{HI},\mathrm{eq}}){\exp }^{-t_{\mathrm{Q}}/t_{\mathrm{eq}}},$$

$$t_{\mathrm{eq}}\approx \frac{1}{{\Gamma }_{\mathrm{UVB}}+{\Gamma }_{\mathrm{QSO}}}.$$

Full text

††slugcomment: Draft Version of

1. Introduction

A major goal of observational cosmology is to understand the epoch of reionization, when the universe transitioned from the cosmic “dark ages” following recombination, into the ionized universe we observe today. Despite much progress in the last decade, there are still many open questions regarding the exact timing and the morphology of the reionization process. Studies of the evolution of the (Lyman-$\alpha$) Ly$\alpha$ absorption features in the spectra of distant quasars are one of the key observational probes of this epoch. Both the steep rise in the Ly$\alpha$ optical depth of the intergalactic medium (IGM) with redshift for $z\gtrsim 5.5$, as well as the increased scatter in the measurements, suggest a qualitative change in the state of the intergalactic medium (IGM), likely resulting from a rapid rise in the volume averaged neutral fraction which is expected to occur during the end stages of reionization (Fan et al., 2006b; Becker et al., 2015). Indeed, the absence of large Gunn-Peterson (GP) troughs (Gunn & Peterson, 1965) in spectra at $z\lesssim 5.5$ indicates the epoch of hydrogen reionization must be completed at that time (McGreer et al., 2015). However, constraining the neutral hydrogen fraction of the IGM at $z\gtrsim 6$ with quasar absorption spectroscopy has proven extremely difficult, because Lyman series transitions are overly sensitive, and saturate already for volume averaged neutral hydrogen fractions $\langle f_{\rm HI}\rangle\gtrsim 10^{-4}$.

In this paper we focus on a different technique that has been applied in previous work to constrain the neutral gas fraction at $z\sim 6$. Near the end of the cosmic reionization epoch, luminous quasars exhibit a region of enhanced transmission immediately blueward of the Ly$\alpha$ emission line in the so-called quasar proximity zone, before the onset of near complete Gunn-Peterson absorption. This enhanced transmission is caused by the radiation from the quasar itself, which ionizes the surrounding IGM (see e.g. Madau & Rees, 2000; Cen & Haiman, 2000; Haiman & Cen, 2001; Wyithe et al., 2005; Bolton & Haehnelt, 2007b, a; Lidz et al., 2007; Bolton et al., 2011; Keating et al., 2015). It has been argued that the evolution of the proximity zone sizes with redshift can constrain the late stages of the reionization epoch (Fan et al., 2006b; Carilli et al., 2010; Venemans et al., 2015).

Haiman & Cen (2001) showed that in a very simplistic picture of reionization in which isolated ionized H II regions expand into a neutral uniform ambient IGM, the (proper) size $R_{\rm ion}$ of the ionized region around the quasar depends on the neutral hydrogen fraction $f_{\rm HI}$, the rate at which ionizing photons are emitted $\dot{N}_{\gamma}$, and the quasar age $t_{\rm Q}$:

[TABLE]

where $n_{\rm H}$ is the hydrogen number density. In eqn. (1) recombinations of the ionized gas inside the H II region are neglected, which should be unimportant on the timescale that the quasar turned on, because the recombination timescale of hydrogen is comparable to the Hubble time. However, there are a number of other effects that are not taken into consideration in this relation, such as overlapping ionized H II regions, large-scale structure effects, and pre-ionization by local galaxies or the clumpiness of the IGM (see also Fan et al., 2006b; Bolton & Haehnelt, 2007b, a; Lidz et al., 2007; Maselli et al., 2007, 2009; Khrykin et al., 2016).

While eqn. (1) provides a reasonable description of the expansion rate of an H II ionization front embedded in a neutral IGM, it does not predict the distribution of residual neutral hydrogen within the ionized bubble. However, the observed size of the proximity zone $R_{p}$ will depend sensitively on the fraction and distribution of the residual neutral hydrogen, because this can cause saturated absorption well before the ionization front is reached (Bolton & Haehnelt, 2007b). Thus the extent of the ionized H II bubble around the quasar can be significantly larger than the observed proximity zone, which is defined to end where the transmitted flux drops below the $10\%$ level (see e.g. Fan et al., 2006b; Carilli et al., 2010) and hence does not necessarily probe the location of the ionization front.

Furthermore Bolton & Haehnelt (2007b) showed that a quasar embedded in a highly ionized IGM can produce a proximity zone that appears qualitatively similar to that of one in a neutral ambient IGM. They show that the observed size of the proximity zone of a quasar within an already ionized surrounding IGM is independent of the neutral gas fraction of the ambient IGM and as such, may not provide insights into the evolution of the neutral hydrogen fraction during the epoch of reionization at all.

However, previous observational studies of the sizes of quasar proximity zones found evidence for a steep decrease in proximity zone size with increasing redshift within the redshift range of $5.7\leq z\leq 6.4$ (Fan et al., 2006b; Willott et al., 2007, 2010; Carilli et al., 2010). Although there is large scatter in the observations, this has been interpreted as a strong evolution of the neutral gas fraction of the IGM assuming that the observed proximity zones trace the extent of the ionized H II region presented in eqn. (1), i.e. $R_{p}\approx R_{\rm ion}$. With the discovery of the current highest redshift quasar $\rm ULASJ1120+0641$ (Mortlock et al., 2011) and the analysis of its proximity zone (Bolton et al., 2011; Simcoe et al., 2012; Venemans et al., 2015; Bosman & Becker, 2015), this steep decrease has become somewhat shallower, but still indicates a strong evolution in proximity zone size with redshift. However, the measurements of the proximity zone sizes contain a number of uncertainties, most importantly the large uncertainties in the systemic redshifts of the quasars and the resulting ambiguity in the beginning of the proximity zone.

The extent of the H II region around the quasar additionally depends on the age of the quasar (eqn. (1)). Even if the IGM is highly ionized prior to the quasar turning on, the gas in the ionized region responds on a finite timescale (Khrykin et al., 2016) and therefore also the observed proximity zone sizes depend on the quasar ages. Thus, assumptions about the quasar ages are required to interpret quasar proximity zones, but the age remains uncertain by several orders of magnitude (Martini, 2004).

To be more precise, we will distinguish between several different timescales that govern the duration of quasar activity. The duty cycle $t_{\rm dc}$ refers to the total time that galaxies shine as active quasars integrated over the age of the universe. However, quasar activity could be episodic, and we refer to the duration $t_{\rm episodic}$ of a single emission episode as the episodic lifetime. The sizes of quasar proximity zones actually depend on the quasar age. If we denote by $t=0$ the time at which light emitted by a quasar just reaches our telescopes on Earth, then the quasar age, which we will henceforth refer to as $t_{\rm Q}$, is defined such that the quasar actually turned on at a time $-t_{\rm Q}$ in the past. Strictly speaking the age $t_{\rm Q}$ is a lower limit on the quasar episodic lifetime $t_{\rm episodic}$, because the quasar episode might indeed continue, which we can only record on Earth if we could conduct observations in the future.

The quasar duty cycle can be constrained by comparing the number density of quasars to their host dark matter halo abundance inferred from their clustering strength (Haiman & Hui, 2001; Martini & Weinberg, 2001; Martini, 2004; White et al., 2008; Furlanetto & Lidz, 2011). But to date this method has yielded only weak constraints on $t_{\rm dc}\sim 10^{6}-10^{9}$ yr owing to uncertainties in how quasars populate dark matter halos (Shen et al., 2009; White et al., 2012; Conroy & White, 2013; Cen & Safarzadeh, 2015). An upper limit on the duty cycle of quasars, $t_{\rm dc}<10^{9}$ yr, is set by the observed evolution of the quasar luminosity function, since the whole quasar population rises and falls over roughly this timescale (e.g. Osmer, 1998).

Moreover, these constraints on the duty cycle of quasars do not give insights into the duration of the individual accretion episodes of the quasar activity. A population of quasars emitting $1000$ individual bursts each with a lifetime of $t_{\rm episodic}\sim 10^{5}$ yr, would be indistinguishable from quasars with steady continuous emission for $t_{\rm episodic}\sim 10^{8}$ yr. A lower limit on the individual quasar bursts of $t_{\rm episodic}\sim 3\times 10^{4}$ yr is based on the argument that quasars need to maintain their ionizing luminosity long enough to explain the observed proximity zones in the Ly$\alpha$ forest (e.g. Bajtlik et al., 1988; Khrykin et al., 2016).

It has been argued, that in order to grow the observed sizes of SMBHs, i.e. $M_{\rm SMBH}\sim 10^{9}-10^{10}M_{\sun}$ at $z\sim 6-7$ (Mortlock et al., 2011; Venemans et al., 2013; De Rosa et al., 2014; Wu et al., 2015), the quasars require very massive initial seeds and need to accrete matter over timescales comparable to the age of the Universe at these high redshifts (Hopkins & Hernquist, 2009; Volonteri, 2010, 2012). It follows that the duty cycle of these quasars needs to be of the order of the Hubble time. If quasar activity is episodic, the episodic bursts need to be very long or the quiescent time in between the bursts needs to be short for the SMBHs to grow to their observed sizes.

The dozens of $z\sim 6$ quasars that have been uncovered over the past decade (see Bañados et al., 2016, for a recent compilation) from wide field surveys results in many new quasars for proximity zone measurements. The statistical power of these data alone motivates revisiting this type of analysis to further understand the resulting constraints on reionization and the quasar emission timescales. In this paper we re-investigate the evolution of the quasar proximity zone sizes at $z\geq 5.77$ with a homogeneous analysis of a significantly enlarged sample of $31$ quasar spectra with higher quality data, including updated and more precise redshift measurements from CO and [C II] line observations, and consistently measured absolute magnitudes. We then compare our analysis with state-of-the-art radiative transfer simulations in order to better understand and interpret our measurements.

This paper is structured as follows: we describe our data set and the reduction pipeline in § 2. In § 3 we specify the methods we use for the continuum normalization and the measurement of the proximity zone size of each quasar. In § 4 we describe a suite of radiative transfer simulations that we compare to our observations. We show our measurements of the proximity zone sizes in § 5 and discuss their evolution with redshift. We highlight three exceptionally small proximity zones in our sample in § 6, and discuss possible explanations for their origin. We summarize and conclude in § 7.

2. High Redshift Quasar Sample

Our initial data set consists of $34$ quasar spectra at $5.77\leq z_{\rm em}\leq 6.54$, observed at optical wavelengths ($4000$ Å - $10000$ Å) with the Echellette Spectrograph and Imager (ESI; Sheinis et al., 2002) at the Keck II Telescope in the years 2001 to 2016. We collected data from the Keck Observatory Archive111https://koa.ipac.caltech.edu/cgi-bin/KOA/nph-KOAlogin and supplemented it with our own observations. The slit widths used range from $0.75"-1.0"$, giving a resolution of $R\approx 4000-5400$. The exposure times vary from $0.3$ h $\lesssim t_{\rm exp}\lesssim 25$ h resulting in median signal-to-noise ratios in the quasar continuum at rest-frame wavelength of $1250$Å -$1280$Å  in the range of $2\lesssim\rm S/N\lesssim 108$ per pixel.

We obtained spectra of the four quasars $\rm PSOJ0226+0302$222also known as PSO J036.5078+03.0498 (Venemans et al., 2015)., $\rm PSOJ0402+2452$333also known as PSO J060.5529+24.8567 (Bañados et al., 2016). , $\rm SDSSJ0100+2802$, and $\rm SDSSJ1137+3549$ with ESI on on January 11th and 12th, 2016. We used a $1"$ slit and exposure times varied between $1-3.25$ hours.

For one object $\rm CFHQSJ0227-0605$ of interest because of its small proximity zone (see § 6), we obtained an additional spectrum with the Low Resolution Imaging Spectrometer (LRIS; Oke et al., 1995) on the Keck I Telescope with an exposure time of $3600$ s. The observations were conducted on September 16th, 2016 using the $600/7500$ grating and a slit width of $1"$, resulting in a resolution of $R\approx 1800$. Note that we only use this LRIS spectrum for follow-up analysis, i.e. searching for metal absorption lines. For the main analysis of the proximity zone of this quasar we use the ESI spectrum, in order to be consistent with the remaining data sample.

The details of all of the observations can be found in Table 1.

2.1. Data Reduction

We reduce all spectra uniformly using the ESIRedux pipeline444http://www2.keck.hawaii.edu/inst/esi/ESIRedux/ developed as part of the XIDL555http://www.ucolick.org/~xavier/IDL/ suite of astronomical routines in the Interactive Data Language (IDL). This pipeline employs standard data reduction techniques which can be summarized as follows: Images are overscan subtracted, flat fielded using a normalized flat field image, and then wavelength calibrated using a wavelength image constructed from afternoon arc lamp calibration images. Objects are identified in the science frames, and then background subtracted using $B$-spline fits (Kelson, 2003; Bernstein et al., 2015) to object free regions of the slit. Object profiles are also fit with $B$-splines, and optimal extraction is performed on the sky-subtracted frames. One-dimensional spectra of overlapping echelle orders are combined to produce a final spectrum for each exposure, and individual exposures are co-added into our final one-dimensional spectra. See (Bochanski et al., 2009) for a more detailed description of the algorithms used.

We further optimized the XIDL ESI pipeline to improve the data reduction for high redshift quasars. The main improvement was to remove the fringing pattern from reddest orders by differencing two images (ideally taken during the same run) with similar exposure times, analogous to the standard difference imaging techniques performed for near-infrared observations. However, this procedure only works on dithered exposures for which the trace of the science object lands at different spatial locations on the slit. Since not every observer dithered their object along the slit it was not possible for us to apply this procedure to $\approx 10\%$ of the exposures that we took from the archive.

We also co-added exposures from different observing runs taken by different PIs, resulting in higher $\rm S/N$ data for some quasars than previous analyses of these objects. We weight each one-dimensional spectrum by its squared signal-to-noise ratio ($\rm S/N^{2}$) that was determined in the quasar continuum region of each spectrum, i.e. at wavelengths longer than the Ly$\alpha$ emission line. In this way spectral regions with low or no transmitted flux, which are common in high redshift quasar spectra, are weighted by the same $\rm S/N$ ratio of regions with more transmitted flux.

2.2. Quasar Properties

Determining precise redshifts for quasars is very challenging due to the broad widths of emission lines, Gunn-Peterson absorption, and offsets between different ionization lines (Gaskell, 1982; Tytler & Fan, 1992; Vanden Berk et al., 2001; Richards et al., 2002; Shen et al., 2016). Most quasars show strong internal motions and winds, which displace many of the emission lines, such as the Ly$\alpha$ line or far-ultraviolet (far-UV) lines, far from the systemic redshift of the host galaxy. Thus the most precise determination of the location of the quasar are low ionization lines such as Mg II lines or, even better, emission lines from the molecular gas reservoir of the host galaxy itself, such as CO or [C II] lines.

We account for uncertainties in the systemic redshift of each quasar depending on the measured emission lines to determine their redshifts. For redshifts determined from the detection of a [C II] line at $158\,\rm\mu$m or a CO line we assign a redshift error of $\Delta v=100$ km/s. For quasars with a redshift measurement from a Mg II line we assume a redshift error of $\Delta v=270$ km/s, in order to account for the dispersion between the redshift of the Mg II line and the redshift of the host galaxy (e.g. Richards et al., 2002; Hennawi et al., 2006; Shen & Ménard, 2012; Venemans et al., 2016; Shen et al., 2016), and for the remaining quasars for which the redshift was determined by the Ly$\alpha$ line or far-UV lines, which usually suffer from large velocity offsets, we assume a redshift error of $\Delta v=1000$ km/s. Note that these uncertainties are chosen to be very conservative (Shen et al., 2016), because the redshift measurements are taken from the literature, measured by various authors possibly using different data and methods. For a quasar at $z=6$ these offsets in velocity result in a distance uncertainty of $\Delta R\approx 0.14$ proper Mpc (pMpc) for quasars with redshift measurements from CO or [C II] lines, $\Delta R\approx 0.39$ pMpc for quasars with redshift measurements from their Mg II line, and $\Delta R\approx 1.43$ pMpc for other rest-frame UV lines.

We take the absolute magnitudes $M_{\rm 1450}$ defined at $\lambda_{\rm rest}=1450$ Å in the rest-frame from Bañados et al. (2016), who determined $M_{\rm 1450}$ for all sources in a consistent way. They assume a power law continuum slope $\alpha_{\nu}=-0.3$ and require the extrapolated $y_{\rm P1}$- or $J$-band magnitudes (at $\lambda_{\rm eff}=9627.7$ Å or $\lambda_{\rm eff}=12444.0$ Å, respectively) to be consistent with their measurements, since most optical quasar spectra at these high redshifts have limited wavelength coverage at $\lambda_{\rm rest}=1450$Å. Table 2 summarizes the properties of all quasars in our sample.

Note that we exclude three quasars ($\rm SDSSJ1048+4637$, $\rm SDSSJ0353+0104$, and $\rm ULASJ0203+0012$) from this analysis, because of their broad absorption line (BAL) features, which can contaminate the proximity zone and make it difficult to determine a precise redshift or the continuum level of these quasars. For a similar reason we exclude the object $\rm SDSSJ0840+5624$ from our analysis of the proximity zones. In this spectrum we found an absorption system with associated Ly$\alpha$ and metal line absorption close to the quasar that attenuates enough flux such that it results in a spuriously small proximity zone. In this case, as in the case of BAL quasars, an absorption system associated with the quasar itself contaminates the proximity zone and thus we exclude this object from our analysis. Note that we do not model such intrinsic absorption line systems associated with the quasar or its host galaxy in the radiative transfer simulations (see § 4), such that any comparison would result in biases if we were not excluding these objects. We further discuss the spectrum of this excluded object $\rm SDSSJ0840+5624$ in Appendix A.

3. Methods

3.1. Quasar Continuum Normalization

We normalize all quasar spectra to unity in a region free of emission lines at $\lambda_{\rm rest}=1280$ Å. Note that for a handful of quasars, that have very low signal-to-noise ratio data at this wavelength, we normalize their spectra to unity at a slightly different wavelength. We then estimate the quasar continua with principal component spectra (PCS) from a principal component analysis (PCA) of low redshift quasar spectra (Suzuki, 2006; Pâris et al., 2011). The idea of the PCA is to represent the continuum spectrum $\ket{q_{\lambda}}$ of a quasar by a reconstructed spectrum that consists of a mean spectrum $\ket{\mu_{\lambda}}$ and a sum of $m$ weighted PCS $\ket{\xi_{\lambda}}$, where the index $\lambda$ denotes the wavelength. Hence

[TABLE]

where the index $i$ refers to the $i$th PCS and $\alpha_{i}$ is its weight. For the majority of our continuum fits, we use the PCS from Pâris et al. (2011) that were derived to characterize quasar continua from $78$ high-quality spectra of bright $z_{\rm em}\sim 3$ quasars in the Sloan Digital Sky Survey (SDSS).

The quasars in our sample are all at very high redshift and thus experience absorption due to the intervening residual neutral hydrogen along the line of sight bluewards of the Ly$\alpha$ emission line, i.e. in the Ly$\alpha$ forest. Thus we estimate the quasar continuum solely from wavelengths redwards of rest-frame Ly$\alpha$ which suffer only from modest absorption by metal lines. Pâris et al. (2011) provide a set of PCS for wavelengths $1215.67$ Å$\leq\lambda_{\rm rest}\leq 1600$ Å, whose coefficients $\bm{\alpha_{\rm red}}$ we estimate by $\chi^{2}$ minimization using the noise vector from the spectra.

We use a projection matrix $\bm{P}$ to transfer the estimated coefficients for the PCS redwards of Ly$\alpha$, $\bm{\alpha_{\rm red}}$, onto coefficients $\bm{\alpha}$ for a set of PCS that cover the entire spectral region between $1020$ Å $\leq\lambda_{\rm rest}\leq 1600$ Å, i.e.

[TABLE]

This projection matrix $\bm{P}$ has been calibrated by Pâris et al. (2011) using the set of PCS for both the red wavelength side only and the whole spectral region covering wavelengths bluewards and redwards of Ly$\alpha$.

Eqn. (2) gives us a model for the continuum of each quasar that includes the (absorbed) blue side of the spectrum. An example of a quasar spectrum from our data set and its continuum model is shown in Fig. 1. Note that in most cases we take five PCS on the red side into account to estimate the continuum model. However, a visual inspection of the continuum models reveals that for four quasars a continuum model with three or seven PCS results in a better fit.

Our PCA modeling of quasar continua is limited by the finite number of principal components and the selection of quasars that were used to create the PCS. For a few of the quasars in our sample, the PCS from Pâris et al. (2011) do not provide an acceptable fit to the continuum, particularly for the fainter quasars in our sample. For these objects we instead use the set of PCS provided by Suzuki (2006), based on fainter lower redshift quasars, which qualitatively provide better fits to this handful of objects. When using the Suzuki (2006) PCS basis, we determine the continuum in the Ly$\alpha$ forest region using the method they advocate. Namely, we determine the continua at blue wavelengths $\lambda_{\rm rest}<1216.67\,{\rm\AA}$ using the best-fit continua to the red-side $\lambda_{\rm rest}>1216.67\,{\rm\AA}$.

The details of the continuum normalization are listed in Table 3 in Appendix B.

3.2. Measuring the Sizes of Quasar Proximity Zones

In order to calculate the proximity zone sizes for each quasar we adopt the standard definition used in the literature (see Fan et al., 2006b; Willott et al., 2007, 2010; Carilli et al., 2010; Venemans et al., 2015). Namely, we take the continuum normalized quasar spectra and smooth them with a $20$Å-wide (observed frame) boxcar function. This smoothing scale corresponds to $\approx 1.0$ pMpc or a $\approx 705$ km/s window at $z=6$. We define the proximity zone size as the distance to the first of three consecutive pixels666We checked whether it changes the measurements of the proximity zone sizes when requiring up to ten pixels to be below the $10\%$ level, but we do not observe a significant difference. on the blue side of the Ly$\alpha$ emission line that show a drop of the smoothed flux below the $10\%$ level. This is demonstrated in Fig. 2 and Fig. 3 for subsets of bright ($-27.5\leq M_{\rm 1450}\leq-26.5$) and faint ($-25.8\leq M_{\rm 1450}\leq-24.5$) quasars from our data sample, respectively. The depicted quasars all cover a similar range in luminosity, but nevertheless show a wide range of proximity zone sizes, between $0.8$ pMpc $\lesssim R_{p}\lesssim 7.1$ pMpc for bright quasars and $0.5$ pMpc $\lesssim R_{p}\lesssim 4.0$ pMpc for fainter quasars. Similar plots for the remaining objects in our quasar sample are shown in Appendix C. All measurements of the proximity zone sizes $R_{p}$ for our data set are listed in Table 2.

4. Radiative Transfer Simulations

To interpret our measurements of $R_{p}$, we run a series of radiative transfer simulations of the effect of quasar ionizing photons on the IGM along the line of sight similar to those performed by Bolton & Haehnelt (2007b). We apply the one-dimensional ionizing radiative transfer code from Davies et al. (2016) to skewers from a $100$ Mpc$/h$ Eulerian hydrodynamical simulation run with the Nyx code (Almgren et al., 2013; Lukić et al., 2015) at $z=6$ with $4096^{3}$ baryonic (Eulerian) grid elements and dark matter particles. The radiative transfer computes the time-dependent ionization and recombination of six species ($e^{-}$, H I, H II, He I, He II, He III) as well as the associated photoionization heating and cooling by various processes including adiabatic cooling due to the expansion of the Universe and inverse Compton cooling off cosmic microwave background (CMB) photons (see Davies et al., 2016, for details). We use $900$ skewers of density, temperature, and peculiar velocity drawn along the $x$, $y$, and $z$ grid axes from the centers of the $150$ most massive dark matter halos in the Nyx simulation, corresponding to halo masses $M_{\rm h}\gtrsim 4\times 10^{11}$ M⊙777In detail we find that starting the skewers from more or less massive halos has a negligible impact on the resulting proximity zone sizes, in agreement with Bolton & Haehnelt (2007b) and Keating et al. (2015).. Given the relatively coarse resolution of the Nyx simulation ($25$ kpc$/h$) and the lack of metal-line cooling, star formation, and feedback processes, we do not expect to resolve or accurately model dense gas inside of galaxies. For this reason we have removed one object in our sample ($\rm SDSSJ0840+5624$; see Appendix A) which exhibits strong associated metal line absorption at the quasar redshift that likely arises from dense gas in the vicinity of the quasar or from its host halo.

Quasar absolute magnitudes were converted to specific luminosity at the hydrogen ionizing edge ($\nu=\nu_{\rm HI}$) using the Lusso et al. (2015) spectral template, and for consistency with previous studies we assume that the spectrum at $\nu>\nu_{\rm HI}$ behaves as a power law888Note that the magnitudes for the quasars in our data sample have been calibrated by Bañados et al. (2016) according to the template from Selsing et al. (2016) instead of Lusso et al. (2015). However, due to the similarities between the two templates at observed wavelengths (redward of Ly$\alpha$), we expect only a very minor inconsistency in the inferred ionizing luminosity. with $L_{\nu}\propto\nu^{-1.7}$. We assume that the quasar turns on abruptly and emits at constant luminosity for its entire age $t_{\rm Q}$, i.e. a so-called “light bulb” model. For the quasar age we assume a fiducial value of $t_{\rm Q}=10^{7.5}$ yr, but we later investigate the dependence of $R_{p}$ on $t_{\rm Q}$ in § 6.3.3. We consider two initial conditions for the ionization state of gas in the simulation: either the gas is initially highly ionized by a uniform ionizing background or the gas is initially completely neutral. In the highly ionized case, we add ionizing radiation due to the ultraviolet background (UVB) to each cell leading to an ionization rate $\Gamma_{\rm UVB}=2\times 10^{-13}$ s*-1*, consistent with observations of Ly$\alpha$ forest opacity at $z\sim 6$ (e.g. Wyithe & Bolton, 2011) and resulting in a neutral fraction of $f_{\rm HI}=1.5\times 10^{-4}$, with a spectrum characteristic of galactic sources (i.e. a sharp cutoff above the He II ionizing edge). For simplicity we assume that $\Gamma_{\rm UVB}$ is constant with redshift, but we note that the resulting proximity zone sizes are insensitive to changes in $\Gamma_{\rm UVB}$ of a factor of a few, a non-trivial point that we will discuss in more detail in future work (Davies et al., in prep.). We compute radiative transfer outputs for quasars of varying luminosity and redshift bracketing the properties of our observed quasar sample in in both ionized and neutral scenarios. We assume here that the overdensity field does not significantly evolve with redshift from our $z=6$ output across the redshift range we study, and simply re-scale physical densities by $(1+z)^{3}$.

Ly$\alpha$ forest spectra are computed by combining the neutral fraction $f_{\rm HI}$ and gas temperatures from the radiative transfer simulation with the velocity field from the Nyx simulation, summing the absorption from each gas element using the efficient Voigt profile approximation of Tepper-García (2006). Finally, to retrieve $R_{p}$ from the Ly$\alpha$ forest spectra, we perform a similar analysis as is applied to the real spectroscopic data: we smooth the spectra with a boxcar filter of $20$ Å in the observed frame, and then locate where the transmitted flux first drops below $10\%$. A few skewers through the simulation box giving mock Ly$\alpha$ forest spectra in a highly ionized IGM are shown in Fig. 4.

In this work we present the general trends of $R_{p}$ with quasar parameters, leaving a more detailed discussion of the structure and evolution of these simulated proximity effect spectra to future work (Davies et al., in prep.).

5. The Redshift Evolution of Proximity Zone Sizes

In this section we will present our measurements of the proximity zone sizes for the ensemble of quasar spectra in our data set. Since these quasars cover a wide range of luminosities, we first put them on a common luminosity scale, in order to facilitate a study of the redshift evolution of proximity zone sizes. In the last part of this section we compare our measurements to previous work.

5.1. Correcting for the Quasar Luminosity

We expect the size of the proximity zone of a quasar to depend on their luminosity. Proximity zone sizes may also evolve with redshift in a way that tracks the evolution of the ionization state of the IGM driven by the underlying evolution in the UVB. To study the evolution of the proximity zone size with redshift, we have to remove the dependency on quasar luminosity, and normalize our measurements to the same absolute magnitude. The approach that has been adopted in the literature, is to re-scale all quasar luminosities to a common value assuming a particular model for the luminosity dependence, in order to remove the large scatter in proximity zone sizes driven by the wide range of quasar luminosities present in the data set. This then enables one to study the redshift evolution of proximity zone sizes. However, the scaling of the proximity zone size with luminosity depends on the physical conditions of the ambient IGM. In this section we study radiative transfer simulations to better understand the expected luminosity scaling, and its dependence on the IGM ionization state. We compute the luminosity scaling in our data set and compare to these simulations. We use these simulations to advocate for an approach to luminosity correct our measurements and search for a redshift evolution.

In general, the luminosity dependence of proximity zone sizes could depend on the ionization state of the surrounding IGM. Specifically eqn. 1 indicates that the size of an ionization front $R_{\rm ion}$ expanding into a neutral IGM scales as $R\propto\dot{N}_{\gamma}^{1/3}$. On the other hand, using an analytical model Bolton & Haehnelt (2007b) showed that the proximity zone size scales as $R_{p}\propto\dot{N}_{\gamma}^{1/2}$, if the IGM is in fact highly ionized by either the quasar itself or the UVB. To better understand this scaling, we turn to our radiative transfer simulations (see § 4), which were run for an ensemble of quasars that span the same range of luminosities as our data set. Note that the redshifts of the simulated quasars and the IGM are set to $z=6$.

The simulated relation between proximity zone size and quasar magnitude for a highly ionized IGM is shown as the gray dashed line in Fig. 5, with the shaded region illustrating the $1\sigma$ scatter of simulated sizes due to cosmic variance alone (i.e. no measurement error). The relation from the simulations in a highly ionized IGM is reasonably well fit by the power-law

[TABLE]

Our simulations suggest that the proximity zone sizes in a highly ionized IGM scale as $R_{p}\propto\dot{N}_{\gamma}^{1/2.35}$, which lies in between the two theoretically expected relations for a neutral or mostly ionized ambient IGM, shown as the blue dashed-dotted and yellow dotted curves, respectively. Note that the normalization for these two analytic curves has been arbitrarily chosen such that all curves intersect at the same point. The simulated relation does not align exactly with the analytically expected relation in a highly ionized IGM, possibly due to heating effects of He II reionization.

Note however, that we obtain a very similar scaling from our radiative transfer simulations, when assuming a neutral ambient IGM:

[TABLE]

This is due to the fact that the measured proximity zone sizes $R_{p}$ do not trace the extent of the ionization front $R_{\rm ion}$ described in eqn. (1) for an assumed quasar age of $t_{\rm Q}\sim 10^{7.5}$ yr, but end much earlier within the ionized H II region around the quasar. Because the IGM is indeed highly ionized in this region, one expects a scaling similar to that of a highly ionized IGM. Hence, it is no surprise to find the two relations from the radiative transfer simulations in a highly ionized ambient IGM as well as in an originally neutral ambient IGM to behave similarly. We will come back to this in more detail in § 5.3.

Given this new found intuition from radiative transfer simulations, it is now interesting to study the luminosity scaling of proximity zone sizes in our data. To this end, we show the dependence of our measured proximity zone sizes on the quasar’s magnitude, which is proportional to the quasar luminosity, in Fig. 5. Measurements have been color coded by the emission redshift of the quasar. The errorbars of the measurements reflect the uncertainty on the proximity zone sizes due to the uncertainty of the quasar redshift (see § 2). We fit a power-law to our measurements and obtain

[TABLE]

as the best fit, shown as the red dashed line. The $1\sigma$ uncertainty of this fit is given by the shaded region, which was determined by bootstrap re-sampling of our measurements with replacement, and repeating the fit $1000$ times. Note that we do not weight our measurements by the redshift errors in this fit. In general, we find reasonable agreement between the scaling relations from the radiative transfer simulations and our data. The data seem to favor a shallower evolution, but this is complicated by the fact that we fit a sample over a large range of redshifts, and also, possibly, that exceptionally small proximity zones (see § 6) are pulling down the fit.

So which relation should we adopt to put the measured proximity zone sizes onto a common luminosity scale? Due to the large body of evidence for the IGM being mostly ionized at $z\sim 6$ (e.g. Wyithe & Bolton, 2011; Calverley et al., 2011; Becker et al., 2015; McGreer et al., 2015) we choose the relation from our radiative transfer simulations assuming a highly ionized IGM (eqn. (4)) to eliminate the dependence on quasar luminosity and normalize our measurements to a common magnitude of $M_{\rm 1450}=-27$. However, our results are only marginally dependent on this choice, since the luminosity scaling for an ionized and neutral IGM are in fact very similar, and the main conclusions of this paper are not influenced by this choice. Thus we re-scale our measured proximity zone sizes with

[TABLE]

Our measurements of the quasar proximity zone sizes are provided in Table 2, with the two rightmost columns showing the measured proximity zone $R_{p}$ in proper Mpc and the corrected proximity zone size $R_{p,\rm corr}$, which are re-scaled to an absolute magnitude of $M_{\rm 1450}=-27$ following eqn. (7).

5.2. The Redshift Evolution of Quasar Proximity Zone Sizes

After correcting our measurements for the quasar luminosity, we can study the evolution of the proximity zones with redshift. In Fig. 6 we show this redshift evolution of our luminosity corrected measurements, color coded by $M_{\rm 1450}$, in both the left and middle panel. In the left panel, the red dashed line shows the best-fit power-law fit to the redshift evolution

[TABLE]

The shaded region indicates the $1\sigma$ uncertainty on the fit determined by bootstrap re-sampling of our measurements with replacement and repeating the fit $1000$ times.

In the next section we use our radiative transfer simulations to better understand the redshift evolution of proximity zones, and show that the shallow redshift evolution that we measure is indeed expected. However, the black dashed and dash-dotted lines in the right panel of Fig. 6 are linear fits to the measurements of quasar proximity zones from previous analyses of similar data sets by Carilli et al. (2010) and Venemans et al. (2015), which show a much steeper trend. We compare our measurement to previous work in § 5.4, and discuss the sources of this discrepancy.

5.3. Understanding the Shallow Redshift Evolution of Quasar Proximity Zone Sizes

We use our radiative transfer simulations to investigate the expected redshift evolution of proximity zone sizes for both a highly ionized ambient IGM and a neutral IGM. Our results are shown in the middle panel of Fig. 6 as the gray dashed line (ionized) and blue dash-dotted line (neutral), respectively. In the case of an ionized ambient IGM the corrected quasar proximity zones follow a power-law

[TABLE]

whereas the evolution in the case of a neutral IGM is well fit by

[TABLE]

The shaded regions around the curves show the scatter about these relations due to cosmic variance for quasars with magnitude $M_{\rm 1450}=-27$, but they slightly over- and underpredict the scatter for higher and lower luminosity quasars, respectively. The two scenarios of an ionized and a neutral IGM both result in a relatively shallow redshift evolution of proximity zones and are thus consistent with our measured redshift evolution shown in the left panel as the red dashed curve, but inconsistent with the steep evolution found in previous work (right panel).

Why should the redshift evolution of quasar proximity zones be so shallow? The end of the proximity zone around quasars is defined as the location where the smoothed transmitted flux drops below the $10\%$ transmission level and a corresponding limiting optical depth of $\tau_{\rm lim}=2.3$ is reached. Bolton & Haehnelt (2007b) determine in eqn. (8) of their paper the ionization rate $\Gamma_{\rm lim}$ that is necessary to produce a neutral gas fraction $f_{\rm lim}$ that results in the required limiting optical depth $\tau_{\rm lim}$. Assuming a highly ionized gas in ionization equilibrium and making reasonable assumptions about the gas temperature of the IGM at $z\sim 6$, one obtains an ionization rate of $\Gamma_{\rm lim}\sim 4\times 10^{12}$ s*-1* at the end of the proximity zone $R_{p}$, which is an order of magnitude larger than ionization rate of the UVB $\Gamma_{\rm UVB}\sim 2\times 10^{-13}$ s*-1* at this redshift (Wyithe & Bolton, 2011). Thus the total ionization rate at the end of the proximity zone $\Gamma_{\rm lim}$, which is the sum of the ionization rate of the background radiation $\Gamma_{\rm UVB}$ and the ionization rate of the quasar itself $\Gamma_{\rm QSO}$, has to be totally dominated by the latter.

As such, in our simulations $R_{p}$ is essentially independent of $\Gamma_{\rm UVB}$ and hence fairly insensitive to the neutral fraction $f_{\rm HI}$ of a highly ionized IGM, provided the quasar has been emitting light for longer than the equilibration timescale $t_{\rm eq}$ of the gas, which denotes the timescale on which the gas reaches ionization equilibrium (see § 6.3.3). The quasar age thus has to be $t_{\rm Q}\gtrsim t_{\rm eq}\sim 1/\Gamma_{\rm QSO}(R_{p})\sim 10^{5}$ yr, in order for the ambient gas to have reached ionization equilibrium.

If the IGM surrounding the quasars is instead very neutral, i.e. $f_{\rm HI}\sim 0.1-1.0$, eqn. (1) for the location of the ionization front indicates that the size of the ionized H II region around the quasar scales as $R_{\rm ion}\propto(t_{\rm Q}/f_{\rm HI})^{1/3}$. This suggests that the proximity zones could be sensitive to the neutral gas fraction (subject to a degeneracy with the quasar age). Indeed for short ages the measured proximity zone $R_{p}$ will trace the expanding ionization front $R_{\rm ion}$ around the quasar, i.e. $R_{p}\approx R_{\rm ion}$, which increases with the quasar age. However, $R_{p}$ will cease to grow further once it reaches a distance given by the $10\%$ transmission level of the flux, according to the definition of $R_{p}$. Thus even as $R_{\rm ion}$ continues to grow with age, the proximity zone size saturates, and will be insensitive to both the quasar age and neutral fraction. This implies that the measured proximity zone size provides a lower limit on the location of the ionization front $R_{\rm ion}$. The maximum size of the proximity zone $R_{p,\rm max}$ is given by the distance at which a $10\%$ flux transmission level corresponding to a limiting optical depth of $\tau_{\rm lim}=2.3$ is reached. The age of the quasar at this distance is $t_{\rm Q}>t_{\rm Q}(R_{p,\rm max})\sim 10^{6}$ yr under the assumption of a homogeneous IGM (see eqn. (13) in Bolton & Haehnelt, 2007b).

Thus for a highly plausible quasar age of $t_{\rm Q}\sim 10^{7.5}$ yr (see e.g. Martini, 2004) the proximity zones $R_{p}$ in a neutral IGM will look the same as for the highly ionized case, because the proximity zones by definition cease to grow once a flux transmission level of $10\%$ is reached, as has been previously pointed out by (Bolton & Haehnelt, 2007b). This demonstrates the insensitivity of $R_{p}$ to the neutral fraction $f_{\rm HI}$ of the IGM for an appropriate choice of quasar age (see also Khrykin et al., 2016, for an analogous argument in the context of He II proximity zones).

Regarding the redshift evolution of the proximity zone sizes, eqn. (11) in Bolton & Haehnelt (2007b) predicts a very shallow scaling of $R_{p}\propto(1+z)^{-2.25}$ due to the density evolution in the universe. The simulations in Fig. 6 reveal a similar scaling, i.e. $R_{p}\propto(1+z)^{-2.45}$ for a neutral ambient IGM and $R_{p}\propto(1+z)^{-1.62}$ for a highly ionized IGM with a quasar age of $t_{\rm Q}\sim 10^{7.5}$ yr. The true scaling does not match exactly the analytically expected scaling due to different heating effects.

Hence we conclude that the redshift evolution of $R_{p}$ is not a very useful probe of the ionization state of the IGM. In a highly ionized IGM due to the definition of proximity zones, the ionization rate at the end of the proximity zone is totally dominated by the ionization rate of the quasar itself and thus essentially independent of the background radiation and the ionization state of the surrounding IGM. In a neutral IGM, $R_{p}$ ceases to grow once a maximum size $R_{p,\rm max}$ corresponding to a limiting optical depth $\tau_{\rm lim}$ at the $10\%$ flux transmission level is reached. The observed shallow redshift evolution is thus perfectly consistent with our models for both a highly ionized and a neutral IGM for a quasar age of $t_{\rm Q}>10^{6}$ yr, in particular for a fiducial quasar age of $t_{\rm Q}\sim 10^{7.5}$ yr. None of our simulations reveal the steep evolution that was measured by previous analyses (Carilli et al., 2010; Venemans et al., 2015). Our models do not predict any large proximity zones of $R_{p}\gtrsim 10$ pMpc that would be consistent with their best-fit curve at $z\approx 5.8$. We will now investigate possible reasons leading to the discrepancy between ours and previous analyses.

5.4. Comparison to Previous Analyses

There are several possible reasons causing the differences between ours and previous analyses that result in the different measured redshift evolution of quasar proximity zones: we use a different scaling with luminosity applied to correct all proximity zone measurements and put them on the same luminosity scale, different quasar spectra999Although many of the objects are overlapping between ours and previous data sets, we co-added different exposures from different observing runs and have thus gained higher quality spectra. , different approaches to fitting continua and measuring the proximity zone sizes, and our measurements also rely on updated luminosities and redshift measurements. Below we discuss each of these differences and their impact on our results in turn.

To address the issue of different luminosity scalings we facilitate a comparison to previous work by Carilli et al. (2010) and Venemans et al. (2015) by taking their measurements of $R_{p}$ and correcting them with our luminosity correction ($R_{p,\rm corr}\propto\dot{N}_{\gamma}^{1/2.35}$ instead of their $R_{p,\rm corr}\propto\dot{N}_{\gamma}^{1/3}$; see eqn. (7)) using updated measurements of $M_{\rm 1450}$ from Bañados et al. (2016) for their objects. In their work they chose a linear relation between proximity zone size and redshift, i.e. $R_{p}\propto z$, and thus we also show a linear fit to their corrected measurements in the right panel of Fig. 6 in order to simplify the comparison. Our best fit slope to the measurements of Carilli et al. (2010) is $m\approx-7.07$ (black dotted line) and $m\approx=-6.28$ (black dash-dotted line) for measurements by Venemans et al. (2015), which is considerably steeper than the slope $m\approx-1.12$ that we obtain when applying a linear fit to our data. Thus we conclude that the difference in luminosity scaling does not cause the discrepancy between the different analyses.

However, both previous analyses contain objects in their data sample that are not included in our sample. Although the data sets from (Carilli et al., 2010) and Venemans et al. (2015) partially overlap with our sample, they also contain distinct objects. In order to verify that differences in the slope of the redshift evolution are not being driven by differences in the quasar samples, we restrict the analysis now to measurements of the $18$ quasars that are common to our sample and those of Carilli et al. (2010) and Venemans et al. (2015). For five of these quasars, our analysis uses updated redshift measurements and thus we correct the measurements of $R_{p}$ from previous work for the difference between the old and the new redshift measurements in order to have consistent redshifts and afterwards use eqn. (7) to obtain $R_{p,\rm corr}$ from their measurements of $R_{p}$. Their corrected measurements of the objects overlapping in both data sets are shown in the right panel of Fig. 6 and the best linear fit for the redshift evolution is shown as the green dashed curve. It still shows a relatively steep slope of $m\approx-5.33$ although slightly shallower than before. Thus we conclude that differences in the data sets are not driving the discrepancy in the results.

Note however, that the differences in the resulting redshift evolution of $R_{p}$ appear to be driven by a handful of objects, where we measure significantly different proximity zone sizes than previous work. In Fig. 7 we compare our measurements of the (uncorrected) $R_{p}$ to those measured by Willott et al. (2007), Willott et al. (2010), Carilli et al. (2010) (who updated the measurements originally performed by Fan et al. (2006b)), and Venemans et al. (2015) for all objects that are overlapping between our data sample and these previous analyses. We have adopted consistent redshifts for all objects here in order to facilitate a one-to-one comparison that is not driven by changes in redshift, i.e. for the objects for which we have updated redshift measurements, we correct the previously analyzed proximity zones for the difference.

The measurements which lie along the gray dashed line indicate agreement between the different analyses. In general, we measure smaller sizes of the proximity zones than previously obtained, but for most objects the measurements agree within $\Delta R_{p}\lesssim 2$ pMpc. Despite the fact the fact that for a significant fraction of these overlapping quasars the analyzed data comes from the same instrument (Keck/ESI), small differences in the measurements can be attributed to different data reduction pipelines and in many cases we have co-added data from different runs to get higher $\rm S/N$ spectra and thus the final quasar spectra might differ in quality. However, there are a few outliers in Fig. 7, for which our measurements differ significantly from previous work, i.e. $\Delta R_{p}>2$ pMpc, which we discuss in further detail.

We measure much smaller proximity zones for three quasars in particular: $\rm SDSSJ0836+0054$ ($R_{p,\rm Carilli}\approx 13.0$ pMpc vs. our measurement $R_{p}=5.06\pm 0.40$ pMpc) at $z=5.810$, $\rm SDSSJ0002+2550$, ($R_{p,\rm Carilli}\approx 11.5$ pMpc vs. our measurement $R_{p}=5.43\pm 1.49$ pMpc) at $z=5.82$, and $\rm CFHQSJ2329-0301$, ($R_{p,\rm Willott}\approx 7.0$ pMpc vs. our measurement $R_{p}=2.45\pm 0.35$ pMpc) at $z=6.417$.

The spectra of the first two objects, $\rm SDSSJ0836+0054$ and $\rm SDSSJ0002+2550$ and our measured proximity zones (black dashed lines) as well as the measurements from Carilli et al. (2010) (blue dashed lines) are shown in Fig. 8. Both analyses use the same redshift measurements and thus no differences due to redshift errors should be causing the discrepancy we see here. The spectrum of $\rm SDSSJ0836+0054$ in the upper panel shows transmitted flux throughout the spectrum with the smoothed flux oscillating around the $10\%$ level. Our spectrum nevertheless shows a significant drop below the $10\%$ level at $R_{p}\approx 5.06$ pMpc. This measured $R_{p}$ does not change significantly when fitting the continuum model with a different set of PCA components.

For $\rm SDSSJ0002+2550$ (lower panel) the size of the proximity zone is more sensitive to the precise placement of the quasar continuum level. A slightly different continuum normalization (i.e. taking a different set of PCA components) would increase the measurement of the proximity zone from $R_{p}\approx 5.43$ pMpc to $R_{p}\approx 8.85$ pMpc, but still falls considerably short of the $R_{p}=11.5$ pMpc previous measurement by Carilli et al. (2010). Another reason for the discrepancy in the $R_{p}$ measurements could be the higher signal-to-noise data in our sample. However, it remains unclear whether the differences in the continuum normalization and the better quality data could cause the whole discrepancy of the two measurements.

The third object that has a significantly different proximity zone measurement is $\rm CFHQSJ2329-0301$, whose proximity zone was determined by Willott et al. (2007) to be $R_{p}\approx 6.3$ pMpc (after updating their measurements with a new redshift measurement it is now $R_{p}\approx 7.0$ pMpc). However, most of this discrepancy between their measurement and ours can be contributed to the fact that for this particular object, Willott et al. (2007) decided after inspecting the spectrum by eye that they take the second drop of the smoothed flux below $10\%$ level as the end of the proximity zone instead of the first drop according to the standard definition of proximity zones (Fan et al., 2006b). They determine the first drop in flux below the $10\%$ transmission level to be at $R_{p}\approx 3.7$ pMpc (with an updated redshift measurements this would be $R_{p}\approx 4.4$ pMpc), which would be much more consistent with our measurement of $R_{p}=2.45\pm 0.35$ pMpc.

Another reason for the discrepancy in the redshift evolution of the proximity zone measurements is caused by the exclusion of the quasar $\rm SDSSJ1335+3533$ at $z=5.9012$ from previous analyses. For this particular object, which has a very small proximity zone of $R_{p}=0.78\pm 0.15$ pMpc, our measured size agrees with the one measured by Carilli et al. (2010). However, this object is a weak emission line quasar and has been somewhat arbitrarily excluded from previous analyses due to this fact (see Carilli et al., 2010; Wyithe & Bolton, 2011). The reason given in the aforementioned papers for is the “fundamentally different nature” of such objects. However Diamond-Stanic et al. (2009) showed that weak emission line quasars do not show significantly different UV continuum slopes apart from the emission lines and thus do not differ in their physical properties. The original reason for Fan et al. (2006b) to exclude this object from their analysis was that the only redshift measurement they had at that time came from the onset of strong Ly$\alpha$ absorption. Given that we now have an accurate redshift measurement, we include $\rm SDSSJ1335+3533$ into our analysis, which causes the slope of the redshift evolution of the proximity zones to become shallower due to its small zone. We will further discuss the implications of small proximity zones in § 6.

Further differences between the different analyses can be attributed to the different continuum fitting methods. Previous analyses applied a power-law fit to the quasar continuum redwards of the Ly$\alpha$ emission line and fitted the Ly$\alpha$ and N V lines with Gaussian curves, whereas we chose to model the continuum of each quasar with two different sets of PCA (i.e. Pâris et al. (2011) or Suzuki (2006)). Although differences in the continuum estimation can change the size of individual proximity zones, we do not expect them to alter the distribution of proximity zone sizes, since the scatter on $R_{p}$ resulting from continuum uncertainties is much smaller than the intrinsic scatter due to density fluctuations (Kramer & Haiman, 2009).

The continuum uncertainties in our analysis arising due to differences between the two sets of PCA are not very significant. The median difference between the proximity zone sizes when estimating the quasar continua with different PCA models is $\langle\Delta R_{p}\rangle=\langle|R_{p,\rm Paris}-R_{p,\rm Suzuki}|\rangle\approx 0.09$ pMpc. The scatter in the distribution of $\Delta R_{p}$ determined from the $16$th and $84$th percentile is $\sigma_{\Delta R_{p}}\approx 0.73$ pMpc.

To conclude, it is not completely clear to us, how other authors obtained some very large proximity zones, particularly for a few objects at $z<6$ that are the main drivers for the steep redshift evolution. Some of the differences between ours and previous analyses and the resulting shallower redshift evolution of the proximity zone sizes can be attributed to a variety of reasons: first, we stick to a rigorous definition of the proximity zone sizes and continue to take the first drop of the smoothed flux below the $10\%$ level, since we adopt the same rigorous treatment for the mock spectra from the radiative transfer simulations. Second, we analyze higher quality data, which can help to determine the drop in flux below the $10\%$ level easier. Third, we include all quasars of the ensemble in our analysis and do not exclude the weak emission line quasar in our sample. Finally, we analyze a larger sample of quasar spectra than previous analyses. Differences in the continuum estimation, updated redshift and magnitude measurements, and the different correction for the quasar luminosity do not have a significant influence on the discrepancy.

6. Exceptionally Small Quasar Proximity Zones

Several of the quasars we studied have particularly small proximity zones, as can be seen from Fig. 6. The proximity zones of two objects, i.e. $\rm CFHQSJ2229+1457$ and $\rm SDSSJ1335+3533$, are $R_{p}<1$ pMpc and $R_{p,\rm corr}\lesssim 1$ pMpc. Additionally, given its extreme brightness, the proximity zone of $\rm SDSSJ0100+2802$, $R_{p}=7.12\pm 0.13$ pMpc and $R_{p,\rm corr}=3.09\pm 0.06$ pMpc, is also exceptionally small. In this section, we discuss the properties of these objects and possible explanations for their small proximity zones.

6.1. Individual Objects with Small Proximity Zones

CFHQS J2229+1457

This relatively faint quasar, $M_{\rm 1450}=-24.78$ (Bañados et al., 2016), was discovered in the Canada-France High-z Quasar Survey (CFHQS) and published by Willott et al. (2009). The measurements of the proximity zone are $R_{p}=0.45\pm 0.14$ pMpc and, when normalized to a magnitude of $M_{\rm 1450}=-27$, $R_{p,\rm corr}=1.07\pm 0.33$ pMpc. The quasar has a precise redshift measurement of $z=6.1517$ from [C II] emission from the host galaxy (Willott et al., 2015) and hence a very small redshift uncertainty $\Delta v=100$ km/s resulting in an uncertainty of the proximity measurement of $\Delta R_{p}\approx 0.14$ pMpc. The measurements of the proximity zone are independent of the choice of PCS used to model the quasar continuum, i.e. continuum uncertainties do not influence our measurements.The top panel of Fig. 3 shows the continuum normalized spectrum and the proximity zone of this object. In the upper panel of Fig. 9 we show the whole spectrum, revealing no BAL features.

SDSS J1335+3533

This quasar, that has been discovered in the Sloan Digital Sky Survey (SDSS) by Fan et al. (2004), has a precise redshift measurement of $z=5.9012$ from CO $(6-5)$ emission from the host galaxy (Wang et al., 2010), and is fairly bright with a magnitude of $M_{\rm 1450}=-26.67$ (Bañados et al., 2016). The measured proximity zone size for this object is $R_{p}=0.78\pm 0.15$ pMpc and the luminosity corrected size of the proximity zone is $R_{p,\rm corr}=0.89\pm 0.17$ pMpc. The redshift uncertainty of $\Delta v=100$ km/s results in an uncertainty of the proximity zone measurement of $\Delta R_{p}\approx 0.15$ pMpc. These measurements of the proximity zone are again independent on continuum modeling uncertainties. The top panel of Fig. 2 shows the continuum normalized spectrum and the proximity zone of this object. The whole spectrum of the quasar shown in the middle panel of Fig. 9 is completely devoid of broad emission lines, as previously noticed from its discovery spectrum by Fan et al. (2004). The spectrum does not show any BAL features.

SDSS J0100+2802

This object is the brightest high redshift quasar known so far (Wu et al., 2015) with an absolute magnitude of $M_{\rm 1450}=-29.14$ (Bañados et al., 2016). Its redshift, $z=6.3258$, has been measured precisely by the detection of the [C II] line (Wang et al., 2016). Given its extreme brightness the quasar exhibits a very small proximity zone of $R_{p}=7.12\pm 0.13$ pMpc ($R_{p,\rm corr}=3.09\pm 0.06$ pMpc, when normalized to a luminosity of $M_{\rm 1450}=-27$) compared to a proximity zone size of $R_{p}=12.0\pm 2.0$ pMpc, that one would expect to see for an object this bright. The redshift uncertainty $\Delta v=100$ km/s results in an uncertainty of the proximity zone measurement of $\Delta R_{p}\approx 0.13$ pMpc. Continuum uncertainties do not play any significant role for these measurements. The lower panel of Fig. 9 shows the spectrum of this object, which does not show any BAL features. The continuum normalized spectrum and its proximity zone is shown in the second to last panel of Fig. 16 in Appendix C.

6.2. How Common are Small Proximity Zones?

In order to quantify the probability of finding objects with such small proximity zones, we compare in Fig. 10 our measured proximity zones $R_{p}$ from a subset of quasar spectra from our data sample to proximity zones from our radiative transfer simulations. We show the distribution of quasar proximity zone sizes for three different luminosity ranges, quasars with magnitudes $-25.8\leq M_{\rm 1450}\leq-24.5$ are shown in the left panel, brighter objects with magnitudes $-27.5\leq M_{\rm 1450}\leq-26.5$ are shown on the middle panel, and the right panel shows the distribution of proximity zones for very bright quasars with $M_{\rm 1450}=-29.14$. The magnitude intervals were chosen, such that each bin has a width of $\sim 1$ dex and includes one of the small proximity zone objects. This results in six quasars in the faintest magnitude range, eleven quasars in the brighter range of magnitudes and only one object at $M_{\rm 1450}\approx-29$ (Wu et al., 2015). Because of this split into magnitude intervals, we are comparing $R_{p}$ instead of $R_{p,\rm corr}$, and model the ranges of quasar luminosities in the simulations. Note that the redshift range covered by the quasars in each bin is very broad, i.e. $5.844\leq z\leq 6.4323$ for the faintest magnitude bin, $5.7722\leq z\leq 6.5412$ for the brighter magnitude range, and the only one bright quasar has a redshift $z\approx 6.3258$. We forward model the redshift evolution of the quasars in the radiative transfer simulations, such that the broad range of luminosities in each bin does not influence our results.

The red histograms show the distribution of measured proximity zones from our data set. In this figure we only consider quasars that have their redshifts determined by the detection of a Mg II, [C II] or CO line, i.e. with small redshift uncertainties. These histograms have been normalized to the number of objects available in our data set within the given magnitude interval.

The gray histograms show the distribution of proximity zone sizes expected from our radiative transfer runs. We simulate multiple realizations of each quasar in our sample matching the redshift, luminosity, and redshift uncertainties (which we add as a Gaussian uncertainty to the quasar position in our mock spectra; see Fig. 4). For the two leftmost panels of Fig. 10 we simulate $100$ skewers for each quasar in our data sample. Thus the gray histograms show the proximity zones of $600$ skewers in the left panel and $1100$ skewers in the middle panel. For the single very bright quasar in the right panel of Fig. 10 we simulate $900$ realizations to improve the statistics. We do not model continuum uncertainties in the simulations, which we expect to have negligible influence on the distribution of proximity zone sizes (see discussion in § 5.4, Kramer & Haiman, 2009). For the radiative transfer runs we assume a highly ionized IGM, which is consistent with empirical constraints on the neutral gas fraction at $z\sim 6$ (Calverley et al., 2011; Wyithe & Bolton, 2011; Becker et al., 2015; McGreer et al., 2015). Following the arguments of § 5.3 we do not expect any qualitative changes when assuming a mostly neutral IGM.

In both the faint and the brighter case (left and middle panel) the bulk of the distributions of measured and simulated proximity zones agrees very well. Note that our simulations do not reproduce the very large proximity zones of $R_{p}\sim 10-13$ pMpc that Carilli et al. (2010) measured previously.

However, the simulations also do a poor job of reproducing the frequency of very small proximity zones that we find. For the fainter quasars (left panel) our simulations indicate that the probability of finding a quasar with $R_{p}\approx 0.45$ pMpc, which is our measurement for $\rm CFHQSJ2229+1457$, in a sample of six quasar spectra is $\approx 3\%$. However, the significance of this small proximity zone is limited by the definition of the proximity zone. The smoothed flux in the spectrum of $\rm CFHQSJ2229+1457$ remains below the $<10\%$ level within a distance of $10$ pMpc to the quasar, whereas a by eye inspection of the simulated small proximity zones reveals that in most spectra the flux increases above $10\%$ again just outside of their proximity zone. If we would adopt a different definition of the proximity zone that would also be sensitive to the length of the GP trough outside of the proximity zone, this object would be an even greater outlier.

In the middle panel showing quasars with magnitudes $-27.5\leq M_{\rm 1450}\leq-26.5$ we have one object, $\rm SDSSJ1335+3533$, with a measured proximity zone $R_{p}\approx 0.78$ pMpc. In the sample of $1100$ simulated quasars with similar redshifts and magnitudes only one of them has a proximity zone size that small. Thus the probability of finding such an object in a sample of eleven quasars is $\approx 1\%$.

In the right panel we have only one object in our data set, $\rm SDSSJ0100+2802$. Its proximity zone is relatively small given its extreme brightness. We estimate the probability of finding a quasar with a proximity zone of $R_{p}\lesssim 7.12$ pMpc to be $\approx 3\%$.

In summary, in all three respective magnitude ranges, the occurrence of small proximity zones is much higher in our data sample than the simulations predict. We will now investigate several possible scenarios that could explain the exceptionally small proximity zones of the three aforementioned objects.

6.3. Possible Explanations for Small Proximity Zones

We address three scenarios that could possibly explain the exceptionally small proximity zones of the three quasars mentioned in the previous subsection: The proximity zones could be prematurely truncated due to associated dense absorbers, such as damped Ly$\alpha$ systems (DLA) or Lyman limit systems (LLS), patches of remaining neutral hydrogen within the IGM could truncate the proximity zones, or the quasars could be very young.

6.3.1 Truncation Due to Damped Ly$\alpha$ Systems or Lyman Limit Systems

One possible explanation for the exceptionally small proximity zone sizes could be the truncation of the zones due to strong absorbers, such as DLAs or LLSs, which could either be intervening or physically associated with the quasar environment (Bañados et al. in prep.). The presence of such self-shielding absorbers could prematurely truncate the proximity zone at small radii by blocking the quasar’s ionizing flux. These optically thick absorption line systems are not correctly modeled by our radiative transfer simulations (see § 4). We model only the gas densities in the IGM, but do not include such dense gas patches that could be coming from galaxy like overdensities resulting in these absorbers. Thus, proximity zones truncated prematurely by dense absorption systems are not reproduced in our simulations and need to be identified and eliminated from any comparison between the data and the simulations.

To this end, we search for signatures of strong absorption line systems in the quasar spectra near the end of their proximity zones by visually inspecting the spectra, and searching for evidence of damping wings which would indicate the presence of a strong absorber, as well as associated ionic metal-line transitions.

In the spectrum of $\rm SDSSJ1335+3533$, the only strong metal line absorber we could find was a Mg II absorber at $z\approx 2.10$. However, we have to note that the signal-to-noise ratio of the spectrum is insufficient to identify very weak metal absorption features.

Searching the ESI spectrum of $\rm CFHQSJ2229+1457$ for metal absorption lines associated with nearby dense absorption systems did not reveal any absorbers due to the very low signal-to-noise data ($\rm S/N\approx 2$ per pixel). To facilitate the search for metal lines, we obtained a higher quality spectrum with the Low Resolution Imaging Spectrometer (LRIS) at the Keck I telescope, which has a lower resolution ($R\sim 1800$) than the ESI spectrum ($R\sim 4000$) but a higher signal-to-noise ratio ($\rm S/N\approx 7$ per pixel). An excerpt of this spectrum is shown in Fig. 11 and reveals an N V doublet associated with an absorber at $z\approx 6.136$, whose H I absorption Voigt profile is shown as the red curve when assuming a column density of $N_{\rm HI}=10^{19}$ cm*-2* and a Doppler parameter of $b=40$ km/s. The absorber is located at a distance of $0.91$ pMpc from the quasar and thus it is highly unlikely that this absorber influences the size of the proximity zone significantly, because the proximity zone ends at a higher redshift of $z\approx 6.144$ at a distance $R_{p}\approx 0.45$ pMpc. Additionally, at the location of the absorber in the spectrum, there is a transmitted flux spike in the Ly$\alpha$ forest partly visible (at $\lambda\approx 8670$ Å), which indicates that the absorber cannot be saturated, implying $N_{\rm HI}\lesssim 10^{14}$ cm*-2*, i.e. its column density would need to be much less than the here assumed $N_{\rm HI}=10^{19}$ cm*-2* in order not to violate the spectrum, and is thus unlikely to truncate the proximity zone. We also do not find evidence for any low-ion absorption lines, such as Si II at $\lambda_{\rm rest}=1260.42$ Å, that should be present in DLAs.

In the spectrum of $\rm SDSSJ0100+2802$ we find a low ionization absorption system close to the object at $z\approx 6.144$ (see Fig. 11), which has previously been identified by Wu et al. (2015). However, this absorber is at such large line-of-sight distance $R\approx 9.94$ pMpc from the quasar that it cannot truncate its proximity zone, that already ends at a distance of $R_{p}=7.12\pm 0.13$ pMpc. We show the H I Voigt absorption profile (red curve) assuming an absorber with $N_{\rm HI}=10^{19}$ cm*-2* and $b=40$ km/s, and its associated lines (Si II at $\lambda_{\rm rest}=1260.42$ Å and $\lambda_{\rm rest}=1304.37$ Å, O I at $\lambda_{\rm rest}=1302.16$ Å, [C II] at $\lambda_{\rm rest}=1334.53$ Å). Note that an absorption system with a higher column density of $N_{\rm HI}\sim 10^{20.3}$ cm*-2* characteristic of DLAs, would be too broad and would be inconsistent with the presence of flux transmission that are observed in the Ly$\beta$ forest at this redshift (see § 6.3.2).

In light of the possibility that proximate absorbers could truncate the proximity zones, it is interesting to estimate the number of those proximate absorbers that one would expect to find around quasars. Therefore, we start with an estimate of intervening absorbers along the line-of-sight to quasars by Songaila & Cowie (2010), who have estimated the number density $\mathop{}\!\mathrm{d}N/\mathop{}\!\mathrm{d}z$ of intervening LLSs in Fig. $5$ and Table $3$ of their paper. They determine the best power-law fit to their measurements with a maximum likelihood analysis to be

[TABLE]

which results in an estimate of $\mathop{}\!\mathrm{d}N/\mathop{}\!\mathrm{d}z\approx 6.99$ at $z=6$.

On the one hand the abundance of LLSs could be enhanced in the proximity of quasars because of the overdense quasar environment (Hennawi et al., 2006; Prochaska et al., 2008, 2013). However, the intense radiation from the quasar will also ionize the dense gas in its surroundings making it less likely to self-shield, possibly lowering the abundance of so-called proximate LLSs (Hennawi & Prochaska, 2007; Prochaska et al., 2008). Prochaska et al. (2010) found empirically that at $z\sim 4$ the number density of proximate LLSs (occurrence within $\Delta v\leq 3000$ km s*-1* of the quasar emission redshift) roughly equals the number density of intervening LLSs (see Fig. $15$ of Prochaska et al., 2010). Assuming the same approximate equality holds at $z\sim 6$, we estimate the probability $p$ of finding a proximate LLS within $1$ pMpc to a quasar at $z=6$ to be $p\approx n(z)\Delta z\approx 0.11$. Given the size of our data set, which consists out of $30$ quasars, we would expect to find a proximate LLS in $\sim 3$ spectra.

Where are the $\sim 3$ proximate LLSs that we expect to find in our data sample? We found only one in the spectrum of the quasar $\rm SDSSJ0840+5624$, which we excluded from our analysis for this reason (see Appendix A). The absence of other proximate LLSs detected in our data sample could be attributed to a few different reasons: the estimate of proximate LLSs by Prochaska et al. (2010) assumes a $\Delta v\leq 3000$ km s*-1* window, which corresponds to a distance of $\approx 4.3$ pMpc to a quasar at $z=6$. We are just interested in the innermost $1$ pMpc to the quasar, where the quasar’s radiation is stronger and could thus have photoionized all possible LLSs (Hennawi & Prochaska, 2007). A second possible reason for the lack of detections of more proximate LLSs in our sample, could be due to the fact that our data sample at $z\sim 6$ is brighter than the one analyzed by Prochaska et al. (2010) at $z\sim 3.5$ and thus the higher radiation could have photoionized more LLSs. Hence, it might be not surprising, that we are only detecting one proximate LLS in our data sample instead of the expected $\sim 3$.

Although we cannot rule out the presence of absorbers with weak metal lines or rare metal-free absorbers (Fumagalli et al., 2011; Simcoe et al., 2012; Cooper et al., 2015; Crighton et al., 2016; Cooke et al., 2017) in all objects, which would require higher quality data, we do not see evidence for proximate LLSs in the three spectra showing small proximity zones, particularly in the spectrum of the very bright quasar $\rm SDSSJ0100+2802$, where we would most likely be able to detect it. Follow up higher $\rm S/N$ and resolution observations are required to confirm the lack of associated dense absorbers in these spectra.

6.3.2 Neutral Islands of Gas in a Patchily Reionized IGM

Another possible explanation for the existence of very small proximity zones could be the truncation of the proximity zones due to islands of very neutral gas in an otherwise highly ionized IGM, which can occur in patchy reionization models (e.g. Lidz et al., 2007; Mesinger, 2010). This scenario implies that reionization is not yet complete at $z\sim 6$ and these quasars lie in (presumably rare) neutral environments that could absorb all incoming flux and thus truncate the proximity zones. Opposite to the previously described scenario where LLSs, i.e. dense gas on galactic scales, truncate the proximity zones, we assume here neutral patches of $\sim 10$ pMpc scales that are associated with reionization and do not have any low redshift counter parts. This explanation requires the presence of large patches of highly neutral gas just outside the proximity zone, which would also absorb flux in the Lyman-series (Ly$\alpha$, Ly$\beta$, etc.) forests. In other words, in a patchy reionization scenarios one expects to see long GP troughs that are line black in both Ly$\alpha$ and Ly$\beta$, beginning right at the edge of the proximity zone.

In Fig. 13 we show continuum normalized spectra of the Ly$\alpha$ and the Ly$\beta$ forest (top and bottom panels, respectively) near the three quasars in question. The black dashed lines in all panels indicate the extent of the quasar proximity zone. The left panels show the Ly$\alpha$ and the Ly$\beta$ forest close to the quasar $\rm CFHQSJ2229+1457$. While we do not detect any transmitted flux in the Ly$\beta$ forest, we do see a few transmission spikes in the Ly$\alpha$ forest at $R\approx 1.1$ pMpc and $R\approx 2.8$ pMpc, indicated by the red dashed lines. The transmission in the Ly$\beta$ forest at the same distances must have been absorbed due to foreground Ly$\alpha$ absorption at $z\approx 5.02$ and $z\approx 4.99$.

The middle panels show the Ly$\alpha$ and Ly$\beta$ forest of $\rm SDSSJ1335+3533$, where we see transmitted flux in the Ly$\alpha$ forest just outside of the proximity zone, i.e. at $R\approx 0.9$ pMpc and $R\approx 1.0$ pMpc, and a bit further away at $R\approx 2.9$ pMpc and $R\approx 3.6$ pMpc. In the Ly$\beta$ forest we can see clear transmission spikes just outside of the proximity zone at $R\approx 1.4$ pMpc, $R\approx 2.8$ pMpc and $R\approx 3.3$ pMpc.

The right panels show excerpts of the Ly$\alpha$ and Ly$\beta$ forest for $\rm SDSSJ0100+2802$. One can clearly see several transmitted flux spikes in the Ly$\alpha$ forest right outside the proximity zone. Also the Ly$\beta$ forest shows a few flux spikes close to the proximity zone, although less prominent ones.

These detections of transmitted flux just outside the proximity zones in the Ly$\alpha$ or Ly$\beta$ forest suggest that if there are patches of neutral gas present close to the quasar, they need to be $\lesssim 1$ pMpc in size, corresponding to the size of the regions that are both line black in the Ly$\alpha$ and Ly$\beta$ forest. Although we cannot completely rule out a patchy reionization scenario with such small islands of neutral gas, our constraints on their sizes already suggest that they would be difficult to distinguish from metal poor LLSs. Deeper spectra (at least for the quasars $\rm CFHQSJ2229+1457$ and $\rm SDSSJ1335+3533$) covering their proximity zones further down the Lyman series could reveal additional Ly$\beta$ or Ly$\gamma$ transmission spikes that would further constrain this scenario, as well as enabling more sensitive searches for metal lines to rule out an LLS scenario.

6.3.3 Young Quasar Age

Finally, another possibility to explain small proximity zones is that the quasars may have turned on only very recently. When a quasar turns on the IGM requires a finite amount of time to adjust to the new higher photoionization rate. We will focus the discussion in this section on a highly ionized surrounding IGM, consistent with optical depth measurements at $z\sim 6$ (e.g. Wyithe & Bolton, 2011; Calverley et al., 2011; Becker et al., 2015; McGreer et al., 2015), although we note similar arguments apply in the neutral case (see Khrykin et al., 2016, for further details). Following (Khrykin et al., 2016)121212Khrykin et al. (2016) wrote down the analogous equation for the He II singly ionized fraction in the context of He II Ly$\alpha$ proximity zones at $z\sim 3$, the exact same set of arguments applies here for neutral hydrogen around $z\sim 6$ quasars. the time evolution of the neutral fraction in the proximity zone is well described by the equation

[TABLE]

where $t_{\rm Q}$ is the quasar age, $x_{\rm HI,0}\approx n_{\rm e}\alpha_{\rm HII}/\Gamma_{\rm UVB}$ is the neutral fraction of the IGM before the quasar turned on, and $x_{\rm HI,eq}\approx n_{\rm e}\alpha_{\rm HII}/(\Gamma_{\rm UVB}+\Gamma_{\rm QSO})$ is the new lower neutral fraction that the IGM will reach once ionization equilibrium is attained. Here $n_{\rm e}$ and $\alpha_{\rm HII}$ are the IGM electron density and recombination coefficient, respectively. The characteristic timescale to reach ionization equilibrium, the equilibration timescale $t_{\rm eq}$, is

[TABLE]

Eqn. (12) thus implies that a quasar has to have been shining for at least $t_{\rm Q}\gtrsim t_{\rm eq}$, in order for the surrounding IGM to have reached ionization equilibrium with the quasar and result in the maximum proximity zone size. Ionizing equilibrium is achieved more rapidly close the quasar due to the stronger radiation field ($\Gamma_{\rm QSO}\propto R^{-2}$). In order to obtain an optical depth of $\tau_{\rm lim}=2.3$ at $z\sim 6$, which corresponds to a transmitted flux level of $10\%$ in the Ly$\alpha$ forest, a value of $\Gamma_{\rm QSO}\sim 10^{-12}$ $\rm s^{-1}$ is implied (see § 5.3), resulting in a typical equilibration time scale of $t_{\rm eq}\sim 3\times 10^{4}$ yr. Full ionization equilibrium is reached after a few equilibration times, or for $t_{\rm Q}\sim 10^{5}$ yr. The region in which ionization equilibrium is reached, and hence also the size of the quasar’s proximity zone, will grow with time until $t_{\rm Q}\approx 3t_{\rm eq}$. Thus, even for a highly ionized IGM, the sizes of the proximity zones are dependent on the quasar age for $t_{\rm Q}\lesssim 10^{5}$ yr, whereas for ages $t_{\rm Q}\gtrsim 10^{5}$ yr the proximity zone sizes are largely independent on the precise quasar age131313A mild increase in proximity zone sizes for long quasar ages, i.e. $t_{\rm Q}\gtrsim 10^{7}$ yr, can be attributed to heating effects due to He II reionization (Bolton et al., 2010, 2012; Davies et al., in prep.).

In Fig. 14 we show the dependence of proximity zone size on quasar age determined from radiative transfer simulations in a highly ionized IGM. Quasars with magnitudes $M_{\rm 1450}=-25.0$ (left) $M_{\rm 1450}=-27.0$ (middle) and $M_{\rm 1450}=-29.0$ (right) are shown, where we have simulated $900$ proximity zones at each magnitude. The black horizontal dashed lines show the measured sizes of the proximity zones of $\rm CFHQSJ2229+1457$ (left panel), $\rm SDSSJ1335+3533$ (middle panel) and $\rm SDSSJ0100+2802$ (right panel).

The distribution of proximity zone sizes in the left and middle panel encompasses values of $R_{p}\approx 1$ pMpc that we have measured only when assuming a short quasar age of $t_{\rm Q}\lesssim 10^{5}$ yr. Hence a plausible explanation for the exceptionally small sizes of the proximity zones could be young ages for these quasars. The proximity zone of the very bright quasar $\rm SDSSJ0100+2802$ (right panel) is not as significant of an outlier as the other two objects, but nevertheless the small size of its zone would be more probable for $t_{\rm Q}\sim 10^{5}$ yr.

At lower redshifts $z\lesssim 4$, studies of quasar clustering imply that the duty cycle of quasar activity is $t_{\rm dc}\sim 10^{9}$ yr (Shen et al., 2007; White et al., 2008)141414Whether this also applies at redshifts $z\sim 6$ still remains to be determined. . However, whereas clustering constraints the duty cycle, proximity zones actually probe a different timescale, which is the duration of quasar emission episodes which could be considerably shorter. Indeed, we actually obtain a lower limit on the episodic lifetime, because the quasar could continue to emit for many years after we observe it on Earth. Nevertheless, for a given episodic lifetime $t_{\rm episodic}$ the probability of measuring an age of $t_{\rm Q}$ is $p=t_{\rm Q}/t_{\rm episodic}$, assuming the simplest “light-bulb” lightcurve for the quasars. If we assume an average episodic lifetime of $t_{\rm episodic}\sim 10^{8}$ yr, the probability of observing a quasar that has only been shining for $\sim 10^{5}$ yr is $p\sim 0.1\%$, or $p\sim 1\%$ for $t_{\rm episodic}\sim 10^{7}$ yr. We have discovered three objects suggesting ages of $t_{\rm Q}\sim 10^{5}\,{\rm yr}$ in a sample of $30$ quasars, i.e. $p\approx 10\%$. To be consistent with finding a few of these small proximity zones, the quasar episodic lifetimes would need to be $t_{\rm episodic}\sim 10^{6}$ yr. This would leave the sizes of the proximity zones of the vast majority of quasars unchanged, which we typically observe much later in their evolution, but could explain the very small zones we find.

7. Summary

In this paper we analyze a sample of $34$ high redshift quasar spectra taken with the ESI instrument on the Keck II telescope. We reduce the spectra in a homogeneous way, and analyze the sizes of their proximity zones for a subset of $30$ quasars which do not exhibit BAL features or have obvious signatures of nearby associated absorption line systems that could prematurely truncate their proximity zones. Our analysis uses updated redshift measurements, fully consistent values for the quasar absolute magnitudes $M_{\rm 1450}$, and carefully fits to the quasar continua based on principal component spectra.

Previous work found a strong evolution of proximity zone sizes with redshift, and it was argued that this provided evidence for rapid evolution of the IGM neutral fraction during the epoch of reionization. We instead find a much shallower redshift evolution, which is however consistent with the evolution predicted by our radiative transfer simulations, irrespective of assumptions about the ionization state of the IGM. The size of the proximity zone ends at a distance corresponding to a flux transmission level of $10\%$ according to the definition of quasar proximity zones (Fan et al., 2006b) and, in a highly ionized IGM, this distance is reached at a location where the ionization rate of the quasar dominates the ionization rate of the UVB by at least an order of magnitude. As such, the size of the quasar proximity zone is essentially independent of the UVB and neutral gas fraction. Assuming a neutral ambient IGM the observed proximity zone sizes $R_{p}$ track the growth of the ionized H II region $R_{\rm ion}$, but then cease to grow when the distance with a $10\%$ flux transmission level is reached. For a highly plausible quasar age of $t_{\rm Q}\sim 10^{7.5}$ yr, $R_{p}$ is thus independent of the ionization state of the IGM for a mostly neutral IGM as well. Thus contrary to previous claims, both the observed shallow redshift evolution and the results from our simulations imply that the redshift evolution of proximity zone size $R_{p}$ does not significantly constrain the evolution of the neutral gas fraction during the epoch of reionization.

Our analysis uncovered three quasars with exceptionally small proximity zones, two of them with $R_{p}\lesssim 1$ pMpc. The occurrence of zones this small is extremely rare in our simulations. We estimate the probability of finding these objects in a sample of $30$ quasar spectra to be $p\approx 1\%$ for $\rm SDSSJ1335+3533$, and $p\approx 3\%$ for $\rm CFHQSJ2229+1457$ and $\rm SDSSJ0100+2802$. We search for evidence of proximate absorption systems such as DLAs or LLSs, or islands of remaining neutral hydrogen associated with patchy reionization, both of which could result in truncation of the proximity zones closer to the quasars. However the absence of strong metal absorption lines or large GP troughs (in both the Ly$\alpha$ and the Ly$\beta$ forest) appears to disfavor these truncation scenarios. Nevertheless we cannot completely rule out the presence of low metallicity or metal-free absorbers due to the low signal-to-noise ratio of the data of $\rm CFHQSJ2229+1457$ and $\rm SDSSJ1335+3533$.

The most plausible explanation for the small proximity zones that we observe is that these quasars are relatively young. It takes the IGM roughly $\sim 10^{5}\,{\rm yr}$ to reach ionization equilibrium with the quasar radiation. Proximity zones of quasars with ages longer than this timescale are essentially independent of the exact age of the quasar. However, for shorter quasar ages the surrounding IGM has not yet reached ionization equilibrium resulting in proximity zones comparable to the small sizes we observe. If the duration of quasar emission episodes is $t_{\rm episodic}\sim 10^{8}$ yr, the detection of these small zones would be very unlikely, i.e. $p\sim 0.1\%$. This discrepancy can be resolved if one assumes a shorter duration of $t_{\rm episodic}\sim 10^{6}$ yr, resulting in a probability of $p\sim 10\%$ for finding these objects with ages of $t_{\rm Q}\sim 10^{5}$ yr.

However, an average episodic lifetime of $t_{\rm episodic}\sim 10^{6}$ yr causes significant tension with the inferred sizes of SMBHs at these redshifts, since the presence of $\sim 10^{9}M_{\odot}$ SMBHs at $z\sim 6$ requires that quasars accrete continuously at around the Eddington limit for nearly the entire Hubble time of $\sim 10^{9}{\rm yr}$ (Volonteri, 2012). Thus, although the proximity zone measurements constrain the episodic lifetime $t_{\rm episodic}$, whereas the sizes of SMBH depend on the integral over all emission episodes, i.e. the duty cycle, the presence of large SMBHs at $z\sim 6$ either requires long episodic lifetimes comparable to the duty cycle $\sim 10^{9}\,{\rm yr}$, or that the quiescent time between episodic emission bursts must be very short, implying that SMBHs grow via many short episodic phases (Schawinski et al., 2015; Mao & Kim, 2016; Davies et al., in prep.).

It would clearly be interesting to uncover more quasars with small proximity zones, particularly if they are indeed young quasars that have only been shining for $t_{\rm Q}\sim 10^{5}\,{\rm yr}$. In the future, we plan to use large quasar surveys, such as Pan-STARRS or SDSS, to identify candidate small-zone quasars and measure accurate redshifts from [C II] or CO emission lines with the Atacama Large Millimeter Array (ALMA), in order to confirm the potentially small sizes of their proximity zones. Deep optical and near-IR follow-up observations will be helpful to rule out a premature truncation due to associated absorption systems. Further theoretical work is also required to investigate the causes and implications of small proximity zones, and their relationship to the distribution of quasar ages, which is a subject we address in future work (Davies et al., in prep.).

Acknowledgment

The authors would like to thank J. Bolton and M. Haehnelt for their questions about the small proximity zones in our sample at the “Dark Ages” meeting in Heidelberg, which provided part of the motivation for this work. Additionally, we wish to thank E. Bañados, X. Fan, B. Venemans, and F. Walter for valuable input and discussion.

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 research has made use of the Keck Observatory Archive (KOA), which is operated by the W. M. Keck Observatory and the NASA Exoplanet Science Institute (NExScI), under contract with the National Aeronautics and Space Administration.

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

ZL was supported by the Scientific Discovery through Advanced Computing (SciDAC) program funded by U.S. Department of Energy Office of Advanced Scientific Computing Research (ASCR) and the Office of High Energy Physics. Calculations presented in this paper used resources of the National Energy Research Scientific Computing Center (NERSC), which is supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231. The authors would like to thank Dmitriy Morozov and Gunther H. Weber for the halo finder work, which was funded by the ASCR project “Scalable Analysis Methods and In Situ Infrastructure for Extreme Scale Knowledge Discovery,” program manager Lucy Nowell.

Appendix A Associated Absorbers in SDSSJ0840+5624

This object has been excluded from our analysis because we found several absorption systems at high redshift associated with the quasar itself. These absorbers cause additional absorption within the proximity zone of the quasar, such that our analysis recovers a spuriously small proximity zone of $R_{p}\approx 0.88$ pMpc. However, inspecting the spectrum by eye reveals that the actual proximity zone of this object extends much further out to $R_{p}\gtrsim 5.0$ pMpc.

We find two N V absorption line systems ($\lambda_{\rm rest}=1238.82$ Å and $\lambda_{\rm rest}=1242.80$ Å) associated with close absorbers at $z\approx 5.835$ and $z\approx 5.825$, which are shown in Fig. 16. We also find another absorption system in the spectrum at $z\approx 5.594$ with associated Si II absorption ($\lambda_{\rm rest}=1260.42$ Å) that falls within the proximity zone. Thus there is a lot of absorption within the proximity zone that cannot be contributed to residual neutral hydrogen only and an exact determination of the proximity zone for this object without contamination from other absorption systems is impossible.

Appendix B Details of the Continuum Modeling of Each Quasar

Table 3 shows the details of modeling the continuum of each quasar spectrum.

Appendix C Zones of Remaining Quasars from our Data Set

Fig. 16 shows the same as Fig. 2 and Fig. 3 for the remaining quasars in our data sample sorted by the redshift of the quasar.