MAMMOTH-Subaru. V. Effects of Cosmic Variance on Lyα Luminosity Functions at z = 2.2–2.3

Cosmic variance introduces significant uncertainties into galaxy number density properties when surveying the high-redshift Universe with a small volume. Such uncertainties produce the field-to-field variance σ g of galaxy numbers in observational astronomy, which significantly affects the luminosity function (LF) measurement of Lyα emitters (LAEs). For most previous Lyα LF studies, σ g is often adopted from predictions by cosmological simulations, but barely confirmed by observations. Measuring cosmic variance requires a huge sample over a large volume, exceeding the capabilities of most astronomical instruments. In this study, we demonstrate an observational approach for measuring the cosmic variance contribution for z ≈ 2.2 Lyα LFs. The LAE candidates are observed using the narrow band and broad band of the Subaru/Hyper Suprime-Cam in eight independent fields, making the total survey area ≃11.62 deg2 and a comoving volume of ≃8.71 × 106 Mpc3. We report a best-fit Schechter function with parameters α = −1.75 (fixed), LLyα*=5.95−0.96+1.22×1042 erg s−1, and ϕLyα*=5.26−1.27+1.65 × 10−4 Mpc−3 for the overall Lyα LFs. After clipping out the regions that may bias the cosmic variance measurements, we calculate σ g by sampling LAEs within multiple pointings on the field image. We investigate the relation between σ g and survey volume V, and fit a simple power-law σg=k×V105Mpc3β . We find best-fit values of −1.399−0.156+0.160 for β and 1.249−0.193+0.213 for k. We compare our measurements with predictions from simulations and find that the cosmic variance of LAEs is likely larger than that of general star-forming galaxies.


