# Modeling Quasar Proximity Zones in a Realistic Cosmological Environment with
a Self-consistent Light Curve

Yihao Zhou,1 Huanqing Chen,2 Tiziana Di Matteo,1 Yueying Ni,3 Rupert A.C.
Croft, 1 and Simeon Bird4
1McWillams Center for Cosmology, Department of Physics, Carnegie Mellon
University, Pittsburgh, PA 15213, USA
2Canadian Institute for Theoretical Astrophysics, University of Toronto,60 St
George St, Toronto, ON M5R 2M8, Canada
3Harvard-Smithsonian Center for Astrophysics, 60 Garden Street, Cambridge, MA
02138, USA
4Department of Physics & Astronomy, University of California, Riverside, 900
University Ave., Riverside, CA 92521, USA E-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

We study quasar proximity zones in a simulation that includes a self-
consistent quasar formation model and realistic IGM environments. The quasar
host halo is $10^{13}\ M_{\mathrm{\odot}}$ at $z=6$, more massive than typical
halos studied in previous work. Between $6<z<7.5$, the quasar luminosity
varies rapidly, with a mean magnitude of $M_{\mathrm{UV,mean}}=-24.8$ and the
fluctuation reaching up to two orders of magnitude. Using this light curve to
post-process the dense environment around the quasar, we find that the
proximity zone size ($R_{\mathrm{p}}$) ranges between $0.5$-5 pMpc. We show
that the light curve variability causes a similar degree of scatter in
$R_{\mathrm{p}}$ as does the density fluctuation, both of which result in a
standard deviation of $\sim 0.3\ \mathrm{pMpc}$. The $R_{\mathrm{p}}$ traces
the light curve fluctuations closely but with a time delay of $\sim 10^{4}\
\mathrm{yr}$, breaking the correspondence between the $R_{\mathrm{p}}$ and the
contemporaneous $M_{\mathrm{UV}}$. This also indicates that we can only infer
quasar activity within the past $\sim 10^{4}$ years instead of the integrated
lifetime from $R_{\mathrm{p}}$ in the later part of cosmic reionization.
Compared with the variable light curve, a constant light curve underestimates
the $R_{\mathrm{p}}$ by 13% at the dim end ($M_{\mathrm{UV}}\sim-23.5$), and
overestimates the $R_{\mathrm{p}}$ by 30% at the bright end
($M_{\mathrm{UV}}\sim-26$). By calculating the $R_{\mathrm{p}}$ generated by a
number of quasars, we show that variable light curves predict a wider
$R_{\mathrm{p}}$ distribution than lightbulb models, and readily explain the
extremely small $R_{\mathrm{p}}$ values that have been observed.

###### keywords:

quasars: supermassive black holes – intergalactic medium – radiative transfer
– galaxies: high-redshift

††pubyear: 2023††pagerange: Modeling Quasar Proximity Zones in a Realistic
Cosmological Environment with a Self-consistent Light Curve–12

##  1 Introduction

A bright quasar at high redshift usually creates a large region, commonly
referred to as a ‘quasar proximity zone’, where the ionizing radiation
contributed from the quasar significantly exceeds the cosmic ionizing
background. Within a quasar proximity zone, the hydrogen neutral fraction is
considerably lower than typical regions in the Universe. At $z>6$, these are
the only regions where we can observe non-zero Lyman $\alpha\ $transmitted
flux (Bajtlik et al., 1988; Cen & Haiman, 2000; Wyithe et al., 2005; Bolton &
Haehnelt, 2007a, b; Lidz et al., 2007). As a result, quasar proximity zones
are unique windows for probing the distant universe.

One key observational measurement related to a quasar proximity zone is its
size, which is traditionally defined in quasar spectra as the distance from
the systematic Lyman $\alpha\ $line center to the first point where the
transmitted flux drops below 10% of the continuum level after being smoothed
by a 20Å top-hat kernel (Fan et al., 2006; Carilli et al., 2010; Eilers et
al., 2017, 2020; Mazzucchelli et al., 2017; Ishimoto et al., 2020). Fan et al.
(2006) compiled the first large sample of $z\gtrsim 6$ quasar spectra and
measured a proximity zone size $R_{\mathrm{p}}\sim 8\ \mathrm{Mpc}$ at $z=6$
for quasars of magnitude $M_{\mathrm{1450}}=27$. Carilli et al. (2010)
analyzed the proximity zone sizes of 27 quasars with more accurate redshift
measurements. In the last decade, the number of high redshift quasar spectra
with well-measured Lyman $\alpha\ $proximity zone sizes has grown considerably
(Reed et al., 2017; Bañados et al., 2018; Matsuoka et al., 2019; Wang et al.,
2019). For example, Eilers et al. (2017) studied quasar proximity zones in the
redshift range $5.77\leqslant z\leqslant 6.54$ with a homogeneous analysis of
34 medium resolution spectra, and Ishimoto et al. (2020) presented
measurements of the proximity zone size for 11 low-luminosity
($-26.16\leqslant M_{\mathrm{1450}}\leqslant-22.83$) quasars at $z\sim 6$. An
unexpected result from these observations obtained in recent years is that
some quasars display very small proximity zones, such as the $R_{\mathrm{p}}\
$$\sim 0.37$ Mpc measured in Eilers et al. (2021). One possible explanation
for several small $R_{\mathrm{p}}$ seen in $z\gtrsim 7$ quasar spectra is that
a large amount of hydrogen at $z\sim 7$ is still neutral; e.g., Miralda-Escudé
& Rees 1998; Mortlock et al. 2011; Bolton et al. 2011; Bosman & Becker 2015a;
Bañados et al. 2018; Davies et al. 2018; Wang et al. 2020; Yang et al. 2020;
Bosman & Becker 2015b; Greig et al. 2017, 2022). However, the population of
small $R_{\mathrm{p}}$ quasars at $z\sim 6$ remains perplexing.

The sizes of quasar proximity zones could provide valuable insights into
quasar activity and cosmic reionization, an epoch when the IGM transitioned
from a mostly neutral state into a mostly ionized state. Several physically
motivated (semi-)analytic models of proximity zones have been proposed, which
lead to scaling relations between $R_{\mathrm{p}}$ and the intrinsic
properties of quasars (e.g., the luminosity and the lifetime) as well as the
surrounding IGM. In a nearly neutral universe, the proximity zone size is
closely related to the size of the quasar ionized bubble, which is sensitive
to the ionization fraction of the local IGM and the total number of ionizing
photons the quasar emits (Cen & Haiman, 2000; Haiman & Cen, 2001):

$R_{\mathrm{ion}}=\left(\frac{3\dot{N}t_{\mathrm{q}}}{4\pi\,n_{\mathrm{H}}\,x_{\mathrm{H\,\textsc{i}}}}\right)^{1/3},$
(1)

where $\dot{N}$ is the emitted ionizing photon rate, $t_{\mathrm{q}}$ is the
quasar lifetime, $n_{\mathrm{H}}$ is the hydrogen number density and
$x_{\mathrm{H\,\textsc{i}}}$ is the neutral hydrogen fraction. On the other
hand, if the quasar is embedded in an already ionized IGM, Bolton & Haehnelt
(2007a) showed that the proximity zone size quickly reaches a maximum
$R_{\mathrm{p}}^{\mathrm{max}}$ which is independent on the neutral fraction:

$\displaystyle R_{\mathrm{p}}^{\mathrm{max}}$
$\displaystyle=\frac{3.14}{\Delta_{\mathrm{lim}}}\
\left(\frac{\dot{N}}{2\times 10^{57}\ \mathrm{s}^{-1}}\right)^{1/2}\
\left(\frac{T}{2\times 10^{4}\ \mathrm{K}}\right)^{0.35}$ (2)
$\displaystyle\times\left(\frac{\tau_{\mathrm{lim}}}{2.3}\right)^{1/2}\,\left(\frac{\alpha^{-1}\left(\alpha+3\right)}{3}\right)^{-1/2}\left(\frac{1+z}{7}\right)^{-9/4}\
\ \ \mathrm{pMpc}.$

Here $\tau_{\mathrm{lim}}=-\ln(F_{\mathrm{lim}})$ is the limiting Lyman
$\alpha\ $optical depth (typically $F_{\mathrm{lim}}$ is set to be $0.1$) and
$\alpha$ denotes the spectral index of the quasar power-law spectrum.
Recently, Davies et al. (2020) developed a semi-analytic model to trace the
time-dependent evolution of the mean proximity zone size, considering both the
photoionization and photoheating effects.

Although (semi-)analytical models delineate the typical trend in proximity
zone size evolution, it is difficult for them to describe the inhomogeneity of
the IGM, and the complex dependence of $R_{\mathrm{p}}$ on quasar emission
history. As a result, many works combine hydrodynamic simulations and post-
processing with radiative transfer (RT) to investigate the formation and the
properties of proximity zones. For example, Lidz et al. (2007) calculated the
mock Lyman $\alpha\ $flux in quasar proximity zones from simulations and
showed that the absorption spectra are sensitive to the level of small-scale
structure in the IGM. They suggested that the proximity zone measurements are
compatible with the fact that $z\sim 6$ quasars reside in a pre-ionized IGM.
Keating et al. (2015) used simulations that included a wide range of host halo
masses ($10^{10-12.5}\ M_{\odot}$) and discovered that $R_{\mathrm{p}}$
depends weakly on the mass of the host halo, but strongly on direction,
leading to a large $R_{\mathrm{p}}$ scatter across lines of sight. Chen &
Gnedin (2021) found that about $1\%\sim 2\%$ of old quasars exhibit
extraordinarily small proximity zone sizes. Their results are based on the
Cosmic Reionization On Computers (CROC) simulations, which can resolve Lyman
Limit Systems (LLSs) because of its relatively high spatial resolution. These
small proximity zones are predominantly due to the occurrence of a damped
Ly$\alpha$ absorber (DLA) or an LLS along the line of sight. They can be
distinguished from the small proximity zones of young quasars because of the
metal contamination.

Most of the modeling of proximity zones using cosmological simulations post-
processed with RT assumes a ‘lightbulb’ model for quasar light curves. In such
a model, the quasar switches on abruptly and then maintains the same
luminosity. However, realistic quasar light curves are much more complex. They
may instead undergo short periods of flickering, with the fluctuation
magnitude as well as the timescale varying continuously (Novak et al., 2011;
Gabor & Bournaud, 2013; King & Nixon, 2015; Schawinski et al., 2015;
Oppenheimer et al., 2018; Shen, 2021). Extreme variability is observed in some
quasars with magnitude changes of $\sim 1\ \mathrm{mag}$ within $\sim 15\
\mathrm{yr}$ (Rumbaugh et al., 2018). Furthermore, lightbulb models can hardly
explain the small proximity zones observed around luminous quasars at $z\sim
6$. Although Eilers et al. (2017, 2021) offer the suggestion that these small
proximity zones are from young quasars that just turned on for
$t_{\mathrm{q}}=10^{4}$ yr, these estimated lifetime values are significantly
smaller than the timescales required by current supermassive black hole (SMBH)
formation models (Volonteri, 2010; Li, 2012; Madau et al., 2014; Volonteri et
al., 2015; Smith et al., 2017; Weinberger et al., 2018; Regan et al., 2019;
Inayoshi et al., 2020).

To model quasar variability, several works have built toy models for
flickering light curves, most of which have fixed duty cycle
$f_{\mathrm{duty}}$ or episodic lifetime $t_{\mathrm{ep}}$ (Šoltinský et al.,
2023; Satyavolu et al., 2023), and predicted a very different proximity zone
evolution. Davies et al. (2020) applied a light curve consisting of different
time-scale variations and computed $R_{\mathrm{p}}$ based on their time-
dependent semi-analytical model. They found that a variable light curve
yielded a $R_{\mathrm{p}}$ distribution strongly skewed towards small values
compared to a typical lightbulb model.

In this work, we go a step further and study a more realistic quasar light
curve. We utilize data from a state-of-the-art cosmological simulation that
incorporates subgrid models of galaxy and black hole formation, and extract
the light curve from the simulated active galactic nuclei (AGNs). We model the
accretion onto AGNs in the same fashion as the BLUETIDES (Feng et al., 2016)
and MassiveBlack I & II (Di Matteo et al., 2012; Khandai et al., 2015)
simulations. They achieve close agreement with observational measurements of
the relationships between SMBHs and their host galaxies (DeGraf et al., 2015;
Ding et al., 2020, 2022), the observed quasar luminosity functions (QLF), and
the correlation length $r_{0}$ of AGN clustering (Khandai et al., 2015).

