UV Background Fluctuations and Three-Point Correlations in the Large Scale Clustering of the Lyman-alpha Forest

TL;DR

This paper predicts the three-point correlation function of the Lyman-alpha forest at redshift 2.3 using simulations, exploring the impact of UV background fluctuations and assessing the detectability of the signal in current and future surveys.

Contribution

It provides the first theoretical predictions for the 3PCF of the Lyα forest considering UV background fluctuations and evaluates its observational detectability.

Findings

The 3PCF shows hierarchical scaling with the square of the 2PCF.

UVB fluctuations can either suppress or boost the 2PCF and 3PCF depending on the mean free path.

Predicted signal-to-noise ratios suggest the 3PCF is detectable in current and upcoming surveys.

Abstract

Using the Ly$\alpha$ mass assignment scheme (LyMAS), we make theoretical predictions for the 3-dimensional 3-point correlation function (3PCF) of the Ly$\alpha$ forest at redshift $z=2.3$. We bootstrap results from the (100 $h^{-1} \mbox{ Mpc}$)$^3$ Horizon hydrodynamic simulation to a (1 $h^{-1}$ Gpc)$^3$ $N$-body simulation, considering both a uniform UV background (UVB) and a fluctuating UVB sourced by quasars with a comoving $n_q \approx 10^{-5}$ $h^3$ Mpc$^{-3}$ placed either in massive halos or randomly. On scales of $10-30$ $h^{-1} \mbox{ Mpc}$, the flux 3PCF displays hierarchical scaling with the square of the 2PCF, but with an unusual value of $Q \equiv \zeta_{123}/(\xi_{12} \xi_{13} + \xi_{12} \xi_{23} + \xi_{13} \xi_{23}) \approx -4.5$ that reflects the low bias of the Ly$\alpha$ forest and the anti-correlation between mass density and transmitted flux. For halo-based quasars and an ionizing photon mean free path of $\lambda = 300$ $h^{-1} \mbox{ Mpc}$ comoving, UVB fluctuations moderately depress the 2PCF and 3PCF, with cancelling effects on $Q$. For $\lambda = 100$ $h^{-1} \mbox{ Mpc}$ or 50 $h^{-1} \mbox{ Mpc}$, UVB fluctuations substantially boost the 2PCF and 3PCF on large scales, shifting the hierarchical ratio to $Q \approx -3$. We scale our simulation results to derive rough estimate of the 3PCF detectability in observational data sets for the redshift range $z=2.1 - 2.6$. At $r = 10$ $h^{-1} \mbox{ Mpc}$ and 20 $h^{-1} \mbox{ Mpc}$, we predict a signal-to-noise (SNR) of $\sim$ 9 and $\sim$ 7, respectively, for both BOSS and eBOSS, and $\sim$ 37 and $\sim$ 25 for DESI. At $r = 40$ $h^{-1} \mbox{ Mpc}$ the predicted SNR is lower by $\sim$ 3$-$5 times. Measuring the flux 3PCF would be a novel test of the conventional paradigm of the Ly$\alpha$ forest and help separate the contributions of UVB fluctuations and density fluctuations to Ly$\alpha$ forest clustering.

Equations

$$\zeta (r_{12},r_{13},r_{23})=Q\left[\xi (r_{12})\xi (r_{23})+\xi (r_{23})\xi (r_{31})+\xi (r_{31})\xi (r_{12})\right],$$

$${\delta }^F(\mathbf{x})\equiv \frac{F(\mathbf{x})}{\bar{F}(z)}-1,$$

$$\xi (\mathbf{r})=⟨{\delta }^F(\mathbf{x}){\delta }^F(\mathbf{x}+\mathbf{r})⟩,$$

$$\xi (r_{12})=⟨{\delta }_1^F{\delta }_2^F⟩,$$

$$\zeta (r_{12},r_{13},\theta )\equiv ⟨{\delta }_1{\delta }_2{\delta }_3⟩,$$

$$Q(r_{12},r_{13},\theta )=\frac{\zeta (r_{12},r_{13},\theta )}{\xi (r_{12})\xi (r_{13})+\xi (r_{12})\xi (r_{23})+\xi (r_{13})\xi (r_{23})}.$$

$$P(F_1,F_2|{\delta }_1,{\delta }_2)=P(F_1|{\delta }_1)P(F_2|{\delta }_2).$$

$$⟨{\delta }_1^F{\delta }_2^F⟩=\frac{⟨F_1{F}_2⟩-⟨F{⟩}^2}{⟨F{⟩}^2},$$

$$⟨F_1{F}_2⟩$$

$$P({\delta }_1,{\delta }_2)d{\delta }_{1}d{\delta }_2$$

$$⟨F_1{F}_2⟩$$

$$P({\delta }_1,{\delta }_2)d{\delta }_{1}d{\delta }_2$$

$$=\iint \overline{F_1}({\delta }_1)\overline{F_2}({\delta }_2)P({\delta }_1,{\delta }_2)d{\delta }_{1}d{\delta }_2,$$

$$\bar{F}(\delta )=\int F\cdot P(F|\delta )dF.$$

$${\Gamma }_i(x,y,z)\propto L\frac{e^{-d_i/\lambda }}{4\pi d_i^2}.$$

$$Q_i(r)$$

$${\zeta }_{123,i}$$

$${\xi }_{12,i}$$

$${\xi }_{13,i}$$

$${\xi }_{23,i}$$

$$Q=\frac{Q_m}{b}+\frac{b_2}{b^2}$$

$${\delta }^F=b\delta +\frac{b_2}{2}{\delta }^2-\frac{1}{2}{b}_2⟨{\delta }^2⟩$$

Full text

1 Introduction

The Ly$\alpha$ forest arises from the low column density ($N_{\mathrm{HI}}\sim 10^{14}\mbox{ cm}^{-2}$) tenuous gas in mildly overdense regions of the intergalactic medium (IGM). Initially thought to stem from discrete gas clouds along the line of sight (Lynds, 1971; Sargent et al., 1980), a combination of cosmological simulations, analytic models, and improved observations in the mid-1990s established the now standard view of the Ly$\alpha$ forest as tracing a smoothly fluctuating and continuous matter distribution (Cen et al., 1994; Zhang et al., 1995; Hernquist et al., 1996; Miralda-Escudé et al., 1996; Bi & Davidsen, 1997; Croft et al., 1997; Rauch et al., 1997), an inhomogeneous version of the classic Gunn-Peterson effect (Gunn & Peterson, 1965). In this standard picture, the absorbing gas is in photoionization equilibrium with the ionizing background radiation, with Ly$\alpha$ optical depth $\tau=-\ln F\propto n_{H}^{2}T^{-0.7}\Gamma^{-1}$, where $F$ is the continuum-normalized transmitted flux, $n_{H}$ is the total hydrogen density, $T$ is the IGM gas temperature, and $\Gamma$ is the hydrogen photoionization rate. The low density gas that fills most of the volume also obeys a power-law temperature-density relation (Katz et al., 1996; Hui & Gnedin, 1997) and approximately traces the underlying dark matter distribution (Croft et al., 1999; Peeples et al., 2010). This allows a quantitative connection between the Ly$\alpha$ forest and the dark matter density field known as the fluctuating Gunn-Peterson approximation (FGPA, Weinberg et al. 1998).

