Mapping the Morphology and Kinematics of a Lyman-alpha-selected Nebula at z=3.15 with MUSE

Recent wide-field integral field spectroscopy has revealed the detailed properties of high redshift Lyman-alpha (Lya) nebulae, most often targeted due to the presence of an active galactic nucleus (AGN). Here, we use VLT/MUSE to resolve the morphology and kinematics of a nebula initially identified due to strong Lya emission at $z\sim3.2$ (LABn06; Nilsson et al. 2006). Our observations reveal a two-lobed Lya nebula, at least $\sim$173 pkpc in diameter, with a light-weighted centroid near a mid-infrared source (within $\approx$17.2 pkpc) that appears to host an obscured AGN. The Lya emission near the AGN is also coincident in velocity with the kinematic center of the nebula, suggesting that the nebula is both morphologically and kinematically centered on the AGN. Compared to AGN-selected Lya nebulae, the surface brightness profile of this nebula follows a typical exponential profile at large radii ($>$25 pkpc), although at small radii, the profile shows an unusual dip at the location of the AGN. The kinematics and asymmetry are similar to, and the CIV and HeII upper limits are consistent with other AGN-powered Lya nebulae. Double-peaked and asymmetric line profiles suggest that Lya resonant scattering may be important in this nebula. These results support the picture of the AGN being responsible for powering a Lya nebula that is oriented roughly in the plane of the sky. Further observations will explore whether the central surface brightness depression is indicative of either an unusual gas or dust distribution or variation in the ionizing output of the AGN over time.


INTRODUCTION
Recent studies have revealed extended diffuse Lyα emission that covers nearly 100% of the sky around high redshift galaxies (Wisotzki et al. 2018;Leclercq et al. 2020). The mechanisms responsible for lighting up the cosmic web via Lyα radiation have been debated for a long time. Processes such as fluorescence powered by Active Galactic Nuclei (AGN; e.g., Cantalupo et al. 2005;Kollmeier et al. 2010;Kimock et al. 2020), resonance scattering of centrally produced Lyα photons (e.g., Verhamme et al. 2006;Steidel et al. 2011), shockheating in powerful galactic winds (e.g., Taniguchi & Shioya 2000;Taniguchi et al. 2001;Mori et al. 2004), and gravitational cooling (e.g., Rosdahl & Blaizot 2012) of collisionally excited HI atoms have all been put forth to explain the observed emission. In some cases multiple powering mechanisms are believed to contribute to the observed emission (e.g. Vanzella et al. 2017;Herenz et al. 2020); in others, one mechanism appears to dominate (e.g. Borisova et al. 2016;Arrigoni Battaia et al. 2018).
Being able to resolve the environments in which Lyα nebulae reside and the kinematics of the emitting gas will provide a better opportunity to pinpoint which powering mechanisms contribute to the overall emission. Narrowband (NB) and broadband imaging (e.g., Fynbo et al. 1999;Steidel et al. 2000;Matsuda et al. 2004;Yang et al. 2009;Prescott et al. 2012a,b) combined with follow-up long slit spectroscopy (e.g., Weidinger et al. 2005;Matsuda et al. 2006;Yang et al. 2009;Zafar et al. 2011;Prescott et al. 2015a) have been useful in constraining the nebula morphologies and in investigating the kinematics along selected dimensions through the emitting gas, but a key difficulty is not being able to fully map the gas dynamics of the entire nebula with respect to the associated galaxy population. Integral field spectrograph (IFS) observations provide a spatially resolved view of this gas (e.g., Borisova et al. 2016;North et al. 2017;Vanzella et al. 2017;Arrigoni Battaia et al. 2018;Herenz et al. 2020) and allow us to search for signatures of either in-flowing or out-flowing material, potentially providing evidence for cold accretion flows or feedback processes, respectively. IFS data can also be used to detect other emission lines in addition to Lyα that can further constrain the dominant powering source.
Using IFS data from the VLT's MUSE (Bacon et al. 2010), we aim to investigate what mechanisms are playing a role in powering a z≈3.2 Lyα nebulae at the center of a long-standing debate (Nilsson et al. 2006, hereafter "LABn06"). The initial VLT/FORS data revealed revealed an isolated nebula that lacked any associated continuum or radio-loud counterparts, and the surface brightness profile of the nebula resembled those seen in the theoretical simulations of gravitationally cooling nebulae at the time. It was thus concluded that the nebular emission of this source was due to gravitational cooling. Subsequently, Prescott et al. (2015b) studied the environment of LABn06 using deep HST CANDELS imaging and 3D-HST grism spectroscopy as well as improved mid-infrared (MIR) photometry and photometric redshifts from Spitzer and Herschel. They identified 6 continuum sources associated with the nebula along with evidence for an obscured AGN in the vicinity, at the location of a nearby MIR source named 'Source 6' (hereafter S6). In this paper, we return to this controversial Lyα nebula with newly obtained MUSE ob-servations in order to investigate the morphology and kinematics of the nebula and what these findings tell us about the powering of the Lyα emission.
In Section 2, we discuss our observations and data reduction, and in Section 3 we describe the 3D moment analysis of our data and the Lyα line profiles across the nebula. In Section 4, we present our results on the morphology, kinematics, and emission line constraints. In Section 5, we discuss the nature of S6 and how this particular system compares with the larger population of Lyα-emitting nebulae. We conclude in Section 6. We assume the standard ΛCDM cosmology, i.e., Ω m = 0.27, Ω Λ = 0.73, h = 0.7, correspondiing to an angular scale of 7.8 physical kiloparsecs/arcsec (pkpc/ ) at z≈3.2.  (UT 2021 Feb 3, 4, 8, 9, 10, 14, & 15). For the 2019 observations, each OB was split into two 1447 second exposures on the nebula, for a total of 4.02 hours of ontarget exposure time. For the 2021 observations, each OB was split into two 1387 second exposures, for a total of 7.71 hours of on-target exposure time. In 2021, the pointing was shifted by 10 to the West relative to that used in 2019 in order to ensure a sufficiently bright guide star for the Slow Guiding System (SGS). The combined 15 OBs resulted in a total of 11.73 hours of on-target exposure time. The pixel scale for MUSE observations is 0.2 arcsec/pixel, and the wavelength dispersion for each integral field unit (IFU) is 1.25Å/pixel. Observations were taken using the nominal wavelength range (4800-9300Å). The resolving power for the instrument ranges from R∼1770 on the blue side to R∼3590 on the red side. Details of these observations, including the seeing for each night, are listed in Table 1.

