Time-dependent behaviour of quasar proximity zones at $z \sim 6$

TL;DR

This study uses radiative transfer simulations and a semi-analytic model to analyze the time-dependent behavior of quasar proximity zones at redshift around 6, providing insights into quasar activity timescales and growth of supermassive black holes.

Contribution

It introduces a semi-analytic model that captures the equilibrium and non-equilibrium behavior of proximity zones and links quasar variability to observed zone sizes.

Findings

Simulated proximity zone sizes match observations well.

Large luminosity variations on timescales less than 10^4 years are unlikely.

Proximity zones can be explained by steady emission models without requiring rapid variability.

Abstract

Since the discovery of $z\sim 6$ quasars two decades ago, studies of their Ly$\alpha$-transparent proximity zones have largely focused on their utility as a probe of cosmic reionization. But even when in a highly ionized intergalactic medium, these zones provide a rich laboratory for determining the timescales that govern quasar activity and the concomitant growth of their supermassive black holes. In this work, we use a suite of 1D radiative transfer simulations of quasar proximity zones to explore their time-dependent behaviour for activity timescales from $\sim10^3$ to $10^8$ years. The sizes of the simulated proximity zones, as quantified by the distance at which the smoothed Ly$\alpha$ transmission drops below 10% (denoted $R_p$), are in excellent agreement with observations, with the exception of a handful of particularly small zones that have been attributed to extremely short $\lesssim 10^4$ year lifetimes. We develop a physically motivated semi-analytic model of proximity zones which captures the bulk of their equilibrium and non-equilibrium behaviour, and use this model to investigate how quasar variability on $\lesssim10^5$ year timescales is imprinted on the distribution of observed proximity zone sizes. We show that large variations in the ionizing luminosity of quasars on timescales of $\lesssim10^4$ years are disfavored based on the good agreement between the observed distribution of $R_p$ and our model prediction based on "lightbulb" (i.e. steady constant emission) light curves.

Equations

$$R_{\mathrm{ion}}={\left(\frac{3{\dot{N}}_{\mathrm{ion}}{t}_{\mathrm{q}}}{4\pi n_{\mathrm{H}}{x}_{\mathrm{HI}}}\right)}^{1/3},$$

$${\tau }_{\mathrm{GP}}=\frac{{\sigma }_{\alpha }c{n}_{\mathrm{H}}{x}_{\mathrm{HI}}}{H(z)},$$

$$R_p=\frac{3.14\mathrm{proper}\mathrm{Mpc}}{{\Delta }_{\lim }}$$

$$\frac{d{n}_{\mathrm{HII}}}{dt}$$

$$\frac{d{n}_{\mathrm{HeII}}}{dt}$$

$$\frac{d{n}_{\mathrm{HeIII}}}{dt}$$

$$n_{\mathrm{HI}}$$

$$n_{\mathrm{HeI}}$$

$$n_e$$

$$\frac{dT}{dt}=\frac{({\gamma }_{\mathrm{ad}}-1)\mu m_{\mathrm{H}}}{k_{\mathrm{B}}\rho }(\mathcal{H}-\Lambda )-2H(z)T-\frac{T}{n}\frac{d{n}_e}{dt},$$

$$R_{\mathrm{ion}}^{\mathrm{He}}={\left(\frac{3{\dot{N}}_{\mathrm{ion}}^{\mathrm{He}\mathrm{II}}{t}_{\mathrm{q}}}{4\pi n_{\mathrm{He}}{x}_{\mathrm{He}\mathrm{II}}}\right)}^{1/3},$$

$${\tau }_{\mathrm{eff}}(r)={\tau }_{\mathrm{b}}{\left(\frac{{\Gamma }_{\mathrm{q}}(r)+{\Gamma }_{\mathrm{b}}}{{\Gamma }_{\mathrm{b}}}\right)}^{-\alpha },$$

$$r_{\mathrm{b}}$$

$${\tau }_{\mathrm{eff}}(r)={\tau }_{\mathrm{b}}{\left[1+{\left(\frac{r}{r_{\mathrm{b}}}\right)}^{-2}\right]}^{-\alpha }.$$

$$R_p=r_{\mathrm{b}}{\left[{\left(\frac{{\tau }_{\mathrm{b}}}{{\tau }_{\lim }}\right)}^{1/\alpha }-1\right]}^{-1/2}.$$

$$R_p\approx r_{\mathrm{b}}({\tau }_{\mathrm{b}}/{\tau }_{\lim }{)}^{-1/2\alpha }\propto r_{\mathrm{b}}{\Gamma }_{\mathrm{b}}^{1/2}\propto {\Gamma }_{\mathrm{b}}^0,$$

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

$$\frac{x_{\mathrm{HI}}(t)}{x_{\mathrm{HI},0}}=1+\left(\frac{x_{\mathrm{HI},\mathrm{eq}}}{x_{\mathrm{HI},0}}-1\right)\left(1-e^{-t/t_{\mathrm{eq}}}\right).$$

$${\left(\frac{{\tau }_{\mathrm{eff}}(t)}{{\tau }_{\mathrm{b}}}\right)}^{1/\alpha }=1+\left(\frac{{\Gamma }_{\mathrm{b}}}{{\Gamma }_{\mathrm{HI}}}-1\right)\left(1-e^{-t/t_{\mathrm{eq}}}\right).$$

$${\tau }_{\mathrm{eff}}(r,t)={\tau }_{\mathrm{b}}{\left[1-\left(\frac{{\Gamma }_{\mathrm{q}}(r)}{{\Gamma }_{\mathrm{q}}(r)+{\Gamma }_{\mathrm{b}}}\right)\left(1-e^{-t[{\Gamma }_{\mathrm{q}}(r)+{\Gamma }_{\mathrm{b}}]}\right)\right]}^{\alpha },$$

$${\tau }_{\mathrm{eff}}(r,t)\approx {\tau }_{\mathrm{b}}{e}^{-\alpha t{\Gamma }_{\mathrm{q}}(r)},$$

$$R_p\approx R_{p,\mathrm{eq}}{\left(\frac{\alpha t_{\mathrm{q}}{\Gamma }_{\mathrm{q}}(R_{p,\mathrm{eq}})}{\ln {\tau }_{\mathrm{b}}-\ln {\tau }_{\lim }}\right)}^{1/2},$$

$$\frac{x_{\mathrm{HI}}(t)}{x_{\mathrm{HI},\mathrm{eq}}^{\mathrm{on}}}=1+\left(\frac{x_{\mathrm{HI},\mathrm{b}}}{x_{\mathrm{HI},\mathrm{eq}}^{\mathrm{on}}}-1\right)\left(1-e^{-t{\Gamma }_{\mathrm{b}}}\right),$$

$$\frac{x_{\mathrm{HI}}(r,t)}{x_{\mathrm{HI}}(r,t_{\mathrm{on}})}=1+\left(\frac{x_{\mathrm{HI},\mathrm{eq}}^{\mathrm{on}}(r)}{x_{\mathrm{HI}}(r,t_{\mathrm{on}})}-1\right)\left(1-e^{-(t-t_{\mathrm{on}})({\Gamma }_{\mathrm{q}}(r)+{\Gamma }_{\mathrm{b}})}\right),$$

$$\frac{x_{\mathrm{HI}}(r,t)}{x_{\mathrm{HI}}(r,t_{\mathrm{off}})}=1+\left(\frac{x_{\mathrm{HI},\mathrm{b}}}{x_{\mathrm{HI}}(r,t_{\mathrm{off}})}-1\right)\left(1-e^{-(t-t_{\mathrm{off}}){\Gamma }_{\mathrm{b}}}\right),$$

