Near-zone sizes and the rest-frame extreme ultraviolet spectral index of the highest redshift quasars

Abstract

The discovery of quasars with redshifts higher than six has prompted a great deal of discussion in the literature regarding the role of quasars, both as sources of reionization, and as probes of the ionization state of the intergalactic medium (IGM). However the extreme ultraviolet (EUV) spectral index cannot be measured directly for high-redshift quasars owing to absorption at frequencies above the Lyman limit, and as a result, studies of the impact of quasars on the IGM during reionization must assume a spectral energy distribution in the EUV based on observations at low redshift, z≲ 1. In this paper we use regions of high Lyα transmission (near-zones) around the highest redshift quasars to measure the quasar EUV spectral index at z∼ 6. We jointly fit the available observations for variation of near-zone size with both redshift and luminosity, and propose that the observed relation provides evidence for an EUV spectral index that varies with absolute magnitude in the high-redshift quasar sample, becoming softer at higher luminosity. Using a large suite of detailed numerical simulations, we find that the typical value of spectral index for a luminous quasar (M₁₄₅₀=−27) at z∼ 6 is constrained to be α= 1.3^+0.4_−0.3 for a specific luminosity of the form L_ν∝ν^−α. This value is sensitive to the background ionization rate at z∼ 6, and we present an updated estimate of Γ= (0.18^+0.18_−0.09) × 10⁻¹² s⁻¹ as part of our analysis. We find the scatter in spectral index among individual quasars to be in the range Δα∼ 0.75–1.25. These values are in agreement with direct observations at low redshift, and indicate that there has been no significant evolution in the EUV spectral index of quasars over 90 per cent of cosmic time.

galaxies: high-redshift, intergalactic medium, quasars: general, cosmology: theory

1 INTRODUCTION

The discovery of distant quasars has allowed detailed absorption studies of the state of the high-redshift intergalactic medium (IGM) at a time when the Universe was less than a billion years old (Fan et al. 2006a; Willott et al. 2007, 2009). Beyond z∼ 6 several of these quasars show a complete Gunn–Peterson trough in their spectra blueward of the Lyα line (White et al. 2003). However, the spectra of these distant quasars also show enhanced Lyα transmission in the region surrounding the quasar, implying the presence of either an H ii region (Cen & Haiman 2000) or of a proximity zone. The interpretation of these transmission regions has been a matter of some debate. Different arguments in favour of a rapidly evolving IGM at z > 6 are based on the properties of the putative H ii regions inferred around the highest redshift quasars (Mesinger, Haiman & Cen 2004; Wyithe & Loeb 2004; Wyithe, Loeb & Carilli 2005; Kramer & Haiman 2009). On the other hand, Bolton & Haehnelt (2007a), Maselli et al. (2007) and Lidz et al. (2007) have demonstrated that the features in high-redshift quasar spectra blueward of the Lyα line could also be produced by a classical proximity zone. In this case, the spectra provide no evidence for a rapidly evolving IGM. Importantly, these, and other analyses, rely on an assumed value for the extreme ultraviolet (EUV) spectral index to convert from the observed luminosity (redward of rest-frame Lyα) to an ionizing flux, which is the quantity of importance for studies of reionization. As a further example, a value for the EUV spectral index must be assumed to estimate the quasar contribution to the reionization of hydrogen (e.g. Srbinovsky & Wyithe 2007).

The spectral energy distributions of luminous quasars are thought to show little evolution out to high redshift (Fan 2006). For example, broad emission line ratios at z∼ 5 have similar values to those observed at low redshift (Hamann & Ferland 1993; Dietrich et al. 2003). Moreover, optical and infrared spectroscopy of some z∼ 6 quasars has indicated a lack of evolution in the optical–UV spectral properties redward of Lyα (Vanden Berk et al. 2001; Pentericci et al. 2003). At higher energies, the optical/IR-to-X-ray flux ratios (e.g. Brandt et al. 2002; Strateva et al. 2005) and X-ray continuum shapes (Vignali et al. 2003) show at most mild evolution from low redshift. However Jiang et al. (2010) have recently reported hot dust-free quasars at z∼ 6 which have no counterparts in the more local Universe. These quasars are thought to be at an early evolutionary stage, and offer the first evidence that quasar properties have evolved since the reionization era.

As yet, the spectral energy distribution blueward of the Lyman limit, which is critical for studies of reionization, is not measured at z∼ 6 since the quasar spectra are subject to complete absorption. In the future, the EUV spectral index could be measured from 21-cm observations of the thickness of ionizing fronts of the H ii regions around high-redshift quasars (Kramer & Haiman 2008). However in the meantime the assumed value of the spectral index is usually based on direct observations of the EUV spectrum of quasars at z≲ 1, where they are not subject to absorption from the IGM.

In this paper, we show that the extent of Lyα transmission regions around high-redshift quasars is sensitive to the EUV spectral index, and so provides a probe of this otherwise hidden portion of the spectrum. Our paper is presented in the following parts. We first summarize quasar near-zones in Section 2, the observed near-zone relation in Section 3 and our associated numerical modelling in Section 4. We then present constraints on the relation between near-zone size, quasar redshift and quasar luminosity (Section 4). In Sections 6 and 7, we present a discussion of the astrophysical interpretation of our results. As part of our analysis, we provide a revised estimate of the ionizing background at z = 6 and jointly constrain the ionizing background and EUV spectral index. Some additional uncertainties in our modelling are discussed in Section 8, and our conclusions are presented in Section 9. In our numerical examples, we adopt the standard set of cosmological parameters (Komatsu et al. 2009), with values of Ω_b = 0.044, Ω_m = 0.24 and Ω_Λ = 0.76 for the matter, baryon and dark energy fractional density respectively, and h = 0.72 for the dimensionless Hubble constant.

2 NEAR-ZONES