Data Reduction
Observations of all 15 OBs were reduced using the MUSE ESOREX (Freudling et al. 2013;Weilbacher et al. 2020, v2.6.2) command line reduction pipeline. Each exposure was processed separately using the pipeline's standard calibration procedure. Initially, each exposure was bias-subtracted (muse bias), dark-subtracted (muse dark), flatfieldcorrected (muse f lat), and wavelength-calibrated (muse wavecal). The master twilight calibration was used to apply a 3D illumination correction to the data cube (muse twilight). All of these calibrations were applied to the science and standard star exposures. The recipe muse scibasic was used to apply these basic calibrations to the science exposures. During this stage, the data were also astrometrically aligned to the MUSE WCS. For the astrometric calibration, a geometry table was used to compute the relative locations of each slice within the IFU. We used the geometry calibration file provided by the ESOREX reduction pipeline, as it has been found to show little variation over time (Weilbacher et al. 2020, see Section 3.6).
The recipe muse standard was then used to create a standard star response file. The standard flux calibration was applied to each science exposure, and sky subtraction was performed in the muse scipost recipe. The sky subtraction step requires a Line Spread Function (LSF) as input to create a model of the sky lines. Although the ESOREX pipeline provides a standard version of the LSF profile for the MUSE instrument, we chose to create our own in order to obtain the most accurate pipeline sky subtraction possible for our data. We computed the spatial offsets between all 30 exposures relative to a reference exposure using the recipe muse exp align, and we determined flux scaling factors for each exposure with respect to our best observation in order to remove the effects of variable observing conditions. We then combined the individual exposures into a final data cube, applying the spatial offsets and flux scaling factors, using the recipe muse exp combine. The exposures were weighted by their integration times for this combination resulting in approximately equal weights. We followed this initial reduction with a single round of self-calibration, which is necessary to remove any low-level instrumental signatures that may not have been removed during the standard calibration steps (Bacon et al. 2017). To do this, we used a white light image of the final MUSE datacube to create a 2D source mask, where only the brightest sources in the field-ofview (FOV) were masked. This 2D mask was then used as an input to perform self-calibration on each individual exposure, and the self-calibrated exposures were then recombined into a final datacube.
The MUSE reduction pipeline sky subtraction step (performed in muse scipost) in version 2.6.2 of the pipeline tends to leave sky residuals in the final datacube. Therefore, we used the ESO recommended postprocessing tool Zurich Atmosphere Purge (ZAP; Soto et al. 2016, v.2.0) to remove any pipeline residuals. We found that optimal sky subtraction was achieved using the ZAP process function with a median level sky subtraction method (zlevel='median') and a continuum filter method for the singular value decomposition (SVD) computation (cftype = 'median'; cfwidthSVD=200). To confirm that we were not over/under-subtracting the sky emission, we used an emission line source at the edge of our FOV for which the [OIII]4959/5007 doublet at 7965.25Å happens to coincide with an optical skyline 1 .
We selected the number of eigenspectra that preserved the expected doublet line ratio for the [OIII] doublet (f 5007 /f 4959 =3; Storey & Zeippen 2000). We note, however, that for the Lyα-emitting nebula at z≈3.2 being investigated here, the Lyα line is on the blue side of the spectrum where concerns about over-subtracted sky lines are minimal.

Flux Calibration Validation
As a check on the flux calibration of our data, we compared our MUSE observations to GOODS-S/CANDELS (Giavalisco et al. 2004;Koekemoer et al. 2011) F775W imaging mosaics produced by the 3D HST Team (Brammer et al. 2012;Skelton et al. 2014;Momcheva et al. 2016). We applied the HST F775W filter transmission curve (including instrumental response) to the MUSE data cube to create a pseudo F775W broadband image. To determine the point-spread-function (PSF) of the MUSE data, we fit two 1D Gaussian models in x and y to the brightest source in the FOV of the filtered MUSE 2D image (a bright, unresolved continuum source located at 03:32:12.644 -27:43:30.11) and measured the full-width-half-maximum (FWHM) of the MUSE PSF to be x=0.89 and y=0.89 . These measurements are consistent with the night-to-night seeing in Table 1, which makes sense given that the source is compact. We convolved the HST mosaic with a PSF kernel matched to the MUSE PSF, and we convolved the MUSE image with a PSF kernel matched to the HST PSF, assuming the HST F775W PSF (0.105 ) derived by Skelton et al. (2014). Finally, we extracted aperture photometry of the source within a 3 radial aperture in both the HST and the MUSE data. The ratio between them provided a flux scaling of f HST /f M U SE =0.98 ± 0.04. We also performed this measurement for the second brightest source in the FOV (an emission line source located at 03:32:13.198 -27:42:40.57) resulting in a flux scale of f HST /f M U SE =1.0 ± 0.07. The weighted average of these two resulting ratios was f HST /f M U SE =0.99 ± 0.03, indicating accurate flux calibration for our final datacube.

Variance Cube Correction
A known issue in the ESOREX MUSE pipeline reduction is that it underestimates the variance cube when resampling the data values from PIXTABLE format to datacube format. This process introduces correlated errors that are not accounted for by the pipeline. Bacon et al. (2017) performed an experiment where they created a test PIXTABLE filled with perfect Gaussian noise (centered on 0.0 with a variance of 1.0), and then resampled this into a final datacube. This resulted in a pixel-to-pixel standard deviation of 0.6, implying the pipeline variance data needed to be multiplied by a correction factor of ( 1 0.6 ) 2 in order to recover the true input variance. This correction factor was obtained using a pixfrac drizzle parameter of 0.8, the default value in the ESOREX MUSE pipeline and the value used for the reduction of our data. Therefore, we applied this correction factor to the variance extension of our final sky-subtracted datacube.