This picture, together with improving cosmological simulations and observational data sets, has turned the Ly$\alpha$ forest into a powerful probe of matter clustering at redshifts $z=2-4$. Early cosmological studies focused on the line-of-sight power spectrum or the one-point probability distribution function (PDF) of the transmitted flux (Croft et al., 1998, 1999; McDonald et al., 2000; Croft et al., 2002), with a large leap in precision enabled by the enormous sample of quasar spectra from the Sloan Digital Sky Survey (SDSS, McDonald et al. 2005, 2006). The Baryon Oscillation Spectroscopic Survey (BOSS, Dawson et al. 2013) of SDSS-III (Eisenstein et al., 2011) transformed Ly$\alpha$ forest cosmology by providing a dense enough grid of sight-lines to enable measurements of 3-d flux auto-correlation functions across sight-lines (Slosar et al., 2011) and precise measurements of cross-correlations between the Ly$\alpha$ forest and other tracers such as damped-Ly$\alpha$ systems and quasars (Font-Ribera et al., 2012, 2013, 2014). These 3-d measurements are especially powerful for cosmology because they enable measurements of the distance-redshift relation and the Hubble expansion via baryon acoustic oscillations (Busca et al., 2013; Slosar et al., 2013; Delubac et al., 2015; Bautista et al., 2017; du Mas des Bourboux et al., 2017). The large and uniform sample of BOSS spectra also enables highly precise measurements of the line-of-sight power spectrum (Palanque-Delabrouille et al., 2013) and flux PDF (Lee et al., 2015). These 1-d statistics from BOSS and from high-resolution spectra probe small scale dark matter physics, neutrino masses, the amplitude of matter correlations, and the thermal state of the IGM (e.g., Bolton et al. 2008; Viel et al. 2013; Bolton et al. 2014; Palanque-Delabrouille et al. 2015; Rossi 2017; Walther et al. 2019; Khaire et al. 2019).

In this paper we present theoretical predictions for the 3-dimensional 3-point correlation function (3PCF), $\zeta(r_{12},r_{13},r_{23})$, of the Ly$\alpha$ forest at $z=2.3$. Here $\zeta\equiv\langle\delta^{F}_{1}\delta^{F}_{2}\delta^{F}_{3}\rangle$ where $\delta^{F}=(F-\bar{F})/\bar{F}$ is the fractional deviation of the transmitted flux at three positions, denoted by the subscripts, that form a triangle with side lengths $r_{ij}$. The 3PCF is the Fourier transform of the bispectrum, just as the 2-point correlation function (2PCF), $\xi(r)$, is the Fourier transform of the power spectrum. A volume average of the 3PCF yields the skewness $\langle\delta_{S}^{3}\rangle$ of the smoothed $\delta$ field just as a volume average of the 2PCF yields the variance $\langle\delta_{S}^{2}\rangle$. Mandelbaum et al. (2003) and Viel et al. (2004) presented numerical and analytic predictions and measurements of the line-of-sight 1-d flux bispectrum, and Zaldarriaga et al. (2001) investigated a correlation between large scale fluctuations and small scale power that is also a form of 1-d bispectrum. To our knowledge, however, ours is the first investigation of the 3-dimensional 3-point flux correlations. We carry this out using a modified form of the Ly$\alpha$ Mass Association Scheme (LyMAS, Peirani et al. 2014; Lochhaas et al. 2016), which bootstraps results from high-resolution hydrodynamic simulations onto large cosmological $N$-body volumes. Our study is motivated by the prospect of measuring 3-point correlations with the large Ly$\alpha$ forest sample expected from the Dark Energy Spectroscopic Instrument (DESI, DESI Collaboration et al. 2016), which will measure $10^{6}-10^{7}$ Ly$\alpha$ forest spectra over 14,000 deg2, as well as the possibility of first detections with existing data from BOSS and its SDSS-IV successor eBOSS (Dawson et al., 2016).

The Gaussian initial conditions predicted by standard inflationary models have a vanishing 3-point function. However, gravitational instability of Gaussian initial conditions generates a non-vanishing 3-point function at second order in perturbation theory, with the scaling

[TABLE]

where $Q$, often referred to as the reduced 3PCF, is a dimensionless quantity of order unity with moderate dependence on the shape of the matter power spectrum and the shape of the triangle (Fry, 1984). The analogous “hierarchical” relation for moments of the smoothed matter density field is $\langle\delta_{m}^{3}\rangle=S_{3}\langle\delta_{m}^{2}\rangle^{2}$ with $S_{3}\approx 3Q$ (Juszkiewicz et al., 1993). A local bias relation $\delta=f(\delta_{m})$ between the matter density contrast and that of a tracer field preserves the hierarchical form of equation (1) at second order but changes the value of $Q$ and its dependence on triangle shape (Fry & Gastanaga, 1993; Fry, 1994; Juszkiewicz et al., 1995).

We will show that the 3PCF of the Ly$\alpha$ forest scales like the square of the 2PCF as in equation (1), but with an unusual value of $Q$ that reflects the low bias factor of the forest flux fluctuations (Slosar et al., 2011). In the nonlinear and strong clustering regime (0.1 $\lesssim$ $r$ $\lesssim$10 $h^{-1}$ Mpc), $Q$ for galaxies has been observed to be constant at $\approx 1.3$ with no clear dependence on triangle shape (Peebles & Groth, 1975; Groth & Peebles, 1977), consistent with $N$-body simulations of the matter distribution (Fry et al., 1993; Matsubara & Suto, 1994; Scoccimarro et al., 1998; Scoccimarro & Frieman, 1999). On larger scales, observations, simulations, and perturbation theory suggests that galaxies do not strictly show a constant $Q$ but exhibit scale and shape dependence (Jing & Börner, 1998; Scoccimarro et al., 1998; Takada & Jain, 2003; McBride et al., 2011; Hoffmann et al., 2018). The BAO feature has been detected in the 3PCF measurements of BOSS galaxies (Slepian et al., 2017), while other studies focus on the galaxy bispectrum (e.g., Tellarini et al., 2016; Desjacques et al., 2018; Gualdi et al., 2019).

Spatial fluctuations of the ionizing ultraviolet background (UVB) and the IGM temperature-density relation can imprint structure on the Ly$\alpha$ forest in addition to the clustering generated by the density and velocity fields. Some level of spatial variation of $\Gamma$ is inevitable because much of the ionizing background at $z=2-4$ comes from relatively rare quasars, and the expected mean free path of ionizing photons is only $\sim 100-600$ comoving $h^{-1}\mbox{ Mpc}$ (Meiksin & White, 2004; Worseck et al., 2014). Fluctuations of the temperature-density relation at this redshift could arise from the residual effects of inhomogeneous He II reionization (Lai et al., 2006; White et al., 2010; McQuinn et al., 2011). These effects complicate the relation between the Ly$\alpha$ forest and the underlying matter density, and they are a source of systematic uncertainty in cosmological interpretation of Ly$\alpha$ forest clustering. Diagnostics of ionizing background fluctuations or temperature fluctuations are valuable both as direct probes of these physical processes and to help control the cosmological systematics.

Early studies of the impact of UVB fluctuations on the forest focused on the column density distribution and correlation function of Ly$\alpha$ absorption lines (Zuo, 1992a, b; Fardall & Shull, 1993). Croft et al. (1999) and Gnedin & Hamilton (2002) studied the effect of UVB spatial variations on the flux power spectrum and the recovered matter power spectrum and found a negligible effect on small scales but a potential effect on large scales. By further including the effect of quasar lifetimes, Croft (2004) found that UVB fluctuations weakly suppress the flux power spectrum at small scales. Meiksin & White (2004) examined similar effects at redshifts $z>5$, where fluctuations are large because of the short photon mean free path. Recent analytical studies by Pontzen (2014) and Gontcho A Gontcho et al. (2014) demonstrated the scale-dependence of a UVB fluctuation imprint on the flux power spectrum and the resultant broadband distortion to the correlation function of the forest. Suarez & Pontzen (2017) extended these studies to include the effect of quasar emission geometry. The impact of temperature fluctuations from inhomogeneous He II reionization is less well explored, but effects are expected to be present (Lai et al., 2006; White et al., 2010; McQuinn et al., 2011).