In view of the difficulty of determining the physical conditions of the quasar environment, as well as the observational challenges associated with measuring the Lyα transmission at high redshift, Fan et al. (2006a) defined a specific radius, hereafter referred to as the near-zone radius, at which the normalized Lyα transmission drops to 10 per cent after being smoothed to a resolution of 20 Å. This value of 10 per cent is arbitrary, and is chosen to be larger than the average Gunn–Peterson transmission in z≃ 6 quasar spectra (<4 per cent). In order to make the measurement unambiguous, the near-zone radius is defined at the point where the transmission first drops below this level, even if it rises back above 10 per cent at larger radii. The measured near-zone radius is therefore dependent on the spectral resolution.

Fan et al. (2006a) found a striking relation between near-zone size and redshift, which they quantify using the expression

1

Here the near-zone size, R_cor, has been corrected for luminosity according to the relation

2

where B = 3 and R₂₇ is the mean size of a near-zone at z = 6 around a quasar of M₁₄₅₀=−27. The value of B = 3 (which is roughly consistent with observations; Carilli et al. 2010) is motivated by the evolution of an H ii region (Haiman 2002):

3

where f_{H i} is the fraction of hydrogen that is neutral, is the rate of ionizing photons produced by the quasar and t is the quasar age.¹ The size of the near-zones is found to increase by a factor of 2 over the redshift range 6.4 > z > 5.7, from which Fan et al. (2006a) inferred the neutral fraction to increase by an order of magnitude over that time.

On the other hand, if the near-zone size is set by resonant absorption in an otherwise ionized IGM rather than by the boundary of an H ii region, then its value is instead approximated by the expression

4

where Δ_lim is the normalized baryon density in the transmitting regions as a fraction of the cosmic mean (τ_lim = 2.3, corresponding to the 10 per cent transmission), and T is the IGM temperature (Bolton & Haehnelt 2007b). The parameter α is the EUV spectral index. Equation (4) differs from equation (3) in three important ways. First, if near-zone size is set by resonant absorption then it is not sensitive to the quasar lifetime. Secondly, resonant absorption results in a power-law index describing the dependence of near-zone size on luminosity with a value of B = 2 rather than B = 3. Thirdly, the near-zone size is sensitive to α, suggesting that near-zone size provides an opportunity to measure this unknown parameter.

The precise definition of near-zone radius suggested by Fan et al. (2006a) allows for quantitative comparison with simulations. Wyithe, Bolton & Haehnelt (2008) and Bolton et al. (2010) have modelled near-zones using hydrodynamical simulations with radiative transfer, combined with a semi-analytic model for the evolving, density-dependent ionizing background. The modelling in the latter work is able to produce both the amplitude and evolution of the near-zone sizes, as well as the statistics of absorption lines within the near-zone, and is based on quasars in a highly ionized IGM. In this regime, it is equation (4) with B = 2 rather than equation (3) with B = 3 that is applicable. This indicates that a varying ionizing background could be responsible for the observed relation, providing an alternative scenario to a strongly evolving neutral fraction. Unfortunately, this scenario implies that the evolution of quasar near-zones cannot necessarily be used to infer an evolution in neutral fraction. Thus, although intriguing, it could be argued that the observed near-zone relation has not yet taught us anything concrete about reionization.

When discussing near-zones, most authors have concentrated on the relation between size and redshift, with measured size corrected for luminosity dependence as part of this process. In most cases, these studies used the scaling of near-zone in proportion to luminosity to the one third power (e.g. Fan et al. 2006a; Carilli et al. 2010). As described above, this value is appropriate if the near-zone is expanding into a neutral IGM. On the other hand, a power of one half (B = 2) is more appropriate if the near-zone corresponds to a radius where resonant absorption in an ionized IGM results in the 10 per cent transmission (Bolton & Haehnelt 2007a). Thus the correct scaling between near-zone radius and quasar luminosity is uncertain. In this work, we therefore treat the power-law index B as a free parameter in a two-dimensional relation describing near-zone radius as a function of quasar luminosity and redshift. We find that current near-zone data place strong constraints on the value of B. Moreover, rather than focus on the interpretation of the parameter A with respect to the end of hydrogen reionization, in this paper we concentrate instead on understanding the implications of the parameters B and R₂₇ for the EUV spectral index blueward of the Lyman limit for z∼ 6 quasars. We show that the near-zone relation offers a probe of the EUV spectrum of quasars at these early epochs, information that is otherwise concealed by the Lyman-limit absorption, and so has previously only been studied for quasars at redshifts of z≲ 1 (Telfer et al. 2002).

3 OBSERVED NEAR-ZONE SAMPLE

Carilli et al. (2010) have assembled a sample of 25 z∼ 6 quasars which have quality rest-frame UV spectra and redshift measurements (Fan et al. 2006a; Jiang et al. 2008). Eight of the sample are detected (Wang et al. 2010) in CO. Another nine have redshifts measured (Jiang et al. 2007; Kurk et al. 2007) from the Mg ii line emission. For the other eight objects, Carilli et al. (2010) adopt redshifts from the relevant discovery papers, which are mainly determined with the Lyα+NV lines (Fan et al. 2006b). The measurements are summarized in table 1 of Carilli et al. (2010), including the quasar absolute AB magnitude at 1450 Å (M₁₄₅₀), quasar redshifts (z) and near-zone size (R_p). Three sources in the original sample (Wang et al. 2010) listed in table 1 of Carilli et al. (2010), J0353+0104, J1044−0125, and J1048+4637, are broad absorption line quasars (Fan et al. 2006b; Jiang et al. 2008), while the source J1335+3533 has a lineless (Fan et al. 2006b) UV spectrum. We exclude these quasars from our analysis owing to the fundamentally different nature of their intrinsic spectra (Fan et al. 2006a; Carilli et al. 2010). Thus the sample consists of 21 z∼ 6 quasars which we use to analyse the near-zone relation. Near-zone sizes are plotted in Fig. 1 as a function of both absolute magnitude and redshift. Typical measurement errors (Carilli et al. 2010) are 1.2, 0.4 and 0.1 Mpc in the estimated near-zone size (R_p) introduced by UV, Mg ii and CO determined redshift uncertainties, respectively. There are clear trends with both luminosity (Carilli et al. 2010) and redshift (Fan et al. 2006a; Carilli et al. 2010).

