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

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_{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_{p}$) ranges between $0.5-5$ pMpc. We show that the light curve variability causes a similar degree of scatter in $R_{p}$ as does the density fluctuation, both of which result in a standard deviation of $\sim 0.3$ pMpc). The $R_{p}$ traces the light curve fluctuations closely but with a time delay of $\sim 10^4\ \mathrm{yr}$, breaking the correspondence between the $R_{p}$ and the contemporaneous $M_{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_{p}$ in the later part of cosmic reionization. Compared with the variable light curve, a constant light curve underestimates the $R_{p}$ by 13% at the dim end ($M_{UV}\sim -23.5$), and overestimates the $R_{p}$ by 30% at the bright end ($M_{UV}\sim -26$). By calculating the $R_{p}$ generated by a number of quasars, we show that variable light curves predict a wider $R_{p}$ distribution than lightbulb models, and readily explain the extremely small $R_{p}$ values that have been observed.


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  > 6, these are the only regions where we can observe non-zero Lyman  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  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. 2017Eilers et al. , 2020;;Mazzucchelli et al. 2017;Ishimoto et al. 2020).Fan et al. (2006) compiled the first large sample of  ≳ 6 quasar spectra and measured a proximity zone size  p ∼ 8 Mpc at  = 6 for quasars of magnitude  1450 = 27.Carilli et al. (2010) analyzed the proximity zone sizes of 27 quasars with more accurate ★ E-mail: yihaoz@andrew.cmu.eduredshift measurements.In the last decade, the number of high redshift quasar spectra with well-measured Lyman  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 ⩽  ⩽ 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⩽  1450 ⩽ −22.83) quasars at  ∼ 6.An unexpected result from these observations obtained in recent years is that some quasars display very small proximity zones, such as the  p ∼ 0.37 Mpc measured in Eilers et al. (2021).One possible explanation for several small  p seen in  ≳ 7 quasar spectra is that a large amount of hydrogen at  ∼ 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. 2017Greig et al. , 2022)).However, the population of small  p quasars at  ∼ 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  p and the intrinsic properties of quasars (e.g., the lumi-nosity 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): where  is the emitted ionizing photon rate,  q is the quasar lifetime,  H is the hydrogen number density and  H 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  max p which is independent on the neutral fraction: pMpc. (2) Here  lim = − ln( lim ) is the limiting Lyman  optical depth (typically  lim is set to be 0.1) and  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  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  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  ∼ 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  ⊙ ) and discovered that  p depends weakly on the mass of the host halo, but strongly on direction, leading to a large  p scatter across lines of sight.Chen & Gnedin (2021) found that about 1% ∼ 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 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 ∼ 1 mag within ∼ 15 yr (Rumbaugh et al. 2018).Furthermore, lightbulb models can hardly explain the small proximity zones observed around luminous quasars at  ∼ 6.Although Eilers et al. (2017Eilers et al. ( , 2021) ) offer the suggestion that these small proximity zones are from young quasars that just turned on for  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  duty or episodic lifetime  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  p based on their time-dependent semi-analytical model.They found that a variable light curve yielded a  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 Massive-Black 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;Ding et al. 2022), the observed quasar luminosity functions (QLF), and the correlation length  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.