In this paper we aim to establish basic theoretical expectations for the flux 3PCF at $z\sim 2.3$ and to investigate how UVB fluctuations affect the flux 2PCF and 3PCF on scales of $\sim 5-50$ $h^{-1}$ Mpc. Because ionizing background fluctuations modulate the Ly$\alpha$ flux with a field that is non-Gaussian and has a different power spectrum than the underlying density field, their impact on the 3PCF could be distinctive. We find that a fluctuating UVB changes the 2PCF and 3PCF of the Ly$\alpha$ forest at all scales to give systematically larger values as the UVB becomes more inhomogeneous. A combination of the 2PCF and 3PCF could then allow better separation between UVB fluctuations and other astrophysical and cosmological parameters. We also use our simulations to give an estimate of the achievable signal-to-noise ratio (neglecting observational noise such as photon noise) of a 3PCF measurement for future and current surveys, in which we predicted a 3PCF detection with a SNR of $\sim 7$ for BOSS and eBOSS and $\sim 25$ for DESI within the redshift range $z=2.1-2.6$.

In §2 we define our notation for the Ly$\alpha$ forest 2PCF and 3PCF measurements. In §3 we explain how we use LyMAS to predict these clustering statistics for a uniform ionizing background and for a fluctuating background sourced by quasars in massive halos or placed at random, with different choices of source volume density and photon mean free path. Section 4 presents our clustering results with uniform and fluctuating UVB and a rough estimate of detectability of the 3PCF. We summarize our findings in §5.

2 Correlation functions

For measurements of the Ly$\alpha$ forest at redshift $z$, we define the flux fluctuations for a pixel at redshift-space position $\mathbf{x}$,

[TABLE]

where $F=e^{-\tau_{\mathrm{Ly}\alpha}}$ is the ratio of the transmitted flux to the quasar continuum and $\bar{F}(z)$ is the mean transmitted flux at redshift $z$. We define the flux 2PCF

[TABLE]

where the average is over all available pixel pairs in a sample of sight-lines with redshift-space separation $\mathbf{r}$. In general, the clustering of the Ly$\alpha$ forest is highly anisotropic due to redshift-space distortions (Slosar et al., 2011). For simplicity, in this paper we will restrict our attention to purely transverse or (in §3.3) nearly transverse pixel separations, so that $\mathbf{r}=|\mathbf{r}|$ refers to the transverse separation of sight-lines. For a set of sight-lines through a simulation, we measure the 2PCF by considering all pairs of sight-lines with transverse separations $r_{12}\rightarrow r_{12}+\Delta r_{12}$ and computing

[TABLE]

where the average includes all transverse pixel pairs along all sight-line pairs.

The formalism for the transverse 3PCF follows similarly to that of the 2PCF:

[TABLE]

where $r_{12}$ is the separation between the first and second line of sight, $r_{13}$ is the separation between the first and third line of sight, and $\theta$ is the angle between the vectors $r_{12}$ and $r_{13}$. The reduced 3PCF, $Q$, can be constructed from the ratio of the 2PCF and 3PCF according to equation (1) as

[TABLE]

3 Simulations and method

3.1 Predicting Ly$\alpha$ forest correlations with LyMAS

Accurately modeling the Ly$\alpha$ forest with hydrodynamic simulations requires resolving the pressure-support scale (Jeans scale) of the diffuse IGM, which is of order $\lambda_{J}\sim 0.25$ $h^{-1}\mbox{ Mpc}$ comoving for a matter overdensity $\delta\approx 10$ (Peeples et al., 2010, eq. 2). Predicting the 3PCF on scales accessible to BOSS and DESI requires simulation volumes of $\sim 1$ Gpc3 or more, and this combination of volume and resolution is impractical with current capabilities. We therefore compute our flux predictions with LyMAS (Peirani et al., 2014), which uses a high-resolution hydrodynamic simulation to compute the conditional PDF, $P(F_{s}|\delta_{s})$, and creates artificial spectra from the density field $\delta_{s}$ of a large volume $N$-body simulation by drawing flux values from $P(F_{s}|\delta_{s})$. Here $F_{s}$ represents the transmitted flux field smoothed in 1-d along the line of sight by the spectral resolution of the survey being modeled, and $\delta_{s}$ represents the matter density field smoothed in 3-d over a scale resolved adequately in the large volume simulation. In the remainder of the paper, we drop the $s$ subscripts and use $F$ and $\delta$ to refer to the smoothed fluxes and matter density contrasts, respectively.

In this paper, as in Lochhaas et al. (2016), we calibrate LyMAS using the Horizon simulation of Dubois et al. (2014) with no AGN feedback (Peirani et al., 2017), and we apply it to a 20483 $N$-body simulation of a (1 $h^{-1}$ Gpc)3 comoving volume that is executed with GADGET2 (Springel, 2005). We adopt line-of-sight Gaussian smoothing of dispersion $\sigma=0.696$ $h^{-1}\mbox{ Mpc}$ comoving, appropriate to BOSS spectral resolution at $z\approx 2.5$, and 3-d Gaussian density smoothing with dispersion $\sigma=0.5$ $h^{-1}$ Mpc. Our simulations use WMAP7 cosmological parameters (Komatsu et al., 2011), where $\Omega_{m}=0.272$, $\Omega_{\Lambda}=0.7284$, $\Omega_{b}=0.045$, $h=0.704$, $\sigma_{8}=0.81$, and $n_{s}=0.967$. We expect that changing to Planck cosmological parameters would have a small impact on our predicted 2PCF and 3PCF but would not qualitatively change our conclusions. For further details, see Lochhaas et al. (2016).

The fundamental assumption of LyMAS is that any correlation between the fluxes arises only from the correlation of the underlying matter distribution. In other words, each draw of the flux value from the conditional PDFs $P(F|\delta)$ is independent, implying

[TABLE]

This approximation breaks down on small scales but becomes more accurate at large separations (Peirani et al., 2014).

Peirani et al. (2014) focused on calculating the flux joint conditional PDFs (Miralda-Escudé et al., 1997) as a model statistic. For calculating flux correlation functions, LyMAS can be simplified. The flux 2PCF can be written generally as

[TABLE]

with

[TABLE]

This expression has no approximations – we can compute $\langle F_{1}F_{2}\rangle$ by integrating over the full joint PDF of the matter density contrasts $\delta_{1}$, $\delta_{2}$ and over the full conditional joint PDF of the fluxes given $\delta_{1}$, $\delta_{2}$. We can now apply the LyMAS ansatz of equation (7) to write

[TABLE]

where the conditional mean flux is

[TABLE]

We therefore obtain the same 2PCF if we deterministically assign fluxes to $N$-body pixels using the conditional mean $\bar{F}(\delta)$ and if we draw from the full conditional $P(F|\delta)$. Averaging over pixel pairs from the simulated density field performs the integral over $P(\delta_{1},\delta_{2})d\delta_{1}d\delta_{2}$ by Monte Carlo integration. A similar argument holds for the 3PCF. We have confirmed numerically that conditional mean fluxes yield the same flux correlation functions as draws from the conditional PDFs, except for the impact of random fluctuations on the latter. Using the conditional mean flux rather than draws from $P(F|\delta)$ has the advantage of producing spectra that are coherent along the line of sight, removing the need for the “percentile field” mapping of Peirani et al. (2014) to create smooth mock spectra. The ‘full LyMAS’ prescription of Peirani et al. (2014) also rescales the Fourier components of the flux field to reproduce the 1-d flux power spectrum of the hydrodynamic simulation, but we omit this step here.

