{'fraction_non_alphanumeric': 0.031675459496255956, 'fraction_numerical': 0.02216643975493533, 'mean_word_length': 4.365711676749229, 'pattern_counts': {'":': 0, '<': 7, '<?xml version=': 0, '>': 6, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Large-scale 3D mapping of the intergalactic medium using the Lyman Alpha Forest

Maps of the large-scale structure of the Universe at redshifts 2-4 can be made with the Lyman-alpha forest which are complementary to low redshift galaxy surveys. We apply the Wiener interpolation method of Caucci et al. to construct three-dimensional maps from sets of Lyman-alpha forest spectra taken from cosmological hydrodynamic simulations. We mimic some current and future quasar redshift surveys (BOSS, eBOSS and MS-DESI) by choosing similar sightline densities. We use these appropriate subsets of the Lyman-alpha absorption sightlines to reconstruct the full three dimensional Lyman-alpha flux field and perform comparisons between the true and the reconstructed fields. We study global statistical properties of the intergalactic medium (IGM) maps with auto-correlation and cross-correlation analysis, slice plots, local peaks and point by point scatter. We find that both the density field and the statistical proper- ties of the IGM are recovered well enough that the resulting IGM maps can be meaningfully considered to represent large-scale maps of the Universe in agreement with Caucci et al., on larger scales and for sparser sightlines than had been tested previously. Quantitatively, for sightline parameters comparable to current and near future surveys the correlation coefficient between true and reconstructed fields is r>0.9 on scales>30 h^-1 Mpc. The properties of the maps are relatively insensitive to the precise form of the covariance matrix used. The final BOSS quasar Lyman-alpha forest sample will allow maps to be made with a resolution of ~ 30 h^-1 Mpc over a volume of ~ 15 h^-3 Gpc^3 between redshifts 1.9 and 2.3.


INTRODUCTION
The structure of the Intergalactic Medium (IGM) can be studied using the Lyman-α forest, the absorption features due to neutral hydrogen seen in quasar spectra (Lynds 1971). Besides quasars, other background sources can be used to probe absorption along the line of sight. These include Lyman-break galaxies (LBGs) at redshifts higher than 2 (Steidel & Hamilton 1992;Steidel et al. 1998;Steidel, Pettini & Adelberger 2001;Adelberger et al. 2005;Lee et al. 2014a), although their much fainter magnitudes makes this significantly more difficult. Gamma Ray Bursts (GRB) can also be used (Totani et al. 2013). At z = 3, each individual spectrum provides one dimensional information for a ∼ 400 h −1 Mpc skewer in length if spectral coverage is obtained for the Lyman-α -Lyman-β region. This 1D information can be used to constrain the matter fluctuation amplitude (Weinberg et al. 1999; Kim et al. 2004) and the temperature-density relation ; Lee et al. 2015). Transverse information can be obtained with pairs of quasars (Dinshaw et al. 1994(Dinshaw et al. , 1995Fang et al. 1995; Petry, Impey & Foltz 1998; Hennawi & Prochaska 2007). Using this trans-verse information, the three dimensional correlation function of the Lyman-α forest on large scales (up to 100 h −1 Mpc) was measured for the first time, using data from the Baryon Oscillation Spectroscopic Survey (BOSS) Data Release 9 (DR9) (Slosar et al. 2011). BOSS has also enabled the first measurement of Baryon Acoustic Oscillations using the Lyα forest on even larger scales (Busca et al. 2012;Slosar et al. 2013;Delubac et al. 2014;Aubourg et al. 2014). Such BAO measurements at z ∼ 2 with the Lyman-α forest are currently the only way to put constraints on dark energy at these redshifts.
If the line of sight density is sufficiently high, inversion methods can be utilized in order to directly make three-dimensional maps of the intergalactic medium from collections of 1D spectra . Caucci et al. (2008) and Lee et al. (2014a,b) have shown that the topological and statistical properties of the IGM can be reconstructed accurately even on small scales that correspond to mean line of sight separations of the order of 1 arcmin. In this paper, we explore how well maps can be made on much larger scales, corresponding to mean quasar angular separations of tens of arc minutes. Our work, focusing on large-scale structure, is therefore complementary to the studies where individual structures like protoclusters and voids in the Lyman-α forest have been identified (Stark et al. 2014(Stark et al. , 2015. In a companion paper (Ozbek et al. 2015) we present maps made from the BOSS final data release (Alam et al. 2015), which contains such a sample of quasars covering a wide area.
The Lyman-α forest provides a means of studying the large scale structure in the redshift range 2 z 3.5. Before recent large-scale structure surveys such as BOSS (Dawson et al. 2013), the sky density of known background quasars in this redshift range over most of the sky was of the order of 1 per square degree (e.g., from the 2dF quasar survey (Miller et al. 2002;Outram et al. 2003;Croom et al. 2004;Miller et al. 2004), and from SDSS I and II (Schneider et al. 2002(Schneider et al. , 2003(Schneider et al. , 2005Richards et al. 2006). Except for some small areas with higher observed densities of objects (e.g., Rollinde et al. 2003) the Lyman-α forest was treated as a collection of discrete 1D individual quasar sightlines.
Recently, however, the increasing number of discovered quasars with suitable redshifts (z > 2) for ground based study has made it possible to correlate information over large scales in three dimensions. BOSS features a high QSO density of ∼ 15 deg −2 (∼ 180, 000 QSOs in the redshift range 2.15 < z < 4 over 10,000 deg 2 ). Each QSO provides Lyman-α forest information along a skewer of length ∼ 400 h −1 Mpc, and the typical mean separation for spectra in BOSS is ∼ 20 comoving h −1 Mpc (Lee et al. 2013). This is what has enabled clustering statistics of the Lyman-α forest to be measured in three dimensions as mentioned above. In the future, one can expect yet higher densities of sightlines and more precise measurements, in view of the fact that even more quasars will be available for analysis (e.g., 45 deg −2 proposed for MS-DESI ), see Table 1). In this paper, we use the quantity nLOS for the areal quasar density in observational surveys, as given in Table 1, first row. For the number of sightlines chosen at random to carry out the interpolation in the simulation cube, we define a new quantity NLOS, whose appropriate values are chosen according to nLOS values from observational surveys, as shown with the red markers in Fig. 1.
At the redshifts relevant to the work we present here (z = 2-4), most of the volume of the IGM is in photoionization equilibrium. The thermal state of the IGM is determined by photoionization heating and adiabatic cooling (Hui & Gnedin 1997). Lya forest traces weakly overdense regions at these redshifts (Bi 1993), and so we can use it to study the large scale structure in this diffuse medium. The IGM density fluctuations follow those of the total matter potential, tracing dark matter at large scales ( 0.1 h −1 Mpc) and respond to pressure smoothing on small scales (< 100 h −1 kpc) (Sargent et al. 1980;Ikeuchi 1986;Rees 1986).
The Lyman-α forest is useful for studying large scale structure at high redshift (2 < z < 4), in a fashion complementary to galaxy surveys at low redshift. The redshift range accessible for surveying large samples of the latter is limited by the surface brightness of galaxies scaling as (1 + z) −4 . The space density of tracer objects in current large-scale maps of the Universe above redshift z = 0.5 declines rapidly (e.g., in the CMASS sample of BOSS DR9 (Nuza et al. 2013), the comoving space density is ∼ 3.6×10 −4 h 3 M pc −3 at ∼ z = 0.5 but falls off rapidly with increasing redshift). Here we define large-scale maps as those which cover an appreciable fraction of the sky area. At redshifts z > 2, quasars are the only tracers that can be used to make such maps, and the space density of objects declines to the order of 10 −6 h 3 M pc −3 (Dawson et al. 2013). With a mean 3D separation of ∼ 100 h −1 Mpc or larger, the shot noise in maps made using quasars is extremely high. Using absorption lines, for example metal lines, where multiple lines are available for each quasar offers a way to increase the space density of tracers (Vikas et al. 2013;Zhu & Ménard 2013). The Lyman-α forest, being a continuous field offers a way of reducing the level of noise in maps even further.
When one has saturation in the Lyman-α forest, the Lyman-β transition should provide information due to its lower cross section, which makes it a potentially better probe at high overdensities (Shull et al. 2000) as high as 10 times the mean density at z = 2−3. Since Lyman-β absorption occurs at a lower rest wavelength (1026 A), the Lyman-α forest overlaps the Lyman-β forest, and so statistical techniques are needed to make use of the Lyman-β information (Dijkstra, Lidz & Hui 2003;Iršič et al. 2013). In principle higher order transitions could be used together with the Lyman-α forest in mapmaking also.
The structure of this paper is as follows: We first describe the properties of the cosmological hydrodynamic simulation we use and our choices of simulated data samples at redshifts 2 and 3 with varying line of sight densities and noise levels in accordance with current and future sky surveys in §2. In §3, we review Wiener filtering and how it is used to carry out the recovery the flux field of the IGM. We look at point to point comparisons of true and reconstructed fields in §4.1, and then provide correlation plots to test the validity of the reconstruction as a function of scale, followed by non-Gaussianity analysis via Kolmogorov-Smirnov tests in the following subsections. In §4.4, we look for local peaks in the true field and the reconstructed field and examine slice images to make a visual comparison between the true and reconstructed fields. Finally, we test the use of an alternative input correlation matrix in the Wiener filtering. We summarize our results in §5 and also provide closing comments.

SIMULATION
In order to evaluate the expected performance of map-making reconstruction on Lyman-α forest data from BOSS and other observational surveys, we make use of a large hydrodynamic cosmological simulation of the ΛCDM model. We use the smoothed particle hydrodynamics code P-GADGET (see Springel 2005;Di Matteo et al. 2012) to evolve a distribution of 2 × 4096 2 = 137 billion particles in a cubical periodic volume of side length 400 h −1 Mpc. The simulation cosmological parameters were h = 0.702, ΩΛ = 0.725, Ωm = 0.275, Ω b = 0.046, ns = 0.968 and σ8 = 0.82. The mass per particle was 1.19×10 7 h −1 M (gas) and 5.92×10 7 h −1 M (dark matter). A gravitational force resolution of 3.25 kpc/h comoving was used. The power spectrum of the simulation initial conditions was taken from CAMB (Lewis, Challinor & Lasenby 2000). The simulation was run with an ultraviolet background radiation field consistent with Haardt & Madau (1995). Cooling and star formation were included. However the latter used a lower density threshold than usual (for example in Springel & Hernquist 2002) so that gas particles are rapidly converted to collisionless gas particles. This was done to speed up execution of the simulation. As a result the stellar properties of galaxies in the simulation are not predicted reliably but this has no significant effect on the diffuse IGM that gives rise to the Lyman-α forest. Black hole formation and feedback from stars were also switched off in the simulation.  . The red markers show the fiducial choices for our work with simulation data. The BOSS quasar catalog indicates that the number of quasars peaks at z ∼ 2.25 and decreases rapidly at higher redshifts.

Data
We use two simulation snapshots at redshifts of z = 2 and z = 3 to generate two sets of Lyman-α spectra using information from the particle distribution (Hernquist et al. 1996). We make spectra set out on a grid with 176 2 = 30, 976 evenly spaced sightlines, resulting in 2.27 h −1 Mpc spacing. This can be compared in the line of sight direction with BOSS pixels of width ∆v = 69.02 km s −1 (Lee et al. 2013), which is ∼ 0.6 h −1 Mpc at z = 3. Each simulation sightline was generated with high resolution, 10,560 pixels, in order to resolve the thermal broadening when computing the optical depth. The spectra were then downsampled (by averaging the transmitted flux over 60 pixels) to 176 pixels. The full set of simulation data sets therefore consist of 176 3 data values each. In order to roughly approximate the noise which will be expected in observational data, we add random uncorrelated Gaussian pixel noise to the data sets with a signal to noise ratio S/N = 1 or 2 per unit simulation pixel in 176 pixels per 400 h −1 Mpc sightline. This is similar to BOSS which has S/N of order unity . If BOSS Lyman-α data were binned to 9.3 h −1 Mpc pixels, the mean S/N ratio would be 2.5. Similarly, the simulation data with S/N = 1, when binned to pixels of the same size results in an S/N ratio of 4. Therefore, we use the S/N = 1 case as a close match for the BOSS noise level, whereas the other noise level, S/N = 2, is  Table 2. Our choices of simulated data sets at redshifts 2 and 3 with different LOS density and noise levels.
given as an example with less noise. A random subset of the 176 2 sightlines was chosen, according to LOS area densities from the experiments BOSS (Alam et al. 2015), eBOSS (Raichoor et al. 2015) and MS-DESI ) (e.g., NLOS ∼ 400 LOS passing through the simulation volume for BOSS) from our simulation box to carry out the reconstruction (see Table 1 and Fig. 1). We define NLOS as the number of lines of sight chosen to reconstruct the entire volume in the simulation box, whereas nLOS (see, e.g., Table  1) denotes the total number of sightlines along the entire redshift range from the observer to quasars. Our data points derived from the simulation are optical depths, (τ = − log e F/F0), where F is the flux received at a certain location in space and F0 is the unabsorbed flux. The data are convolved with the peculiar velocity, therefore we work in redshift space. In (Caucci et al. 2008) the authors worked with the density field directly. We will however work with the flux, and our maps will be reconstructions of the three dimensional flux field, in redshift space. The relation between the gas density and the optical depth is where δ(x) ≈ ρ−ρ ρ is the density contrast, and α and A(z) are redshift dependent factors. We present our results in terms of flux contrast, δF = (F/ F ) − 1, where F is the mean transmitted flux computed from all spectra.

Simulated datasets
We have made 12 simulated data sets with different sightline densities and noise levels. These are summarised in Table 2. Sightline density choices were made to mimic those of current or future observational surveys, as shown in Fig. 1 with red markers. Some data sets have a LOS density that is even lower than that of BOSS (e.g. the data set labelled "z2 N200") but still allow an accurate recovery of the flux field, as we will see in §4. Data sets with higher sightline densities (e.g. "z2 N1000", which is comparable to that of MS-DESI) result in an even better inference of the field. Therefore, as observational surveys find more quasars, even more accurate density maps will be available with the Lyα forest.

RECONSTRUCTION
There are several methods which can be used to interpolate between the sparse absorption skewers in the Lyman-α forest. For example, recent work by Cisewski et al. (2014) used local polynomial smoothing for this purpose. The method we choose in this paper is Wiener filtering, pioneered in this context by Pichon et al. (2001), and used by Caucci et al. (2008), and Lee et al. (2014b) to make maps from simulated data, and by Lee et al. (2014a) to make the first 3 dimensional maps from observations.
We consider the values of the flux contrast in the reconstructed field to be entries in a column vector M, and the values of the flux contrast in the absorption skewer data to be entries in a column vector D. In general the entries of M will represent values on a uniform grid of voxels as we are constructing a map which covers all space within the map boundary. In our simulation tests, they will be covering the cubical simulation volume uniformly. We choose not to make use of the simulation periodic boundary conditions, in order to mimic some aspects of real data. Using Wiener filtering, the reconstructed 3D field M can be inferred from the absorption skewer data D by computing where CMD and CDD are the map-data and data-data covariance matrices and N is the diagonal noise matrix. In the present work, we assume the noise to be uncorrelated, so that the entries of N are inversely proportional to the square root of the number of pixels in each cell. The covariance matrices encode the expected correlation structure of the field. In most of our work we use the following simple form advocated by Pichon et al. (2001) and Caucci et al. (2008), where (x1 − x2) and |x 1⊥ − x 2⊥ | represent the distances between two pixels, parallel and perpendicular to the LOSs respectively, L || and L ⊥ are correlation lengths parallel and perpendicular to the LOSs, while the variance σ 2 is calculated directly from the field. The CDD covariance matrix contains correlation information between the initial data points only (the D array), whereas CMD contains information of the 3D pixel locations of the map to be inferred and the D array. In order to test how well the reconstruction works as a function of line of sight density, we make several different data samples by choosing a subset of our lines of sight at random. The areal density of the sightlines are those that correspond to some current and planned experiments, e.g. BOSS (Dawson et al. 2013) and DESI ) (see Table 1 and Fig. 1).
Our numerical code to carry out the reconstruction splits the simulation volume into "subcubes". The interpolation is then carried out separately for each subcube in parallel with the others and in the final step the results are combined to form the whole reconstructed simulation cube. In order to make the calculations more computationally tractable, we decrease the resolution of the field from 176 3 to 44 3 pixels.
We introduce a buffer volume on the edges of the subcubes, allowing them to overlap, in order to avoid edge artefacts. In our fiducial reconstruction of the simulation we use 64 subcubes overall and a buffer of 40 h −1 Mpc on each side for each subcube. Each subcube therefore has a side length of 180 h −1 Mpc, including the buffer regions. We have tested and checked that adjusting the num-ber of subcubes or changing the number of pixels does not significantly alter the results.
The code used does not take into account the periodic boundary conditions of the whole simulation box, in order to approximate the situation which will occur for observational data. This means that the reconstruction will be less accurate near the edges of the cubical simulation volume. For this reason, when choosing slice images to compare real and reconstructed fields, we choose the middle planes of the cube rather than the edges. We find by visual inspection that there is no significant difference in the quality of reconstructions when increasing the separation from the box edge by greater than 50 h −1 Mpc. We repeat some statistical calculations after truncating the cube by 50 h −1 Mpc from each edge in order to test the importance of edge artefacts.
The resolution of the maps is determined by the mean separation between quasar lines of sight: For NLOS = 200, dLOS is equal to 28.28 h −1 Mpc. For a study of the map resolution as a function of exposure time, the reader is referred to Lee et al. 2014b have shown that typical values for the correlation lengths L || and L ⊥ should be of the order of dLOS in order to avoid numerical instabilities in the matrix inversion leading to spurious structures. We smooth both true and reconstructed fields with an isotropic Gaussian filter with a standard deviation σS = 1.4 dLOS , in the latter case after carrying out the reconstruction.

ANALYSIS
After the reconstruction of the field we analyse the relationship between the true and reconstructed fields, bearing in mind that both have been smoothed, as stated above.

Scatter Plots
We first show point by point scatter plots in Fig. 2. We plot the reconstructed flux contrast, δrecon against the true flux contrast, δorig. The results of the point to point comparison of the fields are summarised in Table 3. Throughout the paper, we use "original field" (or "true field") with the meaning that we keep all of the Lyman-α skewers in the cube, while "recovered field" or "reconstructed field" means the flux field inferred with the given LOS density with the quasars located at redshift 2 or 3. In the top left panel of Fig. 2, we show the results for the the NLOS = 200 dataset with no noise. We can see the slope of the relation between δrecon and δorig is biased (this was also found by Lee et al. 2014a), in the sense that the recovered field has more contrast than the original field. After fitting a linear regression we find that the slope is 1.73, whereas the y-intercept is consistent with zero. This bias depends on the interplay between the Wiener filter smoothing scales and fluctuations in the field that are missed in the sparse sampling. The bias is larger when the number density of sightlines is low (compare the top left panel of Fig. 2 which has a fitted slope of 1.73 and a sightline density 5 times less than the bottom left panel, which has a fitted slope of 1.35). Any correction for this bias is likely to be empirical, and therefore in the rest of our analysis we apply the simplest correction, by renormalizing the δrecon according to the slope of the regression. In order to quantitatively test the quality of the reconstruction, we compute the error by calculating the ratio of the root mean square (RMS) of the pixel by pixel difference to the RMS of the true field (δorig) which only includes 95 per cent of the true pixels (±2σ from the mean, which is 0), therefore avoiding outlier points: From Fig. 2, we can see from the top row of panels that the addition of noise to the input field does affect the reconstruction. The RMS error (after bias correction) is 20.4 per cent for the noiseless case and 37.0 per cent and 57.3 per cent for the cases with S/N=2 and S/N=1 respectively in the pixels of size 9.30 h −1 Mpc used in our analysis. We remind the reader that the BOSS Lyα forest data when rebinned in this way has a mean S/N ratio of 2.5. Therefore, we estimate that the RMS error for the reconstructed BOSS data will be 36.7 per cent.
We also provide the Pearson coefficient (r) as a measure of the linear correlation between δorig and δrecon. Total positive correlation between the original field and the reconstruction would correspond to r = 1, while no correlation would be r = 0 and total negative correlation would be r = −1. At redshift z = 2 with NLOS = 200, r = 0.783, and it is even higher when the number of sightlines is increased, as expected.
Increasing the number of sightlines, as shown in the bottom left panel of Fig. 2 allows the smoothing scale to be reduced and the resolution of finer features in the flux contrast field. The RMS flux contrast fluctuations increase to 0.0193 for this sample (z2 N1000, see Table 3) and the percentage error on the reconstruction stays approximately the same as the lower resolution reconstruction in the top left panel.
Finally we show results for the higher redshift, z = 3 in the bottom right panel of Fig. 2. For the sample z3 N200, we use the same number of sightlines, as the top left (z = 2) panel, but the RMS accuracy of the reconstruction is lower by a factor of 0.75. The quality of the reconstruction is low for the sample z3 N60, as can be seen in Fig. 10. However, the percentage error in Table 3 is comparable to those from other samples due to the fact that the dynamic range in the true field is low because of the much higher σS value for that sample. The flux contrast fluctuations are larger at the higher redshift, because of the greater overall optical depth in the Lyα forest, but this does not translate into a better reconstruction.

Cross Correlation
In the previous section we have seen how the true and reconstructed fields compare on a point-by-point basis for one specific value of the smoothing scale, given by σ = √ 2 dLOS . We can also study how similar the fields are at different scales. Instead of using different filter scales, we make use of the correlation functions of the field, measuring the correlation between points separated by a distance r. We compute the auto-correlation of the true field, the cross correlation and the standardized cross correlation between the true and the recovered fields. In general, the correlation function of two fields 1 and 2 is defined by: For the auto-correlation, 1 and 2 represent the same fields. The standardized cross correlation is defined by: where ξ12 is the cross-correlatin and ξ11 and ξ22 are the autocorrelations. C12 enables us to quantify the accuracy of the reconstruction as a function of scale, with a value of unity indicating perfect fidelity. It should be noted that we do not expect good agreement for scales smaller than the fiducial smoothing length, the natural resolution of the map. According to Equation 4, the fields with NLOS = 200 are smoothed with σS = 39.6 h −1 Mpc, while those with NLOS = 1000 are smoothed with σS = 17.7 h −1 Mpc.
If we look at the case of low sampling density, NLOS = 200, we see that C12 can be as low as 0.45 on scales which are approximately 1.5-2 times the smoothing scale. Repeating the same analysis for unsmoothed fields with this LOS density does give better agreement at low scales -clearly, some structure is erased due to smoothing. For larger values of NLOS, on the other hand, we get good agreement despite the smoothing (Fig. 3). After smoothing both fields with the fiducial smoothing length, the recovery improves with distance until very large distances, where it decreases again. This decrease is due to edge effects, as on large scales, much of the volume is close to an edge. In panel (d) of Fig. 3 we restrict the measurements of ξ12 to the volume of the simulation cube left after eliminating all regions within 50 h −1 Mpc of an edge. Comparing with panel (c) of Fig. 3, we can see that the agreement on small scales has also been improved, showing the positive effects of getting rid of the edge artefacts. In a large survey such as BOSS which spans a contiguous volume of several gigaparsecs, most of the volume will be much further from an edge than 50 h −1 Mpc, so that edge effects should be a small issue. When the edge effects have been removed, we can see that the C12 measurement is close to 1 for all scales greater than the smoothing filter scale at z = 2, indicating essentially perfect statistical agreement. At redshift z = 3, C12 is never greater than 0.8, which may indicate that the less evolved structures at higher redshift make accurate reconstruction more difficult.

Non-Gaussianity and the Flux Probability Distribution
The density field probed by the Lyα forest is expected to be in the mildly non-linear regime. When smoothed on large scales, which we necessarily must do in order to construct our interpolated maps, we expect that the flux probability distribution should be quite close to Gaussian. Indeed this Gaussian assumption underlies the reconstruction carried out with the Wiener filter in Equation 3. It is therefore of interest to compare the reconstructed and true flux probability density functions with each other and with a normal distribution.
In order to look for deviations in the flux pdfs, we use the Kolmogorov-Smirnov (KS) test. We first compute the mean and the standard deviation σ of a Gaussian fit to the true field. We construct a cumulative probability distribution from this and compare it to the cumulative probability distributions of the true field and the reconstructions, for various NLOS values, redshifts and levels of noise.
With the KS test, we compute quantitative measures of the similarity of the flux pdfs to normal distributions. The test statistic (D value) is the maximum of the difference in the cumulative distribution functions of the particular field being tested and the Gaussian. The closer this value is to 0, the more likely it is that the data sets are have been drawn from the same distribution. Furthermore, the p value, which is computed from the test statistic, represents the significance level threshold below which the null hypothesis (that the data sets come from the same distribution) will be accepted.
When computing the flux pdfs, we would like the data points to be as independent as possible, and so our data points should at least be separated by distances greater than the smoothing scale, because smoothing would correlate the measurements. Because of this, we downsample each data set, picking only a 5 3 grid of values (i.e., data points separated by 80 h −1 Mpc in each direction in the 400 h −1 Mpc volume). Our KS test results are shown in Table 4, for different line of sight densities, redshifts and signal to noise ratios.
In all cases for the true field we find high p values and low D values, which means that the original field was approximately Gaussian to start with. The recovered field also shows the same property as well in Table 4, for recovery from data samples with no added noise. When noise was added, however, the reconstructed maps became significantly non-Gaussian, with the p value decreasing as the signal to noise ratio decreased. Furthermore, the fact that the D and p values between samples with different LOS density are significantly different can be attributed to the fact that the smoothing filter size depends on the LOS density itself. As we have mentioned above, the noisiest data, which has S/N=1, is significantly worse than the majority of BOSS data, for example, but the effect of noise in changing the pdf shape of the reconstructed field should still be borne in mind in an analysis of observed data.

Peaks In the Density Field
Searching for local maxima in the reconstructed flux density field offers one means of defining objects and finding them. Such peaks are likely to correspond to the locations of forming clusters or superclusters of galaxies. The properties of these density maxima can be used to constrain the cosmological model (Bardeen et al. 1986;Croft & Gaztanaga 1998;De & Croft 2007, 2010. It is therefore of interest to compare the peaks of the reconstructed flux density field with those in the true flux density field in the simulation. We search for density peaks in the three dimensional volume of the simulation. Density and flux are inversely related, therefore  Table 4. Kolmogorov-Smirnov test results giving the probability (p value) of the flux pdfs of the real and the reconstructed data being drawn from a Gaussian distribution. We show results for different numbers of quasar sightlines through our simulation volume, N LOS , redshifts and signal to noise ratios.
we identify a simulation 3D pixel as a local peak if its flux value is the smallest amongst the 26 neighbouring 3D pixels surrounding it. As expected, we find that the number density of local peaks is strongly dependent on the smoothing filter size. We find that for a filter size of 39.6 h −1 Mpc, appropriate for NLOS = 200, we find 9 local peaks in the simulation volume at z = 2 (Fig. 4a), and for a filter size of 17.8 h −1 Mpc, appropriate for NLOS = 1000, we find 87 peaks (Fig. 4c), where both of these figures are for the true field. Furthermore, z = 2 and z = 3 data samples have very similar true peak locations for the same filter size. When noise is added to mimic observational data, we discover 8 local peaks for the real simulation flux field, while the reconstructed field contains 10 local peaks for the data set z2 N200 SN1. In Fig. 4 we show a comparison between the true field local peaks and those of the reconstructed field. The number of peaks in both cases match exactly for the three different combinations of NLOS and redshift shown. The positions of the peaks are visually a reasonable match, with better agreement for z = 2 than z = 3. The structures traced out in the plot with NLOS = 1000 by the reconstructed peaks do seem visually to trace out those in the real peaks.
There is not a one to one correspondence however. We can quantify the level of agreement by counting the number of peaks in the true field which have a peak in the reconstructed field within one smoothing length. This is 33.3 per cent for z = 2, NLOS = 200. The expectation from randomly positioned peaks with the same number density (computed using 1000 Monte Carlo trials) is 11.5 per cent. This means that the reconstruction is a factor of 2.89 better than random. The equivalent measures for z = 3, NLOS = 200 and z = 2, NLOS = 1000 are 11.1 per cent and 32.1 per cent peaks within 1 smoothing length respectively and factors of 1.13 and 10.4 better than random.
This type of analysis could be potentially extended to look at the structures that are enclosed within isodensity contours. This would reveal the morphology of the IGM on large scales. At the scales we are probing here (smoothing scales > 10 h −1 Mpc), sheet and filament-like topologies are relatively difficult to see (as we shall see in our visualisations in the next section). On smaller scales, these characteristics are readily apparent in simulated maps (e.g. Pichon et al. 2001) and the first maps made from observational data with these techniques (Lee et al. 2014a). The most straightforward cosmological constraints will come from the peak density, and the reconstruction technique does very well: We get perfect agreement between the real and reconstructed fields for noiseless comparison, and for the noisy case (with z2 N200 SN1), the number of peaks agrees within 25 per cent.

Slice Image Comparisons
We now turn to a visual comparison of the structures in the real and reconstructed maps. The three-dimensional datacubes have z axes oriented parallel to the line of sight, and x and y axes perpendicular to it. The sampling of pixels in a mock data set is therefore different depending on the plot orientation, and this could influence the recovery of structure. We therefore show two orientations for each plot, one in the y − z plane (an "x" slice) and one in the x − y plane (a "z" slice). In our plots we show the flux contrast in a slice of thickness one grid cell. As our volumes are 44 cells on a side, this corresponds to a thickness 400/43 = 9 h −1 Mpc. We have also seen in Section 4.1 that there is a bias in the reconstructed field which leads it to have higher contrast. Changes in NLOS and adjusting the correlation lengths do not alter this and so we follow (Lee et al. 2014a) in applying a bias correction before visualizing the fields.  In Figs 5 through 11 we present image slices through the simulation volume. The images show the flux contrast δF = (F/ F )− 1, which means that low values correspond to high values of the matter density. The red colour shows these higher density regions and white those of lower density. The image slices are taken from the centre of the cube in directions parallel and perpendicular to the LOSs. As explained in Section 3, the motivation for choosing to display slices through the centre of the cube is because we are not using information about the periodic boundary conditions in the simulation when carrying out the reconstruction. The edges of each image slice will therefore give an idea of how well the reconstruction would succeed at the edges of a survey volume. We have checked other random slices and verified that the reconstruction re- covers the general features of the original field, even when close to edges of the simulation volume. In Fig. 5 we can see the results for our lowest density of sightlines (we have NLOS = 200), at redshift z = 2. We can see that the general morphology of the field is recognizably similar in the true and reconstructed maps. In detail, the maps have some differences, but the maxima, minima and their gross shapes are fairly well reproduced, and one could therefore expect that observational data from the BOSS survey (which has approximately this number density of quasars) would yield visually quite accurate maps of the large scale structure, at least when smoothed on the relevant filter scale (a filter of 39.6 h −1 Mpc was used here.).
The top and bottom rows of Fig. 5 show results for slices parallel and perpendicular to the line of sight. We see no obvious difference in the fidelity of reconstruction for each, and there is no obvious sign of the discrete sampling of the field by pixels and sightlines (which is different for the top and bottom rows). General features of the field are recovered well, especially for mildly dense regions.
Figs 6 and 7 demonstrate the effect of adding uncorrelated Gaussian noise to the flux field in order to better mimic observational data, and how the fidelity of the reconstruction changes when noise is introduced. Noise levels (S/N = 1 or 2) indicate the amount of noise for a simulation pixel (∼ 0.76 h −1 Mpc) wide. Since our pixels are rebinned to ∼ 9 h −1 Mpc, the added noise is reduced by a factor of 3.4. Hence, the difference between the true fields before and after adding noise is small. However, the reconstruc- tion is sensitive to the amount of noise, therefore the fidelity of the noisy reconstruction is noticeably worse, especially for overdense and underdense regions. We increase the density of the sightlines in Figs 8 and 9. As a result, the quality of the reconstruction is visually better. Although the reconstruction code does not take into account the periodic boundary conditions of the simulation, the fields are comparable even at the edges. This is likely due to lower smoothing levels, as the smoothing level scales inversely with NLOS.
Fluctuations in the flux field are greater at higher redshifts. As Figs 10 and 11 clearly show, this results in a decrease in the fidelity of the reconstruction. It is obvious that at redshift z = 3, the LOS density of NLOS = 60 is not good enough to yield a comparable reconstructed field. For observational data, we naturally expect a better map at z = 2 than at z = 3, as the LOS density is higher at z = 2. In this study, although the LOS density is set to be the same at both redshifts, we get a better recovery of the field at z = 2.
Furthermore, in Fig. 13, we show one-dimensional visual comparisons along four lines of sight chosen at random using the data set z2 N 200, whose source density matches the areal LOS density of BOSS, and observe that the recovered skewers capture general features of the original ones.
The recovery is accurate for scales larger than ∼ 1.4 dLOS , as found in (Caucci et al. 2008), especially for mildly dense regions (standardized correlation plots). Due to the isotropic nature of the recovery and the smoothing, we do not notice any significant statistical difference between the directions parallel to and perpendicular to LOSs. As NLOS increases, naturally, the recovery gets substantially better. This means that with future experiments like eBOSS and MS-DESI, which have higher areal density of LOSs, a very accurate large map of the IGM can be generated. Adding noise to our data (pixel by pixel) and carrying out the recovery is an important step in order to better simulate real data from experiments. It is clearly seen from the figures that adding noise makes the recovery of overdense and underdense regions significantly worse. Furthermore, our results with the data set at z = 3 are significantly worse than the other data set at z = 2. and (d) show reconstructed field slices. The smoothed reconstructed field recovers the general features of the simulation, although the quality of the reconstruction is lower than that with the z = 2 data set.

Using Observed Correlations in the Covariance Matrix
Along with Caucci et al. 2008, andLee et al., 2014, we use a simple Gaussian form for the correlation function which appears in the Wiener interpolation covariance matrix (Equation 3). This is motivated by simplicity, and the fact that it is well behaved numerically at large separations. One might expect covariance matrices computed from the actual correlation functions of the field to give more accurate reconstruction results, however, and we now test this.
In (Slosar et al. 2011), the three-dimensional correlation function of the absorption in the Lyman-α forest was measured for the first time. The measurement was extended to greater than 100 h −1 Mpc scales by (Busca et al. 2013) and (Slosar et al. 2013). We use this measurement of correlation function to construct a correlation matrix instead of the Gaussian covariances we have used (Equation 3).
The correlation function measured from the observational Lyα forest data is anisotropic because of redshift distortions. We construct the correlation matrix not from the observational data results of (Slosar et al. 2011), but from the linear theory CDM model consistent with the data. This redshift space model fit is given by Equations 4.5 -4.13 of (Slosar et al. 2011). We use these equations, along with the linear theory correlation function from Section 2, and the following parameters: bias factor b = 0.2, and redshift distortion factor β = 1.5 to compute ξF (r ⊥ , r ), the flux correlation function for line of sight separation r and transverse separation r ⊥ . The Wiener covariance (replacing Equation 3) is then given by . Slice plots from the middle of the cube are shown above from a data sample at redshift z = 2 with N LOS = 1000. We use the original Gaussian correlation function as well as a correlation function measured from observations to recover the field. The top row shows results for the middle of the cube, perpendicular to LOSs. The bottom row has plots in the middle of the cube, parallel to LOSs. The recovery of the field using the Gaussian correlation function yields better results than using the CDM correlation function obtained from observations.
where r = (x1 − x2) and r ⊥ = |x 1⊥ − x 2⊥ |. After reconstructing the simulation field using the CDM fit to the (Slosar et al. 2011) results in the covariance matrix, we compare the results to our fiducial reconstruction technique (Fig. 12). We find that the recovery of the field with the fiducial (Gaussian) correlation functions yields slightly better results than with the CDM correlation function. Although the dynamic range with the Gaussian correlation seems to be slightly higher, general features of the original field are recovered better. For example, for the data set z2 N1000, instead of our original result of the RMS percentage error 17.2, we find 20.3 with the CDM correlation function, which is significantly worse. It is worth noting that the actual correlation function in the simulations is probably not the same as it is estimated by Slosar et al. (2011).

SUMMARY AND CONCLUSIONS
Using Wiener interpolation, we reconstruct the entire simulation box with a subset of the Lyman-alpha absorption skewers chosen randomly. This subset of skewers, NLOS, sets a natural resolution of our maps. The number of the skewers chosen at random is decided by matching it with the areal LOS density of current and future spectroscopic surveys such as BOSS and MS-DESI. Using the Lya forest with this method, one can make maps of the large scale structure at high redshifts (2 < z < 3.5).
The standardized cross correlation plot (Fig. 3, panel(c)) indi-cates that the reconstruction is much better at z = 2 than at z = 3 using BOSS areal LOS densities. Naturally, the fidelity of the reconstruction is better as NLOS is increased. Truncating the cube 50 h −1 Mpc from each edge, in order to remove the edge artefacts resulting from periodic boundary conditions of the simulation and to better mimic observational data, yields significantly better reconstruction (Fig. 3, panel(d)). We find that the data set at z = 2 yields clearly better results than the one at z = 3 for the simulation, even with the same NLOS. This is most easily understood in terms of the growth of structure through gravitational instability between z = 3 and z = 2. For observational surveys, in view of the fact that the areal LOS density is also much greater at z = 2 than at z = 3, one naturally expects that the large scale structure map will be significantly better at lower redshifts.
The overall bias seen in point to point flux values in real and reconstructed fields is an issue which does not have an easy explanation. Adjusting NLOS, the correlation lengths and the buffer length does not change the situation, but using an empirical bias correction allows the fields to be well-reconstructed.
In the high redshift range covered by the Lya forest, the IGM density field is expected to be in the mildly non-linear regime, therefore, we look for non-Gaussianity in the probability density functions of our reconstructed maps. While its behaviour is nearly Gaussian for noiseless data samples, it becomes less Gaussian as the noise level is increased, as Table 4 indicates.
Since the smoothing levels used in this study are greater than 10 h −1 Mpc, it is not possible to see the filamentary structure in  the IGM topology. We remind the reader that there are no widefield galaxy surveys that can detect the topology of the IGM at z > 2, as it is increasingly expensive to detect galaxies at higher redshifts to reach a high source density, even with 8-10 m telescopes (Le Fèvre et al. 2013). However, searching for local peaks allows us to discover the potential locations of the superclusters, which can be cross-correlated with galaxy surveys. Both the number of the local peaks and their locations are reproduced reasonably well by the interpolation.
Slice images allow a visual comparison between the original and the reconstructed fields. General features of the flux field are well reproduced by the interpolation, especially for mildly overdense regions. Using a correlation matrix derived from a CDM-fit to observational data instead of the simple Gaussian correlation matrix used in our fiducial Wiener filtering leads to a slightly worse recovery of the field.
As an improvement, the isotropic smoothing of both fields can be altered, as one does not necessarily expect the same statistics parallel and perpendicular to the LOSs. Furthermore, as future sur-veys like eBOSS and MS-DESI discover more quasars, the fidelity of the large scale structure maps will improve. For example, in order to reach resolutions in the sub 10 h −1 Mpc regime at z ∼ 2 to study the IGM filamentary structure, QSO densities of over 100 deg −2 will be necessary.