Dark Galaxy Candidates at Redshift ~3.5 Detected with MUSE

Recent theoretical models suggest that the early phase of galaxy formation could involve an epoch when galaxies are gas-rich but inefficient at forming stars: a"dark galaxy"phase. Here, we report the results of our MUSE (Multi Unit Spectroscopic Explorer) survey for dark galaxies fluorescently illuminated by quasars at $z>3$. Compared to previous studies which are based on deep narrow-band (NB) imaging, our integral field survey provides a nearly uniform sensitivity coverage over a large volume in redshift space around the quasars as well as full spectral information at each location. Thanks to these unique features, we are able to build control samples at large redshift distances from the quasars using the same data taken under the same conditions. By comparing the rest-frame equivalent width (EW$_{0}$) distributions of the Ly$\alpha$ sources detected in proximity to the quasars and in control samples, we detect a clear correlation between the locations of high EW$_{0}$ objects and the quasars. This correlation is not seen in other properties such as Ly$\alpha$ luminosities or volume overdensities, suggesting the possible fluorescent nature of at least some of these objects. Among these, we find 6 sources without continuum counterparts and EW$_{0}$ limits larger than $240\,\mathrm{\AA}$ that are the best candidates for dark galaxies in our survey at $z>3.5$. The volume densities and properties, including inferred gas masses and star formation efficiencies, of these dark galaxy candidates are similar to previously detected candidates at $z\approx2.4$ in NB surveys. Moreover, if the most distant of these are fluorescently illuminated by the quasar, our results also provide a lower limit of $t=60$ Myr on the quasar lifetime.


INTRODUCTION
Despite a great deal of progress in defining the demographics of galaxies at high redshift (z > 3), our knowledge about the fuel for the formation of the first stars, i.e. the cold gas (T 10 4 K) surrounding the galaxies, is still limited. In addition, due to small sample sizes and technical limitations of the current facilities (Fumagalli et al. 2014), both how this gas forms the large−scale structure of the Universe, the Intergalactic Medium (IGM), and how it keeps star formation active over time are unclear processes (Cantalupo et al. 2012, hereafter C12).
It is well recognized that the densest and most filamentary parts of the IGM play a key role in the formation and evolution of galaxies (Meiksin 2009 and references therein). Recent observations have raised our awareness of the nature of the IGM and CGM (Circum Galactic Medium), thanks to both the absorption (Giavalisco et al. 2011;Turner et al. 2014) and emission (e.g., Wisotzki et al. 2016;Borisova et al. 2016b) signatures of hydrogen at several scales and in different environments, from quasars (QSOs) to radio galaxies (e.g., Cantalupo et al. 2014;Swinbank et al. 2015;Cantalupo 2017).
Theoretical models have suggested the existence of a primordial phase − almost optically dark − of galaxy formation in which there were gas−rich and residing in low−mass halos (e. g., Dekel et al. 2009;Krumholz & Dekel 2012;Kuhlen et al. 2013, amongst others) with very low star formation efficiencies (SFEs < 10 −11 yr −1 ). This less efficient star formation phase of the IGM gas at high redshift could be due to the metal−free gas present in the environment at that epoch, or to the H 2 self−regulation effect  or even to a reduced CGM cooling rate (Cantalupo 2010).
Different approaches have been taken to further investigate this dark phase of galaxy formation in the literature. The different methods that have been used in the past to try to detect the "starless" IGM gas, i. e. just before a considerable star formation occurs, are: (i) HI absorption systems along the line−of−sight to bright background sources (QSOs) at high redshift (e.g., Fumagalli et al. 2011;Prochaska et al. 2013a,b;Lee et al. 2014;Johnson et al. 2015, among others) using one−dimensional data. This method does not help to discern between real isolated dark clouds or gas reservoirs within/around galaxies without the additional information on spatial extent that comes from the emission of the neutral gas.
(ii) HI−21 cm direct imaging (e.g., Giovanelli et al. 2005;Gavazzi et al. 2008). This approach is ob-servationally limited to the dark clouds detected in the local Universe because this line is too weak to be detected at high redshift using current ground−based telescopes.
(iii) Fluorescent emission induced by the cosmic ultraviolet background (UVB), as proposed by the pioneering works of Hogan & Weymann (1987) and Gould & Weinberg (1996). This radiation is produced by ionized gas that recombines and emits fluorescent HI Lyα 1 photons (Cantalupo et al. 2005). The main drawback of this method is the intrinsic faintness of the UVB emission that would imply a Lyα surface brightness of SB ∼ 10 −20 erg s −1 cm −2 arcsec −2 (Rauch et al. 2008), which makes the detection with current facilities very challenging.
(iv) QSO−induced fluorescent Lyα emission can locally boost the signal from dense and otherwise dark gas clouds by orders of magnitude (Haiman & Rees 2001;Cantalupo et al. 2005;Kollmeier et al. 2010, C12) acting as a flashlight on its surroundings. Notwithstanding the complex interpretation of the physics behind the Lyα fluorescence (e.g., Fynbo et al. 2003;Francis & Bland-Hawthorn 2004;Cantalupo et al. 2007;Rauch et al. 2008;Hennawi & Prochaska 2013, among others), and thanks to the support of the 3D radiative transfer models, this seems to be the most promising observational approach and forms the basis of the present investigation using MUSE.
Despite the predictions by several numerical simulations and observational efforts with 8−10 meter class telescopes, in most of the studies conducted so far, the proto−galactic phase preceding the first spark of SF has been poorly constrained. The most convincing observational evidence for this dark phase at high redshift are the objects presented in C12. Using the fluorescent emission induced by the QSO UM287 at redshift 2.4 they detected in a 20hr deep narrow−band (NB) image with VLT−FORS 12 dense and compact gas−rich emitters, named "Dark Galaxies" (DG hereafter), with no detected continuum (stellar) counterpart. The rest−frame equivalent widths, EW 0 , > 240Å of these DG cannot be easily explained by normal star forming regions (Salpeter stellar initial mass function, Malhotra & Rhoads 2002;Charlot & Fall 1993). There are several limitations in the methodology employed in C12. For instance, it required a custom−made NB filter centered on the QSO redshift. This implies that (i) the estimation of the QSO redshift must be very precise; (ii) the results have to take into account possible filter losses and (iii) the candidates need to be confirmed with spectroscopic data. Another limitation concerns the comparison of their results with previous works, because their control samples can be affected by the different observational strategies of the "blank-field" surveys in the literature.
Therefore, the challenging question that we would like to consider here regarding the nature of the dark galaxies is: Do DGs exist at higher (z > 2.4) redshifts and what can be learned from their redshift evolution?
In order to answer this question, we will use (1) an alternative approach of searching for the fluorescent Lyα emission produced by bright QSOs at z > 3 as well as (2) the advantages of an Integral Field Unit (IFU) like the MUSE instrument (Bacon et al. 2010) with the final aim at investigating how the IGM gas is converted into stars. The homogeneous data quality and large wavelength range, translating into a large cosmological volume, offered by the MUSE data presents an unparalleled opportunity for this kind of study, including an analysis of this process with bi−dimensional information. With the help of the third dimension, i. e. the wavelength information missing from NB surveys, we have direct spectroscopic confirmation and also the possibility to explore the presence of other emission lines (e. g., [C iv] λλ 1548,1550 and [He ii] λ 1640 ). More importantly, the use of Integral Field Spectroscopy (IFS) provides the ability to build control samples with essentially the same instrumental and observational conditions, as well as data reduction and analysis techniques, with respect to the main dataset. One drawback will be the relatively small MUSE field of view (MUSE FoV 1 × 1 ) with respect to previous NB images by C12 (VLT−FORS FoV ∼ 7 × 7 ) in exploring the fluorescent volume around the QSO. Indeed, based on the C12 work, we expect to find only 1 or 2 DGs per MUSE field around each QSO. For this reason, in this paper we are combining medium-deep MUSE observations (> 9 hours total exposure time per field) obtained on 6 different fields containing bright quasars.
The QSOs photoionize the surrounding gas, boosting the faint Lyα fluorescent glow expected from the cold gas by a factor of 100−1000 (within a distance of about 10 comoving Mpc) with respect to fluorescence due to the UVB only. Uncertainties include the variable luminosities of the QSOs, the uncertain UV continuum (Lusso et al. 2015), the QSO opening angle (Trainor & Steidel 2013) and further complexities related to the resonant nature of the Lyα line.
Here, we present the MUSE detection of 11 high EW 0 (> 240Å) objects within six medium-deep (> 9 hours) fields at z > 3, of which 8 of these intriguing objects are possible DG candidates fluorescently illuminated by the QSOs. In addition, we present the discovery of a (control sample) population of ∼ 200 Lyα emitters (LAEs) detected in the same fields.
The paper is organized as follows. In § 2 we describe the sample providing details of the MUSE observations, data reduction and post processing. In § 3 we present the systematic analysis of both continuum detected and undetected Lyα emitters within the six MUSE fields. Our results are presented in § 4 and we discuss our findings in § 5. The summary and the conclusions are presented in § 6. Finally, we publish the catalog of LAEs in the Appendix.
Throughout the paper we adopt a flat ΛCDM cosmology with Wilkinson Microwave Anisotropy Probe 9 cosmological parameters of Ω Λ = 0.714, Ω M = 0.286 and h = 0.693 (Hinshaw et al. 2013), corresponding to ∼ 7.5 kpc/ at redshift ∼ 3. We use vacuum wavelengths for the spectral analysis and all magnitudes are in units of the AB system (Oke & Gunn 1983).