Creating an Emission Line Cube
To detect emission line sources in the datacube, we used the Line Source Detection and Cataloguing tool (LSDcat; Herenz & Wisotzki 2017). The tool performs continuum source subtraction and then smooths the data both spatially (1.19") and spectrally (300 km/s). The smoothing kernel widths used for this step were chosen to match the maximum seeing during the observations and the spectral extent of the target Lyα emission line, derived using visual inspection of the datacube. The maximum extent of the final datacube is 81.58 × 61.19 ×4601.25Å; the portion with full exposure time depth spans dimensions of 60.76 × 61.19 ×4400.0Å.

Masking
In order to construct a 3D mask for the analysis of this nebula, we followed the prescription detailed in Arrigoni Battaia et al. (2019). Specifically, for each wavelength slice in a signal-to-noise ratio (SNR) cube, we identified voxels with SNR≥2. The wavelength slice with the largest connected region of SNR≥2 voxels was 5048.48Å, which we therefore used to define the systemic velocity of the Lyα emission (z Lyα = 3.15). After identifying the wavelength slice containing the largest connected area, we stepped through the cube in the direction of both increasing and decreasing wavelength, attaching new segments to the 3D mask only if there was at least one SNR≥2 voxel that overlapped spatially with the 3D mask defined in previous layers. In doing this, we obtained a final 3D mask composed of 30,983 voxels, extending across 31.25Å, and centered on 5048.48Å. The full extent of the connected SNR≥2 region, measured as the maximum projected edge-to-edge distance of the collapsed mask, is 172.6 pkpc.
We display a voxel density map in Figure 1 panel (a), showing the number of voxels along a given spaxel with SNR≥2. In Figure 1 panel (b), we show the maximum SNR along each spaxel. The locations of the 6 associated continuum sources and S6 are indicated (Prescott et al. 2015b).

MAPPING THE Lyα EMISSION
3.1. Adaptive Narrowband Lyα Image To visualize the Lyα nebula, we created an adaptive narrowband image following Herenz et al. (2020). Voxels with SNR≥2 are summed within a given spaxel, and spaxels that contain no voxels above this SNR threshold are summed along 5Å centered on the central wavelength (5048.48Å). The resulting adaptive narrowband image is shown in Figure 2. We note that this optimallyextracted narrowband image enhances the SNR of the object in a given spaxel because the filter width over which voxels are summed varies from spaxel to spaxel. This differs from standard narrowband imaging in which the filter width is fixed for all spatial positions across a source (see Borisova et al. 2016, Appendix A, for a comparison). We use this adaptive narrowband map to calculate the light-weighted centroid of the Lyα emission and indicate this on all maps presented in this paper. We note that the light-weighted centroid is unchanged if we instead use the MUSE datacube to create a standard narrowband image of the Lyα emission with a fixed width of 31.25Å.

Moment Analysis and Line Profile Shape
We calculated the 0th, 1st, and 2nd moments of our data in order to investigate the surface brightness profile and kinematics of the emitting nebula. Only voxels within the 3D mask with SNR≥2 are included in these moment calculations. The 0th moment is shown in Figure 3 panel (a); this is identical to Figure 2, but contains only spaxels within the 3D mask. Particular regions of interest, selected as the brightest knots of Lyα emission Figure 2. The adaptive narrowband map. Spaxels with SNR≥2 are summed over all wavelengths contained in the 3D mask; spaxels with SNR less than the cutoff are summed over 5Å centered on the central wavelength (5048.48Å). We include the positions of the 6 associated galaxies and S6 identified by Prescott et al. (2015b), along with the lightweighted centroid of the Lyα emission calculated from the adaptive narrowband map (open circle). The dashed circle represents the FWHM (1.72 ) of the 4.6 µm Spitzer/IRAC data used to identify S6.
across the nebula, are indicated by markers A-D, along with the location of S6. In Figure 3 panel (b) we show the Lyα line profiles for each region of interest, mea-sured within an aperture diameter of 12 native MUSE pixels (2.4 arcsec).
The 1st and 2nd moments, corresponding to the velocity centroid and standard deviation, are shown in Figure 4, panels (a) and (b), respectively. Many previous analyses of Lyα-emitting halos have chosen to display a FWHM map as a measure of the line width. This is done by assuming a Gaussian emission line profile, measuring the standard deviation from line center, and multiplying this by ≈2.355, the conversion factor between standard deviation and FWHM for a Gaussian. As seen in Figure 3, the assumption of a Gaussian profile is not necessarily valid for Lyα emission lines; the profiles often show multiple peaks or enhanced tails, likely due to the resonant nature of Lyα. Thus, we instead chose to show the standard deviation directly as an "apparent velocity dispersion" map, without converting to FWHM or making any assumptions about line shape (Herenz et al. 2020).

Morphology and Size
From the maps presented in Figure 1 and 2, we can see multiple regions of significant Lyα emission, but the brightest parts of the nebula are not coincident with the associated galaxies in the region. Instead, LABn06 consists of a two-lobed structure centered on the position of S6, which lies within 2.2 (17.2 pkpc) of the lightweighted centroid of the nebula. Assuming the source powering the emission lies near the light-weighted centroid of the nebula, the two lobes could be tracing out a bipolar structure with a wide opening angle.
We quantify the diameter of the two-lobe structure using the distance from region D to region A (as labeled in Figure 3) -a distance of 7.8 (61 pkpc). Some low signal-to-noise emission appears to extend over to the West of the light-weighted centroid, connecting the brightest Northern peaks (A and B) to the three associated galaxies near region C. The vertical dashed line in Figure 3 marks a dip in the Lyα surface brightness distribution. We will refer to the material west of this vertical dashed line as the "western lobe" throughout the remainder of the paper. The maximum source size down to the limit of our data, as defined in Section 2.6, is 172.6 pkpc.

Surface Brightness Profile
It is interesting to investigate how the surface brightness profile of LABn06 compares to other Lyα-emitting sources in the literature. In Figure 5 panel (a), we show the circularly-averaged, radial surface brightness profile for the nebula centered on the location of S6. For com-parison, we also show the measured profiles for a stacked sample of Lyman Break Galaxies ), a sample of Lyα halos surrounding 2 radio-loud and 16 radio-quiet Type I quasars (Borisova et al. 2016) and four halos around Type II AGN (den Brok et al. 2020). All profiles displayed are redshift-dimming-corrected to z=3 by the factor ( 1+z 4 ) 4 . The most striking finding is that the inner region of LABn06 appears distinct from all of the other profiles, showing a dip in the surface brightness at the location of S6. Outside of 25 pkpc, the profile resembles that of other nebulae. To make a quantitative comparison, we fit both a power law and an exponential profile to LABn06, only including surface brightness values at radial distances greater than 25 pkpc (the distance of peak A) from S6. The nebula is best fit by an exponential with a scale radius of 24.2 pkpc (χ 2 =6.5). We compare LABn06 to a subset of sources well fit by exponential profiles in Figure 5 panel (b), which emphasizes the distinct dip in the center of LABn06 compared to other sources.