Figure 1

Measured sizes R_p (not scaled for quasar luminosity) plotted against M₁₄₅₀ (left) and z (right), respectively, for observed (black dots) and simulated (triangles) near-zones. Measured near-zone sizes are also plotted for CFHQS J1509−1749, CFFQS J2329−0301 and CFHQS J0210−0456 (squares). The simulated data have been computed at the best fit in α₀ and α₁ (see text for details).

We note that quasar near-zone sizes have also been measured (Willott et al. 2007, 2010) for CFHQS J1509−1749, CFFQS J2329−0301 and CFHQS J0210−0456. Following Carilli et al. (2010) we do not include these near-zones, which are drawn from a different data set, in our analysis. However as may be seen from Fig. 1, these near-zones lie on the (Carilli et al. 2010) correlation, and we have checked that their addition does not significantly alter our results.

4 MODELLING OF NEAR-ZONES

Our near-zone simulations combine a semi-analytical model for the evolving, density-dependent photoionization rate in the biased regions surrounding quasars, with a radiative transfer implementation, and realistic density distributions drawn from a high-resolution cosmological hydrodynamical simulation. These simulations are discussed in detail elsewhere (Bolton et al. 2010), and their description is not reproduced here. However for completeness we show the evolution of the background photoionization rate for our fiducial simulations in the left-hand panel of Fig. 2. The grey curve shows the photoionization rate in the mean IGM, while the solid, dotted and dashed black curves show the mean and 1σ range of the photoionization rate in the biased regions within 5 Mpc of a 10¹³ M_⊙ halo. Importantly, we construct the absorption spectra using an ionizing background computed as a function of proper time along the trajectory of a photon emitted by the quasar, rather than at the proper time of the quasar. This effect is appropriate if considering spectra at the end of the reionization era, when the ionizing background can evolve significantly during the light travel time across a quasar near-zone. The right-hand panel of Fig. 2 illustrates the background photoionization rate as a function of radius in our model (for a quasar observed at z = 6).

Figure 2

The fiducial model for the ionizing background. Left: the evolution of the ionizing background (Wyithe et al. 2008). The grey curve shows the background in the mean IGM, while the solid, dotted and dashed black curves show the mean and 1σ range of ionizing background in the biased regions within 5 Mpc of a 10¹³ M_⊙ quasar host halo. The data points show our estimate of the ionizing background at z = 5 and 6 (error bars represent the 68 per cent range). The point at z = 4 has been reproduced from Bolton & Haehnelt (2007a). Right: the ionizing background as a function of radius in our fiducial model (for a quasar observed at z = 6). The ionizing background has been computed as a function of proper time along the trajectory of a photon emitted by the quasar, rather than at the proper time of the quasar. This results in the slow rise of the background level at large radii.

For each of the model near-zone sets constructed, the sizes of the simulated near-zones were calculated for a sample of synthetic quasar spectra with the same redshifts and quasar absolute magnitudes as the Carilli et al. (2010) compilation. The synthetic spectra are designed to match the spectral resolution (R∼ 2500) and signal-to-noise ratio (S/N) of the observational data (⁠ per 3.5 Å pixel⁻¹) as closely as possible. For the fiducial simulations, an EUV spectral index of α = 1.5 (where dL/dν∝ν^−α) was assumed to calculate the corresponding ionizing photon rate (⁠⁠) from the observed M₁₄₅₀. Additionally, a temperature at mean density of T = 23600⁺⁵⁰⁰⁰₋₆₉₀₀ K has been measured (Bolton et al. 2010) within the near-zone around SDSS J0818+1722. Our fiducial simulations are consistent with this constraint, and we further assume this as the best-fitting temperature when modelling near-zones throughout our analysis. This level of complexity is important as the near-zone size (Fan et al. 2006a) is measured based on the first resolution element where the transmission in the spectrum drops below 10 per cent, meaning that the Lyα forest must be simulated in detail in order for the modelling to be considered reliable. Examples of some of the synthetic quasar spectra used in our analysis are displayed in Fig. 3.

Figure 3

Example synthetic near-zone spectra. In both panels, the dotted curve shows the spectrum after being smoothed to a resolution of 20 Å in the observed frame, while the vertical dashed line marks the measured near-zone size. This corresponds to the position where the smoothed flux first drops below a transmission of 10 per cent. Top: quasar at z = 5.8918 with M₁₄₅₀=−26.09. Bottom: quasar at z = 6.247 with M₁₄₅₀=−26.67.

5 NEAR-ZONE RELATION CONSTRAINTS

As noted in the introduction, the appropriate value of B describing the power-law relation between near-zone size and luminosity is uncertain. Therefore, rather than correct for luminosity by assuming a fixed B and fit equation (1) to determine the evolution of near-zone size with redshift as has been done previously, in our paper we instead model the observed near-zone plane using free parameters to describe both the evolution with redshift (A) and the evolution with absolute magnitude (B),

5

In this section, we use equation (5) to quantify constraints from the near-zone plane at z∼ 6.

We find the parameter sets that describe the near-zone plane by fitting equation (5) to both the observed near-zone relation and to our simulated data sets, and calculating the resulting values of χ². As part of this process, we introduce intrinsic scatter into equation (5) so that the reduced χ² of the best-fitting model is unity. We find a smaller intrinsic scatter (σ_sim = 1.27 ± 0.04 Mpc) in the simulated data, indicating that sightline variation is not responsible for all of the observed scatter (σ_obs = 1.86 ± 0.06). The black and grey contours in Fig. 4 show the resulting projections of Δχ²=χ²−χ²_min, where χ_min is the minimum value, on to the parameter spaces (A, B) and (A, R₂₇) for the observed and fiducial simulated near-zones, respectively. We note that our model has three parameters, but that we are interested in the constraints on these individually. The natural confidence regions of choice are the projections of the three-dimensional regions defined by fixed Δχ² on to the one-dimensional spaces of interest. These projected values of Δχ² are distributed as a chi-squared distribution with 1 degree of freedom (e.g. Press et al. 1992). The projection of the contours in Fig. 4 (as well as latter figures) on to the axis representing a particular parameter therefore describes confidence intervals of 68, 84, 91 and 97 per cent, respectively.