The rest of this paper is organized as follows. In Section 2, we briefly
summarize the basic features of the cosmological simulations used in this work
and then illustrate the method used to compute proximity zone size evolution
with RT. In Section 3, we present our results, and discuss our key results in
Section 4. We conclude with Section 5.

##  2 Methodology

###  2.1 Constrained realization simulation

Figure 1: Illustration of the host environment of the quasar prior to the RT
post-processing. The first two rows are the snapshots of the gas density
fields, with the quasar positioned at the center of each panel, at $z=7.5$
(upper) and $z=6.0$ (middle). For regions with high density, the color hue is
set by the hydrogen neutral fraction $x_{\mathrm{H\,\textsc{i}}}$. The sidebar
transitions from red to blue, representing the progression from ionized to
neutral states. The left two panels display cubical boxes with a comoving side
length of $10\ h^{-1}\mathrm{Mpc}$, and the two right panels depict the
zoomed-in central regions of diameter $1\ h^{-1}\mathrm{Mpc}$. The blue,
yellow, green, and red straight lines in the left panels show the directions
of 4 lines of sight used in this work. In the bottom panel, we show the gas
density contrast $\Delta_{\mathrm{g}}$ ($=\rho/\bar{\rho}$) along all 48
directions at $z=7.5$ in grey, with the four examples in the upper panels
highlighted in blue, yellow, green, and red, respectively.  Figure 2: The gas
density contrast ($\Delta_{\mathrm{g}}$) PDF for the lines of sight from the
CR simulation used in this work (blue), and those from the CROC simulation
(orange) used in Chen & Gnedin (2021). Each pixel is $4\ \mathrm{pkpc}$ in
length, and all the data are drawn from $0.1\sim 2\ \mathrm{pMpc}$ regions
from the quasar in the $z=6.5$ snapshot.

To simulate the proximity zones around a long-lived luminous quasar, we use
one of the cosmological simulations with constrained realizations (CR)
described in Ni et al. (2021). Early quasars are extremely rare, and
generally, only large volume cosmological simulations (with box size $\gtrsim
500\ \mathrm{cMpc}$) can capture them. However, an alternative method is to
use the CR technique to model these rare objects. The CR technique (Hoffman &
Ribak, 1991; van de Weygaert & Bertschinger, 1996) imposes various user-
specified constraints on the Gaussian random field of the initial conditions
(ICs). Within a relatively small box, the CR technique allows us to model a
rare high-density peak in the primordial density field that will form a
massive halo hosting a bright quasar at high redshift ($z\sim 7$). Among the
CR simulations carried out in Ni et al. (2022), we select the one that
contains the quasar with the largest instantaneous bolometric luminosity,
since we focus on the brightest quasars in this study. The quasar’s black hole
mass is $M_{\mathrm{BH}}\sim 10^{9}\ M_{\odot}$ at $z=6$.

The CR simulation has a box size of $20\,h^{-1}\ \mathrm{cMpc}$ per side. The
ICs are generated using the most recent implementation of
gaussianCR111https://github.com/yueyingn/gaussianCR, and the simulation is run
to $z=6$ with the massively parallel cosmological smoothed-particle
hydrodynamics (SPH) simulation software MP-Gadget (Feng et al., 2018). The
peak height of our quasar-hosting halo is $5\,\sigma_{0}$ with a peculiar
velocity $\boldsymbol{v}=\boldsymbol{0}\,\text{km}\,\text{s}^{-1}$, where
$\sigma_{0}$ is the density variance after convolution with a Gaussian kernel
of size $R_{\mathrm{G}}=1\ h^{-1}\,\text{Mpc}$. For full details of the
underlying formalism of the CR simulation, we refer the reader to van de
Weygaert & Bertschinger (1996); Ni et al. (2021).

The box contains $2\times 352^{3}$ particles, with a mass resolution
$M_{\mathrm{DM}}=1.2\times 10^{7}M_{\odot}\,h^{-1}$ for the dark matter
particles and $M_{\mathrm{gas}}=2.4\times 10^{7}M_{\odot}\,h^{-1}$ for the gas
particles in the ICs. Star particles have a mass of $M_{*}=6\times
10^{5}M_{\odot}\,h^{-1}$. The gravitational smoothing length is
$1.8\,h^{-1}\text{ckpc}$ for both dark matter and gas particles. The
cosmological model is consistent with the nine-year _Wilkinson Microwave
Anisotropy Probe_ (WMAP) data (Hinshaw et al., 2013) ($\Omega_{0}=0.2814$,
$\Omega_{\Lambda}=0.7186$, $\Omega_{\mathrm{b}}=0.0464$, $\sigma_{8}=0.82$,
$h=0.697$, $n_{s}=0.971$).

###  2.2 Physics: hydrodynamics and sub-grid modelling

Most of the model physics implemented in the CR simulation is the same as in
BLUETIDES (Feng et al., 2016). The gravity is solved with the treePM approach.
The pressure-entropy formulation of SPH is adopted to solve the Euler
equations (Read et al., 2010; Hopkins, 2013). The density estimator uses a
quintic kernel to reduce the noise in the SPH density and gradient estimation
(Liu & Liu, 2010). A range of sub-grid models are applied to simulate galaxy
and black hole formation and associated feedback processes. Radiative cooling
from both primordial gas (Katz et al., 1996) and metals (Vogelsberger et al.,
2014) is considered. Star formation is implemented based on the multiphase
star formation model (Springel & Hernquist, 2003), but incorporating several
effects described in Vogelsberger et al. (2013). The formation of molecular
hydrogen is computed according to the prescription of Krumholz & Gnedin
(2011), and its effect on star formation at low metallicities is considered.
Type II supernova wind feedback is included, using the same model as in the
Illustris simulation (Nelson et al., 2015; Okamoto et al., 2010). Their wind
speeds are assumed to be proportional to the local one-dimensional dark matter
velocity dispersion $\sigma_{\mathrm{DM}}$:
$v_{\mathrm{w}}=\kappa_{\mathrm{w}}\,\sigma_{\mathrm{DM}}$, where
$v_{\mathrm{w}}$ is the wind speed, and the dimensionless parameter
$\kappa_{\mathrm{w}}=3.7$ (Vogelsberger et al., 2013).

Black hole growth and AGN feedback are modeled in the same way as in the
MassiveBlack I & II simulations, based on the black hole sub-grid model
developed in Springel et al. (2005) and Di Matteo et al. (2005). Black holes
are seeded with an initial seed mass of $M_{\mathrm{seed}}=5\times
10^{5}M_{\odot}\,h^{-1}$ in halos with a mass larger than $5\times
10^{10}M_{\odot}\,h^{-1}$. Note that our choice of seed mass is close to that
expected from direct collapse scenarios (e.g., Latif et al. 2013; Schleicher
et al. 2013; Ferrara et al. 2014), but our seeding scheme makes no direct
assumption for the black hole seed formation mechanism.

The gas accretion rate of the black hole is given by the Bondi-Hoyle rate
(Bondi & Hoyle, 1944):

$\dot{M}_{\mathrm{B}}=4\pi\,G^{2}\,M_{\mathrm{BH}}^{2}\,\rho_{\mathrm{BH}}\left(c_{\mathrm{s}}^{2}+v^{2}_{\mathrm{vel}}\right)^{-3/2},$
(3)

where $c_{\mathrm{s}}$ is the local sound speed, $\rho_{\mathrm{BH}}$ is the
gas density around the quasar, and $v_{\mathrm{vel}}$ is the velocity of the
black hole relative to the surrounding gas. Super-Eddington accretion is
allowed with an upper limit of twice the Eddington accretion rate
$\dot{M}_{\mathrm{Edd}}$. Therefore the black hole accretion rate
$\dot{M}_{\mathrm{BH}}$ is determined by
$\dot{M}_{\mathrm{BH}}=\min\left(\dot{M}_{\mathrm{B}},2\dot{M}_{\mathrm{Edd}}\right)$.
With a radiative efficiency $\eta=0.1$ (Shakura & Sunyaev, 1973), the black
hole radiates with a bolometric luminosity $L_{\mathrm{bol}}$ proportional to
the accretion rate: $L_{\mathrm{bol}}=\eta\,\dot{M}_{\mathrm{BH}}c^{2}$. Five
percent of the radiated energy is thermally coupled to the gas residing within
twice the radius of the SPH smoothing kernel of the black hole particle, which
is typically about $1\%\sim 3\%$ of the virial radius of the halo.

The patchy reionization model (Battaglia et al., 2013) is not included in the
CR simulation because of the small box, instead the reionization is assumed to
occur instantaneously. The applied global ionization history is consistent
with that in BLUETIDES, and reionization is almost completed at redshift $z=8$
(Fig. 2 in Feng et al. (2016)). Consequently, the gas in the snapshots we use
herein ($z\leqslant 7.5$) is originally highly ionized.

Figure 3: The mass evolution (blue curve; using right y-axis) and the ionizing
photon number emitted per second by the accretion disk (black curve; left
y-axis) of the brightest quasar in the simulation. The red solid line
represents the photon number rate corresponding to twice the Eddington limit
(upper limit set by our simulation). The yellow curve denotes the averaged
photon number rate $\dot{N}_{\mathrm{mean}}$ computed over a time kernel of
$5\ \mathrm{Myr}$. The upper axis indicates the age of the Universe at a given
redshift. Figure 4: The statistics for the quasar light curves with
$6\leqslant z\leqslant 7.5$. Left panel: the PDF of the ionizing photon rate
contrast $\dot{N}/\dot{N}_{\mathrm{mean}}$, where $\dot{N}_{\mathrm{mean}}$ is
the mean of the ionizing photon number rate $\dot{N}$ emitted by the quasar.
Right panel: the dimensionless power spectrum $kP(k)$ for the time evolution
of $\dot{N}/\dot{N}_{\mathrm{mean}}$.

###  2.3 Line of sight gas densities

The projected gas density fields around the quasar at $z=7.5$ and $z=6.0$ are
shown in the first two rows of Fig. 1, with the quasar located at the center
of each panel. Regions with high density are color-coded by the hydrogen
ionization fraction $x_{\mathrm{H\,\textsc{i}}}$, from red to blue indicating
ionized ($x_{\mathrm{H\,\textsc{i}}}\lesssim 10^{-4}$) to neutral
($x_{\mathrm{H\,\textsc{i}}}\sim 0.1$), as shown by the color bars. The left
two panels display the regions $10\ h^{-1}\mathrm{Mpc}$ in width centered on
the quasar, while the right two panels depict the central zoomed-in regions of
width $1\ h^{-1}\mathrm{Mpc}$.

We use HEALPY222https://github.com/healpy/healpy to cast 48 evenly spaced
lines of sight, starting from the position of the quasar, and employ the SPH
formalism to calculate the gas properties $f(\boldsymbol{x})$ (e.g., density,
velocity, ionization fraction) at position $\boldsymbol{x}$ on the line of
sight (Liu & Liu, 2010):

$\left<f(\boldsymbol{x})\right>=\sum_{j=1}^{N}\frac{m_{j}}{\rho_{j}}f\left(\boldsymbol{x}_{j}\right)\,W\left(\boldsymbol{x}-\boldsymbol{x}_{j},q\right),$
(4)

where $\sum_{j}$ is the sum over all the neighboring gas particles within the
smoothing length $q$, and $W$ is the quintic density kernel. $m_{j}$,
$\rho_{j}$, and $\boldsymbol{x}_{j}$ are the mass, density, and position of
each particle, respectively. Taking advantage of the periodic boundary
conditions of the simulation box, we extend each line of sight to a length of
$40\,h^{-1}\text{Mpc}$. The sightlines are drawn significant off-axes (not
parallel to the $x$, $y$ or $z$ axes), ensuring that none of them travel
through the massive halo again. The spatial resolution is set to be $30$
comoving kpc, equivalent to $\sim 4$ kpc in proper units. We indicate the
directions of 4 lines of sight in the left panels in Fig. 1 (blue, yellow,
green, and red straight lines). In the bottom panel, we show the gas density
contrast $\Delta_{\mathrm{g}}(=\rho/\bar{\rho}$; where $\bar{\rho}$ is the
mean gas density of the Universe) at $z=7.5$ for all the 48 lines of sight.
The four examples shown in the upper panels are also included, represented by
correspondingly colored lines. These $\Delta_{\mathrm{g}}$ profiles enable us
to see that the density peaks in different lines of sight are located at
different radii, and that the density fluctuations along each direction span
about two orders of magnitude.