Asymmetry
From Figures 1, 2, and 3 panel (a), it is obvious that LABn06 is not completely circularly symmetric, instead showing a bipolar structure. In an attempt to quantify the level of asymmetry, we apply both a moment-based method and a Fourier Decomposition method, as detailed in den Brok et al. (2020).
For the moment-based method, we calculate the second-order moment of the spatial pixels with respect to the known location of S6, including all pixels contained by the collapsed 3D mask. We do not include flux-weighting of the pixels in order to emphasize how diffuse emission is affecting the asymmetry. These second-order moments are then used to quantify the circular asymmetry of the nebula via the dimensionless parameter α (den Brok et al. 2020): where M xx , M yy , and M xy are the second-order image moments. A value of α = 1 corresponds to a circularly symmetric nebula while α < 1 corresponds to a less circularly symmetric distribution of emission. For LABn06, we find α = 0.60, indicating a degree of asymmetry comparable to other Type II AGN (α =0.4-0.7) but outside the 25 -75 percentile range for Type I AGN (α =0.65-0.85) (den Brok et al. 2020).
The Fourier Decomposition method quantifies the nebula asymmetry as a function of distance from the known location of S6. Following den Brok et al. (2020), we first reproject the (x,y) pairs of the surface brightness pixels onto a grid of (r,θ) values. We then decompose SB(r, θ) into its Fourier series and solve for the Fourier coefficients a k (r) and b k (r). These coefficients can then be combined into a single coefficient: A circularly symmetric distribution will be dominated by the 0th order coefficients (a 0 , b 0 ). By contrast, an elliptical profile centered on S6 will have a significant contribution from 2nd order coefficients (a 2 , b 2 ), and an off-center elliptical profile will have a significant contribution from the 1st order coefficients, (a 1 , b 1 ). Therefore, the ratio c k (r) c0(r) quantifies how much the profile diverges from a circularly symmetric one at a given radial distance.
In Figure 6, we display ( c1(r) c0(r) ) 2 +( c2(r) c0(r) ) 2 for our nebula along with profiles of Type I and Type II AGN samples (den Brok et al. 2020, private communication). We find that the nebula shows higher asymmetry at very small radii (r < 10 pkpc), although this region is somewhat compromised by the seeing limit of the data. LABn06 is fairly symmetric at r = 10 − 50 pkpc, but begins to deviate from circular symmetry at larger distances from S6 (r > 50 pkpc).

Kinematics
Overall, the velocity field (1st moment) and apparent velocity dispersion (2nd moment) maps in Figure 4 give the impression of relatively quiescent Lyα-emitting gas in LABn06. Most of the velocity values lie in the range V≈[-500,250] km s −1 , where 0 km s −1 corresponds to z Lyα = 3.15. A hint of ordered motion is apparent as a gradient between positive velocity values (V≈250 km s −1 ; pink) to the South-East of S6 and negative velocity values (V≈-500 km s −1 ; blue) to the North-West of S6. The one exception is the steep velocity gradient seen across the western lobe of the nebula, where the velocities range from -450 to 250 km s −1 over a distance of about ∼5 (∼40 pkpc). The western lobe is in the vicinity of three of the associated galaxies, including the brightest one and the only one that is spectroscopically confirmed (Prescott et al. 2015b). This steep velocity gradient could be evidence that the gas in the western lobe is kinematically distinct from the gas centered on S6 and is actually linked with the 3 associated galaxies in the western lobe.
The apparent velocity dispersion map shows a range of line widths from σ v ∼200 km s −1 to σ v ∼500 km s −1 . The regions of highest apparent velocity dispersion seem to be concentrated mostly to the North of S6's MIR position as well as in the western lobe of the nebula. However, this is likely driven by the presence of double peaked line profiles in these regions (Section 3.2).