Figure 4

Projections of χ² for observed (black) and simulated (grey) near-zone sets assuming the fiducial evolving model. We have introduced intrinsic scatter (σ_obs = 1.86 Mpc and σ_sim = 1.27 Mpc, respectively) so that the reduced χ² is unity for the best-fitting model. The contours represent Δχ² = 1, 2, 3 and 5 relative to the corresponding minimum value of χ². The projection of these contours on to the axis representing a particular parameter describe confidence intervals of 68, 84, 91 and 97 per cent, respectively. The best-fitting values for the observed relation (shown as large dots) are A = 8.0 ± 2.5, B = 3.3^+1.7_−0.8 and R₂₇ = 7.75 ± 0.35, while the best-fitting values describing the simulations are A = 9.8 ± 2.0, B = 1.9 ± 0.3 and R₂₇ = 8.45 ± 0.35.

While the simulations reproduce the size evolution with redshift, the luminosity dependence is not well described, with B∼ 2 for the simulations (equation 4) but B∼ 3 for the observed sample. There are two explanations for this disagreement between the data and fiducial model which we now discuss in turn.

6 NEAR-ZONE RELATION DETERMINED BY H ii REGIONS

First, the observed value of B = 3 is expected if the quasars were surrounded by H ii regions embedded in a neutral IGM (Cen & Haiman 2000). Naively, this could be interpreted as evidence for a significantly neutral IGM. To test this, we calculate the expected scatter in the near-zone relation for comparison with observation. In contrast to the case for a proximity zone in an optically thin IGM, the size of an H ii region depends on the age of the quasar at the time of observation. As a result, there is scatter in the observed size even for fixed redshift and quasar luminosity. We can estimate the scatter due to the age of the quasars by noting that if there is equal probability of observing a quasar at any age during its lifetime, then the probability distribution for the H ii region radius at fixed luminosity L, hydrogen neutral fraction and redshift z, given a quasar lifetime t_max, is

6

where R_max is the size reached after t_max. Given this distribution, the mean is 〈R_p〉 = 3/4 R_max, and the variance is σ²_R≡〈R²_p〉−〈R_p〉²=〈R_p〉²/15. Observationally, since we have 〈R_p〉∼R₂₇∼ 8 Mpc, the scatter in H ii region radius owing to random quasar age is σ_R∼ 2 Mpc. Thus the scatter in observed quasar age accounts for all of the observed scatter (σ_obs = 1.86 ± 0.06) in the near-zone sizes. However in the case of H ii regions, we would also expect additional scatter in the size owing to other quantities like the quasar lifetime and the spectral index (we return to this point below), as well as scatter due to inhomogeneities in the density field and ionization structure along the different quasar sightlines. In the case of the latter, the sizes of H ii regions near the end of reionization (Furlanetto, Zaldarriaga & Hernquist 2004) are thought to vary over the range 1–2 proper Mpc, which would lead to a component of scatter in addition to the ∼2 Mpc expected from the quasar age. Thus, the total scatter in an H ii region generated near-zone relation would be in excess of σ_R∼ 3 Mpc. We therefore infer that although quasar H ii regions would lead to a near-zone relation with the appropriate value of B∼ 3, the observed relation between near-zone size and magnitude would be too tight in this case. This scatter-based constraint would be alleviated if the near-zone sizes corresponded to resonant absorption within an H ii region (the scatter would be smaller since the size is independent of lifetime in this case). However, in this alternative case where resonant absorption sets the near-zone size, our modelling suggests that B∼ 2 rather than B∼ 3 should be observed. As a result, H ii regions cannot provide a viable physical explanation for the observational data (see also Maselli, Ferrara & Gallerani 2009).

7 CONSTRAINING THE EUV SPECTRAL INDEX

An alternative scenario to explain the observed relation between near-zone size and magnitude is provided by an EUV spectral index which is luminosity-dependent, so that at brighter absolute magnitudes the ionizing flux is smaller. To quantify this point, we first note that the observed quantities are the near-zone size and the absolute UV magnitude. Given this UV magnitude, the ionization properties of the surrounding medium are dependent on the ionizing luminosity in photons per second (⁠⁠), and the spectral index blueward of the Lyman break. The relation between and α is

7

where is the luminosity at the Lyman limit ν_Ly, h is Planck's constant and x_max is the ratio between the frequency where photons cease to contribute to ionization and the Lyman limit (we take x_max = 100, though this choice does not affect our results). To estimate the value of for the quasars in the Carilli et al. (2010) sample, the observed M₁₄₅₀ has previously (Bolton et al. 2010) been combined with an assumed spectral index of α = 1.5. Our approach to investigate whether a variable spectral index is able to explain the discrepancy in B is therefore to modify the value of absolute magnitude that corresponds to the ionization rate at the near-zone radius measured in the simulations. This allows us to investigate the near-zone relations that result from a range of spectral indices without recalculating the radiative transfer simulations. Beginning with equations (4) and (7) and assuming a spectrum of the form (e.g. Vanden Berk et al. 2001; Telfer et al. 2002; Pentericci et al. 2003):

8

we adjust the absolute magnitude corresponding to the model quasars (at fixed near-zone radius) for a variation in spectral index using

9

where M_1450,1.5 is the absolute magnitude used in the fiducial model. Here we note that simulations performed with α = 0.5, 1.0, 1.5, 2.0 and 2.5 yield near-zone temperatures of T∼ 2.81 × 10⁴, ∼2.34 × 10⁴, ∼2.06 × 10⁴, ∼1.89 × 10⁴ and ∼1.76 × 10⁴ K, respectively. These values are fit (to within 1 per cent) by

10

This dependence links the temperature and spectral index in equation (4), and contributes the last product in equation (9).

To quantify the possible dependence of spectral index on luminosity, we parametrize the value of α as a function of M₁₄₅₀ using