One of the unique properties of our simulation is the constrained initial
conditions. It is thus useful to compare our sightlines with one of the
simulations studied in Chen & Gnedin (2021), which models a more commonly
occurring type of region without constrained conditions. Chen & Gnedin (2021)
studied sightlines drawn from the B40E CROC simulation, where the box size is
$40\,h^{-1}\ \mathrm{cMpc}$ on each side. The CROC simulation is run with the
Adaptive Refinement Tree code (Kravtsov, 1999; Kravtsov et al., 2002; Rudd et
al., 2008), with a base resolution of $39\,h^{-1}\mathrm{ckpc}$ and a peak
resolution of 100 pc. All the sightlines are drawn from halos with dark matter
mass larger than $1.5\times 10^{11}\ M_{\odot}$. Another major difference is
that the CROC simulation models reionization by star particles self-
consistently, resulting in a volume-weighted neutral fraction
$\left<x_{\mathrm{H\,\textsc{i}}}\right>_{\rm v}=0.13$ at $z=7.33$ and
$\left<x_{\mathrm{H\,\textsc{i}}}\right>_{\rm v}<6.7\,\times 10^{-4}$ after
$z=6.7$.

In Figure 2, we compare the density contrast of pixels, each 4 pkpc in length,
in the same range from 0.1-2 pMpc between our CR simulation and CROC at
$z=6.5$. We display the density contrast Probability Density Function (PDF)
for the lines of sight drawn from the CR simulation (blue curve) and that for
the CROC (orange curve) in Fig. 2. It can be seen that our lines of sight have
many more pixels with high density: the fraction of points with
$\Delta_{\mathrm{g}}>100$ is one order of magnitude larger than that in CROC.
This difference is probably because the quasar host halo in the CR simulation,
which reaches a mass of $\sim 10^{13}\ M_{\odot}$ at $z=6$, is much more
massive than the halos selected in Chen & Gnedin (2021), whose halo masses are
$\gtrsim 1.5\times 10^{11}M_{\odot}$. In fact, the halo chosen in this work is
more massive than almost all the halos in previous simulations of proximity
zones, for example, $M_{\mathrm{h}}=2.5\times 10^{12}\ h^{-1}M_{\odot}$ in
Keating et al. (2015), $M_{\mathrm{h}}\gtrsim 10^{11.5}\ M_{\odot}$ in Davies
et al. (2020), and $M_{\mathrm{h}}<10^{12}\ M_{\odot}$ in Satyavolu et al.
(2023). We test how this difference in the surrounding density field affects
the resultant $R_{\mathrm{p}}$ in Appendix A using a constant lightbulb model.
Considering that black holes were added by hand at the center of the halo in
most previous studies, our lines of sight, which are self-consistently drawn
from the host halo of the black hole particle, reflect the environment of the
high-redshift quasar more realistically.

###  2.4 Quasar light curve

To convert the bolometric luminosity $L_{\mathrm{bol}}$ computed in Section
2.2 to the UV ionizing photon number emitted by the quasar per second
($\dot{N}$), we follow the standard procedure (see, also Chen & Gnedin, 2021)
and use a power-law spectral energy distribution (SED) from 1450 Å to 912 Å:
$L_{\nu}\propto\nu^{-\alpha}$ with a spectral index $\alpha=1.5$, normalized
by $L_{\mathrm{bol}}$. This leads to $M_{\mathrm{UV}}$ upon applying the
appropriate bolometric correction (Fontanot et al., 2012):

$M_{\mathrm{UV}}=-2.5\log_{10}\frac{L_{\mathrm{bol}}}{f_{\mathrm{B}}\mu_{\mathrm{B}}}+\Delta_{\mathrm{B,\,UV}}+34.1,$
(5)

where $f_{\mathrm{B}}=10.2$, $\mu_{\mathrm{B}}=6.7\times
10^{14}\,\mathrm{Hz}$, and $\Delta_{\mathrm{B,\,UV}}=-0.48$. We assume the
escape fraction for the quasar is $f_{\mathrm{esc}}=100\%$, consistent with
the large escape fractions inferred from observations (Eilers et al., 2021;
Stevans et al., 2014; Worseck et al., 2014). The total ionizing photon rate is
then given by $\dot{N}=\int_{13.6\,\mathrm{eV}}^{\infty}L_{\nu}/h\nu\,d\nu$,
which translates $M_{\mathrm{UV}}=-26.66$ to $\dot{N}=1\times
10^{57}\,\text{s}^{-1}$, and $M_{\mathrm{UV}}=-27$ to $\dot{N}=1.36\times
10^{57}\,\text{s}^{-1}$.

We compute $\dot{N}$ for the quasar and show the light curve (black curve) in
Fig. 3, compared with $\dot{N}_{\mathrm{Edd}}$ (red curve), which represents
the photon number rate corresponding to $2\dot{M}_{\mathrm{Edd}}$, the upper
limit of the accretion rate. The yellow curve indicates the photon number rate
averaged using a $5\ \mathrm{Myr}$ top-hat kernel ($\dot{N}_{\mathrm{mean}}$).
In this work, we focus on the light curve within $6<z\leqslant 7.5$, a period
during which the average quasar luminosity ($\dot{N}_{\mathrm{mean}}$) reaches
a plateau. The quasar has a mean magnitude of $M_{\mathrm{UV,mean}}=-24.8$ in
this redshift range, which is comparable to currently observed quasars.

The quasar exhibits significant variation in the light curve as opposed to
maintaining a fixed $\dot{N}$. We display the PDF (left panel) and the
dimensionless power spectrum $kP(k)$ (right panel) of
$\dot{N}/\dot{N}_{\mathrm{mean}}$ in Fig. 4, which shows that the quasar
experiences variation in $\dot{N}$ spanning over two orders of magnitude. The
characteristic fluctuation timescale, which is indicated by the peak of the
power spectrum, is around $t=1$ Myr.

###  2.5 Radiative transfer code

To interpret the effect of quasar radiation on the surrounding IGM, we carry
out one-dimensional RT in the manner of Chen & Gnedin (2021). This is
implemented by postprocessing the CR simulation. The RT code solves the time-
dependent ionization and recombination of H i, He i, He ii including the
effect of quasar photoionizing radiation and the cosmic ionizing background.
Temperature evolution is also calculated considering recombination cooling,
collisional ionization cooling, collisional excitation cooling, Bremsstrahlung
cooling, and inverse Compton cooling, as well as the expansion of the
Universe. One improvement in the code compared to previous work (Bolton &
Haehnelt, 2007a; Davies et al., 2020) is the implementation of an adaptive
prediction-correction scheme, which is motivated by the vastly different
temporal behavior of gas at different distances from the quasar.

The quasar light curve from the simulation is passed to the first cell, and
the neutral fraction H i, He i, He ii and temperature are evolved with an
adaptive scheme. At each adaptive time step, the transmitted ionizing spectrum
is passed to the next cell as the incidental spectrum to evolve the next cell.
These operations are executed iteratively for consecutive cells along the line
of sight. For a more complete explanation of the RT code, we direct readers to
Chen & Gnedin (2021).

In order to compute the Lyman $\alpha\ $flux spectra, we convolve the
absorption contributed by all cells with the approximate Voigt profile
proposed by Tepper-García (2006). This profile is calculated using the
hydrogen neutral fraction $x_{\mathrm{H\,\textsc{i}}}$ and gas temperature
output from the RT simulation, in conjunction with the velocity and the
density field from the CR simulation.

##  3 Results

Figure 5:  The spectra computed by the RT post-processing code for one line of
sight drawn from the snapshot at $z=7$. Panels from top to bottom show the gas
density contrast $\Delta_{\mathrm{g}}$ ($=\rho/\bar{\rho}$), the Lyman
$\alpha\ $transmitted flux, the gas temperature, and the fraction of H i/He
i/He ii. The same color scheme is used in both the flux and temperature panels
to illustrate the outputs from the RT simulation at $t_{\mathrm{evol}}=1\
\mathrm{Myr}$ (blue), $5\ \text{Myr}$ (orange), $15\text{Myr}$ (green), $20\
\text{Myr}$ (red), and $25\ \text{Myr}$ (purple). In the second row, a subplot
provides an overview of the light curve along with the $\dot{N}$ for the
selected $t_{\mathrm{evol}}$ marked with colored stars, whose values are
enumerated in the legend. The solid curves represent the spectra smoothed by a
20 Å top-hat kernel, with the dashed curves giving the original flux profiles,
and the horizontal dotted line demonstrates the 10% flux threshold. In the
third panel, the black dashed curve depicts the temperature at
$t_{\mathrm{evol}}=0\ \text{Myr}$ (i.e., no quasar radiation). In the bottom
row, H i/He i/He ii fraction (solid blue/orange/green curves) at
$t_{\mathrm{evol}}=5\,\text{Myr}$ are compared with the background ionization
fractions(dashed curves).  Figure 6:  The evolution of the proximity zone size
$R_{\mathrm{p}}$ with a time resolution of $\Delta t_{\mathrm{evol}}=0.25\
\mathrm{Myr}$. The top 6 panels depict the $R_{\mathrm{p}}$ (blue curve; left
y-axis) for 6 lines of sight, compared with the quasar light curve $\dot{N}$
(red curves, right y-axis). The vertical grey dash lines denote the redshift
of available snapshots $z_{\mathrm{snap,i}}$, and the upper axis indicates the
age of the Universe at a given redshift. The bottom panel shows the average
proximity zone $\left<R_{\mathrm{p}}\right>$ (blue curve) and the 16-84th
percentile scatter $\sigma_{\mathrm{R}_{\mathrm{p}}}$ (blue shaded area)
across 48 directions for the specific redshift. The standard deviation of
$\left<R_{\mathrm{p}}\right>$, which represents the scatter caused by the
light curve variation, is $0.33\ \mathrm{Mpc}$. And the mean of
$\sigma_{\mathrm{R}_{\mathrm{p}}}$ for the entire redshift range, indicating
the influence of density fluctuations, is $0.28\ \mathrm{Mpc}$.

###  3.1 Evolution of $R_{\mathrm{p}}$

Figure 7: The dependence of proximity zone size $R_{\mathrm{p}}$ on quasar
instantaneous magnitude $M_{\mathrm{UV}}$. Solid curves in the lower left
panel show median $R_{\mathrm{p}}$ as a function of $M_{\mathrm{UV}}$ for
$6.0\leqslant z\leqslant 7.5$ (blue), which is the entire redshift range for
this study, and $6.0\leqslant z\leqslant 6.5$ (red). The blue dotted curve
indicates the best power-law fit $R_{\mathrm{p}}-M_{\mathrm{UV}}$ scaling
relation: $R_{\mathrm{p}}\propto L^{\mathrm{0.32}}$. Observational
measurements from Eilers et al. (2017), Eilers et al. (2020) and Ishimoto et
al. (2020) are displayed by the black and green dots, which are made for
quasars at $6.0\leqslant z\lesssim 6.5$. The minimum and maximum simulated
values for $6.0\leqslant z\lesssim 6.5$ are depicted with the red dashed
curves for comparison. The upper and the right panels present the marginal
distribution of $M_{\mathrm{UV}}$ and $R_{\mathrm{p}}$, respectively. The blue
dash-dot lines represent the mean values of the $M_{\mathrm{UV}}$ and the
$R_{\mathrm{p}}$ for $6.0\leqslant z\leqslant 7.5$:
$M_{\mathrm{UV,mean}}=-24.8$, $R_{\mathrm{p,mean}}=1.37$ pMpc.