INTRODUCTION
The cosmic structure formation arises from fluctuations during the inflationary epoch (Guth 1981) and introduces inhomogeneity to the Universe except on large scales.Due to the presence of the scale-dependent inhomogeneity, observations on the number density and Luminosity Function (LF) within a small survey volume have an inevitable level of uncertainty (e.g.Szapudi & Colombi 1996;Trenti & Stiavelli 2008;Robertson 2010a), known as the cosmic variance.Note cosmic variance is referred to as a sample variance due to different survey sizes within our Universe, instead of arising from an ensemble across multiple Universes.This source of error dominates uncertainties when the survey volume is smaller than the typical clustering scale, or at high redshift where galaxies tend to be more clustered than dark matter (Kauffmann et al. 1999;Coil et al. 2004).Therefore, cosmic variance is significant in observational studies of luminosity functions (LF) for high redshift galaxies, where the number density is expressed as a function of the emission line luminosity (For a detailed approach, see Gronke et al. 2015).
Although many works have been done for LAE surveys at z ∼ 2 − 3 (e.g.Matthee et al. 2017;Sobral et al. 2018;Spinoso et al. 2020), significant disagreements existed between different studies (e.g.Ciardullo et al. (2012) and e.g.Konno et al. (2016)), and it has been suggested by many studies (e.g.Sobral et al. 2015Sobral et al. , 2017) ) that strong cosmic variance can be a possible explanation.Thus, quantifying the cosmic variance is required to understand the formation of the cosmic structures.
The fluctuations in the galaxy number densities caused by cosmic variance are sometimes referred to as the field-to-field variations and are often denoted as σ g where g stands for galaxies.The value of this σ g can vary with galaxy types, redshift, survey volumes, and sometimes the survey geometry (e.g.Newman & Davis 2002;Moster et al. 2010;Robertson 2010b).Theoretically, the field-to-field variation for a certain type of galaxy with a similar survey geometry (e.g.'pencil beam' shape for narrowband imaging surveys) at a given redshift only depends on the survey volume.In general, measuring the field-to-field variance requires a vast sampling, however, due to the relatively small sample size, it is difficult to estimate the cosmic variance with a sufficient level of accuracy for previous studies.Therefore, the common approach is to use the field-to-field variation predicted from cosmological simulations (e.g.Newman & Davis 2002;Trenti & Stiavelli 2008;Moster et al. 2010).Nevertheless, such a field-to-field variance has not been confirmed observationally.
In this study, we provide an observational approach for investigating the cosmic variance through Lyα LFs of z ∼ 2.2 LAEs.We have developed a simple relation between σ g and the survey volume for z ∼ 2.2 narrowband imaging surveys.The data are obtained from the Subaru/Hyper Suprime-Cam (HSC, Miyazaki et al. (2018)) with a wide Field-of-View (FoV) of 1.5 deg in diameter, which gives us a large survey volume available for Monte-Carlo sampling.Eight fields targeting the candidates of the MApping the Most Massive Overdensity Through Hydrogen (MAMMOTH) project (Cai et al. 2016(Cai et al. , 2017b,a) ,a) are used, with each field covered by one HSC FoV.
This paper is organized as follows: We summarised the information on observational data and target fields, the data processing, and the selection criteria of LAEs in Section 2. A detailed explanation for computing the Lyα LFs is described in Section 3. In Section 4, we demonstrate the method and results for the field-to-field variance estimations using Lyα LFs.Finally, a summary of this study is listed in Section 5. Throughout this paper, we adopt the AB magnitude system (Oke & Gunn 1983) and a flat ΛCDM cosmology with Ω m = 0.3, Ω Λ = 0.7, and h = 0.70.
We obtained our observations using the gigantic mosaic CCD camera HSC, with a wide Field-of-View (FoV) of 1.5 deg in diameter, and a pixel scale of 0.168".Narrowband imagings are carried out using two NB filters, NB387 (λ c = 3862 Å, FWHM = 56 Å) and NB400 (λ 0 = 4003 Å, FWHM = 92 Å).The corresponding redshifts of detecting Lyα emissions are z = 2.18 ± 0.02 and z = 2.29 ± 0.04 respectively.We also use the gband to evaluate the continuum of the detected objects. Figure 1 shows the transmission curve of these 3 filters used, which takes the transmittance accounting in CCD quantum efficiency, dewar window, and the Primary Focus Unit and the reflectivity of the Prime Mirror into account.The imaging data are reduced with the HSC pipeline, hscP ipe (Bosch et al. 2018;Aihara et al. 2019), see Liang et al. (2021) and Zhang et al. in prep.for details of data reduction.Due to the poor quality of the J0210 NB387 data, NB387 and g-band data of J0210 field are reduced with two different versions of hscP ipe, which results in a slight difference in the final image depth, as shown in table 1.The final catalog of J0210 is a combination of the two.and 6 are the 5σ limiting magnitudes measured within a 1.7" (2.5" for J0210 image) aperture for the final stacked NB387/NB400 image and the PSF-matched g-band, respectively.Data of J0210 field is reduced using two versions of hscP ipe, whose depths are also indicated in the first two rows.

Photometric processing
Object detection and the following photometry are performed using the SExtractor Bertin & Arnouts (1996).The g-band and NB387/400 images are PSF matched by convolving a proper Gaussian kernel in each field (Liang et al. (2021), Zhang et al. in prep.).The detection threshold is set as 15 continuous pixels over the 1.2σ sky background.We also apply masks to the pixels with low S/N signals or saturated by bright stars and artifacts during object detection.In addition, the sky background root-mean-square map is used as the weighting map in SExtractor, to minimize the influence of the fluctuations in the image depth for each field.The aperture diameters for photometry are 15 pixels (∼ 2.5") for J0210 field and 10 pixels (∼ 1.7") for the other 7 fields.
We estimate the total magnitude using AUTO-MAG, which uses the automatically determined elliptical aperture for Kron photometry by SExtractor Bernardeau & Kofman (1995).Both AUTO-MAG and APER-MAG output by SExtractor are used in the selection crite- ria described in the next section.Furthermore, we also replace the g-band magnitudes with the corresponding 2σ limiting magnitudes, for the objects with a g-band fainter than the 2σ limit.The 5σ limiting magnitudes within the 1.7" (2.5" for J0210) aperture of the final stacked g-band and NB387/400 images are listed in table 1.The depths for the NB400 images are about 1 mag deeper than that for the NB387 images in general, the seeing and the depth of J0210 are relatively poorer than all other fields.

LAE Selection
We use the narrowband color excesses in the broadband to select out our LAE samples, this approach is widely used in previous studies (e.g.Guaita et al. 2010;Ouchi et al. 2008;Konno et al. 2016).Most of these studies use multi-wavelength data and thus multiple color excesses to define their selection criteria (e.g.Uband), while we only use one broadband (HSC g-band) to estimate the continuum.However, Liang et al. (2021) have proved that it is sufficient enough for the z∼ 2.2 LAE selection.
We assume the LAE spectrum has a Gaussian-like redshifted Lyα emission with a rest Equivalent Width (EW) of 20 Å, and a continuum that follows a simple power law of f λ = λ β .We also take IGM absorption (Madau 1995;Madau & Dickinson 2014) into account when calculating the observed magnitude (Inoue et al. 2014).In addition, we find that the value of g -NB calculated from APER-MAG is not necessarily the same as that from AUTO-MAG, as shown in figure 2 (The differences between g-NB values calculated from AUTO-MAG and that from APER-MAG are not zero for every object).As a result, for some objects with APER-MAG values passing the first three criteria, its g -NB value calculated from AUTO-MAG can still be smaller than the corresponding EW 0 = 20 Å values, or even smaller than 0. Since we use AUTO-MAG to estimate the total magnitude, this could be a problem when calculating the Lyα luminosity.To address this issue, we also apply the criteria for assessing the color excess and color error to the AUTO-MAG, so that (g -NB) AUTO is also above the required values.
For NB387 data, the selection criteria are defined as: 20.5 < NB387 APER NB lim,5σ For NB400 data, the selection criteria are defined as: Here, The first line defines the detection limits using APER-MAG, where the lower limits are set to avoid saturations, and the upper limits are set as the corresponding 5σ limiting magnitudes to ensure the reliability of the source detections.The color excesses are equivalent to setting EW 0 > 20 Å, which are ∼ 0.3 for NB387 filter and ∼ 0.4 for NB400 filter.3σ(NB) corresponds to the color error, and is defined following Shibuya et al. (2018): which is designed to reject false selections of faint contaminants passing the criteria due to statistical fluctuations.We plot the LAE selections for each field in figure 2.Moreover, the NB387 data are re-selected from the original catalog of Liang et al. (2021), where a 2σ color error and APER-MAG-only criteria are used.We make the selection criteria consistent across two filters in this study to investigate the field-to-field variations.For J0210, two versions of data reductions yield two different image depths, thus both criteria are overplotted for clarification.Furthermore, the g-band magnitudes are replaced as the corresponding 2σ limiting magnitude, for any object with a g-band fainter than the respective 2σ limit.The g -NB values for these objects are shown as a lower limit in the diagrams.