We calibrate the conditional mean flux $\bar{F}(\delta)$ using the $z=2.3$ output of the Horizon-noAGN simulation (Dubois et al., 2014; Peirani et al., 2017), a (100 $h^{-1}\mbox{ Mpc}$)3 comoving volume simulated using RAMSES (Teyssier, 2002) in which the initially uniform grid is adaptively refined down to 1 proper kpc at all times, then sampled on a 2562 grid of sight-lines with the box $z$-axis taken as the line of sight. For our uniform UVB simulation, we choose an HI photoionization rate $\Gamma_{0}$ that yields a mean flux $\bar{F}(z)=0.80$ averaged over all sight-lines, in agreement with observational estimates (Faucher-Giguère et al., 2008; Becker et al., 2013). We also calibrate $\bar{F}(\delta)$ for other choices of the ionizing background intensity $\Gamma$, sampling values of ln($\Gamma_{0}/\Gamma$) = $-1.5$ to 1.5 with a separation of 0.1. To do so we rescale the optical depth $\tau=-\mathrm{ln}F$ of the full resolution Horizon-noAGN spectra by $\Gamma_{0}/\Gamma$, then apply the 0.696 $h^{-1}\mbox{ Mpc}$ line-of-sight smoothing to these rescaled spectra. This method assumes that the neutral hydrogen density is inversely proportional to $\Gamma$, which is an accurate approximation for the diffuse, highly photoionized gas that produces the Ly$\alpha$ forest (Rauch et al., 1997; Peeples et al., 2010). We tabulate $\bar{F}(\delta)$ at values of log10($1+\delta)$ from $-1.275$ to 1.695 in steps of 0.01, yielding a 2-d lookup table from which we can interpolate to find $\bar{F}$ at any value of $\delta$ and $\Gamma_{0}/\Gamma$ within the range studied. For both the Horizon-noAGN simulation and the (1 $h^{-1}$ Gpc)3, 20483 $N$-body simulation, the redshift-space dark matter density field is smoothed with a 3-d Gaussian of dispersion 0.5 $h^{-1}\mbox{ Mpc}$ as described by Peirani et al. (2014). Figure 1 shows $\bar{F}(\delta)$ for a subset of our ln($\Gamma_{0}/\Gamma$) values.

When applied to the dark matter density field of the calibrating hydrodynamic simulation, LyMAS reproduces the 2PCF of the full hydro spectra well but not perfectly, with the largest deviations arising for separations that are elongated along the line of sight (Peirani et al. 2014, Figure 20). Perturbation theory treatments of the Ly$\alpha$ forest consider separate bias factors associated with the density contrast and the line of sight velocity gradient $\eta$ (McDonald, 2003; Seljak, 2012), and it may be possible to improve LyMAS by calibrating fluxes conditioned on both $\delta$ and $\eta$. We leave such an investigation to future work and for this paper note that our predicted amplitudes of the 2PCF and 3PCF could be inaccurate at the $20-30\%$ level based on the comparisons in Peirani et al. (2014) and our investigations with the Horizon simulation. Unfortunately, the hydrodynamic simulation volume is itself too small to characterize this inaccuracy with precision.

We expect that our qualitative conclusions about the dependence of $Q$ on scale and triangle shape and the influence of UVB fluctuations on $\zeta,\xi$ and $Q$ to hold despite this quantitative uncertainty. LyMAS should be considerably more accurate than calculations based on applying the fluctuating Gunn-Peterson approximation (FGPA) to a large $N$-body simulation (e.g., Slosar et al. 2009), which would effectively replace the curves in Figure 1 with linear relations between $-\mathrm{ln}F$ and log10($1+\delta)$. The essential problem with the FGPA for large volume simulations is that the tight relation between optical depth and matter density holds at the Jeans scale of the diffuse IGM, but does not hold between the smoothed matter density contrast and the smoothed Ly$\alpha$ forest spectrum (see Peirani et al. 2014, Figure 4). We therefore regard LyMAS as the most promising method to make predictions for non-linear 3-d structure in the Ly$\alpha$ forest at the $20-50$ $h^{-1}\mbox{ Mpc}$ scales probed well by BOSS and DESI, since full hydrodynamic simulations of the requisite resolution and volume remain impractical.

To create simulated Ly$\alpha$ forest spectra with a fluctuating UVB, we first compute the quantity ln($\Gamma_{0}/\Gamma$) on a uniform 3-d grid of 20 $h^{-1}\mbox{ Mpc}$ spacing in the 1 $h^{-1}$ Gpc simulation cube using the method described below in §3.2, where $\Gamma_{0}$ represents the mean photoionization rate averaged over all points in the grid. At each pixel along each spectrum, we compute ln($\Gamma_{0}/\Gamma$) by linear interpolation among the surrounding grid points, then assign the value of $F$ by linear interpolation on our 2-d table of $\bar{F}(\mathrm{log}_{10}(1+\delta),\mathrm{ln}(\Gamma_{0}/\Gamma))$. We apply a final multiplicative scaling of all $\Gamma$ values such that the mean flux along all spectra is again $\bar{F}=0.8$.

3.2 Implementing a fluctuating ionizing background

To obtain a fluctuating radiation field, we assume quasars as our ionizing sources and place them either randomly in the box or in a random subset of massive DM halos. For the clustered quasar population, we use the DM halos identified by Lochhaas et al. (2016) using a friends-of-friends algorithm (Davis et al., 1985). We place quasars in halos with $M_{h}\geq$ $3\times 10^{12}$ $M_{\odot}$, consistent with the host halo mass inferred from the clustering of BOSS quasars (Font-Ribera et al., 2013; Eftekharzadeh et al., 2015). This mass cut selects $\sim$ 97,000 halos in the (1 $h^{-1}$ Gpc)3 simulation volume. For our fiducial fluctuating UVB model, we adopt a quasar duty cycle of 10%, i.e., we randomly select 10% of these halos to represent active quasars at $z=2.3$. This random selection results in 9606 quasars in the box, which is a comoving volume density of $n_{q}\approx 10^{-5}$ $h^{3}$ Mpc*-3*.

Comparing the clustering results for randomly placed quasars and quasars in massive halos allows us to separate the impact of shot noise and quasar clustering (see Gontcho A Gontcho et al. (2014) and Pontzen (2014) for analytic discussion). In both cases we have simplified reality by assigning all quasar sources the same luminosity rather than drawing from a luminosity function. For randomly distributed quasars of constant luminosity $L_{q}$ and mean volume density $n_{q}$, the mean and variance of the total luminosity emitted in a volume $V$ are $Vn_{q}L_{q}$ and $Vn_{q}L_{q}^{2}$, respectively, because the variance in quasar number for a Poisson distribution is equal to the mean. For randomly distributed quasars drawn from a luminosity function $\phi(L)$, the mean and variance are $V\cdot\int_{0}^{\infty}L\phi(L)dL$ and $V\cdot\int_{0}^{\infty}L^{2}\phi(L)dL$. Taking the quasar luminosity function of Kulkarni et al. (2018) at $z=2.5$, a double power-law with $\phi_{*}\approx 10^{-6}$ ($h^{-1}$ Mpc)3, $\alpha\approx-4$, $\beta\approx-1.75$ (see their Figure 4), we find an rms fractional fluctuation of 0.292 $(V\phi_{*})^{-1/2}$, for $V$ in comoving ($h^{-1}$ Mpc)3, which is equal to that of a constant $L_{q}$ population of volume density $n_{q}=1.17\times 10^{-5}$ ($h^{-1}$ Mpc)3. Our fiducial case of $n_{q}\approx 10^{-5}$ ($h^{-1}$ Mpc)3 should therefore be representative of the UVB fluctuations expected from the observed quasar population at this redshift.

We also vary the space densities for random and clustered quasar populations by a factor of eight higher and lower to map out the dependence of the 2PCF and 3PCF on the UVB emissivity fluctuations. The contribution of galaxies to the UVB at this redshift ($z=2.3$) is uncertain, but it could potentially be non-negligible (Haardt & Madau, 2012; Khaire & Srianand, 2019). If galaxies make a large contribution to the UVB, then the UVB would be smoother than our $n_{q}=8\times 10^{-5}$ ($h^{-1}$ Mpc)3 case, since the shot noise would be lower and the clustering bias of galaxies is weaker than that of quasars.