We present the post-processed spectra obtained from the RT code for a single
line of sight at $z=7$ in Fig. 5. The gas density contrast
$\Delta_{\mathrm{g}}$ is shown in the first row and the temperature prior to
quasar activation (i.e., $t_{\mathrm{evol}}=0\,\text{Myr}$) is shown in the
third row with the black dashed line. The transmitted flux and IGM temperature
at $t_{\mathrm{evol}}=1\,\text{Myr}$ (blue), $5\,\text{Myr}$ (orange),
$15\,\text{Myr}$ (green), $20\,\text{Myr}$ (red), and $25\,\text{Myr}$
(purple) are depicted in the second and third rows, respectively, with the
corresponding $\dot{N}$ values enumerated in the legend of the second panel. A
subplot in the second row provides an overview of the light curve variation
along with the selected $\dot{N}$ marked by stars. The traditional definition
of the proximity zone edge is the point where the Lyman $\alpha\ $flux first
drops below 10% after being smoothed by a 20Å top-hat window. We show the
smoothed flux (solid curves) as well as the flux threshold 10% (horizontal
grey dashed line) in the second panel, and the $R_{\mathrm{p}}$ is indicated
by the intersection of the horizontal line and the solid curves. It is
noteworthy that the higher the stars marked in the light curve, i.e., larger
instantaneous $\dot{N}$, the higher the corresponding flux is. This implies
that the levels of the Lyman $\alpha\ $flux, and consequently the proximity
zone size $R_{\mathrm{p}}$, are primarily determined by $\dot{N}$ while
showing no correlation with $t_{\mathrm{evol}}$.

The profiles of H i/He i/He ii fractions (solid blue/orange/green curves) at
$t_{\mathrm{evol}}=5\,\text{Myr}$ are displayed in the bottom row, compared
with the background profile (dashed lines). In the original CR simulation,
helium is predominantly found in the He ii state. As the quasar’s radiation
ionizes the He ii, a ‘He ii proximity zone’ emerges, and the energy injected
from the He ii reionization heats the surrounding IGM, which is known as the
‘thermal proximity effect’ (Bolton et al., 2010, 2012; Meiksin et al., 2010).
This phenomenon also enhances the size of the Lyman $\alpha\ $proximity zone
since the H i fraction is partially temperature-dependent (Davies et al.,
2020). As observed in the third panel of Fig. 5, the heated region, unlike
$R_{\mathrm{p}}$, continues to expand monotonically with the increasing
$t_{\mathrm{evol}}$, suggesting its potential application in the estimation of
quasar lifetimes (see Section 4.1).

In order to fully exploit the 21 snapshots within the redshift range:
$6<z_{\mathrm{snap,i}}\leqslant 7.5\ (1\leqslant i\leqslant 21;\
z_{\mathrm{snap,i+1}}<z_{\mathrm{snap,i}})$, we conduct an RT calculation on
the $i\mathrm{th}$ snapshot for a period of

$t_{\mathrm{evol,i}}=t_{\mathrm{age}}(z_{\mathrm{snap,i+1}})-t_{\mathrm{age}}(z_{\mathrm{snap,i}}),$
(6)

where $t_{\mathrm{age}}(z)$ is the age of the universe at redshift $z$,
$z_{\mathrm{snap,i+1}}$ is the redshift of the subsequent snapshot and
$z_{\mathrm{snap,22}}=6$. When we concatenate the next snapshot at
$z_{\mathrm{snap,i}}\ (i>1)$, we draw the density and velocity along the lines
of sight from the new snapshot, which has not been post-processed, while using
the temperature and the ionization fraction of H i/He i/He ii output by the RT
code from the previous snapshot. We then concatenate these evolution segments,
each reflecting the $R_{\mathrm{p}}$ during
$t_{\mathrm{age}}(z_{\mathrm{snap,i}})\sim
t_{\mathrm{age}}(z_{\mathrm{snap,i+1}})$ interval. Such a stitching operation
is implemented across all the lines of sight, each maintaining a fixed
direction throughout the entire redshift range. We estimate the uncertainty in
the resultant $R_{\mathrm{p}}$ caused by the breaking of the continuous
evolution of the IGM caused by this procedure in Appendix B. We find it tiny
compared to the $R_{\mathrm{p}}$ scatter caused by the underlying density
fluctuations.

We capture the variations in proximity zone sizes across 48 directions with a
variable light curve for $\sim 240\ \mathrm{Myr}$ after the quasar’s
activation. Fig. 6 presents the evolution of $R_{\mathrm{p}}$ with a time
resolution of $\Delta t_{\mathrm{evol}}=0.25\ \mathrm{Myr}$ for 6 sample
sightlines (blue curves in the upper six panels) as well as an average value
for 48 directions (blue curve in the bottom row). As a comparison, we also
plot the quasar light curves $\dot{N}$ using the red curves. The blue shaded
area in the bottom panel depicts the 16-84th percentile scatter among
different directions. The vertical grey dash lines indicate the redshifts of
available snapshots $z_{\mathrm{snap,i}}$. Since the reionization is virtually
completed at $z=7.5$, the scatter between different lines of sight is entirely
attributable to the underlying density field at a specific point in time (Lidz
et al., 2006). On the other hand, the fluctuations in the quasar light curve
contribute significantly to the variability in $R_{\mathrm{p}}$ for an
individual direction within short time frames. The extent of this variability
in $R_{\mathrm{p}}$ depends on the specific density field. This is evidenced
by a comparison of the proximity zone evolution between sightlines 6 and 7
(the second and third rows in Fig. 6): the same light curve fluctuations
result in changes of $\Delta R_{\mathrm{p}}\sim 1.7\ \mathrm{Mpc}$ for
sightline 7 within $0.25\ \mathrm{Myr}$, while for sightline 6, the changes
are nearly all less than $\Delta R_{\mathrm{p}}\sim 0.5\ \mathrm{Mpc}$.

We can quantify the influence of the underlying density fluctuations by
computing the standard deviation of proximity zone sizes
($\sigma_{R_{\mathrm{p}}}$) across sightlines for the same redshift. We find
the mean of this standard deviation ($\left<\sigma_{R_{\mathrm{p}}}\right>$)
is $0.28\ \mathrm{Mpc}$ across the entire redshift range. On the other hand,
to show the influence of light curve variability, we compute the mean of
$R_{\mathrm{p}}\ $$(\left<\mathrm{R}_{\mathrm{p}}\right>)$ across all
sightlines at a specific redshift, and then calculate the standard deviation
of these mean $R_{\mathrm{p}}\ $values for different redshifts
($\sigma_{\left<\mathrm{R}_{\mathrm{p}}\right>}$), which is $0.33\
\mathrm{Mpc}$. This illustrates that compared to the density fluctuations, the
variations in the light curve have a similar influence, or slightly larger, on
the scatter of $R_{\mathrm{p}}$ values.

Figure 8: The averaged evolution of the proximity zone size $R_{\mathrm{p}}$
(blue curve; left y-axis) compared with the quasar light curve $\dot{N}$ (red
curve; right y-axis) for $5\ \mathrm{Myr}$ following
$z_{\mathrm{snap}}=7.160$. The blue shaded region represents the 16-84th
percentile scatter. The subplots zoom in to $80\ \mathrm{kyr}$ time spans
around four peaks in the light curve (red dotted lines in subplots) and their
corresponding $R_{\mathrm{p}}$ peaks (blue dotted lines in subplots), which
are denoted by the black shaded rectangles. In the subplots, we label the time
lag between the light curve peaks and the $R_{\mathrm{p}}$ peaks, which are
$\sim 10^{4}\ \mathrm{yr}$.

###  3.2 $R_{\mathrm{p}}\,-M_{\mathrm{UV}}$ scaling relation

In Fig. 7, the lower left panel displays the proximity zone size
$R_{\mathrm{p}}$ as a function of the quasar’s instantaneous magnitude
$M_{\mathrm{UV}}$. The median $R_{\mathrm{p}}$ across the entire redshift
range ($6.0\leqslant z\leqslant 7.5$) for 48 lines of sight is represented by
the blue solid curve. The best power-law fit $R_{\mathrm{p}}-M_{\mathrm{UV}}$
scaling relation is found to be $\log
R_{\mathrm{p}}\propto-0.13\,M_{\mathrm{uv}}$, i.e., $R_{\mathrm{p}}\propto
L^{0.32}$ (blue dotted curve). The upper and the right panels show the
marginal distribution of $M_{\mathrm{UV}}$ and $R_{\mathrm{p}}$. The blue
dash-dot lines represent the mean values of $M_{\mathrm{UV}}$ and
$R_{\mathrm{p}}$: $M_{\mathrm{UV,mean}}=-24.8$, $R_{\mathrm{p,mean}}=1.37$
pMpc. The $M_{\mathrm{UV}}$ histogram shows that the quasar is relatively
faint ($M_{\mathrm{UV}}\gtrsim-25$) over most of its lifetime, and only
becomes luminous ($M_{\mathrm{UV}}\lesssim-26$) occasionally.

For comparison with observational results, in Fig.7 we plot the observed
proximity zone size for quasars in the range $6\leqslant z\lesssim 6.5$
measured by Eilers et al. (2017, 2020) (black dots) and Ishimoto et al. (2020)
(green dots). We also show the minimum/maximum and median values yielded by
our simulation within the same redshift range using red dashed lines and a
solid line, respectively. Our simulation reproduces the $R_{\mathrm{p}}$ range
for most of the observed quasars with $M_{\mathrm{UV}}>-26$. The inability to
yield $R_{\mathrm{p}}<0.5\ \mathrm{Mpc}$ could be due to the failure to
resolve LLSs or DLAs. On the bright end, our simulation predicts smaller
$R_{\mathrm{p}}$ than some of the observational measurements. This probably
stems from the quasar’s low average luminosity, which we discuss in more
detail in Section 4.2.

Our derived $R_{\mathrm{p}}-M_{\mathrm{UV}}$ scaling relation exhibits a
flatter trend compared to the $R_{\mathrm{p}}\propto L^{0.5}$ predicted by
Bolton & Haehnelt (2007a) for an idealized ionized IGM based on a semi-
analytical model (see equation 2). One possible reason for the disparity is
that the $M_{\mathrm{UV}}$ for our simulated quasar is non-uniformly
distributed across a broad redshift range, as shown in the upper panel of Fig.
7. The $R_{\mathrm{p}}-M_{\mathrm{UV}}$ scaling relation is redshift-
dependent, which can be seen in Fig. 11, where the lower redshift environment
typically produces more extensive proximity zones. According to our
simulation, the optimal fit is $\log R_{\mathrm{p}}\approx\log(3.20\
\mathrm{Mpc})+0.36\left[-0.4(M_{\mathrm{UV}}+27)\right]$ at $6.0\leqslant
z\leqslant 6.5$ and $\log R_{\mathrm{p}}\approx\log(2.41\
\mathrm{Mpc})+0.41\left[-0.4(M_{\mathrm{UV}}+27)\right]$ at $7.3\leqslant
z\leqslant 7.5$. The slope remains shallower than that in Bolton & Haehnelt
(2007a), even when the data is constrained within a narrower redshift band.
Such a weaker dependence of $R_{\mathrm{p}}$ on the instantaneous magnitude
probably originates from the variation in the light curve, which breaks the
correspondence between the $R_{\mathrm{p}}$ and the contemporaneous
$M_{\mathrm{UV}}$ as we discuss in Section 4.1.

##  4 Discussion

###  4.1 Response of $R_{\mathrm{p}}\ $to variable light curve

Several recent studies have focused on constraining quasar lifetimes using
proximity zone sizes under the assumption of a lightbulb model (Morey et al.,
2021; Khrykin et al., 2021). In this section, we briefly discuss the influence
of a variable light curve on quasar lifetime estimation.

We start by exploring the response behavior of $R_{\mathrm{p}}$ to the quasar
light curve. We plot a $5\ \mathrm{Myr}$ duration in the evolution of
$R_{\mathrm{p}}$ starting from $z_{\mathrm{snap}}=7.160$ in Fig. 8, with the
blue curve representing the average value for 48 lines of sight and the shaded
area showing the 16-84th percentile scatter. The subplots zoom in to $80\
\mathrm{kyr}$ time spans around four peaks in the light curve (red dotted
lines in subplots) and their corresponding $R_{\mathrm{p}}$ peaks (blue dotted
lines in subplots). We label the time lags between the light curve peaks and
the $R_{\mathrm{p}}$ peaks in each subplot, which are $12,\ 9,\ 12,\ 14\ $
kyr, respectively. They are comparable to the hydrogen equilibrium time at the
edge of the proximity zone:
$t_{\mathrm{eq}}^{\mathrm{H\,\textsc{i}}}=1/\Gamma^{\mathrm{H\,\textsc{i}}}\sim
10^{4}$ yr, where $\Gamma^{\mathrm{H\,\textsc{i}}}$ is the photoionization
rate of hydrogen (Bolton & Haehnelt, 2007b; Eilers et al., 2018; Davies et
al., 2020). This illustrates that $R_{\mathrm{p}}$ traces the fluctuations in
the light curve closely but with a short delay of $\sim 10^{4}$ yr, which
breaks the correspondence between $R_{\mathrm{p}}$ and the contemporaneous
$M_{\mathrm{UV}}$.