$$x_{\mathrm{HI}}(t)=x_{\mathrm{HI},0}{e}^{-{\int }_0^t\Gamma (t^{\prime })d{t}^{\prime }}+{\alpha }_{\mathrm{HII}}{n}_e{\int }_0^t{e}^{-{\int }_{t^{\prime }}^t\Gamma (t^{\prime \prime })d{t}^{\prime \prime }}d{t}^{\prime },$$

$$\displaystyle\tau_{\rm eff}(r,t)=\tau_{\rm eff}(r,0)\Big{[}e^{-\int_{0}^{t}\Gamma_{\rm HI}(r,t^{\prime})dt^{\prime}}$$

$$\displaystyle+\Gamma_{\rm HI}(r,0)\int_{0}^{t}e^{-\int_{t^{\prime}}^{t}\Gamma_{\rm HI}(r,t^{\prime\prime})dt^{\prime\prime}}dt^{\prime}\Big{]}^{\alpha},$$

$$P(R_p^{\mathrm{alt}}=r_i)=\left(\prod _{r_j<r_i}P({\tau }_{\mathrm{eff},j}<{\tau }_{\lim })\right)\times P({\tau }_{\mathrm{eff},i}>{\tau }_{\lim }),$$

$$P(R_p^{\mathrm{alt}}=r_i)$$

Full text

††pubyear: 2019††pagerange: Time-dependent behaviour of quasar proximity zones at $z\sim 6$–A

1 Introduction

The “proximity effect" is the well-known tendency for decreased Ly$\alpha$ forest absorption in the intergalactic medium (IGM) close to luminous quasars along the line of sight (e.g. Bajtlik et al. 1988). This decrease in absorption is due to the strong ionizing radiation field produced by the quasar itself, which can greatly outshine the intergalactic ionizing background at small physical separations.

A related effect was predicted for quasars observed during the epoch of reionization (e.g. Shapiro & Giroux 1987; Madau & Rees 2000; Cen & Haiman 2000, and see Hogan et al. 1997 for a He II analogy), wherein the cumulative ionizing photon output from a luminous quasar carves out a transparent H II region along the line of sight. The size of this ionized region, in the absence of recombinations, is given by

[TABLE]

where $\dot{N}_{\rm ion}$ is the emission rate of ionizing photons, $t_{\rm q}$ is the age of the quasar, $n_{\rm H}$ is the (average) number density of hydrogen atoms, and $x_{\rm HI}$ is the neutral fraction of the IGM. If luminous quasars shine for $>10^{6}$ years, $R_{\rm ion}$ can reach scales of several proper Mpc, corresponding to a few thousand km s*-1* along the line of sight in the $z\gtrsim 6$ Ly$\alpha$ forest.

Observations of transparent proximity zones around quasars at $z\ga 6$ terminated by opaque Gunn-Peterson troughs (Becker et al., 2001; White et al., 2003) appeared to confirm the picture of quasar-ionized bubbles in a neutral IGM, and led to claims of constraints on the reionization history of the Universe (Wyithe & Loeb, 2004) and measurement of the evolving residual IGM neutral hydrogen fraction after reionization was complete (Fan et al., 2006; Carilli et al., 2010). The properties of the first large sample of $z\ga 6$ quasar spectra compiled by Fan et al. (2006) led them to define the size of the proximity zone, denoted by $R_{p}$, as the first location where the Ly$\alpha$ forest spectrum drops below 10% transmitted flux when smoothed to 20 Å in the observed frame, roughly corresponding to a physical scale of 1 proper Mpc at $z\sim 6$.

These transparent proximity zones observed at $z\ga 6$ can be interpreted in the context of the ionized bubbles expanding around quasars during the epoch of reionization, as described above. In this model, the IGM interior to $R_{\rm ion}$ is assumed to be transparent while the external IGM is fully opaque, thus $R_{p}$ is identified with the ionization front radius $R_{\rm ion}$. The scaling of $R_{p}$ with various parameters can then be read directly from equation (1): $R_{p}\propto t_{\rm q}^{1/3}$, $R_{p}\propto\dot{N}_{\rm ion}^{1/3}$, and $R_{p}\propto x_{\rm HI}^{-1/3}$. The $\dot{N}_{\rm ion}$ scaling in particular has been used to “correct" measurements of $R_{p}$ onto a fixed quasar luminosity scale (Fan et al., 2006; Carilli et al., 2010; Venemans et al., 2015).

However, as noted by Bolton & Haehnelt (2007), the Ly$\alpha$ forest at $z\sim 6$ becomes opaque at neutral fractions as low as $x_{\rm HI}\sim 10^{-4}$, and thus the observed proximity zone (as defined by $R_{p}$) may be truncated well before any actual ionization front (i.e. at distances less than $R_{\rm ion}$). Indeed, Bolton & Haehnelt (2007) found that if the IGM is already highly ionized, consistent with measurements of the mean transmitted flux in the IGM at $z\sim 6$, the proximity zone sizes will be roughly $R_{p}\sim 5$–$10$ proper Mpc, in agreement with the existing observations which had been previously interpreted in terms of the location of $R_{\rm ion}$ in a neutral IGM. In this highly ionized regime, $R_{p}$ is no longer directly connected to reionization. Instead, $R_{p}$ reflects the distance at which the enhanced ionizing flux from the quasar (combined with the ionizing background radiation) becomes weak enough for Ly$\alpha$ transmission to fall below 10%.

In the Bolton & Haehnelt (2007) model, the Ly$\alpha$ transmission of the IGM is approximated by the Gunn-Peterson (GP) optical depth (e.g. Weinberg et al. 1997),

[TABLE]

where $\sigma_{\alpha}$ is the Ly$\alpha$ scattering cross-section, $c$ is the speed of light, and $H(z)$ is the Hubble parameter. Allowing $x_{\rm HI}$ to vary along the line of sight due to the enhanced photoionization rate of the quasar, the observed proximity zone size, assuming a transmission threshold of 10%, is given by

[TABLE]

where $T$ is the temperature of the IGM and $\Delta_{\rm lim}$ is an “effective" IGM overdensity which represents an unknown renormalization from the GP optical depth to the true effective optical depth of the Ly$\alpha$ forest implicit in the expression. In this case, $R_{p}$ scales as $\dot{N}_{\rm ion}^{1/2}$, but is otherwise independent of the residual IGM neutral fraction as the quasar is assumed to dominate over the ionizing background. Note that the IGM is assumed to be in ionization equilibrium, so $R_{p}$ is also independent of the quasar lifetime. In this model, however, it is difficult to explain the rapid evolution in $R_{p}$ measured by Fan et al. (2006) and Carilli et al. (2010) (see also Venemans et al. 2015).

For quasars at $z\sim 6$, it is quite likely that the IGM around quasar hosts is ionized before the quasar turns on. The neutral fraction of the Universe at $z\sim 6$ is $\la 0.2$ (McGreer et al., 2015) and the standard inside-out model of reionization (Furlanetto et al., 2004) predicts that the large-scale environment of massive dark matter halos was ionized relatively early (Alvarez & Abel, 2007; Lidz et al., 2007). It is worth noting, however, that the ionized IGM regime for $R_{p}$ will still apply to a partially neutral IGM if the pre-existing ionized region is significantly larger than what $R_{p}$ would have been in the absence of neutral gas (Bolton & Haehnelt, 2007; Maselli et al., 2009). In the following we will assume that, for the purposes of analyzing quasar proximity zones via $R_{p}$, the Universe is fully reionized.