The other critical parameter controlling UVB fluctuations is the mean free path $\lambda$ of ionizing photons. A smaller $\lambda$ implies that the ionizing flux at a given location comes from a smaller number of sources and is therefore subject to larger fluctuations. The mean free path is challenging to estimate observationally because absorption is dominated by systems with $\tau\sim 1$ at the Lyman limit, and these systems are relatively rare ($\sim 1$ per quasar sight-line) and their column densities are difficult to measure because their Ly$\alpha$ absorption is saturated. O’Meara et al. (2013) find $\lambda\sim 570$ $h^{-1}\mbox{ Mpc}$ at $z=2.44$, and Fumagalli et al. (2013) find $\lambda\sim 300$ $h^{-1}\mbox{ Mpc}$ at $z=3.0$ (see Worseck et al. (2014) for a broader compilation). For our calculations, we consider $\lambda$ = 300, 100, and 50 $h^{-1}\mbox{ Mpc}$. Our $\lambda$ = 300 $h^{-1}\mbox{ Mpc}$ is closest to (but larger than) than observational estimates near $z=2.3$, while the smaller values help illustrate behavior with stronger UVB fluctuations, which is useful for intuitive understanding and may be relevant at higher redshifts. It would be useful to have results for a still larger value of $\lambda$, but even our 1 $h^{-1}$ Gpc box is not large enough to do this.

We assume that quasars are radiating isotropically at a constant luminosity $L$, so that the photoionization rate from quasar $i$ located a distance $d_{i}$ away from a point $(x,y,z)$, including periodic boundary conditions, is given by

[TABLE]

The value of $L$ is fixed implicitly by choosing the mean ionization rate to yield $\bar{F}=0.8$ averaged over all sight-lines. We do not account for clustering of absorbers in the same large scale structure that hosts the quasars and the Ly$\alpha$ forest, as this would require a much more complex radiative transfer calculation. In the analytic treatment of Gontcho A Gontcho et al. (2014) and Pontzen (2014), the impact of absorbers is roughly equivalent to modifying the quasar bias factor, so results with clustered Lyman limit absorption might be intermediate between our clustered and random quasar cases. However, a fully non-linear calculation with clustered absorption remains a goal for future work.

Figure 2 shows the distribution of quasars in massive halos for a slice in the 1 $h^{-1}$ Gpc box for our fiducial 10% duty cycle. The UVB flux at a given location is dominated by the nearest number of quasars, $N_{q}\sim(4\pi/3)\lambda^{3}n_{q}\sim 1131(n_{q}/10^{-5})(\lambda/300)^{3}$. Figure 3 shows the combined effect of density and UVB fluctuations on the transmitted flux for a random selected sight-line through the box. As expected, the fractional flux variations ln($\Gamma/\Gamma_{0}$) become much larger for the shorter $\lambda$ values. Although the structure of the Ly$\alpha$ forest spectrum is imprinted principally by the density fluctuations, it is modulated by the UVB fluctuations. Near $z$=800 $h^{-1}\mbox{ Mpc}$, a large scale overdensity is also a location of a concentration of quasars and thus a peak in the UVB intensity. The Ly$\alpha$ forest absorption is therefore reduced relative to the uniform background case (see zoom panel), more so for the shortest $\lambda$. With clustered ionizing background sources, density and UVB fluctuations tend to have opposite impact on the Ly$\alpha$ forest absorption. However, even for $\lambda=50$ $h^{-1}\mbox{ Mpc}$, the scale of UVB fluctuations is much larger than that of the density fluctuations that produce order unity Ly$\alpha$ flux variations.

Figure 4 shows the PDF of transmitted flux from all sight-lines through the box. This PDF is remarkably insensitive to the presence of UVB fluctuations. However, we will show that these fluctuations have a significant impact on the flux 2PCF and 3PCF.

3.3 Calculating the Ly$\alpha$ forest clustering

There are 65,536 sight-lines in our 1 $h^{-3}$ Gpc3 $N$-body box, with a minimum sight-line separation of $ds$ = 3.91 $h^{-1}$ Mpc and each spectrum consisting of 4096 pixels. We currently only correlate sight-lines and pixels at the same redshifts (or $z$ positions, i.e. the planes of the triplets are perpendicular to the line of sight), with sight-line separations up to a maximum of 60 $h^{-1}$ Mpc.

We select triplets with roughly equal side lengths, $r\sim r_{12}\sim r_{13}$ and for three different triangle opening angles $\theta$ = 90*∘, 60∘, and 20∘, each with an angle margin of $\pm$ 5∘*. Recall that $\theta$ is defined as the angle between the vectors $r_{12}$ and $r_{13}$. We choose $r_{13}$ separations spanning from 0.8$r_{12}$ to 1.2$r_{12}$, where $r_{12}=ds,2ds,3ds...15ds=3.91,7.82,11.73...58.65$ $h^{-1}\mbox{ Mpc}$ . For each $r_{12}$, we iteratively use every sight-line in the box as a primary sight-line, then randomly select one of the four possible second sight-line located $r_{12}$ away, and finally select all possible third sight-lines to complete the triplet within the $r_{13}$ and angle ranges. Since sight-lines are repeated, the error bars of the 2- and 3-point correlation functions at various scales are correlated. Figure 5 shows an example of a triplet configuration for each $\theta$, and Figure 6 shows the total number of triplets in our box as a function of separation.

We calculate the 2PCF, 3PCF, and $Q$ of the triplets according to Equations (1), (3) and (4), in bins of 4 $h^{-1}$ Mpc. The final $Q$ values are obtained using all sight-lines in the entire box. To estimate the error bars, we divide the triplets into nine subvolumes and calculate the 2PCF and 3PCF using all sight-lines in each subvolume. The subvolumes are divided in $X$ and $Y$, but not in $Z$, so each subvolume is essentially a long narrow rectangular prism. The reduced 3PCF for a subvolume $i$, $Q_{i}$, is obtained accordingly using the respective correlation functions.

[TABLE]

where $N_{i}(r_{b})$ is the number of triplets in each subvolume $i$ in bin $r_{b}$. We estimate our error bars as the standard deviation in the correlation functions among the subvolumes divided by the square root of the number of subvolumes. They therefore represent our estimate of the uncertainty in the theoretical prediction from the full (1 $h^{-1}$ Gpc)3 simulation volume. We discuss the source of statistical uncertainty in our predictions below, especially in Appendix A, concluding that it is dominated by cosmic variance of large scale structure within our survey volume.

4 Results

4.1 Reduced 3PCF for a uniform ionizing background

We show the reduced 3PCF in a uniform ionizing background as a function of $r$ for all three triangle shapes we investigated in Figure 7. The $Q$ value has little dependence on the triangle shape and remains approximately constant at $\sim-4$ despite the 3PCF changing by more than two orders of magnitude (see Figure 8). Compared to galaxies, which have a positive and small $Q_{m}$ ($\approx 1.3$) (Peebles & Groth, 1975; Groth & Peebles, 1977), the value of $Q$ for the Ly$\alpha$ forest is negative and large. The negative value of $Q$ (and the 3PCF) arises because the forest is in absorption, so that high density produces low flux. The large amplitude of $Q$ reflects the low bias factor of the forest. With a local bias model of the forest at second order,

[TABLE]

we assume the forest flux fluctuation is related to the DM overdensity by

[TABLE]

(Fry & Gastanaga, 1993; Fry, 1994; Juszkiewicz et al., 1995). We get $Q=-7.6$ for $b_{2}=0$ when adopting $Q_{m}=1.3$ and $b=-0.17$ (Slosar et al., 2011); reproducing our simulation value of $Q\sim 4$ requires $b_{2}\approx 0.1$. Thus small values of $b$ and $b_{2}$ for the forest naturally give rise to a large value of $Q$.