Previous studies have extensively discussed the implications of quasar
proximity zone size for their ‘lifetime,’ using a lightbulb model. This model
describes a quasar turning on suddenly, with its luminosity remaining constant
thereafter. Our simulated light curve features some episodes where the
luminosity increases nearly by a factor of $\times 4$ within $10^{3}$ years
(e.g., at 830 kyr in the first zoom-in panel and at 3945 kyr in the last zoom-
in panel, as shown in Fig. 8). Therefore, these sudden jumps in luminosity can
be viewed as the beginning of an ‘episode’, a term previous studies have used
to describe the quasar’s episodic lifetime (Eilers et al., 2017, 2021).
However, we note that the light curve varies rapidly and seldom behaves like a
lightbulb for more than a few times $10^{3}$ years. By the time $10^{4}$ yr
have passed, denoted as the typical $R_{\mathrm{p}}$ delay time, the quasar’s
luminosity has already changed significantly. Moreover, there are many periods
during which the quasar luminosity evolves slowly (e.g., as seen in the second
and third zoom-in panels of Fig. 8), making the ‘episodic lifetime’ ill-
defined.

In the latter half of reionization, the integrated lifetime of the quasar can
hardly be measured solely based on the size of the proximity zone. The
$R_{\mathrm{p}}$ value only informs us about the quasar’s luminosity within a
span of $10^{4}$ years (as seen in Fig. 8), while the integrated lifetime of
our quasar has been hundreds of million years. To measure this total duration
for which the quasar has been shining, one can use observable associated with
the thermal states of the IGM around the quasar, like the He II proximity
zone, as it has a longer response time $t_{\mathrm{eq}}\sim 10^{6}\
\mathrm{yr}$ (Worseck et al., 2021; Khrykin et al., 2017, 2021; Šoltinský et
al., 2023; Chen et al., 2023).

Figure 9:  Left panel: the evolution of the mean proximity zone sizes
$R_{\mathrm{p}}$ with the variable light curve (blue curves) or a lightbulb
model fixed at $M_{\mathrm{UV,mean}}=-24.72$ (yellow curves). The quasar is
assumed to turn on at $z=6.162$ and last for $30\ \mathrm{Myr}$. The blue
shaded region and the yellow dashed lines indicate the 68% scatter for the two
$R_{\mathrm{p}}$ groups. Middle panel: two-dimensional distribution of
$R_{\mathrm{p}}$ and instantaneous $M_{\mathrm{UV}}$ for
$10<t_{\mathrm{evol}}<30$ Myr. The blue pixels represent the $R_{\mathrm{p}}$
for the variable light curve, whose median/68%/95% confidence regions are
shown by the red solid/dashed/dotted curves. The dots depict the
$R_{\mathrm{p}}$ generated by the lightbulb models fixing the magnitude at
$M_{\mathrm{UV}}=-23.5,\ -24,\ -24.5,\ -25,\ -25.5,\ -26,\ -26.5,\ -27$ (black
dots) and $M_{\mathrm{UV,mean}}=-24.74$ (yellow dot). The thick/thin errorbars
correspond to 68%/95% confidence intervals. The subplot displays the
difference between the median $R_{\mathrm{p}}$ produced by the lightbulb
models and the variable light curve. The horizontal dotted line represents
where the two median values are the same, and the vertical dotted line shows
$M_{\mathrm{UV,mean}}$. Right panel: one-dimensional $R_{\mathrm{p}}$
distribution generated by the part of the variable light curve with
$-26.5<M_{\mathrm{UV}}<-25.5$ (red curve), $-24<M_{\mathrm{UV}}<-23$ (blue
curve), and by a lightbulb model with $M_{\mathrm{UV}}=-26$ (red shaded area),
$M_{\mathrm{UV}}=-23.5$ (blue shaded area) for $10\
\mathrm{Myr}<t_{\mathrm{evol}}<30\ \mathrm{Myr}$.

###  4.2 Influence of light curve variation

In this section, we discuss how the rapid fluctuation in the light curve in
our simulation affects the resultant $R_{\mathrm{p}}$, and compare it with the
proximity zone sizes yielded by a light curve that remains constant over an
extended period of time.

We excerpt a portion of the variable light curve from the simulation at
$z=6.162$ with a time span of $30\ \mathrm{Myr}$, short enough that the
density evolution is negligible. This excerpted light curve has a mean
magnitude of $M_{\mathrm{UV,mean}}=-24.72$ and a mean ionizing photon rate of
$\dot{N}_{\mathrm{mean}}=1.67\times 10^{56}\ \mathrm{s}^{-1}$. In contrast to
this variable light curve, we construct another light curve with a constant
flux of the same $\dot{N}_{\mathrm{mean}}$ (‘lightbulb’ model). We evolve the
same set of $48$ sightlines with these two light curves for $30$ Myr. In the
left panel of Fig. 9, we display the mean $R_{\mathrm{p}}$ at different
$t_{\mathrm{evol}}$ of the variable light curve (blue curve) and the lightbulb
model (solid yellow curve). The blue shaded region and the yellow dashed lines
indicate the 68% scatter for the two $R_{\mathrm{p}}$ groups. For the
lightbulb model, $R_{\mathrm{p}}$ remains nearly unchanged after rapid growth
during the first $\sim 1\ \mathrm{Myr}$, with only a slight decline owing to
the Universe cooling, which aligns with the results of previous studies
(Davies et al., 2020; Eilers et al., 2021). In the middle panel, we show the
$R_{\mathrm{p}}$ distributions as a function of instantaneous magnitude
$M_{\mathrm{UV}}$ for $10<t_{\mathrm{evol}}<30$ Myr, during which the
lightbulb model reaches a stable stage and produces a similar $R_{\mathrm{p}}\
$. The blue pixels represent the $R_{\mathrm{p}}$ for the variable light
curve, whose median, 68%, and 95% scatter regions are shown by the red solid,
dashed, and dotted curves, respectively. The $R_{\mathrm{p}}$ values generated
by the lightbulb models are shown as dots, whose thick and thin errorbars
correspond to the 68% and 95% scatter, respectively. With a similar mean
proximity zone size $\left<R_{\mathrm{p}}\right>\sim 1.6\ \mathrm{pMpc}$
across the entire magnitude range, the variable light curve yields a scatter
($\sigma_{\mathrm{R_{\mathrm{p}}}}=0.46\ \mathrm{pMpc}$) 28% larger than that
of the $\dot{N}_{\mathrm{mean}}$ lightbulb
($\sigma_{\mathrm{R_{\mathrm{p}}}}=0.36\ \mathrm{pMpc}$). The
$\sigma_{\mathrm{R_{\mathrm{p}}}}$ values for the variable light curve
encompass contributions from both light curve variation and underlying density
field fluctuation. On the other hand, the $\sigma_{\mathrm{R_{\mathrm{p}}}}$
for the lightbulb model is almost totally attributed to the density
differences. In addition to the lightbulb model fixed at
$M_{\mathrm{UV,mean}}$ (yellow dot), we also simulate the lightbulb with
different magnitudes (black dots): $M_{\mathrm{UV}}=-23.5,\ -24,\ -24.5,\
-25,\ -25.5,\ -26,\ -26.5,\ -27$ in the middle panel. It can be seen from the
error bars that the influence of the density fluctuation for a lightbulb is
strongly correlated with $M_{\mathrm{UV}}$, and the high luminosity magnifies
the variance between directions, i.e., brighter $M_{\mathrm{UV}}$ leads to an
increase in $\sigma_{R_{\mathrm{p}}}$.

An important feature illustrated by the middle panel of Fig. 9 is that the
median $R_{\mathrm{p}}$ from the variable light curve at a specific magnitude
coincides with the lightbulb model only around $M_{\mathrm{UV,mean}}=-24.72$,
while it tends to yield smaller $R_{\mathrm{p}}$ compared to the lightbulb
when $M_{\mathrm{UV}}<M_{\mathrm{UV,mean}}$, and conversely, larger
$R_{\mathrm{p}}$ when $M_{\mathrm{UV}}>M_{\mathrm{UV,mean}}$. This is more
clearly demonstrated in the subplot, where we show the difference between the
median $R_{\mathrm{p}}$ produced by the lightbulb models and the variable
light curve. The horizontal dotted line represents where the two median values
are the same, and the vertical dotted line shows $M_{\mathrm{UV,mean}}$. More
specifically, in the right panel of Fig. 9 we compare the one-dimensional
$R_{\mathrm{p}}$ distributions within a narrow $M_{\mathrm{UV}}$ bin for $10\
\mathrm{Myr}<t_{\mathrm{evol}}<30\ \mathrm{Myr}$. On the bright end
($M_{\mathrm{UV}}\sim-26$), the lightbulb model predicts a median
$R_{\mathrm{p}}$ 30% larger than that produced by the variable light curve
($2.94$ pMpc versus $2.27$ pMpc). While on the dim end
($M_{\mathrm{UV}}\sim-23.5$), the lightbulb model yields a median
$R_{\mathrm{p}}$ 13% smaller than that of the variable light curve ($1.03$
pMpc versus $1.19$ pMpc). However, these two models give similar scatter in
$R_{\mathrm{p}}$ for this quasar. The standard deviations of $R_{\mathrm{p}}$
for both the variable light curve and the lightbulb model are $\sim 0.6$ pMpc
around $M_{\mathrm{UV}}=-26$, and $\sim 0.3$ pMpc around
$M_{\mathrm{UV}}=-23.5$.

By building a toy model of the fluctuating light curves, Davies et al. (2020)
noticed that the $R_{\mathrm{p}}$ simulated based on variable light curves
skewed towards smaller values as opposed to a lightbulb fixed at a relatively
high luminosity. A similar bias on the bright end emerges in our simulation,
while our computation herein further demonstrates that the discrepancy between
the lightbulb model and the variable light curve is contingent upon the
specific magnitude bin. Such a discrepancy occurs because $R_{\mathrm{p}}$ is
governed by the entire light curve within the most recent $\sim 10^{4}\
\mathrm{yr}$, rather than the contemporaneous instantaneous luminosity. As
depicted by the PDF in Fig. 4, the distribution of $\dot{N}$ is essentially
Gaussian centering around $\dot{N}_{\mathrm{mean}}$, which implies that the
$\dot{N}$ value $10^{4}\ \mathrm{yr}$ preceding a bright or a dim point in the
light curve is probably close to $\dot{N}_{\mathrm{mean}}$, producing a
$R_{\mathrm{p}}$ close to that given by a lightbulb model fixed at
$\dot{N}_{\mathrm{UV,mean}}$. Furthermore, the higher luminosities correspond
to more significant variations since the variation amplitude generally equals
$M_{\mathrm{UV,mean}}-M_{\mathrm{UV}}$. Large variations result in more
remarkable discrepancies (see Appendix C), which explains the more substantial
shift at smaller magnitudes observed in the middle panel of Fig. 9.

Therefore, for an individual quasar whose light curve persistently fluctuates
around a certain value, its proximity zone size displays a shallow evolution
with instantaneous magnitude, and diverts from the lightbulb model in a
$M_{\mathrm{UV}}$-dependent way. Such divergence accounts for the difference
between our predicted $R_{\mathrm{p}}$ and the observational measurements at
$M_{\mathrm{UV}}<-26$ shown in Fig. 7: our simulated light curve generally has
a lower luminosity, which makes the $R_{\mathrm{p}}$ in this magnitude range
smaller; while the observed quasars with $M_{\mathrm{UV}}<-26$ probably have
larger overall luminosity, and so generate large proximity zones.

Figure 10: The $R_{\mathrm{p}}$ distribution with
$-26.5<M_{\mathrm{UV}}<-25.5$ produced by the variable light curves (red solid
curve), and the lightbulb models with $10<t_{\mathrm{evol}}<30$ Myr (red
shaded area), compared with the observed $R_{\mathrm{p}}$ for the quasars
within the same magnitude bin measured in Eilers et al. (2017, 2020); Ishimoto
et al. (2020).