Detection completeness
The detection completeness f det is estimated by Monte Carlo simulations with the GALSIM package (Rowe et al. 2015), following the procedure described in previous studies (e.g.Konno et al. 2016Konno et al. , 2018;;Itoh et al. 2018).Mock galaxies are constructed according to 14 specified NB387/400 magnitude bins with centers ranging from 20.25 to 26.75 with a bin width ∆m = 0.5.Source detection and photometry of the mocks are performed using SExtractor 2.25.0 (Bertin & Arnouts 1996), in the same manner as our sample selection (See Section 2).The detection completeness f det is defined as the ratio between the total number of mocks detected and that of mocks created (i.e.500 for each bin), as a discrete function of the NB387/400 magnitudes.We repeat the process 10 times for each field with varying random seeds and obtain the final value by taking the average.The results are displayed in figure 3. Values of f det are typically above 90% for the bright-end sources and fall rapidly to zero beyond the 5σ limiting magnitudes.

Contamination
Possible contaminants that can pass our selection criteria are mainly considered to be: [OII] λ3727 emitters at z ∼ 0.04, CIII] λ1909 emitters at z ∼ 1, and CIV λ1549 emitters at z ∼ 1.5.For the [OII] emitters, the corresponding survey volume is about three orders of magnitude smaller than that of LAEs at z ∼ 2.2, together with the [OII] LFs at z ∼ 0.1 investigated by Ciardullo et al. (2013), we conclude that the amount of low-z [OII] emitters is negligible for our sample.For CIII] and CIV contaminants, the rest-frame equivalent width would be larger than 30 Å, much larger than that of star-forming galaxies.To investigate the portion of these CIII] and CIV contaminants in our sample, we perform a spectroscopic follow-up for our LAEs.
We present an estimate of the overall contamination rate using spectroscopic observations with Magellan/IMACS on September 29 and 30th, 2022.The Magellan/IMACS is set in the multi-slit spectroscopy mode with the f/2 camera, with a field of view of 12 arcsec in radius.Slits have a width of 1.2 arcsec and a length of 8.0 arcsec.We observe one pointing in each of J0210 and J0222 fields, and two pointings in the J0240 field.The on-source exposure times are 7500s and 6000s for NB387 and NB400 LAEs, respectively.
The spectra cover wavelength between ∼ 3600 and 5700 Å with a spectral resolution of ≈ 7 Å.Since the spectral resolution is not small enough to resolve the doublet of low-z [OII] contaminants, we only use the spectra with ≥ two emission lines when estimating the contamination rate.After data reduction, we obtain 120 spectra with detected emission lines at the expected wavelength of Lyα.22 spectra are confirmed LAEs at z = 2.2 with ≥ 2 emission lines, and 2 spectra are foreground objects.Therefore, we estimated that contamination rate is about 2/(2+22)=8%.The narrowband magnitudes of the 24 objects with ≥ 2 emission lines are 21.3-25.3mag.