Recently, Eilers et al. (2017, henceforth E17) compiled a large number of high-redshift quasar spectra from the Keck Observatory Archive and new Keck observations to perform a homogeneous measurement of $R_{p}$ with higher signal-to-noise spectra, consistently defined $M_{1450}$, and principal component analysis continuum estimation. From their sample of 31 quasars covering $5.77\leq z\leq 6.54$, they measured a very shallow evolution of $R_{p}$ with redshift, inconsistent with previous works but consistent with radiative transfer simulations of quasars in the post-reionization IGM. They also discovered a handful of quasars with much smaller proximity zones than predicted, and suggested that these small zones could be due to non-equilibrium ionization of the IGM close to the quasar, requiring the duration of current quasar activity to be shorter than $\sim 10^{4-5}$ years (see also Eilers et al. 2018). Additional evidence for non-equilibrium effects in quasar proximity zones has also recently been measured in the He II Ly$\alpha$ forest (Khrykin et al., 2019). However, this short timescale behaviour was only studied in the context of a standard “lightbulb" model of quasar activity, wherein the quasar switched on at some time in the past and has maintained the same luminosity since then.

In this work, we expand upon previous explorations of $z\ga 6$ quasar proximity zones from a theoretical perspective, in light of the observations in E17. In § 2, we describe our 1D radiative transfer simulations and summarize the physical and observational picture of quasar proximity zones as modeled by skewers through a large-volume high-resolution hydrodynamical simulation. In § 3, using the behaviour of the Ly$\alpha$ forest in our hydrodynamical simulation as a guide, we build a new semi-analytic model for the Ly$\alpha$ transmission profile in quasar proximity zones in a highly ionized IGM. In § 4 we explore how quasar variability on timescales less than $10^{5}$ years can be imprinted on the proximity zone profile. Finally, we conclude in § 5.

We assume a $\Lambda$CDM cosmology with $\Omega_{m}=0.3$, $\Omega_{\Lambda}=0.7$, $\Omega_{\rm b}=0.047$, and $h=0.685$, consistent with the CMB constraints from Planck Collaboration et al. (2016).

2 Radiative Transfer modelling of Proximity Zones

Our radiative transfer modelling of quasar proximity zones closely follows the methods described in Davies et al. (2016) and E17, which we summarize briefly below.

2.1 Hydrodynamical simulation

We use hydrodynamical simulations skewers drawn from a Nyx Eulerian grid hydrodynamical simulation (Almgren et al., 2013; Lukić et al., 2015) 100 Mpc$/h$ on a side with $4096^{3}$ dark matter particles and baryon grid elements. We draw $900$ skewers from the $z=5.5$, $6.0$, and $6.5$ outputs of the simulation starting from the centers of the $150$ most massive dark matter halos at each redshift ($M_{h}\ga 10^{11.5}$ M⊙). While the Eulerian grid does not resolve the very dense regions ($\la r_{\rm vir}$) associated with the small-scale environment ($R<R_{\rm vir}$) of these massive halos particularly well, the lower-density IGM that gives rise to transmission in the Ly$\alpha$ forest on large scales is well-converged (J. Oñorbe, priv. comm.), and the rarity of proximate neutral absorbers in lower redshift quasars (e.g. Prochaska et al. 2008) suggests that the material on such scales does not usually affect the proximity zone (although notable exceptions exist, see D’Odorico et al. 2018; Bañados et al. 2019).

The reionization and thermal history of the hydrodynamical simulation was set by the (Haardt & Madau, 2012) model for the evolution of the ionizing background, leading to reionization at $z>10$ with an associated heat input of $\sim 10^{4}$ K (Lukić et al., 2015; Oñorbe et al., 2017). The resulting thermal state of the IGM at $z=6$ is well-described by a power law relationship between temperature ($T$) and overdensity ($\Delta\equiv\rho/\bar{\rho}$) with $T=T_{0}\Delta^{\gamma-1}$, $T_{0}\sim 7500$ K, and $\gamma\sim 1.5$. The thermal state of the ambient IGM sets the normalization of the relationship between the ionization rate and the mean transmission through the Ly$\alpha$ forest, which we discuss later in § 3.1.

2.2 1D radiative transfer

To compute the effect of quasar ionizing radiation on the IGM, we perform 1D radiative transfer (RT) along the simulations skewers described above. We use the 1D RT code developed in Davies et al. (2016), originally based on the method described in Bolton & Haehnelt (2007), with some minor adjustments as described in E17. We summarize the basic model below, but the code is described in much more detail in Davies et al. (2016).

The code solves the following time-dependent equations for ionized H and He species,

[TABLE]

where $n_{i}$ are the number densities of species $i$, $\Gamma_{i}$ are the photoionization rates of species $i$, $\Gamma^{e}_{i}$ are the collisional ionization rates from Theuns et al. (1998), and $\alpha^{\rm A}_{i}$ are the Case A recombination rates of species $i$ from Hui & Gnedin (1997). The photoionization rates include secondary ionizations by energetic photoelectrons (Shull & van Steenberg, 1985), for which we use the results of Furlanetto & Johnson Stoever (2010). The quasar SED is assumed to follow the Lusso et al. (2015) radio-quiet SED, allowing one to convert from $M_{1450}$ to the specific luminosity $L_{\nu}$ at the hydrogen ionizing edge, beyond which we assume a power law spectrum with $L_{\nu}\propto\nu^{-1.7}$. The neutral species and the electron density are recovered by the closing conditions

[TABLE]

Finally, the gas temperature $T$ is computed by solving

[TABLE]

where $\gamma_{\rm ad}=5/3$ is the adiabatic index, $\mu$ is the mean molecular weight, $\rho$ is the gas density, $\mathcal{H}$ is the heating rate, $\Lambda$ is the cooling rate, and $n$ is the total number density of all species. The first term represents the balance between the various gas-phase heating and cooling rates, the second term represents the adiabatic cooling due to the expansion of the Universe, and the last term evenly distributes thermal energy between species as the number of particles changes.

For simplicity we assume that $z$ is constant with time in every gas cell during the radiative transfer calculation, however in an attempt to produce slightly more accurate transmission profiles on large scales, we allow the redshift of each cell – and the physical gas densities, assuming $\rho\propto(1+z)^{3}$ – to vary as a function of distance from the quasar according to the cosmological line increment $dl/dz=c/[(1+z)H(z)]$. In practice this has a very modest effect on the mock spectra.

In Figure 1, we show the physical properties and Ly$\alpha$ transmission spectrum of a simulated skewer for a $z=6$ quasar with $M_{1450}=-27$ which has been on for $t_{\rm q}=10^{7.5}$ years, the fiducial quasar age assumed in E17. The excess Ly$\alpha$ transmission due to the proximity effect is clearly evident out to distances of nearly 10 proper Mpc, with excess heat from He II ionization visible out to $\sim 6$ proper Mpc.

2.3 The proximity zone size, $R_{p}$

We measure $R_{p}$ in the RT simulations analogously to the observations. The Ly$\alpha$ transmission spectra are smoothed by a 20 Å top hat filter, and then the first location where the transmission falls below 10% (starting from $R=0$) is recorded as $R_{p}$. The blue curve in the top panel of Figure 1 shows the smoothed spectrum and the location of $R_{p}$. For this particular simulated skewer, $R_{p}\approx 6.0$ proper Mpc, which is typical for quasars of this luminosity at $z\sim 6$ (E17).

In Figure 2 we show how $R_{p}$ evolves as a function of time for three different quasar luminosities. Two distinct regimes of $R_{p}$ growth are evident: an initial rapid rise for $t_{\rm q}\lesssim 10^{5}$ yr, and a modest increase at $10^{7}\lesssim t_{\rm q}\lesssim 10^{8}$ yr. As discussed in E17 and Eilers et al. (2018), the first of these increases is due to the finite response time of hydrogen to the ionizing photons from the quasar. As described in McQuinn (2009), and discussed more in § 3, the IGM responds to changes in ionizing radiation over a characteristic timescale of $t_{\rm eq}\sim\Gamma^{-1}$, which for hydrogen in $z\sim 6$ quasar proximity zones corresponds to $\sim 3\times 10^{4}$ years (E17). The second increase is due to He II photoionization heating, which we discuss further below.