SAMPLE AND OBSERVATIONS
Our observations were carried out with MUSE, the second generation IFU mounted on the Very Large Telescope (VLT) at the Nasmyth B focus of the Yepun (Unit Telescope 4) in Paranal, Chile. MUSE has uniquely powerful performance: relatively large FoV (in wide−field mode, WFM, 1 × 1 ) combined with the excellent spatial sampling (0.2 ) and spectral resolutions (R from ∼ 1750 to ∼ 3500) over the wide optical wavelength window (from 4650Å to 9300Å) and high throughput (35 % at 7500Å).

Sample
The six medium-deep fields at z > 3 analyzed in this study were observed between September 2014 and April 2016. We use SExtractor (Bertin & Arnouts 1996) and QFitsView 2 on the NB images centered at 7000Å to measure the seeing (mean full width half maximum, FWHM) on the final combined datacubes. We perform a Gaussian fit to the brightest point sources in each frame. From this, we obtain an average seeing across all frames better than 0.85Å. Most of the observations were carried out under clear or photometric conditions. From the quality assessment of the final combined MUSE datacubes, we obtain a mean (over the six fields) 3σ flux continuum limit in a 1 diameter aperture of 28.5 AB mag whereas 30.5 AB mag represents the mean sensitivity value for the Lyα flux detection (see Section 3.4 for details on how these sensitivities were computed). Table 1 summarizes the measured properties for each field. Their short individual descriptions are provided in the next section. The composite pseudo−color images constructed from the MUSE datacube combining the broad V−, R− and I−band images are shown in Figures 1 and 2. We decided to split our sample into two sub-samples by the redshifts of the targeted QSOs in the respective fields, since our observations target six fields with a difference in the QSO redshift of ∆ z ≈ 0.7 (maximum). Such difference can be important in terms of both cosmological surface brightness dimming (Tolman 1930(Tolman , 1934) that scales as (1+z) 4 , as well as in terms of the explored physical volume. Throughout the paper we will use the term "lower redshift sample" to refer to the fields at z < 3.2 and "high redshift" to refer to those at z > 3.7.

Notes on individual Fields
Low redshift sample − Q0422-3837 or Bulb Nebula: This is the lowest redshift field, z = 3.094, within our sample. Differently from other fields, this observation was targeting a known Lyα nebula around a galaxy that is ∼ 19 comoving Mpc (cMpc) from a bright QSO. It had been discovered through NB imaging (Borisova et al. 2016a).
Our MUSE observations revealed a previously unknown Type II AGN at its center, α (J2000) =04 : 22 : 01.5 and δ (J2000) =-38 : 37 : 19 . It was observed during 20 hours with MUSE. The size of the Point Spread Function (PSF) measured on the final datacube at 7000Å and based on different point sources is 0.7 (the best seeing in our sample). This field is present in both the GALEX (Seibert et al. 2012) and Spitzer (Capak et al. 2012) catalogs, but to our knowledge nothing remarkable about this field has been previously published. The name Bulb comes from the appearance of the Lyα nebula around this AGN in the NB survey that will be presented in a forthcoming paper (Cantalupo et al. in prep.). The RGB synthetic image is shown in the left panel of Fig.  1, where the position of the AGN is marked with the red cross. − Q2321+0135 or Hammerhead Nebula: The second field in our low redshift sample is centered on a radio quiet (RQ) 3 QSO at z = 3.199, α (J2000) =23 : 21 : 14.7 and δ (J2000) =+01 : 35 : 54, and it is presented in the right panel of Fig. 1. This QSO was first spectroscopically discovered in Lyα emission by Schmidt et al. (1987) and was also observed in the Sloan Digital Sky Survey (SDSS, York et al. 2000) with a subsequent follow up by the Baryon Oscillation Spectroscopic Survey (BOSS, Pâris et al. 2012). These authors confirmed the detection of the C IV λ 1550 line, which is likewise detected in the MUSE integrated spectra. The PSF measured at 7000Å is 0.76 . Similarly to the Bulb case, a huge Lyα nebula around this QSO was discovered in NB imaging (Borisova 2016). More details on the Hammerhead will be provided in Marino et al. (in prep.).
− Q1317-0507: Q1317-0507 is a RQ QSO at α (J2000) =13 : 20 : 30.0 and δ (J2000) =-05 : 23 : 35, at z = 3.7. Despite the poor photometric data available in the literature, this QSO has good spectral coverage with UVES. The original time exposure was 10 hours but due to a satellite passing by during one observation, we simply rejected one exposure (15 minutes). The RGB image of this field is shown in the top−right panel of Fig. 2 and the PSF measured is 0.74 .
− Q1621-0042: This RQ QSO, α (J2000) =16 : 21 : 16.9 and δ (J2000) =-00 : 42 : 50, with z = 3.7, is part of the SDSS-DR7 quasar catalog by Schneider et al. (2010). Due to the availability of panchromatic photometric observations together with UVES spectra, this is one of the metal rich QSO used to probe the time evolution of the C IV absorbers (Cooksey et al. 2013). The PSF for the 35 combined exposures (i. e. 8.45 hours, we had to exclude one problematic exposure due to its offset shifts) is 0.77 .
− Q2000-330: The highest redshift field and the only radio loud (RL) QSO within our sample is located at α (J2000) =20 : 03 : 24.0 and δ (J2000) = -32 : 51 : 44 with z = 3.783. The high resolution spectrum of this QSO was taken with the HIRES (High Resolution Echelle Spectrometer, Vogt et al. 1994) instrument and it is part of the KODIAQ survey (O'Meara et al. 2015) along with several other investigations mainly focused on characterizing the CGM. It was observed with MUSE during 10 hours with a PA of 30 • . The PSF in the final datacube has a Gaussian FWHM of 0.84 at λ = 7000Å.

Data Reduction and Post-Processing
The reduction of all 65 hr MUSE data was performed using some of the standard recipes from the latest version of the ESO MUSE Data Reduction Software (DRS, pipeline version 1.6, Weilbacher 2015), complemented with the CubExtractor package (CubEx hereafter, version 1.6; Cantalupo, in prep.) developed to optimally improve the flat−fielding correction and the sky subtraction steps for our specific science case. After retrieving the raw data for each night, we first created the master calibration files using the MUSE pipeline, i.e. the master−bias, the master−flat, the twilight and illumination correction, and wavelength calibration files. Using the DRS routine MUSE scibasic we then processed each individual science exposure, both standard stars and QSO fields, applying the master calibration correction with the recommended parameters. For the illumination correction step we always used the lamp flat−field and the twilight frames closest in time to each individual observation. All these instrumental signatures are removed for each IFU (24 in total), and as output this recipe gives the pre−reduced pixel tables for every IFU exposure. Next we use the MUSE scipost routine to create the individual datacubes, by merging the pixel tables from all IFUs of each exposure. During this step, we also performed the flux calibration using the response curve and telluric absorption correction from one spectrophotometric standard observed during the same night. In addition, scipost applies the geometry and astrometry tables available for each run to the science frames and performs a re−sampling (drizzle algorithm that maximizes the pixel fraction used) onto a 3D grid in order to construct the final datacube. Due to the fact that our observing strategy for each field included a 4 × 90 • rotation pattern with small (< 1 ) offsets, the automatic correction for the absolute astrometry obtained with the pipeline is a source of some uncertainty. For this reason, a double check of the pipeline astrometry correction was required and in the case of clear residual offsets we followed a more classic SExtractor approach in order to correct these offsets.
Once the pipeline−level datacubes were registered, we performed the post−processing using the routines CubeFix, CubeAdd2Mask, CubeSharp and CubeCombine within the CubEx package (Cantalupo, in prep.), since we are interested in reaching very faint surface brightness levels. In particular, using CubeFix we were able to remove the typical checker−board pattern that is seen after the standard data reduction with the pipeline. This is achieved because we self−calibrate each individual exposure at the level of the IFU, slice−by−slice and vertical stacks using the sky−continuum and the sky−lines as "flat sources" together with an iterative masking of any possible continuum sources. Thanks to the CubeFix flat−fielding correction, we were able to reduce the residuals to less than 0.1% of the sky level. Afterwards, we visually inspected the white−light (WL) images created from each CubeFixed−datacube. In those cases where the edges of the individual IFU slices were still visible, or if there is a bright satellite trail or even a problematic channel, we performed a manual masking using CubeAdd2Mask.
Then, we performed a local and flux−conserving sky subtraction on the CubeFixed−CubeMasked− datacube using the CubeSharp routine. This empirical correction takes into account the sky line spread function (LSF) shifts and the variation across the MUSE FoV, conserving the flux and minimizing the residuals. Both CubeFix and CubeSharp were performed twice in order to minimize the contamination from possible unmasked sources when the illumination correction was applied.
Lastly, the CubeFixed-CubeSharped datacubes were combined with a 3σ clipping using both mean and median statistics with the CubeCombine routine. In the case of our analysis, we use the mean−combined datacubes. We also refer the reader to Borisova et al. (2016b); Fumagalli et al. (2016Fumagalli et al. ( , 2017; North et al. (2017) for further details and additional applications of these reduction procedures.

ANALYSIS
The goal of our analysis is to detect Lyα emitters within our sample. In this section, we describe our systematic search and classification of the Lyα emission candidates detected in the MUSE datacubes. In order to perform a consistent comparison between the 6 MUSE fields, we emphasize that the same methodology has been applied to all the MUSE fields as detailed below.

PSF and Continuum Subtraction
Ideally, since we are interested in emission line objects around the QSO (in both spectral and spatial directions) and not in QSO Lyα nebula, removing the nuclear contribution of the quasar should not be necessary for the detection of faint and compact targets. Nonetheless, we decided to perform a PSF subtraction to ensure minimum contamination from the QSO PSF in our LAEs detection by making use of the empirical PSF subtraction of the CubePSFSub routine (part of the CubeExtractor package). Using an averaged−sigma−clipping algorithm, CubePSFSub constructs and rescales the QSO PSF using the NB images created for each wavelength layer, giving excellent results on large scales around the QSO (Borisova et al. 2016b). The next step was the subtraction of the brightest foreground continuum sources within our fields that were carefully removed using CubeBKGSub. This routine estimates the continuum voxel−by−voxel 4 on the basis of a median−filtering performed on the spectrum, which is integrated in 50Å bins and smoothed with a median filter radius of 3 pixels. This allows us to avoid any prominent line features and also to reduce the computational time. Some residuals are still visible in the output datacube, but this has a minimal impact on the extraction procedure of our LAEs considering that we are masking all the bright continuum sources detected from the WL image.

Detection and extraction of Lyα emitters
One of the most important advantages of the IFS is that we can explore the same spatial area over a wide spectral range. To exploit the full capabilities of our MUSE data, our strategy to detect Lyα emitters within our sample was to build three different sub−cubes from each datacube with the same spectral width 200Å (or 160 spectral pixels). The on−source datacube is centered on the QSO Lyα wavelength. Two control sample sub−cubes adjacent to the on−source datacube were extracted on the blue and red sides. For practical reasons, they have the same spectral width as the on−source sub−cube. This choice of the spectral width is justified in terms of the maximum volume (10 cMpc, Trainor & Steidel 2013) where the signature of the fluorescent emission can be detected.
In total we extracted 6 on−source datacubes. These are represented with green symbols in Figures 4, 5 and 7. We also extracted a total of 12 control sample datacubes represented with blue and red colors in the same figures. As mentioned above, the difference in redshifts between our fields corresponds to slightly different analyzed volumes along the spectral direction, because of the constant area coverage. These distances span a range from 36 physical Mpc (pMpc) at redshift < 3.2 to 27 pMpc at redshift > 3.7. We blindly implemented three−dimensional source detection on the 18 reduced and post−processed datacubes using CubExtractor with the same threshold parameters.
Aside from the routines described above, the main purpose of the CubExtractor software is the 3D automatic extraction of sources based on a novel approach used in computer science vision to detect connected re-  gions in binary digital images (see Cantalupo in prep.). The algorithm uses subsets of connected components uniquely labeled on a user−defined property basis, i. e. connected−labelling−component (Shapiro & Stockman 2001). Specifically, we first smooth (with a radius of 0.4 ) both the science and variance datacubes only in the spatial directions for each wavelength layer. Then we require that all detected objects fulfill three conditions: a minimum of 40 connected voxels above a signal−to−noise ratio (SNR) threshold of 3.5 (after the re−scaling factor accounting for the propagated variance is applied) along with a SNR measured on the 1D extracted spectrum above 4.5.
Since the extraction process is based on the noise, estimating the noise correctly is a crucial ingredient of our selection criteria. Since the MUSE pipeline variance tends to be an underestimate of this noise (see Sect.3 in Bacon et al. 2015), we use the propagated variance datacube computed by CubEx that takes the noise sources introduced by both the MUSE pipeline and the CubEx post-processing steps into account. The propagated variance is used to calculate the re−scaling factor applied to each wavelength layer, which in the most extreme case is ≈ 1.95. We also carefully mask the brightest and extended continuum sources detected in the WL image of each datacube 5 , as well as possible skyline residuals to minimize possible artificial detections.
As a result of the three−dimensional segmentation map, we obtain a full catalog of all the line emitters automatically detected in each MUSE field for the on−source and control sample datacubes.

Classification of the Lyα emitters
Although we extensively tested our selection criteria, visual inspection is necessary to remove possible spurious detections of LAEs, such as possible contaminants from [O ii] λλ 3726,3729, [O iii] λ5007 and AGN emitters that were able to pass through the previous masking. Therefore, for each object in our catalog, we tabulated both spectral and photometric information. Specifically, we visually checked the extracted 1D spectrum, where, accounting for different redshift solutions, we were able to distinguish pure Lyα and other emitters by identifying the most prominent emission and absorption line features.
Regarding the photometric properties, from the MUSE datacubes we produce (1) the optimally−extracted (OE), (2) the classical pseudo-NB, (3) continuum and (4) WL images centered on each candidate with a typical size of 30 × 30 . The OE images are constructed by combining all voxels along the wavelength direction that are within the corresponding 3D mask of each detected object from the PSF− and continuum−subtracted MUSE datacubes. This image can be interpreted as a pseudo−NB with a spectral width optimized for the SNR of the candidate (see also Appendix A in Borisova et al. 2016b for a detailed comparison of the OE with the pseudo−NB images).
As we will discuss in the next section, the choice of the continuum image is very critical, especially because based on this image we define a line emitter to be continuum (or not) detected. The ideal case would be the availability of Hubble Space Telescope (HST) images, but these are not available for these fields. We can however take advantage of our IFU datacubes and build the broad−band continuum image. Hence, our approach was to create three continuum images by coadding different spectral ranges. For each field, we considered the spectral layers redward of their QSO Lyα emission and the continuum images were created combining 800 (1000Å), 1600 (2000Å) and all (∼ 3000Å) the wavelength layers in the red part of the datacube. Finally, due to the limited coverage on the blue side of the QSO Lyα emission, we conservatively assumed a continuum slope of β=-2 (f λ ∝ λ β in wavelength space, Meurer et al. 1999) to take in account the shape of the continuum. Then we performed a global statistic of these continuum images while masking the sources in each field. Our final selection of the best continuum image was the deepest one of the three. From the tests performed on our data, the 2000Å continuum image turns out to be the deepest, because its width represents the best spectral compromise able to minimize the contribution of the sky−residual layers. For the sake of completeness, we also checked the classical pseudo−NB and white light images. The results of our classification are summarized in Table 2

Estimation of our detection limits
In order to compute the minimum flux for which we would not be able to detect any candidates, we determine our detection limits for both the continuum and Lyα emission line using the standard deviation (std ) of 100 random locations for each field in our sample. The std is calculated on the continuum and pseudo−NB Lyα images where we mask out all the bright sources with special attention to the scattered light and halos of bright foreground stars. We explored successively larger apertures, with radii from 0.2 to 2 (including the PSF radius) and a 3σ clipping algorithm. We also compare these values with the results from pixel−by−pixel statistics, i.e. the theoretical photon count noise variance, to measure the level of systematics resulting from the sky and continuum subtraction.  The continuum flux is computed as the maximum between the fluxes measured in 9 adjacent PSF size apertures, i.e. it will be always positive. c The rest−frame EWs were determined using the PSF− aperture approach, see Eq. 4. d The gas masses are computed using Eq. 8 in C12. e This measurement is relatively high due to the position of this target in the edge of the FoV.
the continuum and 10 −19 erg s −1 cm −2 arcsec −2 for the Lyα emission. Fig. 3 shows the distribution of the Lyα fluxes and luminosities of the selected LAE candidates.

RESULTS
In this section we present our sample of ∼ 200 LAEs detected in our MUSE datacubes in proximity of quasars and in the control regions within a total volume of ∼ 90 physical Mpc 3 . In particular, we will focus on the LAE Lyα luminosities and equivalent width and their distribution in function of distance from the quasars. The overall properties of the sample is presented in Table 5 in Appendix A.

Lyα flux estimations
Given the recent findings of the extended and diffuse nature of the Lyα emission from LAEs (Wisotzki et al. 2016, Leclercq et al. 2017, measuring reliable Lyα fluxes is not a trivial task because it might depend on both the methodology and available data. In our analysis, the Lyα fluxes were accurately computed from the curve−of−growth analysis (C.o.G. following Drake et al. 2016;Wisotzki et al. 2016) performed on the pseudo NB image centered on the QSO Lyα wavelength with a width of 200Å. By collapsing the corresponding spectral channels of the on−source datacube and assuming the CubeEx coordinates for each target, the Lyα C.o.G. was computed using the fluxes extracted from concentric circular annuli of increasing radii (in steps of 0.2 ) up to 4 . This results in a reasonable value for the characterization of compact objects and their possible extended emission. The total Lyα flux of each object was then determined from the integrated value out to the radius where the surface brightness within a 0.2 annulus is equal to or less than zero. Using the C.o.G. approach, we are able to recover LAEs as faint as 10 −19 erg s −1 cm −2 . Fig. 3 shows the distribution of the Lyα fluxes and luminosities for each low redshift field (Bulb in orange and Hammerhead in blue), for the high redshift fields (in green) and for the full sample (in purple). Although there are definitely uncertainties and limitation in our calculations of Lyα fluxes, we stress that we have used exactly the same method for both the main and the control sample.

The distribution of the Lyα equivalent width
The equivalent width (EW) is a quantitative way of describing the strength of spectral features, both in emission and absorption, compared to the continuum emission. Physically, EWs depend on the initial mass function (IMF) and of the gas metallicity from which stars form, as well as being a useful diagnostic to understand what kind of mechanisms are triggering and sustaining the star formation (e.g., Schaerer 2002Schaerer , 2003. Similar to the Lyα fluxes, the Lyα EW estimation is not unique, and it is very sensitive to the methodology used as well as to the data available for the EW measurements.
In general, we compute the EW as the following ratio: where the numerator corresponds to the Lyα flux. Flux Lyα is computed from the C.o.G. analysis and it is in units of erg s −1 cm −2 . The denominator is the continuum flux density measured in the MUSE continuum image (centered at λ ∼ 6000Å) and extrapolated to the wavelength of the line, assuming that the monochromatic fluxes f ν of all objects are flat in the frequency space. The unit in this case is erg s −1 cm −2Å −1 . However, as explained below, we will use different estimates of these fluxes depending on the nature of the analyzed object. The rest−frame EW(Lyα), EW 0 (Lyα), is: The redshift used in the above equation is defined as the flux centroid of the three−dimensional segmentation mask associated with each detected object. Stellar population synthesis models predict that in the case of continuously star−forming galaxies, the EW 0 (Lyα) produced by Population II stars (hereafter PopII stars) cannot be higher than 240Å except in very extreme cases (Charlot & Fall 1993;Schaerer 2002). EW 0 (Lyα) values above this value may in principle be expected for metal−free PopIII stellar systems (Schaerer 2003;Raiter et al. 2010) and/or Dark Galaxies (C12).
In order to compute the EW 0 (Lyα) of our targets, we decide to follow two different approaches depending on the detection (or not) of our LAE in the continuum image. First, in order to establish if our LAE is detected in the continuum, we measure the continuum flux of our target as the maximum value obtained from the measured continuum flux in 9 different and contiguous positions around the central coordinates of the targets within an aperture with radius equal to the PSF size. This method takes into account possible off-  sets between the spatial peak of the Lyα emission and the stellar continuum (note that the PSF values, listed in Table 1, are all larger than the offsets proposed in Shibuya et al. 2014). Second, if the continuum flux of the target within the PSF size aperture, F Cont @ PSF , is higher than 3 times the standard deviation, std, of the continuum image (3σ Cont , i.e. the local noise, see Section 3.4 for a detailed explanation on how we computed this value) the LAE is considered detected in the continuum.
In the case of F Cont @ PSF < 3 σ Cont our LAE is considered continuum undetected. Of the 186 LAEs selected in our sample, 54% were undetected in the continuum.
In the 4th and 5th columns of Table 2 this statistic is provided for each field. In the case of the continuum detected (CD) LAEs, we used the matched−aperture approach as in C12 and the EW 0 (Lyα) is computed as follows: where Flux Lyα (R) is the Lyα flux within the radius R derived from the C.o.G. analysis, σ(R) is the std of the continuum scaled to the same R apertures and Flux Cont (R) is the continuum flux measured in the same aperture as the Lyα flux. We also masked the contribution of the visible bright continuum objects that were contaminating the measurements extracted from the target aperture, as well as possible contamination from fainter foreground objects inside the aperture.
For those LAEs undetected in the continuum image (CU), we used the PSF−aperture approach and EW 0 (Lyα) is obtained via: (4) where the Flux Lyα (R) is derived as in the case of the CD LAEs and here the continuum flux is computed using the one in the PSF aperture plus 1σ. This method proposed by Feldman & Cousins (1998) ensures an upper limit for the continuum estimation, if the flux in the PSF aperture is positive, otherwise the continuum flux is at least 1σ. This upper limit in the continuum will yield a lower limit in the estimation of the EW 0 .
Despite the complexity and the limitations in estimating the EW 0 , we would like to stress here that we are more interested in the relative distribution of the EW 0 values around the QSOs rather than their absolute values. Similarly to any other measured properties of the Lyα emitters in our sample, we have used exactly the same methods to estimate the EW 0 independent of the position of the object relative to the quasar redshift, both in the main and in the control sample.
In Figures 4 and 5, we present the EW 0 (Lyα) distribution as a function of the redshift difference (spectral distance) from the QSO for the low and the high redshift samples, respectively. The vertical yellow shaded area represents the position of the QSO, while the grey lines indicate the masked position of OH skylines. The CD LAEs are plotted with diamond symbols while the arrows symbolize the lower limit EW 0 (Lyα) estimations for the CU LAEs. Green colors represent the LAEs detected in the on−source (QSO) samples, while the blue and the red ones indicate the control samples. The horizontal dashed line at 240Å denotes the EW 0 (Lyα) limit expected for "normal" star−forming galaxies.
In all MUSE high−z fields we found a higher occurrence of objects with EW 0 (Lyα)> 240Å closer to the QSOs rather than in the control samples. In order to quantify the observed overdensity of high EW 0 (Lyα) objects around the QSO, i. e. in the on−source sample with respect to the control samples, we looked at the EW 0 (Lyα) cumulative distribution. In the left hand panel of Figure 6, the green line indicates the EW 0 (Lyα) cumulative distribution of all (CD and CU) LAEs detected around the QSO. The cyan line denotes the de-tections in the control samples. In the right hand panel, we plot the same but for the continuum undetected LAEs. It is clear in both cases that for EW 0 (Lyα) > 240Å the number of LAEs in the on−source samples is higher. We quantified the probability that the on−source and the control samples are drawn from the same parent population using two non−parametric statistical tests; the Anderson−Darling (AD) test, which is more sensitive to the tails of the distribution, and the Kolmogorov−Smirnov (KS) test, which is more sensitive to the center of the distribution. We decided to use both tests due to our moderate sample size and the fact that the difference between the two samples is more prominent for EW 0 (Lyα) > 240Å. The resulting p-values are lower than 0.007 in both tests. Specifically, in the case of CD and CU LAEs (left panel of Figure 6), we obtained p KS = 0.007 and p AD = 0.005, whereas in the case of the CU LAEs alone (right panel of Figure 6) the p-values are 0.001 in both KS and AD tests. Such low p-values allow us to reject the null hypothesis that the two samples belong to the same population, hence the on−source and the control samples are statistically different.
In the left panel of Figure 7, we combined the EW 0 (Lyα) distribution of the high redshift fields in order to highlight that most of the high EW detections are located closer in redshift to the QSO in the on−source sample. This difference is not due to luminosity effects: If we analyze the Lyα luminosities of these LAEs, plotted in the right panel of Figure 7 as empty diamonds, the distribution of our data does not suggest any significant difference in the luminosities of the LAEs in the on−source with respect to the ones in the control samples.
Similarly, this excess of high EW objects is not connected to an apparent enhancement in the number density of LAE in proximity of the quasars with respect to the control fields, as shown in Fig. 8 where we plot the distribution of the LAEs as a function of the distance from the central ionizing source (AGN in the case of the Bulb field and QSOs for the others). With the exception of the Bulb field, which hosts a lower luminosity AGN, we do not find evidence for an overdensity of LAEs around any of the MUSE QSOs, although the statistical sample is small. Our result is in agreement with the recent findings of Uchiyama et al. (2017) using a sample of ∼ 150 QSOs and of Kikuta et al. (2017) using ∼ 300 LAEs in different environments.
We will discuss in Section 5 the implication of these results in light of our search for dark galaxies candidates fluorescently illuminated by the quasars. Figure 9. Dark Galaxy candidates detected in the MUSE z < 3.2 fields. − Left: the MUSE spectrum within a wavelength range highlighting the observed Lyα emission. The spectrum has been smoothed with a 2 pixel gaussian filter. − Middle: The MUSE Lyα pseudo narrow−band image is shown. The position of the candidate is marked by the red circle. The image was smoothed using a 2 pixels gaussian kernel and the Lyα flux is shown in z−scale. − Right: Continuum broad−band image obtained from the MUSE datacube. We applied a gaussian smoothing with a 2 pixels radius. The continuum flux is plotted with a z−scale stretch between ± 5 σ. In each panel North is up and East is left. Plate scale is 0.2 /pix.

High EW 0 sources
As shown in the previous section, 11 of the 200 LAEs in the total volume explored in this study, including the control samples, present a lower limit on their EW 0 (Lyα) larger than 240Å (arrows in Figs. 4 and 5 above the purple horizontal dashed line). We have demonstrated that these high EW 0 objects tend to be more frequent in proximity of the quasars and in our high redshift sample. In particular, 6 of these are detected in our on−source sub−cubes around the 4 high redsfhit quasars, representing about 25% of the total detected LAEs (24) in this volume. This value is significantly larger than the corresponding fraction in the control samples for the high redshift quasars (about 4%) and for the two fields at low redshift.
In total, 8 high EW 0 objects are present in the on−source samples, i.e. within 10 3 km/s from the quasars (AGN in the case of the Bulb). In Figs. 9 and 10, we show the spectra and postage stamps of these 8 high EW 0 objects detected in the low and high redshift samples, respectively. In particular, for each target, the left panel illustrates a zoom−in of the MUSE spectrum around the detected Lyα emission line while the central and right panels specifically show the Lyα pseudo NB and continuum images obtained from the MUSE datacubes. The position of each object is indicated with a red circle. Their Lyα emission appear compact, similarly to their analogues detected at z ≈ 2.4 by C12. Coordinates, derived photometric and spectral properties, as well as EW 0 lower limits are reported in Table  3.
The Lyα line profiles of these sources is typically asymmetric and in two cases, highlighted in Fig. 11, the emission appears double−peaked. Since the shape of the Lyα profile may be sensitive to the gas kinematics, HI geometry and dust content, our plan is to further investigate these two double−peaked high EW 0 sources as well as the ∼ 60 double−peaked LAEs in our total sample with the help of radiative transfer models in a separate paper.
The main properties of the 3 high EW 0 sources in our control samples are summarized in Table 4 and their postage stamps are shown in Fig. 12 (in the Appendix). We note that these objects do not show any other prominent lines in their spectra. When we considered the 3D extension, i. e. spatial and spectral pixels detected above a threshold, we do not find any significant difference be- tween the 8 objects near to the AGN/QSO and these three high EW 0 objects.

DISCUSSION
The most prominent and characteristic feature of quasar fluorescent illumination is a boost in the EW 0 (Lyα) of LAEs, leading to: (i) a higher frequency of objects without continuum counterparts and (ii) EW 0 limits above 240Å with respect to "blank−fields" (e.g., Cantalupo et al. 2005Cantalupo et al. , 2007. Because the measurement of EWs 0 relies on different methodologies in the literature and because of the different observational techniques and instruments, a proper comparison between the EW 0 of LAEs detected in "quasar−fields" and "blank−fields" has been difficult in previous surveys. Thanks to the new MUSE Integral Field Spectrograph, we were able to obtain a homogeneous sample of Lyα emitting sources around 6 AGN/QSO at z > 3.2 and we were able to build control samples using the same data, the same data reduction and analysis techniques. As expected in the case of fluorescent illumination, we detected an overall excess of high EW 0 sources in proximity of the quasars with respect to the control samples (Figs. 4 and 5). We stress again that, despite the uncertainties and limitation on the measurement of absolute values or limits for the EW 0 , we have used exactly the same methods for our estimates for each source independent of its distance from the quasar. The excess of high EW 0 sources is more prominent in the four quasar fields at z ∼ 3.7. The field−to−field variations could be possibly due to the relatively small MUSE FoV and limited volume probed around each individual quasar. However, they could also suggest intrinsic differences in the quasar properties, such as, e.g., opening angle or age.
In any case, as demonstrated in Section 4.2, the EW 0 distribution in the combined sample around the quasars (on−source) is statistically different than the EW 0 in the control samples at a high significance level.
Is there any other mechanism intrinsic to the sources that would enhance the EW 0 (Lyα) in proximity of quasars without the need for fluorescent "illumination"? High values of EW 0 , if intrinsic, may be due to younger stellar population, different IMF or lower metallicities (see, e.g., Charlot & Fall 1993;Malhotra & Rhoads 2002;Schaerer 2002;Krumholz & Dekel 2012;Orsi et al. 2012). In order for these processes to produce an excess of high EW 0 sources in proximity of the quasar, a relation between the quasar environment and intrinsic galaxy properties would be required. We have explored if the Lyα luminosity and the number density of galaxies are different in proximity of the quasar, possibly indicating a different "environment" but we have found no statistically different results between the on−source and the control samples with respect to these quantities. Moreover, the compact Lyα morphology and the isolated nature of our high EW 0 objects do not suggest any possible effects due to merger activities, although our spatial resolution and the lack of HST imaging would not allow us to detect interactions below scales of a few kpc. While we cannot categorically rule out such a possibility, we see no reason to favour it.
In contrast, the high luminosities of our quasars, the demonstrated existence of "quasar proximity effect" in absorption (at least along our line-of-sight Carswell et al. 1982;Dall'Aglio et al. 2008;Calverley et al. 2011), and the detection of bright Lyα nebulae around these quasars (Borisova et al. 2016b, Marino et al., in prep.) showing that quasars are illuminating their surroundings, all suggest that quasar fluorescence is the most likely explanation for the excess of the compact high EW 0 sources correlated with the quasar redshift in our survey. In this case, the 8 high EW 0 sources without detectable continuum counterparts and EW 0 limits larger than 240Å are the best candidates for Dark Galaxies fluorescently illuminated by the quasars in our survey. The number densities, luminosities and morphologies of these sources are very similar to their 12 analogues detected by C12 at z ≈ 2.4 using NB imaging around a single bright quasars.
How many of these sources have intrinsically high EW 0 without the need of fluorescent "illumination"? Let us consider the fraction of high EW 0 sources in our on−source and control sample at different redshifts. The combined high redshift sample has 25% high EW object on−source and only about 5% in the control sample suggesting that about 1 to 2 of the 6 high redshift LAE with EW 0 limits above 240Å could be objects with intrinsically high EW 0 . Our fraction of 5% of high EW 0 away from quasars at z ∼ 3.6 is consistent with other studies, despite the different methodologies to measure the EW 0 . For instance, Hashimoto et al. (2017, submitted) measured a fraction of about 3% of high EW 0 object in the Hubble Ultra Deep Field (HUDF) using deep MUSE Lyα datacubes and the deepest HST continuum measurements available to date. Despite the small number statistics, this suggest that a significant fraction of the 8 sources with EW 0 limits above 240Å and without continuum counterparts in our survey are strong candidates for Dark Galaxies detected at z > 3.
From the luminosities of these sources and following the approach of C12, we can estimate their total gas masses and star formation efficiencies. In particular, using equation 8 in C12, we estimate gas masses spanning a range between M gas ∼ 0.2 and 6 × 10 9 M similarly to the DG candidates in C12.
To estimate their limit on the star formation rate (SFR), we use the limit on their continuum magnitude (28.8 AB mag extracted from a 1 diameter aperture from the stacked continuum image, see Fig. 13, re−centered at the position of the DG candidates) and convert this value into a SFR using Otí-Floranes & Mas-Hesse (2010) and assuming: i) a Salpeter initial mass function (IMF), ii) a color excess E(B−V)=0 and, iii) an extended burst of 250 Myr. The constraint achieved for the SFR is 0.02 M yr −1 , which yields a star formation efficiency SFE (=SFR/M gas ) of 2.13 × 10 −11 yr −1 indicating that, similar to their analogues at z ≈ 2.4 (C12), our DG are very inefficient at forming stars.
Finally, the distribution of boosted fluorescent Lyα emitters can be also used to constraint the QSO lifetimes (e.g., Cantalupo et al. 2007;Trainor & Steidel 2012;Borisova et al. 2016b). Assuming that our DG candidates are fluorescently illuminated by the QSO, we used a simple geometrical model presented in Borisova et al. 2016a in order to constrain how long the QSO was shining on these proto−clouds of neutral gas, i. e. the quasar life time t Q . Considering the most distant DG candidate within our sample and taking the mean error in the systemic redshift into account, we obtain a distance of 8.7 physical Mpc that corresponds to t Q ∼ 60 Myr. Our estimate is compatible with the results obtained for different QSOs at redshift ∼ 3 analyzed in previous studies (Borisova et al. 2016a).

SUMMARY AND CONCLUSIONS
Making use of medium−deep (∼ 10 hr) MUSE IFU GTO observations around five bright QSO and one Type-II AGN, we have searched for fluorescently illuminated dark galaxies at z > 3.2 among Lyα emitters in proximity of the quasars. Differently than previous surveys based on narrow-band imaging (e.g., C12) and therefore restricted to a fixed volume, we have been able to build control samples at large distances from the quasars using the same data, the same data reduction and analysis techniques thanks to the new capability of the MUSE instrument.
Within a volume of 90 physical Mpc 3 , including the control sample regions, we have identified ∼ 200 line emitters using the automatic source extraction software CubExtractor (Cantalupo, in prep.) complemented with visual analysis. After inspecting their spectral properties in the large wavelength range provided by MUSE, we have found that 186 of these sources are Lyα emitters (LAEs) between redshifts 3.1 and 4.0. We estimated their EW 0 (Lyα) in a homogenous way among the main and the control samples using two different approaches depending whether the sources are detected or not in the continuum. Among all LAEs, we found 11 objects with EW 0 (Lyα) lower limits larger than 240 A -the theoretical limit for galaxies with PopII stellar population (Malhotra & Rhoads 2002;Charlot & Fall 1993). The analysis of the EW 0 (Lyα) distribution reveals that these high EW 0 LAEs tend to preferentially reside within ∼ 10 4 km/s from the quasar systemic redshift. In particular, 6 of the 8 LAEs with EW 0 (Lyα)> 240Å in our high redshift sample lie in close proximity of the quasars. These sources represent about 25% of all LAEs detected within a velocity distance of ∼ 10 4 km/s from the high redshift quasars in our sample. This fraction is significantly higher than the corresponding value in the control samples (4%).
This excess of high EW 0 sources correlated with distance from the quasar is totally consistent with the expectations from quasar fluorescent illumination (e.g., Cantalupo et al. 2007;Trainor & Steidel 2012;Borisova et al. 2016b, C12). Alternative scenarios would require a tight link between distance from the quasars and intrinsic galaxy properties. However, the lack of any correlation between number density of detected LAEs and their luminosities with the distance from the quasars (consistently with other surveys) do not offer currently any support to these alternative scenarios.
In the fluorescent case, the 8 LAEs with EW 0 (Lyα)> 240Å and without continuum counterpart located in close proximity of the quasars represent the best candidates so far for Dark Galaxies at z > 3. Their properties, such as their number densities, compact morphology, luminosities, derived gas masses (∼ 10 9 M ) and star formation efficiencies (SFE < 2.13 × 10 −11 yr −1 ) are remarkably similar to their analogues detected at z ≈ 2.4 with NB imaging by C12.
Although our current sample is limited, this study demonstrate the potential of MUSE observations for the robust detection and characterisation of Dark Galaxies candidates fluorescently illuminated by quasars at z > 3. Compared to NB imaging, the main limitation given by the relatively small MUSE FoV is compensated by the large wavelength range (offering the opportunity to build robust control samples), the immediate spectroscopic confirmation, and the lack of filter (and slit) losses. Every quasar field observed with MUSE will therefore offer the potential to discover new Dark Galaxy candidates and provide crucial information on the early and dark phases of galaxy formation. grant 336736-CALENDS. JB acknowledges support by Fundação para a Ciência e a Tecnologia (FCT) through national funds (UID/FIS/04434/2013) and Investigador FCT contract IF/01654/2014/CP1215/CT0003., and by FEDER through COMPETE2020 (POCI-01-0145-FEDER-007672). TC acknowledges support of the ANR FOGHAR (ANR-13-BS05-0010-02), the OCEVU Labex (ANR-11-LABX-0060) and the A*MIDEX project (ANR-11-IDEX-0001-02) funded by the "Investisse-