Lyα luminosity functions
We derive the Lyα LFs in the same manner as many previous studies (e.g.Ouchi et al. 2008;Itoh et al. 2018).We calculate Lyα luminosities using the total magnitude of NB387/400 and g-band, estimated by the AUTO-MAG.We calculate the Lyα line flux (f line ) and the rest-frame UV continuum flux (f c ) from the NB387/400 and g-band magnitudes (m NB and m BB ) using the following equation: Here, T NB and T BB are the transmittance curve in frequency space of the NB and BB filters, respectively.We take the same assumption as Itoh et al. (2018) that f line is a δ-function and f c is a constant.
The volume number densities φ(logL) in a Lyα luminosity bin of [logL − dlogL/2, logL + dlogL/2] is defined as following: where the sum is taken over all LAEs in the specified bin.f i,det is the detection completeness for the corresponding NB387/NB400 magnitude of the ith LAE, and V eff is the effective survey volume assuming a top-hat filter transmittance.We set the bin width to be dlogL = 0.1dex, which is the same as that of Konno et al. (2016).Note that the actual NB387 and NB400 filter transmission curves are not perfect top-hat.Sobral et al. (2017) investigated the incompleteness due to this effect and found the correction factor range from 0.97 in the faintest bin to 1.3 in the brightest bin.This study mainly focuses on the faint end and excludes bright end as described below, therefore, we take the assumption that the filter transmittance shapes of NB387 and NB400 filters are top-hat.
The left two panels of figure 4 show LFs derived from 8 fields, as well as the overall LFs of NB387 and NB400.The error bars include the Poisson uncertainties calculated from Gehrels (1986), as well as the field-to-field variance estimated from our measurement in section 4. The values are 0.063 for NB387 and 0.034 for NB400 (See section 4.2 for details).We include the field-tofield variance in this section.The low luminosity point has been excluded for low completeness.Note that we choose the value of logL > 42.3 for all fields for keeping the fitting process consistent.Using the minimum χ 2 fitting, the LFs are fitted to a Schechter function (Schechter 1976), defined as: Here, Schechter parameters α, L * , and φ * characterize the faint-end slope, luminosity, and density respectively (Schechter 1976).Konno et al. (2016) surveyed the Lya LFs for z∼2.2 LAEs using the NB387 data from Subaru/Suprime Cam (Miyazaki et al. 2014) in five fields and fitted three Schechter parameters simultaneously as free parameters.They reported a faint-end slope α of −1.75 +0.10 −0.09 , since this value is reasonably wellconstrained (Konno et al. 2016), we fix our α at this value.The goal of this study is to investigate the fieldto-field variation, and thus having α fixed at the same value has the benefit of reducing differences in the expected galaxy abundance for all fittings, so that the statistical fluctuations become significant.The right panels of figure 4 show the 68% and 90% confidence levels, as well as the best-fit values of the other two parameters L * and φ * .The best-fitting results are demonstrated in table 2. We combine NB387 and NB400 data to compute the overall LF for all 8 fields (See red points and line in figure 5).Due to the redshift difference between NB387 and NB400, we argue that the difference in the luminosity function between NB387 and NB400 is due to field-tofield variation rather than redshift evolution.
An excess in number densities at the bright-end with logL 43.1 is presented in the left panels of figure 4. One possible explanation is the magnification bias (e.g.Mason et al. 2015), the luminosity of several LAEs is raised due to gravitational lensing by foreground galaxies.However, past studies (e.g.Konno et al. 2016;Sobral et al. 2017) suggest that the number density excess at the bright end is due to type-I AGN rather than magnification bias.Zhang et al. (2021) have also confirmed this conclusion using the HETDEX spectroscopic surveys.Konno et al. (2016) stated that, due to the presence of significantly smaller error bars at the faint end, the humps at bright end do not affect the fitting results significantly.We attempted to fit a Schechter function to our overall LF both with (42.3 <logL< 44.1) and without the bright-end (42.3<logL<43.1), and obtained two sets of best-fit parameters of (α = −1.75(fixed), L * = 6.23 ).We thus conclude that the bright-end hump does not affect our fitting significantly.We chose to perform curve-fitting without the bright-end hump to make the following cosmic variance investigation as accurate as possible.This results in a logL Lyα coverage range of 42.3 − 43.1.
Further, we also fit a power law log 10 φ = A × log 10 (L) + B according to Sobral et al. (2017) to the bright-end LF with logL > 43.1, and obtain best-fit parameters A = −1.16+0.15  −0.22 and B = 46.04+9.73 −6.31 .We show this power-law fitting using the red dashed line in figure 5.

Comparison with previous studies
We show the comparison between our Lyα LFs and that of previous studies in figure 5. Results from blankfield spectroscopic surveys (Blanc et al. 2011;Cassata et al. 2011;Ciardullo et al. 2014) are plotted with black lines with different line styles, while the results from narrowband imaging surveys (Hayes et al. 2010;Ciardullo et al. 2012;Konno et al. 2016;Sobral et al. 2017) are plotted with dashed lines with different colors.We also plot Lyα LF from the HETDEX spectroscopic survey by Zhang et al. (2021)    2.0 < z < 3.  2017) carried out deep narrowband imaging with a custom-built NB392 filter for the Isaac Newton Telescope's Wide Field Camera (INT/WFC), they derived Lyα LFs from z ∼ 2.23 LAEs detected in the COSMOS and UDS fields.The Lyα EW 0 criteria of color excess in narrowband in these narrowband imaging surveys are all EW 0 =20 Å, which are consistent with our selections.The best-fit Schechter parameters of the literature mentioned above and that of our studies are summarised in table 3.
In figure 5, our results show a small offset with Konno et al. (2016).Apart from the result of Hayes et al. (2010) which subjects to a huge uncertainty owing to small statistics, all Lyα LFs from the narrowband imaging surveys yield similar of α and L * , but small deviations in φ * .