Constrained realization simulation
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 ≳ 500 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 ( ∼ 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  BH ∼ 10 9  ⊙ at  = 6.
The CR simulation has a box size of 20 ℎ −1 cMpc per side.The ICs are generated using the most recent implementation of gaussianCR 1 , and the simulation is run to  = 6 with the massively parallel cosmological smoothed-particle hydrodynamics (SPH) simulation software MP-Gadget (Feng et al. 2018).The peak height of 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  = 7.5 (upper) and  = 6.0 (middle).For regions with high density, the color hue is set by the hydrogen neutral fraction  H 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 ℎ −1 Mpc, and the two right panels depict the zoomed-in central regions of diameter 1 ℎ −1 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 Δ g (= / ρ) along all 48 directions at  = 7.5 in grey, with the four examples in the upper panels highlighted in blue, yellow, green, and red, respectively.The gas density contrast (Δ 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 pkpc in length, and all the data are drawn from 0.1 ∼ 2 pMpc regions from the quasar in the  = 6.5 snapshot.
our quasar-hosting halo is 5  0 with a peculiar velocity  = 0 km s −1 , where  0 is the density variance after convolution with a Gaussian kernel of size  G = 1 ℎ −1 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).

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  DM :  w =  w  DM , where  w is the wind speed, and the dimensionless parameter  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 subgrid model developed in Springel et al. (2005) and Di Matteo et al. (2005).Black holes are seeded with an initial seed mass of  seed = 5 × 10 5  ⊙ ℎ −1 in halos with a mass larger than 5 × 10 10  ⊙ ℎ −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): where  s is the local sound speed,  BH is the gas density around the quasar, and  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  Edd .Therefore the black hole accretion rate  BH is determined by With a radiative efficiency  = 0.1 (Shakura & Sunyaev 1973), the black hole radiates with a bolometric luminosity  bol proportional to the accretion rate:  bol =   BH 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% ∼ 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  = 8 (Fig. 2 in Feng et al. (2016)).Consequently, the gas in the snapshots we use herein ( ⩽ 7.5) is originally highly ionized.

Line of sight gas densities
The projected gas density fields around the quasar at  = 7.5 and  = 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  H i , from red to blue indicating ionized ( H i ≲ 10 −4 ) to neutral ( H i ∼ 0.1), as shown by the color bars.The left two panels display the regions 10 ℎ −1 Mpc in width centered on the quasar, while the right two panels depict the central zoomed-in regions of width 1 ℎ −1 Mpc.
We use HEALPY 2 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  () (e.g., density, velocity, ionization fraction) at position  on the line of sight (Liu & Liu 2010): where  is the sum over all the neighboring gas particles within the smoothing length , and  is the quintic density kernel.  ,   , and   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 ℎ −1 Mpc.The sightlines are drawn significant off-axes (not parallel to the ,  or  axes), ensuring that none of them travel through the massive halo again.  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 Δ g (= / ρ; where ρ is the mean gas density of the Universe) at  = 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 Δ 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 ℎ −1 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 ℎ −1 ckpc and a peak resolution of 100 pc.All the sightlines are drawn from halos with dark matter mass larger than 1.5 × 10 11  ⊙ .Another major difference is that the CROC simulation models reionization by star particles self-consistently, resulting in a volume-weighted neutral fraction ⟨ H i ⟩ v = 0.13 at  = 7.33 and ⟨ H i ⟩ v < 6.7 × 10 −4 after  = 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  = 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 Δ 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 ∼ 10 13  ⊙ at  = 6, is much more massive than the halos selected in Chen & Gnedin (2021), whose halo masses are ≳ 1.5×10 11  ⊙ .In fact, the halo chosen in this work is more massive than almost all the halos in previous simulations of proximity zones, for example,  h = 2.5 × 10 12 ℎ −1  ⊙ in Keating et al. (2015),  h ≳ 10 11.5  ⊙ in Davies et al. (2020), and  h < 10 12  ⊙ in Satyavolu et al. (2023).We test how this difference in the surrounding density field affects the resultant  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.
We compute  for the quasar and show the light curve (black curve) in Fig. 3, compared with  Edd (red curve), which represents the photon number rate corresponding to 2  Edd , the upper limit of the accretion rate.The yellow curve indicates the photon number rate averaged using a 5 Myr top-hat kernel (  mean ).In this work, we focus on the light curve within 6 <  ⩽ 7.5, a period during which the average quasar luminosity (  mean ) reaches a plateau.The quasar has a mean magnitude of  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 .We display the PDF (left panel) and the dimensionless power spectrum  () (right panel) of /  mean in Fig. 4, which shows that the quasar experiences variation in  spanning over two orders of magnitude.The characteristic fluctuation timescale, which is indicated by the peak of the power spectrum, is around  = 1 Myr.

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  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  H i and gas temperature output from the RT simulation, in conjunction with the velocity and the density field from the CR simulation.