###  4.3 Implications of the observed $R_{\mathrm{p}}$ distribution

Our simulation shows that one single quasar can vary significantly over its
lifetime, with a scatter in luminosity spanning approximately two orders of
magnitude. If we define it according to its mean luminosity
$M_{\mathrm{UV}}^{\star}\ $, our quasar is a relatively faint one
($M_{\mathrm{UV}}\ $=-24.8) during the period $z=7.5\sim 6$. However, it still
has a $16\%$ chance to be caught in a relatively bright phase with
$M_{\mathrm{UV}}\ $<-25.5. In such a bright phase, the distribution of
$R_{\mathrm{p}}\ $is shifted towards the shorter end compared to the case of
constant luminosity (the right panel of Fig. 9). This has profound
implications for interpreting the $R_{\mathrm{p}}\ $distribution at given
observed magnitude $M_{\mathrm{UV}}\ $($P(R_{\mathrm{p}}|M_{\mathrm{UV}})$).

The distribution of $P(R_{\mathrm{p}}|M_{\mathrm{UV}})$ can be formulated as
the following conditional distribution:

$P(R_{\mathrm{p}}|M_{\mathrm{UV}})=\frac{P(M_{\mathrm{UV}}^{\star})P(M_{\mathrm{UV}}|M_{\mathrm{UV}}^{\star})P(R_{\mathrm{p}}|M_{\mathrm{UV}},M_{\mathrm{UV}}^{\star})}{P(M_{\mathrm{UV}}^{\star})P(M_{\mathrm{UV}}|M_{\mathrm{UV}}^{\star})},$

where $P(M_{\mathrm{UV}}^{\star})$ is the probability function of quasars with
a certain mean magnitude $M_{\mathrm{UV}}^{\star}\ $,
$P(M_{\mathrm{UV}}|M_{\mathrm{UV}}^{\star})$ is the probability that a quasar
of $M_{\mathrm{UV}}^{\star}\ $is observed with magnitude $M_{\mathrm{UV}}\ $,
and $P(R_{\mathrm{p}}|M_{\mathrm{UV}},M_{\mathrm{UV}}^{\star})$ is the
probability that the quasar of $M_{\mathrm{UV}}^{\star}\ $at the observed
magnitude $M_{\mathrm{UV}}\ $displays a proximity zone size of
$R_{\mathrm{p}}\ $.

To calculate such a distribution, certain assumptions need to be made. The
first term $P(M_{\mathrm{UV}}^{\star})$ is similar to the quasar luminosity
function, but it is for the mean magnitude $M_{\mathrm{UV}}^{\star}\ $instead
of the observed luminosity. Measuring such a $M_{\mathrm{UV}}^{\star}\
$directly is challenging. For simplicity, here we assume that
$P(M_{\mathrm{UV}}^{\star})$ is equal to the observed quasar luminosity
function (QLF) measured by Matsuoka et al. (2018) for quasars at $z=6$:

$P(M_{\mathrm{UV}}^{\star})\propto\left[10^{-0.156\left(M_{\mathrm{UV}}^{\star}+25.30\right)}+10^{-0.716\left(M_{\mathrm{UV}}^{\star}+25.30\right)}\right]^{-1}.$
(7)

To estimate $P(M_{\mathrm{UV}}|M_{\mathrm{UV}}^{\star})$, we need to know the
PDF of the light curve. Motivated by the luminosity PDF of our simulated
quasar (upper panel of Fig. 7), we assume that the light curve has a Gaussian
distribution centered at $M_{\mathrm{UV}}^{\star}$:

$P(M_{\mathrm{UV}}|M_{\mathrm{UV}}^{\star})=\frac{1}{\sigma\sqrt{2\pi}}\,\exp\left[\frac{1}{2}\left(\frac{M_{\mathrm{UV}}-M_{\mathrm{UV}}^{\star}}{\sigma}\right)^{2}\right].$
(8)

with a fixed scatter $\sigma=0.7$ for all the light curves.

Finally, to model $P(R_{\mathrm{p}}|M_{\mathrm{UV}},M_{\mathrm{UV}}^{\star})$,
we make the following assumptions: (1)
$P(R_{\mathrm{p}}|M_{\mathrm{UV}},M_{\mathrm{UV}}^{\star})$ has the same shape
as the $R_{\mathrm{p}}$ distribution generated by the constant light curve
with magnitude fixed at $M_{\mathrm{UV}}$, which we label as
$P_{\mathrm{lb},M_{\mathrm{UV}}}(R_{\mathrm{p}})$, but is shifted towards a
different mean $R_{\mathrm{p}}$ with the amount $B$. (2) the value of $B$ is
only determined by $M_{\mathrm{UV}}-M_{\mathrm{UV}}^{\star}$ and is
independent of the specific values of $M_{\mathrm{UV}}^{\star}$ and
$M_{\mathrm{UV}}$. Hence, the probability of $R_{\mathrm{p}}$ for a given
variable light curve is

$P(R_{\mathrm{p}}|M_{\mathrm{UV}},M_{\mathrm{UV}}^{\star})=P_{\mathrm{lb},M_{\mathrm{UV}}}(R_{\mathrm{p}}+B).$
(9)

The first assumption is motivated by the right panel of Fig. 9, which shows
that for an individual quasar, the $R_{\mathrm{p}}$ distributions generated by
the variable light curve and the lightbulb have similar shapes: they have
roughly the same scatter but different mean $R_{\mathrm{p}}$. To formulate
$B(M_{\mathrm{UV}}-M_{\mathrm{UV}}^{\star})$, we use our results in Section.
4.2 as a guideline, i.e., we linearly interpolate the difference between the
median $R_{\mathrm{p}}$ of the lightbulb models and the variable light curve
as a function of the magnitude difference
$M_{\mathrm{UV}}-M_{\mathrm{UV}}^{\star}$, which is depicted by the subplot in
the middle panel of Fig. 9.

With the formulas and assumptions stated above, we use Markov Chain Monte
Carlo (MCMC) to generate $N=1000$ quasar samples to create the final
distribution $P(R_{\mathrm{p}}|M_{\mathrm{UV}})$. We present the
$R_{\mathrm{p}}$ distribution with $-26.5<M_{\mathrm{UV}}<-25.5$ calculated by
this model in Fig. 10 (red solid curve). As a comparison, we plot the combined
$R_{\mathrm{p}}$ distribution produced by the lightbulb models for
$10<t_{\mathrm{evol}}<30\ \mathrm{Myr}$ fixed at $M_{\mathrm{UV}}=-25.5,\
-26,\ -26.5$ (red shaded area). We consider three lightbulb models rather than
only the one with $M_{\mathrm{UV}}=-26$ to show the whole $R_{\mathrm{p}}$
range reached by the lightbulb in this $M_{\mathrm{UV}}$ bin, which spans the
minimum $R_{\mathrm{p}}$, reached when $M_{\mathrm{UV}}=-25.5$, to the maximum
$R_{\mathrm{p}}$, for $M_{\mathrm{UV}}=-26.5$. These three lightbulb models
are sampled based on the QLF (equation 7). We also show the $R_{\mathrm{p}}$
values for the observed quasars with $-26.5<M_{\mathrm{UV}}<-25.5$ measured by
Eilers et al. (2017, 2020); Ishimoto et al. (2020) in Fig. 10 (black curve).
It is evident that although a variable light curve produces a similar maximum
value ($R_{\mathrm{p,max}}\sim 5.5$ Mpc) as the lightbulb model, it can also
result in much smaller proximity zones ($R_{\mathrm{p,min}}\sim 0.2\
\mathrm{Mpc}$). Therefore, it readily explains the observations with small
proximity zones. These extremely small $R_{\mathrm{p}}$ values are generated
by the quasar with low averaged luminosity (i.e., large
$M_{\mathrm{UV}}^{\star}$). The considerable difference between the observed
instantaneous $M_{\mathrm{UV}}$, which is $\sim-26$ in this case, and
$M_{\mathrm{UV}}^{\star}$ moves the $R_{\mathrm{p}}$ distribution away from
the distribution without light curve variation, as seen in the right panel of
Fig. 9. Additionally, due to the relative abundance of faint
$M_{\mathrm{UV}}^{\star}$ quasars over the bright ones, the overall
$R_{\mathrm{p}}$ distribution skews towards smaller values.

Since the calculation for the variable light curve and the lightbulb models
are based on the same groups of lines of sight, and because the lightbulb
model accounts for all the scatter introduced by underlying density
fluctuations, the wider scatter shown by the red curve in Fig. 10 can only be
attributed to light curve variability. This underscores the necessity of
considering light curve variability when investigating quasar proximity zones.

Note that Davies et al. (2020) made a similar comparison in their Fig. 16
between the observed $R_{\mathrm{p}}\ $distribution and predicted
$R_{\mathrm{p}}\ $from different quasar light curve models. They concluded
that the lightbulb model with long episodic quasar lifetime ($\geq 1$ Myr)
leads to a $R_{\mathrm{p}}$ distribution consistent with the observed
$R_{\mathrm{p}}\ $distribution, while their toy quasar light curve with
variation does not. On the contrary, our simulated quasar light curve results
in a $R_{\mathrm{p}}\ $distribution that skews only slightly to the smaller
end compared to the lightbulb model (red line versus transparent red shaded
histograms in Fig. 10). We conduct a Kolmogorov–Smirnov (K-S) test and find
the K-S statistic to be $D=0.23$ and $p$-value$=0.54$ when comparing the
observed $R_{\mathrm{p}}\ $distribution and that from the lightbulb model. On
the other hand, we find $D=0.34$ and $p$-value$=0.12$ between the observed
$R_{\mathrm{p}}\ $distribution and the one from our variable light curve.
Therefore, both the lightbulb model and our simulated variable light curve are
compatible with the observed $R_{\mathrm{p}}\ $distribution. We reach this
different conclusion from Davies et al. (2020) because (1) their toy quasar
light curve is constructed to have large variability at very small time scales
($10^{2}$ yr and $10^{4}$ yr), while our simulated quasar light curve has low
power at such small timescales (see Fig. 4); (2) we consider the combined
$R_{\mathrm{p}}$ distribution generated by a group of variable quasars with
different mean magnitude, rather than using an individual quasar; (3) we use
an observational sample from Eilers et al. (2017, 2020); Ishimoto et al.
(2020) while Davies et al. (2020) only have the data from Eilers et al.
(2017).

As discussed in this subsection, the distribution of $R_{\mathrm{p}}\ $for a
given observed $M_{\mathrm{UV}}\ $bin is influenced by both the underlying
quasar mean-luminosity function and quasar variability, which includes both
the luminosity PDF and the light curve power spectrum. Measuring
$R_{\mathrm{p}}\ $can therefore provide constraints on these critical
properties of the first quasars. However, it is important to note that the
current observed $R_{\mathrm{p}}\ $distribution may be incomplete, and the
selection function is complex. As a result, comparisons between models and
observed $R_{\mathrm{p}}\ $should be interpreted cautiously. Obtaining a large
complete sample of $R_{\mathrm{p}}\ $for reionization-era quasars could
significantly enhance our understanding of quasar variability. From a modeling
perspective, future work will explore how different light curve power spectra
affect the $R_{\mathrm{p}}\ $distribution.

##  5 Conclusions

In this work, we study the proximity zone around a high-redshift quasar
through RT post-processing the lines of sight from a cosmological simulation
with constrained initial conditions. This constrained realization creates a
quasar host halo of $M_{h}=10^{13}\ M_{\odot}$ at $z=6$, more massive than
most halos studied in previous simulations. The simulation also includes
galaxy and black hole formation models, resulting in a variable quasar light
curve, which is more realistic than the widely used lightbulb model. The
simulated light curve exhibits extreme variability, with the changes in
luminosity spanning up to two orders of magnitude around the average value.