The shaded regions in Figure 2 show the central 68% scatter of $R_{p}$ between different IGM skewers at fixed $t_{\rm q}$. This scatter arises from fluctuations in the density field causing the (smoothed) transmitted flux to drop below 10% at different locations in each spectrum. We show the full distributions of $R_{p}$ for three representative quasar ages in Figure 3, demonstrating that as $R_{p}$ increases, the scatter in $R_{p}$ increases as well. The distributions of $R_{p}$ are largely Gaussian in shape, with weak tails to lower and higher values at early and late times, respectively.

Figure 4 shows the dependence of $R_{p}$ on quasar luminosity in our simulations compared to the measurements from E17. The majority of quasar spectra analyzed by E17 have $R_{p}$ consistent with $t_{\rm q}\ga 10^{6}$ yr, with a few notable exceptions that may indicate much younger ages (see also Eilers et al. 2018). Common theoretical luminosity scalings of $R_{p}$ are shown by the dashed ($R_{p}\propto\dot{N}_{\rm ion}^{1/2}$) and dotted curves ($R_{p}\propto\dot{N}_{\rm ion}^{1/3}$), anchored to the $t_{\rm q}=10^{6}$ yr simulations at $M_{1450}=-27$. The $\dot{N}_{\rm ion}^{1/2}$ scaling arises from the assumption that $R_{p}$ occurs at a fixed quasar ionizing flux (Bolton & Haehnelt, 2007), while the $\dot{N}_{\rm ion}^{1/3}$ scaling comes from balancing the number of emitted ionizing photons with the number of hydrogen atoms in the IGM (equation 1, Cen & Haiman 2000). The scaling in our RT simulations of a highly ionized IGM roughly agree with the $\dot{N}_{\rm ion}^{1/2}$, but subtly disagree due to the effect of IGM density fluctuations and He II reionization heating, which we discuss further below.

Figure 5 shows how $R_{p}$ evolves with redshift with a fixed UV background intensity. The redshift evolution of $R_{p}$ is then driven by cosmological density evolution due to the expansion of the Universe as well as structure formation. We improve upon the redshift evolution of the simulations presented in E17 by showing RT simulations using hydrodynamical simulation outputs at $z=5.5$ and $z=6.5$ rather than simply re-scaling the $z=6.0$ density skewers by $(1+z)^{3}$. The evolution of $R_{p}$ in the simulations is found to roughly follow $R_{p}\propto(1+z)^{-3.2}$, which is fairly shallow over the redshift range covered by E17, consistent with the non-evolution implied by their measurements. In the top panel of Figure 5 we show the “raw" $R_{p}$ measurements from E17 compared to our simulations as a function of $M_{1450}$, with a fixed $t_{\rm q}=10^{6}$ yr. In the bottom panel, we instead show varying $t_{\rm q}$ at fixed $M_{1450}=-27.0$ in the simulations, compared to the “corrected" $R_{p}$ values from E17, $R_{p,{\rm corr}}$, which have been rescaled so as to remove the dependence on quasar luminosity. While Fan et al. (2006) and Carilli et al. (2010) made these corrections assuming $R_{p}\propto\dot{N}_{\rm ion}^{1/3}$, E17 instead used $R_{p}\propto\dot{N}_{\rm ion}^{1/2.35}$, a relationship derived from our RT simulations. However, it is worth noting that this relationship was derived for $t_{\rm q}=10^{7.5}$ yr, but we find that different quasar ages result in subtly different dependences. For example, for $t_{\rm q}=10^{6}$ yr we find $R_{p}\propto\dot{N}_{\rm ion}^{1/2.2}$ is a better fit to the behaviour of the median $R_{p}$.

In general, the measured properties of quasar proximity zones – as traced by $R_{p}$ – appear to be in good agreement with our RT simulations. This agreement is not trivial, as it does not represent any kind of arbitrary tuning of, e.g., the IGM density field or of the quasar SEDs to produce more or fewer ionizing photons.

2.4 Impact of He II reionization heating

Prior to the quasar turning on, the helium in the simulation is almost entirely in the He II ionization state. The reionization of He II by the quasar heats the gas by $\sim 10^{4}$ K within a radius given by the He II analogy of equation (1),

[TABLE]

where $\dot{N}_{\rm ion}^{\rm He\,II}=4^{-1.7}\dot{N}_{\rm ion}$ is the number of He II-ionizing photons emitted by the quasar for our choice of quasar SED, $n_{\rm He}$ is the helium number density, and $x_{\rm He\,II}$ is the He II fraction (which should be close to unity at $z>5$). This excess heat, known as the thermal proximity effect (Meiksin et al., 2010; Bolton et al., 2010, 2012; Khrykin et al., 2017), reduces the H I fraction of the IGM somewhat ($x_{\rm HI}\propto\Gamma_{\rm HI}^{-1}T^{-0.7}$), and thus leads to extra transmission at $R<R_{\rm ion}^{\rm He}$. The increase in Ly$\alpha$ transmission, and corresponding increase in $R_{p}$ at late times is then due to $R_{\rm ion}^{\rm He}$ reaching the unheated $R_{p}$.

We show the effect of this helium heating via stacked Ly$\alpha$ transmission spectra for a progression of $t_{\rm q}$ in Figure 6. The evolution from $t_{\rm q}=10^{4}$ (purple) to $10^{6}$ years (red) is predominantly due to the non-equilibrium ionization of hydrogen at $10^{4}$ years discussed above – after $10^{6}$ years, the hydrogen is in photoionization equilibrium, and so the profile should remain static at later times. However, we can see that after $10^{8}$ years (blue), there is substantially more transmission at very large radii. The impact of the helium heating can be seen by the extended transmission excess in the $10^{8}$ year stack relative to the $10^{6}$ year stack, where the dashed vertical lines show that the typical He 3 region (computed via equation 12) increases in size from $\sim 2$ to $\sim 9$ proper Mpc (the same effect has been seen in simulations of the He II proximity effect, Khrykin et al. 2016).

This “helium bump" does not manifest as a sharp transition at $R_{\rm ion}^{\rm He}$ in the RT simulation stacks for a few different reasons. First, the inhomogeneous density field along the simulation skewers introduces scatter in the location of $R_{\rm ion}^{\rm He}$ from sightline to sightline. Second, the He 3 ionization front itself is somewhat broad, on the scale of $\sim 1$ proper Mpc (e.g. Khrykin et al. 2016). Finally, there is considerable heating of the IGM beyond the He 3 front from photons at X-ray energies ($\gtrsim 500$ eV) that leads to an extended tail of heated gas at large radii (e.g. McQuinn 2012; Davies et al. 2016).

3 A New Semi-Analytic Model of Proximity Zones

It is worthwhile understanding where the $R_{p}$ scalings in the RT simulations come from to build intuition into how to interpret the observed proximity zones. In this section we build upon the analytic model of $R_{p}$ in the ionized IGM introduced by Bolton & Haehnelt (2007) in two important ways. First, instead of assuming that the transmission through the IGM can be described by the GP optical depth, we calibrate a relation between the hydrogen photoionization rate and the effective optical depth $\tau_{\mathrm{eff}}$ using our hydrodynamical simulation described in § 2.1. Second, given the exceptionally small proximity zones in E17 and the time dependencies seen in the RT simulations (§ 2.3), we introduce non-equilibrium effects to explore how they could manifest in proximity zone measurements.

3.1 Equilibrium model

A natural location one might predict for $R_{p}$ is the point where the mean transmission of the Ly$\alpha$ forest in the enhanced ionizing radiation field of the quasar is equal to 10% (Maselli et al., 2009). It is not immediately obvious that this location is where the first crossing below 10% transmission should be in individual spectra, but as we will show later, it appears to be a reasonable approximation. Similar to Maselli et al. (2009)111Contrary to Maselli et al. (2009), however, because we assume a post-reionization IGM we do not consider this $R_{p}$ to be a “maximum” value, but rather a characteristic one., we define the mean transmission profile in terms of the effective optical depth $\tau_{\mathrm{eff}}\equiv-\ln{\langle F\rangle}$ as determined from our hydrodynamical simulation instead of $\tau_{\rm GP}$, because $\tau_{\mathrm{eff}}$ explicitly accounts for the effect of density and velocity structure of gas in the IGM on the Ly$\alpha$ transmission (on average).