11

where α₀ represents the mean spectral slope of quasars with M₁₄₅₀=−27 and α₁ = dα/dM₁₄₅₀ describes the variation of spectral index with magnitude. Note that since α is a function of M₁₄₅₀, we evaluate it not at the observed M₁₄₅₀, but instead use an absolute magnitude M′₁₄₅₀ that is adjusted to account for the average index α₀ at M₁₄₅₀=−27,

12

We use these equations to constrain the likelihood of parameter combinations (α₀, α₁) through convolution of the resulting joint probability distribution for the parameters A, B and R₂₇ that describe simulated near-zone relations, with the corresponding relation for the observations:

where

and

13

Here we have assumed flat prior probability distributions for each of the parameters A, B and R₂₇. As part of this procedure, we calculate the intrinsic scatter for each set of parameters (α₀, α₁), and impose the condition that this scatter be smaller than the observed scatter by weighting the calculated likelihood with

14

Here the value of 0.06 dex represents the uncertainty in σ_obs.

Contours of likelihood for the values of α₀ and α₁ based on our fiducial model are plotted in the lower left panel of Fig. 5. We find that the preferred value of α₀ = 1.9 ± 0.3 is greater than the assumed α = 1.5. In addition, since the fiducial simulations underpredict B (Fig. 4), the value of α₁ is found to be greater than zero with high confidence, so that more luminous quasars have softer spectra. This is consistent with the results of Strateva et al. (2005),who found that brighter quasars are relatively less luminous in the X-ray, and is also expected from theoretical modelling of the quasar spectrum (Wandel & Petrosian 1988).

Figure 5

Upper panels: projections of χ² for observed (black) and simulated (grey) near-zone sets for the best-fitting simulations (in α₀ and α₁) assuming the fiducial evolving ionizing background. We have introduced intrinsic scatter so that the reduced χ² is unity for the best-fitting model. The contours represent Δχ² = 1, 2, 3 and 5 relative to the corresponding minimum value of χ². The best-fitting values (shown as large dots) describing the simulations are A = 9.9 ± 1.8, B = 3.1^+0.4_−0.3 and R₂₇ = 7.6 ± 0.2, which can be compared with the values for the observed relation of A = 8.0 ± 2.5, B = 3.3^+1.7_−0.8 and R₂₇ = 7.75 ± 0.35. Lower left panel: contours of likelihood for the values of α₀ and α₁ (at 61 and 14 per cent of the peak likelihood). The values of 61 per cent and 14 per cent are chosen to equal the contour height corresponding to Δχ² = 1 and Δχ² = 4, and therefore the projection of these contours on to the axis for a particular parameter represents the 1σ and 2σ ranges, respectively. The most likely values (shown as a large dot) are α₀ = 1.9 ± 0.3 and α₁ = 1.25^+1.50_−0.75. Lower right panel: contours of scatter in spectral index α (thick solid and thin dotted). The point with error bars shows the corresponding parameters measured (Telfer et al. 2002) at low redshift. For comparison, the dashed contour repeats the 60 per cent contour of likelihood for the values of α₀ and α₁.

The corresponding contours (grey curves) of Δχ² in the parameter spaces (A, B) and (A, R₂₇) for the simulated near-zone relations computed using the most likely values of α₀ and α₁ are shown in the upper panels of Fig. 5, along with contours for the observed near-zone sample (black curves). The best-fitting simulations reproduce the observed values of average size R₂₇ and index B. The parameter A, which is sensitive to the evolution of the ionizing background, is also consistent with observations (Fan et al. 2006a; Bolton & Haehnelt 2007b). The near-zone sizes for the best-fitting model are plotted in Fig. 1 (triangles).

7.1 Intrinsic scatter in the spectral index

As mentioned previously in Section 5, the intrinsic scatter in the simulated near-zones (σ_sim = 1.25 for the best-fitting model) is smaller than the scatter σ_obs = 1.86 in the observed near-zone relation. Since the numerical simulations already include line-of-sight density variations in the IGM, the missing scatter,

15

could be provided by scatter in the spectral index among individual quasars, which is known to be significant at lower redshift (Telfer et al. 2002). Using the relation for the missing scatter in R_p, we obtain the magnitude of scatter in the spectral index that is necessary to simulate the observed near-zone scatter:

16

Contours of this scatter are shown in the right-hand panel of Fig. 5. Values for scatter in spectral index ranging between Δα∼ 0.75 and ∼1.25 correspond to the range of parameters α₀ and α₁, which lead to models for the near-zone relation that are in good agreement with the data.

7.2 Constant ionizing background models

Our simulations of quasar near-zones assume a semi-analytic model for the ionizing background (see Wyithe et al. 2008). It is therefore important to ask whether the results in this paper pertaining to the dependence of near-zone size on luminosity (i.e. the constraints on B) and the resulting conclusions regarding the quasar spectral index α are sensitive to this model. To address this issue, we have repeated our modelling of quasar near-zones assuming an ionizing background that is independent of redshift. We first take the z = 6 value from our fiducial model; Fig. 6 shows contours describing the resulting constraints (grey curves). The results are very similar to the fiducial case, with the exception of A, which shows that the redshift evolution of near-zone size in the constant background model is smaller than observed. This is consistent with the previous inference that the trend of near-zone size with redshift is being driven by the rising intensity of the ionizing background (Wyithe et al. 2008) at z∼ 6. We have also tested sensitivity to the ionizing background amplitude by repeating our analysis for 16 different values ranging between 1/100 and times the fiducial model. In the left-hand panel of Fig. 7, we show the resulting range for α₀ as a function of the background photoionization rate ⁠. This figure illustrates that results for α depend on the assumed value of Γ₁₂. As a result, we next constrain Γ₁₂ following the work of Bolton & Haehnelt (2007b).

Figure 6