4.2 Impact of ionizing radiation fluctuations and source clustering

We compare the 2PCF and 3PCF of the Ly$\alpha$ forest in different fluctuating UV backgrounds and with different triplet configurations of sight-lines. Figure 8 and 9 show the results for quasars found in massive halos and for randomly-distributed quasars, respectively. UVB fluctuation changes the clustering at all scales and produces increased signal as $\lambda$ gets smaller. Although the flux PDF remains unchanged with UVB fluctuations (Figure 4), the 2PCF and 3PCF are clearly changed.

For a clustered quasar source population with $\lambda=300$, the 2PCF of the forest is moderately suppressed by a factor of $\sim 1.2-2$, whereas the 2PCF from unclustered quasar sources is similar to that for a uniform background at small scales but slightly enhanced at large scales. In the unclustered source cases, flux variations are a source of additional large scale structure in the forest. For halo-based quasars on the other hand, high density regions also have higher UVB, and the cancellation suppresses clustering overall.

For $\lambda=100$, the UVB fluctuations are larger in amplitude, and for unclustered quasar sources the flux 2PCF and 3PCF are enhanced significantly at all scales. With clustered sources there is again partial cancellation, but at large scales the 2PCF is now enhanced relative to the uniform UVB instead of suppressed. For $\lambda=50$, the flux 2PCF and 3PCF are dramatically enhanced at all scales; in this case the large scale clustering of the forest is dominated by UVB variations rather than gas density fluctuations. One might hope that the transition to a different clustering origin would lead to a sharp departure from the hierarchical relation of the 3PCF and 2PCF, but the $\lambda=50$ and $\lambda=100$ cases still have approximately constant $Q$ out to $r=30$ $h^{-1}\mbox{ Mpc}$, with a moderate reduction from $|Q|\approx-4.5$ from the smooth UVB case to $|Q|\approx-3$.

There could be significant changes in $Q$ at larger scales, but the error bars from our finite simulation volume are too large to tell. In all cases the triangle shape has only moderate impact on $Q$, though for $\theta=20^{\circ}$ the error bars at large $r$ are reduced because $r_{23}$ remains relatively small, so $\zeta_{123}$ and $\xi_{23}$ are larger and better measured. For this triangle shape we see only moderate scale dependence of $Q$ out to $r=60$ $h^{-1}\mbox{ Mpc}$, and the value of $|Q|$ is lower for the $\lambda=50$ and $\lambda=100$ cases relative to the $\lambda=300$ and uniform cases.

Figure 10 compares the impact of clustered vs. unclustered radiation sources more directly, for triplets with $\theta=60^{\circ}$. Quasars in massive halos tend to produce weaker forest clustering than randomly-placed quasars for all $\lambda$ values. Because hierarchical behavior is a “special” consequence of gravitational instability of Gaussian initial conditions, we anticipated that we might see sharp scale-dependence of $Q$ associated with the scale of $\lambda$. However, this is not evident within our errors. Larger simulation boxes are needed to further test this conjecture. As discussed further in Appendix A, our results are stable against different random realizations of quasar distribution. Rather than being limited by the number of sight-lines or triplets in our box, our measurements are limited by the variance due to large scale structures. We therefore need larger simulation volumes rather than more complete sampling of this simulation.

4.3 Impact of shot noise

Shot noise from the rarity of the ionizing sources can affect and complicate interpretation of Ly$\alpha$ forest clustering measurements. We investigate the impact of shot noise by changing the number of quasars in the different fluctuating backgrounds by a factor of eight from the fiducial $N_{q}=10,000$. We again assume either a halo-based or random quasar distribution. For halo-based quasars, we use the same set of DM halos and the same mass cut of $M_{h}\geq$ $3\times 10^{12}$ for the quasar hosts as before, but vary the quasar duty cycle to 80% and 1.25%. For random quasar distributions, we specify the desired numbers of quasars exactly (either higher or lower by a factor of eight from the fiducial) and randomly assign their $(X,Y,Z)$ positions in the box. We follow the same steps listed in §3.1 to generate the resultant Ly$\alpha$ forest fluxes using the new UVB grids, making sure to renormalize the new fluxes to the observed mean flux of 0.8. The resultant flux histograms (1-point PDFs) are nearly unchanged for different quasar densities, similar to Figure 4.

Figure 11 shows the 2PCF, 3PCF, and $Q$ measurements for clustered quasar populations at the three volume densities for an equilateral triplet configuration in a $\lambda=300$ $h^{-1}\mbox{ Mpc}$ ionizing background. Shot noise adds broadband power on all scales, giving rise to the largest clustering signal for the lowest $n_{q}$ (green points) while being suppressed in the largest $n_{q}$ (blue points). Values of $Q$ are not sensitive to $n_{q}$, at least relative to our error bars.

Figure 12 shows the $\lambda=50$ and $\lambda=100$ cases with different $n_{q}$. The clustering signals increase for shorter $\lambda$ or lower $n_{q}$ as expected. Lines with the same color show $(n_{q},\lambda)$ combinations chosen to have the same $N_{q}=4/3\pi\lambda^{3}n_{q}$; blue lines have $N_{q}\sim 5$ and green lines have $N_{q}\sim 42$. The value of $N_{q}$ clearly separates these combinations, though it does not fully determine the forest clustering. The case of low quasar density $n_{q}=1.25\times 10^{-8}$ and $\lambda=100$ is the one combination that shows a very different value of $Q$ at $\sim-1.5$. One can see a suggestion of reduced $|Q|$ for low $n_{q}$ and $\lambda=300$ in Figure 12. Randomly placed quasars follow the same behavior and trend, except with larger correlation function values and smaller $|Q|$ amplitudes than clustered quasars, reflecting what we see in Figure 10.

4.4 Detectability of the 3PCF

Slosar et al. (2011) made the first detection of the 3-d 2PCF redshift-space distortion in the Ly$\alpha$ forest. With subsequent data from BOSS, the 2PCF has been measured with increasing precision to constrain the BAO peak (Busca et al., 2013; Slosar et al., 2013; Delubac et al., 2015; Bautista et al., 2017; du Mas des Bourboux et al., 2017; Blomqvist et al., 2019; de Sainte Agathe et al., 2019). To our knowledge, no 3PCF measurement of the forest has been made. We attempted to measure the 3PCF using recent data from the CLAMATO 3D Ly$\alpha$ forest tomography survey (Lee et al., 2018). Although we obtained a clear 2PCF signal, the 3PCF measurement is consistent with noise. We can therefore ask whether the 3PCF should be detectable in current and future surveys, e.g., in BOSS, eBOSS, and DESI.

To answer this question, we estimate the expected signal-to-noise ratio (SNR) of a 3PCF detection using specifications that approximate these current spectroscopic surveys. The volume probed by these surveys is much larger than our 1 ($h^{-1}$ Gpc)3 box. For example, the BOSS survey volume for the redshift range $z=2.1-2.6$ over 104 deg2 corresponds to $\sim 21$ $h^{-3}$ Gpc3. However, the BOSS sampling density is far lower; given 114,600 quasars between $2.1<z<2.6$ distributed over an area of 9376 deg*-2* (BOSS DR12, Alam et al., 2015), its sampling density is $\sim$ 12 quasars deg*-2*, corresponding to a comoving surface density of $\Sigma=2.6\times 10^{-3}$ $h^{2}$ Mpc*-2* at $z=2.3$. Our previous analyses use 2562 sight-lines through our box, which gives a comoving surface density of $\Sigma=6.55\times 10^{-2}$ $h^{2}$ Mpc*-2*.