By concatenating $R_{\mathrm{p}}$ evolution segments from 21 snapshots
covering $6.0<z\leqslant 7.5$, we capture the variations in proximity zone
sizes with a variable light curve for $\sim 240\ \mathrm{Myr}$ after the
quasar’s activation in a highly ionized IGM. The resultant $R_{\mathrm{p}}$
ranges between $0.5-5\ \mathrm{pMpc}$. We demonstrate that variations in the
light curve contribute an additional scatter, which is separate from the
scatter induced by density variations. The standard deviation in
$R_{\mathrm{p}}$ values caused by each of these effects are approximately
$\sigma(R_{\mathrm{p}})\sim 0.3\ \mathrm{pMpc}$.

Our simulation suggests that in a pre-ionized IGM, the evolution of
$R_{\mathrm{p}}$ traces the variations in the light curve closely with a short
time delay of $\sim 10^{4}\ \mathrm{yr}$. This time lag breaks the
correspondence between the $R_{\mathrm{p}}$ and the contemporaneous
$M_{\mathrm{UV}}$. The $R_{\mathrm{p}}$ is heavily influenced by the magnitude
about $10^{4}$ yr previously, whose difference from the observed
$M_{\mathrm{UV}}$ is uncertain and could be significant.This indicates that
$R_{\mathrm{p}}$ can only be used to infer the quasar episodic lifetime at
best, and does not inform us of the integrated quasar lifetime.

By analyzing the $R_{\mathrm{p}}$ distribution for specific $M_{\mathrm{UV}}$
values, we show that for an individual quasar with a fluctuating light curve,
its proximity zone size increases weakly with brighter instantaneous
magnitude, and diverts from the lightbulb model in a
$M_{\mathrm{UV}}$-dependent way. Compared to the variable light curve, the
lightbulb model underestimates $R_{\mathrm{p}}$ by 13% at the dim end
($M_{\mathrm{UV}}\sim-23.5$), and overestimates the $R_{\mathrm{p}}$ by 30% at
the bright end ($M_{\mathrm{UV}}\sim-26$).

We computed the distribution of $R_{\mathrm{p}}$ based on a set of quasars
sampled from a QLF and found that light curve variability leads to a broad
distribution of $R_{\mathrm{p}}$ at given observed magnitude. Notably,
variable light curves contribute to a group of instantaneously bright quasars
with extremely small proximity zones. These small $R_{\mathrm{p}}$ can hardly
be explained if the quasar light curve stays constant. This shows that it is
necessary to consider the details of light curve variability when
investigating quasar proximity zones.

## Acknowledgements

The authors thank Hy Trac and Nianyi Chen for helpful discussions. HC thanks
the support by the Natural Sciences and Engineering Research Council of Canada
(NSERC), funding reference #DIS-2022-568580. SB acknowledges the funding
support by NASA-80NSSC22K1897. TDM and RAAC acknowledge funding from the NSF
AI Institute: Physics of the Future, NSF PHY-2020295, NASA ATP NNX17AK56G, and
NASA ATP 80NSSC18K101. TDM acknowledges additional support from NASA ATP
19-ATP19-0084, and NASA ATP 80NSSC20K0519.

## Data Availability

The data underlying this article will be shared on reasonable request to the
corresponding author.

## References

  * Bañados et al. (2018) Bañados E., et al., 2018, Nature, 553, 473
  * Bajtlik et al. (1988) Bajtlik S., Duncan R. C., Ostriker J. P., 1988, ApJ, 327, 570
  * Battaglia et al. (2013) Battaglia N., Trac H., Cen R., Loeb A., 2013, ApJ, 776, 81
  * Bolton & Haehnelt (2007a) Bolton J. S., Haehnelt M. G., 2007a, MNRAS, 374, 493
  * Bolton & Haehnelt (2007b) Bolton J. S., Haehnelt M. G., 2007b, MNRAS, 381, L35
  * Bolton et al. (2010) Bolton J. S., Becker G. D., Wyithe J. S. B., Haehnelt M. G., Sargent W. L. W., 2010, MNRAS, 406, 612
  * Bolton et al. (2011) Bolton J. S., Haehnelt M. G., Warren S. J., Hewett P. C., Mortlock D. J., Venemans B. P., McMahon R. G., Simpson C., 2011, MNRAS, 416, L70
  * Bolton et al. (2012) Bolton J. S., Becker G. D., Raskutti S., Wyithe J. S. B., Haehnelt M. G., Sargent W. L. W., 2012, MNRAS, 419, 2880
  * Bondi & Hoyle (1944) Bondi H., Hoyle F., 1944, MNRAS, 104, 273
  * Bosman & Becker (2015a) Bosman S. E. I., Becker G. D., 2015a, MNRAS, 452, 1105
  * Bosman & Becker (2015b) Bosman S. E. I., Becker G. D., 2015b, MNRAS, 452, 1105
  * Carilli et al. (2010) Carilli C. L., et al., 2010, ApJ, 714, 834
  * Cen & Haiman (2000) Cen R., Haiman Z., 2000, ApJ, 542, L75
  * Chen & Gnedin (2021) Chen H., Gnedin N. Y., 2021, ApJ, 911, 60
  * Chen et al. (2023) Chen H., Croft R. A. C., Gnedin N. Y., 2023, MNRAS, 519, 5931
  * Davies et al. (2018) Davies F. B., et al., 2018, ApJ, 864, 142
  * Davies et al. (2020) Davies F. B., Hennawi J. F., Eilers A.-C., 2020, MNRAS, 493, 1330
  * DeGraf et al. (2015) DeGraf C., Di Matteo T., Treu T., Feng Y., Woo J. H., Park D., 2015, MNRAS, 454, 913
  * Di Matteo et al. (2005) Di Matteo T., Springel V., Hernquist L., 2005, Nature, 433, 604
  * Di Matteo et al. (2012) Di Matteo T., Khandai N., DeGraf C., Feng Y., Croft R. A. C., Lopez J., Springel V., 2012, ApJ, 745, L29
  * Ding et al. (2020) Ding X., Treu T., Silverman J. D., Bhowmick A. K., Menci N., Di Matteo T., 2020, ApJ, 896, 159
  * Ding et al. (2022) Ding X., et al., 2022, The Astrophysical Journal, 933, 132
  * Eilers et al. (2017) Eilers A.-C., Davies F. B., Hennawi J. F., Prochaska J. X., Lukić Z., Mazzucchelli C., 2017, ApJ, 840, 24
  * Eilers et al. (2018) Eilers A.-C., Hennawi J. F., Davies F. B., 2018, ApJ, 867, 30
  * Eilers et al. (2020) Eilers A.-C., et al., 2020, ApJ, 900, 37
  * Eilers et al. (2021) Eilers A.-C., Hennawi J. F., Davies F. B., Simcoe R. A., 2021, ApJ, 917, 38
  * Fan et al. (2006) Fan X., et al., 2006, AJ, 132, 117
  * Feng et al. (2016) Feng Y., Di-Matteo T., Croft R. A., Bird S., Battaglia N., Wilkins S., 2016, MNRAS, 455, 2778
  * Feng et al. (2018) Feng Y., Bird S., Anderson L., Font-Ribera A., Pedersen C., 2018, Mp-Gadget/Mp-Gadget: A Tag For Getting A Doi, Zenodo, doi:10.5281/zenodo.1451799
  * Ferrara et al. (2014) Ferrara A., Salvadori S., Yue B., Schleicher D., 2014, MNRAS, 443, 2410
  * Fontanot et al. (2012) Fontanot F., Cristiani S., Vanzella E., 2012, MNRAS, 425, 1413
  * Gabor & Bournaud (2013) Gabor J. M., Bournaud F., 2013, MNRAS, 434, 606
  * Greig et al. (2017) Greig B., Mesinger A., Haiman Z., Simcoe R. A., 2017, MNRAS, 466, 4239
  * Greig et al. (2022) Greig B., Mesinger A., Davies F. B., Wang F., Yang J., Hennawi J. F., 2022, MNRAS, 512, 5390
  * Haiman & Cen (2001) Haiman Z., Cen R., 2001, in Umemura M., Susa H., eds, Astronomical Society of the Pacific Conference Series Vol. 222, The Physics of Galaxy Formation. p. 101
  * Hinshaw et al. (2013) Hinshaw G., et al., 2013, ApJS, 208, 19
  * Hoffman & Ribak (1991) Hoffman Y., Ribak E., 1991, ApJ, 380, L5
  * Hopkins (2013) Hopkins P. F., 2013, MNRAS, 428, 2840
  * Inayoshi et al. (2020) Inayoshi K., Visbal E., Haiman Z., 2020, ARA&A, 58, 27
  * Ishimoto et al. (2020) Ishimoto R., et al., 2020, ApJ, 903, 60
  * Katz et al. (1996) Katz N., Weinberg D. H., Hernquist L., 1996, ApJS, 105, 19
  * Keating et al. (2015) Keating L. C., Haehnelt M. G., Cantalupo S., Puchwein E., 2015, MNRAS, 454, 681
  * Khandai et al. (2015) Khandai N., Di Matteo T., Croft R., Wilkins S., Feng Y., Tucker E., DeGraf C., Liu M.-S., 2015, MNRAS, 450, 1349
  * Khrykin et al. (2017) Khrykin I. S., Hennawi J. F., McQuinn M., 2017, ApJ, 838, 96
  * Khrykin et al. (2021) Khrykin I. S., Hennawi J. F., Worseck G., Davies F. B., 2021, MNRAS, 505, 649
  * King & Nixon (2015) King A., Nixon C., 2015, MNRAS, 453, L46
  * Kravtsov (1999) Kravtsov A. V., 1999, PhD thesis, NEW MEXICO STATE UNIVERSITY
  * Kravtsov et al. (2002) Kravtsov A. V., Klypin A., Hoffman Y., 2002, ApJ, 571, 563
  * Krumholz & Gnedin (2011) Krumholz M. R., Gnedin N. Y., 2011, ApJ, 729, 36
  * Latif et al. (2013) Latif M. A., Schleicher D. R. G., Schmidt W., Niemeyer J. C., 2013, MNRAS, 436, 2989
  * Li (2012) Li L.-X., 2012, MNRAS, 424, 1461
  * Lidz et al. (2006) Lidz A., Oh S. P., Furlanetto S. R., 2006, ApJ, 639, L47
  * Lidz et al. (2007) Lidz A., McQuinn M., Zaldarriaga M., Hernquist L., Dutta S., 2007, ApJ, 670, 39
  * Liu & Liu (2010) Liu M., Liu G., 2010, Archives of Computational Methods in Engineering, 17, 25
  * Madau et al. (2014) Madau P., Haardt F., Dotti M., 2014, ApJ, 784, L38
  * Matsuoka et al. (2018) Matsuoka Y., et al., 2018, ApJ, 869, 150
  * Matsuoka et al. (2019) Matsuoka Y., et al., 2019, ApJ, 883, 183
  * Mazzucchelli et al. (2017) Mazzucchelli C., et al., 2017, ApJ, 849, 91
  * Meiksin et al. (2010) Meiksin A., Tittley E. R., Brown C. K., 2010, MNRAS, 401, 77
  * Miralda-Escudé & Rees (1998) Miralda-Escudé J., Rees M. J., 1998, ApJ, 497, 21
  * Morey et al. (2021) Morey K. A., Eilers A.-C., Davies F. B., Hennawi J. F., Simcoe R. A., 2021, ApJ, 921, 88
  * Mortlock et al. (2011) Mortlock D. J., et al., 2011, Nature, 474, 616
  * Nelson et al. (2015) Nelson D., et al., 2015, Astronomy and Computing, 13, 12
  * Ni et al. (2021) Ni Y., Di Matteo T., Feng Y., 2021, MNRAS, 509, 3043
  * Ni et al. (2022) Ni Y., et al., 2022, MNRAS, 513, 670
  * Novak et al. (2011) Novak G. S., Ostriker J. P., Ciotti L., 2011, ApJ, 737, 26
  * Okamoto et al. (2010) Okamoto T., Frenk C. S., Jenkins A., Theuns T., 2010, MNRAS, 406, 208
  * Oppenheimer et al. (2018) Oppenheimer B. D., Segers M., Schaye J., Richings A. J., Crain R. A., 2018, MNRAS, 474, 4740
  * Read et al. (2010) Read J. I., Hayfield T., Agertz O., 2010, MNRAS, 405, 1513
  * Reed et al. (2017) Reed S. L., et al., 2017, MNRAS, 468, 4702
  * Regan et al. (2019) Regan J. A., Downes T. P., Volonteri M., Beckmann R., Lupi A., Trebitsch M., Dubois Y., 2019, MNRAS, 486, 3892
  * Rudd et al. (2008) Rudd D. H., Zentner A. R., Kravtsov A. V., 2008, ApJ, 672, 19
  * Rumbaugh et al. (2018) Rumbaugh N., et al., 2018, ApJ, 854, 160
  * Satyavolu et al. (2023) Satyavolu S., Kulkarni G., Keating L. C., Haehnelt M. G., 2023, MNRAS, 521, 3108
  * Schawinski et al. (2015) Schawinski K., Koss M., Berney S., Sartori L. F., 2015, MNRAS, 451, 2517
  * Schleicher et al. (2013) Schleicher D. R. G., Palla F., Ferrara A., Galli D., Latif M., 2013, A&A, 558, A59
  * Shakura & Sunyaev (1973) Shakura N. I., Sunyaev R. A., 1973, A&A, 24, 337
  * Shen (2021) Shen Y., 2021, ApJ, 921, 70
  * Smith et al. (2017) Smith A., Bromm V., Loeb A., 2017, Astronomy and Geophysics, 58, 3.22
  * Springel & Hernquist (2003) Springel V., Hernquist L., 2003, MNRAS, 339, 289
  * Springel et al. (2005) Springel V., Di Matteo T., Hernquist L., 2005, MNRAS, 361, 776
  * Stevans et al. (2014) Stevans M. L., Shull J. M., Danforth C. W., Tilton E. M., 2014, ApJ, 794, 75
  * Tepper-García (2006) Tepper-García T., 2006, MNRAS, 369, 2025
  * Vogelsberger et al. (2013) Vogelsberger M., Genel S., Sijacki D., Torrey P., Springel V., Hernquist L., 2013, MNRAS, 436, 3031
  * Vogelsberger et al. (2014) Vogelsberger M., et al., 2014, MNRAS, 444, 1518
  * Volonteri (2010) Volonteri M., 2010, A&ARv, 18, 279
  * Volonteri et al. (2015) Volonteri M., Silk J., Dubus G., 2015, ApJ, 804, 148
  * Wang et al. (2019) Wang F., et al., 2019, ApJ, 884, 30
  * Wang et al. (2020) Wang F., et al., 2020, ApJ, 896, 23
  * Weinberger et al. (2018) Weinberger R., et al., 2018, MNRAS, 479, 4056
  * Worseck et al. (2014) Worseck G., et al., 2014, MNRAS, 445, 1745
  * Worseck et al. (2021) Worseck G., Khrykin I. S., Hennawi J. F., Prochaska J. X., Farina E. P., 2021, MNRAS, 505, 5084
  * Wyithe et al. (2005) Wyithe J. S. B., Loeb A., Carilli C., 2005, ApJ, 628, 575
  * Yang et al. (2020) Yang J., et al., 2020, ApJ, 897, L14
  * Šoltinský et al. (2023) Šoltinský T., Bolton J. S., Molaro M., Hatch N., Haehnelt M. G., Keating L. C., Kulkarni G., Puchwein E., 2023, MNRAS, 519, 3027
  * van de Weygaert & Bertschinger (1996) van de Weygaert R., Bertschinger E., 1996, MNRAS, 281, 84