Evolution of 𝑅 p
We present the post-processed spectra obtained from the RT code for a single line of sight at  = 7 in Fig. 5.The gas density contrast Δ g is shown in the first row and the temperature prior to quasar activation (i.e.,  evol = 0 Myr) is shown in the third row with the black dashed line.The transmitted flux and IGM temperature at  evol = 1 Myr (blue), 5 Myr (orange), 15 Myr (green), 20 Myr (red), and 25 Myr (purple) are depicted in the second and third rows, respectively, with the corresponding  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  marked by stars.The traditional definition of the proximity zone edge is the point where the Lyman  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  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 , the higher the corresponding flux is.This implies that the levels of the Lyman  flux, and consequently the proximity zone size  p , are primarily determined by  while showing no correlation with  evol .
The profiles of H i/He i/He ii fractions (solid blue/orange/green curves) at  evol = 5 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(Bolton et al. , 2012;;Meiksin et al. 2010).This phenomenon also enhances the size of the Lyman  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  p , continues to expand monotonically with the increasing  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: In the third panel, the black dashed curve depicts the temperature at  evol = 0 Myr (i.e., no quasar radiation).In the bottom row, H i/He i/He ii fraction (solid blue/orange/green curves) at  evol = 5 Myr are compared with the background ionization fractions(dashed curves).
6 <  snap,i ⩽ 7.5 (1 ⩽  ⩽ 21;  snap,i+1 <  snap,i ), we conduct an RT calculation on the th snapshot for a period of where  age () is the age of the universe at redshift ,  snap,i+1 is the redshift of the subsequent snapshot and  snap,22 = 6.When we concatenate the next snapshot at  snap,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  p during  age ( snap,i ) ∼  age ( 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  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  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 ∼ 240 Myr after the quasar's activation.Fig. 6 presents the evolution of  p with a time resolution of Δ evol = 0.25 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  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  snap,i .Since the reionization is virtually completed at  = 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  p for an individual direction within short time frames.The extent of this variability in  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 Δ p ∼ 1.7 Mpc for sightline 7 within 0.25 Myr, while for sightline 6, the changes are nearly all less than Δ p ∼ 0.5 Mpc.
We can quantify the influence of the underlying density fluctuations by computing the standard deviation of proximity zone sizes (  p ) across sightlines for the same redshift.We find the mean of this standard deviation (   p ) is 0.28 Mpc across the entire redshift range.On the other hand, to show the influence of light curve variability, we compute the mean of  p ( R p ) across all sightlines at a  2020) are displayed by the black and green dots, which are made for quasars at 6.0 ⩽  ≲ 6.5.The minimum and maximum simulated values for 6.0 ⩽  ≲ 6.5 are depicted with the red dashed curves for comparison.The upper and the right panels present the marginal distribution of  UV and  p , respectively.The blue dash-dot lines represent the mean values of the  UV and the  p for 6.0 ⩽  ⩽ 7.5:  UV,mean = −24.8, p,mean = 1.37 pMpc.
specific redshift, and then calculate the standard deviation of these mean  p values for different redshifts ( ⟨Rp⟩ ), which is 0.33 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  p values.