In Figure 7 we show the relationship between $\tau_{\mathrm{eff}}$ and the photoionization rate $\Gamma_{\rm HI}$ measured by computing the mean transmission through random skewers in our hydrodynamical simulation with different $\Gamma_{\rm HI}$. The dotted line shows the effective optical depth corresponding to the 10% transmission threshold, $\tau_{\rm lim}=-\ln(0.1)=2.3$. We see that $\tau_{\mathrm{eff}}=\tau_{\rm lim}$ at $\Gamma_{\rm HI}\approx 1.3\times 10^{-12}$ s*-1*. The behaviour of $\tau_{\mathrm{eff}}$ as a function of $\Gamma_{\rm HI}$ in the regime of interest ($\tau_{\mathrm{eff}}\sim 2$–$6$) is well-described by a power law, $\tau_{\mathrm{eff}}\propto\Gamma_{\rm HI}^{-\alpha}$ with $\alpha\approx 0.55$222The exact fit is $\tau_{\mathrm{eff}}=5.678\,(\Gamma_{\rm HI}/[2.5\times 10^{-13}\,{\rm s}^{-1}])^{0.5486}$.. Note that this quantitative relationship between $\tau_{\mathrm{eff}}$ and $\Gamma_{\rm HI}$ is sensitive to the thermal state of the IGM in the simulation (e.g. Becker & Bolton 2013), so the absolute scaling of Ly$\alpha$ transmission inside the proximity zone – both in this analytic description and in the radiative transfer simulations – will vary depending on how the heat injection by reionization is modeled.

Using this simple power-law relationship, and neglecting the small-scale density enhancement close to the quasar due to its massive host halo, we can estimate $\tau_{\mathrm{eff}}$ as a function of line-of-sight distance $r$ from a quasar,

[TABLE]

where $\Gamma_{\rm q}$ is the ionization rate from the quasar radiation alone, and $\tau_{\rm b}(\Gamma_{\rm b},z)$ and $\Gamma_{\rm b}$ are the effective optical depth and ionization rate of the IGM in the absence of the quasar (i.e. the background opacity and ionization rate). Assuming that the mean free path of ionizing photons emitted by the quasar is large compared to our region of interest, the ionization rate from the quasar can be written as $\Gamma_{\rm q}(r)=\Gamma_{\rm b}(r/r_{\rm b})^{-2}$ (see the lower panel of Figure 1), where $r_{\rm b}$ is the distance at which the ionization rate from the quasar is equal to the average value in the ambient IGM. For our assumed quasar SED, $r_{\rm b}$ is given by

[TABLE]

where $\dot{N}_{\rm ion}=1.73\times 10^{57}$ s*-1* corresponds to a quasar with $M_{1450}=-27$ and $\Gamma_{\rm b}=2.5\times 10^{-13}$ is consistent with constraints from the Ly$\alpha$ forest (Davies et al., 2018a). The effective optical depth profile along the line of sight to the quasar can then be written as333Calverley et al. (2011) derived a similar expression with $\alpha=1$ to describe the profile of the GP optical depth.

[TABLE]

We show example $\tau_{\mathrm{eff}}$ profiles for different quasar luminosities in Figure 8. The proximity zone size $R_{p}$ is then found by solving for $r$ at the limiting optical depth $\tau_{\rm lim}$,

[TABLE]

From equation (16), we can derive the expected scalings of $R_{p}$ with various parameters. The scaling of $R_{p}$ with quasar ionizing photon output $\dot{N}_{\rm ion}$ is built into the definition of $r_{\rm b}$; at fixed $\Gamma_{\rm b}$, we have $r_{\rm b}\propto\dot{N}_{\rm ion}^{1/2}$, so we recover $R_{p}\propto\dot{N}_{\rm ion}^{1/2}$ as in BH07. In Figure 9 we compare the predicted $R_{p}$ as a function of $M_{1450}$ (black) to the RT simulations at $10^{6}$ (red) and $10^{8}$ (blue) years. As mentioned previously, $\dot{N}_{\rm ion}^{1/2}$ is a slightly steeper dependence than what we find in the RT simulations, due to the unaccounted-for effects of density fluctuations and He II reionization heating (§ 2.4).

Given the way we have written equation (16), the scaling with $\Gamma_{\rm b}$ is somewhat more complicated, as we must also account for the change in the background optical depth $\tau_{\rm b}\propto\Gamma_{\rm b}^{-\alpha}$. In the limit where the background opacity is much larger than the $R_{p}$ threshold, $(\tau_{\rm b}/\tau_{\rm lim})^{1/\alpha}\gg 1$, we have

[TABLE]

where in the last step we use the relation $r_{\rm b}\propto\Gamma_{\rm b}^{-1/2}$ at fixed $\dot{N}_{\rm ion}$ as mentioned above. Thus, if the IGM prior to the quasar turning on is highly opaque, $R_{p}$ should be insensitive to $\Gamma_{\rm b}$ (E17). As $\tau_{\rm b}$ approaches $\tau_{\rm lim}$, however, the predicted $R_{p}$ increases asymptotically, whereas in the RT simulations, IGM density fluctuations can still cause the transmission to drop below 10% even if the mean transmitted flux is greater than 10%. In Figure 10, we compare the predicted $R_{p}$ vs $\Gamma_{\rm b}$ relationship from equation (16) to the RT simulations. The dependence of the analytic $R_{p}$ on $\Gamma_{\rm b}$ steepens greatly as $\tau_{\rm b}$ approaches $\tau_{\rm lim}$, while the RT simulations show a much weaker trend. We also show a dotted curve in Figure 10 representing the interpretation of $R_{p}$ evolution in the context of $R_{\rm ion}$ by Fan et al. (2006) and Carilli et al. (2010) – it is clear that this does not describe the dependence particularly well.

The evolution of $R_{p}$ with redshift depends on the cosmological and dynamical evolution of the IGM density field. We find that, between the different redshift outputs of our hydrodynamical simulation ($z=5.5,6.0,6.5$), $\tau_{\rm b}$ at fixed $\Gamma_{\rm b}$ evolves roughly as $(1+z)^{3.5}$. In the absence of $\Gamma_{\rm b}$ evolution, and assuming that $\tau_{\rm b}\gg\tau_{\rm lim}$, we then have $R_{p}\propto(1+z)^{-3.5/(2\alpha)}\sim(1+z)^{-3.2}$, practically identical to evolution measured from the RT simulations. This evolution is somewhat faster than the $\propto(1+z)^{2.25}$ relationship from Bolton & Haehnelt (2007) (equation 1) due to their neglect of the evolution of cosmological structure, and certainly faster than the $\propto(1+z)$ relationship predicted for a fully neutral IGM (derived from the $n_{\rm H}^{1/3}$ dependence in equation 1).

Let us now consider a fiducial $z=6$ quasar with $M_{1450}=-27$ in an IGM with $\Gamma_{\rm b}=2.5\times 10^{-13}$ (Davies et al., 2018a), corresponding to $\tau_{\rm b}=5.7$ and $r_{\rm b}=11.3$ proper Mpc. From equation (16) we then expect $R_{p}=5.5$ proper Mpc, similar to the observations and to the RT simulations with $t_{\rm q}=10^{8}$ years (see also Figure 9). However, the equilibrium model in equation (16) does not include a prescription for He II reionization heating, and thus it should instead be compared to the RT simulations with $t_{\rm q}\sim 10^{6}$ years, which instead have $R_{p}=4.5$ proper Mpc. Thus, the analytic approach moderately overestimates the normalization of $R_{p}$. This overprediction comes about due to the first-crossing definition of $R_{p}$ – there is considerable path length at shorter distances where IGM fluctuations can drop the transmitted flux below 10%. In the Appendix we further explore the effect of IGM fluctuations on proximity zone sizes, and show that including a prescription IGM fluctuations does indeed make the dependence on $\dot{N}_{\rm ion}$ shallower.

