Large-scale clustering of Lyman-alpha emission intensity from SDSS/BOSS

(Abridged) We detect the large-scale structure of Lya emission in the Universe at redshifts z=2-3.5 by measuring the cross-correlation of Lya surface brightness with quasars in SDSS/BOSS. We use a million spectra targeting Luminous Red Galaxies at z<0.8, after subtracting a best fit model galaxy spectrum from each one, as an estimate of the high-redshift Lya surface brightness. The quasar-Lya emission cross-correlation we detect has a shape consistent with a LambdaCDM model with Omega_M =0.30^+0.10-0.07. The predicted amplitude of this cross-correlation is proportional to the product of the mean Lya surface brightness,, the amplitude of mass fluctuations, and the quasar and Lya emission bias factors. Using known values, we infer(b_alpha/3) = (3.9 +/- 0.9) x 10^-21 erg/s cm^-2 A^-1 arcsec^-2, where b_alpha is the Lya emission bias factor. If the dominant sources of Lya emission are star forming galaxies, we infer rho_SFR = (0.28 +/- 0.07) (3/b_alpha) /yr/Mpc^3 at z=2-3.5. For b_alpha=3, this value is a factor of 21-35 above previous estimates from individually detected Lya emitters, although consistent with the total rho_SFR derived from dust-corrected, continuum UV surveys. 97% of the Lya emission in the Universe at these redshifts is therefore undetected in previous surveys of Lya emitters. Our measurement is much greater than seen from stacking analyses of faint halos surrounding previously detected Lya emitters, but we speculate that it arises from similar Lya halos surrounding all luminous star-forming galaxies. We also detect redshift space anisotropy of the quasar-Lya emission cross-correlation, finding evidence at the 3.0 sigma level that it is radially elongated, consistent with distortions caused by radiative-transfer effects (Zheng et al. (2011)). Our measurements represent the first application of the intensity mapping technique to optical observations.