Upper Limits on NV, CIV, HeII, & CIII]
AGN are a common powering mechanism for Lyα nebulae. To test the AGN-powering hypothesis in   Figure 5. Circularly averaged surface brightness profiles all corrected to z=3. We only show profiles for distances greater than 10 pkpc (1.2 ) as smaller radii are compromised by the seeing. (a) Comparisons between LABn06 (pink) and other Type II AGN observed with MUSE (green, red, yellow, and blue solid lines; den Brok et al. 2020). The profile for LABn06 is best fit by an exponential profile with a scale radius of 24.2 pkpc (dashed pink line). We also plot the median of a sample of Type I radio-quiet and radio-loud quasars (Borisova et al. 2016, black solid line); the 10th and 90th percentiles of the Type I AGN distribution is given by the gray shaded region. We also show the Lyα surface brightness profile derived using a stack of Lyman Break Galaxies (Steidel et al. 2011, dashed black). (b) Surface brightness profiles for the subset of sources well-fit by exponential profiles, normalized to 25 pkpc. Figure 6. Fourier decomposition of the LABn06 surface brightness map as a function of distance from the location of S6 (purple). The dependent variable combines the first three fourier coefficients of the decomposed surface brightness distribution, where the 0th order coefficient corresponds to a circularly symmetric distribution and coefficients 1 and 2 quantify deviations from circular symmetry, as discussed in the text. For comparison, we plot the same profile for a sample of 4 Type II AGN (cyan) and a median profile for a sample of 19 Type I AGN (den Brok et al. 2020). The dashed black lines bounding the gray shaded region indicate the 25th and 75th percentiles around the Type I median. LABn06 shows a hint of asymmetry at small radii (r < 10 pkpc), and a rising asymmetry profile at large radii (r > 50 pkpc), approaching what is seen in the Type II AGN sample.
LABn06, we searched the data for additional rest-frame ultraviolet (UV) emission lines that would be indicative of an AGN (i.e., NVλ1240, CIVλ1550, HeIIλ1640, CIII]λ1909, Feltre et al. 2016). We inspected the datacube in regions contained by the 3D mask at the locations where one would expect to find these lines given a source redshift of z Lyα = 3.15. After aligning the center of the 3D mask created in Section 2.6 to the expected observed wavelength for each line, we summed over the wavelengths included in the 3D mask (a region of 31.25Å across) to create a flattened 2D image. These images are displayed in Figure 7 with several higher SNR peaks indicated by black circles. Peaks outside the lowest surface brightness contour of the nebula all appear to be linked with foreground and background HST sources not associated with this source. Interestingly, within the Lyα nebula the closest source to the peaks in the NVλ1240 and CIII]λ1909 images is S6. We measure the flux of these peaks within a 1.2 (6 pixel) radial aperture centered on the peak and find f N V = 3.05×10 −18 erg s −1 cm −2 (SNR= 3.1) and f CIII] = 3.62×10 −18 erg s −1 cm −2 (SNR= 2.4). We estimated the uncertainty on these measurements by laying down 10000 random apertures across the images and taking the standard deviation of this distribution. We also measured upper limits on all four of these lines following the scanning method of Borisova et al. (2016, Section 4.5) to allow for the possibility of kinematic offsets between Lyα and other emission lines. For each line, we used the previously defined 31.25Å-wide (≈1180-1860 km s −1 ) 3D mask (Section 2.6) to scan a ±3000km s −1 (≈100-160Å) wide region centered on the expected observed wavelength of the line assuming a source redshift of z Lyα = 3.15. Beginning with the lower wavelength edge of the mask located -3000 km/s from the expected line center, we summed all of the voxels contained by the mask. We then shifted the entire mask by one spectral unit (1.25Å), summed over the 3D mask, and repeated this procedure until the upper wavelength edge of the mask was +3000 km/s from the expected line center. Calculating the rms of these summed values provides us with a conservative measurement of the upper limit on the flux of these lines within the Lyα nebula. These results for the full nebula are summarized in Table 2. To make sure that using the entire 31.25Å mask was not overly conservative, we repeated this scanning method with only the central 15Å (≈570-890 km s −1 ) and 7.5Å (≈280-450 km s −1 ) of the mask. In doing this, we found the same SNR peaks as indicated by black circles in Figure 7.
To allow for the possibility that any NV, CIV, HeII, or CIII] emission is confined only to the brightest regions of the Lyα nebula, we also modified the previous method to incorporate an aperture-based approach. In repeating the scanning process described above, we focused on the specific regions of interest (A-D, S6) shown in Figure 3. That is, we summed over the 3D mask but only included voxels within a 1.2 arcsec radial aperture centered on the regions of interest. We then measured the rms of these summed values to derive an upper limit for the NV, CIV, HeII, and CIII] flux within each region A-D as well as the region centered on S6. These results are summarized in Table 2.
In Figure 8, we show how the line ratio upper limits from the extended nebula in LABn06 compare to other Lyα nebulae in the literature on a HeII/Lyα vs. CIV/Lyα diagnostic diagram. These results are consistent with the line ratio measurements and upper limits seen in other Lyα nebulae, most of which have been identified around AGN. We note that many of the other supposedly AGN powered LAN targeted with MUSE (e.g. Borisova et al. 2016;Arrigoni Battaia et al. 2019) do not show these expected AGN emission lines yet their Lyα halos are comparable in size and even more luminous than LABn06. While not conclusive, the UV emission line ratios found for LABn06 are consistent with the possibility of AGN-powering in this system.
In Figure 9, we also show the spectrum extracted from 1.2 and 2.4 diameter regions of the nebula centered on the location of S6, as indicated by its MIR detection. The Lyα emission line along with other UV metal emission lines typically indicative of AGN powering are indicated by dashed black lines and cyan labels. No other emission lines besides Lyα show a significant detection in the MUSE data. While it would be beneficial to compare the UV metal emission upper limits obtained for S6 to those seen in typically faint or obscured AGN, the samples of AGN at z∼3 for which measurements of the UV metal emission exist are 2-3 orders of magnitude brighter at 22µm than S6. Further, these samples were selected to have bright CIVλ1550 and Lyαλ1216 (Alexandroff et al. 2013) or as some of the most optically bright AGN at z∼3 (Borisova et al. 2016). This suggests the need for further study into the UV emission line properties of the obscured AGN population at fainter infrared magnitudes.

DISCUSSION
Our MUSE coverage of LABn06 represents one of the deepest datacubes targeting a large Lyα nebula. In what follows, we discuss the nature of S6, the evidence for AGN-powering and resonant scattering in this system, and how the morphology and kinematics of LABn06 compare to other AGN-powered nebulae. We end by Figure 7. Flattened 2D images centered on the expected wavelength for each UV metal line assuming a source redshift of zLyα = 3.15. Contours of the Lyα emission are plotted at arbitrary levels. We also indicate S6 (purple star) and the locations of some of the largest peaks in these images (black circles; 1.2 radius apertures). S6 seems to lie close to a peak in the NVλ1240 image and one of the peaks in the CIII]λ1909 image.   revisiting the earlier arguments made for gravitational cooling in this system in light of our new data.