Measuring the field-to-field variations
In previous studies (e.g.Konno et al. 2016;Ouchi et al. 2008), the field-to-field variation is characterized by σ g , which represents the expected dispersion in galaxy counts owing to the cosmic variance.The expected galaxy counts N can be obtained by integrating LFs over the effective survey volume probed by a field (e.g.Robertson 2010b).Therefore, we make a series of this value N by assigning pointing on each field image, with the coverage smaller than one HSC field of view.For each pointing assigned, a Lyα LF can be derived from a subsample consisting of LAEs within that pointing.By integrating the best-fit Schechter functions, the galaxy abundance N of each pointing can be obtained: Here, logL lim corresponds to the integration limit for Lyα luminosities.We adopt the value from Konno et al. (2016) logL lim = 41.41 erg s −1 for all fittings.Moreover, logL lim has also been discussed in Lyα luminosity densities (e.g.Konno et al. 2016;Ouchi et al. 2008), which characterize the mean Lyα luminosity per unit volume of the pointing.Different choices of logL lim can result in systematic errors in the integral, thus shifting the position of mean galaxy abundances, but leaving the distribution unaffected.
To keep the calculation process of galaxy abundances N consistent for each pointing, we apply the same fitting procedure for each pointing as follows: fixing the faintend slope α at -1.75, as mentioned in section 3.3, and excluding the data points with logL < 42.3, which are beyond the completeness limit, and the data points with logL > 43.1, contaminated by AGNs.The lower limit for integrating LFs is set to be logL lim = 41.41.In addition, since we aim to investigate the field-to-field variation itself, we only include Poisson uncertainties in the error bars when fitting the LFs for this section.
Furthermore, MAMMOTH fields are selected by targeting overdensities.Overdensities could bias the fieldto-field variations which are supposed to measure in random fields.To reduce the bias on the cosmic variance measurements caused by these overdensities, we clip out the regions pre-selected by MAMMOTH on each field image.Among our 8 fields, six (J0210, J0222, J0924, J1419, J1133, and J1349) are selected targeting coherently strong Lyα absorption systems (CoSLAs) (See Liang et al. 2021;Cai et al. 2016), two (J0240 and J0755) are selected targeting at grouping quasars (See Cai et al. in prep.), we refer these as the MAMMOTH targets.The MAMMOTH pre-selected regions are defined to cover these MAMMOTH targets, with an aperture size of 15 comoving Mpc, which is the typical scale of protoclusters at z ∼ 2.2 (see e.g., Chiang et al. 2013).These regions are treated as masks and the LAEs within these regions are excluded when sampling the pointings.
The pointings are assigned with circular shapes and fixed volume coverages.We scale up the pointing areas assigned in the NB387 field images by a factor of depth NB400 /depth NB387 , so that the volumes probed by pointings in NB387 images is the same as those in NB400 images.The centers are created uniformly within the field of view for each field.We calculate the galaxy number density n i for the i-th pointing as follows: where φ i (logL) is the best-fit Schechter function for the i-th pointing.In such a way, the set {n i } can be considered to have the same effective survey volume V eff .There are two reasons we quantify galaxy numbers by integrating the LFs instead of directly using the cumulative numbers.One is to investigate the effect of cosmic variance on LF itself, the other is that this approach can take the completeness correction into account.
Log-normal distributions provide an excellent fit for cosmological density distribution function from CDM non-linear dynamics (Bernardeau & Kofman 1995).The skewness behaviour towards zero can effectively imitate the dark matter fluctuations at the high-mass end (Coles & Jones 1991).Yang et al. (2010) investigated Lyα nebulae in CDFS, CDFN, and COSMOS fields, and reported a strong field-to-field variance.They quantified the field-to-field variance by assuming a log-normal distribution for their Lyα blob counts.Here we adopted a similar model, which assumes the galaxy number density n follows a log-normal distribution n ∼ Log-N(n, σ 2 LN ), where σ 2 LN is the log-normal variance and is related to the actual variance by σ 2 v = exp(σ 2 LN ) − 1.Since the dispersion of N can be estimated as the field-to-field variance plus the Poisson variance (i.e.σ 2 g N 2 + N , see e.g.Robertson (2010b)), we then estimate σ g using the relation: where the last term corresponds to the Poisson variance.n and V eff correspond to the mean values of galaxy number density and effective survey volume respectively.We fit each set of {n i } to this log-normal distribution: where φ(logn) is defined as the number counts of logn within a specified bin divided by the total number, per bin width.We also take the error ranges of n for each pointings into account, and obtain φ(logn) weighted by the reciprocal of these error ranges.C is an arbitrary normalization factor.We fit the three parameters C, logn and σ 2 LN in equation 10 simultaneously, and obtain a best-fit Gaussian mean logn field for each field.
To eliminate the diversity in the Gaussian mean between different fields, we take the best-fit n field from fittings for individual fields, and compute a quantity log( n n ) field for each field.This quantity should have the same mean value of 0 for all fields, we thus fit the following equation to the overall distributions of log( n n ) combining multiple fields: where σ 2 LN and the arbitrary normalizing factor C are treated as free parameters for fitting.The lower and upper limit for the errors of log( n n ) is calculated as log n lower nupper and log nupper n lower respectively.Furthermore, we scale back each φ(log n n ) field according to the normalizing factor C field obtained in equation 10 when combining {log( n n ) i } field for different fields.