3.2 Non-equilibrium model

The connection between $\tau_{\mathrm{eff}}$ and $\Gamma_{\rm HI}$ assumed in equation (16) relies on the assumption of ionization equilibrium – more accurately, the connection should be described in terms of $\tau_{\mathrm{eff}}$ and $x_{\rm HI}\propto\Gamma_{\rm HI}^{-1}$. When the quasar turns on, $x_{\rm HI}$ in the vicinity of the quasar evolves over a characteristic timescale $t_{\rm eq}=\Gamma_{\rm HI}^{-1}$ (ignoring recombinations). This behaviour is captured by the solution to equation (4) assuming a constant ionization rate (e.g. McQuinn 2009; Khrykin et al. 2016),

[TABLE]

where $x_{\rm HI,eq}$ and $x_{\rm HI,0}$ are the equilibrium ($t=\infty$) and initial ($t=0$) neutral fractions, respectively. We can then re-write the above equation in terms of relative neutral fraction,

[TABLE]

We know that $\tau_{\mathrm{eff}}\propto\Gamma_{\rm HI}^{-\alpha}$ and $x_{\rm HI,eq}\propto\Gamma_{\rm HI}^{-1}$, so assuming that the IGM is initially in equilibrium with the UV background, we have

[TABLE]

Finally, substituting $t_{\rm eq}=1/\Gamma_{\rm HI}$ and $\Gamma_{\rm HI}=\Gamma_{\rm q}(r)+\Gamma_{\rm b}$, and assuming $\Gamma_{\rm HI,0}=\Gamma_{\rm b}$, we have

[TABLE]

which asymptotes to equation (15) when $t\gg t_{\rm eq}$ as expected.

In Figure 11, we show the $\tau_{\mathrm{eff}}(r,t)$ (top) and $F=e^{-\tau_{\mathrm{eff}}}$ (bottom) profiles for a series of quasar ages, assuming $M_{1450}=-27$ and $\Gamma_{\rm b}=2.5\times 10^{-13}$ s*-1*. The innermost regions equilibrate first, because $t_{\rm eq}\sim\Gamma_{\rm q}^{-1}\sim r^{2}$, and the entire proximity zone (within the equilibrium $R_{p}$) has reached equilibrium within $\sim 10^{5}$ years. As discussed above in the context of the RT simulations, any link between the small proximity zones seen in E17 and non-equilibrium effects thus requires quasar ages $\la 10^{5}$ years. Figure 11 also shows the RT simulation stacks for $10^{4}$ (purple), $10^{6}$ (red) and $10^{8}$ years (blue). By comparing the RT simulation stacks to the analytic model, we can identify physical features that are unaccounted for in equation (21). At very small radii ($<1$ proper Mpc) the RT simulations show additional absorption due to dense gas associated with the massive quasar host halos (§ 2.1). On larger scales, we see deviations to low optical depth (high transmission) due to He II reionization heating, which propagates outward according to equation (12) (see discussion in § 2.4). Importantly, at distances where we do expect that the analytic model accounts for most of the physics, i.e. $\gtrsim 1$ proper Mpc at $10^{4}$ years and $\gtrsim 2$ proper Mpc at $10^{6}$ years, the agreement with the RT simulations is quite good.

Equation (21) does not have an analytic solution for $R_{p}$, although at radii where $\Gamma_{\rm q}\gg\Gamma_{\rm b}$ it simplifies to

[TABLE]

from which we can then derive the following expression for the very early time evolution of $R_{p}$,

[TABLE]

where $R_{p,{\rm eq}}$ is $R_{p}$ given by equation (16) and $\Gamma_{\rm q}(R_{p,{\rm eq}})\approx 10^{-12}$ s*-1*. Thus at early times when $R_{p}$ is small, i.e. where $\Gamma_{\rm q}$ is very large, we have $R_{p}\propto t^{1/2}$. In Figure 12 we compare this approximate behaviour, shown by the dot-dashed curve, to the RT simulations. While there is general agreement between the two for $t_{\rm q}\lesssim 10^{4}$ yr, it is clear that this approximation fails spectacularly at all later times.

To better predict the evolution of $R_{p}$ on longer timescales, we measure $R_{p}$ from the full non-equilibrium $\tau_{\mathrm{eff}}(r,t)$ profiles (equation 21), as shown in Figure 11, by numerically determining the distance where $\tau_{\mathrm{eff}}(r,t)=\tau_{\rm lim}=2.3$ after smoothing the profiles with a 20 Å ($\approx 1$ proper Mpc) top-hat filter. The smoothing primarily affects the transmission profile at very small radii $R\lesssim 1$ proper Mpc where the proximity zone is “contaminated" by the unabsorbed continuum of the quasar redward of rest-frame Ly$\alpha$. We show the resulting time evolution of $R_{p}$ as the dashed black curve in Figure 12. While the non-equilibrium model qualitatively matches the early-time evolution of $R_{p}$, similar to equation (16) it consistently overpredicts $R_{p}$ compared to the RT simulations. As mentioned in the previous section, and explored further in the Appendix, this overprediction is due to the model not accounting for IGM density fluctuations, which can lead to fluctuations in transmission below 10% when the mean transmission is higher than 10%. We find that a good match between the two at $t_{\rm q}<10^{6}$ yr can be obtained by instead measuring $R_{p}$ in the analytic model at $\tau_{\rm lim}=1.9$ ($F=0.15$). In the rest of this work we will apply this modified $\tau_{\rm lim}$ “fudge factor" to predict the time-dependent behaviour of $R_{p}$ more accurately.

As noted in § 2, the evolution of $R_{p}$ in the RT simulations shows a $\sim 25$% increase at $t_{\rm q}\sim 10^{7-8}$ years due to the reionization of He II by the quasar. We can approximately include this effect in our analytic model by decreasing $\tau_{\mathrm{eff}}$ within the expected radius of the He 3 bubble from equation (12). As a toy model we approximate the amount of helium heating as a doubling of the IGM temperature (corresponding to $\sim 10^{4}$ K of heat input, similar to expectations from RT simulations, e.g. McQuinn et al. 2009; Khrykin et al. 2016), and since $x_{\rm HI}\propto T^{-0.7}$, we have $\tau_{\mathrm{eff}}\propto T^{-0.7\alpha}$, so the resulting effect is a decrease in $\tau_{\mathrm{eff}}$ within $R_{\rm ion}^{\rm He}$ by a factor of $\approx 0.77$. The purple dotted curve in Figure 12 shows the resulting “helium bump," which roughly mimics the behaviour seen in the RT simulations. In detail, the effect of He II reionization heating in the RT simulations is somewhat more extended in time than predicted by this toy model due to the finite width of the ionization front and the extended X-ray heating at large distances discussed in § 2.4.

4 Implications of non-equilibrium proximity zone behaviour

In the previous section, we derived an analytic approximation to the non-equilibrium behaviour of quasar proximity zones. In this section we use this new tool to show that quasar light curves which vary on timescales comparable to $t_{\rm eq}$ can significantly modify the properties of quasar proximity zones.

4.1 Blinking lightbulb model