The Nature of S6
The nature of the MIR source (S6) is still not obvious. While previous investigation of the MIR SED for this source found it was best-fit by the Mrk 231 Type 2 AGN template, the FIR 70-500 µm fluxes were reminiscent of cold dust heated by star formation (Prescott et al. 2015b). At the same time, Prescott et al. (2015b) found that the MIR fluxes of S6 placed it solidly within the AGN color selection boxes of Donley et al. (2012). However, it is true that MIR AGN color-selection can be contaminated by higher redshift (z>2) star-forming galaxies (Barmby et al. 2006;Park et al. 2010;Donley et al. 2008). In these cases, power-law selection is an efficient way of cleanly isolating obscured AGN. While typical star-forming galaxies tend to display a dip in emission in the MIR, AGN-heated dust is expected to produce a monotonically rising SED across the IRAC bands (Rees et al. 1969;Neugebauer et al. 1979;Assef et al. 2010). For example, while a few contaminating high redshift starburst galaxies fell into the AGN selection box of Donley et al. (2012) due to their measured MIR colors, the AGN nature of these sources could be ruled out due to the non-power-law shape of their IRAC SED. Taking a similar approach, inspection of the IRAC fluxes of S6 shows that its SED rises monotonically between 3.6 and 8µm. Following the procedure in Donley et al. (2008, Section 2.2), we find that its IRAC SED is well-fit by a power law, f ν ∝ ν α , with α = −1.54 ± 0.17 (χ 2 = 1.83 P χ = 0.40), consistent with obscured AGN power-law selection criteria.
Despite the evidence in the MIR for the AGN nature of S6, no counterpart X-ray detection was found in deep 1Ms, 2Ms, or 4Ms exposures (Prescott et al. 2015b) or even in more recent, 7Ms data (Luo et al. 2017). Assuming S6 does contain an AGN, the lack of X-ray emission prompts us to ask whether the AGN is heavily obscured or whether it is a ramped-down source. Previous work showed that it is relatively common for a heavily obscured AGN to be undetected in X-ray observations with t<2Ms (Donley et al. 2007;LaMassa et al. 2019;Pouliasis et al. 2020), however, it is unclear what fraction of power-law selected sources are expected to be undetected at the limit of the 7Ms data (2.7×10 −17 erg s −1 cm −2 ; Luo et al. 2017). Updated studies on the X-ray properties of power-law selected AGN are needed in order to determine whether the lack of an X-ray detection in the 7Ms data can rule out the presence of an obscured AGN in S6. A rampeddown AGN could also explain the lack of X-ray emission.
In this case, the fact that we still observe AGN-powered MIR emission could indicate that the AGN ramp-down took place relatively recently.
In deeper radio observations that now exist for this field (Miller et al. 2013), we find an unresolved 5σ 1.4GHz detection that is only 0.3 away from the location of S6. We can use this unresolved radio detection to investigate how S6 compares to the well-known FIRradio correlation for star-forming galaxies (FRC; Helou et al. 1988;Yun et al. 2001;Magnelli et al. 2015). Using a FIR SED template of the star-forming galaxy Arp 220 (Polletta et al. 2007), shifted to the observed redshift of z=3.15 and scaled to the observed 250µm flux density of S6, we estimated the FIR (42-122 µm) luminosity of S6 to be 7.69×10 45 erg s −1 . We K-corrected the radio luminosity assuming a power law spectral index of -0.71 for S ν ∝ ν α (Read et al. 2018, Equation 1) and estimated the FRC parameter to be q F IR =2.04±0.08. This q F IR value places S6 within the star-forming-galaxy locus of the FRC, a region also populated by radio-quiet seyfert galaxies (Morić et al. 2010;Sargent et al. 2010;Padovani et al. 2011;Del Moro et al. 2013). This indicates that S6 does not have an excess of radio emission, making it consistent with what is expected for radio-quiet AGN and star forming galaxies.
Based on the results of the MIR AGN-selection techniques (IRAC colors and power law slope), the lack of evidence for excess radio emission, and the evidence for star formation-heated dust in the FIR, we conclude that S6 is likely a radio-quiet, obscured AGN/star-forming galaxy composite in which the AGN is either heavily obscured or recently ramped down. In what follows, we compare the LABn06/S6 system with other AGNpowered nebulae and refer to S6 as the obscured AGN.

Morphologically & Kinematically Centered on the AGN
We have confirmed that this highly debated Lyα nebulae is larger than originally reported (at least ≈61 pkpc in size from regions A to D, with fainter emission extending out to 173 pkpc), and is spatially offset from the brightest optically detectable galaxies in the region. Previously, Prescott et al. (2015b) provided evidence that an obscured AGN was buried ≈30 pkpc from the peak of the emission. In our new MUSE observations, we found that the Lyα emission at the location of the obscured AGN is coincident in velocity with the kinematic centroid of the nebula (Figure 3, panel (b)), and that this obscured AGN lies ≈17 pkpc from the light-weighted centroid of the nebula, placing the AGN roughly at the center of a two-lobed Lyα-emitting structure (Figures 1 and 2). The spatial coincidence may be even closer if the gas in the western lobe is distinct from LABn06 as suggested by the steep velocity gradient in this region. To explore this possibility, we perform the following exercise: we remove the gas in the western lobe from our calculation of the light-weighted centroid of the Lyα emission by masking out all values West of the surface brightness dip (vertical dashed line) in Figure 3 panel (a). In this case, we find that the light-weighted centroid lies even closer, within 0.3 (2.3 pkpc) of the AGN. Taken together, the MUSE data show that the AGN is likely located at both the spatial and kinematic center of the Lyα nebula.
While the spatial proximity of the obscured AGN to the light-weighted centroid of the Lyα emission could be a coincidence, we think this scenario is much less likely. The obscured AGN was first identified using standard Spitzer/IRAC MIR color selection (e.g., Lacy et al. 2004Lacy et al. , 2007Donley et al. 2012). Using the compilation of Mendez et al. (2013, their Figure 2) and extrapolating the Donley et al. (2012) curve to the flux density of the obscured AGN in LABn06 (≈1.57×10 −29 erg s −1 cm −2 Hz −1 ), we estimate that the expected surface density of such obscured AGN is ≈10 3 -10 4 degree −2 dex −1 . In LABn06, the obscured AGN is located within a radius of 2.2 (0.3 ) of the light-weighted centroid of the nebula, corresponding to a circular area of ≈1×10 −6 degrees 2 (≈2×10 −8 degrees 2 ). Therefore, by multiplying this area by the expected surface density, we estimate that there is a 0.1-1% (0.002-0.02%) chance of an obscured AGN (as selected by Donley et al. (2012) color-cuts) showing up by chance within 2.2 (0.3 ) of the light-weighted centroid of LABn06.
Our deep MUSE data provide indirect evidence for AGN-powering in LABn06, with the nebula being morphologically and kinematically centered on an obscured AGN, and with derived upper limits on HeII and CIV that are consistent with AGN-selected systems. In what follows, we explore how LABn06 compares to other AGN-powered nebulae, both in terms of morphology and kinematics.