Upper panels: projections of χ² for observed (black) and simulated (grey) near-zone sets for the best-fitting simulations (in α₀ and α₁) assuming the non-evolving ionizing background. We have introduced intrinsic scatter so that the reduced χ² is unity for the best-fitting model. The contours represent Δχ² = 1, 2, 3 and 5 relative to the corresponding minimum value for χ². The best-fitting values (shown as large dots) describing the simulations are A = 4.75 ± 2.25, B = 3.2^+0.9_−0.5 and R₂₇ = 7.8 ± 0.3, which can be compared with the values for the observed relation of A = 8.0 ± 2.5, B = 3.3^+1.7_−0.8 and R₂₇ = 7.75 ± 0.35. Lower left panel: contours of likelihood for the values of α₀ and α₁ (at 61 and 14 per cent of the peak likelihood). The most likely values (shown as a large dot) are α₀ = 1.9 ± 0.4 and α₁ = 1.0^+1.8_−0.9. Lower right panel: contours of scatter in spectral index α. The point with error bars shows the corresponding parameters measured (Telfer et al. 2002) at low redshift.

Figure 7

Left-hand panel: the best-fitting value for α₀ as a function of Γ₁₂ (solid curve), together with the 68 per cent range (dashed curves). Central panel: the joint probability distribution for Γ₁₂ and α₀ (the most likely value is shown as a large dot). The prior probability for Γ₁₂ is assumed to be flat in the logarithm (dP_prior/dΓ₁₂∝ 1/Γ₁₂). Right-hand panel: the joint probability distribution for α₀ and α₁ (the most likely value is shown as a large dot). In the central and right-hand panels, the contours represent levels that are at 61 and 14 per cent of the peak likelihood. The final estimates for the parameters describing the EUV spectral index (taken from the projection of probabilities on to a single parameter axis) are α₀ = 1.3^+0.4_−0.3, with α₁ = 0.9^+1.3_−0.7.

7.2.1 The ionizing background at z = 6

Fan et al. (2006a) present eight values of transmission (with uncertainty ⁠) measured in redshift intervals of Δz = 0.1 centred on redshifts in the range 5.9 ≤z≤ 6.25, and 14 values centred in the redshift range 4.9 ≤z≤ 5.15. These values are listed in Table 1. From these values we calculate the mean value of transmission and ⁠, respectively. To estimate this uncertainty, we have calculated the standard error on the mean transmission via a bootstrap analysis which includes the uncertainties on individual parameters. The values of and are consistent with the values quoted previously (Bolton & Haehnelt 2007a) at z = 6 and 5, respectively, although note that this study uses slightly wider redshift bins and does not use the independent data of Songaila (2004). However our estimates of the uncertainty are smaller, particularly at z = 6 where the previously quoted error (Bolton & Haehnelt 2007a) is a factor of ∼2 larger. This difference can be attributed to the fact that rather than estimate the standard error on the mean transmission, Bolton & Haehnelt (2007a) used an interquartile range as an estimate for the error on the transmission values along various lines of sight.

Table 1

Table of transmissions and uncertainty in Δz = 0.1 regions from (Fan et al. 2006a).

z
4.92	0.1276	0.0011
4.95	0.1139	0.0022
4.98	0.1002	0.0050
5.01	0.1268	0.0117
5.04	0.1567	0.0073
5.06	0.1765	0.0190
5.06	0.1285	0.0042
5.06	0.1509	0.0224
5.07	0.0751	0.0011
5.08	0.0530	0.0025
5.10	0.0898	0.0020
5.11	0.1650	0.0030
5.13	0.1293	0.0064
5.13	0.1243	0.0050
5.90	0.0108	0.0033
5.93	0.0125	0.0022
5.95	0.0038	0.0005
5.95	0.0060	0.0010
6.08	−0.0071	0.0020
6.10	0.0012	0.0010
6.10	0.0051	0.0005
6.25	0.0015	0.0005

To estimate the corresponding background photoionization rate, we next combine previously published (Bolton & Haehnelt 2007a) data with a Monte Carlo method. We first determine the resulting distribution of effective optical depths from our transmission estimates, finding values of τ_eff = 5.5^+0.6_−0.4 and 2.1^+0.08_−0.07 at z = 6 and 5, respectively. Here and below, we have quoted the median and the 68 per cent range. We note that the value of τ_eff is not defined in the small number of cases where the Monte Carlo trial value of is negative. In such cases, we set and note that our results are insensitive to the exact value, provided it is well below the observational limits in the Gunn–Peterson troughs. Given the distribution of effective optical depths, we determine the corresponding probability distribution for the normalized photoionization rate, ⁠. To estimate the uncertainties, we have utilized scaling relations derived from hydrodynamical simulations (see table 4 of Bolton & Haehnelt 2007a) for Γ₁₂, as functions of τ_eff, the IGM temperature [which we take to be log ₁₀(T/10⁴ K) = 0 ± 0.3], and the slope of the temperature–density relation. We also utilize scaling relations as a function of the cosmological parameters.

At z = 5 we find Γ₁₂ = 0.47^+0.3_−0.2, while at z = 6 we obtain Γ₁₂ = 0.18^+0.18_−0.09. We note that the smaller uncertainty in τ_eff relative to the previous estimate (Bolton & Haehnelt 2007a) allows us to quote a probability density for the measured ionizing background (with a non-zero value preferred at 2σ) rather than an upper limit. The revised estimates of the ionizing background at z = 5 and 6 are presented in Fig. 2.

7.2.2 Joint probability forΓ₁₂andα₀

We use the constraints for Γ₁₂, together with the likelihoods for α₀ and α₁ as a function of Γ₁₂, to generate a joint probability distribution for Γ₁₂ and α₀ (central panel of Fig. 7). Similarly, we obtain the joint probability for α₀ and α₁ (right-hand panel of Fig. 7). This procedure yields an estimate of α₀ = 1.3^+0.4_−0.3, with α₁ = 0.9^+1.3_−0.7. These constraints represent the primary measurement of our paper.

7.3 Comparison with spectral index at low redshift