For our SNR estimate, we assume an observed surface density of 10 quasars deg*-2* at $z=2.3$, or $\Sigma=2\times 10^{-3}$ $h^{2}$ Mpc*-2*, roughly comparable to that of BOSS. This translates to 2000 sight-lines for our 1 $h^{-1}$ Gpc box. A typical Ly$\alpha$ forest region spanning the range from the quasar’s Ly$\alpha$ to Ly$\beta$ emission lines is $\sim 300$ $h^{-1}\mbox{ Mpc}$ long, three times shorter than our simulation sight-lines. As we analyzed the entire sight-lines, this means we effectively have 6000 sight-lines or a 3$\times$ larger effective volume. We estimate the noise for BOSS, eBOSS, and DESI by scaling the noise from our box by $1/N^{1/2}_{\mathrm{sightline}}$. This scaling is appropriate because errors in widely separated regions (i.e., larger than our simulation box) should be independent.

After selecting 2000 random sight-lines through our box, we consider loosely equilaterial triangle configurations with a fractional width of $w=\Delta r/r$ = 0.2 at three separations of $r=10,20$, and 40 $h^{-1}$ Mpc. For instance, for $r$ = 10 $h^{-1}$ Mpc, we consider sight-lines that are between 8 $h^{-1}$ Mpc and 12 $h^{-1}$ Mpc away from the primary sight-line. The same fractional width is also applied when we correlate pixels from the sight-line triplets, such that we are not limited to strictly face-on pixel triplets. In other words, we correlate pixels in sight-line $i$ with pixels that are located between $z_{i}-wr$ and $z_{i}+wr$ in the other two sight-lines.

We ran 50 realizations in which we chose 2000 random sight-lines from our uniform UVB box and measure the 3PCF using all sight-line and pixel triples that satisfy the shape criterion mentioned above. These 50 realizations have a mean of 15, 271, and 333 sight-line triplets at $r=10,20$, and 40 $h^{-1}$ Mpc, respectively. We take the dispersion among these 50 realizations to represent the rms error of $\zeta_{123}(r)$ expected for 6000 forest spectra with a surface density of 10 deg*-2* at $z=2.3$. We compute the SNR as the ratio of the mean $\zeta_{123}$ to this dispersion. This calculation implicitly assumes that the statistical errors at this quasar surface density are dominated by sparse sampling of the available structure and not by cosmic variance of the structure itself. We believe this assumption is justified, but we have not demonstrated it.

Another source of noise in observational data is photon noise in the spectra. The photon noise per pixel can be reduced by smoothing the spectra, though this also reduces the number of independent pixel triplets. To help assess this issue, we compute the SNR for our full resolution spectra, for spectra that are boxcar-smoothed over 4 or 16 pixels (each roughly comparable to a BOSS pixel), and for spectra binned over 4 or 16 pixels, which therefore have a factor of 4 or 16 smaller pixel count. To isolate the effect of smoothing from that of number of pixel triplets, we also calculate the SNR after simply choosing every 4th or 16th pixel from that full resolution spectra.

Figure 13 shows our results. For 6000 sight-lines we expect a SNR of $\sim 2$ at $r=10$ $h^{-1}\mbox{ Mpc}$, $\sim 1.5$ at $r=20$ $h^{-1}\mbox{ Mpc}$, and $\sim 0.4$ at $r=40$ $h^{-1}\mbox{ Mpc}$. For the 114,600 sight-lines within the redshift range $z=2.1-2.6$ found in BOSS DR 12 (Alam et al., 2015), we predict a SNR higher by $(114,600/6000)^{1/2}\approx 4.4$, so roughly 8.7, 6.6, and 1.75 at these separations. The eBOSS survey has a comparable surface density of 13.8 quasars deg*-2* within the redshift range $z=2.1-2.6$. For the 129,975 sight-lines within the above redshift range found in eBOSS DR 14 quasar catalog (Pâris et al., 2018), we predict a very similar SNR of roughly 9.3, 7, and 1.9 for these separations. Remarkably, smoothing, binning, or sampling with a scale of 4 or 16 pixels has essentially no impact on the SNR. This is encouraging as it implies that one could reduce photon noise by smoothing up to 16 pixels affecting the SNR. One can therefore either bin or average a larger number of (individually noisier) pixel triplets to reduce photon noise in the observed spectra. Binning has the additional attraction of reducing the computational demands of the 3-point measurement without loss of sensitivity.

We next consider the effect of changing the quasar surface density on the SNR by repeating our simulations for 4000 and 8000 sight-lines, corresponding to increasing the quasar surface density by two (giving 20 quasars deg*-2*) and four (giving 40 quasars deg*-2*). As a comparison, the DESI survey will have a surface density of 50 quasars deg*-2* at $z>2.1$. For the redshift range $z=2.1-2.6$, its surface density is estimated to be 30 quasars deg*-2* (based on Figure 3.17 of DESI Collaboration et al., 2016), resulting in 420,000 total sight-lines.

Since we again use the entire 1 $h^{-1}$ Gpc path length, our 4000 and 8000 unique sight-lines effectively result in 12,000 and 24,000 usable Ly$\alpha$ forest spectra. We show the results in Figure 14. The trend with smoothing the spectra is similar as before, so we only show the comparison at full resolution for brevity. For 8000 sight-lines in our box we expect a SNR of $\sim$ 9, 6, and 2 at $r=10,20$, and 40 $h^{-1}\mbox{ Mpc}$, respectively. To extrapolate our results to DESI, we scale our SNR by $(420,000/24,000)^{1/2}\approx 4$ to get an expected SNR of 37, 25, and 8 at these three separations.

5 Conclusion

The standard picture of the Ly$\alpha$ forest is one where the low-density gas in the IGM remains in photoionization equilibrium with the ionizing background and obeys a tight temperature-density relation. Density fluctuations are thought to dominate the structure of the forest, making the Ly$\alpha$ forest a tracer of large-scale structure and a powerful cosmological tool. In principle however, flux fluctuations can arise from other sources such as spatial variations in the IGM mean temperature and the ionizing background.

We have used LyMAS (Peirani et al., 2014; Lochhaas et al., 2016) to make the first predictions of the 3-d 3PCF fo the Ly$\alpha$ forest, $\delta_{123}=\langle\delta^{F}_{1}\delta^{F}_{2}\delta^{F}_{3}\rangle$, bootstrapping results from the (100 $h^{-1}\mbox{ Mpc}$)3 Horizon hydrodynamic simulation (Dubois et al., 2014) into a (1 $h^{-1}$ Gpc)3 DM-only simulation. We introduce a simplified “conditional mean” formulation of LyMAS, which yields the same results for flux correlation functions as the original “conditional PDF” formulation but makes it much easier to implement effects of a fluctuating UVB. To derive a fluctuating radiation field, we assume quasars as our ionizing sources with various radiation mean free paths, and we either randomly distribute them in space or place them in massive DM halos. For our three-point clustering measurements, we focus on triangle configurations with $r_{12}\sim r_{13}\leq 60$ $h^{-1}$ Mpc, and with opening angles $\theta=90^{\circ},60^{\circ},$ and 20∘.

The predicted 3PCF of the Ly$\alpha$ forest approximately follows the hierarchical behavior expected for matter clustering from Gaussian initial conditions: $\zeta_{123}=Q(\xi_{12}\xi_{13}+\xi_{12}\xi_{23}+\xi_{13}\xi_{23})$, with $Q$ approximately independent of scale. For a uniform UVB, we find $Q\approx-4$ to $-5$ on scales of $10-30$ $h^{-1}\mbox{ Mpc}$, with a weak dependence on triangle size and shape. The large value of $|Q|$ (compared to $|Q|\approx 1$ for matter) likely arises from the low bias factor of the forest, while the negative sign arises because higher densities produce lower fluxes. Even with a (1 $h^{-1}$ Gpc)3 simulation volume, our predictions become noisy on scales larger than $r\sim 30$ $h^{-1}\mbox{ Mpc}$.