Morphology compared to other AGN-powered Nebulae
In Borisova et al. (2016), the authors targeted bright quasars using MUSE and found extended Lyα emission around every single source. These presumably AGNpowered Lyα nebulae span a range of sizes (100-300 pkpc) and luminosities (10 43.3 -10 44.6 erg s −1 ). Compared to this sample, LABn06, with a full extent of ≈173 pkpc and a Lyα luminosity of ≈10 43 erg s −1 , seems to be comparable in size, although only half as luminous. Although rigorous comparisons should use size measure-ments down to the same surface brightness threshold (Wisotzki et al. 2018;Leclercq et al. 2020), this qualitative comparison suggests that LABn06 is quite extended for its low surface brightness and likely falls below the correlation implied for other AGN-powered Lyα nebulae at this redshift (Christensen et al. 2006).
To make a more robust comparison of the nebula size, we use the best-fit scale radius of the circularly averaged surface brightness profile. LABn06 is best described by an exponential with a scale radius of 24.2 pkpc. In Figure 5 panel (b), we plot the profiles for LABn06 and eight other sources best-fit by an exponential profile, all normalized to 25 pkpc. LABn06 appears most similar to the Lyman Break Galaxy stack of Steidel et al. (2011), both in terms of scale radius (25.2 pkpc) and in terms of profile shape at distances between 25-60 pkpc. In the outer regions, LABn06 also resembles other AGNpowered nebulae observed with MUSE in terms of profile shape, although it is a bit smaller in terms of scale radius. LABn06's scale radius is at the lower end of the range spanned by the 5/19 sources from Borisova et al. (2016, r scale = 20.69, 29.71, 38.88, 40.12, 55.66 pkpc) that were also best-fit by an exponential profile.
LABn06 shows a similar level of asymmetry compared to other AGN-powered nebulae, particularly those around Type II AGN. The dimensionless asymmetry parameter (α) calculated for LABn06 is roughly at the median α value of the population of Type II AGN investigated by den Brok et al. (2020). The asymmetry profiles derived using Fourier decomposition in Figure 6 revealed that while LABn06 appears relatively symmetric at intermediate radii, it exceeds the typical asymmetry of Type I AGN at both small (<10 pkpc) and large (>50 pkpc) distances from the AGN. This fact combined with the bipolar morphology, and the relatively mild velocities are all consistent with LABn06 being only moderately inclined along our line-of-sight, potentially giving us an "out-of-the-cone" view of the system.
Despite the similarities to other AGN-powered systems, at small radii LABn06 is unique in showing a pronounced hole around the position of the obscured AGN. This could be due to high central HI column density and/or dust obscuration preventing Lyα photons from escaping the nebula in the vicinity of the obscured AGN. Alternatively, it could be due to a lack of HI gas in the central region. Finally, it is possible that this dip reflects a recent ramp down in AGN power roughly 10 5 years ago, in which case, we are seeing an ionization echo in the outskirts of LABn06.
More speculatively, if the AGN in LABn06 has recently ramped down, the fact that one lobe of the nebula is brighter than the other might suggest that we could be viewing a late-day version of the Lyα nebula in Weidinger et al. (2005), which showed a biconical structure inclined at an angle around a z∼3 quasar. If the central engine shuts off, the light traveling to the observer from the backside cone would be delayed relative to that from the front side cone. This might create enough of a time window that we would see a bicone that is not symmetrically illuminated. If this interpretation is correct for LABn06, it would imply that the northern lobe in this system is further away from us along the line-of-sight than the southern lobe.
While to our knowledge this is the first case of a Lyα nebula with a central hole at the location of the AGN, we note that having an AGN offset from the brightest regions of Lyα emission is not unprecedented. Other examples of offset AGN powering extended Lyα nebulae are known to exist (Kurk et al. 2002;Prescott et al. 2009Prescott et al. , 2012b. In addition, hydrodynamical simulations of Lyα nebulae reveal that the brightest regions can be offset substantially from the location of star forming galaxies if AGN are included as a powering mechanism (Kimock et al. 2020).

Kinematics compared to other AGN-powered Nebulae
AGN-powered nebulae in the literature show a diversity of kinematics. For example, Borisova et al. (2016) found that most of the nebulae showed relatively chaotic velocity fields, without any coherent gradient in velocity. However, the two largest nebulae presented clear velocity shear, in one case along the major axis of the nebula but in the other, surprisingly, along the minor axis. Similarly, Prescott et al. (2015a) found signs of large scale rotation as indicated by Lyα, CIV, HeII, and CIII] in the Lyα nebula PRG1 at z∼1.7. Zafar et al. (2011) identified evidence for ordered motion in the Lyα nebula around the binary quasar Q0151+048 at z∼1.9, and Arrigoni Battaia et al. (2019) identified an obvious velocity gradient centered on two LAE and two AGN in an extended Lyα nebula at z∼3.2. Finally, Herenz et al. (2020) found clear signs of a velocity gradient perpendicular to the principal axis of LAB1 at z∼3. In this context, LABn06 is well within the range of other AGNpowered nebulae, with line kinematics that are somewhat chaotic, but with some signs of coherent velocity shear along the major axis, as seen in Figure 4.
In terms of velocity dispersion, other Lyα nebulae in the literature show a wide range of values. For their sample of z∼3 Lyα-emitting nebulae, Borisova et al. (2016) found velocity dispersions of σ v =150-450 km s −1 , although the three largest sources and the two radioloud sources displayed a peak in velocity dispersion near the location of their AGN. Spectroscopically derived line profiles for these nebulae show a range of both narrow and broad lines with single, asymmetric shapes or, in a few sources, hints of double peaks. Similarly in LAB1, Herenz et al. (2020) found a range of velocity dispersions (σ v =200-550 km s −1 ), although the highest dispersion was found in the region occupied by identified sources. The line profile derived from this region of highest dispersion appeared to be double peaked, although another double peaked profile was found in a region devoid of sources.
In LABn06, we see varying Lyα line profiles across the nebula ( Figure 3) and apparent velocity dispersions of σ v ≈350-450 km s −1 . The apparent velocity dispersions do reach σ v ≈500 km s −1 to the North of the obscured AGN's MIR position. However, from Figure 5 panel (b) we can see that regions of highest apparent velocity dispersion tend to show double peaked profiles. Thus, these high σ v values likely reflect more complicated line profile shapes rather than more chaotic kinematics. Complicated, double-peaked line profiles are expected from Lyα radiative transfer calculations, as Lyα photons produced at line center need to be either red-or blue-shifted in order to escape the core of the profile (Verhamme et al. 2006). While it is difficult to draw strong conclusions in the case of LABn06, the double peaked profiles suggest that resonant scattering of Lyα photons is important in this system, and that in some regions additional kinematics may be responsible for producing asymmetric line profiles. The conclusion that resonant scattering appears to be important and that the obscured AGN is the dominant power source together make a clear prediction that this nebula should show rising Lyα polarization as a function of radius from the AGN location (e.g., Rybicki & Loeb 1999;Hayes et al. 2011). It would therefore be interesting to obtain imaging polarimetry or spectro-polarimetric data to confirm both of these claims.
It is possible that the gas emitting from the biconical structure might have been part of, or impinged upon by, a galactic scale wind driven by massive stars (e.g., Wilman et al. 2005) or an AGN-powered radio jet (e.g., Swinbank et al. 2015). In modeling the evolution of Lyα nebulae in galaxy formation simulations, Kimock et al. (2020) found that ∼ 50% of the gas in the halo of a simulated nebula has crossed the virial radius of the central galaxy one time by z = 2, meaning it was likely part of an outflow (their Figure 17). While the stellar wind scenario will probably be less likely on such large scales, deeper low-frequency radio imaging of this region would be needed to rule out the possibility of a past interaction between a now fading radio jet and the emitting nebula.