Telfer et al. (2002) measured the UV spectral index for quasars at z∼ 1–2. By analysing the composite spectrum, they found a mean value of 〈α〉 = 1.76 ± 0.12. (Telfer et al. 2002) also measured α for a subsample of individual quasars, allowing them to estimate the variation α₁ = 0.3 ± 0.1 with luminosity,² as well as the intrinsic scatter Δα = 0.76. These estimates are plotted alongside our model constraints in Figs 5–8. We find that the parameter combination α₀ = 1.76 ± 0.12 and α₁ = 0.3 ± 0.1 is consistent with our inference for these parameters at z∼ 6. Moreover, we find that our inferred scatter in α at z∼ 6 agrees well with the Δα = 0.76 found by (Telfer et al. 2002). Thus, assuming our fiducial model of the ionizing background, there is no evidence for evolution in the EUV slope of quasars with redshift out to z∼ 6.

Figure 8

Figure showing the systematic uncertainties on the joint probability distribution for α₀ and α₁. In each case the contours represent levels that are at 61 and 14 per cent of the peak likelihood respectively, and the large dots represent the most likely values. Left-hand panel: the dark contours show constraints that have been computed based on a model with α = 0.5. The best-fitting values for this model are α₀ = 2.4^+0.6_−0.4 and α₁ = 1.0^+2.7_−0.9. Central panel: the dark contours show the constraints in a case where the continuum has been deliberately placed at a level that is 5 per cent low. The best-fitting values for this model are α₀ = 2.15^+0.6_−0.4 and α₁ = 0.6^+1.8_−0.8. Right-hand panel: the dark contours show the constraints in a case where self-shielded absorbers have been added to regions with density contrast Δ > 16. The best-fitting values for this model are α₀ = 1.6 ± 0.4 and α₁ = 0.6^+1.3_−0.7. In each case, the model assumes an unevolving ionizing background, for which the light contours in each panel show the constraints. The best-fitting values for the fiducial model are α₀ = 1.9 ± 0.4 and α₁ = 1.0^+1.8_−0.9. The overlap between the projection of the 61 per cent contours and the best-fitting values for our fiducial and ‘modified’ models shows that the differences between the constraints on α₀ and α₁ are less than the 1σ range for each effect considered.

8 UNCERTAINTIES

In addition to those already discussed, we have checked several other sources of uncertainty. We have recalculated our constraints including the assumption of a fiducial model with α = 0.5 rather than α = 1.5, a model where the quasar continuum has been placed systematically high by 5 per cent on the data, and a model that includes an approximate treatment of Lyman limit systems. Each of these effects is found to have a relatively small effect on the derived parameters for the spectral index relative to the quoted uncertainty, as detailed below. Specifically, our results show that for each effect considered, the difference between the best-fitting value for our fiducial model and the best-fitting value for the ‘modified’ model is less than the 1σ range both in α₀ and α₁.

8.1 Near-zone temperature

There is an uncertainty in the absolute temperature of the IGM within the quasar near-zone, as well as scatter in the temperature from one near-zone to another. These uncertainties will impact on our constraints for the EUV spectral index. First, any scatter in temperature among different near-zones will contribute to scatter in the observed near-zone sizes. However this scatter will be degenerate with scatter in the spectral index (Δα). Thus the potential presence of scatter in temperature does not introduce uncertainty into our estimated mean value of α₀ or α₁. However, it does imply that our quoted Δα should be considered an upper limit.

On the other hand, uncertainty in the absolute value of the temperature within the near-zones could introduce a more serious error. The temperature of the IGM in the quasar near-zone has two independent contributions which add linearly, and hence the near-zone temperature has two components of uncertainty. First, the mean IGM has an uncertainty in temperature resulting from heating during hydrogen reionization. This component of uncertainty in the temperature due to hydrogen reionization is already included (Bolton & Haehnelt 2007a) in the calculation of the ionizing background.

However in addition the reionization of He ii by the quasar may heat the IGM in the near-zone to well above this mean IGM temperature (Bolton, Oh & Furlanetto 2009). In our modelling, we have assumed a temperature (Bolton et al. 2010) that is consistent with observations of a single z∼ 6 quasar (SDSS J0818+1722), and that the quasar is responsible for reionization of He ii within the near-zone. This latter assumption means the temperature increase from He ii reionization by the quasar is directly related to the spectral hardness of the quasar, and hence to α, which is the quantity we are trying to measure. As a result, the assumption of an incorrect temperature could bias the inferred value of α in our analysis. A hardening of the spectral index leads both to higher temperatures within the near-zone and to an enhanced photoionization rate. By assuming a fixed temperature for all modelling, our analysis would include only the latter enhancement, whereas both ionization and temperature influence the optical depth. We have therefore accounted for this effect via the scaling in equation (9) based on a calibration of near-zone temperature with spectral index estimated from additional radiative transfer simulations. Furthermore, to check the importance of the temperature effects, and at the same time check that our assumed fiducial value of α does not influence the results, we have repeated the constraints on α₀ and α₁ based on the fiducial unevolving ionizing background using a simulated quasar with a value of α = 0.5 (and correspondingly hotter temperature). As shown in the left-hand panel of Fig. 8, we find that the assumed fiducial α (and corresponding fiducial temperature) does not strongly influence the constraints.

8.2 Continuum placement

In addition to errors in the redshifts of the quasars, uncertainty in continuum placement may also lead to an error in the near-zone radius R_p. This error corresponds to an uncertainty (Fan et al. 2006a) of Δτ∼ 0.05, or equivalently to a relative uncertainty in the transmission of per cent. At the near-zone radius per cent, so that the uncertainty in transmission is ⁠. Inspection of fig. 3 in Carilli et al. (2010) shows that when averaged over the near-zone sample, the gradient of transmission with radius is ⁠. Hence the uncertainty in the near-zone radius that is introduced by the uncertainty in continuum placement is Mpc. Thus, continuum placement is subdominant in terms of contribution to the uncertainty in near-zone size. We note that the uncertainty could be larger for individual quasars. However, this would introduce additional scatter into the near-zone relation which would be degenerate with the other sources of scatter. To check this conclusion, we have again repeated the constraints on α₀ and α₁ based on an unevolving ionizing background using simulated quasars for which the continuum has been systematically placed 5 per cent low. As shown in the central panel of Fig. 8, we find that the continuum placement has a systematic effect on near-zone size and hence α₀. However the shift is much smaller than the statistical uncertainty. We thus conclude that our results are robust to a systematic uncertainty in the continuum at this level.