𝑅 p − 𝑀 UV scaling relation
In Fig. 7, the lower left panel displays the proximity zone size  p as a function of the quasar's instantaneous magnitude  UV .The median  p across the entire redshift range (6.0 ⩽  ⩽ 7.5) for 48 lines of sight is represented by the blue solid curve.The best power-law fit  p −  UV scaling relation is found to be log  p ∝ −0.13  uv , i.e.,  p ∝  0.32 (blue dotted curve).The upper and the right panels show the marginal distribution of  UV and  p .The blue dash-dot lines represent the mean values of  UV and  p :  UV,mean = −24.8, p,mean = 1.37 pMpc.The  UV histogram shows that the quasar is relatively faint ( UV ≳ −25) over most of its lifetime, and only becomes luminous ( UV ≲ −26) occasionally.
For comparison with observational results, in Fig. 7 we plot the observed proximity zone size for quasars in the range 6 ⩽  ≲ 6.5 measured by Eilers et al. (2017Eilers et al. ( , 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  p range for most of the observed quasars with  UV > −26.The inability to yield  p < 0.5 Mpc could be due to the failure to resolve LLSs or DLAs.On the bright end, our simulation predicts smaller  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  p −  UV scaling relation exhibits a flatter trend compared to the  p ∝  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  UV for our simulated quasar is non-uniformly distributed across a broad redshift range, as shown in the upper panel of Fig. 7.The  p −  UV scaling relation is redshift-dependent, which can be seen in Fig. A1, where the lower redshift environment typically produces more extensive proximity zones.According to our simulation, the optimal fit is log  p ≈ log(3.20 Mpc) + 0.36 [−0.4(UV + 27)] at 6.0 ⩽  ⩽ 6.5 and log  p ≈ log(2.41Mpc)+0.41[−0.4(UV + 27)] at 7.3 ⩽  ⩽ 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  p on the instantaneous magnitude probably originates from the variation in the light curve, which breaks the correspondence between the  p and the contemporaneous  UV as we discuss in Section 4.1.

Response of 𝑅 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  p to the quasar light curve.We plot a 5 Myr duration in the evolution of  p starting from  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 kyr time spans around four peaks in the light curve (red dotted lines in subplots) and their corresponding  p peaks (blue dotted lines in subplots).We label the time lags between the light curve peaks and the  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:  H i eq = 1/Γ H i ∼ 10 4 yr, where Γ H i is the photoionization rate of hydrogen (Bolton & Haehnelt 2007b;Eilers et al. 2018;Davies et al. 2020).This illustrates that  p traces the fluctuations in the light curve closely but with a short delay of ∼ 10 4 yr, which breaks the correspondence between  p and the contemporaneous  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 ×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(Eilers et al. , 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  p delay time, the quasar's luminosity has already changed significantly.Moreover, there are many periods during which the quasar luminosity evolves In the subplots, we label the time lag between the light curve peaks and the  p peaks, which are ∼ 10 4 yr.
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  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  eq ∼ 10 6 yr (Worseck et al. 2021;Khrykin et al. 2017Khrykin et al. , 2021;;Šoltinský et al. 2023;Chen et al. 2023).

Influence of light curve variation
In this section, we discuss how the rapid fluctuation in the light curve in our simulation affects the resultant  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  = 6.162 with a time span of 30 Myr, short enough that the density evolution is negligible.This excerpted light curve has a mean magnitude of  UV,mean = −24.72 and a mean ionizing photon rate of  mean = 1.67 × 10 56 s −1 .In contrast to this variable light curve, we construct another light curve with a constant flux of the same  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  p at different  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  p groups.For the lightbulb model,  p remains nearly unchanged after rapid growth during the first ∼ 1 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  p distributions as a function of instantaneous magnitude  UV for 10 <  evol < 30 Myr, during which the lightbulb model reaches a stable stage and produces a similar  p .The blue pixels represent the  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  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  p ∼ 1.6 pMpc across the entire magnitude range, the variable light curve yields a scatter ( R p = 0.46 pMpc) 28% larger than that of the  mean lightbulb ( R p = 0.36 pMpc).The  R p values for the variable light curve encompass contributions from both light curve variation and underlying density field fluctuation.On the other hand, the  R p for the lightbulb model is almost totally attributed to the density differences.In addition to the lightbulb model fixed at  UV,mean (yellow dot), we also simulate the lightbulb with different magnitudes (black dots):  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  UV , and the high luminosity magnifies the variance between directions, i.e., brighter  UV leads to an increase in   p .
An important feature illustrated by the middle panel of Fig. 9 is that the median  p from the variable light curve at a specific magnitude coincides with the lightbulb model only around  UV,mean = −24.72,while it tends to yield smaller  p compared to the lightbulb when  UV <  UV,mean , and conversely, larger  p when  UV >  UV,mean .This is more clearly demonstrated in the subplot, where we show the difference between the median  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  UV,mean .More specifically, in the right panel of Fig. 9 we compare the one-dimensional  p distributions within a narrow  UV bin for 10 Myr <  evol < 30 Myr.On the bright end ( UV ∼ −26), the lightbulb model predicts a median  p 30% larger than that produced by the variable light curve (2.94 pMpc versus 2.27 pMpc).While on the dim end ( UV ∼ −23.5), the lightbulb model yields a median  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  p for this quasar.The standard deviations of  p for both the variable light curve and the lightbulb model are ∼ 0.6 pMpc around  UV = −26, and ∼ 0.  2020) noticed that the  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  p is governed by the entire light curve within the most recent ∼ 10 4 yr, rather than the contemporaneous instantaneous luminosity.As depicted by the PDF in Fig. 4, the distribution of  is essentially Gaussian centering around  mean , which implies that the  value 10 4 yr preceding a bright or a dim point in the light curve is probably close to  mean , producing a  p close to that given by a lightbulb model fixed at  UV,mean .Furthermore, the higher luminosities correspond to more significant variations since the variation amplitude generally equals  UV,mean −  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  UV -dependent way.Such divergence accounts for the difference between our predicted  p and the observational measurements at  UV < −26 shown in Fig. 7: our simulated light curve generally has a lower luminosity, which makes the  p in this magnitude range smaller; while the observed quasars with  UV < −26 probably have larger overall luminosity, and so generate large proximity zones.

Implications of the observed 𝑅 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  ★ UV , our quasar is a relatively faint one ( UV =-24.8) during the period  = 7.5 ∼ 6.However, it still has a 16% chance to be caught in a relatively bright phase with  UV <-25.5.In such a bright phase, the distribution of  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  p distribution at given observed magnitude  UV (( p | UV )).
The distribution of ( p | UV ) can be formulated as the following conditional distribution: , where ( is the probability that the quasar of  ★ UV at the observed magnitude  UV displays a proximity zone size of  p . To calculate such a distribution, certain assumptions need to be made.The first term ( ★ UV ) is similar to the quasar luminosity function, but it is for the mean magnitude  ★ UV instead of the observed luminosity.Measuring such a  ★ UV directly is challenging.For simplicity, here we assume that ( ★ UV ) is equal to the observed quasar luminosity function (QLF) measured by Matsuoka et al. (2018) To estimate ( UV | ★ UV ), 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  ★ UV : with a fixed scatter  = 0.7 for all the light curves.Finally, to model ( p | UV ,  ★ UV ), we make the following assumptions: (1) ( p | UV ,  ★ UV ) has the same shape as the  p distribution generated by the constant light curve with magnitude fixed at  UV , which we label as  lb,  UV ( p ), but is shifted towards a different mean  p with the amount .(2) the value of  is only determined by  UV −  ★ UV and is independent of the specific values of  ★ UV and  UV .Hence, the probability of  p for a given variable light curve is The first assumption is motivated by the right panel of Fig. 9, which shows that for an individual quasar, the  p distributions generated by the variable light curve and the lightbulb have similar shapes: they have roughly the same scatter but different mean  p .To formulate ( UV −  ★ UV ), we use our results in Section.4.2 as a guideline, i.e., we linearly interpolate the difference between the median  p of the lightbulb models and the variable light curve as a function of the magnitude difference  UV −  ★ UV , 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  = 1000 quasar samples to create the final distribution ( p | UV ).We present the  p distribution with −26.5 <  UV < −25.5 calculated by this model in Fig. 10 (red solid curve).As a comparison, we plot the combined  p distribution produced by the lightbulb models for 10 <  evol < 30 Myr fixed at  UV = −25.5,−26, −26.5 (red shaded area).We consider three lightbulb models rather than only the one with  UV = −26 to show the whole  p range reached by the lightbulb in this  UV bin, which spans the minimum  p , reached when  UV = −25.5, to the maximum  p , for  UV = −26.5.These three lightbulb models are sampled based on the QLF (equation 7).We also show the  p values for the observed quasars with −26.5 <  UV < −25.5 measured by Eilers et al. (2017Eilers et al. ( , 2020)); Ishimoto et al. (2020) in Fig. 10 (black curve).It is evident that although a variable light curve produces a similar maximum value ( p,max ∼ 5.5 Mpc) as the lightbulb model, it can also result in much smaller proximity zones ( p,min ∼ 0.2 Mpc).Therefore, it readily explains the observations with small proximity zones.These extremely small  p values are generated by the quasar with low averaged luminosity (i.e., large  ★ UV ).The considerable difference between the observed instantaneous  UV , which is ∼ −26 in this case, and  ★ UV moves the  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  ★ UV quasars over the bright ones, the overall  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  p distribution and predicted  p from different quasar light curve models.They concluded that the lightbulb model with long episodic quasar lifetime (≥ 1 Myr) leads to a  p distribution consistent with the observed  p distribution, while their toy quasar light curve with variation does not.On the contrary, our simulated quasar light curve results in a  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  = 0.23 and -value= 0.54 when comparing the observed  p distribution and that from the lightbulb model.On the other hand, we find  = 0.34 and -value= 0.12 between the observed  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  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  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. (2017Eilers et al. ( , 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  p for a given observed  UV bin is influenced by both the underlying quasar meanluminosity function and quasar variability, which includes both the luminosity PDF and the light curve power spectrum.Measuring  p can therefore provide constraints on these critical properties of the first quasars.However, it is important to note that the current observed  p distribution may be incomplete, and the selection function is complex.As a result, comparisons between models and observed  p should be interpreted cautiously.Obtaining a large complete sample of  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  p distribution.

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  ℎ = 10 13  ⊙ at  = 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 ex-hibits extreme variability, with the changes in luminosity spanning up to two orders of magnitude around the average value.
By concatenating  p evolution segments from 21 snapshots covering 6.0 <  ⩽ 7.5, we capture the variations in proximity zone sizes with a variable light curve for ∼ 240 Myr after the quasar's activation in a highly ionized IGM.The resultant  p ranges between 0.5 − 5 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  p values caused by each of these effects are approximately ( p ) ∼ 0.3 pMpc.
Our simulation suggests that in a pre-ionized IGM, the evolution of  p traces the variations in the light curve closely with a short time delay of ∼ 10 4 yr.This time lag breaks the correspondence between the  p and the contemporaneous  UV .The  p is heavily influenced by the magnitude about 10 4 yr previously, whose difference from the observed  UV is uncertain and could be significant.This indicates that  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  p distribution for specific  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  UV -dependent way.Compared to the variable light curve, the lightbulb model underestimates  p by 13% at the dim end ( UV ∼ −23.5), and overestimates the  p by 30% at the bright end ( UV ∼ −26).
We computed the distribution of  p based on a set of quasars sampled from a QLF and found that light curve variability leads to a broad distribution of  p at given observed magnitude.Notably, variable light curves contribute to a group of instantaneously bright quasars with extremely small proximity zones.These small  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.
Figure2.The gas density contrast (Δ g ) PDF for the lines of sight from the CR simulation used in this work (blue), and those from the CROC simulation (orange) used inChen & Gnedin (2021).Each pixel is 4 pkpc in length, and all the data are drawn from 0.1 ∼ 2 pMpc regions from the quasar in the  = 6.5 snapshot.

Figure 3 .Figure 4 .
Figure3.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  mean computed over a time kernel of 5 Myr.The upper axis indicates the age of the Universe at a given redshift.

Figure 5 .
Figure5.The spectra computed by the RT post-processing code for one line of sight drawn from the snapshot at  = 7. Panels from top to bottom show the gas density contrast Δ g (= / ρ), the Lyman  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  evol = 1 Myr (blue), 5 Myr (orange), 15Myr (green), 20 Myr (red), and 25 Myr (purple).In the second row, a subplot provides an overview of the light curve along with the  for the selected  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  evol = 0 Myr (i.e., no quasar radiation).In the bottom row, H i/He i/He ii fraction (solid blue/orange/green curves) at  evol = 5 Myr are compared with the background ionization fractions(dashed curves).

Figure 6 .
Figure 6.The evolution of the proximity zone size  p with a time resolution of Δ evol = 0.25 Myr.The top 6 panels depict the  p (blue curve; left y-axis) for 6 lines of sight, compared with the quasar light curve  (red curves, right y-axis).The vertical grey dash lines denote the redshift of available snapshots  snap,i , and the upper axis indicates the age of the Universe at a given redshift.The bottom panel shows the average proximity zone  p (blue curve) and the 16-84th percentile scatter  Rp (blue shaded area) across 48 directions for the specific redshift.The standard deviation of  p , which represents the scatter caused by the light curve variation, is 0.33 Mpc.And the mean of  Rp for the entire redshift range, indicating the influence of density fluctuations, is 0.28 Mpc.

Figure 7 .
Figure 7.The dependence of proximity zone size  p on quasar instantaneous magnitude  UV .Solid curves in the lower left panel show median  p as a function of  UV for 6.0 ⩽  ⩽ 7.5 (blue), which is the entire redshift range for this study, and 6.0 ⩽  ⩽ 6.5 (red).The blue dotted curve indicates the best power-law fit  p −  UV scaling relation:  p ∝  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 ⩽  ≲ 6.5.The minimum and maximum simulated values for 6.0 ⩽  ≲ 6.5 are depicted with the red dashed curves for comparison.The upper and the right panels present the marginal distribution of  UV and  p , respectively.The blue dash-dot lines represent the mean values of the  UV and the  p for 6.0 ⩽  ⩽ 7.5:  UV,mean = −24.8, p,mean = 1.37 pMpc.

Figure 8 .
Figure8.The averaged evolution of the proximity zone size  p (blue curve; left y-axis) compared with the quasar light curve  (red curve; right y-axis) for 5 Myr following  snap = 7.160.The blue shaded region represents the 16-84th percentile scatter.The subplots zoom in to 80 kyr time spans around four peaks in the light curve (red dotted lines in subplots) and their corresponding  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  p peaks, which are ∼ 10 4 yr.

Figure 9 .
Figure 9. Left panel: the evolution of the mean proximity zone sizes  p with the variable light curve (blue curves) or a lightbulb model fixed at  UV,mean = −24.72 (yellow curves).The quasar is assumed to turn on at  = 6.162 and last for 30 Myr.The blue shaded region and the yellow dashed lines indicate the 68% scatter for the two  p groups.Middle panel: two-dimensional distribution of  p and instantaneous  UV for 10 <  evol < 30 Myr.The blue pixels represent the  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  p generated by the lightbulb models fixing the magnitude at  UV = −23.5,−24, −24.5, −25, −25.5, −26, −26.5, −27 (black dots) and  UV,mean = −24.74(yellow dot).The thick/thin errorbars correspond to 68%/95% confidence intervals.The subplot displays the difference between the median  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  UV,mean .Right panel: one-dimensional  p distribution generated by the part of the variable light curve with −26.5 <  UV < −25.5 (red curve), −24 <  UV < −23 (blue curve), and by a lightbulb model with  UV = −26 (red shaded area),  UV = −23.5 (blue shaded area) for 10 Myr <  evol < 30 Myr.

Figure 10 .
Figure10.The  p distribution with −26.5 <  UV < −25.5 produced by the variable light curves (red solid curve), and the lightbulb models with 10 <  evol < 30 Myr (red shaded area), compared with the observed  p for the quasars within the same magnitude bin measured inEilers et al. (2017Eilers et al. ( , 2020));Ishimoto et al. (2020).