Volume dependence of σ g
In the sections above, we described our approach to obtain σ g by randomly assigning pointings with fixed volume coverages for 8 fields, we investigate the effects on σ g by changing the effective volume of pointings in this section.
To obtain the relation between σ g and V eff , we need to keep V eff similar for every pointing during each investigation.In the sections above, we randomly create our pointing centers across the whole circular HSC field of view and clip out the MAMMOTH-specified regions.However, this results in a relatively large scatter in the effective survey volumes V eff , since pointings at the edge lack full coverage of the field image and pointings overlapping with masked and clipped regions have smaller V eff than expected.To deal with this effect, we first generate pointing centers within a confined circle with a radius defined as, the radius of HSC field of view minus the aperture size set for the pointings.In this way, we ensure that every pointings fully cover the field image.Next, we repeat the process described in section 4.1 with various sizes of pointings, and obtain multiple pairs of (V eff , log( n n )), with each pair corresponding to a pointing assigned.We then make bins according to the V eff distribution.For each V eff bin, we calculate the mean value of the effective volumes V eff and σ g by fitting the Gaussian to log( n n ) of the pointings whose V eff lies inside the bin.
We plot σ g against V eff and show the results in figure 6.The range of V eff is limited by the survey depths and effective survey areas of the field images.A too small V eff will cause the pointings to pick up a small number of LAEs or even zero LAEs, which leads to a poor Schechter fit of LFs.On the other hand, a large V eff value that covers a whole field image will create a delta function for its distribution, which is meaningless for σ g investigation.The smallest size assigned for the pointings is 0.2 deg 2 for NB400 and 0.34 deg 2 for NB387 (The depths of the two filters are different and the areas are changed accordingly to keep the volume constant), the corresponding mean number of LAEs per pointing is ∼70.The largest size assigned for the pointings is 0.5 deg 2 for NB400 and 0.84 deg 2 for NB387, the corresponding mean number of LAEs per pointing is ∼170.This results in a volume coverage between 1.86 × 10 5 M pc −3 and 4.61 × 10 5 M pc −3 for the pointings (Although, due to the overlapping between pointings and masked or clipped regions, the actual values of V eff are smaller than expected in general).Detailed information on each V eff bins is shown in table 4. Table 4.The information of each V eff bin.From left to right, the columns are: mean value of V eff inside the bin, number of pointings within the V eff bin, best-fit σ 2 LN of the log-normal fitting, and the resulting σg, respectively.Note that number of pointings in the last bin is relatively smaller than those in other bins, this results in a larger error of the σg calculated, which can also be seen in figure 6.
In previous studies (e.g.Ouchi et al. 2008;Konno et al. 2016;Robertson 2010b), σ g is also defined by the following: where b g is the bias parameter of galaxies, defined as the ratio between the galaxy and matter correlation functions ξ gg /ξ mm (Robertson 2010b), and σ 2 DM (R, z) is the dark matter fluctuation in a sphere with radius R at redshift z.This suggests that σ g only depends on the survey volume at a given redshift at large scales.σ DM can be calculated as following (e.g.Newman & Davis 2002;Moster et al. 2010): where ξ is the two-point correlation function of galaxies, which can be treated as a power-law ξ(r) = (r 0 /r) γ in previous studies (e.g.Somerville et al. 2004;Moster et al. 2010).This makes equation 13 analytically solvable into a closed form: σ 2 DM = C(r 0 /r) γ (Somerville et al. 2004) Since the bias parameter b g only depends on the redshift, the volume-dependency of σ g comes from σ DM only.In this study, we assume a simple power-law σ 2 DM ∝ V β , and thus fit our σ g and V eff,pointing to a simple relation σ g = k × ( V eff 10 5 Mpc 3 ) β with some constant β and k.We find a best-fit value of −1.209 +0.106 −0.106 for β and 0.986 +0.108 −0.100 for k.This is shown as the black dashed line in figure 6.
Using this relation, we find a σ g of 0.063 +0.026 −0.018 for NB387 with survey volume 9.72 × 10 5 Mpc 3 , and 0.034 +0.016 −0.006 for NB400 with survey volume 16.34 × 10 5 Mpc 3 .We use these two values as the field-to-field variance when reporting LF results in section 3 (See section 3.3).