8.3 Lyman limit systems

Finally, a potential shortcoming of our modelling is the assumption that the IGM is optically thin, which results in the absence of absorbers that are self-shielded with respect to the ionizing background. To address this, we have tested the sensitivity of near-zone size by adding a population of self-shielded Lyman-limit systems. We achieve this by setting all gas with density contrast Δ > 16 (corresponding to the density contrast of a clump that will self-shield if immersed in our fiducial ionizing background; Furlanetto & Oh 2005) to be neutral prior to the quasar turning on (Bolton & Haehnelt 2007a). We find that their inclusion does not alter the predicted distribution of near-zone sizes in most cases, and that as shown in Fig. 8, the presence of self-shielded systems should not significantly modify our measurement of α₀ relative to other uncertainties.

A related question concerns an alternative explanation for the observed B = 3 due to an enhanced number density of absorbing systems in the biased environments surrounding luminous quasars. Such an enhancement could preferentially retard the observed near-zone sizes, leading to an increase in near-zone size with luminosity that is slower than expected in the optically thin regime. However we have previously (Bolton & Haehnelt 2007a) investigated the dependence of near-zone size on quasar host-halo mass, finding a negligible effect. This is because the absorbers that set the near-zone size are located between 30 and 100 comoving Mpc from the quasar. Thus, the environmental dependence of absorber abundance is unlikely to explain the observed B = 3.

9 CONCLUSION

The discovery of a strong correlation of the sizes of near-zones around high-redshift quasars and their redshift in the range 5.7 ≲z≲ 6.4 (Fan et al. 2006a; Carilli et al. 2010) has prompted a series of studies aimed at understanding its implications for the tail end of reionization (Fan et al. 2006a; Bolton & Haehnelt 2007a; Wyithe et al. 2008; Bolton et al. 2010; Carilli et al. 2010). In this paper, we have noted that the amplitude of the near-zone radius relation is sensitive to the assumed EUV spectral index of the quasar. We therefore performed a large suite of numerical simulations of the near-zone spectra using radiative transfer through a hydrodynamical model of the IGM in the quasar environment, and find that the observed near-zone relation can be used to constrain the UV spectral index of z∼ 6 quasars. We find that the typical value of spectral index for luminous (M₁₄₅₀=−27) quasars at z∼ 6 is α = 1.3^+0.4_−0.3 (for a specific luminosity of the form L_ν∝ν^−α), where the uncertainty includes both the statistical uncertainty in the near-zone relation as well as uncertainty in the ionizing background. Based on comparison of the scatter in our simulations and the scatter in the observed near-zone relation, we infer that the scatter in the spectral index is Δα∼ 0.75–1.25.

As part of our analysis, we present an updated estimate of the background ionization rate at z∼ 6 of Γ = (0.18^+0.18_−0.09) × 10⁻¹² s⁻¹. While the derived value α depends on the ionizing background at z = 6, we find that it is quite insensitive to the evolution in Γ₁₂ with redshift. We have also investigated a series of other potential systematic uncertainties, including the value of α assumed in the simple model, the level of continuum placement in the data and the possible effect of Lyman limit systems. Each of these effects is found to have a relatively small effect on the derived parameters for the spectral index relative to the quoted uncertainty.

(Carilli et al. 2010) noted that the data are consistent with near-zone sizes that scale with quasar luminosity as R_p∝L^1/B, where B = 3, and that the scatter in the near-zone radius–redshift relation is reduced if this scaling is applied. However, as pointed out by Bolton & Haehnelt (2007a), the value of the power-law index B describing the scaling of near-zone radius with luminosity that is physically appropriate depends on whether near-zone sizes are set by H ii region boundaries (B = 3) or resonant absorption in an ionized IGM (B = 2). Motivated by the possibility of using this second parameter in the near-zone relation to gain a better understanding of its meaning, in this paper we jointly fit the available observations for variation of near-zone size with both redshift and luminosity. From this analysis we find that the observations have an index of B∼ 3, with an uncertainty that excludes a value of B = 2 with high confidence. In contrast to the observations, we find that our fiducial numerical simulations predict a value of B = 2.

We have discussed two possibilities for this difference. First our simulations are conducted in a highly ionized IGM. It is possible that although these simulations are able to describe the size and evolution of the near-zones (Bolton et al. 2010), they do not represent physical conditions of the IGM at z∼ 6, which could instead contain hydrogen that is significantly neutral as previously argued in Wyithe & Loeb (2004) and Wyithe et al. (2005). Although this would lead to the observed B = 3, we do not believe that this is the case because the observed relation is too tight to be explained by H ii regions, the sizes of which contain a range of additional sources of scatter. Instead, we argue that the observed value of B∼ 3 is evidence for an UV spectral index that varies with absolute magnitude in the high-redshift quasar sample, becoming softer at higher luminosity. This finding is in agreement with theoretical modelling of quasar spectra (Wandel & Petrosian 1988).

The results of this paper provide the first constraints on the properties of the EUV spectral index of the highest redshift quasars, which are not observable owing to absorption blueward of the Lyman limit. The values we find are in agreement with direct observations at low redshift, and therefore indicate that there is no evidence for evolution in the EUV properties of quasars over most of cosmic time. The new constraints will aid all studies of reionization involving high-redshift quasars, which have previously relied on an assumed value of the mean EUV spectral index measured at lower redshift.

1

This calculation assumes clumping and recombinations to be unimportant, that the density is at the cosmic mean over the large scales being considered, and that the quasar lifetime is much less than the age of the Universe.

2

Note that the evolution of α with luminosity is dominated by radio loud quasars, which formed half the total sample of quasars for which individual measurements of α were made by (Telfer et al. 2002).

This work was supported in part by the Australian Research Council. We thank Chris Carilli and Martin Haehnelt for helpful comments.