of redshifts (e.g., Hu & McMahon 1996, Keel et al. 1999, Fujita et al. 2003, Cowie et al. 2010. Another potentially useful technique is that of intensity mapping (e.g., Carilli 2011, Peterson & Suarez 2012, Pullen et al. 2014, which seeks to map the large-scale structure using one emission line or more (see e.g., Wyithe & Morales 2007, Visbal & Loeb 2010, without resolving individual sources (such as galaxies or gas clouds). By measuring this structure, one is sensitive to all clustered emission, without the observational biases which arise from source detection and luminosity measurement (such as detection limits, determination of backgrounds and finite aperture size). In this paper, we seek to perform the first cosmological measurement of intensity mapping in the Lyα line, using a large dataset of spectra from the Sloan Digital Sky Survey III (SDSS-III, Eisenstein et al. 2011) Baryon Oscillation Spectroscopic Survey (BOSS, Dawson et al. 2013). We use spectra that were targeted at massive galaxies at z < 0.8. After subtracting best fit model galaxy spectra, we expect that any high redshift Lyα emitters that are within the fiber aperture result in a residual flux, present also in sky fibers. Even if not detectable as individual sources, we can search for large-scale structure in this emission by determining its spatial cross-correlation function with the positions of BOSS quasars, which are tracers of structure with a known bias factor (e.g., White et al. 2012) at redshifts z > 2, where the Lyα emission line is in the optical part of the spectrum.
Following the early prediction by Partridge & Peebles (1967) that galaxies should be detectable at high redshift from their Lyα emission line, many surveys have been designed to detect individual galaxies as sources of Lyα emission. These include narrow band imaging (e.g., Ouchi et al. 2003, Gronwall et al. 2007, serendipitous slit spectroscopy (e.g., Cassata et al. 2011) and integral field spectroscopy (e.g., van Bruekelen et al. 2005, Blanc et al. 2011). These techniques have resulted in the compilation of catalogs of several hundred to a few thousand Lyα emitting galaxies from redshifts z ∼ 2.1 to the redshifts associated with the end of reionization. These samples have been used to show that Lyα emitters of line luminosity Lα = 10 42 erg s −1 found at z = 3 have space densities of ∼ 10 −3 Mpc −3 (e.g., Gawiser et al. 2007, Cassata et al. 2015 and are therefore expected to be the progenitors of L * galaxies at redshift z = 0. The clustering of these galaxies has been measured on scales of up to 10 Mpc by Gauita et al. (2010) and Gawiser et al. (2007), who find that at redshifts z = 2 − 3 they have a bias factor with respect to the underlying matter (in CDM models) of b ∼ 1.5 − 2. Integrating the luminosity functions of Lyα emitting galaxies, assuming a power-law extrapolation for the faint end slope, has revealed that the comoving volume emissivity of Lyα photons declines significantly from z ∼ 6 to z ∼ 2 (Cassata et al. 2011, Gronwall et al. 2007, Ouchi et al. 2008. This behaviour can be compared to the opposite evolution in redshift of galaxies measured in optically thin parts of the rest-frame UV spectrum (or using the Hα line, e.g., Hayes et al. 2011). This comparison has been used to infer (e.g., by Hayes et al. 2011 andCassata et al. 2011) that the escape fraction of Lyα photons produced in star forming regions has significantly decreased from z = 6 to z = 2.
Because Lyα photons have a high cross-section for scattering off neutral hydrogen, extended Lyα emission is ex-pected to be common in many environments. For example, Lyα radiation from star forming regions in galaxies should undergo hundreds or thousands of scatterings in gas in any circumgalactic medium before finally escaping or else being absorbed by dust. The existence of a general fluorescent emission from the intergalactic medium was also hypothesised by Hogan & Weymann (1987) and Gould & Weinberg (1996). Theoretical work applying line radiative transfer on gas distribution in cosmological hydrodynamic simulations has made predictions for Lyα emission around galaxies and quasars (e.g., Cantalupo et al. 2005, Laursen et al. 2007, Kollmeier et al. 2010, Zheng et al. 2011a, as well as metal line emission (Bertone & Schaye 2012). These studies have resulted in predictions of extended Lyα halos around galaxies with sizes of hundreds of kpc, with a strong dependence of their properties on environment that can lead to new effects on galaxy clustering (Zheng et al. 2011a).
Observational evidence for extended emission includes the discovery and characterization of the so-called Lyα "blobs" (Steidel et al. 2000). Deep spectroscopic searches for diffuse Lyα emission have been completed by Rauch et al. (2008), finding faint Lyα emitting galaxies. Stacking of spectra of damped Lyα absorbers in quasars has also produced measurements of residual Lyα emission (Rahmani et al. 2010, Noterdaeme et al. 2014. Recently, Martin et al. (2014a,b) published the first results from the Cosmic Web Imager, an integral field spectrograph designed to map low surface brightness emission, detecting Lyα emission from filamentary structures around a z = 2.8 quasar as long as 250-400 proper kpc. Diffuse Lyα halos around high redshift galaxies have been found to be ubiquitous by Steidel et al. (2011) and Matsuda et al. (2012). Momose et al. (2014) have assembled several samples of up to 3600 Lyα emitters from Subaru narrowband imaging at a range of redshifts from z = 2.2 to z = 6.6 and, after controlling for atmospheric and instrumental artifacts, they detect diffuse extended Lyα halos with exponential scale lengths of ∼ 5 − 10 kpc from z = 2.2 − 5.7. The large scale studies in our paper are an alternative, complementary observational strategy to these earlier studies, which involve deep integrations over small fields of view.
All of these sources should be clustered on large scales and should contribute to the mean Lyα emission intensity in the Universe. This mean emission is detectable if it cross-correlates as expected with other tracers of largescale structure that we can observe at the same redshift. We shall use the quasars found by BOSS at z > 2 to correlate with Lyα emission in this work. This clustered Lyα emission is extremely faint, but as we shall demonstrate it can be detected with BOSS thanks to the enormous number of spectra that are observed. While large-scale clustering measurements cannot easily allow separation of the signal into various sources, we may expect faint Lyα emitting galaxies to dominate over quasars due to their much larger number density. If this is the case, then the mean Lyα emission intensity clustered with quasars can be used as a measure of the global star formation rate (e.g., Cassata et al. 2011), times the mean (luminosity weighted) bias factor of the distribution of these galaxies. Our measurement of Lyα emission will therefore be useful as a probe of star formation which takes into account all sources of Lyα emission, reaching to arbitrarily faint luminosities and surface brightnesses from extended halos to faint galaxies.
The structure of this paper is as follows. Section 2 describes the data samples we use in our work, which include the galaxy and sky spectra from SDSS DR10, along with the quasar catalog from that data release. We present our measurement of quasar-Lyα emission correlations in Section 3, including the evolution with redshift and clustering parallel and transverse to the line of sight. Section 4 describes our tests involving fitting and subtraction of emission lines. In Section 5, we convert our determination of the Lyα surface brightness into a star formation rate density and compare to other measurements. In Section 6, we summarise our results and in Section 7 discuss them further. There are also three appendices to the paper, A-C; in these we measure stray light contamination, determine a large-scale surface brightness correction, and perform some sample tests of our results.

DATA SAMPLES
This study makes use of data from the SDSS BOSS survey Data Release 10 (DR10, Ahn et al. 2014), including quasar position and redshift data, galaxy spectra and sky fiber spectra. The SDSS camera and telescope are described in Gunn et al. (1998) and Gunn et al. (2006), respectively. Full information on the SDSS/BOSS spectrographs can be found in Smee et al. (2013). The wavelength coverage of the spectrograph is from λ =3560Å to 10400Å the resolving power is R ∼ 1400 for the range λ = 3800Å− 4900Å, and is kept above R = 1000 for the remainder of the wavelength range. The fibers have a diameter of 120 µm, corresponding to 2 arcsec. in angle. We restrict the redshift range of data we use in our analysis to 2.0 < z < 3.5, due to the spectrograph cutoff at low redshift and the limited number of observed quasars at high redshift.

Spectra
The 987482 galaxy spectra in our sample are of targeted LRGs which are within redshifts z ∼ 0.15 and z ∼ 0.7. The redshift range of the original targets is not important to our study, as for each spectrum we make use only of the pixels for which the Lyα emission line lies within the redshift range specified above ( 2.0 < z < 3.5). In observed wavelength units this is 3647Å to 5470Å. We also make use of 146065 sky fiber spectra.
The main BOSS LRG program consists of two galaxy target samples (see Dawson et al. 2013), designated CMASS (for constant mass) and LOWZ (for low-redshift). The LOWZ galaxy sample is composed of massive red galaxies spanning the redshift range 0.15 < ∼ z < ∼ 0.4. The CMASS galaxy sample is composed of massive galaxies spanning the redshift range 0.4 < ∼ z < ∼ 0.7. Both samples are color-selected to provide near-uniform sampling over the combined volume. The faintest galaxies are at r = 19.5 for LOWZ and i = 19.9 for CMASS. Colors and magnitudes for the galaxy selection cuts are corrected for Galactic extinction using Schlegel et al. (1998) dust maps. We do not differentiate between CMASS and LOWZ samples in our analysis.
The spectroscopic measurement pipeline for BOSS is described in detail in Bolton et al. (2012). The most important data products that are used in the present analysis are: (a) Wavelength-calibrated, sky-subtracted, fluxcalibrated, and co-added object spectra, which have been rebinned onto a uniform baseline of ∆ log 10 λ = 10 −4 (about 69 km s −1 pixel −1 ). (b) Statistical error-estimate vectors for each spectrum (expressed as inverse variance) incorporating contributions from photon noise, CCD read noise, and skysubtraction error. (c) Mask vectors for each spectrum.

Data preparation
For each of the LRG spectra, we subtract the best fit model spectrum provided by the pipeline. This template model spectrum (see Bolton et al. 2012 for details) is computed using least-squares minimization comparison of each galaxy spectrum to a full range of galaxy templates. A range of redshifts is explored, with trial redshifts spaced every pixel. At each redshift the spectrum is fit with an error-weighted leastsquares linear combination of redshifted template eigenspectra in combination with a low-order polynomial. The polynomial terms absorb Galactic extinction, intrinsic extinction, and residual spectrophotometric calibration errors (typically at the 10% level) that are not fully spanned by the eigenspectra; there are three polynomial degrees of freedom for galaxies. The template basis sets are derived from restframe principal-component analyses (PCA) of training samples of galaxies, and have four degrees of freedom (eigenspectra). After subtraction of the best fittting template spectrum, we compute the average residual spectrum of all galaxies. This is displayed in Figure 1 where the horizontal axis is labelled in units of the redshift of the Lyα line. Figure 1 also presents the mean sky spectrum and the sky-subtracted sky spectrum.
We can see that the residual surface brightness per unit wavelength (hereafter shortened to "surface brightness"-we use this term to refer to the quantity measured throughout the paper, which is most precisely the flux density per unit solid angle per unit wavelength) in the galaxy fibers is within ±10 −19 erg s −1 cm −2Å −1 arcsec −2 for most of the redshift range. There are however significant excursions corresponding to features including the zero redshift Calcium H and K lines (at 3969 and 3934Å) and a strong Mercury G line from streetlamps (at 4358Å). In our analysis we subtract the mean residual surface brightness, from all spectra before cross-correlating them as we are only interested in the fluctuations in the Lyα surface brightness. In order to reduce noise, we also mask two regions corresponding to large features in the residual surface brightness, 40Å and 30Å windows centered on wavelengths 3900Å and 4357Å respectively (corresponding to redshifts z = 2.21 and 2.58).
Comparing the sky fiber and galaxy fiber residual spectra, we can see that there are differences at the ∼ 10 −19 erg s −1 level over much of the spectra. We attribute these to galaxy surface brightness that was not subtracted perfectly by the galaxy model. In our cross-correlation technique for measuring clustering in the Lyα emission we necessarily subtract the mean surface brightness, therefore residual fluctuations seen in Figure 1 are not problematic except for the noise they contribute. : the average sky fiber surface brightness in all 146065 sky fiber spectra. Cyan line: the average sky fiber surface brightness with model sky subtracted from all sky fiber spectra. In all curves, the prominent emission line at wavelength λ = 4358Åis due to terrestrial airglow from Mercury streetlamps (the Hg G-line).

Quasars
We use quasars from the SDSS/BOSS DR10 catalogue (Ahn et al. 2014). The quasar target sample included both colorselected candidates and known quasars (Bovy et al. 2011;Kirkpatrick et al. 2011;Ross et al. 2012). The candidate quasar spectra were all visually inspected and redshift estimates computed using a principal component analysis (see Pâris et al. 2012 for the details of the procedure as applied to DR9 quasars). We select the 130812 quasars in the DR10 dataset that have redshifts in the range 2.0 < z < 3.5. Because the galaxy pixels cover this redshift range uniformly, the central redshift of our measurements is the mean redshift of these quasars, z = 2.55.

QUASAR-Lyα EMISSION CROSS-CORRELATION
Before computing the quasar-Lyα emission crosscorrelation, we first split the sample of galaxy spectra into 100 subsamples of approximately equal sky area based on contiguous groupings of plates. We then convert the galaxy spectrum pixels and the quasar angular positions and redshifts into comoving Cartesian coordinates using a flat cosmological model with matter density Ωm = 0.315, consistent with the Planck, Ade (2014) results (cosmological constant density ΩΛ = 0.685). This fiducial model is used throughout the paper. We compute the quasar-Lyα emission surface brightness cross-correlation, ξqα(r) , using a sum over all quasar-galaxy spectrum pixel pairs separated by r within a certain bin: where N (r) is the number of pixels in the bin centered on quasar-pixel distance r, and ∆µ,ri = µri − µ(z) is the residual surface brightness in the spectrum at pixel i for the bin r. Note here that we have a different list of pixels labeled as i for each bin in the separation r between a pixel and a quasar, which has Lyα surface brightness µri. The residual flux at each pixel is obtained by subtracting the mean at each redshift, µ(z) . We weight each pixel by wri = 1/σ 2 ri , where σ 2 ri is the pipeline estimate of the inverse variance of the flux at each pixel. We first present our results as a function of only the modulus of the quasar-pixel separation r in comoving h −1 Mpc, in 20 bins logarithmically spaced between r = 0.5h −1 Mpc and r = 150h −1 Mpc. In Section 3.3 we will also examine redshift space anisotropies in the correlation function ξqα by considering bins in the parallel and perpendicular components of r, using the same formulation.
When evaluating equation 1, a possible signficant systematic error is caused by stray light from the quasars themselves contaminating spectra of nearby galaxies. This occurs because the light from the various fibers is dispersed onto a single CCD, so that extraction of each spectrum along one dimension (Bolton et al. 2012) may include light from adjacent fibers of bright sources. We see strong evidence of this stray light from quasars in galaxy spectra when the quasar and galaxy spectra are four fibers apart or fewer, in the list of fibers as they are ordered in the CCD. The effect is discussed in detail in Appendix A. When the galaxy and quasar spectra in a pair are more than 4 fibers apart, we see no evidence for this contamination, and the results are statistically consistent with using only pairs of quasars and galaxies on different plates (see Appendix A). In order to safely eliminate this stray quasar light when computing the flux cross-correlation with equation 1, we therefore apply the constraint that the quasar and galaxy fibers must be at least six fibers apart.
There is also the possibility that some clustering in the plane of the sky is generated by effects (e.g., galactic obscuration) which modulate both Lyα surface brightness and quasar target selection. Appendix B presents measurements of ξqα for quasar-pixel pairs which are close together on the sky (i.e., in the transverse separation) but widely separated along the line of sight. This measurement enables us to quantify how much clustering could be caused by effects such as Galactic obscuration and to compute a ξqα correction term to be subtracted from our fiducial clustering result. We also measure ξqα for pairs which are close in the line-of-sight separation but widely separated on the sky. This latter measurement constrains how much spurious clustering is caused by large-scale variations in the line-ofsight direction, for example redshift evolution in the efficiency of galaxy subtraction, or flux calibration errors with wavelength that may be associated with sky lines. We apply the corrections to ξqα from Appendix B to our analysis below and in the other sections of the paper. We discuss the small-scale anisotropy of ξqα in Section 3.3 below. For now we note that application of the correction factors described above changes the amplitude and shape of our measured ξqα by less than 1σ in all cases.
We perform the pairwise computation of equation 1 for each of our 100 subsamples, and then compute the mean and standard deviation of ξqα(r) using a Jackknife estimator. The Jackknife estimator is also used to compute the covariance matrix of ξqα(r) : where ξ qα,k (ri) is the cross-correlation in bin i for Jackknife sample k, ξqα(ri) is the cross-correlation for bin i for the full dataset, and the number of Jackknife samples is M = 100.

Fiducial result
We show ξqα(r) for our fiducial sample (which is the entire dataset over the redshift range 2.0 < z < 3.5) in Figure  2. The mean Lyα redshift of the galaxy pixels in this sample is z = 2.71 and of the quasars z = 2.55. Because the galaxy pixel distribution is uniform in redshift, and because quasar-pixel pairs with small separations contribute most to the clustering signal, we adopt the effective mean redshift of our fiducial measurement to be z = 2.55. The crosscorrelation function is in units of the surface brightness of Lyα emission, and its amplitude is directly proportional to that surface brightness. The ξqα(r) points reveal that there is significant measurable large-scale structure present in the Lyα emission, on scales from 1 to ∼ 15h −1 Mpc. Figure 2 also displays a linear CDM fit to the cross-correlation function, which is consistent with observational results on scales from 1 − 100h −1 Mpc. We turn to this fit in the next subsection.

Model fit
If the Lyα emission clustering is due to a linearly biased version of the density field, then a model for the isotropically averaged quasar-Lyα cross-correlation ξqα(r) is as follows: where µα is the mean surface brightness of Lyα emission, bq and bα are the quasar and Lyα emission linear bias factors, and ξ(r) is the linear ΛCDM correlation function. It is important to note that bα is the bias factor for Lyα surface brightness fluctuations, and is different in definition from the usually quoted bias factor of Lyα emitters, bLAE (e.g., as measured by Gawiser et al. 2007 andGuaita et al. 2010). The bias factor bLAE reflects the relation between the fluctuations δn in the number density n of Lyα emitters and that in the matter density δ, where δ = (ρ − ρ )/ ρ and ρ is the matter density field. The factor bα in Equation 3 relates fluctuations δµ in the Lyα surface brightness µ to matter fluctuations according to In the absence of radiative transfer effect (Zheng et al.  Equation 1). The points represent results for the fiducial sample that covers redshift range 2.0 < z < 3.5. The error bars have been calculated using a jackknife estimator and 100 subsamples of the data. The smooth curve is a best fit linear CDM correlation function (see Section 3.2). The top panel shows the ξqα(r) results with a log y-axis scale, and the bottom panel displays rξqα(r) on a linear scale in order to allow points which are negative to be visible.
2011a; see Section 6.4), the Lyα surface brightness µ is proportional to the Lyα luminosity density ρL of the underlying star-forming galaxy population. The fluctuations δL of the latter can be characterized by the bias factor bL, and we have bα = bL. As bL reflects weighting by luminosity rather than by number, it is likely to be significantly higher than bLAE, because higher luminosity emitters tend to be more strongly clustered. We will return to this topic in Section 5. Radiative transfer effect leads to a modification in the relation bα = bL. In a simple model, we have bα = bL + α1 with α1 a positive number (see Section 6.4).
Overall, we expect bα to be substantially higher than bLAE.
In Equation 3, f β is a constant enhancement to the correlation function on linear scales of the form that is caused by peculiar velocity redshift-space distortions (Kaiser 1987). We use the linear CDM transfer function of Lewis et al. (2000) to compute ξ(r). In our computations we choose to vary the shape of the correlation function by changing Ωm, the matter density, in the context of the ΛCDM model, keep-ing the other parameters which influence the shape (such as h, and the baryon density Ω b ) fixed. Parameters have been reported for the best fit ΛCDM model to the Planck Satellite data by Ade et al. 2014. We assume Planck values for Ω b = 0.049 and the spectral index ns = 0.9603 but set h = 0.7. Note that we are merely using Ωm to parametrize the shape of the correlation function to see if it is consistent with other observations, and are not presenting our results for Ωm as properly marginalized measurements of that parameter.
The other free parameter is the amplitude, bqbαf β µ . We assume in all cases that the underlying amplitude of mass fluctuations σ8(z = 0) = 0.83, and therefore that σ8(z = 2.55) = 0.294, again consistent with Ade et al. 2014.
We fit our model to the data in Figure 2 by varying these two parameters. The χ 2 value is given by where the sum is over the N = 20 bins, ξ obs qα (ri) is the observed cross-correlation measured in bin i, ξ mod qα (ri) is the model prediction for bin i, and Cij is the covariance matrix computed using our 100 Jackknife samples in Equation 2.
The best fit values and one sigma error bars are as follows: and Ωm = 0.296 +0.103 −0.071 . The shape parameter Ωm is consistent with the best fit value from the Planck satellite results (Ade et al. 2014, Ωm = 0.30). The 1, 2 and 3 σ confidence contours in these parameters considered together (i.e. ∆χ 2 =2.3, 6.17 and 11.8) are displayed along with the best fit values in Figure 3.
When carrying out our linear fit to the isotropically averaged correlation function we are neglecting non-linear effects and redshift measurement errors which are expected to suppress the correlation function on small scales (Davis and Peebles 1983). We are also not including in our model of Equation 3 the one-halo term (see e.g., Cooray and Sheth 2002) representing structure inside virialized halos, which would cause the correlation function to be increased at small scales. An estimate of the size of the redshift-space distortion effects is given below. There is some evidence for suppression of redshift-space distortions of the form that is expected from radiative transfer effects (Zheng et al. 2011a). The effect of the suppression terms in the cross-correlation function will be of opposite sign to the boost due to a onehalo term, so that the partial cancellation will further lower the impact of these effects. Because of these interactions, we leave more detailed nonlinear and halo modelling of the correlation function to future work.

Clustering transverse and parallel to line of sight
We also compute the quasar-Lyα cross-correlation as a function of r and r ⊥ , the quasar-pixel pair separation along and across the line of sight, shown in Figure 4, on a linear scale.
We can see that the contours are relatively symmetric about the r = 0 axis and somewhat stretched along the r direction. Font-Ribera et al. (2013) found a redshift offset between quasars and the Lyα forest of δz = −160 km s −1 , due to the quasar catalogue redshifts being on average too small by this amount. This correlation resulted in the quasar-Lyα forest cross-correlation being shifted upwards by this amount. The precision of the quasar-Lyα emission crosscorrelation in our paper is smaller, but visual inspection of Figure 4 reveals that the position of the centroid of the cross-correlation is consistent with a small upward shift of this magnitude (1.6h −1 Mpc at this redshift).
On the scales where the cross-correlation is easy to discern (r < ∼ 20h −1 Mpc), there is no sign of compression due to linear infall (Kaiser 1987). The redshift-space distortions should be less prominent compared to the Lyα forest because of the expected higher bias for the Lyα emission. In reality there appears to be stretching along the line of sight (we quantify this below), which might be due to a combination of quasar redshift errors, the intrinsic velocity dispersion of quasars in their host halos, or the intrinsic velocity dispersion of the sources of Lyα emission.
Another source of apparent clustering anisotropy could be the radiative transfer effects predicted by Zheng et al. (2011a). It was shown by these authors, using cosmological radiation hydrodynamic simulations, that Lyα radiative transfer has a strong environmental dependence which can cause the apparent spatial distribution of Lyα emission to become anisotropic with respect to the line-of-sight direction. Density fluctuations along the line-of-sight direction are found to preferentially emit the Lyα radiation in that direction in overdense regions, mainly because of the effect of peculiar velocity gradients on the Lyα radiative transfer. This causes a suppression of the line-of-sight fluctuation, which can be modeled similarly to the Kaiser effect (also caused by the peculiar velocity gradient), even though the sign of the effect is opposite.

Fitting redshift-space distortions
In order to approximately quantify the level of distortion in Figure 4 and its statistical significance, we have investigated fitting a redshift-space distortion model to the ξqα (r ⊥ , r ) data. To compute the model for ξqα (r ⊥ , r ), we first assume the linear ΛCDM correlation function shape used in Equation 3 and then use a model for peculiar velocities to distort it in redshift space. Our pecular velocity model includes standard linear infall for large scale flows (Kaiser 1987) and a small scale random velocity dispersion (e.g., Davis and Peebles 1983).
The parameterization of the model for linear infall allows for stretching (outflow) as well as squashing (infall) along the line of sight. Although gravitational processeses are expected to result in infall, as mentioned above Zheng et al. (2011b) have shown that radiative transfer effects on the anisotropy of clustering can be approximately parametrized with the same model. We do this here, allowing a net linear outflow measured by the model to be interpreted as the radiative transfer effect.
The effects of coherent flows on the correlation function in linear theory were presented by Hamilton (1992). We use the formulation of Hawkins (2003), with modifications to make it appropriate for the case of cross-correlation functions. This modification involves the use of the two bias factors bq and bα from Equation 3 to compute redshiftspace distortion factors βq = Ωm(z = 2.55) 0.6 /bq and βα = Ωm(z = 2.55) 0.6 /bα. The linearly distorted quasar-Lyα emission cross-correlation function is then given by where µ = r /r, and ξ0(s) = 1 + 1 3 (βq + βα) + 1 5 βqβα ξ(r) (10) with Here ξ(r) is the linear ΛCDM correlation function of equation (3). We use these relations to create a model ξ ′ qα (r ⊥ , r ) which we convolve with the distribution function of random pairwise motions, f (v), to produce the final model ξqα(r ⊥ , r ): The random velocity dispersion we use is an exponential where σ is the pairwise velocity dispersion of quasars and Lyα emission, which we assume to be independent of pair separation.
There is a strong degeneracy between βq and βα, but the bias factor for BOSS quasars is reasonably well measured. The first BOSS measurement, of bq = 3.6 ± 0.6, was made using the quasar autocorrelation function by White et al. (2012). Font-Ribera et al. (2013), find an even more precise value of bq = 3.64 +0.13 −0.15 from the cross-correlation of BOSS quasars with the Lyα forest. Because of this we set bq = 3.64, giving βq = 0.27. The free parameters in our distortion model are therefore βα, σ and an amplitude parameter bqbα µα .
We set the parameter governing the shape in our model, Ωm, equal to the Planck value, Ωm = 0.30. We compute ξqα (r ⊥ , r ) for our model for a grid of values of varying βα, σ, and bqbα µ . We then compare our model to the observed ξqα (r ⊥ , r ) from Figure 4, performing a χ 2 fit to all points within r = 40h −1 Mpc. We again use jackknife error bars computed from 100 subsamples, but because of difficulties with noisy matrix inversion, we do not use the offdiagonal (C i =j ) terms when inverting the covariance matrix. We marginalize over the amplitude parameter bqbα µα , and show the confidence contours (for the two remaining degrees of freedom) in Figure 5.
We have allowed the βα parameter to be negative in our fit not because we believe that outflow of Lyα emission is likely around quasars but, as mentioned above, because this allows quantifying the stretching along the line of sight seen in Figure 4 and its significance. Our best fit values with 1 σ error bars are βα = −0.76 ± 0.36 and σ = 490 ± 300 km s −1 . From these values and by observing the contours in Figure   5, we infer a detection of anisotropies in the quasar-Lyα emission correlation function that is opposite in sign to that expected from peculiar velocity flows due to gravitational evolution: our constraint on βα is 2.1σ from βα = 0. We interpret this result as indicating that there are strong nongravitational effects on the Lyα emission causing the elongation of the cross-correlation contours along the line of sight extending to large separations; the good fit obtained with the redshift-space distortion model with negative β needs to be understood in this case as a coincidence, since the model is not physically correct.
If we limit the fit to points with 20 h −1 Mpc, we find central values for βα and σ that are consistent with those for our fiducial fit, As might be expected, the error bars are larger, however, by a factor of 40% for σ and 65% for βα.
If we assume a Lyα emission bias factor bα = 3, corresponding to highly biased star forming galaxies, (see Section 5 for further discussion of this value), then βα = 0.32 at redshift z = 2.55. For the expected value of σ, we can use as a guide the results of Font-Ribera et al. (2013) who constrained σ < 370 km s −1 at the 1 σ confidence level from the redshift-space quasar-Lyα forest cross-correlation function. We can see from Figure 5 that although the σ measurement is consistent with Font-Ribera et al. , βα and σ considered jointly disagree at the 2.5σ level from our measurement.
We discuss further in Section 6 how our measurements can be interpreted as a detection of clustering anisotropies due to radiative transfer effects (Zheng et al. 2011b). In this case, elongation along the line of sight is expected, which can explain the effective measurement of a negative βα parameter.
The cross-correlation contours of the best fit model are plotted in Figure 4 (right panel), which also reveals stretching along the line-of-sight. The χ 2 value for the fit is 610 measured from 400 bins, with 3 free parameters, a reduced χ 2 of 1.5. The fit is therefore not good, and the discrepancy arises in large part in the central region, where the model has lower surface brightness than the observations. This result may be a sign that adding a one-halo term to the correlation function would provide a better fit, and should be addressed in future work with a larger data sample.

Large-scale tests
The cross-correlation across and along the line of sight over the whole spectrum offers a way to test whether the detected signal is reasonable. One can search for any significant crosscorrelation signal if a different wavelength other than Lyα is used, which would indicate either contributions from other emission lines, or that other effects are causing the crosscorrelation signal meaasured. An equivalent approach is to extend the line-of-sight range of the cross-correlation to large distances. Figure 6, extends the contours shown in Figure 4 to much larger scales. The positive signal seen in Figure 4 is the most significant feature, centered at r = 0, r ⊥ = 0. This is a good indication that Lyα emission is the dominant contribution to our signal. Signal from lines at longer (shorter) wavelengths than Lyα would appear at positive (negative) values of r .
The second most prominent feature, at r ⊥ = 0, r ∼ 60h −1 Mpc, is significant at the ∼ 1.5σ level (the pixel at the center of the feature is 1.5 σ from the zero level, using  the Jackknife error bars from Section 3.4), and so is consistent with noise. Strong lines that might be an issue, such as Lyman-β and Carbon-IV, are very far away (at r = −490 and +640 h −1 Mpc respectively) and so are not a concern. Si-III would appear at r = −22h −1 Mpc if it was present.

Evolution with redshift
The redshift coverage of our data sample (z = 2−3.5) is sufficient that we can separate it into different bins in redshift and search for evolution. We do this for four different redshift bins where each bin contains one quarter of the quasar data. The bin boundaries and mean quasar redshifts of each bin are given in Table 1.
The quasar-Lyα emission cross-correlation results are shown in the four panels of Figure 7. The global CDM model fit to the full sample averaged over all redshifts is indicated as the dash-dot line in every panel. Although the results for the redshift bins are relatively noisy, as expected, they are broadly similar to the global result in shape and amplitude and show no clear evidence for any redshift evolution. Figure 8 presents the best fit CDM shape and amplitude parameters from Section 3.2, along with confidence contours. The fitting was again performed using jackknife error bars. Figure 8 demonstrates that the fiducial model results all lie within the 2σ confidence contours for the different redshift bins. Assuming that the shape of the cross-correlation function remains fixed in comoving coordinates (as it would do if governed by linear biasing), we can search for changes in the amplitude of clustering and the mean Lyα surface brightness as a function of redshift. We set Ωm = 0.30 (the CDM shape determined by the Planck results; Ade et al. 2014) and then determine the best fit amplitude parameters, bqbαf β µ , at each redshift bin. The result, shown in Figure 9, indicates a decreasing cross-correlation amplitude with redshift, although the errors are large and a horizontal line would not be an unreasonable fit "by eye". To express this quantitively, we have carried out a χ 2 fit to the function log (bqbαf β µ )= a + bz, finding a slope parameter b = −0.40 ± 0.20, meaning that the hint of redshift evolution is significant at the 2.0σ level only. The values for the fit parameters bqbαf β µ for the different redshift bins (which were used to plot Figure 8 are listed in Table 1. We have also looked at redshift evolution of the CDM shape governed by the parameter Ωm. Within the assumption of linear biasing, the shape should remain constant with redshift. The results are examined in Figure 10, where the results are indeed consistent with a constant Ωm, within the uncertainties. The Planck value (Ωm = 0.30) is also shown and is consistent with our results.

EMISSION LINES
Given that the sources of the clustered Lyα emission seen in Section 3 could be discrete objects such as star forming galaxies, we must investigate see if individual emission lines can be detected in our spectra and whether discrete detectable lines can account for the signal. The BOSS spectra have relatively short integrations on a small (2.5 m) telescope, and cannot be expected to compete in individual detections with other surveys such as that described in Cassata et al. (2011): we can only detect the most luminous objects. However, our cross-correlation technique enables us to find the mean total surface brightness, which includes all emission-line objects no matter how faint they may be. The difference between the cross-correlation signal with and without individually identified lines therefore enables us to discover what fraction of the Lyα surface brightness lies below our line detection limit.
We note that our line detection procedure is less sophisticated than that in the BELLS survey (Brownstein et al. 2012), which used line detections of galaxies behind LRGs to find gravitational lenses. In particular, we are not seeking confirmed detections of objects (which requires multiple emission lines) and we do not deal with interlopers, except statistically.

Line fitting
For each LRG spectrum (see Section 2.1), we subtract the best fit galaxy spectrum model, as we do in our fiducial cross-correlation analysis. We then fit lines to this residual Figure 7. The quasar-Lyα emission cross-correlation function ξqα(r) for different redshift ranges, which correspond to one quarter of the full dataset each. The solid lines indicate the best fitting CDM model at each redshift bin and the dash-dot line represents the best fit to the full dataset (shown in Figure 3). Both the shape and amplitude of the model fit are consistent at the 2σ level with no evolution over the redshift range (this is shown in Figures 9 and 10). flux, centering our fitted line profile on the center of each spectrum pixel, one at a time. In this first stage, each spectrum pixel is therefore the center of a best fit line profile -we remove overlapping lines later. For each of a grid of values of amplitude A and line width σ (between 0 and 20Å) we compute the χ 2 value of a positive Gaussian emission line profile G(x), where the profile has the form G(x) = A exp −(x 2 /2σ 2 ) and x is the separation between the line center and the pixels we include in our fit. We use pixels in our fit that are in a 40Å region centered on the line center, excluding masked regions as in Section 2.1.1. Once the best fitting values of A and σ are found, we estimate the significance of each fitted line from the χ 2 difference between the line fit and a flat interval with zero flux. After fitting to all the pixels we pass through the list and eliminate overlapping lines, removing the lower significance line when there is an overlap.
We find 5200 lines with a nominal significance of 5σ (∆χ 2 > 30.1 for the 2 degrees of freedom fitted), and 1.6×10 6 lines with a nominal significance of 3σ (∆χ 2 > 11.8) There are 1.3 × 10 9 pixels in the search regions of the spectra. The detected lines are constrained to be at least 40Å (37 pixels) apart, but for rare lines this should not change . Likelihood contours for power spectrum parameters Ωm (which we are using to parametrize the CDM shape) and bqbαf β µ (the amplitude) for different redshift ranges. The filled dot shows the best fitting values for that redshift range and the open circles the best fitting values for the full dataset (shown in Figure 3). the random expectation. One would therefore expect to find approximately 370 and 1.8×10 6 5σ and 3σ lines respectively from positive noise fluctuations alone.
For the 3σ lines, the fact that we have detected fewer lines than even pure noise fluctuations predict is likely to be a sign that the fluctuations do not exactly obey a Gaussian noise model. An additional complication is that the noise estimate from the standard data pipeline which we have used has been shown to be underestimated by up to 16% for the relevant wavelengths (Palanque-Delabrouille et al. 2013). The detection of more 5σ lines than randomly expected is likely to indicate that there are false detections arising from unsubtracted features in the galaxy spectra, sky lines we have not accounted for, and other systematics. There is also the possibility of interloping [OII] emission lines from lower redshift galaxies (see Noterdaeme et al. 2010, Menard et al. 2011. For our wavelength coverage of 3800-5500Å, this interloper emission would arise from between z = 0.02 and z = 0.48. Without a significant additional effort to remove false detections and interlopers, our dataset is not useful for computing the luminosity function of Lyα lines. Instead, we turn to the statistical cross-correlation to test for the fraction of these lines which are really Lyα emission lines in our redshift interval.

Cross-correlation
We subtract the flux in the lines detected in Section 4.1 from each LRG spectrum, and then recalculate the crosscorrelation of quasars and Lyα emission (Equation 1). The results for ξqα(r) (again computing the error bars using a jackknife estimator and 100 subsamples) are shown in Figure 11, using our two thresholds on the significance of the removed lines, 5σ and 3σ. We can see that in each case, the clustering signal is still visible and the shape traces that of a CDM curve. This shows that most of the surface brightness of Lyα emission is not accounted for by these lines. As expected, most of the lines are due to noise features. By subtracting these lines, however, we are also subtracting any possible real lines, and the change in the amplitude of ξqα(r) is a measure of the fraction of surface brightness that is actually contributed by strong lines.
The shape and amplitude fitting parameters (Ωm and bqbαf β µ ) for these two cases (> 5σ and > 3σ lines subtracted) are shown in Figure 12. The shape parameter Ωm is very similar in the two cases and almost the same as in the fiducial case. The amplitude is lower, as would be expected for subtraction of some real lines, but the fiducial result (with no line subtraction) lies well within the 1 σ error contour of both of the panels in Figure 12, implying that the contribution to the cross-correlation from emission lines that are detected and removed is not statistically significant. Quantitatively, this can be seen by considering that we find the amplitude parameter for the > 5σ case to be bqbαf β µ = 3.18 +0.39 −0.41 × 10 −20 erg s −1 cm −2Å−1 arcsec −2 , and for the 3σ case to be bqbαf β µ = 2.89 +0.43 −0.37 × 10 −20 erg s −1 cm −2 A −1 arcsec −2 . The amplitude parameter is therefore 4±12% and 13 +11 −13 % lower than the fiducial case for the > 5σ and > 3σ line removal cases, but both of these are consistent with zero within the errors. The analysis is therefore consistent with our line fitting having found no true Lyα emission lines at all.
In order to relate the significance levels to line luminosity, we have computed the luminosity from the surface brightness for each line (bearing in mind that our measurements are restricted to a 1 arcsec radius fiber aperture). We find that the median luminosity of the > 5σ lines is L = 9.0 × 10 42 erg s −1 and the > 3σ lines have a median luminosity L = 1.9 × 10 42 erg s −1 . We can compare these luminosities measured with some published values from Lyα emitter surveys. The flux limit of the Guaita et al. (2010) data sample was 2 × 10 −17 erg s −1 cm −2 (emission line flux) at z=2.1. This corresponds to a Lyα luminosity of 5 × 10 41 erg s −1 . For the Gawiser et al. (2007) sample at z=3.1, the line flux limit was 1.5 × 10 −17 erg s −1 cm −2 , corresponding to a line luminosity of 1.3 × 10 42 erg s −1 .
In our calculation above we have seen that at the 1 σ confidence level, 13% of the Lyα cross-correlation signal in our dataset could be due to lines with a median luminosity (in a 1 arcsec radius aperture) 1.9 × 10 42 ergs −1 . This is similar to the values of Gawiser et al. and Guaita et al. , although given our small aperture, the intrinsic luminosity of our fitted line emitters will be even higher. A small fraction of the bqbαf β µ value we are seeing could therefore be contributed by emitters similar to those in these two surveys. As noted above, the error bar on this fraction is large, and our result is consistent with zero contribution from such lines.
There are various possibilities for the nature of the ma- Figure 11. The quasar-Lyα emission cross-correlation function ξqα(r) (see Equation 1), as in Figure 2, but computed after subtracting emission lines that are apparently detected in the spectra at the 5σ significance level in panel (a), and the 3σ level in panel (b). The smooth curve is the best fit linear CDM correlation function (see Section 3.2) and the dash-dotted line is the best fit CDM curve for the fiducial sample (i.e., before subtracting the apparent emission lines).
jor contributor to the Lyα cross-correlation signal. The first is fainter lines than those seen in Lyα emitter surveys. This is unlikely because bα is luminosity weighted, so low luminosity lines contribute much less to the signal than higher luminosity lines. The second is high luminosity, low surface brightness emission. This emission is much more difficult to detect, and cannot be seen by our line search algorithm, which is sensitive to high surface brightness, narrow lines. It is also unlikely to have been seen in previous surveys (we return to this in Section 6). This type of emission would also be highly biased, and so if present would contribute strongly to the bqbαf β µ measurement. We discuss the various possibilities in more detail in Section 6.

STAR FORMATION RATE DENSITY
If the Lyα surface brightness we are seeing is produced by star forming galaxies, we can convert the mean Lyα surface brightness into a star formation rate density (SFRD) at the mean redshift z = 2.55 of our observations. Traditionally (e.g., Gronwall et al. 2007) narrow band surveys have been used to detect Lyα emitters, compute their luminosity function, and integrate it to compute a mean Lyα luminosity density ǫα before converting it into a star formation rate, using a relationship such as: (Cassata et al. 2011), where Lα is the Lyα luminosity. The conversion factor is based on a stellar population with a Salpeter IMF and with no correction for effects like dust and escape fraction, which is accurate to within a factor of a few for a range of population age, high mass cutoff of stars, and metallicity (Leitherer et al. 1999). This method assumes that the surveys of Lyα emitters Figure 13. The star formation rate density (ρ SFR ) inferred from our measurement of the mean Lyα surface brightness in the Universe between z = 2−3.5 (see Section 5) is shown as the red point with solid line error bars, assuming that the linear bias factor for Lyα emission is bα = 3, a reasonable value for the luminosityweighted clustering of star forming forming galaxies (see Section 5). The true value of bα is unknown, so this data point should be scaled by 3/bα. Other data values plotted with open (black) symbols are from published ρ SFR values which used UV estimators. The solid (blue) points show estimates of ρ SFR computed from the luminosity functions of surveys for Lyα emitters. The references are given in Section 5. The shaded area represents the range of dust corrected UV estimates compiled by Bouwens et al. (2010).
are able to capture all the radiation from young stars. However, these surveys can only detect the high surface brightness portion of sources within a small angular aperture, and may be missing much of the Lyα line intensity when it is scattered far out into the galaxy halo. In our case, the statistical cross-correlation technique we are using should not be affected by any threshold in Lyα surface brightness. We should therefore be able to compute the total star formation rate density from our measurement. One large uncertainty is absorption due to dust, which is known to significantly affect UV continuum and line estimators of star formation. Recall that our measurement is of the quantity bqbαf β µ , so to compute the Lyα surface brightness we need to have independent knowledge of bq, bα and f β . For the quasar bias factor we use the value measured for BOSS quasars by Font-Ribera et al. (2013), bq = 3.64 +0.13 −0.15 . We recall that the bias factor bα is related to a luminosity-weighted bias factor bL, from the definitions in Equations 5-6 (bα = bL in the absence of radiative transfer effect), and bL is different from the number weighted bias factor of Lyα emitters, bLAE. To understand the difference in the values of the two bias factors, we start with the following simple model for the Lyα emission.
If we assume that there are N (M ) galaxies per dark matter halo of mass M and that Lyα emission comes from galaxies in halos above a mass limit Mmin, then the spatial bias can be computed as follows (e.g. Berlind & Weinberg 2002): where b h (M ) is the bias factor for halos of mass M and dn/dM is the halo mass function. For the luminosity weighted bias factor, we have where L(M ) is the average Lyα luminosity in halos of mass M .
To estimate the luminosity-weighted bias factor bL, we need to know the relation between luminosity and halo mass. Assuming that L(M ) ∝ M p , for Mmin = 10 11 M⊙ we obtain bL = 2.38, 4.40, and 6.84 for the values p = 1, 2, and 3, respectively. The differences between bLAE and bL are therefore substantial, the latter being usually much larger. To proceed further, we need to consider the likely relation between L and M . The Lyα luminosity is related to the star formation rate in galaxies. Along the star-forming sequence, the star formation rate is inferred to be approximately proportional to the stellar mass (e.g. Daddi et al. 2007;Pannella et al. 2009;Lee et al. 2011). We therefore first make use of the relationship between halo mass and stellar mass found with abundance matching by Moster et al. (2010). It is in the form of a softened broken power law, and at z ∼ 2.55 the low-mass-end (high-mass-end) slope of the L-M relation is ∼ 2.5 (∼ 0.6) with a transition mass around 10 12 M⊙. Using this relation in Equation 19 yields bL = 2.82 for Mmin = 10 11 M⊙. In fact, the result is insensitive to Mmin (bL = 2.81 for Mmin = 10 9 M⊙), given the steep L-M relation (so low-mass halos are weighted less). If we modify the high-mass-end slope to ∼ 1 to approximately account for the luminosity contribution from satellite galaxies, we obtain bL = 3.15.
There are various uncertainties involved in deriving bL with this simple model. First, the slope of the star formation rate versus the stellar mass relation for the star-forming galaxy sequence can be slightly different from unity. We find that a 10% deviation from unity in the above slope leads to a ∼ 5% change in the value of bL. Second, we assume that the Lyα luminosity is proportional to the star formation rate. The way these two quantities track each other may vary as a function of star formation rate if, for example, the escape fraction of Lyα photons varies. Another possibility is that a large fraction of the Lyα emission comes from previously undetected sources or low surface brightness halos around galaxies. These factors will change the L-M relation and Figure 14. Probability distribution functions (PDFs) of Lyα luminosity density (solid curves) and halo-bias-weighted Lyα luminosity density (dashed curves) for our fiducial model. The PDFs are computed as proportional to the average Lyα luminosity or halo-bias-weighted Lyα luminosity in halos of mass M multiplied by the differential halo mass function dn/d log M . The blue curves use only the Lyα luminosity of central galaxies, and the red curves include contributions from satellite galaxies. All curves have been normalized to unity at their respective maxima. See the text for further details concerning the model. therefore the derived value of bL. Even with the above uncertainties, it is likely that bL is around 3.
With the above model, we can compute the contribution to the Lyα luminosity density from halos of different masses, which is simply proportional to the average Lyα luminosity in halos of mass M times the differential halo mass function dn/d log M . The solid blue curve in Figure 14 uses central galaxy Lyα luminosity only, which peaks around log(M/M⊙) = 12.25. Including the contribution from the satellite galaxies shifts the peak slightly to a higher mass, around log(M/M⊙) = 12.45, as shown by the solid red curve. The curve gives the probability density of a random Lyα photon to come from a halo of mass M . The Lyα emission detected through the cross-correlation technique probes the halo-bias-weighted luminosity density distribution. The dashed curves in Figure 14 show these probability distributions. With the satellite contribution included, the curve peaks around log(M/M⊙) = 12.6. Taking the full width at half maximum of the curve, the fiducial model implies that the signal in the cross-correlation should mainly come from Lyα emission in halos of mass (1-20)×10 12 M⊙.
Finally, as mentioned in Section 3.2, bL is not the same as bα once the radiative transfer effect is taken into account. A simple model shows that bα = bL + α1 with α1 a positive number (see Section 6.4). The value of α1 is not readily known without detailed radiative transfer modeling. Overall, we expect bα to be larger than ∼3. We choose to parameterize derived quantities in terms of (3/bα).
An additional uncertain factor to consider is the modification to clustering caused by redshift-space distortions. This is embodied in the f β parameter. We have seen in Section 3.4 that measurements of anisotropies in clustering give a measurement of βα = Ω 0.6 m (z = 2.55)/bα = −0.76 ± 0.36. This negative value of βα is of the form expected to be caused by radiative transfer effects on clustering (Zheng et al. 2011a) and is opposite in sign to the usual Kaiser (1987) peculiar velocity redshift-space distortions. Neverthless, this redshift-space distortion model was shown in Section 3.4 to give a reasonable fit to the data and one can use this to compute the factor f β as f β = 1 + 1 3 (βq + βα) + 1 5 (βqβα) (from equation [10]). If we do this we find that the value of f β = 0.80±0.15, which we take as a reasonable estimate of the reduction of the monopole term due the spreading the correlation along the line of sight. Even though gravitational evolution is not the physically correct model for interpreting our observations owing to the negative value we obtain for βα, a model with redshift-space distortion plus radiative transfer effect does seem to work reasonably well here (see Section 6.4).
(20) We convert this into a comoving Lyα luminosity density ǫα using where c is the speed of light and λα = 1216Å. We find the value ǫα = 3.1 × 10 41 (3/bα) erg s −1 Mpc −3 . We then use Equation 17 to convert this into a measurement of the star formation rate density ρSFR(z = 2.55) = (0.28 ± 0.07)(3/bα) M⊙yr −1 Mpc −3 . As mentioned before, the conversion depends on the assumption about the underlying stellar population. A younger population and lower metallicity would lead to a lower inferred SFR, which could be an important effect for interpreting our results. Keeping this possibility in mind, we proceed with the discussion by using the above result from the commonly adopted conversion factor. We plot this result in Figure 13 as the red point, for the chosen value of bα = 3. We note that the true value for the parameter bα is not well determined, and that the ρSFR datapoint should be scaled up and down by the factor 3/bα. Our discussion above suggests that bα is likely to be larger than 3, so that it is more probable for ρSFR to be scaled downwards than upwards. Figure 13 also shows various estimates of ρSFR from Lyα emitter surveys as well as UV continuum estimates of the star formation rate. We can see that our measurement is about 30 times higher (for bα = 3) than the Lyα emitter based measurements of Gronwall et al. (2007), Ouchi et al. (2008) or Cassata et al. (2011). Note that these Lyα emitter based measurements result from a direct integration over the observed Lyα luminosity function without corrections for any possible dust effect.
A complication which adds substantial uncertainty is dust absorption, which may affect Lyα and continuum ra-diation differently. One estimator of the level of dust extinction in the continuum is the rest-frame UV continuum slope, β, which specifies how the flux density of a galaxy varies with wavelength (i.e., f λ ∝ λ β ) in the UV continuum region (∼ 1300Å to ∼ 3500Å) of its spectrum. If an intrinsic β-dust extinction relationship is assumed, (usually that measured empirically from z ∼ 0 galaxies by Meurer et al. 1999), one can use observations of β for high redshift galaxies to infer a dust-corrected UV luminosity and then star formation rate. This has been done by several authors, including Bouwens et al. (2009Bouwens et al. ( ,2010. In Figure 13 we show as a blue band the compilation of dust-corrected UV star formation densities from Bouwens et al. (2010) computed using this technique. According to Bouwens et al. (2010), the dust correction for a limiting luminosity of 0.3 L * z=3 is 6.0±2.5 at redshift z = 2.5. Support for the validity of these corrections comes from the agreement of dust-corrected UV values with ρSFR estimated from infrared observations (see the recent review by Madau & Dickinson 2014).
As for the effect of dust on the Lyα radiation, Cassata et al. (2011) have speculated that it could be even stronger than proposed by Bouwens et al. (2010), as the Lyα emitter inferred ρSFR is less than 20% of the non-dust corrected UV continuum value (as can be seen in Figure 13). One reason which favors this interpretation is the fact that resonantly scattered radiation does have to cover a longer path length than continuum radiation before it leaves a galaxy so that it could be more vulnerable to dust extinction. On other hand, there is one well known mechanism (Neufeld 1991) which could preferentially enhance the escape of Lyα radiation: in a clumpy medium of dusty clouds, continuum (UV) photons are absorbed as soon as their path crosses an optically thick dust cloud, whereas Lyα photons can bounce off the cloud surfaces and find their way through the clouds to escape, leading to a lower extinction for Lyα than for continuum photons if the intercloud medium is sufficiently devoid of dust. The anisotropic escape of Lyα radiation (Zheng & Wallace 2014) caused by, for example, a bipolar galactic wind, can also help make Lyα photons follow the path of lower extinction optical depth to escape, while UV continuum photons are emitted isotropically and on average experience more extinction. From Figure 13, it appears that some mechanism of this sort is needed if we are to explain our results. We discuss these issues further in Section 6.
In conclusion, the rather surprising result seen in Figure 13 is that the fiducial value of the Lyα surface brightness from our measurement is consistent with all Lyα photons produced in stars at z = 2.55 escaping from their host galaxies and being detected. The dustcorrected results of Bouwens et al. (2010) imply ρSFR = 0.19 +0.08 −0.06 M⊙yr −1 Mpc −3 at z = 2.55, and from our measurement, ρSFR = (0.28 ± 0.07)(3/bα) M⊙yr −1 Mpc −3 . This means that bα > 3.0 is needed for our measurement not to imply detection of more Lyα photons than are actually produced at more than the 1 − σ level. Our simple model for bα does satisfy this limit. We note that the intensity mapping technique we use in this paper will detect Lyα photons which are scattered into our sightline from arbitrarily large distances from the emitting galaxy, and at arbitrarily low surface brightness. One can characterise our measurement as being consistent with a "total escape fraction" of Lyα photons from star-forming galaxies of 100 %. This to-tal escape fraction could be much higher from the "detected escape fraction" measured from traditional surveys of Lyα emitters, which have surface brightness limitations. We also discuss this in more detail in Section 6.

DISCUSSION
We can frame further discussion of our results in terms of the following questions: (i) Can the observed Lyα surface brightness be explained by known Lyα emitters? We have seen in Section 5 that the simple answer appears to be no.
(ii) Can the observed Lyα surface brightness be explained by faint Lyα emitters below the threshold of published surveys? Ouchi et al. (2008) have shown that changing the extrapolated luminosity function faint end slope α from −1.0 to −2.0 changes the total integrated Lyα luminosity density they infer by only 20%. The answer to this question is therefore also no.
(iii) Can the observed Lyα surface brightness be explained by extended halos around the known Lyα emitters? We show in Section 6.3.5 below that this is not the case, and that extended halos around known Lyα emitters, while adding to the mean surface brightness, fall short of accounting for our results by an order of magnitude .
(iv) Can the observed Lyα surface brightness be explained by star-forming galaxies, when these are estimated based on the extinction-corrected UV continuum SFR density? We have seen in Section 5 that there appears to be just enough star formation per unit volume in the Universe at z = 2.5 that if most of it led to escaped Lyα emission, this could explain what we are seeing. We further discuss the implications of this below in Sections 6.2, 6.5 and 7.2.
(v) What is the contribution to our meaasurement of sources of Lyα emission beyond star-forming galaxies? We address this in Section 6.3 below, showing that other contributions are likely to be very small.

Potential systematic errors in the measurement
Our measurement of Lyα intensity clustering relies on statistical cross-correlation techniques applied to a large sample of spectra with relatively low signal to noise, all of which were targeted at bright foreground galaxies that we have removed in post-processing. There is therefore ample scope for small instrumental or other effects to influence the signal we measure. Bearing this in mind, we have carried out a range of tests, detailed in Appendices A-C, to make as certain as possible that the signal is real. Most importantly, these include tests of our methods for eliminating contaminating light from neighboring fibers, which does have a strong effect. We have also tested the effect of eliminating quasars clustered with the quasars we are using in our cross-correlation, and we have checked the dependence of the signal on the luminosity of the galaxy fiber target and on the quasar luminosity. For these latter tests, we have found no significant effect and therefore conclude that the signal is real, to the extent we have been able to ascertain.
There remains the possibility that some previously unknown effect, instrumental or otherwise, is responsible for the cross-correlation signal. This clustering signal would have to have a shape and amplitude consistent with ΛCDM, yield a star formation rate density consistent with dustcorrected UV estimates, have clustering distortions of the form predicted by Lyα radiative transfer effects, and pass all the tests mentioned above. This would constitute an unlucky coincidence; in particular, contamination by quasar light seems very unlikely once we have eliminated the effect of neighboring fibers in the BOSS camera and we have tested the absence of a dependence on the quasar measured flux. Galactic dust absorption or other systematics on the spectral continua also cannot induce the enhanced emission that is measured at the inferred Lyα wavelength at the redshift of each observed quasar. However, the reader should bear in mind this possibility until future independent work is able to confirm our measurement.
We note that the effect of gravitational lensing on our results should be zero, even though we use spectra of bright galaxies, because our Lyα measurements are of surface brightness, which is conserved under lensing. If on the other hand we were detecting Lyα emitters in the fibers and obtaining their luminosity function, this would be subject to the well known magnification bias (e.g., Turner 1980). In our case, we are computing the cross-correlation function of the surface brightness measured from all fibers with quasars, and the magnification of Lyα sources cannot change the amplitude of the cross-correlation. Dust associated with the foreground galaxies might reduce the Lyα emission coming from higher redshift, but this could only further increase the inferred brightness of the Lyα background.

Star-forming galaxies and the photon budget
We have found in Section 5 that our detected signal of crosscorrelation of Lyα surface brightness with quasars implies a brightness for the mean Lyα photon background given by equation (20). This at the same time implies an emissivity of Lyα radiation of ǫα = 3.1 × 10 41 (3/bα) erg s −1 Mpc −3 . This emissivity can be reexpressed in terms of the rate at which Lyα photons must have been emitted for each baryon in the universe at the mean redshift of our observation, z = 2.55. Using the comoving number density of baryons n b = 2.5 × 10 −7 cm −3 , and an expansion rate at z = 2.55 of H(z = 2.55) = 261 km s −1 Mpc −1 (using the parameters Ω b h 2 = 0.0221, H0 = 68 km s −1 Mpc −1 , and Ωm = 0.315, consistent with the most recent determinations from Planck in Ade et al. 2014), we find the following result: The first, most simple assumption we make is that these photons are mostly originating from star formation in galaxies. The Lyα photons created by recombinations in the HII regions produced around massive stars can then be scattered out to gaseous halos surrounding galaxies, from which they give rise to the background we detect in the quasar-Lyα emission cross-correlation. As discussed in Section 5, this implies a very large star formation rate at z = 2.55. Equation (17) can also be recast in terms of the number of Lyα photons emitted for each baryon that forms stars, nα/n b : Comparing to equation (22), we see that this implies that, for bα = 3, about 10% of all the baryons in the universe would have to turn into stars if the star formation rate is maintained over the age of the universe at z = 2.55, (2/3)H −1 (z = 2.55) = 2.5 × 10 9 years. The difficulty with this very high star formation rate is twofold: estimates of the total fraction of baryons in the form of stars at present are near 6% (Fukugita & Peebles 2014; Shapley 2011), and as described in the previous section, the total star formation rate at the redshift of our measurement can reach this value only for the maximum estimates of dust absorption, which would imply that while the UV continuum has to be absorbed by factors of ∼ 5, the Lyα photons would have to emerge suffering little dust absorption. A first possible solution to the problem of this extremely high inferred star formation rate is to modify the Initial Mass Function (IMF) of the stellar population that is assumed in deriving the Lyα photons emitted per baryon in equation (22) from population synthesis models. If the slope of the IMF is flatter in the high mass range of 10 − 100 M⊙, then stars above ∼ 20 M⊙, which dominate the production of ionizing photons, increase their abundance compared to ∼ 10 M⊙ stars, which dominate the observed UV continuum. Most of the ionizing photons can be absorbed in HII regions, but most of the Lyα photons produced by subsequent recombinations can escape, and their production rate can be increased compared to that inferred from the UV continuum observations. If the IMF stays flat down to lower masses, that can also greatly reduce the total star formation rate that is implied, as well as the stellar mass that is derived for the present universe which is measured from the old stellar population dominating the present luminosity of galaxies. If this flat IMF occurs particularly in massive galaxies with high metallicity, then the UV continuum observed at z = 2.55 can be further reduced due to the suppression of blue horizontal branch stars, and the luminosity-weighted bias factor of the Lyα emission can be further increased. A top-heavy IMF during the epoch when most stars were formed implies a large increase in the production of heavy elements, but this may be consistent with observations that show relatively high metallicities in the intracluster medium and in massive galaxies (Renzini & Andreon 2014).
For the rest of this section we assume a standard IMF with the ratio of Lyα photons produced per baryon forming stars in Equation (22). We next enumerate additional sources that might contribute to this background other than star-forming galaxies, estimating if any of these contributions might be substantial. Finally, we discuss the effects of Lyα radiative transfer on the quasar-Lyα background crosscorrelation we detect.

Other sources of Lyα emission beyond star-forming galaxies
We shall discuss here four possible contributions to the quasar -Lyα emission cross-correlation not arising from star-forming galaxies clustered around the quasars: (1) Scattering of the quasar Lyα broad emission line by the Lyα forest.
(2) Fluorescence of the ionizing radiation from the quasar.
(3) Scattering of the cosmic UV background by the Lyα forest in the overdense IGM around quasars. (4) Fluorescence of the general cosmic ionizing background by the overdense IGM around quasars. (5) Lyα cooling radiation from radiative dissipation of gas in halos that are correlated with quasars.

Scattering of the quasar Lyα broad emission line
The average observed flux of our sample of BOSS quasars within the central ∼ 2000 km s −1 of the Lyα broad emission line is close to fα ∼ 10 −16 erg s −1 cm −2Å −1 . In general, the mean fraction of light that is found to be absorbed by the Lyα forest at z = 2.55 is 1−F ≃ 0.2 (F is the mean transmitted fraction; e.g., Faucher-Giguère et al. 2008). At a characteristic impact parameter of ∼ 10h −1 Mpc (inside which our cross-correlation signal is strongest), corresponding to an angular separation θ ∼ 500 arcsec, the surface brightness of the scattered radiation should be fα(1−F )/(ψπθ 2 ), where ψ is a dimensionless number that depends on the geometry of the scattering gas around the quasar and the shape of the Lyα emission line, and has a value ψ ≃ 4 for a uniform gas density and a flat emission line profile. This yields a surface brightness ∼ 7 × 10 −24 erg s −1 cm −2Å −1 arcsec −2 , more than two orders of magnitude lower than our measured excess surface brightness at an impact parameter of 10h −1 Mpc from a quasar, as shown in Figure 2. Scattered light from the quasar broad Lyα emission line is therefore a negligible contribution to our detected Lyα background. This is consistent with our test in Appendix C.3 showing no dependence of the cross-correlation amplitude on the quasar luminosity.

Fluorescence of the quasar ionizing radiation
An excess of Lyα emission around the quasar may also arise from fluorescence of the ionizing radiation from the quasar. The ionizing radiation is absorbed in intergalactic absorption systems, and about half of the energy is reemitted as Lyα photons when the recombinations that maintain ionization equilibrium take place. These Lyα photons should be predominantly emitted from systems with column densities NHI ∼ 10 17 cm −2 , for which the Lyman limit optical depth is of order unity.
The ratio of this fluorescent emission to the scattered Lyα radiation from the quasar can be estimated as follows. The characteristic equivalent width of the Lyα emission line of quasars is cWα/λα ∼ 2 × 10 4 km s −1 (see, e.g., the composite spectrum of BOSS quasars in Figure 16). Assuming a continuum slope bluewards of the Lyα line fν ∝ ν −1.5 for the average quasar, we find that the ratio of the number of ionizing photons to Lyα photons emitted by a quasar is (3/4) 1.5 (2/3)λα/Wα ≃ 7. Around two thirds of these absorptions of ionizing photons result in a Lyα photon (the rest end up producing two-photon emission from the 2s state of the hydrogen atom), whereas only ∼ 20% of the Lyα photons on the blue half of the quasar emission line are scattered in the Lyα forest at z = 2.55. The ratio of the total number of fluorescent to scattered Lyα photons from a quasar is therefore 7(2/3)/0.2/0.5 ∼ 40.
To estimate the ratio of the contribution to the excess Lyα surface brightness near a quasar from the two processes, we also need to take into account the different pathlengths over which the photons are absorbed or scattered. Whereas Lyα photons in the blue half of the quasar emission line will be scattered over a pathlength corresponding to the typical half-width of the Lyα emission line, ℓα ∼ 5000 km s −1 , the ionizing photons will be absorbed over the longer mean free path between Lyman limit absorbers at z = 2.55, which is ℓLL ∼ 30000 km s −1 (Prochaska et al. 2014). The surface brightness produced near the quasar is proportional to the total number of Lyα photons produced divided by the pathlength over which they are emitted. Therefore, we conclude that the fluorescent emission from quasars should be 40/6 ∼ 7 times brighter in surface brightness than the scattered light from their Lyα broad emission line.
Even this factor of 7, however, raises the estimated surface brightness contribution from fluorescence at an impact parameter r = 10h −1 Mpc to ξqα(r) ≃ 5 × 10 −23 erg s −1 cm −2Å −1 arcsec −2 , which is still a factor of ∼ 30 below our measured value at this impact parameter from Figure 2.
We note here that this Lyα fluorescent emission from the quasar would not be of uniform intensity, but would predominantly arise from Lyman limit system absorbers in the quasar vicinity (see, e.g., the simulations in Kollmeier et al. 2010). This does not matter for our measurement, which is sensitive only to the mean surface brightness. In any case, our discussion suggests that the main difficulty for detecting this component of fluorescent emission caused by individual quasars probably lies in distinguishing it from the Lyα emission scattered in gas halos surrounding star-forming galaxies that are clustered with the quasar, which is what we believe dominates the production of the Lyα photons in our signal.
We also point out here that the quasar fluorescence contribution is affected by the anisotropy of the quasar emission and the time-delay between our observation of the quasar and that of the fluorescent light. Typically, if quasars are highly anisotropic and variable, the quasar emission in a direction away from the line of sight, or at ∼ 10 7 years ago corresponding to the time when the quasar was illuminating an absorber at an impact parameter r ∼ 10h −1 Mpc, would be systematically lower than the one we observe at present because of the flux-limited selection of quasar catalogues. These effects are likely to be important in any detailed modelling of the fluorescent emission due to the quasar.

Fluorescence of the cosmic ionizing background
Fluorescence from the mean cosmic ionizing background can also contribute to the quasar-Lyα cross-correlation we measure. In this case, if we assume the intensity of ionizing photons from distant sources is uniform around the quasar, the effect on the cross-correlation would arise only from the overdensity of Lyman limit absorbers near the quasar, because a uniform Lyα brightness does not contribute to our detected signal. We consider a value of the photoionization rate at z = 2.5 of Γ = 10 −12 s −1 , which corresponds to a proper intensity per unit wavelength of the ionizing background of i λ ≃ 3 × 10 −20 erg s −1 cm −2Å arcsec −2 at the ionization edge λ = 912Å (e.g., Faucher-Giguère et al. 2008).
This ionizing background is being absorbed over a mean free path λi ≃ cH −1 /10 at z = 2.5 (Prochaska et al. 2014). We assume that the ionizing background has a spectral slope −β in frequency (i.e., i λ ∝ λ β−2 below λ = 912Å). Each ionizing photon produces a Lyα photon with probability close to 2/3, so the intensity of the Lyα background that is produced is The factor (3/4) 2 arises from the ratio of the Lyα to the Lyman limit wavelength (one factor for the energy of each photon, and one for the fact that the intensity is per unit of wavelength). This approximates λi as being independent of frequency, which is valid for optically thick absorbing systems (for optically thin systems, β in the denominator needs to be replaced by β + 3, and the reality should be intermediate).
Using β = 2, we obtain iα ≃ 5×10 −20 erg s −1 cm −2Å −1 arcsec −2 , at z = 2.55. When observed at the present time, the intensity of this Lyα background is reduced to This mean intensity now needs to be multiplied by the crosscorrelation function of Lyman limit absorbers and quasars, to obtain the contribution to our measured cross-correlation. The value of the mass autocorrelation at impact parameter r = 10h −1 Mpc and z = 2.55 is ξm ≃ 0.05, which needs to be multiplied by the bias factor of quasars and the bias factor of Lyman limit absorbers. The bias factor for the BOSS quasars is known, bq ≃ 3.5, and that of the Lyman limit systems is unlikely to be larger than the value for DLAs, bD ≃ 2 (Font- Ribera et al. 2012), implying that f β ∼ 1.3 and bqbLLf β ξm(r = 10h −1 Mpc) is ∼ 0.45. This yields a value of the contribution to the quasar-Lyα cross-correlation of ∼ 4.5 × 10 −23 erg s −1 cm −2Å −1 arcsec −2 , a factor of ∼ 30 below our measured value from Figure 2. We note that this estimate is consistent with a cruder calculation which uses the ionizing background observationally inferred by Fontanot et al. (2014) from the comoving space density of quasars and star forming galaxies. After exploring the likely parameter space of limiting magnitudes and escape fractions, Fontanot et al. (2014) find that a central value for the ionizing background comoving emissivity is about 3 × 10 50 photons s −1 Mpc −3 at z ∼ 2.55. If the emissivity of the ionizing photons that are converted to Lyα photons is at a similar level, with a conversion efficiency of 2/3, the fluorescent Lyα emissivity is then 2 × 10 50 photons s −1 Mpc −3 . This is 3.2 × 10 39 erg s −1 Mpc −3 , ∼ 100 times lower than the value inferred from our results.
We therefore conclude that fluorescence from the cosmic ionizing background also cannot be responsible for the large Lyα brightness that we detect correlated with quasars.

Scattered Lyα radiation from the radiation background
Just as scattering from the quasar by the Lyα forest can in principle produce some contribution to this Lyα light around quasars, one may think that the general ultraviolet continuum background from distant galaxies and quasars could also give rise to an excess of Lyα photons near quasars because of the overdensity of the Lyα forest that scatters this background radiation. However, this effect actually cancels out for our observation. High redshift galaxies that are behind the quasar should show the Lyα forest absorption, reducing their ultraviolet flux, and this Lyα absorption will be enhanced because of the overdensity surrounding the quasar. The background sources (which are of course too faint to be detected individually, but always contribute to our total background) are therefore fainter near quasars compared to any random fields. This is exactly compensated by the scattered Lyα background from these same sources. For this reason, this contribution to the Lyα background could only be detected if the individual sources behind the quasar were individually detected and subtracted out, before evaluating the Lyα background intensity. In addition to Lyα photons scattered by the Lyα forest, the background would also have a contribution from background photons reaching the Lyγ wavelength when passing near the quasar and being downscattered to Lyα . This contribution is also very small, and is also nearly cancelled unless the population of ultraviolet sources creating the background evolves very fast over the redshift interval corresponding to the ratio of the Lyγ to the Lyα wavelength. (note that the Lyβ forest does not contribute to this diffuse emission by intergalactic gas because Lyβ absorptions can only end in a 2-photon emission from the 2s atomic state of hydrogen.)

Low surface brightness Lyα halos around galaxies
We have described in Section 5 how Lyα emission from galaxies could account for our measured result. Due to the magnitude of our signal, this Lyα emission does not appear to originate in a straightforward way from Lyα emitters seen in narrow band surveys at these redshift.
The Lyα luminosity from observed Lyα emitters is usually defined by pixels above certain surface brightness threshold or measured within an aperture (e.g., 2 ′′ ; Ouchi et al. 2008). In any case, the surface brightness threshold is well above the sky noise and higher than the level associated with our measurement in this paper. An extended Lyα emitting halo below sky noise level, resulting from Lyα photons scattered by neutral gas in the circumgalactic and intergalactic media and clustered unresolved Lyα sources, is predicted to exist around a star-forming galaxy based on radiative transfer modeling (Zheng et al. 2011b). Such diffuse Lyα emitting halos are detected from stacking analysis (Steidel et al. 2011;Matsuda et al. 2012;Momose et al. 2014).
We use the stacked Lyα surface brightness profiles and the fits to the diffuse Lyα halos in Momose et al. (2014) to estimate the contribution of these halos to Lyα emission in our detection. The luminosity inside the aperture of radius 2 ′′ roughly corresponds to the Lyα luminosity from the Lyα emitter survey. We find that at both z = 2.2 and z = 3.1, the diffuse Lyα emission outside of 2 ′′ is about one third of the luminosity inside the aperture. Therefore, the observed diffuse Lyα halos, regardless of their origin, may only contribute an additional one third to the Lyα emission from Lyα emitter surveys.
This implies that our hypothesized extended emission from scattering halos that may account for our measurement would have to arise in many more galaxies than those detected in Lyα -emission surveys.

Cooling radiation
Cooling radiation from gas in galactic halos can produce Lyα emission. A rough estimate of the cooling radiation can be made if we assume that cooling and star formation reach a steady state on average in a galaxy (i.e., 1 M⊙ of gas cools per year to feed 1 M⊙ yr −1 star formation).
For gas initially at virial temperature T dissipating the energy through cooling (as suggested by cosmological simulations; e.g., Fardal et al. 2001;Faucher-Giguère et al. 2010;Goerdt et al. 2010), we can estimate the corresponding cooling luminosity as follows: Let us assume that a fraction fc of all baryons in the Universe fall into halos of virial temperature T , and they dissipate all their energy by emitting Lyα photons. Initially the baryons fall into the halo and are shock heated to temperature T at a radius rv ≃ σt/6 in the halo, where t is the age of the Universe, and in the end they have to reach a radius rg by dissipating their energy through Lyα emission, where rg is the half-radius of the galaxy.
The circular velocity vc = √ 2σ is assumed to be flat, independent of radius. Then, the potential difference from rv to rg is φ = v 2 c log(rv/rc). The energy to be dissipated per baryon is therefore mpφ. The relation between v 2 c and the temperature T is σ 2 = kT /µ, where µ = 0.6mp for the fully ionized mixture of hydrogen and helium from Big Bang nucleosynthesis. So, the total energy dissipated per baryon in the universe is and the number of Lyα photons emitted per baryon is ǫα = ǫ β /(10.2eV) = 3.3fc kT 1.2 × 10 5 K log(rv/rc).
Using T = 3×10 6 K and rv/rc = 20 we get ǫα = 240fc. With plausible values of fc ∼ 0.2, this could amount to nearly 15% of the emission we observed, and more if the bias of this Lyα emission is high. Cooling radiation is therefore the alternative source of Lyα emission which comes closest to explaining our results, being plausibly less than an order of magnitude below our measurements. Cooling radiation could contribute in a substantial way to the observed Lyα emission.

Radiative transfer effect
The above estimates show that the Lyα emission relevant to our clustering measurements is likely to be dominated by contributions from star-forming galaxies. Regardless of the origin, as long as Lyα photons are scattered by neutral hydrogen, we expect them to be affected by a radiative transfer effect. The quasar-Lyα cross-correlation function we measure suggests that this effect is detected, shown as elongated contours along the line of sight on scales as large as ∼ 20h −1 Mpc. This feature cannot solely originate from the dispersion of the relative peculiar veclocity between quasars and galaxies (including redshift errors), but is consistent with the predicted main radiative effect on the clustering of Lyα emitters in Zheng et al. (2011a).
The radiative transfer effect predicted in Zheng et al. (2011a) is a result of the anisotropic emission of scattered Lyα photons from the anisotropic distribution of neutral gas density and velocity (hence anisotropic distribution of Lyα scattering optical depth). In particular, Lyα photons preferentially escape along the direction for which the neutral gas has the largest peculiar velocity gradient. Therefore, the Lyα emission from less overdense regions along the lineof-sight direction would appear enhanced, suppresing the line-of-sight fluctuation. This dependece on the line-of-sight peculiar velocity gradient is analogous to that in the Kaiser effect but of opposite sign, leading to the elongated correlation contours. For all the possible origins of Lyα emission discussed above, this radiative transfer effect should be at work as long as Lyα photons interact with neutral gas on large scales.
Our measurement can be used to constrain the parameters relevant to the radiative transfer effect. Following the simple model presented in Zheng et al. (2011a), the realspace overdensity δα (the Lyα surface brightness fluctuations here) can be related to the matter linear overdensty field δ and peculiar velocity field as where bL is the Lyα luminosity weighted bias of the underlying galaxy population and vz the line-of-sight peculiar velocity. The coefficient α1 represents a combined effect of the dependence of the Lyα radiative transfer (i.e., the Lyα escape emission) on the density and the transverse peculiar velocity gradient, and α2 specifies the impact of the lineof-sight peculiar velocity gradient (see Zheng et al. 2011a, in particular their Appendix A). Both coefficients are expected to be positive. In redshift space, we also need to add the Kaiser effect, and in terms of the Fourier component of the Lyα surface brightness fluctuations, we have The βα parameter (as constrained in § 3.4) is The factor "1" comes from the Kaiser effect. The bias factor bα in Equation 3 is bL +α1. With βα = −0.76±0.36, we have α2 = 1 + (2.35 ± 1.11)(bα/3) (N.B. this is the Cv coefficient defined in Wyithe & Dijkstra 2011). That is, the radiative transfer effect (indicated by a positive α2) shows up or has been detected at a level of ∼ 3.0σ for the fiducial value of bα = 3. Overall, the quasar-Lyα cross-correlation provides tentative evidence to a new clustering effect caused by Lyα radiative transfer. A better measurement with a larger data set and a more detailed modeling will help understand this effect and use it to constrain the neutral gas distribution.

Escape fraction and detected fraction
Assuming now that our fiducial value of the star formation rate density is correctly inferred from the quasar-Lyα cross-correlation, we conclude that it is consistent with the whole dust corrected star formation rate in Bouwens et al. (2010). At face value, this indicates that dust has little effect in reducing the Lyα emission, i.e., almost 100% of the Lyα photons produced from star formation escape. Clumpy dust clouds (Neufeld 1991) and anisotropic Lyα escape (Zheng & Wallace 2014) may be possible mechanisms for explaining the high escape fraction of Lyα photons, together with the much lower one for continuum UV photons.
We emphasize that the escape fraction is not the same as the fraction that is detected in Lyα emitter surveys. The latter, the detected fraction, comes only from the central, high surface brightness part of the Lyα emission from galaxies, and not from the extended Lyα halos of low surface brightness (Zheng et al. 2011b;Steidel et al. 2011;Matsuda et al. 2012;Momose et al. 2014). The detected fraction for Lyα emission at z = 2 to 3, inferred from comparing the Lyα luminosity density and Hα or H β luminosity density, is about 5% (e.g. Hayes et al. 2011;Ciardullo et al. 2014). This fraction is consistent with the ratio between the star formation rate density inferred from Lyα emitters and that from quasar-Lyα cross-correlation in this work (see the blue and red points in Figure 13).

Summary
We have carried out a cross-correlation analysis of residual light in SDSS/BOSS galaxy spectra and SDSS/BOSS quasars at redshifts from z = 2.0 − 3.5. We have used the Lyα emission line (which is at wavelengths λ = 3647 − 5471 A at these redshifts) to trace structure in the crosscorrelation. Our main findings are as follows: (1) We measure large-scale structure in the quasar-Lyα emission cross-correlation at a mean redshift z = 2.55 at the 8 σ level. The cross-correlation function shape is consistent with the linear ΛCDM correlation function.
(2) Looking at the clustering as a function of separation across and along the line of sight we see evidence at the 3.0σ level for distortions of clustering of the form expected to be caused by radiative transfer effects.
(3) We detect clustering independently in 4 subsamples at different redshifts, finding that the shape of the crosscorrelation function is consistent with the fiducial sample. The amplitude of the cross-correlation increases by a factor of 3 between z = 3.5 and z = 2.0, although this detection of evolution is marginal, being consistent with no evolution at the 2.0σ level.
(4) Although the spectra are too shallow to allow making a good catalogue of emission lines, we are able to weakly constrain the contribution of emission lines to our signal statistically by fitting lines, subtracting them and remeasuring the quasar-Lyα cross-correlation function. We find that lines with luminosities (measured in our 1 arcsec radius aperture) of LLyα > 8 × 10 41 erg s −1 may contribute 13 +11 −13 % of the quasar-Lyα cross-correlation amplitude at the relevant redshifts, but this contribution is still consistent with zero.
(5) In one of our sample tests, we measure the crosscorrelation to be independent of quasar luminosity. This is evidence that the large-scale clustering of Lyα surface brightness we measure arises mostly from Lyα emission associated with star formation, and not from any systematic error associated with the quasar light. (6) We estimate the plausible contribution to the quasar-Lyα surface brightness we measure from a variety of physical processes alternative to star formation galaxies that follow the same large-scale structure as quasars, such as fluorescence or scattering of the quasar emission or the ionizing background, and we find these contributions to be almost all negligible. Only cooling radiation may contribute in a substantial fashion (perhaps at the 15% level) to the Lyα surface brightness. (7) Using measurements of clustering from SDSS/BOSS quasars, we convert the measured amplitude of the crosscorrelation function bqbαf β µ to the product of mean Lyα sky brightness at z = 2.55 and linear bias factor of Lyα emission, finding µα (bα/4) = (3.9 ± 0.9) × 10 −21 erg s −1 cm −2Å−1 arcsec −2 . Assuming that this Lyα surface brightness is due to star formation (see points (5) and (6) above), we convert our measured value to the mean star formation rate density in the Universe at redshift z = 2.55, finding ρSFR = (0.28 ± 0.07)(3/bα) M⊙yr −1 Mpc −3 . This is consistent with dust-corrected UV continuum based estimates of star formation, but more than an order of magnitude higher than previous estimates of the SFRD from surveys of individual Lyα emitters.

Implications
We conclude that the high intensity of the Lyα background at z ≃ 2.5 that is derived from our measurement can only be reasonably produced by Lyα -emitting galaxies that are clustered with quasars. If our measurement is confirmed, the consequences for our understanding of galaxy formation and evolution are dramatic. The Lyα emission directly observed in Lyα emitting galaxies so far at these redshifts contributes only 0.01 M⊙ yr −1 Mpc −3 to the mean star formation rate density. Extended halos that have been seen around these Lyα -emitting galaxies detected in surveys can only contribute an additional ∼ 30% to this value. We have detected 21 − 35 (±1σ range) times more Lyα photons than the Lyα emitter surveys (with the uncertainty due to the factor (3/bα), and we have argued that most of these photons also arise from star formation. This amount of star formation represents most of the massive stars that are estimated to have formed in the universe and to have generated the present heavy elements.
To contribute to our measurement, these Lyα photons cannot be absorbed by dust before escaping the galaxies. The question that remains then is how this Lyα radiation could have been missed by previous observations. Putting forward and testing a detailed scenario is beyond the scope of this paper, but we speculate that all star forming galaxies at these redshifts, even if they do not have any detectable Lyα emission line in their central parts, are surrounded by low surface brightness halos that nevertheless have a high total integrated Lyα luminosity. These low surface brightness halos could allow the bulk of Lyα photons to be below the detectable levels in narrow band Lyα emitter surveys. As much of the star formation in the Universe at redshifts z = 2 − 4 occurs in massive galaxies, the implication is that a large fraction of the Lyα emission we detect is from giant low surface brightness halos around massive, bright, starforming galaxies, which absorb most of their continuum photons to reradiate the energy in the infrared, and yet let their Lyα photons escape.

The future
We note that the sky area covered by the million SDSS/BOSS (for our purposes randomly placed) fibers in the current study is ∼ 3 × 10 6 square arc seconds. This is approximately 1/200, 000 of the entire sky area, showing that Lyα intensity mapping holds an enormous promise as a probe of structure in the Universe. In addition to the quasar-Lyα emission cross-correlation employed in the present paper, one can imagine carrying out Lyα forest-Lyα emission cross-correlations and Lyα emission autocorrelation measurements. Correlations of Lyβ absorption and Lyα emission, and vice-versa would avoid common mode systematics in the fluxing of spectrographs and may reduce the possibility of related error. As the Lyα signal is distributed on large-scales, a way to efficiently carry out intensity mapping surveys (for example for baryonic acoustic oscillation experiments) may be to use integral field spectroscopy with relatively low angular (∼ 10 arcsec) resolution on a widefield (∼ few deg.) telescope. If bright point sources could be masked such an instrument could in principle capture a dataset containing 5 orders of magnitude more information on large-scale clustering in the Universe than our present study.
of Tokyo, University of Utah, Vanderbilt University, University of Virginia University of Washington, and Yale University.
As mentioned in Section 3 a potentially very important systematic is the contamination of the galaxy spectra by quasar light. This is particularly relevant for our cross-correlation measurement because we are searching for light from sources that are clustered with quasars, and light from the quasar itself could mimic this if not treated carefully. The issue arises because the 1000 fibers that orginate from the plate in the telescope focal plane are fed into (each red and blue) spectrograph and dispersed onto a single CCD. 1000 spectra are therefore extracted from a 4k×4K chip and some light from neighbouring fibers can end up in CCD columns that are centered on other fibers.

A.1: Measured Quasar-Lyman-α surface brightness correlation
We can test for this contaminating light by use of two facts about the observational setup. The fibers are arranged by fiber number (from 1-1000) on the CCD, so that for each spectrum we can measure light from a certain number of fibers away. Quasars and galaxies that we use in our crosscorrelation will also often be measured on different plates, so that there is no possibility for this contamination. We can therefore also test our result by restricting the calculation of the quasar-Lyα cross-correlation to pairs of quasars and galaxies lying on different plates.
In the first test, we have computed the quasar-Lyα cross-correlation of Equation 1 restricting ourselves to only quasar-galaxy spectrum pairs separated by specific values of ∆ fiber , the difference in fiber number of the two spectra. Spectra with ∆ fiber = 1, are adjacent on the CCD, for example, and have the highest potential for cross-contamination of light. We reduce the quasar-Lyα cross-correlation for each value of ∆ fiber by averaging the ξqα(r) results over a range from r = 5 − 50h −1 Mpc. Our conclusions about stray light contamination are not dependent on the range picked.
The results are shown in Figure 15, where we show ξqα (r = 5 − 50h −1 Mpc) against ∆ fiber . In the plot the value of ξqα (r = 5 − 50h −1 Mpc) for our fiducial sample (Section 3.1), which uses all pairs with ∆ fiber ≥ 6 is shown as a horizontal line. We can see that for ∆ fiber = 1, the mean surface brightness inferred from the cross-correlation is over 20 times the fiducial value. It remains significantly higher for all ∆ fiber ≤ 4. This is an indication that even when quasar and galaxy spectra are separated by 3 other spectra that light from the quasar is able to leak into the galaxy spectrum and contaminate it. Of course the particular region of the spectrum we are looking at (close to the Lyα emission line at the redshift of the quasar) is the one in which the quasar is extremely bright, and one would not expect it to contaminate other parts of galaxy spectrum as much. Nevertheless, for our project, this is the important region of the spectrum, and we therefore must apply a cut on ∆ fiber . Figure 15. A test of stray light from quasars in nearby fibers contaminating the galaxy spectra. We show the quasar-Lyα emission surface brightness cross-correlation averaged over quasarpixel pair separations of r = 5−50h −1 Mpc plotted against ∆ fiber , the difference in fiber number between the quasar and the galaxy spectra. A value of ∆ fiber = 1 means that the quasar and galaxy were in adjacent fibers and their spectra where dispersed next to each other in the CCD. The error bars are jackknife errors computed using 100 subsamples of the data in each case. The horizontal solid line is the quasar-Lyα emission surface brightness cross-correlation averaged over quasar-pixel pair separations of r = 5 − 50h −1 Mpc for our fiducial computation (see Section 3.1), which uses information from all quasar-pixel pairs 6 fibers apart or greater. The dotted line (which lies close to the solid line) is the equivalent result but only using quasar-pixel pairs which are on different plates. Figure 15, we have chosen that ∆ fiber ≥ 6 in our analysis.

Based on
A related test is to compute ξqα (r = 5 − 50h −1 Mpc) for quasars and galaxies that are on different plates (and therefore not dispersed at the same time onto the CCD). We show the results as a dashed line in Figure 15, which is indistinguishable from the results from our fiducial analysis. This is good evidence that our cut ∆ fiber ≥ 6 is sufficient to eliminate stray quasar light.
Yet another related test is to eliminate all galaxy fibers which have a quasar within a certain number of fibers. This is different from eliminating quasar-pixel pairs on a case-bycase basis because it will also eliminate any potential contamination which could come from quasars being clustered in space with other quasars. In the fiducial case, the contaminating quasars will have been eliminated directly, but those which are clustered with them could still contaminate the neigbouring galaxy spectra. One good reason to believe that this is not occuring to a detectable degree is that any contamination of this type would be a convolution of the redshift space autocorrelation function of quasars and the contaminating surface brightness from the quasar lya line, which is highly asymmetric and extremely elongated along the line of sight (see below). This does not appear to describe the measured signal (e.g., Figure 4).
Nevertheless, we have carried out the test, which elim- Figure 16. The mean surface brightness from all DR10 quasar spectra as a function of wavelength. We shift all spectra to their rest wavelength and stack them with equal weight to make the curve.
inates 50% of the galaxy fibers from the data sample. We find an amplitude bqbαf β µ = 1.98 +0.66 −0.65 × 10 −20 erg s −1 cm −2Å−1 arcsec −2 . and shape Ωm = 0.69 +0.71 −0.36 . These are consistent at about the 2 σ with the measurement from the fiducial sample. We note that the detection level of the signal is only about 3 σ overall, compared to ∼ 9σ for the fiducial measurement. This can be explained by the fact that eliminating all galaxies which have a quasar within 5 fibers will disproportionately affect the number of close quasar-pixel pairs. By directly counting, we find the number of quasarpixel pairs with separations below 40 h −1 Mpc has fallen by a factor of 3.5 rather than the factor of 2 expected if pairs were drawn uniformly from all fiber separations.
We have also carried out another, similar test, which eliminates a smaller fraction of the dataset, but which should have the same effect. For this test we remove from our list of galaxy pixels all pixels which have more than 1 quasar within r = 50h −1 Mpc. In this way, we eliminate all potential cross-contamination from quasars clustered with the target quasars in the cross-correlation, but at a much reduced cost (doing this only eliminates 0.3 % of the galaxy spectrum pixels). After fitting the cross-correlation of this sample, we find an amplitude bqbαf β µ = 3.04 +0.37 −0.45 × 10 −20 erg s −1 cm −2Å−1 arcsec −2 . and shape Ωm = 0.33 +0.13 −0.10 . This is consistent with the fiducial result within the error bars. It is slightly lower (by less than 1 σ), which should not be surprising, as we are presumably preferentially removing some pixels in overdense regions. Overall these two test results are a good sign that significant cross-contamination from quasars clustered with the target quasar is not ocurring.

A.2 Modelling the stray light contamination
If the excess surface brightness seen above in galaxy spectra which are close on the CCD to quasar spectra is indeed due to cross-talk from quasar light, we would expect the contamination to have a quasar-like spectrum. To examine The stack is centered on the redshift of the Lyα line at redshift z = 2.55 and the fiducial cosmology has been used to convertÅ to comoving h −1 Mpc. The stack is plotted as a function of r and is independent of r ⊥ . Bottom panel: The measured quasar-Lyα emission correlation for quasar-galaxy pixel pairs separated by 1 fiber (∆ fiber = 1) as a function of r and r ⊥ . In both panels a Gaussian filter with σ = 4h −1 Mpc was used to smooth the image (as in Figure 4). this, we first make a stacked spectrum of the DR10 quasar sample, by averaging all the spectra together in the quasar rest frame with unit weight. This is shown in Figure 16. We can see from that figure that the quasar Lyα line is extremely broad, with a FWHM of ∼ 50Å in the rest frame. There is also the noticeable emission feature due to NV on the red wing of the line.
We next take this stacked spectrum, and convert the wavelength units into comoving h −1 Mpc at z = 2.55, the mean redshift of our measurements. If the quasar light is truly a contaminant, we expect the strength of the contam-ination to depend on the difference between fiber numbers (∆ fiber , as defined above) and not on the actual physical separation between quasar and pixel across the line of sight (r ⊥ ). Plotting the quasar-Lyα cross-correlation for particular ∆ fiber values as a function of r ⊥ is therefore a useful check that light contamination is occuring as we believe. The prediction for the contamination, taken directly from the quasar spectrum is shown in the top panel of Figure 17. In that Figure we have scaled the overall amplitude by a factor of 2 × 10 −3 (see below). We can see that the large width of the Lyα line, and its asymmetry due to the NV feature are both apparent.
The data with the strongest contamination should be from galaxy pixel-quasar pairs which are 1 fiber apart (∆ fiber =1). We therefore plot the quasar-Lyα correlation for just those pixel-quasar pairs in the bottom panel of Figure  17. It is immediately apparent that there is a stripe across the middle of the plot corresponding to the contaminating quasar Lyα emission line. The contamination is clearly asymmetric, extending further to positive values of r than negative, and overall is visually consistent with the prediction based on the stacked quasar spectrum (top panel of Figure 17). Because the quasar-pixel pairs that are exactly 1 fiber apart only comprise a very small fraction of the data, we expect there to be lots of noise in the plot, particularly for large values of r ⊥ . This latter is because close pairs of quasars and galaxies on the sky are more likely to be close in ∆ fiber . There is neverthless enough range that we can examine by eye whether the contamination depends on r ⊥ , with the answer appearing to be no. We examine this more quantitively below.
When comparing the top and bottom panels of Figure  17 we note that the amplitude of the stacked quasar spectrum is scaled by a factor of 2 × 10 −3 , which is the level required to match the observations in a χ 2 fit (see below). This means that 0.2 % of the quasar light is scattered into neighboring spectra (i.e. spectra with ∆ fiber = 1).
If there is a signal that is actually coming from quasar-Lyα emission correlations in the Universe, one would expect that to be superimposed on top of any contamination from stray light. We have therefore carried out a simple joint fit of ξqα with a sum of the fiducial theoretical model from Figure 4 (right panel) scaled by an amplitude factor aCDM and the stacked quasar contamination scaled by a factor acontam. The fit is therefore as follows: where ξ stack (r ⊥ , r ) (multiplied by 2 × 10 −3 ) is shown in Figure 16 (top panel). We separate the observational ξqα results by ∆ fiber and then carrying out a χ 2 fit of the form given by Equation 31 to determine the best amplitude parameters acontam and aCDM as well as the confidence limits on those parameters. We carry out the fit for quasar-pixel pair separations r ≤ 50h −1 Mpc, although our conclusions are insensitive to this value.
The results are plotted in Figure 18. For the ∆ fiber = 1 dataset, we find a total χ 2 for the best fit of 1526 for 988 datapoints, which shows that the light contamination which dominates the fit is fairly well modelled by the quasar stack. It is perhaps not surprising the fit to the contamination is not perfect, as the contamination model is extremely simple. We have plotted the 1,2 and 3 σ confidence contours on the parameters aCDM and acontam in Figure 18. The central value of acontam is 1.75 × 10 −3 . The 1 σ confidence interval on aCDM is consistent with unity, and also with zero. The results from ∆ fiber = 1 are therefore consistent with the sum of significant contamination from quasar light and a quasar-Lyα correlation at the level predicted by the CDM model fit to the fiducial dataset.
We have carried out the same fitting to datasets with different restricted values of ∆ fiber in the other panels of Figure 18. We can see that for cases with ∆ fiber ≤ 4 there is a significant detection of a contamination contribution. The value of acontam is lower for ∆ fiber 2 − 4 than for ∆ fiber = 1 at the factor of 10 − 4 level. In all cases the measurements are also consistent with an additional quasar-Lyα correlation at the level of aCDM = 1. For ∆ fiber = 5 − 10 and ∆ fiber ≥ 6, the contamination is not detectable.
The results of the fitting shown in Figure 18 are a good sign that the signal and contamination are behaving in a particular way which can be accounted for by excising pairs of pixels from close fibers. It is useful however to examine the results for individual datapoints in more detail. Of partic-ular interest is the measured quasar-Lyα surface brightness as a function of r ⊥ . As mentioned above, we expect contamination at the CCD level for a given fiber separation to be independent of r ⊥ . We have taken the measured ξqα(r ⊥ , r ) values for different datasets limited by ∆ fiber and averaged them between r = ±50h −1 Mpc .We plot this average, ξqα as a function of r ⊥ in Figure 19. In each panel we also plot a dashed horizonal line showing a level of contamination that is independent of r ⊥ , a dotted line showing the level of the signal from quasar-Lyα clustering (we plot ξqα (r ⊥ ) computed from the model shown in Figure4, right panel), and a solid line which is the sum of the two.
We can see that in the ∆ fiber = 1 panel there is the strong contamination signal which is consistent with being independent of r ⊥ . The small relative contribution of the actual signal makes little difference in this panel. In the ∆ fiber = 2 and ∆ fiber = 3 panels the level of contamination is lower (note that the y-axes are different in the different panels), but it still dominates over the expected signal, with again no sign of a dependence on r ⊥ . For the fiducial panel (bottom right), we can see that the measured ξqα (r ⊥ does have significant dependence on r ⊥ , and this has a similar form and amplitude to the best fit model to the data in Section 3.2. There is no sign of a contaminating (independent of r ⊥ ) component to this measurement.
Having carried out these tests on contamination from the Lyα emission line from nearby quasars, it is pretty clear that the signal we are seeing in our fiducial dataset cannot be caused by straightforward scattered light coming from nearby fibers. There remains the possibility however that there is an additional contaminating component from some other mechanism that has a dependence on r ⊥ . This is difficult to imagine, but one can construct tests for this based on other spectral features than the Lyα line. In the next two subsections we do this, first centered on the quasar CIV emission line and then using stars and the Hα line.

A.3 Tests with CIV
As a test of scattered light contamination, we carry out a different correlation between galaxy spectrum pixels and quasar positions. The only difference from our standard analysis (Section 3) is that we compute the cross-correlation at the wavelength of CIV at the redshift of the quasar. We use a rest wavelength of 1550Å for CIV (close to the center of CIV doublet). If our interpretation of the contamination from nearby quasars from the previous subsection is correct, we would expect to see strong contamination from the CIV line for low ∆ fiber and then no signal ξqCIV for large ∆ fiber (because there should be no strong intergalactic CIV line emission).

Figure 20.
A null test using Carbon IV-quasar crosscorrelations. Here we show ξ qCIV measured from BOSS data averaged between r = ±50h −1 Mpc and plotted as function of r ⊥ . We show results for 2 different values of ∆ fiber as indicated in each panel. The points with error bars show the observational results. To aid comparison with previous plots, the prediction of the theoretical model from Equation 11 for the quasar-Lyα cross correlation is shown as a dotted line. A dashed line indicates scattered light contamination of the sort discussed in Section 3.2 which is independent of r ⊥ .
As with Figure 19, we plot ξqα (r ⊥ ), in Figure 20. We show results for ∆ fiber = 1 and for ∆ fiber ≥ 6. In the case of ∆ fiber = 1 there is strong contamination from the quasar, which does not depend on r ⊥ . The amplitude of the contamination is approximately three times smaller than that from Lyα , which is the ratio one would expect from the relative strengths of the two lines in Figure 16. For the ∆ fiber ≥ 6 we see that there is no evidence for contamination, and, although the error bars are relatively large, there is no sign of any ξqCIV signal. We have plotted the fiducial ξqα model on both panels with a dotted line. That the contamination behaves the same way as the Lyα results in the previous subsection but that there is no ξqCIV signal is further evidence for the reality of the ξqα measurement.

A.4 Tests with stars
A further test of stray light contamination which we can carry out is to cross-correlate the positions of stars with the Figure 21. A null test using Hα -star cross-correlations. Here we show ξ sH α measured from BOSS data averaged between r = ±50h −1 Mpc and plotted as function of r ⊥ . We show results for 2 different values of ∆ fiber as indicated in each panel. The points with error bars show the observational results. To aid comparison with previous plots, the prediction of the theoretical model from Equation 11 for the quasar-Lyα cross correlation is shown as a dotted line. A dashed line indicates scattered light contamination of the sort discussed in Section 3.2 which is independent of r ⊥ . galaxy pixels. We use all stars in the SDSS DR10 catalog (Ahn et al. 2014), which is 171612 objects in total. In this case we choose the Hα line, with rest wavelength 6562.8 A and cross-correlate star positions with galaxy pixels centered on this wavelength. We measure the resulting star-Hα cross-correlation, ξsHα (r ⊥ is for two values of ∆ fiber , 1 and ≥ 6. In this test, all the stars are obviously at the same rest wavelength, which effectively means that in this measurement we are stacking galaxy spectra together. This test is therefore likely to be more sensitive to the CIV test to residual artifacts in photocalibration of the spectrometer (see e.g. Figure 4 of Busca et al. 2013). Nevertheless it is useful to see whether any contamination has an r ⊥ dependence.
We plot the results for ξsHα (r ⊥ in Figure 21, where we can see that there is again a sign of stellar light contamination when ∆ fiber = 1 and that this contamination appears to be independent of r ⊥ . For ∆ fiber ≥ 6, the measurement is consistent with zero, again a sign that the cut in ∆ fiber removes scattered light contamination.
The results from the tests in Appendix A therefore suggest that although light contamination from quasars can be readily seen in the data, it can be eliminated to a high degree by excluding close pairs of fibers. When this is done, the resulting measured quasar-Lyα surface brightness correlation seems to be real. If instead it is produced by some as yet unknown systematic this effect would have to reproduce the r ⊥ and r dependence expected from a cosmological signal, be measurable only in the Lyα wavelength range and not around other emission lines, and be independent of the position of the spectra on the detector.
Having searched for and as far as we can tell eliminated the effects of strong light contamination, we now turn to other potential systematic effects in the next 2 appendices.

APPENDIX B: LARGE-SCALE SURFACE BRIGHTNESS CORRECTION
A potential systematic error on our measurement is any obscuration (such as that produced by Galactic dust) that may affect both quasar selection and Lyα surface brightness, thus producing a spurious cross-correlation between the two on the plane of the sky. One can correct for such a cross-correlation by searching for an ξqα (r ⊥ , r ) signal which is non-zero for large values of r, where no physical (3D) correlation would be expected. This correlation which will be a function of r ⊥ only (in the parallel line-of-sight approximation) can then be subtracted from the ξqα signal. We have done this, computing the surface brightness correction µα given by ξqα(r , r ⊥ )dr (32) where xa = 400h −1 Mpc and x b = 80h −1 Mpc and ξqα (r ⊥ , r ) is the qso-Lyα emission cross-correlation. The value of µα(r ⊥ ) is shown in the bottom panel of figure 22. We can see that there is no coherent structure, as µα(r ⊥ ) becomes both positive and negative as r ⊥ varies. We have also plotted the best fit linear CDM model (as a function of r) for the qso-Lyα cross-correlation (from Section 3.2), and we can see that this dominates over the µα(r ⊥ ) correction on small scales r ⊥ < 40h −1 Mpc, as we would hope.
There is also another analogous correction, but one that applies in the orthogonal direction. Any systematic effect which affects the Lyα emission in the line-of-sight direction (including redshift evolution in the Lyα surface brightness), and which is correlated with evolution in the quasar population with redshift (or at least evolution in the efficiency of quasar selection with redshift) could produce a spurious qso-Lyα emission cross-correlation. We compute how ξqα varies as a function of r for large r ⊥ values where there should be minimal physical clustering. In this case the the surface brightness correction µα is given by where we use the same values of xa and x b as in Equation 32. The results for µα(r ) are shown in the top panel of Figure 22. We can see that in this case there is a significant trend in the surface brightness correction with µα(r ) gradually decreasing from positive values of ∼ 4 × 10 −22 erg s −1 cm −2Å−1 arcsec −2 for r = −150h −1 Mpc to ∼ −10 −23 erg s −1 cm −2Å−1 arcsec −2 . for r = +150h −1 Mpc. As with the previous (r ⊥ ) correction, the CDM model dominates on scales r ⊥ < 40h −1 Mpc, where most of the clustering signal is located. In our analyses in the main body of the paper, we have applied these µα(r ) and µα(r ⊥ ) corrections to the computation of the qso-Lyα emission cross-correlation ξqα . We have done this on a quasar-pixel pair basis, computing r and r ⊥ for each pair, and then subtracting the appropriate value of µα. In all cases, the effect of the correction is small, as we would expect given the much large amplitude of the clustering signal compared to the surface brightness corrections that can be seen in Figure 22. For example, in our fiducial case, without applying the large-scale surfacebrightness corrections, we find a best fitting values of the shape and amplitude parameters of Section 3.2 of bqbαf β µ = 3.5 erg s −1 cm −2Å−1 arcsec −2 . and Ωm = 0.26, which represent differences of 0.5σ and 0.5σ in the parameters respectively.

APPENDIX C: SAMPLE TESTS
In this appendix we report on three consistency tests of our cross-correlation results. One test uses sky fibers to search for Lyα emission instead of galaxy fibers, and the second looks at how the luminosity of the orginally targeted  galaxy affects the quasar-Lyα cross-correlation. The third test investigates how quasar luminosity affects the crosscorrelation.

C.1 Sky fibers
Ideally one would like to be able to use fibers positioned on random areas of the sky to carry out intensity mapping, without needing to worry about subtracting a foreground galaxy. This approach is being carried out for example by HETDEX (Blanc et al. 2011). In our case, however, there are a number of such random fibers which were obtained for each plate, to use in sky-subtraction. The number of fibers available for use to use is 146065, approximately 15% as many as there are galaxy fibers in our fiducial dataset. We have taken these sky fibers and carried out the quasar-Lyα cross-correlation of Equation 1. The calculation was the same as in the fiducial case, including subtraction of a background from fiber ( Figure 1).
We show the results for ξqα(r) for the sky fibers in Figure 24. We can see that the datapoints appear to be consistent with the trend delineated by the best fit to the fiducial computation (also shown, as a dash-dotted line), although the measurement is much noisier. We fit the same CDM linear correlation function as was carried out in Section 3.2 and show the best fit parameters and confidence contours in Figure 23. We find bqbαf β µ = 2.2 +1.0 −1.0 × 10 −20 erg s −1 cm −2Å−1 arcsec −2 and Ωm = 0.52 +0.48 −0.36 . Here we can see that there that the confidence contours are very wide, and that the clustering signal is present at only the 2σ confidence level. The results are also consistent with the results from the fiducial sample at the 1.5 σ level. At present therefore, this shows that the use of galaxy and sky fibers in computing the cross-correlation is approximately equivalent, although the uncertainty is obviously large.
In the future, it would arguably be best to use randomly distributed fibers to carry out the cross-correlation, to avoid any potential selection biases. In the current sample case these biases are below the level of detectability, but one can imagine two types of bias, related to the fact that galaxy fibers are selected in regions of above average sky brightness and sky fibers selected to be in regions of below average sky brightness. Both effects could in principle bias the measurement of Lyα surface brightness in opposite directions. We now carry out a further test to constrain this effect, by splitting the spectra into two halves bases on target LRG luminosity.

C.2 Sample split by galaxy luminosity
We divide the galaxy spectrum sample into two halves, based on the measured SDSS r band luminosity (no k-correction was applied). One half consists of galaxies above the median luminosity of the whole dataset (5.2 × 10 40 erg s −1 ), and the other half those below it. The median r band luminosities of the halves are 6.9 × 10 40 erg s −1 and 3.8 × 10 40 erg s −1 We note that the LRG galaxy surface brightnesses measured from BOSS spectra are of order 10 −17 erg s −1 cm −2 A −1 arcsec −2 , approximately 2 − 3 orders of magnitude brighter than the mean Lyα surface brightness µα that we have measured in Section 5. We therefore expect that it is unlikely that regions of excess background Lyα surface brightness could have caused certain LRGs to be preferentially selected and therefore bias our measurements.
We measure the quasar-Lyα cross correlation ξqα(r) for the two samples of spectra and show the results in Figure  25. Both subsamples show evidence of clustering that is consistent with the CDM model shape and the fiducial amplitude. This can be seen in a quantitative manner from Figure  26, where we show the fit parameters. We find bqbαf β µ =3.1 +0.6 −0.7 × 10 −20 erg s −1 cm −2Å−1 arcsec −2 , for the bright galaxy spectra and bqbαf β µ =(3.8 ± 0.5) × 10 −20 erg s −1 cm −2Å−1 arcsec −2 , for the faint galaxy spectra. Both samples are consistent with the amplitude of the fiducial result at the 1σ level (the fiducial sample fractional error bar is Figure 25. Quasar-Lyα cross-correlation ξqα(r) for two subsamples of spectra originally targeted at LRGs with (a) luminosity above the median value and (b) below the median value. The solid curve shows the best fit CDM model and the dash-dotted line the CDM fit to the fiducial sample (all LRGs).

bar +12%
−11% for fixed Ωm = 0.30). This indicates that any bias present in the mean surface brightness µα due to fiber target selection is likely to be at the ∼ 10% level or less.

C.3 Sample split by quasar luminosity
A potential source of Lyα emission clustered with quasars is recombination radiation from dense IGM systems illuminated by the quasars themselves (e.g., Kollmeier et al. 2010). We have made an estimate of the amplitude of this signal in Section 6 and find it likely that it is much smaller than the signal from star forming galaxies clustered with the quasar. One way of testing this directly is by splitting the sample into high and low luminosity quasars and then measuring the ξqα(r) signal for each. We note that this test can also function as a diagnostic for stray light contamination from quasars, as any stray light should be correlated with quasar luminosity.
In this section we do this, measuring ξqα(r) the subsample of quasars with SDSS r band luminosity above the Figure 26. CDM amplitude and shape parameters fitted to the Quasar-Lyα cross-correlation ξqα(r) for two subsamples of spectra originally targeted at LRGs with (a) luminosity above the median value and (b) below the median value. The points shows the best fit values of the amplitude (bqbαf β µ ) and shape (Ωm) and the contours show the 1, 2 and 3σ confidence intervals. The open circles show the best fit values of the fit parameters for the fiducial sample. median, and with the subsample with luminosities below the median. The median luminosity of the bright subsample is 3.45 times the median luminosity of the faint sample, and the median redshifts of the quasars in each are z = 2.66 and z = 2.44 respectively.
We can see in Figures 27 that the ξqα(r) results for two subsamples do not look very different. In particular if the quasar luminosities were causing significant fluorescent emission from nearby IGM material one might expect the brighter subsample to exhibit a steeper ξqα(r) on small scales, which is not the case. From Figure 28 we can see that the fit parameters for the bright subsample are within 1σ of the fiducial result and the faint subsample within 2σ. The amplitude of parameter (bqbαf β µ ) (for Planck Ωm = 0.30) is actually slightly lower for the bright subsample, being (2.3 ± 0.6) × 10 −20 erg s −1 cm −2Å−1 arcsec −2 , whereas for the faint subsample it is 3.9 +0.7 −0.6 × 10 −20 erg s −1 cm −2Å−1 arcsec −1 . Given that the quasar subsamples are different in luminosities by a factor of 3.45 one would expect there to be a significant difference the ξqα(r) results in the oppo- Figure 27. Quasar-Lyα cross-correlation ξqα(r) for two subsamples of quasar data with (a) luminosity above the median value and (b) below the median value. The solid curve shows the best fit CDM model and the dash-dotted line the CDM fit to the fiducial sample (all quasars).
site direction if quasar properties were significantly affecting the large-scale Lyα intensity around them. As mentioned above, this lack of dependence on quasar luminosity is also strong evidence that stray light from quasars scattered in the spectrograph is not contaminating our measured crosscorrelation signal.

C.4 Other: sample split by fiber position in spectrograph
It is possible to imagine other sample tests based on various observational parameters, which could help to identify sources of systematic error. One such test is to split the sample into two subsets based on position of fibers in the spectrograph. Of the 1000 fibers, those labelled with numbers 100-400 and 600-900 are positioned towards the center of the spectrograph and the rest are positioned more towards the edge. The reason for such a test is that the optical quality of the instrument cameras degrades towards the edge, Figure 28. CDM amplitude and shape parameters fitted to the Quasar-Lyα cross-correlation ξqα(r) for two subsamples of quasars with (a) luminosity above the median value and (b) below the median value. The points shows the best fit values of the amplitude (bqbαf β µ ) and shape (Ωm) and the contours show the 1, 2 and 3σ confidence intervals. The open circles show the best fit values of the fit parameters for the fiducial sample. and this might affect the signal to noise ratio and/or cause other problems.
We carry out this split, and find that the average residual flux after subtraction of the model LRG spectra is somewhat different for the two subsamples. The residual flux for the fiducial sample was plotted in Figure 1, where we saw that it fluctuated between approximately ±3 × 10 20 erg s −1 cm −2Å −1 arcsec −2 for most of the relevant redshift range. In Figure 29 we plot for each subsample (edge and center) ∆(SB residual ), the residual flux for the subsample minus the residual flux for the fiducial sample. We can see that there are indeed differences-the residual flux in fibers close to the edge of the cameras is systematically higher than the fiducial sample (and the center of camera sample) by ∼ 10 −20 erg s −1 cm −2Å −1 arcsec −2 over the wavelength range ∼ 4000 − 5000Å. As this covers the important Lyα redshift range z = 2.3−3.0, this could indeed have consequences for our measurements.
Using the appropriate mean residual flux for each subsample, we compute the quasar-Lyα emission crosscorrelation ξqα . It should be noted that the subsamples are Figure 29. The difference in average residual flux for two subsamples of LRG spectra selected to be close to the center of of the spectrograph camera (blue line) or edge of spectrograph camera (red line). The residual flux is that left after subtraction of the best fit galaxy model from each spectrum, and the difference plotted here is the difference between each subsample and the fiducial result plotted as a black line in Figure 1. not quite the same size (there are 400 "edge" fibers per plate vs 600 "center" fibers). We fit the usual CDM model parameters and find bqbαf β µ =3.9 +0.5 −0.5 × 10 −20 erg s −1 cm −2 A −1 arcsec −2 , and Ωm = 0.31 +0.14 −0.10 for the central fibers and bqbαf β µ =(1.5 ± 0.7) × 10 −20 erg s −1 cm −2Å−1 arcsec −2 , and Ωm = 0.44 +0.96 −0.33 for the edge fibers. The central fiber results are within 1 σ of the fiducial results, but the edge fiber parameters have larger errors and are 2.5 σ below the fiducial results. We have also tried averaging the ξqα results from the two subsamples before fitting and find fit parameters consistent at the 1σ level with the fiducial results.
This test appears to show that the degradation of the optical quality of the spectrograph cameras close to the edge of their field of view does affect our ability to measure ξqα . Given this degradation, it is reasonable to expect that the error bars for the edge sample are underestimated and that this may account for the 2.5 σ difference with the fiducial sample. The fibers close to the center of camera sample appears to be responsible for most of our ability to determine ξqα . This information should be useful for future measurements of this type, but we do not use it to effect any changes to the fiducial analysis in this paper.