Comparison with predictions from simulation
In previous studies investigating Lyα LFs (e.g.Ouchi et al. 2008;Konno et al. 2016;Robertson 2010b), cosmic variance is estimated from predictions of the ΛCDM model and N-body simulations.Using this approach, Moster et al. ( 2010) have derived a model to calculate σ g for five surveys (UDF, GOODS, GEMS, AEGIS, COS-MOS), as a function of mean redshift z, redshift bin size ∆z, and the stellar mass of the galaxy population m * .We apply their model with z = 2.2, ∆z = 0.08 (consistent with our data), and 8.5 < log(m * /M ) < 9.0 (typical stellar mass for LAEs, see e.g.Ouchi et al. 2020), and obtain the theoretical values of σ g predicted for 5 surveys at z = 2.16 − 2.24.We over-plot these σ g against their corresponding survey volume as blue squares in figure 7. Since the corresponding survey volume of the COSMOS survey is ∼ 18.8×10 5 M pc 3 , which exceeds the range of our measurement, we exclude this data point.We also over-plot the cosmic variance adopted by Konno et al. (2016) as a blue square in figure 7, which is also estimated using the same method.Measurements on σ g taken by our study and the fitting model described in section 4.2 are plotted as red dots and lines in figure 7. It is shown that our measurements of σ g are significantly larger than the theoretical values (about 2-4 times).This suggests that the cosmic variance of LAEs might be different from that of gen-eral star-forming galaxies.This conclusion is consistent with Ito et al. (2021)'s results, which showed that the bias factor of LAEs might be strongly influenced by the HI distribution, and thus might be different from that of general galaxies.
Our results also suggest that the cosmic variances used by previous Lyα LF studies are likely underestimated.A possible cause is that previous approaches use correlation functions (e.g.Robertson 2010a) to estimate the galaxy bias, which only depends on the coordinates of detected galaxies, and then uses the galaxy bias and σ DM to calculate σ g .On the other hand, our measurements are carried out by integrating the LFs which take Lyα luminosity and completeness corrections into account.As a result, introducing Lyα luminosity and completeness correction may cause a higher cosmic variance that is likely more realistic for LF studies.
The underestimated cosmic variance can also explain the inconsistencies on LF fittings between different narrowband imaging surveys (See e.g.Konno et al. 2016;Sobral et al. 2017).For instance, Konno et al. (2016) adopt a value of 0.099±0.017for σ g using predictions by simulation and ΛCDM, while our model gives a value of 0.335 +0.074 −0.061 for σ g under the same survey volume, which is significantly larger.This underestimation can further enlarge the error range of Konno et al. (2016)'s LFs, and thus provide another reason for the offset between Konno et al. (2016)'s LFs and ours presented in figure 5. Similar arguments can be applied to other LF studies, and thus explain the inconsistencies in the LF fittings by different surveys.Although, given that the cosmic variances are only added to the error bars during LF calculation, this underestimation will likely leave the major conclusions unchanged, only enlarging their error ranges.

SUMMARY
We have observed 8 overdense fields targeting the MAMMOTH candidates, with a total survey area of ∼ 11.63 deg 2 .We have carried out a narrowband imaging survey using HSC with two narrowband filters NB387 and NB400 and investigated the Lyα LFs of LAEs selected using the narrowband color excess.The results are summarised below.
1. We fit a Schechter function to our LFs by having the faint-end slope α fixed at -1.75, which is adopted from Konno et al. (2016).2. Using the Lyα LFs of 200 LAE subsamples per field within the pointings created, we have investigated the field-to-field variation that arises from the cosmic variance.After clipping out the MAM-MOTH pre-selected regions to reduce the bias, we create circular pointings on the field images with a fixed volume.We then calculate the LAE number densities for these pointings, by integrating the Lyα LFs with an integration limit of logL lim = 41.41 erg s −1 .We fit a log-normal function to the resulting distribution and compute σ g .
tronomy, the University of Hawaii, the Pan-STARRS Project Office, the Max-Planck Society and its participating institutes, the Max Planck Institute for Astronomy, Heidelberg and the Max Planck Institute for Extraterrestrial Physics, Garching, The Johns Hopkins University, Durham University, the University of Edinburgh, Queen's University Belfast, the Harvard-Smithsonian Center for Astrophysics, the Las Cumbres Observatory Global Telescope Network Incorporated, the National Central University of Taiwan, the Space Telescope Science Institute, the National Aeronautics and Space Administration under Grant No.