Gravitational Cooling, Revisited
LABn06 has often been cited as an example of gravitational cooling in Lyα-emitting nebulae, but additional data continues to call this into question. Most of the Lyα emission seen in LABn06 is offset substantially from star-forming galaxies in the region and is relatively clumpy within a two-lobed structure. By contrast, Lyα nebulae that are lit up via cooling radiation in hydrodynamical simulations are centered on the region of densest star formation (Rosdahl & Blaizot 2012), and appear smoother and less concentrated than Lyα nebulae that host AGN (Kimock et al. 2020). The presence of a "hole" in the center of the Lyα nebula is particularly at odds with theoretical expectations for gravitational cooling. Furthermore, gravitational cooling nebulae should be associated with the most massive halos (Ao et al. 2020), but even in cases massive enough to produce gravitational cooling radiation, Lyα emission driven by AGN fluorescence is still expected to dominate over the lower surface brightness cooling emission (Kimock et al. 2020).

CONCLUSION
The vast majority of Lyα nebulae targeted with widefield IFS instruments to date have been selected due to the presence of a known AGN. Here, we used the MUSE instrument aboard the VLT's 8.2m telescope to simultaneously image and obtain spectroscopic coverage of a Lyα-selected nebula at z ∼ 3.2, finding that: • The source originally identified as Source 6 is likely a radio-quiet, obscured AGN/star-forming composite galaxy.
• The location of the obscured AGN lies within 17 pkpc of the light-weighted centroid of the Lyα emission, and the gas surrounding the AGN is at the kinematic center of the nebula. The nebula therefore appears to be morphologically and kinematically centered on the obscured AGN.
• The full extent of LABn06 is ≈173 pkpc, measured down to a surface brightness of 4.0×10 −19 erg s −1 cm −2 arcsec −2 . At large radii ( 25 kpc), the Lyα nebula surface brightness profile is wellfit with an exponential, resembling the profiles of other AGN-powered nebulae and stacked Lymanbreak galaxies. At small radii, however, LABn06 shows an unusual dip in its surface brightness profile, with the brightest peak offset spatially from the location of the AGN.
• LABn06 shows a level of asymmetry in between the nebulae seen around Type I and Type II AGN.
In terms of second-order moments, LABn06 is quantified as being circularly asymmetric, similar to Type II AGN, and in terms of a Fourier decomposition approach, it shows circular symmetry at moderate radial distances but a relatively asymmetric profile compared to the majority of Type I AGN at larger radial distances.
• The kinematics of the nebula show signs of coherent motion and a range of apparent line widths up to σ v ≈500 km s −1 , similar to what has been observed around other AGN-powered nebulae observed with MUSE.
• The Lyα emission lines at different locations across the nebula show double peaks and asymmetric profiles, suggesting that resonant scattering is playing a role in the system.
• AGN-powering cannot be conclusively demonstrated based on high ionization lines, but the constraints on the CIV/Lyα and HeII/Lyα ratios measured within the nebula are consistent with measurements from other AGN-powered Lyα nebulae.
These results are consistent with the obscured AGN being largely responsible for powering the extended Lyα emission in LABn06. Further observations will be needed to understand the central dip in the surface brightness profile, and whether this indicates an unusual distribution of gas or dust in this system or timevariability in the ionizing output from the AGN. Additionally, obtaining Lyα polarimetry to constrain the effects of Lyα photon scattering and low-frequency radio observations to rule out the presence of past radio-jet interactions with the gas will greatly benefit our understanding of this source.