In § 3 we assumed that the quasar turns on once and has constant luminosity – also known as the “lightbulb" model for quasar light curves. However, quasars may instead undergo short periods of low and high accretion rates (Novak et al., 2011; Schawinski et al., 2010; Schawinski et al., 2015), with flux-limited observations naturally drawing preferentially from the latter. Here we use the analytic model to explore the effect of a simple “blinking lightbulb" toy model where the quasar turns on and off with a fixed duty cycle444The “duty cycle” of quasars is often defined in a cosmological context, i.e. as the fraction of the Hubble time that quasars are active. Here we use the term in a more simple manner to describe the fraction of time that the quasar is luminous after it initially turns on.. For simplicity we neglect heating due to He II reionization, which introduces time dependence on much longer timescales than the equilibration time555Any excursions to low $R_{p}$ will likely cause it to intersect with $R_{\rm ion}^{\rm He}$, even at relatively short quasar ages, so in principle the heating should always be included. However, here we wish to focus on the effects of quasar variability alone..

First let us consider what happens to the proximity zone when a quasar turns off (although, due to the absence of the quasar, it is not observable in this state). For a quasar turning off after an extended period of on-time (i.e. $t_{\rm on}\gg t_{\rm eq}$), the neutral fraction will relax from its equilibrium state as

[TABLE]

where $x_{\rm HI,b}$ is the equilibrium neutral fraction of the IGM in the absence of the quasar, and $x_{\rm HI,eq}^{\rm on}$ is the equilibrium neutral fraction when the quasar was turned on (i.e. the initial state of the gas in this case). The neutral fraction then returns to its original (quasar-free) state within a few $t_{\rm eq}\sim\Gamma_{\rm b}^{-1}$. Subsequent turn-on after an extended off-time $t_{\rm off}\gg t_{\rm eq}$ would then result in evolution identical to equation (18).

The situation is more complicated if the quasar is turning on and off on timescales comparable to $t_{\rm eq}$. In this case, the neutral fraction evolution can be solved as a chain of time-dependent “on" and "off" solutions, with the initial conditions for one set by the final state of the other. In the “on" state, the neutral fraction profile evolves according to equation (19),

[TABLE]

where $t_{\rm on}$ is the time at which the quasar turned on. In the “off" state, we have instead

[TABLE]

where $t_{\rm off}$ is the time at which the quasar turned off. The effective optical depth profiles, and $R_{p}$, can then be computed via the procedure described in § 3.2.

In the top panel of Figure 13 we show examples of $R_{p}$ evolution computed in this fashion for such “blinking lightbulb" quasars with duty cycle $f_{\rm duty}=0.5$ and varying episodic lifetime $t_{\rm ep}$, where the “on" state (when the quasar can be observed) is shown by solid curves and the “off" state (when the quasar is invisible) is shown by dotted curves. For episodes (and off times) longer than the equilibration time, shown by the black curve, we see roughly identical behaviour to the original lightbulb case. The blue curve shows that as the episodes get shorter, the proximity zone never has enough time to fully equilibrate, and so the maximum $R_{p}$ starts to become truncated. If the episodic lifetime is substantially smaller than $t_{\rm eq,off}$ (i.e. the orange curve), the region of enhanced transmission does not have time to relax completely, and $R_{p}$ oscillates around a value corresponding to the equilibrium size for $\dot{N}_{\rm ion}\times f_{\rm duty}$. In the bottom panel of Figure 13 we show quasars with fixed episodic lifetime $t_{\rm ep}=4\times 10^{4}$ yr and varying duty cycle. Large duty cycles, shown by the grey curve, result in a narrow range of possible $R_{p}$ values, in particular not allowing $R_{p}$ to be very small. Smaller duty cycles, shown by the purple and green curves, regularly return to small $R_{p}$ as the IGM is allowed to return closer to its original ionization state.

If observations of quasar proximity zones sample such flickering lightcurves, then the observed distribution of $R_{p}$ will reflect samples from the “on" curves (solid) in Figure 13 plus scatter due to IGM fluctuations. We approximate IGM fluctuations by assuming that $R_{p}$ is drawn from a Gaussian centered on the typical value with scatter that varies as $\sigma(R_{p})=A+BR_{p}^{C}$, where we fit $A=0.155$, $B=0.091$, and $C=1.478$ to the scatter in the RT simulations at $t_{\rm q}<2\times 10^{5}$ years for a $M_{1450}=-27$ quasar. In Figure 14 we show samples from the varying episodic lifetime (left panel) and varying duty cycle (right panel) lightcurves shown in Figure 13. Existing samples of quasar spectra should be able to distinguish between these distributions, and thus constrain both the lifetime and duty cycle of luminous quasar activity.

4.2 General lightcurves

For more complicated quasar lightcurves, equation (4) can be solved more generally,

[TABLE]

where $\Gamma(t)$ is the total H I ionization rate. We can then write the non-equilibrium effective optical depth profile as

[TABLE]

where $\Gamma_{\rm HI}(r,t)$ is the background ionization rate plus the variable contribution from the quasar, and we have assumed ionization equilibrium at $t=0$, i.e. $x_{\rm HI,0}=\alpha_{\mathrm{HII}}n_{e}/\Gamma_{\rm HI}(r,0)$.

In Figure 15 we show a toy quasar lightcurve and its corresponding $R_{p}$ evolution by computing $\tau_{\rm eff}(r,t)$ profiles following equation (28) and numerically determining $R_{p}$. The toy lightcurve, shown in the top panel, consists of short timescale ($\Delta t\sim 100$ yr) and long timescale ($\Delta t\sim 10^{4}$ yr) variations over multiple orders of magnitude in luminosity. The duration of the long timescale variations was chosen to be comparable to the equilibration timescale at the equilibrium $R_{p}$, i.e. $\sim 1/\Gamma_{\rm HI}(R_{p,{\rm eq}})\sim 2.5\times 10^{4}$ yr. The resulting $R_{p}$ evolution, shown by the black curve in the bottom panel, reflects mostly the long timescale evolution, as the ionization equilibration time is $\sim 1/\Gamma_{\rm HI}(R_{p})\sim 2.5\times 10^{4}$ yr. This strong variability decouples the observed luminosity of the quasar from $R_{p}$, as demonstrated by the substantial disagreement with the grey curve in the bottom panel which shows the predicted equilibrium $R_{p}$ from equation (16). The blue curve in the bottom panel of Figure 15 also shows the evolution of a lightbulb lightcurve quasar with $M_{1450}=-27$ for comparison. Despite the variable lightcurve frequently reaching brighter magnitudes, $R_{p}$ fails to reach the equilibrium size for the lightbulb case.

A variable lightcurve, such as the toy model presented in Figure 15, will thus affect the distribution of observed $R_{p}$ at a fixed $M_{1450}$. The black curve in Figure 16 shows the resulting distribution of $R_{p}$ when observed during times when $-27.5<M_{1450}<-26.5$, with approximate IGM density fluctuations added as in § 4.1. The distribution is strongly skewed towards small $R_{p}$, with only rare excursions to $R_{p}>5$ proper Mpc. The blue curve shows the corresponding distribution of $R_{p}$ for the lightbulb model, drawn from the RT simulations (§ 2.2). The orange curve shows the observed distribution of $R_{p}$ from E17 for $z\sim 6$ quasars with $-27.5<M_{1450}<-26.5$ and systemic redshifts coming from either sub-millimeter observations or the Mg 2 emission line (i.e., accurate enough that their uncertainties are subdominant compared to IGM fluctuations). The observed $R_{p}$ distribution is clearly incompatible with the toy lightcurve model, but consistent with the $10^{6}$ year lightbulb model (and as shown in E17, a $10^{7.5}$ year lightbulb model is consistent with the observations as well).

5 Discussion & Conclusion