For a fluctuating UVB, we consider three values of quasar volume density $n_{q}$ and three values of the ionizing photon mean free path $\lambda=300$ $h^{-1}\mbox{ Mpc}$, 100 $h^{-1}\mbox{ Mpc}$, and 50 $h^{-1}\mbox{ Mpc}$ (comoving). Even the longest of these $\lambda$ is shorter than the observationally inferred value at $z=2.3$, but we would need larger simulation volumes to model larger $\lambda$. Our shorter $\lambda$ values amplify UVB fluctuations to a level expected at higher redshifts; the Worseck et al. (2014) estimate corresponds to a comoving $\lambda\approx 350$ $h^{-1}\mbox{ Mpc}$ at $z=3$ and 130 $h^{-1}\mbox{ Mpc}$ at $z=4$.

For randomly-placed quasars, UVB fluctuations boost the 2PCF and 3PCF on all scales, with larger enhancements for smaller $\lambda$ or smaller $n_{q}$ as expected. With $\lambda=300$ $h^{-1}\mbox{ Mpc}$ the enhancements are small. With $\lambda\leq 100$ $h^{-1}\mbox{ Mpc}$ the large scale enhancements are a factor of two or more, making UVB fluctuations the dominant source of large scale flux correlations. The value of $Q$ remains approximately constant on the scales we can reliably measure, with a somewhat smaller $|Q|\sim 3$ for the smaller $\lambda$ values. For halo-based quasars, a fluctuating UVB with $\lambda=300$ $h^{-1}\mbox{ Mpc}$ depresses the 2PCF and 3PCF on all scales relative to a uniform UVB because overdense regions have a higher average UVB that counteracts the higher average IGM density. For $\lambda=100$ $h^{-1}\mbox{ Mpc}$ or $\lambda=50$ $h^{-1}\mbox{ Mpc}$ , the 2PCF and 3PCF are higher than those for a uniform UVB but lower than those for randomly placed quasars. The value of $|Q|$ is nearly unchanged for $\lambda=300$ $h^{-1}\mbox{ Mpc}$ , and it is again moderately reduced (to $|Q|\sim 3$) for $\lambda=50$ or 100 $h^{-1}\mbox{ Mpc}$ . For $\lambda=300$ $h^{-1}\mbox{ Mpc}$ and halo-based quasars, raising $n_{q}$ from 10*-5* ($h^{-1}\mbox{ Mpc}$)3 to 8 $\times$ 10*-5* ($h^{-1}\mbox{ Mpc}$)3 further depresses the 2PCF and 3PCF, while lowering $n_{q}$ to 1.25 $\times$ 10*-6* ($h^{-1}\mbox{ Mpc}$)3 strongly boosts the predicted correlation functions and reduces $|Q|$ to $\sim 2$.

Because hierarchical behavior of the 2PCF and 3PCF is a distinctive prediction of gravitational instability and Gaussian initial conditions, we hoped that we might see a marked transition to a scale-dependent $|Q|$ on scales where UVB fluctuations become a significant driver of flux correlations. However, we do not see a clear sign of such a transition in our results. Unfortunately our simulation volume is too small to yield precise 3PCF measurements on scales $\geq\lambda$. Confirming or refuting the conjecture of a scale-dependent $|Q|$ from UVB fluctuations will require further studies with larger simulation volumes.

Finally, we derive a rough estimate of the detectability of the 3PCF in data sets such as BOSS, eBOSS, and DESI. We reduce the sight-line density to values comparable to these surveys, consider loosely equilateral triangle configurations that are approximately transverse to the line of sight, and assume that the SNR will scale as $N_{q}^{-1/2}$, where $N_{q}$ is the number of quasar sight-lines. In the absence of observational noise, we estimate SNR $\sim$ 7 and $\sim$ 9 for a BOSS- and eBOSS-like data set at $r=10$ $h^{-1}\mbox{ Mpc}$ and 20 $h^{-1}\mbox{ Mpc}$ , increasing to $\sim$ 37 and $\sim$ 25 for DESI. At $r=40$ $h^{-1}\mbox{ Mpc}$ the predicted SNR is lower by a factor of $\sim$ 3$-$5. Smoothing or binning the spectra over a scale of 16 BOSS-like pixels barely alters the SNR of the 3PCF measurement, which should simplify observational analyses.

Higher-order moments of large-scale structure contain richer and more complex information than two-point statistics alone. Bispectrum-like measurements along individual lines of sight already show promise as tests of the gravitational instability paradigm for the Ly$\alpha$ forest and constraints on other sources of structure (Zaldarriaga et al., 2001; Fang & White, 2004). Dense, wide-area spectroscopic surveys such as BOSS, eBOSS, and DESI offer the prospect of measuring 3-point correlations of the Ly$\alpha$ forest in 3-dimensional redshift space. These measurements can provide new diagnostics of non-gravitational physics affecting the Ly$\alpha$ forest, and reproducing higher-order statistics will allow more confident use of Ly$\alpha$  forest BAO as a probe of dark energy.

Acknowledgment

We thank J. Blaizot, J. Devriendt, Y. Dubois, and C. Pichon for their collaborative contributions to the Horizon hydrodynamic simulation used in this work. We also thank Cassandra Lochhaas for providing her list of dark matter halos and Debopam Som for providing the number of eBOSS Ly$\alpha$ quasars between $z=2.1$ and $z=2.6$.

This work was supported by the U.S. Department of Energy, Office of Science, Office of High Energy Physics under Award Number DE-SC-0011726 and NSF grant AST-1516997. It was also supported by collaborative visits funded by the Cosmology and Astroparticle Student and Postdoc Exchange Network (CASPEN). This work has been done in part within the Labex ILP (reference ANR-10-LABX-63) part of the Idex SUPER, and received financial state aid managed by the Agence Nationale de la Recherche, as part of the programme Investissements d’avenir under the reference ANR-11-IDEX-0004-02. TS was funded by Conacyt (Consejo Nacional de Ciencia y Technologia) and UCL (University College London). WZ was supported by the Beatrice and Vincent Tremaine Fellowship.

This research made use of Astropy, a community-developed core Python package for Astronomy (Astropy Collaboration, 2013).

Appendix A

We check to make sure that our measurements of the correlation functions are stable against different realizations of the quasar distributions. Figure 15 shows a comparison of the $Q$ amplitude for three extra realizations of halo-based and randomly distributed quasars for the fiducial quasar volume density $n_{q}=10^{-5}$ $h^{3}$ Mpc*-3*. The different realizations display the same overall trend in the reduced 3PCF values in the different UVB, with variations falling within the error bars.

We next investigate if our clustering measurements are limited by statistics or cosmic variance. We divide our 1 $h^{-1}$ Gpc box with a smooth UV background into nine equal subvolumes and assign each sight-line triplet to a random and the correct subvolume. Assigning triplets to random subvolumes in principle lets us average out the variance due to large scale structure. The correct subvolume assignment is done based on the position of the primary sight-line regardless of the positions of the second and third sight-lines. The random subvolume assignment results in a uniform number of triplets in each subvolume.

We then compare the errors bars of the correlation functions measured from triplets using the correct and random subvolume assignment, which is shown in Figure 16. We obtain overall larger fractional errors for the 2PCF, 3PCF, and $Q$ when the sight-line triplets are distributed correctly compared to when they are distributed randomly. This is especially true at increasingly large scales, $r\gtrsim$ 10 $h^{-1}$ Mpc. This suggests that our clustering measurements are limited by variance due to large scale structure, rather than by the sampling of these structures from the available sight-lines in our box.

We also compare the error bars on the correlation functions from varying the number of sight-lines in our box, in which we use all sight-lines, half of all sight-lines, and a quarter of all sight-lines. As shown in Figure 17, the fractional errors of the correlation functions are approximately the same whether we use all or a quarter of the available sight-lines. This further suggests that we are not limited by the number of our triplets.