##  Appendix A $R_{\mathrm{p}}$ for the lightbulb model

Figure 11: Based on a lightbulb model fixed at $\dot{N}=1\times 10^{57}\
\mathrm{s}^{-1}$, the evolution of the proximity zone size $R_{\mathrm{p}}$ as
a function of the quasar turn-on redshift for $t_{\mathrm{evol}}=0.03\
\mathrm{Myr}$ (blue), $1\ \mathrm{Myr}$ (orange), and $30\ \mathrm{Myr}$
(green). The dots are results from the CR simulation, whose errorbars
represent the 16-84th scatter, and the solid curves are the median
$R_{\mathrm{p}}$ modeled by Chen & Gnedin (2021) with the CROC simulation.

In order to compare to previous work, we simulate proximity zone sizes for a
lightbulb model with $\dot{N}=10^{57}\ \mathrm{s}^{-1}$ based on the CR
simulation used in this paper. For each line of sight, we turn on the quasars
at each snapshot redshift $z_{\mathrm{snap}}$ and evolve the radiation for
$30\ \mathrm{Myr}$. In Fig. 11, we show the measured $R_{\mathrm{p}}$
distribution after $t_{\mathrm{evol}}=0.03\ \mathrm{Myr}$ (blue),
$t_{\mathrm{evol}}=1\ \mathrm{Myr}$ (orange), and $t_{\mathrm{evol}}=30\
\mathrm{Myr}$ (green) from the CR simulation using the dots, with the errorbar
indicating the 16-84th percentiles of scatter. We also portray the median
$R_{\mathrm{p}}$ modeled by Chen & Gnedin (2021) with the solid curves.

Fig.11 demonstrates that our simulation generates similar $R_{\mathrm{p}}$ to
Chen & Gnedin (2021) at $t_{\mathrm{evol}}=0.03\ \mathrm{Myr}$ and
$t_{\mathrm{evol}}=30\ \mathrm{Myr}$, while predicting smaller proximity zones
at $t_{\mathrm{evol}}=1\ \mathrm{Myr}$. This difference in results probably
stems from the more massive quasar host halo ($\sim 10^{13}\ M_{\odot}$ rather
than $\gtrsim 1.5\times 10^{11}\ M_{\odot}$) we considered in our simulation,
as shown in Fig. 2. With a given ionization fraction, the proximity zones in
the high-density environment grow slower while reaching a similar maximum
$R_{\mathrm{p}}$ with those in the low-density environment (see Fig. 7 in
Keating et al. (2015)). At high redshift $z>7$, our results seem to evolve
faster than those produced by Chen & Gnedin (2021) at $t_{\mathrm{evol}}=0.03$
Myr. This is because the CR simulation has completed reionization at $z=8$
(see Section 2.2), while CROC still has a large volume of neutral hydrogen in
the IGM around $z=7.3$ ($x_{\mathrm{H\,\textsc{i}}}=0.13$).

##  Appendix B Influence of background evolution

In this section, we explore the uncertainty in the simulated $R_{\mathrm{p}}$
caused by the stitching procedure applied in Section 3.1. When we calculate
the evolution of $R_{\mathrm{p}}$, we only update the properties of the IGM
(including the density and velocity field) at $z_{\mathrm{snap}}$, which
breaks the continuous evolution of IGM and might make our results differ from
reality.

To estimate the uncertainty caused by this method of background updating, we
use the semi-analytical model proposed in Davies et al. (2020) and calculate
the largest change in $R_{\mathrm{p}}$ when the background is updated, which
roughly equals to the largest error caused by the stitching procedure if we
assume the simulated $R_{\mathrm{p}}$ at newly-updated $z_{\mathrm{snap}}$ is
accurate. Taking into account that the IGM environments adopted in Davies et
al. (2020) are different from ours (such as the density, ionization fraction),
we first re-scale the semi-analytical model for our simulation:

$\displaystyle R_{\mathrm{p}}$
$\displaystyle=3.1544\left(\frac{\Gamma_{\mathrm{bkg}}^{\mathrm{H\,\textsc{i}}}}{2.455\times
10^{-13}\,\text{s}^{-1}}\,\right)^{-1/2}\left[\left(\frac{\tau_{\mathrm{bkg}}}{2.3}\right)^{1/\alpha}-1\right]^{-1/2}$
(10) $\displaystyle\times\left(\frac{\dot{N}}{1.3197\times
10^{56}\,\text{s}^{-1}}\right)^{1/2}\left(\frac{1+z}{7.5}\right)^{-3/2}\
\mathrm{Mpc},$

where
$\tau_{\mathrm{bkg}}=5.678\left(\Gamma^{\mathrm{H\,\textsc{i}}}_{\mathrm{bkg}}/\left[2.5\times
10^{-13}\,\mathrm{s}^{-1}\right]\right)^{-\alpha}$ and $\alpha=0.5486$. In
this model, all of the influences of the background are encoded in
$\Gamma_{\mathrm{bkg}}^{\mathrm{H\,\textsc{i}}}$, which is the collisional
ionization rate of hydrogen caused by the background radiation. We calculate
the value of $\Gamma_{\mathrm{bkg}}^{\mathrm{H\,\textsc{i}}}$ in the same way
as Chen & Gnedin (2021) (see their equation 3), which is determined by the
density, temperature and the ionization fraction of the IGM without the quasar
radiation. The largest time gap between two consecutive snapshots is $25\
\mathrm{Myr}$, which exists between $z_{\mathrm{snap}}=7.500$ and
$z_{\mathrm{snap}}=7.305$, corresponding to a change in
$\Gamma_{\mathrm{bkg}}^{\mathrm{H\,\textsc{i}}}$ of $\lesssim 2\times
10^{-14}\ \mathrm{s}^{-1}$. Based on equation 10, such change in
$\Gamma_{\mathrm{bkg}}^{\mathrm{H\,\textsc{i}}}$ results in an uncertainty of
$\Delta R_{\mathrm{p}}\sim 0.05\ \mathrm{Mpc}$. Compared with the
$R_{\mathrm{p}}$ scatter stemming from the density field fluctuations and the
light curve variation, both of which are $\sim 0.3$ Mpc as mentioned in
Section 3.1, such uncertainty is insignificant. Hence, the background
evolution between two adjacent snapshots has negligible effects on our
simulation.

##  Appendix C The Influence of Variation Amplitudes in Quasar Light Curves

Here we briefly discuss how the variation amplitude in the light curve
influences the amount of discrepancy between the variable light curve and the
lightbulb models. We build toy models for three light curves varying around
$M_{\mathrm{UV}}=-25.5$ with different variation amplitudes: $\Delta
M_{\mathrm{UV}}=\pm 1$, $\Delta M_{\mathrm{UV}}=\pm 2$, and $\Delta
M_{\mathrm{UV}}=\pm 3$, respectively. These light curves have zigzag shapes,
whose magnitudes change linearly with time from $-25.5+|\Delta
M_{\mathrm{UV}}|$ to $-25.5-|\Delta M_{\mathrm{UV}}|$, and then drop linearly
to $-25.5+|\Delta M_{\mathrm{UV}}|$ to complete a period. They all have a
fluctuation period of $0.1\ \mathrm{Myr}$ and we evolve them for $0.5$ Myr. We
calculate the $R_{\mathrm{p}}$ generated by these light curves on the lines of
sight from $z_{\mathrm{snap}}=6.162$ snapshot with a time resolution of
$2\times 10^{3}\ \mathrm{yr}$. By comparing the resultant $R_{\mathrm{p}}$
distributions, we are allowed to see the influence of variation amplitudes.

In Fig. 12, we depict the one-dimensional histogram for the $R_{\mathrm{p}}$
sampling from the variable light curves within $-26.5<M_{\mathrm{UV}}<-25.5$
(solid curves) as well as those produced by a lightbulb model fixed at
$M_{\mathrm{UV}}=-26$ (red shaded area). The mean $R_{\mathrm{p}}$ for the
lightbulb model is $2.93$ pMpc. While for the variable light curves with
$\Delta M_{\mathrm{UV}}=\pm 1,\ \pm 2,\ \pm 3$, the mean $R_{\mathrm{p}}$ are
$2.15$ pMpc, $2.19$ pMpc, and $2.32$ pMpc, respectively. All the variable
light curves generate an $R_{\mathrm{p}}$ distribution peaking at smaller
values. Furthermore, larger variation magnitudes result in a slightly larger
shift from the lightbulb model.

Figure 12:  The $R_{\mathrm{p}}$ distribution generated by the variable light
curves with $-26.5<M_{\mathrm{UV}}<-25.5$ (solid curves) and by a lightbulb
model fixed at $M_{\mathrm{UV}}=-26$ (red shade area). The three variable
light curves have a fixed fluctuation period of $10^{5}$ yr, and vary around
$M_{\mathrm{UV}}=-25.5$ with the variation amplitudes of $\Delta
M_{\mathrm{UV}}=\pm 1$ (blue), $\pm 2$ (orange), $\pm 3$ (green),
respectively. All the light curves evolve $0.5\ \mathrm{Myr}$. The
$R_{\mathrm{p}}$ are computed on the $z_{\mathrm{snap}}=6.162$ snapshot and
are sampled with a time resolution of $2\times 10^{3}\ \mathrm{yr}$.