In this work we presented a suite of RT simulations of quasar proximity zones at $z=6$ to investigate the time-dependent properties of the effective proximity zone size $R_{p}$. The early time ($t_{\rm q}\lesssim 10^{5}$ yr) behaviour of $R_{p}$ is dominated by equilibration of the IGM to the newly elevated ionization rate close to the quasar (E17; Eilers et al. 2018), while the late time ($t_{\rm q}\gtrsim 10^{7}$ yr) behaviour of $R_{p}$ demonstrates a sensitivity to the thermal proximity effect from He II reionization. Leveraging the connection between the ionization rate and effective optical depth of the Ly$\alpha$ forest from a hydrodynamical simulation, we constructed a novel semi-analytic model for the proximity zone to better understand its dependencies on quasar age, luminosity, UV background, and redshift. By basing our model on the effective optical depth of the IGM rather than the GP optical depth, our model provides a quantitative picture for $R_{p}$ that only suffers from the lack of IGM density fluctuations, which we explored in the Appendix. Our model also considers non-equilibrium ionization inside post-reionization proximity zones for the first time666The first time for hydrogen, that is, see Khrykin et al. (2019) for discussion of the much longer equilibration timescale for He II proximity zones at $z\lesssim 4$., in light of the discovery of apparently-young quasars with exceptionally small proximity zones (E17; Eilers et al. 2018). Finally, we employed this model as a tool to efficiently explore how measurements of proximity zone properties can reflect quasar variability on timescales comparable to the equilibration time, or even longer if the “helium bump" can be detected.

Our predictions for the effect of quasar variability on $R_{p}$ in § 4 suggest that the ensemble properties of high-redshift quasar proximity zones can be used to constrain the lifetime and duty cycle of luminous quasar events. Generically, strong quasar variability on timescales comparable to $t_{\rm eq}\sim 1/\Gamma_{\rm HI}(R_{p})\sim 2.5\times 10^{4}$ yr can “decouple" the present-day observed quasar luminosity from the size of the proximity zone as measured by $R_{p}$, as shown by the discrepancy between the black and grey curves in Figure 15. The good agreement between the majority of $R_{p}$ measurements at $z\sim 6$ and the model predictions for $t_{\rm q}\sim 10^{6-8}$ years (E17) should rule out such strong variability – however, the existence of a handful of quasars with very small $R_{p}$ (E17; Eilers et al. 2018) suggests that strong variability on timescales $\lesssim 10^{4}$ years does in fact exist, but must be relatively rare.

We note that the quantitative statements in this work regarding the absolute scale of $R_{p}$ as a function of observed quasar properties depend on our choice of hydrodynamical simulation in two ways. First, as mentioned in § 3 the state of the hydrodynamical simulation at $z\sim 6$ reflects the assumption of a uniform reionization event at $z>10$ with minimal heat injection, which is inconsistent with observations of reionization occurring at much lower redshift (e.g. Planck Collaboration et al. 2018; Davies et al. 2018b) and predictions of inhomogeneous heating of the IGM by $\sim 2\times 10^{4}$ K (e.g. D’Aloisio et al. 2018; Oñorbe et al. 2018). In future work we will investigate how quasar proximity zones in a more realistic hydrodynamical simulation incorporating these features differs from the model presented here. Second, given the small volume of our simulation $(100\,\rm{Mpc}/h)^{3}$ relative to the $\sim 1\,\rm{Gpc}^{-3}$ space density of luminous $z\sim 6$ quasars (e.g. Jiang et al. 2016), we may be underestimating the effect of the local quasar environment on the structure of the proximity zone, although the dependence with halo mass has been shown to be fairly weak (Bolton & Haehnelt, 2007; Keating et al., 2015).

Future samples of high-quality $z\sim 6$ quasar spectra with precise systemic redshift estimates from current observational programs will soon increase the number of $z\sim 6$ proximity zone measurements by a factor of a few. Studying these proximity zones in the context of the non-equilibrium phenomenology developed here will allow for powerful constraints on the lightcurves of accreting supermassive black holes, especially when combined with additional IGM-based probes of quasar activity at lower redshift that are sensitive to much longer timescales (Khrykin et al., 2016; Khrykin et al., 2017, 2019; Schmidt et al., 2017, 2018b, 2018a). While we have focused on $R_{p}$ in this work, in principle the timescale at which the proximity zone is sensitive to variability differs along the line of sight as $t_{\rm eq}\propto r^{2}$. By investigating the entire transmission profile of an ensemble of $z\sim 6$ quasars, it should then be possible to constrain the power spectrum of quasar variability across a wide range of timescales that will not be directly observable for millennia.

Acknowledgements

We thank Z. Lukić for providing access to the Nyx hydrodynamical simulation, and we thank M. McQuinn for helpful discussions during which he derived equation (27) for us. FBD acknowledges support from the Space Telescope Science Institute, which is operated by AURA for NASA, through the grant HST-AR-15014.

Appendix A Analytic model for scatter in proximity zone sizes

The analytic models discussed in § 3 make predictions for $R_{p}$ from the mean transmission profiles in the proximity zone. In reality, however, there is significant scatter due to IGM density fluctuations, and $R_{p}$ is defined as the first time that the transmission falls below 10%. Here we introduce a method to estimate the effect of IGM scatter on proximity zone sizes – however, for simplicity, we assume that the observed spectrum is additionally re-binned to the size of the smoothing scale777This alternative definition already exists in the literature, e.g. Willott et al. (2007), but most works are not clear about whether their smoothing procedure also includes re-binning of the spectrum., 20 Å in the observed frame or $\approx 1$ proper Mpc at $z=6$, to avoid having to model strong correlations from pixel to pixel. We denote this alternative $R_{p}$ definition as $R_{p}^{\rm alt}$.

Assuming no correlations between pixels, the probability for $R_{p}^{\rm alt}=r_{i}$ for any (binned) pixel $i$ can be written as

[TABLE]

where $\tau_{{\rm eff},j}$ is the effective optical depth for pixel $j$ at distance $r_{j}$. That is, the probability that any given distance $r_{i}$ is equal to $R_{p}^{\rm alt}$ is the probability that all smaller radii have $\tau_{\rm eff}<\tau_{\rm lim}$ times the probability that $\tau_{\rm eff}(r_{i})>\tau_{\rm lim}$.

We estimate $P(\tau_{\rm eff})$ by calibrating the strength of IGM density fluctuations on the smoothing scale from our Nyx hydrodynamical simulation, similar to the prescription for $\tau_{\mathrm{eff}}(\Gamma_{\rm HI})$ calibrated in § 3.1. We find that the distribution of $\tau_{\mathrm{eff}}$ on 1 proper Mpc scales, when the mean transmission is 10%, very closely follows a lognormal distribution (e.g. Becker et al. 2007) with $\sigma_{\rm ln\tau}\approx 0.3$. We can then rewrite equation (29) as

[TABLE]

where $\bar{\tau}_{{\rm eff},j}$ is the mean effective optical depth expected in pixel $j$ (i.e. from the mean transmitted flux profile as described in the previous sections). Note that equation (30) is a very rough approximation that does not take into account correlations in the IGM density field or the gradient in $\tau_{\rm eff}$ inside of each 1 proper Mpc pixel. Nevertheless, it should provide a baseline expectation for the fluctuations in $R_{p}^{\rm alt}$.

We compare this simple model for proximity zone scatter to the RT simulations in Figure 17. The median $R_{p}^{\rm alt}$, shown as the black solid curve and shaded region, almost exactly reproduces the dependence of $R_{p}^{\rm alt}$ with $M_{1450}$ from the RT simulations, shown in red. For very luminous quasars, the cumulative probability for $\tau_{\rm eff}>\tau_{\rm lim}$ in regions where $\bar{\tau}_{\rm eff}<\tau_{\rm lim}$ is enough to shift the median $R_{p}$ to lower values. Perhaps more interestingly, the scatter in $R_{p}^{\rm alt}$ is also well-reproduced by equation (30), suggesting that the variation in $R_{p}^{\rm alt}$ between sightlines (and thus likely $R_{p}$ as well) is indeed dominated by the IGM scatter on the smoothing scale.