Figure 1 .
Figure 1.Transmission curve of the NB387 (blue), NB400 (yellow) and HSC-g (green) band.The curves represent the total transmittance accounting in CCD quantum efficiency, dewar window, the Primary Focus Unit, and the reflectivity of the Prime Mirror.

Figure 2 .
Figure 2. g -NB vs. NB diagram of APER-MAG for LAE selections in each field.The color bar shows the the difference between the g -NB values calculated from AUTO-MAG and that from APER-MAG.The first three selection criteria for the APER-MAGs are plotted as black dashed lines.For J0210, two versions of data reductions yield two different image depths, thus both criteria are overplotted for clarification.Furthermore, the g-band magnitudes are replaced as the corresponding 2σ limiting magnitude, for any object with a g-band fainter than the respective 2σ limit.The g -NB values for these objects are shown as a lower limit in the diagrams.

Figure 3 .
Figure 3. Detection completeness, f det , of NB387 (top) and NB400 (bottom) images.The completeness is calculated for each NB387/400 magnitude bin with width ∆m = 0.5 mag.The red pentagon, blue square, yellow triangle, and green diamond represent the detection completeness of J0210, J0222, J0924, and J1419 fields on the top panel, and J0240, J0755, J1133, and J1349 fields in the bottom panel, respectively.Colored dashed lines are the corresponding 5σ limiting magnitudes for each field in both panels.

Figure 4 .
Figure 4. Lyα LFs and error contours for NB387 & NB400 fields and their corresponding overall results, the points are shifted slightly for clarification.Left panels: LF of each field are shown as the grey data points with different formats specified in the legends, and the overall LFs for 4 NB387 fields (blue) and 4 NB400 fields (yellow) and their corresponding Schechter (Schechter 1976) fits are plotted by solid points and lines.The fitting is done by fixing α at a value concluded by Konno et al. (2016), -1.75.During the fitting, we rule out the points with logL 42.3 which is beyond the completeness correction limit, and the points with logL 43.1 which are contaminated by AGNs.These data points are indicated as hollow circles.Right panels: 1σ and 2σ error contour of fittings to L * and φ * for 8 fields, NB387 overall and NB400 overall LFs.

Figure 5 .
Figure 5.Comparison of our Lyα LFs with previous studies of LAEs at z ∼ 2. The red filled circles, red solid line, and red dashed line are the LF data points, the best-fit Schechter function, and the best-fit power law function derived from our overall sample respectively.The black dotted, dashed and solid lines represent the best-fit Schechter functions obtained by spectroscopic surveys of Blanc et al. (2011); Cassata et al. (2011) and Ciardullo et al. (2014), at redshift range 1.9 < z < 3.8, 2.0 < z < 3.2, and 1.90 < z < 2.35 respectively.The gold, cyan, green, and magenta hollow circles and dashed lines are the data points and the best-fit Schechter functions from narrowband imaging surveys by Hayes et al. (2010); Ciardullo et al. (2012); Konno et al. (2016) and Sobral et al. (2017) respectively.We show the lines within the luminosity ranges limited by corresponding studies.The black hollow triangles are LFs computed by Zhang et al. (2021), which confirmed that the bright-end hump is mainly caused by AGN.

Figure 6 .
Figure 6.Log-normal fits of the overall {( n n )i} distribution for 8 fields, with varying mean effective volume per bin V eff .The bottom panel shows the distributions of {log( n n )i}, and the corresponding best-fit Gaussian function (equation 11).Each V eff bin is indicated by different colors shown in the color bar.The top panel shows σg obtained by each investigation plotted against V eff , we fit a simple power-law σg = k×( V eff 10 5 Mpc 3 ) β , and demonstrate it using a black dashed line.The best-fit β and k has a value of −1.209 +0.106−0.106and 0.986+0.108−0.100 respectively.The vertical and horizontal error bars show the 1-σ error ranges of σg adopted from the Gaussian fitting and the coverages of each V eff bin respectively.Only vertical error bars are taken into account when performing the power-law fitting.

Figure 7 .
Figure 7.Comparison between σg measured by this study (red dots and line) and predictions from simulations (blue squares).The red dots are measurements taken by altering the effective volume covered by each pointing.The red solid line is the best-fit power law described in section 4.2, and the red dashed lines correspond to the 1-σ error range of the fitting parameters.The blue squares are calculated using the approach by Moster et al. (2010).

Table 2 .
The best-fit Schechter parameters for the Lyα LFs of individual fields and the overall data of NB387, NB400, and all 8 fields.All fittings are performed within a logL range of 42.3-43.1 erg s −1 and the faint-end slopes are fixed at -1.75.