3D tomography of the giant Ly α nebulae of z ≈ 3–5 radio-loud AGN

Ly α emission nebulae are ubiquitous around high-redshift galaxies and are tracers of the gaseous environment on scales out to ≳ 100 pkpc (proper kiloparsec). High-redshift radio galaxies (HzRGs, type-2 radio-loud quasars) host large scale nebulae observed in the ionised gas di ff er from those seen in other types of high-redshift quasars. In this work, we exploit MUSE observations of Ly α nebulae around eight HzRGs (2 . 92 < z < 4 . 51). All the HzRGs have large scale Ly α emission nebulae with seven of them extended over 100 pkpc at the observed surface brightness limit ( ∼ 2 − 20 × 10 − 19 ergs − 1 cm − 2 arcsec − 2 ). Because the emission line profiles are significantly a ff ected by neutral hydrogen absorbers across the entire nebulae extent, we perform an absorption correction to infer maps of the intrinsic Ly α surface brightness, central velocity and velocity width, all at the last scattering surface of the observed Ly α photons. We find: (i) The intrinsic surface brightness radial profiles of our sample can be described by an inner exponential profile and a power law in the low luminosity extended part; (ii) our HzRGs have higher surface brightness and more asymmetric nebulae than both radio-loud and radio-quiet type-1 quasars; (iii) intrinsic nebula kinematics of four HzRGs show evidence of jet-driven outflows but we find no general trends for the whole sample; (iv) a relation between the maximum spatial extent of the Ly α nebula and the projected distance between the AGN and the centroids of the Ly α nebula; (v) an alignment between radio jet position angles and the Ly α nebula morphology. All of these findings support a scenario in which the orientation of the AGN has an impact on the observed nebular morphologies and resonant scattering may a ff ect the shape of the surface brightness profiles, nebular kinematics and relations between the observed Ly α morphologies. Furthermore, we find evidence showing that the outskirts of the ionised gas nebulae may be ‘contaminated’ by Ly α photons from nearby emission halos and that the radio jet a ff ects the morphology and kinematics of the nebulae. Overall, this work provides results which allow us to compare Ly α nebulae around various classes of quasars at and beyond Cosmic Noon ( z ∼ 3).


Introduction
Being the most abundant element in the Universe, hydrogen (especially the cold gas, i.e. neutral hydrogen atom and molecular hydrogen, H 2 ) is the building block of the baryonic Universe.Studying H 2 directly is difficult due to lack of prominent transition lines.It is often probed using low-J CO transitions as a proxy which unfortunately results in added uncertainties, for example in the conversion factor (e.g.Bolatto et al. 2013).In contrast, neutral atomic hydrogen can be easily ionised (E H 0 = 13.6 eV) and cascade with line emissions being produced.The H i Lyαλ1216 (Lyα hereafter) line is the most promi-nent one among them.For high-redshift galaxies, it is a commonly targeted emission line which can easily be observed in the optical to near-infrared bands (e.g.Hu & McMahon 1996;Cowie & Hu 1998;Shimasaku et al. 2006;Dawson et al. 2007;Leclercq et al. 2017;Wisotzki et al. 2018;Umehata et al. 2019;Ono et al. 2021;Ouchi et al. 2020, and reference therein).Lyα emission can be detected on a range of spatial scales, for example at interstellar medium (ISM) to circumgalactic medium (CGM, Tumlinson et al. 2017) scales and even beyond the viral radius of the central object out to intergalactic medium (IGM) scales (e.g.Cantalupo et al. 2014;Cai et al. 2019;Ouchi et al. 2020).However, it is non-trivial to identify the origin of Lyα emission (e.g.due to resonant nature of Lyα emission and various potential ionising source acting at once), which is essential to understanding the physics of the emitting gas observed on different scales and around various types of objects (Dijkstra 2019;Ouchi et al. 2020).This is further complicated when AGN (active galactic nuclei) are present.
Active galaxies hosting AGN, especially the ones with quasar level activities (bolometric luminosity, L bol ≳ 10 45 erg s −1 ), at high-redshift are known to host Lyα nebulae on scales of a few 100 kpc (e.g.Heckman et al. 1991a;Basu-Zych & Scharf 2004;Weidinger et al. 2004Weidinger et al. , 2005;;Dey et al. 2005;Prescott et al. 2015;Cantalupo et al. 2014;Arrigoni Battaia et al. 2016;Borisova et al. 2016;Cai et al. 2019;Arrigoni Battaia et al. 2019).The central powerful AGN act as a main ionising mechanism for the surrounding gas which is responsible for the detection of these extended Lyα nebulae (as predicted by theoretical works, e.g.Costa et al. 2022).In addition, the diffuse emission from galaxies nearby to the AGN host can also contribute to the overall profile observed of the central target (e.g.Byrohl et al. 2021).In some of the giant nebulae, it is natural to find various mechanisms functioning at different scales and positions (e.g.Vernet et al. 2017).Therefore, despite leaving internal physics entangled, Lyα acts as a simpler tool for detecting gaseous environment through out cosmic time.
Before wide field integral field spectrographs (IFS) became available, narrow-band imaging and long slit spectroscopy provided effective methods to detect diffuse Lyα nebulae (e.g.Steidel et al. 2000;Francis et al. 2001;Matsuda et al. 2004;Saito et al. 2006;Yang et al. 2009Yang et al. , 2010;;Cantalupo et al. 2012Cantalupo et al. , 2014;;Hennawi & Prochaska 2013;Prescott et al. 2015;Arrigoni Battaia et al. 2016).However, these observations have been limited by uncertainties in the systemic redshift measurements and limited spatial coverage, respectively.Integral field unit observations (e.g.Multi-Unit Spectroscopic Explorer or MUSE/VLT and Keck Cosmic Web Imager or KCWI/Keck) allow us to measure the extent of the nebulae together with the information of their dynamics.Numerous works of Lyα nebulae around quasars report (10s of kpc to over 100 kpc) extended emission across a large range of redshifts (z ∼ 2 to z ∼ 6.3) and quasar types (e.g.radio-quiet and radio-loud type-1, radio-quiet type-2 and extremely red quasar, ERQ, Christensen et al. 2006;Borisova et al. 2016;Arrigoni Battaia et al. 2019;Cai et al. 2019;Farina et al. 2019;den Brok et al. 2020;Fossati et al. 2021;Mackenzie et al. 2021;Lau et al. 2022;Vayner et al. 2023;Zhang et al. 2023).This diversity in nebula properties suggest a range of driving mechanisms, dependencies on orientation, and demonstrate that well-selected samples are needed.Despite the effort has been done on this topic, a link between the aforementioned types and type-2 radio-loud quasars on CGM scale is missing.
Among the high-redshift quasar population, high-redshift radio galaxies (HzRGs) are a unique sample despite smaller in number (see Miley & De Breuck 2008, as a review).They host type-2 quasars and have powerful radio jets.They have been shown to reside in dense protocluster environments (Venemans et al. 2007;Wylezalek et al. 2013Wylezalek et al. , 2014;;Noirot et al. 2016Noirot et al. , 2018) ) which may evolve to modern galaxy clusters.HzRGs were among the first sources where giant Lyα nebulae were discovered (∼ 10 44 erg s −1 , ≳ 100 kpc, e.g.Hippelein & Meisenheimer 1993;van Ojik et al. 1996van Ojik et al. , 1997;;Reuland et al. 2003;Villar-Martín et al. 2006, 2007b) and observed with the previous generation of IFU instruments (e.g.Adam et al. 1997).The Lyα nebulae of HzRGs have been found to have two distinctive parts, namely the high surface brightness kinematically disturbed inner part and the quiescent low surface brightness extended outer nebula (e.g.Villar-Martín et al. 2002, 2003, 2007a).The spatial separation of these two parts seem to be consistent with the extent of the radio jets (e.g.Villar-Martín et al. 2003) suggesting that the jet plays a role in disturbing the inner part.Specifically, there is evidence that the Lyα nebulae around HzRGs are related to jet-driven outflows (Humphrey et al. 2006) while some of the quiescent gas may be related to infalling material (Humphrey et al. 2007).AGN photoionisation is likely the main mechanism of exciting these nebulae (e.g.Villar-Martín et al. 2002, 2003;Morais et al. 2017), but ionisation by fast shocks might also play a role (e.g.Bicknell et al. 2000;Morais et al. 2017).Polarisation measurements show that some of the Lyα emission in HzRGs is scattered (Humphrey et al. 2013).Despite these works, however, a comparison of the nebulae of HzRGs and other quasar samples has yet to be performed which is the motivation of this work.
The Lyα nebulae of HzRGs are known to be partially absorbed by neutral hydrogen (H i absorbers, e.g.Rottgering et al. 1995;van Ojik et al. 1997;Jarvis et al. 2003;Wilman et al. 2004;Humphrey et al. 2008;Kolwa et al. 2019).The absorbing gas is found to be extended on galaxy-wides scales and likely related to outflowing gas from the host galaxy (e.g.Binette et al. 2000;Swinbank et al. 2015;Silva et al. 2018a;Wang et al. 2021b).The correction of these absorption is only possible through spectral observation.Without careful treatment, a considerable amount (a factor of ≳ 5) of flux will be missed, and inaccurate conclusions will be derived.Alternatively, some absorption trough features might potentially be explained by radiative transfer effects (Dijkstra 2014; Gronke et al. 2015;Gronke & Dijkstra 2016;Gronke et al. 2016).Although it is interesting to compare the different treatments of the observed Lyα spectra, it is beyond the scope of this work.
There was also clear observational evidence that the morphology of the continuum and line emission regions of HzRGs are aligned with jet direction (e.g.Chambers et al. 1987;Pentericci et al. 1999;Miley et al. 2004;Zirm et al. 2005;Duncan et al. 2022) on relatively smaller scale (several kpc to 10s of kpc).Molecular gas detected around HzRGs was reported to be distributed along the jet within and outside the hot spot which may suggest several scenarii (e.g.jet-driven outflow, jet-induced gas cooling and jet propagating into dense molecular gas medium, Emonts et al. 2014;Gullberg et al. 2016;Falkendal et al. 2021).On Mpc scale, West (1991) found that the radio jet often points towards nearby galaxies.Eales (1992) proposed a model explaining the alignment effect, suggesting that the high-redshift radio emission is often detected when the jet travels close to the major axis of surrounding asymmetrically distributed gas.With the advanced IFS observation and 100s kpc gas tracer of Lyα, we are able to probe the intrinsic (i.e.corrected for absorption) gaseous nebula around HzRGs in this work, test its distribution with respect to the radio jets and seek evidence following these pioneering works.
In this paper, we utilise the power of MUSE IFU to fully map the Lyα emission nebulae of a HzRGs sample over a redshift range of 2.92−4.51and initiate a comparison with type-1 quasars and study of CGM-scale environments.We introduce our sample of HzRGs, the MUSE observations and data reduction in Sect. 2. We present how we measure the maximum extent of the nebulae in Sect.3.1 and summarise the spectral fitting procedure in Sect.3.2.We then present the results of surface brightness, kinematics and morphology in Sect. 4 followed by a discussion in Sect. 5. Finally, we conclude in Sect.6.In this paper, we assume a flat ΛCDM cosmology with H 0 = 70 km s −1 Mpc −1 and Ω m = 0.3.Following this cosmology, 1 arcsec ≃ 6.6 − 7.7 pkpc for our sample redshifts.Throughout the paper, pkpc stands for proper kiloparsec and ckpc represents comoving kiloparsec, ckpc= (1 + z)pkpc.In this paper, we use 'intrinsic' to refer to the absorption corrected Lyα emission.

Sample selection
The 8 HzRGs at 2.92< z <4.51 (Table 1) that we investigate in this paper were selected to (i) be at z > 2.9 for Lyα to be covered by MUSE ; (ii) have a known extended bright Lyα (> 10 ′′ ) emission nebula; and (iii) be at DEC < 25 • to be observable by ground-based telescopes in the southern hemisphere.This sample also has a wealth of high quality supporting data obtained by our team, including deep Spitzer/IRAC and Spitzer/MIPS 24 µm imaging, and Herschel/SPIRE detections (Seymour et al. 2007;De Breuck et al. 2010).ALMA Band 3 or 4 data are also available for the sample targeting dust continuum and molecular lines (Falkendal et al. 2019;Kolwa et al. 2023).Being identified as radio galaxies, the radio observations (e.g.VLA, Carilli et al. 1997) provide information on the jet morphology and polarisation.Based on these supporting data sets, we have estimates of the total stellar mass of the host galaxies (several 10 11 M ⊙ for all targets, De Breuck et al. 2010) and the star formation rates ranging from uppers limit of < 84 M ⊙ yr −1 to constraints of 626 M ⊙ yr −1 (Falkendal et al. 2019).

AGN bolometric luminosity estimation
To put the HzRGs into context with other quasar species, we plan to link our Lyα nebulae to literature works based on AGN bolometric luminosity.There are different methods for estimating the bolometric luminosity of AGN, L bol, AGN , for example through scaling of the far-IR AGN-heated dust luminosity (e.g.Drouart et al. 2014), scaling the IR flux density (e.g.f 3.45 µm which is used for type-1 quasars, Lau et al. 2022) and through [O iii] emission (which can be affected by star formation and/or shocks Reyes et al. 2008;Allen et al. 2008).However, there is a large uncertainty between the values derived through these different methods which makes it non-trivial to directly compare the L bol, AGN of type-1s and type-2s.For instance, the estimates for type-2 AGN are affected by obscuration by the dusty torus assuming the AGN unification model (e.g.Antonucci 1993).Accounting for this by applying an extinction correction factor would lead to a large uncertainty (e.g.Drouart et al. 2012) if we were to use the same method for type-1s to estimate the L bol, AGN for our sample.We report that the L bol, AGN estimated for our sample using those different methods varies from 10 45.9 to 10 48.5 erg s −1 .Given this large uncertainty, we find it is unreasonable to draw further conclusions from the comparison of L bol, AGN between type-1s and our HzRGs.However, it is worthwhile to report this estimation procedure and the resulted inconsistency under different assumptions.A systematic study of the L bol, AGN is beyond the scope of this work and may be done more thoroughly through multi-wavelength approach.

Jet kinematics
To distinguish between the approaching and receding sides of the jet, we use the kinematics information from [O iii] as a proxy which is often used for studying quasar outflow (e.g.Veilleux et al. 2005;Zakamska et al. 2016;Nesvadba et al. 2017a,b;Vayner et al. 2021).5 out of 8 of our sample targets have been observed by SINFONI from which the [O iii] velocity shifts are available (Nesvadba et al. 2007(Nesvadba et al. , 2008(Nesvadba et al. , 2017a)).For MRC0943-242 and TN J1338-1942, we use the radio hot spot polarisation information as indicator where the more depolarised indicates the far side (receding) of the jet (Carilli et al. 1997;Pentericci et al. 2000).These are also consistent with the tentative [O ii] velocity gradient of TN J1338-1942 found in Nesvadba et al. (2017a) (also He ii kinematics in Kolwa et al. 2023) and MRC0943-242 He iiλ1640Å (He ii) kinematics in Kolwa et al. (2019). For 4C+04.11, Parijskij et al. (2014) gives the jet kinematics based on high-resolution radio polarisation.We note here that the reported approaching and receding directions based on the current observations should be treated with caution.The polarisation of the radio lobes could especially be affected by the intervening ionised structures.We also quantified the size of the jets by calculating the angular distance between the jet hot spots on either side to the AGN position (presented in Appendix D).

MUSE observations
In this work, we analyse data from MUSE integral field spectrograph (Bacon et al. 2010;Bacon et al. 2014) mounted on the ESO Very Large Telescopes (VLT) Yepun (UT4).All observations were carried-out in Wide-Field Mode (WFM) offering a 1×1 arcmin 2 field of view and spatial sampling of 0.2 arcsec pixel −1 .MUSE provides two sets of wavelength coverage: a nominal range (N, 480−930 nm) and an extended range (E, 465−930 nm) without using of the adaptive optics (AO).For observations carried in AO mode, the wavelength coverage of 582−597 nm is excluded due to the Na Notch filter.The MUSE spectrograph has the spectral sampling of 0.125 nm pixel −1 and resolving power of 1750−3750 for 465−930 nm which corresponds to ∆v ∼ 171 − 90 km s −1 .
The observations of our sample were carried mostly in service mode under the program IDs 094.B-0699, 096.B-0752 and 097.B-0323 (PI: J. Vernet).For MRC 0943-242, we also include the data of MUSE commissioning observation under the program ID 60.A-9100(A) (e.g.Gullberg et al. 2016).The extended wavelength coverage was employed for MRC 0943-242, the lowest redshift sample target, to cover its Lyα emission (L Lyα, obs = 4769 Å).We use the MUSE commissioning and science verification data of TN J1338-1942 under the program IDs 60.A9100(B) and 60.A-9318(A) (e.g.Swinbank et al. 2015).For 4C+03.24, we adopt the data released from the MUSE WFM-AO commissioning observations under the program ID 60.A-9100(G).The information of the observations of our sample, in the order of redshift, is summarised in Table 1.For each object, observations consist of 1 (4C+03.24)to 6 (TN J 1338-1942) observing blocks (OBs).Within each OB, the 2 or 3 exposures of 20−30 minutes were slightly dithered (with a < 1 ′′ amplitude pattern) and rotated by 90 degrees from each other.

Data processing
The reduction of the raw MUSE data are carried out following the standard procedure using the MUSE pipeline (Weilbacher et al. 2020, version 2.8.4) executed by EsoRex (ESO Recipe Execution Tool; ESO CPL Development Team 2015).For studying the extended Lyα nebulae to the faintest edge, we reduce the data following the optimised procedure developed in our pilot study of 4C+04.11(Wang et al. 2021b).We first reduce each exposure Notes. (†) The redshifts are determined from the He iiλ1640Å or [C i](1-0) emission line (Kolwa et al. 2023).For MRC0316-257 which is not included in Kolwa et al. (2023), we reported its z sys in this paper based on He ii fit from our MUSE data (Appendix E).
(*) The seeing reported here is determined from the fitted 2D Moffat FWHM (full width at half maximum) of a star in the white-light image (5000 − 9000 Å) produced from the combined cube.We note that the stars used are red in color, i.e. the image quality in the Lyα wavelength should in general be worse than the reported seeing (e.g.larger by 10 to 20%). (a) The seeing is determined from a star in the overlapping region of the two pointings. (b) There is no available star in the FoV.The seeing is determined from the fit of the most point-like source.
individually with the standard pipeline doing the sky-line subtraction and then using ZAP (Zurich Atmosphere Purge, Soto et al. 2016) to remove the sky-line residuals (see below details regarding the ZAP execution).We then combine all exposures to the final data cube using MPDAF Cubelist.combine(MUSE Python Data Analysis Framework Bacon et al. 2016).We correct the astrometry of the final combined cubes using star positions from the available Gaia EDR3 catalogue (Early Data Release 3, Gaia Collaboration et al. 2021).Two sources had no GAIA star within the MUSE field-of-view (FoV).For TN J0121+1320 we use the SDSS DR16 (16th Data Release, Ahumada et al. 2020) catalogue instead.For MRC 0316-257, we use Gaia EDR3 to first correct the astrometry of the HST/ACS F814W image and then matched the MUSE cube to the HST image.
Using ZAP directly for sky-line residual removal without applying masks may remove faint narrow Lyα line emission at the outskirt of our sample.Since the Lyα nebulae in our sample extend much further beyond the continuum emission regime of the host galaxy and become narrower in line width (e.g.Villar-Martín et al. 2003;Humphrey et al. 2007) such that they are mistakenly treated as sky-line residuals and removed.To alleviate this problem (Soto et al. 2016), for each source, we (i) generate a first version of the combined data cube without masks in the ZAP step; (ii) construct a Lyα mask that covers most line-emission region1 ; (iii) re-run ZAP using this Lyα mask on individual cubes for each exposures; (iv) combine the newly obtained individual cubes to the final version data cube with MPDAF.
We also correct for small residual (mostly) negative background level offsets probably due to a slight over-subtraction of the sky continuum in previous steps.To do so, we (i) extract a median spectrum from an r ≃ 10 ′′ circular aperture around the radio galaxy masking all continuum sources falling in the aperture; (ii) mask the Lyα line emission wavelength range and strong sky-lines (> 10 16 erg s −1 cm −2 Å −1 arcsec −2 , Hanuschik 2003) for this median spectrum; (iii) fit a 6th-order polynomial to this masked spectrum; (iv) subtract this solution from the whole cube.
Finally, to correct for the known underestimation of the variance in the standard pipeline reduction (see Weilbacher et al. 2020), variance scaling is implemented as described in Wang et al. (2021b).Specifically, we scale the variance extension propagated by the pipeline based on the scale factor calculated in source-free regions using the variance estimated from the data extension.

Lyα nebulae extent and tessellation
To systematically study the Lyα nebulae of our HzRGs sample, we first need to determine all the voxels (volume pixel) containing usable Ly-alpha signal (Sect.3.1.1)and bin the data to a sufficient S/N using a tessellation technique (Sect.3.1.2)before fitting the emission feature described in Sect.3.2.

Maximum extent of the nebulae
To select the Lyα signal with optimised sensitivity and capture the very low surface brightness structures of the nebulae, we used our own version of the adaptive smoothing algorithm described in Martin et al. (2014) (see also Vernet et al. 2017, for an application to one of the sources in our sample).We first smooth the data cube in the wavelength direction by averaging n λ neighbouring pixels.Then for each wavelength plane, the algorithm iteratively smoothes spatially with a growing gaussian kernel selecting pixels passing a given SNR threshold (T S/N ) and leaving to the next iterations only spaxels below this S/N threshold, until a maximum smoothing radius is reached (σ max ).The spaxels not selected by the end of the iterative process are masked out.To further clean the smoothed data cube from spurious noise features and make sure that a proper line fitting can be made, we mask spatial positions selected by the adaptive smoothing algorithm in less than n c consecutive wavelength bin.Notes. (†) Position of the AGN and/or host galaxy. (*) The surface brightness limit is the 2σ limit extracted from continuum-source-and Lyα-free regions in a narrow band image.The narrow band image is collapsed from v 05 to v 95 (see Section 4.1).The cgs unit is erg s −1 cm −2 arcsec −2 .
( ‡) Intrinsic Lyα luminosity (i.e.corrected for absorption) of the nebula.It is integrated over the entire area selected in Sect.3.1 and multiplied with the luminosity distance using the cosmological parameters (Sect.1). (§) Observed Lyα luminosity of the nebula integrated over the entire area from v 05 to v 95 (see text).
To determine the optimal combination of the 4 parameters (n λ , T S/N , σ max and n c ), we explore a range of possible combinations and select the set that is most sensitive to the extended low-surface brightness emission while at the same time minimising the number of detached 'island-like' structures (see Appendix A.1 for details).We note that the maximum nebulae extents selected by this method are similar to the results from previous studies of individual targets by different procedure (TNJ1338-1942 from Swinbank et al. 2015) or pure manual selection (MRC0316-257 in Vernet et al. 2017).We then manually clean up this map for the few remaining isolated island-like regions with further checking spectra extracted from these regions.This clean-up is accompanied by signal checking through spectrum extraction and only affects low S/N regions (Appendix A.1). Thus, the bulk of the detection map remains unchanged.This resulting detection map defines the pixels that we consider as part of the nebula and that we use in the analysis in this paper (see also Appendix A.1).

Tessellation procedure
In order to increase the S/N to a level that allows fitting of the Lyα line, especially close to the detection limit at the periphery of the nebulae, we tessellate the Lyα detection map.To construct the tessellation map, we firstly use a S/N map based on a narrowband image (∼ 15 Å wide) extracted around the Lyα emission peak.We implement a two-step Voronoi binning (Cappellari & Copin 2003) procedure which optimises the performance for both high S/N and low S/N regions by tessellating individually on these two parts.Specifically, the two-step procedure uses different target S/N for inner and outer regions.In this way, we can avoid large size tiles at the low S/N (outer) regions which may unnecessarily smear spatial resolution by imposing too high target S/N.We then combine the tessellated regions from the twostep process into one map.We emphasise that the tessellation is a trade-off between spatial resolution and S/N.The main goal of the work is to study the extend Lyα nebulae to the detection limit.This can only be achieved by sacrificing the spatial information.The details of this tessellation process are described in Appendix A.2, and in A.3 we present the resulting the maps.

Lyα absorption modeling
In this work we treat the Lyα emission system of HzRGs as an idealised case where several assumptions have been made prior to the analysis: (i) the radio galaxies reside in giant reservoirs of neutral hydrogen (∼100s kpc); (ii) the neutral hydrogen is rather diffuse with large covering factor; (iii) the geometry of the giant reservoirs is unknown but can be highly asymmetrical due to the influence of the radio jet.Under these assumptions, it is natural that we observe absorption effecting the Lyα profiles.Indeed, such absorption troughs are observed in our Lyα spectra (see Fig. 1) that need to be accounted for when drawing conclusions about the intrinsic emission line flux and higher moment measurements.Specifically, high resolution spectrocopy using the UltraViolet and Echelle Spectrograph (UVES) on the VLT exists for seven out of eight of the targets in our sample (Jarvis et al. 2003;Wilman et al. 2004, Ritter et al. in prep.).These spectra with ∼ 30 higher resolution than MUSE display sharp edges which is fully consistent with a well-defined absorption profile rather than radiative transfer effects.We note that the the term 'radiative transfer' used in the paper refers to the process where Lyα photons are scattered in frequency (wavelength) but are still captured in the spectrum (i.e.not 'lost' in the observer's line of sight).In our assumptions, contrary to this, the photons are 'lost' either due to being scattered outside our line of sight or being absorbed by dust and remitted at longer wavelength.We therefore adopt the technique used in our pilot study (Wang et al. 2021b, and equations therein) and fit the spectra using a combination of Gaussian emission line profiles and Voigt absorption troughs (e.g.Tepper-García 2006, 2007;Krogager 2018).This procedure has also been implemented successfully in the literature for fitting the Lyα line emission in HzRGs (e.g.Swinbank et al. 2015;Silva et al. 2018b;Kolwa et al. 2019).
The known degeneracy between the H i column density and Doppler parameter in our fits (e.g.Silva et al. 2018a) does not affect the reconstructed intrinsic emission which is the focus of this work.We show the 'Master Lyα spectrum' extracted from a central 1 ′′ aperture in Fig. 1a which presents how the intrinsic profile compares to the observed spectrum (see Sect. 4.1 for details).1, Kolwa et al. 2023).(b) Intrinsic Lyα surface brightness map.The flux in each tile is the integrated flux of the line emission corrected for absorption, i.e. total flux of the one or two Gaussians, see Sect.3.2.The light blue circle shows the aperture where the Master spectrum is extracted from.Green triangles mark the positions of the radio lobes.We place a green bar linking the triangles on TN J0121+1320 to indicate the unresolved state of its radio emission.The length of the bar represents the linear size of the 3σ contour along the east-west direction.The white hatched regions are the ones where the flux uncertainty is higher than 50% of the fitted intrinsic flux.The white bar indicates the 50 pkpc at the redshift of the radio galaxy.The unit of the surface brightness is 10 −16 erg s −1 cm −2 arcsec −2 .We apply the same color scale for all targets.(c) v 50 map of the intrinsic Lyα nebula.The zero velocity used for each target is determined by the systemic redshift (Table 1).Green contours show the morphology of the radio jet in arbitrary values.The green cross mark the AGN position (Table 2).We emphasise that the term 'intrinsic Lyα emission' throughout the work refers to the nebula Lyα emission corrected for intervening absorbers.The absorption troughs seen on the spectra (Fig. 1a) are due to the Lyα emission being absorbed by these neutral hydrogen gas clouds or shells along the line of sight.Under the aforementioned assumptions, a natural consequence is that the absorbers must be distributed across the whole projected extension of the nebula.The fact that we mostly ob-serve these features continuously across the extent of the nebulae in most HzRGs indeed indicates they are coherent intergalacticscale structures.This can be found in Fig. 2, 3 and 4 where similar absorption features are seen in the selected spectra at larger distance (10s of kpc) away from AGN.Similar maps of the remaining sources are shown in Appendix C which are the ones have been previously published (Swinbank et al. 2015;Gullberg et al. 2016;Vernet et al. 2017;Falkendal et al. 2021;Wang et   2021b).Our approach is a common interpretation in studies of HzRGs.Conversely, such absorbers are not often seen in the Lyα nebulae of other quasars.This reinforces the interpretation that strong (radio-mode) feedback on intergalactic scales is needed to create such 'shells' of H i material.The use of a Gaussian as underlying intrinsic emission profile is supported both by observational and modeling works (e.g.Arrigoni Battaia et al. 2019;Chang et al. 2022).This could be a result of prior radiative transfer effects of Lyα (e.g.local scattering or scattering from the broad line region of the AGN Verhamme et al. 2006;Gronke & Dijkstra 2016;Gronke et al. 2016;Li et al. 2022).The radiative transfer modelling requires assumptions about the composition and geometry (and kinematics) of the gas near the AGN which is not the focus of this paper.Hence, we just assume the Gaussian shape of the Lyα (which could be due to the radiative transfer effects) and correct for the absorption along the line of sight to reconstruct the intrinsic emission on CGM scales.Incorporating radiative transfer calculations into the study of HzRGs Lyα nebulae is beyond the scope of this current work.Further developments of theoretical works are required (e.g.adding jet and resolving shells in simulations), and our dataset would be well suited for such studies.We therefore stress that the presented results are only valid for the stated assumptions that absorption rather than radiative transfer is primarily responsible for the line profiles.We will discuss the limitation of this treatment in Sect.5.1.We note that the O v]λλ1213.8,1218.2 (O v]) line underneath the Lyα can affect the obtained flux especially in the nuclear region where the ionisation parameter (and metallicity) is higher (Humphrey 2019).In our pilot study (Wang et al. 2021b) of 4C+04.11,we found the contribution from O v] is negligible.Hence, we do not further include O v] in our line fitting.We leave the inspection to future work when data of metal lines (e.g.N vλλ1238, 1243 which is found to be related with O v]) and high resolution spectra are analysed.

Fitting procedure
To reconstruct the intrinsic Lyα emission across the nebula, we fit each spectrum in each tessellation bin (see Sect. 3.1.2)following the procedure described in Sect.3.2.1.We take into account the physical connection between neighbouring tiles by using the fit results of a previous connected bin as the starting parameters for the next bin (see Appendix A for details on the ordering).We determine the number of absorbers based on the Master spectrum where the S/N is the highest (Fig. 1a).We then use that same number of absorbers across the nebula, where the centroid, column density and Doppler parameter of absorbers are fitted in a given range (Appendix B).This assumption is supported by the profile shapes at the largest spatial extents (see Fig. 2, 3, and 4 and also Appendix C).We note that the number of absorbers selected here may be incomplete but this has minor effects on the results of this paper: (i) the absorbers that impact most the intrinsic flux (i.e.spatially extended ∼ 10 ′′ and having higher column density and/or larger Doppler parameter) are included; (ii) absorbers that seem to be 'superfluous' at the wings have only minor effects on the reconstructed flux where S/N is low (Fig. 2, 3 and 4).Future work using high spectral resolution data will address these issues also taking into account that some of these absorbers have counterparts in metal lines covered by the MUSE data (e.g.N vλλ1238, 1243 and C ivλλ1548, 1551 Kolwa et al. 2019;Wang et al. 2021b).We perform the fit in each bin using both one and two Gaussian emission line components and we choose the solution that minimises the reduced χ 2 .The fit is done using a Least-squares method followed by a Markov chain Monte Carlo sampling (MCMC, emcee Foreman-Mackey et al. 2013).The uncertainties we report are either the direct output of the 1σ error by the MCMC or the propagated 1σ error.A detailed description of the fitting procedure is provided in Ap-pendix B. We reiterate that we do not report any further parameters on the absorption features which will be analysed in future work in combination with higher spectral resolution data (e.g.Jarvis et al. 2003;Wilman et al. 2004;Kolwa et al. 2019, Ritter et al. in prep.).We present the results of this procedure for all of our sources in Fig. 1.

Intrinsic mapping
In this section, we present the intrinsic maps (i.e.corrected for absorption) constructed following the fitting procedure described in Sect.3.2.For each sources we show the Master spectrum together with its best fit in Fig. 1a as an example (Sect.3.2.1).We also show the non-resonant He ii spectrum extracted in the same aperture (green histogram, not continuum subtracted) which is used for systemic redshift (∆v = 0 km s −1 , Table 1) determination.We note that there is no He ii detected at the AGN position for 4C+03.24(Sect.4.2.4).In addition, to illustrate how fitting procedure works spatially (Sect.3.2.2), the The intrinsic Lyα surface brightness maps are shown in Fig. 1b on the same flux scale.Regions with larger fitting uncertainties (≳ 50%) that should be treated with caution are indicated by the overlaid hatched tiles.We report the total intrinsic Lyα luminosities (L Lyα,int ) of the nebulae in Table 2 and their maximum linear extent, d max , in Table 3. Down to the surface brightness limit (Table 2), seven of our nebulae are extended over 100 pkpc with the largest being ∼ 347 pkpc (MRC0316-257).TNJ0121+1320 is the only target with nebula < 100 pkpc (∼ 72 pkpc).The total intrinsic surface brightness (L Lyα,int ) of the nebulae ranges from 2 to 29 × 10 44 erg s −1 .
To characterise the kinematic information of the intrinsic nebulae that are fitted with one or two Gaussians, we use a set of non-parametric emission line measurements (see e.g.Liu et al. 2013) derived from the cumulative line flux as a function of velocity Φ(v) defined as: where f (v ′ ) is the flux density at v ′ .The often used v 50 is the velocity where the cumulative flux reaches 50% of the total integrated value, Φ(v 50 ) = 0.5Φ(∞).The v 05 , v 10 , v 90 and v 95 are defined similarly.The line width measurement, W 80 , defined in this context is W 80 = v 90 −v 10 .In case of single Gaussian fits, W 80 is directly related to the FWHM and v 50 is the Gaussian centroid.The non-parametric velocity shift (v 50 ) and line width (W 80 ) of the nebulae are shown respectively in panels (c) and (d) of Fig. 1.The v 50 maps do not show clear trend on larger scale (i.e.beyond the jet hot spots) for the whole sample.This is foreseeable given that (i) Lyα is a resonant line which is sensitive to scattering (i.e. it will not necessarily show the bulk velocity of the gas), and we only observe the last scattering surface; (ii) the size of the tile far from the centre is larger which could smear out potential velocity structures; (iii) the line emissions on several 10s of pkpc could trace the inflowing gas (or other gas components not governed by the host galaxy and/or kinematically related to the quasar outflow, e.g.Vernet et al. 2017).Within the extent of the radio jets, 3 targets (MRC0943-242, MRC0316-257 and TN J1338-1942) show tentative velocity gradients consistent with the jet kinematics (Sect.2.1.3).For the line width maps, W 80 , 3 targets (4C+03.24,4C+19.71and TN J1338-1942) show a trend with the line being broader near the centre and becoming narrower outwards.There are some tiles on the periphery of the nebulae, for example the southwest tile of 4C+04.11,displaying larger W 80 values (≳ 2500 km s −2 ).Except 4C+19.71,all targets show a line width of ∼ 800 − 2500 km s −1 .For 4C+19.71, due to the strong 5577 Å sky-line located close with the observed Lyα peak wavelength, its line width should be treated as lower  1998, Sect. 4.2.4).The yellow shaded regions show the wavelength range excluded due to the 5577 Å sky-line.limit (≳ 600km s −1 ) especially for the tiles in the outskirts of the nebula.
We note that the non-parametric measurements used in this mapping are based on intrinsic (= absorption-corrected) line profiles which are determined through model fitting, same as Wang et al. (2021b).In Appendix C, we present the maps of observed surface brightness and flux ratio as supplementary material.

Circularly averaged surface brightness radial profiles
In this section, we present the surface brightness radial profile of the eight Lyα nebulae.In order to compare our HzRGs to other quasar samples, we extract the surface brightness profile centred around the AGN in circular annuli.The annuli over which the profiles extracted are shown in Fig. D.1.We compute the surface brightness in each annulus as the mean of the surface brightness of each contributing spaxel weighted by the fraction of the spaxel area covered by the annulus.Table D.2 lists the extracted intrinsic profile values.
Figure 5 shows the radial profiles after correction for cosmological dimming and in comoving units for observed (upper panel) and intrinsic (lower panel, corrected for absorption) maps.The dashed lines in the upper panel represent the comparison quasar samples or single targets (ERQ and 2 radio-loud quasars from Lau et al. 2022;Vayner et al. 2023, respectively) ).Those two are on the higher surface brightness end of the comparison samples, and we will examine them quantitatively along with both the intrinsic and observed HzRGs profiles.We note that Vayner et al. (2023) fitted the Lyα absorbers from the spatially integrated 1D spectrum and found ∼ 10 13.5 cm −2 for the column densities.For Farina et al. (2019), there is not much evidence of absorption.There-  2) for each target in the same color.Lower panel: Intrinsic radial profile in thicker solid lines.The shaded regions around each profiles indicates the uncertainty range of the surface brightness from fitting.In this panel, we show again the same profiles of F19 and V22 7C as in the upper panel for comparison.Our HzRGs are extended further with higher surface brightnesses (or flattening in some sources) at larger radii (∼ 300 ckpc) compared to other samples.fore, the comparison is legitimate.The best fit profiles to the observed Lyα nebulae of radio loud quasars in Arrigoni Battaia Article number, page 11 of 43 Fig. 6.Surface brightness radial profiles for approaching (blue squares) and receding (red circles) directions along the jet axis.The dotted curves in corresponding colors show the exponential+power law fits for the two directional profiles.We also include the fits for the circularly averaged profile in solid magenta lines.In each panel, the magenta shaded region mark again the same uncertainty range for the intrinsic surface brightness profile as Fig. 5.The solid green curve is the normalised radial profile of a star extracted up to 2 ′′ (the one in the FoV of MRC0316-257 is extracted from a round galaxy due to no available star) showing the PSF (Table 1).The vertical dashed lines indicate the distances of the jet hot spots in corresponding colors.The profile along the receding side of the jet is brighter than along approaching side for most sources within the extent of the jets except 4C+03.24and 4C+04.11.This may indicate different gas density distribution (see Sect. 4.2.2).We also identify flatting of the profile at ≳ 100 ckpc for MRC0943-242, MRC0316-257 and TNJ1338-1942 which may related to nearby companions (see Sect. 5.4).et al. ( 2019) are included in both panels of Fig. 5 which can be used as a reference between the two panels.
Except for the ERQ from Lau et al. (2022) (which is also highly obscured), type-1 radial profiles are dominated by direct emission from bright AGN point source in the inner regions (∼50 ckpc or ∼10 pkpc).Hence, due to point spread function (PSF) subtraction, the inner-most radius covered in the comparison samples is limited to ∼50 ckpc in most cases (except Vayner et al. 2023).At larger radii, the contamination by the PSF should be negligible.Of the three single target profiles, V22 7C (7C 1354+2552 from Vayner et al. 2023) has the highest surface brightness.At a radius of ∼50 ckpc, the intrinsic surface brightness of our HzRG sample has a factor of 0.5 − 7 compared with V22 7C (7 of our targets are brighter).This source then shows a faster drop off compared with HzRGs.At the faint end corresponding to ∼ 300 ckpc (except TN J0121+1320), the HzRGs have a factor of 7 − 100 higher surface brightness than V22 7C.The profile of Farina et al. (2019) shows the highest surface brightness among the comparison samples.At ∼50 ckpc, the intrinsic HzRG profiles are still a factor of 1.1 − 15 (or 4 − 40 at ∼ 400 ckpc, except TN J0121+1320) brighter than the 75th percentile of Farina et al. (2019).These indicate that our eight observed HzRG have some of the brightest known Lyα nebulae (Sect.2.1).We note that the jet compression is also known to result in high Lyα nebula luminosity (e.g.Heckman et al. 1991a,b).Compared to quasars with similarly deep observations (i.e.avoiding the surface brightness detection limit), our HzRG sample generally maintains a high surface brightness out to larger radii (5 out of 8 > 500 ckpc).We note again that the detected extent of our nebulae will have similar range even if adopting other detection methods than the ones in this paper (Sect.3.1).For example, Gullberg et al. ( 2016) reported the simi-lar extend Lyα nebula in MRC0943-242 with less exposure time.Vernet et al. (2017) detected the nebula of MRC0316-257 > 700 ckpc based on visual detection.Swinbank et al. (2015) found the > 500 ckpc nebula of TN J1338-1942 with (or even without) a simpler binning algorithm.Hence, we are sure that the detection of the ≳ 500 ckpc nebulae in our sample based on our method is robust.However, we do caution that this sample is not representative since they are selected to have bright and extended Lyα emission.The profiles of MRC0943-242, MRC 0316-257 and TN J1338-1942 show a flattening at r AGN > 200 ckpc.For the comparison samples, their profiles drop off monotonically and drop below detection limit at radii smaller than our HzRGs.The lowest surface brightness of HzRG intrinsic profiles is ∼ 1 × 10 −15 erg s −1 cm −2 arcsec −2 (MRC0316-257, corrected for cosmological dimming) which is higher than the faint end of the quasar samples by a factor of 5 − 40 (not at similar comoving distances).These indicate that we are observing some of the most extend Lyα nebulae, in two cases (MRC0316-257 and TN J1338-1942) even extending beyond the field of view of MUSE.By simply comparing our intrinsic profiles to the exponential and power law fits of Arrigoni Battaia et al. (2019), we find that the inner part of HzRGs profiles are exponential-like (especially MRC0943-242 and TN J1338-1942) while extended parts show power law decline.We note, however, that the exponential part is affected by seeing smearing.
If we do not correct for the Lyα absorption and instead measure at the observed radial profiles (Fig. 5 upper panel), 5 of the HzRGs are still brighter than the comparison samples, but by a lower factor of ∼ 2−4 (∼ 2−6) at radii of ∼400 (∼50) ckpc compared to the 75th percentile of Farina et al. (2019).Comparing the results from intrinsic and observed profiles, this suggest that Fig. 7. Directional v 50 and W 80 profiles for approaching (blue squares) and receding (red circles) sides along the jet axis extracted from the intrinsic maps (i.e.corrected for absorption, Fig. 1cd).The vertical dashed lines indicate the distances of the jet hot spots (blue for approaching, red for receding, Appendix D).We note that the radio emission of TN J0121+1320 is unresolved.The grey shaded regions are >5 arcsec from the host galaxy.The data points in the shaded regions should be treated with caution given the large tile size may smear kinematic structures.The horizontal black dotted line in the v 50 panel marks the 0 km s −1 derived from systemic redshift.The dashed horizontal black line in the v 50 panel of 4C+03.24indicates the velocity shift of Hβ redshift (z Hβ ≃ 3.566, −1100 km s −1 with respect to the systemic redshift) (Nesvadba et al. 2017b) with respect to its systemic used in this paper (Table .1, see text).We note that the range of the x-axis is customised for each target and that the W 80 is shown in logarithmic scale.We show a zoom-in view of the central part of the v 50 profiles of MRC 0316-257 and TN J1338-1942 in the insets.For MRC 0943-242, MRC 0316-257, TN J0205+2242 and TN J1338-1942, we mark the fit to the v 50 profiles within the jet hot spots (vertical lines) in green lines to guide the eye of the evidence of nebula velocity gradient following jet kinematics.In general, there is no clear evidence of a trend in bulk motion identified for the whole sample.For some targets (4C+03.24,4C+19.71and TN J1338-1942), W 80 decreases with increasing radial distance which may indicate that the jet is disturbing the gas.
the quasar samples may miss a non-negligible amount of flux (≳ 5) due to uncorrected for absorption.
The radial surface brightness profiles of the comparison samples are extracted from a fixed velocity or wavelength range, for example ±2000 km s −1 in Cai et al. (2019), 30 Å in Arrigoni Battaia et al. ( 2019) and ±500 km s −1 in Farina et al. (2019).Considering the redshift difference between these samples, the integration range adopted are consistent.For our study, particularly for the observed radial profile, our extraction is based on the v 05 and v 95 which are determined based on intrinsic fitting (Sect.4.1).In this way, we can minimise the uncertainties coming from the observed line width difference, for example between the emission lines in the vicinity of the host and outskirts of the nebula.Our velocity range (v 05 and v 95 ) used is basically the value of W 90 which has the range of ∼ 800 − 2700 km s −1 for all tiles of all targets.Nevertheless, we conduct a check by extracting observed circular radial profiles through integration of 30 Å around the systemic redshift of our targets for comparison.The results vary by ∼ ±10% in each annulus to the profiles in Fig. 5 (from v 05 and v 95 ), especially for emissions at > 50 ckpc where the line width is narrower comparing to the centre.The 30 Å extracted profile could be 40% less than the v 05 − v 95 extraction in the centre regions (wider line width) for high-redshift targets where the fixed wavelength range in observed frame corresponds to a narrower rest frame range.Hence, to alleviate this problem brought by the difference in line width and redshift range spanned by our sample, we keep the v 05 -to-v 95 extraction.As for the intrinsic radial profile, it is redundant to integrate from a narrower range when we can have the direct fit results for the integrated Lyα line.We show the flux ratio between the intrinsic and observed maps in the same velocity range in Appendix C which can be used as a proxy for scaling between the two profiles.Therefore, the different wavelength (velocity) ranges used when extracting radial profiles for our study and comparison quasar samples will not bring additional discrepancy besides the relatively large surface brightness value in our sample.

Directional surface brightness profiles
Since the shapes of our Lyα nebulae are asymmetric (Sect.4.3), the radial profiles extracted in Sect.4.2.1 smear out directiondependent features.For instance, several HzRG Lyα nebulae display features aligned with their radio jets, such as having higher line width and elongated morphology along the jet axis (e.g.van Ojik et al. 1997;Villar-Martín et al. 2003;Miley et al. 2004;Zirm et al. 2005;Humphrey et al. 2007;Morais et al. 2017).Hence, in this section, we study the radial profile of the intrinsic Lyα emission along the direction of the radio jet which could exert a kinematical and/or electromagnetic influence on the surrounding gas.Due to the limited S/N, we split our nebulae into two half parts (approaching and receding, Appendix D) along the jet direction and extract the surface brightness profile in each direction using the same annuli as Sect.4.2.1.Figure 6 shows these directional profiles.We also show the position of the jet hot spot for the receding (red) and approaching (blue) side with vertical dashed lines.2Qualitatively speaking, the surface brightness on the receding side is higher than on the approaching side within the radio jet extent for most of our sources (except 4C+03.24and 4C+04.11).In three sources, the receding jet hot spot is closer to the AGN: MRC0316-257, TN J0205+2242 and TN J1338-1942.This result was first reported by McCarthy et al. (1991) where the authors found that the line emission is brighter in HzGRs on the side with shorter radio jet.They interpreted this as a large scale asymmetry in the density of gas on either side of the nucleus: the denser gas absorbs more ionising radiation resulting in brighter emission lines, while the radio jet is more contained as it travels more slowly through the denser medium.

Fitting the surface brightness profiles
To quantify the shape of the profiles, we fit the circularly averaged intrinsic profile and two directional intrinsic profiles with a piecewise function split into an exponential for the inner part and power law for the outer part.This can be mathematically represented by where r h is the scale length of the exponential profile, r b is the distance at which the inner and outer profiles separate and C e and C p are normalisation parameters for exponential and power law profiles, respectively (C p = C e exp(−r b /r h )/r α b ).The determination of the piecewise function is motivated by previous studies of quasar Lyα nebula (e.g.Arrigoni Battaia et al. 2019;Cai et al. 2019;den Brok et al. 2020) which fit the profile by either power law or exponential.We also test to fit our profiles use only one of the two functions.The single profile, however, cannot fit some targets well.For example, the reduced χ 2 are high (> 20) for MRC0943-242, MRC0316-257 and TN J1338-1942 with the single-function fit.We therefore fit all of the profiles with the piecewise function for consistency.Figure 6 shows the fits and Table D.3 presents the fitted parameters.
For most of our targets, the two directional surface brightness profiles are similar to the circularly averaged profile.One exception is the approaching side of MRC0316-257 which has ∼ 1 dex lower than the receding side.This could be an extreme case of uneven Lyα emitting which may trace the different gas distribution.In Fig. 6, we also show the distance of the jet hot spots on both directions (Appendix D).There is no correlation between the distance of the jet hot spot and r b (nor r h ).As Sect.4.2.1 described, our HzRGs are high in surface brightness (large C e ); the reasons for this include (i) our sample is composed of HzRGs with bright Lyα emission, (ii) our profiles are absorption corrected, (iii) the quasar surface brightness is extracted from a fixed wavelength range (Sect.4.2.1).The exponential shape is also seen in other quasar samples, for example Arrigoni Battaia et al. (2019), Farina et al. (2019), den Brok et al. (2020) and Lau et al. (2022).The r h values derived for our sample are mostly < 20 pkpc (Table D.3) which is consistent with the quasars.This suggests a similarity between the central (high surface brightness) part of HzRGs to other quasars (type-1 radio-loud and radio-quiet, type-2 radio-quiet), despite the high surface brightness in our sample.We note that the PSF-subtraction of quasar samples and resolution effects will impact the inner part to the profile.We will further discuss the power law declining (flattening) part of our nebula in Sect. 5 combining the information from nebular morphology (Sect.4.3).

Radial profiles of kinematic tracers
It is of interest to study how the nebula kinematics changes radially which may offer evidence of outflow and/or inflow (e.g.Humphrey et al. 2007;Swinbank et al. 2015;Vernet et al. 2017).We stress that it is beyond the scope of this work to separate different Lyα kinematics emission components (e.g.systemic and outflow) which will be inspected through high resolution spectroscopic data.Hence, we adopt the v 50 and W 80 parameters to measure the overall kinematics of the line emitting gas (Sect.4.1).We caution that the kinematics derived in this way may be biased, for example if there are several gas components with different kinematics but on similar flux levels.
In Fig. 7, we show the directional radial profiles of v 50 (Fig. 1c) and W 80 (Fig. 1d), respectively.The profiles are extracted in a similar way as the directional surface brightness in Sect.4.2.2 by splitting the map into two halves (approaching and receding).The v 50 (W 80 ) value shown at each radius is averaged in the corresponding annulus.Within the extent of the radio jet hot spots (vertical dashed lines), MRC0943-242, MRC0316-257, TN J0205+2242 and TN J1338+1942 show evidence of jetdriven outflows (e.g.Nesvadba et al. 2008Nesvadba et al. , 2017a) ) if we ignore the absolute v 50 value but focus on the relative gradient.That is to say the velocity shift at the approaching side is higher than the receding side.For these four targets, we overplot a solid green line to show the fit of the velocity radial profile within the radio jet extent in Fig. 7.The same velocity gradient is also identified in He ii for MRC0943-242 (Kolwa et al. 2019), MRC0316-257 (Appendix E) and TN J1338+1942 (Swinbank et al. 2015).There is no other evidence from the v 50 of ordered gas bulk motion for the overall sample.This further suggests that Lyα is an unreliable tracer of kinematics at least on 10s to 100s pkpc scale in HzRGs.We note that the tessellation implemented, especially for tiles with larger size (∼ 5 arcsec 2 ) which are usually located in the low S/N region away from the host galaxy, may smear out potential kinematic features.One possible consequence of combining different kinematic components is broadening of the line width.This may be the case for the receding side of MRC0316-257 and both sides of TNJ0121+1320 and 4C+04.11.In general, the W 80 does not show an increasing trend toward larger r AGN .However, if the line width decreases intrinsically away from the AGN, this will counteract the broadening which makes it difficult to check the impact of smearing.Therefore, we mark the regions with r AGN > 5 ′′ on the kinematic radial profile using grey shade to flag the possible high uncertainty in Fig. 7.
If we assume that the bulk of the gas resides in the potential well of the radio galaxy, we expect to see the Lyα emission gas centred around systemic velocity, at least in the vicinity of the AGN.Offsets of the v 50 levels at r AGN ∼ 0 ckpc from 0 km s −1 (based on systemic redshift, Table 1) are identified for most of our targets which may be due to the aforementioned bias from different kinematic components and scattering of Lyα photons.The most noticeable case is 4C+03.24which has an offset of ∼ 900 km s −1 .We note that its systemic redshift (Table 1) is based on [C i](1-0) emission (Kolwa et al. 2023) due to lack of He ii from the AGN position in the MUSE data (Fig. 1a).This offset can be eased if we use the redshift of Hβ, z Hβ ≃ 3.566, reported by Nesvadba et al. (2017b) as zero velocity.It is marked in black dashed horizontal line in the v 50 panel of 4C+03.24 in Fig. 7.This corresponds to −1100 km s −1 with respective to the systemic velocity shift used in this paper.We caution that, however, the Hβ was not exclusively extracted at the AGN position (radio core).There is also a known jet-gas interaction in the south of the AGN (see bend of the radio jet contours in Fig. 1b and also van Ojik et al. 1996).From the K−band image (van Breugel et al. 1998), we can find a second continuum emission peak in the south.The position of this emission peak is marked by the red square in Fig. 4. Given these pieces of evidence, we propose that there is a companion galaxy at z Hβ ≃ 3.566 in the south of our radio galaxy (z = 3.5828).If there is (or was) an interaction between these two galaxies, the companion may have deprived gas from the AGN resulting in a gas poor AGN host (no detection of He ii and less molecular gas detected at the AGN position, Kolwa et al. 2023).The companion then becomes sufficiently massive and gas-rich to deflect the jet.Therefore, the Lyα nebula of 4C+03.24 may trace the CGM of the companion galaxy.Scheduled JWST data (Wang et al. 2021a) will offer a clearer view of this particular situation.
For the W 80 radial profiles in Fig. 7, we can first identify that most of the HzRGs have high W 80 even at larger radii (∼ 1500 km s −1 ).The exception is 4C+19.71whose measurement is affected by sky-line residuals (Sect.4.1).The W 80 reported here is similar to FWHM (especially at large radii, > 100 ckpc or ≳ 22 pkpc) where most of our fit are done with a single Gaussian (Sect.3.2).In 4C+03.24,4C+19.71and TN J1338-1942, we can see a clear radial decrease of W 80 along both directions.This may be related to results found in Villar-Martín et al. (2003) who observed a Lyα FWHM drop off at distance beyond the extent of the radio jets in a sample of HzRGs (including MRC0943-242) using deep Keck long slit spectroscopy.In our study, however, we firstly do not find such a decrease in all targets.The grey shaded regions (r AGN > 5 ′′ ) should be treated with caution.We note that the FWHM in Villar-Martín et al. (2003) was derived without correction for absorption.Since we see a high spatial coverage of absorbers (Fig. 2, 3 and 4), the correction indeed helps with recovering the close-to-intrinsic gas kinematics at large radii.In Fig. D.2, we show the position-velocity diagram extracted based on our observed and intrinsic surface brightness maps along the jet as a direct comparison with the long slit spectroscopic study.Although it resembles the detection of Villar-Martín et al. (2003) at first glance, we note that this is due to the tessellation and the contrast between high surface brightness and low surface brightness part.
By considering both the radial profiles from v 50 and W 80 , we can generally find that the profiles within the jet extents have behave differently compared to the profiles outside the jet extent.This again suggests the jet is disturbing or entraining the Lyα emitting gas.There are unclear signs of kinematics other than outflows or inflows seen mostly at larger radial distance ∼ 300 ckpc.For example, MRC0943-242 stays relatively flat (for both v 50 and W 80 ), while MRC0316-257 has a decrease in W 80 followed by an increase beyond the jet extent on the receding side.We reiterate that in this analysis we do not distinguish between (potential) different kinematics components by using v 50 and W 80 to quantify the overall velocity shift and line width.This may bring bias of the measured values.Additionally, the measured kinematics farther from the AGN are averaged from larger annulus (e.g.projected area of ∼ 4 × 10 4 pkpc 2 at ∼ 60 pkpc or ∼ 300 ckpc) which will bring another bias.We point out again that we use grey shade to mark the data > 5 ′′ from the center which has larger tile size.Nevertheless, we note that the detected W 80 of ∼ 10 3 km s −1 (and abrupt velocity shift) at large radii (∼ 300 ckpc) in some of the profiles could be caused by the fact that the detected Lyα emission is dominated by emission halos of nearby companions (e.g.Byrohl et al. 2021).

Morphology of the nebulae
The nebula morphology is related to the ionising sources, gas dynamics and galaxy environment (Byrohl et al. 2021;Costa et al. 2022;Nelson et al. 2016).Especially when the Lyα nebulae (in our sample) can probe the CGM gas beyond 100 pkpc.By visual inspection, we observe that the shape of our Lyα nebulae are asymmetric (e.g.Fig. 1).In this section, we study quantitatively the nebula morphology.We first focus on the whole nebula by introducing the morphology quantification measurements (ellipticity, nebula orientation and offset between nebula centroid and AGN position) in Sect.4.3.1 and compare with other samples in Sect.4.3.2.Then, in Sect.4.3.3,we study how the nebula asym-  Notes. (†) Intrinsic flux-weighted (e weight ), observed flux-weighted (e weight, obs ) and unweighted elliptical asymmetric measurements.The fluxweighted ellipticity is sensitive to the morphology of the high surface brightness part of the nebula.The unweighted ellipticity quantifies the morphology of the whole nebula.We note that for e → 0, the nebula is closer to round shape and vice versa. (*) Absolute difference between fluxweighted (unweighted) nebula position angle, θ weight (θ unweight ), and radio axis position angle, PA radio . (‡) Offset between the intrinsic flux-weighted centroid of Lyα nebula and AGN position. (§) Maximum linear extent spanned by the Lyα nebula.
metry changes with radial distance for individual targets.We also report the detected morphology correlations between different measurements (i.e.offset between AGN and nebula centroid position, nebula ellipticity and nebula linear size) in Sect.4.3.4.These shed light on how the central quasar and nearby companions can affect the observed nebula morphology.Finally, in Sect.4.3.5, we show the non-random oration of jet axis and its relation with the elongated direction of nebula which hints at the CGM gas distribution.

Morphology quantification measurements
To quantify the asymmetry, we introduce a set of morphology measurements.Arrigoni Battaia et al. (2019) used flux-weighted asymmetry measurements for the Lyα nebulae which is sensitive to the high surface brightness part.In other works (e.g.den Brok et al. 2020), an unweighted asymmetry measurement was adopted for better studying the extended structure of the nebulae (sensitive to the morphology of the whole nebula).To better characterise the morphology of our HzRGs and perform comparison with other samples, we analyse the asymmetry with both the flux-weighted and flux-unweighted methods.First, we follow Arrigoni Battaia et al. (2019) and calculate the flux-weighted asymmetry, α weight (see Arrigoni Battaia et al.  2021)).We mark the median e unweight for type-1s (0.69) and type-2s (0.80) in solid and dashed horizontal lines, separately.The e weight is sensitive to the morphology of the high surface brightness part of the nebula while the e unweight quantifies the morphology of the whole nebula.We note that for e → 0, the nebula is closer to round shape and vice versa.At the bottom right, we show the median uncertainty of the intrinsic L Lyα for our sample in logarithmic scale, 0.04.The ellipticity for our sample are higher compared to the other quasars for both high surface brightness part and whole nebula.There is no clear evidence that the nebula ellipticity correlates with luminosity.
2019, for the definition).This quantifies the asymmetry of the nebula in two perpendicular directions.Together with the asym-metry measurement, we also obtain the flux-weighted position angle θ weight , which we use as an indicator for the elongation direction of the nebula after converting it to the same reference system as the radio jet axis (i.e.angle measured east from north).The flux-weighted nebula centroid (centre of the nebula) is also computed.We note that the intrinsic flux and its uncertainty are used as weight to measure these three parameters and to calculate the corresponding uncertainties, respectively.We also derive the asymmetry measurement weighted by observed flux for comparison.Second, to compare with flux-unweighted asymmetry reported for other quasar samples (e.g.Borisova et al. 2016;den Brok et al. 2020), we calculate the α unweight following den Brok et al. (2020).In this context, we also derive θ unweight (fluxunweighted position angle) to examine the jet-nebula relation with respect to the entire nebula.Fig. 8 visualises the weighted (nebula centroids and θ weight ) and unweighted (θ unweight ) parameters on the eight nebulae.
In Lau et al. (2022), the authors compared morphology of different quasar samples.Following their comparison, we convert the aforementioned asymmetry measurement (for both fluxweighted and unweighted), α, to an intuitive elliptical asymmetry measurement (or ellipticity) e = √ 1 − α 2 .For e weight → 0, the nebula is closer to round shape and vice versa.Table 3 re 2019) samples, our HzRGs are consistent in asymmetry measurements and on the higher end of their distribution (type-1 median e weight ≈ 0.72, dashed horizontal in Fig. 9a).The e weight of radio-loud type-1s from (Arrigoni Battaia et al. 2019) have a median of 0.69 which is even lower than the value of all type-1 targets (Arrigoni Battaia et al. 2019;Cai et al. 2019).This indicates that the radio emission in type-1 does not disturb the gaseous nebula as in our HzRGs at least along the plane of the sky.This further suggests that orientation is a critical factor (Sect. 5.3).By comparing the intrinsic flux-weighted and observed flux-weighted elliptical asymmetry measurements, we find the e weight can vary significantly (e.g.MRC 0316-257 and 4C+04.11).For MRC0316-257, we can already identify its asymmetric morphology through visual checking (Fig. 1).There is also a large error bar associated with the intrinsic fluxweighted ellipticity.Hence, the morphology of MRC0316-257 is more towards the asymmetric end.The large difference between its intrinsic and observed e weight could be due to the absorption correction elevates the flux difference between the high and low S/N regions thus gives more weight to the central nebula.As for 4C+04.11,this could be due to the potential over-correction of the absorption in the low S/N regions given we use nine absorbers across the nebula (Sect.3.2, 5.1).However, we point again that the absorption is necessary to reconstruct the intrinsic flux given that the absorption features across the nebula were observed (e.g.Swinbank et al. 2015;Wang et al. 2021b).

Comparison of nebula asymmetry with other quasar samples
Our HzRGs have a median e unweight ≈ 0.70 which is lower than the measurement of type-2s and relatively consistent with the type-1s (median e unweight ≈ 0.80 and ≈ 0.69, respectively, Borisova et al. 2016;den Brok et al. 2020;Sanderson et al. 2021;Mackenzie et al. 2021).The ellipicity of the comparison type-2s is calculated from five sources which may not be representative.We note that there is a large scatter in the e unweight measured for our sample with four clustered around the type-2s and other four below the median e unweight of type-1s.This result could be biased by the implemented analysis methods (i.e.detection map and tessellation, Sect.3.1): (i) the construction of the detection map may smooth out the nebula asymmetry at larger radii (lower S/N); (ii) the tessellation results in lager bin sizes along the direction perpendicular to the radio jet (lower surface brightness and S/N, Sect.5.5) which shapes the nebula morphology to be more round.This brings more effects of the quantification of the entire nebula.Hence, the HzRGs nebulae will likely have larger e unweight (i.e. more asymmetric) than the quantified value.In spite of that, we can conclude that at least half of our sample, together with the type-2s, distribute at the most asymmetric end in terms of the whole nebula.The rest has the possibility to be skewed to higher e unweight .Our HzRGs have diverse properties and are not necessarily representative for the entire HzRGs population.Inspection for individual targets is required.
As for the L Lyα , our HzRGs host the most luminous Lyα nebula compared with other quasar types.This is due to: (i) sample selection3 ; (ii) absorption correction; (iii) quasar PSF subtraction of comparison samples (Sect 4.2.3).
We further compare our HzRGs to the ERQ from Lau et al. (2022).In Fig. 9, both the e weight (0.44) and e unweight (0.69) of L22 are lower than the measurements from HzRGs.Its L Lyα is also lower than our HzRGs by ∼ 2 orders of magnitude.This confirms that the nebula of this EQR is type-1 like but highly obscured (e.g.Lau et al. 2022).Since it is the only source of this type, we do not further discuss it.

Asymmetry radial profile
To further quantify the morphology of individual nebula as a function of distance from the AGN, we follow den Brok et al. ( 2020) and decompose our HzRGs intrinsic surface brightness using Fourier coefficients, a k (r) and b k (r): (3)  where SB Lyα (r, θ) is the surface brightness at given polar coordinate (r, θ).The coefficients a k (r) and b k (r) are defined as: The detailed description can be found in den Brok et al. (2020).
Then we combine a k (r) and b k (r): We measure the radial dependence of the asymmetry (i.e.how much it deviate from a circular shape) of the nebulae by using the ratio between kth and 0th modes, c k (r)/c 0 (r).In practice, only the first three modes are used since the higher mode coefficients are much smaller (den Brok et al. 2020).Fig. 10 shows the radial profiles of this asymmetry measurement (represented by (c 1 /c 0 ) 2 + (c 2 /c 0 ) 2 on the y-axis) for our HzRGs.
The larger the value, the more asymmetric the morphology shift from circular shape at a given radial distance.We also include type-2 measurements from den Brok et al. ( 2020 ) is also shown.Our HzRGs generally show an increase of the asymmetry as a function of radial distance (some have a smaller secondary peak at ≲ 50 ckpc) and a steep rise to > 1.5 (7 out of 8) at ∼ 250 ckpc.The most noticeable exception is MRC0316-257 which has a secondary peak (∼ 0.8) at ∼ 50 ckpc (∼ 14 pkpc).This flux-weighted measurement confirms the large difference we observed in the high surface brightness part of the directional profiles of MRC0316-257 in Fig. 6.As we already stated in Sect.4.2, the detected extent of our HzRGs (≳ 400 ckpc) are larger than the comparison quasars (∼ 300 ckpc).If we limit the comparison to the largest extent (∼ 250 ckpc) reached by the type-1s from (Borisova et al. 2016), we  3.We mark the obvious outliers in the two distributions in corresponding colors.This shows that most of our HzRGs have their Lyα nebula elongated along the jet axis.
find that our HzRGs are similar in radial asymmetry measurements.Specifically, both have a relatively flat profile to ∼ 200 ckpc which followed by a shallow rise to the value of ∼ 0.5.However, in Sect.4.3.2(Fig. 9) we find that the HzRGs ellipticity is larger than type-1s.Together with the radial asymmetry profile in this section, we can conclude that the asymmetry of the nebulae associated with HzRGs and type-1s are different due to structures at larger radial distance (> 250 ckpc) from AGN.This is also suggested by den Brok et al. (2020).Although we caution the large W 80 at large distance in Sect 4.2.4 and mark the region with r > 5 arcsec in grey in Fig. 7, we can still find that most targets have at least W 80 ∼ 10 3 km s −1 .This may indicate that these structures are likely disturbed and not 'quiescent' gas in the Cosmic Web (Sect.5).From the comparison with type-2s (den Brok et al. 2020;Sanderson et al. 2021, which also have a steep rise), we find that the radial asymmetry of our HzRGs rises at larger radial distance (≳ 250 ckpc) and reaches higher asymmetry value (> 1.5 compared with ∼ 1.4 of type-2s).At least one of the radial asymmetry measurement (Cdfs 15, with the furthest extent r AGN ∼ 300 ckpc) from den Brok et al. ( 2020) shows an indication of continuous rising up to the detection limit (1 hour integration with MUSE); we cannot rule out the possibility that this target may show a similar trend as the HzRGs with deeper observations.

Morphological correlations
In Fig. 11, we compare the e weight (ellipticity that is sensitive to high surface brightness nebula) and d max (maximum nebula extent) against d AGN−neb (offset between AGN and nebula flux-weighted centroid).For a consistent comparison for the unweighted measurements, we should use the unweighted offset (centroid) and the unweighted ellipticity.We note that the calculation of unweighted ellipticities, e unweight , assumes that the centre of nebula corresponds to the AGN position (den Brok et al. 2020).It would involve large uncertainties if we calculated the flux-unweighted nebula centroids (i.e. the geometric center of the nebula) which are entirely depended on the spaxel distribution from our detection maps (Sect.3.1.1).Hence, we report only the correlation between flux-weighted ellipticity and offset, and use it as a proxy of the nebula even though they are more sensitive to the high surface brightness part.The (r-value, p-value) of Spearman's correlation coefficients are (0.89, 0.007) and (0.88, 0.004) for the e weight vs. d AGN−neb and d max vs. d AGN−neb relations.We note that due to the large uncertainties of e weight for MRC0316-257, we exclude it from the correlation measurement.From the r-values and p-values, we can conclude that e weight and d max are both correlated with d AGN−neb .This suggests that the more asymmetric and more extended nebulae have larger offsets between AGN position and nebular centroid.We also include the Arrigoni Battaia et al. (2019) 2019) sample and the radio-loud targets are (0.4,0.003) and (0.4, 0.1), respectively, suggesting a moderate correlation in their whole sample but not in their radio-loud target.This indicates that the radio-loud type-1s are different from our radio-loud type-2s HzRGs.We will revisit these correlations in Sect. 5.

Jet-nebula position angle distribution
We present the distribution of the |θ − PA radio | (Sect.4.3.1,Table 3) in Fig. 13.Both the flux-weighted and unweighted measurements are shown which are sensitive to high-surface brightness parts and the entire nebulae, respectively.Fig. 13 shows that most HzRGs have |θ − PA radio | < 30 • .A similar result was reported by Heckman et al. (1991b).This observation is consistent with the scenario proposed by Heckman et al. (1991b,a) that the compression of the gas by the radio jet results in brighter emission along the jet.We will discuss the indications behind the alignments further in 5.5.The exception is 4C+04.11which has both large |θ weight −PA radio | (81.4 • ) and |θ unweight −PA radio | (45.3 • ), indicating different conditions that in the other targets (such as inflows, Sect.5.6).In TN J0121+1320 we find a flux-unweighted angle difference of 72.1 • .Given its round Lyα nebula morphology and compact radio emission, there is a large uncertainty in angle difference measurement and we do not discuss this source separately.If we assume the observed angle difference is equally distributed from 0 − 90 • , the probability for us to find 7 (or 6) targets in a sample of 8 with |θ − PA radio | < 30 • is only 0.2% (1.7%) using bimodal distribution.

Absorption correction and radiative transfer
We stress that the analysis of our sampled Lyα nebulae in this paper is under the idealised assumptions stated in Sect.3.2.1.Specifically, we interpret the trough features in the Lyα spectra as absorption features by the neutral hydrogen gas surrounding the radio galaxy.In this section, we discuss the limitations and possible caveats of this treatment for reconstructing the intrinsic Lyα flux along with the possible effects brought by radiative transfer.
Under our assumptions, the absorbing gas is located in between the observer and the last scattering surface of Lyα photons along the line of sight.The intrinsic Lyα flux reconstructed under our treatment is thus assumed to be the one after radiative transfer processes have shaped the Lyα nebula, and is approximated by a Gaussian profile.We also assume the absorbed Lyα photons in the absorbers to be 'lost' rather than continuing their resonant scattering path within the nebula; this may happen when photons are absorbed by dust and re-emitted as infrared thermal emission, or preferentially scattered away from the line of sight due to a particular geometry (see Sect. 3.2.1).The latter may occur because photons originating from the backside of the galaxy have a higher chance to be absorbed by dust when transiting the host galaxy (Liu et al. 2013).These absorbers can be interpreted as intervening low column density gaseous shells (e.g.van Ojik et al. 1997;Swinbank et al. 2015;Kolwa et al. 2019).The fact that most of the absorbers are located in the blue wing of the Lyα profile as well as the fact the trough of the main absorber is often seen across several tens of kpc scales supports this idea (see Fig. 1a, and TNJ0205+2242 in Fig. 2, see also TNJ1338-1942  The situation may be different when we encounter embedded absorbers which are supposed to be closer to the source of Lyα photons (i.e.central AGN in our case).This leaves the trough (or 'main absorber') at around the systemic redshift of the AGN as predicted by radiative transfer studies (e.g.Verhamme et al. 2006).Even though the scales probed in those simulations are different (CGM scale versus sub-kpc), this might be the case in some of our targets, for example MRC0943-242 and TNJ0121+1320 (Fig. 3 and C.1).Therefore, one consequence of our absorption treatment (with 'lost' photons) is that we may double-count Lyα photons that are resonantly scattered to the wing and/or other directions by redundantly adding more when correcting for absorption.We implemented a simple test for checking the double-counting effect, which takes advantage of the IFU nature where slit or spaxel loss effects are compensated by photons resonantly scattered in the neighbouring spaxels.Specifically, we summed all individual spectra into one and reconstruct the intrinsic flux for this single spectrum ( f Lyα,full FOV ); this value removes the IFU information, but should be a good measure of the total value emanating from the nebula in the observer's direction.This is then compared to the sum of individually constructed intrinsic flux in each tile (Fig. 1b, Σ f Lyα,tile fit ).Fig. 12 shows the result, where a value of 1.0 is expected if our treatment is fully accurate.Instead, we observe a median value of 0.71 for our sample, indicating that there may be double-counting of photons.However, several effects may also contribute to this.First, our assumption of a Gaussian shape for the intrinsic profile may not be accurate, as prior radiative transfer effects may have created troughs and boosted the line wings (Verhamme et al. 2006).A future paper presenting high spectral resolution UVES observations of our sample will help to better model the profiles (Ritter et al. in prep.).Second, the assumption that all absorbers extend across the entire nebula may also be oversimplified.While many absorbers are indeed seen on 10s (or 100s) of kpc scales for our targets (Fig. 2, 3 and 4 and also Appendix C), there may be exceptions.
The Lyα nebula properties of the type-1 (and radio-quiet type-2) sources in the comparison samples are derived without correcting for absorption.They are still good comparison samples given that not many absorption features are detected in them (e.g.Arrigoni Battaia et al. 2019).This fact alone is already an indication that there may be an environmental differences between the nebulae of HzRGs and other quasars.Alternatively, this difference may related to intrinsic differences between different AGN types.
Overall, it is likely that both the absorption and radiative transfer effect are working together shaping the the Lyα profile.Our analysis assumes that the absorption correction is the dominant factor.

Lyα nebula and AGN unification model
The unification model of AGN (e.g.Antonucci 1993) proposed that the fundamental difference between type-1 and type-2 is due to the different orientation of the obscuring dusty torus (ionisation cone).Specifically, the ionisation cone of type-1s is more aligned with the observed line of sight than in type-2s.
Within this unification model, we assume that the power of AGN is on a similar level for type-1s and type-2s and that their gaseous nebulae therefore have similar distributions and masses.If we further assume that the ionising photons are primarily produced by the AGN, then the nebulae should have similar morphologies and luminosities.In this picture, the Lyα nebulae are elongated along the direction of the ionisation cone and have a 'rugby-ball' shape.For type-1s, the ionisation cone would be pointing towards the observer resulting in a rounder nebula morphology.For type-2s, the ionisation cone would aligned along the plane of the sky resulting in the observed elliptical morphology.Evidence for such a scenario was indeed reported in He et al. ( 2018) using ionized gas nebulae but for local AGN and on small scales (sub-kpc to kpc).
Using both the e weight (sensitive to the high surface brightness nebula) and e unweight (whole nebula) quantifying the ellipticity of the nebula, we find that the HzRGs (and other radio quiet type-2s) are more asymmetric than the type-1s in Sect.4.3.1.This observation agrees with the AGN orientation unification model.The orientation of type-1s probably still vary in a range which causes the large range of the ellipticities.The AGN luminosity and dark matter halo mass (gas distribution) can also play an important role in shaping the observed morphology (Fig. 9).
With the jet axis indicating the direction of the ionisation cone (Drouart et al. 2012) and the alignment between the jet axis and the Lyα nebula (at least in the high surface brightness part where the photons of AGN are presumably dominating the ionisation, Sect.4.3.5 and Fig. 13), we argue that the AGN orientation between type-1s and type-2s (HzRGs) can explain the observed morphological difference.The evidence for this claim comes mostly from the observed ellipticity distribution (Fig. 9): the e weight (median 0.69, Sect.4.3.2) for radio-loud type-1s (Arrigoni Battaia et al. 2019) are lower than our HzRGs (median 0.78), which is consistent with the jet (ionisation cone) pointing more towards observers in type-1s.By checking the available radio maps of these radio-loud type-1s (e.g.Fig. B1 in Arrigoni Battaia et al. 2019), we indeed find that they have compact radio morphology (i.e.not the two-side jet like HzRGs) suggesting that the jets are aligned more perpendicular to the plane of the sky.We note that TNJ0121+1320 also has a compact radio morphology (i.e.not having two-side jet) and the lowest e weight (0.66) and e unweight (0.48) in our sample (more consistent with type-1 results, see Fig. 1 and 9) which again fits the unification model.
A similar explanation was also suggested by den Brok et al.  2020) found more symmetrical morphologies for their type-2s at small radial distances < 30 pkpc that are more similar to type-1s.den Brok et al. (2020) suggested that additional ionising sources other than the AGN (e.g.star forming processes) may contribute to this observation and smear out the differences between type-1 and type-2 at such small radii.As for our HzRGs, the jet-gas interaction at ≲ 20 − 30 pkpc (∼ 100 ckpc) could be a reason for the observed high e weight shown in Fig. 9a (compared to type-1s) and the reason of the small rise (showing higher asymmetry compared to type-2s) in the radial asymmetry measurement at ∼ 60 ckpc in Fig. 10.

The role of Lyα resonant scattering
As presented in Sect.4.1 and 4.3, Lyα nebulae around our HzRGs are extended in size (≳ 150 pkpc) and are asymmetric in shape.Interestingly, there is a correlation between the maximum nebula extent d max (e weight : ellipticity measurement that is sensitive to high surface brightness nebula) and the offset between AGN position and nebulae's flux-weighted centroid d AGN−neb (Sect. 4.3.4).This correlation may also reflect our findings regarding the surface brightness and kinematic radial profiles in Sect.4.2.1 and 4.2.2, such as the exponential shape of the surface brightness in the inner nebula and high W 80 value.The resonant nature of Lyα photons may be related to this observation.
Villar- Martin et al. (1996) first proposed the idea of resonant scattering being related to the observed extended nebula emission around HzRGs.In simulations, Costa et al. (2022) found that the scattering of Lyα photons, regardless of the powering source, could result in an offset between the luminosity centroid of the nebula and the quasar position (d AGN−neb ).This offset can vary depending on the line of sight and can reach ∼ 15 pkpc.Such offsets are consistent with what we find (Table 3).The authors also found that scattering can shape the nebula into a more asymmetric morphology (e weight → 1).This depends on the gas distribution of neutral hydrogen and observed line of sight as described in Costa et al. (2022).Given the common case that gas density is the highest on galaxy scales, and that we target type-2 AGN, we expect the Lyα photons to scatter out to larger projected distances rather than travelling directly towards the observer.Such a scenario may explain our observed correlation in Fig. 11 between e weight and d AGN−neb .Specifically, when the scattering is more efficient, we may observe the nebula to be more asymmetric and more extended, at least in the high surface brightness regime.
The inner part (luminous) of the surface brightness radial profiles have an exponential shape (Sect.4.2.1).This gradual change in surface brightness and profile slope at smaller radii is suggested to be due to scattering (Costa et al. 2022).The high W 80 values out to radius of ∼ 50 pkpc (or ∼ 230 ckpc, Sect.4.2.4) could also be related with efficient scattering (Fig. B1 in Costa et al. 2022).The velocity shift of the nebulae can also be complicated due to scattering at the similar radial distance (Fig. B1 in Costa et al. 2022) which is the case of our v 50 (Fig. 1c,Sect. 4.1).As shown by the stellar radial profiles (Fig. 6), the impact of seeing cannot be fully responsible for the exponential shape at smaller radii.For the observed kinematics, the jet (which is not included in the simulations) may also play a role here.We overlay the jet hot spot distances on the radial profiles (Fig. 6 and 7).Qualitatively, the v 50 and W 80 show different behaviours within and outside the extent of the jet hot spots at least for some targets (e.g.MRC0316-257).We note that the distances marked by the vertical lines are measured from the brightest radio emission po-sitions, i.e. the full extent of jet is even larger (Appendix D).The decrease of W 80 at large radii (≳ 100 pkpc) observed in some targets (e.g.TN J1338-1942) is, beyond the scope of the Costa et al. (2022) simulations but could be related to the (unvirialised) cosmic web.
In the framework where scattering significantly impacting the observed Lyα properties, Costa et al. (2022) furthermore predicts that Lyα nebulae of edge-on AGN should have lower surface brightness, larger d AGN−neb , more asymmetric shape and flatter surface brightness profiles in the centre.Most of the predictions agree with our observations except for lower surface brightnesses, which could be due to our selection criteria and/or jet impact.We note that the analysis of Costa et al. (2022) has been done without correcting for absorption which is reason of the discrepancy.
The orientation (Sect.5.2) may be the reason for the moderate d max − d AGN−neb relation seen in the sample of (Arrigoni Battaia et al. 2019) (i.e. the orientation spans a large range within the sample).The larger the viewing angle (i.e. the more edge-on we are observing the AGN, assuming unification model, Antonucci 1993), the more extended the nebula is expected to be, because both sides of the nebula will be observed.The nebula is expected to become more asymmetric and have larger d AGN−neb as Costa et al. (2022) predicted.The study of type-2 quasars (e.g.den Brok et al. 2020;Sanderson et al. 2021) also found that the difference in nebula morphology when comparing to type-1s is related to AGN orientation.Interestingly, for our HzRGs, we do find a relatively strong correlation between d max and d AGN−neb despite the expectation that most HzRG are observed under similar, large viewing angles, (i.e.close to edge-on).All of our targets have clear two-sided radio lobe morphology (except TN J0121+1320) and none of them shows clear signs of broad-line emission or significant contamination by a bright point-like source in the centre (see also Drouart et al. 2012, and Sect. 2.1.3).
Resonant scattering has the potential to shape the observed inner parts of the radial profiles (surface brightness and kinematics) and the morphology correlations.Byrohl et al. (2021) studied Lyα emission halos and their relation with environmental factors using TNG50 simulation.They found a flattening in the Lyα halo surface brightness radial profile at large radial distances (≳ 30 pkpc).The authors attributed this to the contribution of scattered photons from nearby halos (companions).In Sect. 4.2 (and 4.2.3),we show that the profiles of three of our targets (MRC0943-242, MRC0316-257 and TN J1338-1942) also have an obvious flattening at large radii (Fig. 6).Coincidentally, these targets are known to live in dense environments (on Mpc scale, e.g.Venemans et al. 2007).All of this indicates that our nebulae are 'contaminated' by Lyα halos associated with nearby companions (e.g.Gullberg et al. 2016).In addition, the observed surface brightness radial profile of MRC0316-257 (Fig. 5) at large radii, shows a decline followed of a rise, and is a factor of five brighter than the 2σ surface brightness limit (16 for the intrinsic) which indicate an extra source of Lyα photons.The filamentary Lyα emitting comic structures are observed on Mpc scales (e.g.Umehata et al. 2019;Bacon et al. 2021).Bacon et al. (2021) discussed the possibility that the diffuse Ly-alpha emission is powered by low mass galaxies not directly detectable in the deep (140 hours) MUSE observation.

Large scale environment: nearby Lyα emission halos
These results provide additional evidence that the detected nebula connects with the emission halo of companions.If this is the case, we should revisit the d max − d AGN−neb in Sect. 4.3.4.When a quasar resides in a dense environment, its apparent size will likely be impacted by neighboring Lyα nebulae.This 'contamination' from nearby companions will likely be more related to the large-scale structure and independent of the orientation of the central (quasar) halo.In other words, the hidden parameters behind the d max − d AGN−neb (e weight − d AGN−neb , Sect. 4.3.4,Fig. 11) relations may be related to the distribution of the companion emission halos.Contamination from neighboring halos will result in the centroid of the nebula being offset from the AGN position and cause a more asymmetric nebula morphology.The nearby companions can be the disturbing sources resulting in the ∼ 1000 km, s −1 line width seen at ∼ 100 pkpc (W 80 , Fig. 7).
The radial asymmetry profiles of our targets have a larger value than type-1s and type-2s (Fig. 10) at larger radii, which is unlikely to be entirely due to the AGN orientation (Sect. 4.3.3).This can now be explained by the impact and contamination of nearby companions at 100s of pkcp.We point out that the type-2 from Sanderson et al. (2021) has a projected size of 173 pkpc, a nebula centre offsets from the AGN of ∼ 17 pkpc and high ellipticity (0.8).This source is also found to have nearby companions (∼ 60 pkpc).The difference of the radial asymmetry profiles (Fig. 10) and surface brightness profiles between this source and our HzRGs may be related to the jet (or large scale gas distribution, Sect.5.5).
The stellar masses of the galaxies studied in Byrohl et al. (2021) are in the range of 8.5 < log 10 (M ⋆ /M ⊙ ) < 10 which is approximately 1 − 2 orders of magnitude lower than the stellar masses of our HzRGs, implying lower dark matter halo masses, as well.Such lower masses may explain the difference in the level of Lyα surface brightness where the flattening of the profiles is observed: in the simulations by Byrohl et al. (2021) the flattening happens at a surface brightness level of ∼ 10 −20 erg s −1 cm −2 arcsec −2 (e.g.their Fig.7) while we observe it to happen at a level of ∼ 10 −18 erg s −1 cm −2 arcsec −2 ).The dark matter halo mass difference can also explain the different radial distance at which the companion emission dominates (∼ 30 pkpc in Byrohl et al. 2021, versus ∼ 40 pkpc, Table D.3).We note that we use the r b , the radius at which surface brightness radial profile changes from exponential to power law (Sect.4.2.3), as the distance where nearby halos start to significantly impact the surface brightness.
The ≳ 100 pkpc extents of the Lyα and the dense environments of HzRGs (Wylezalek et al. 2013) suggest that nearby halos may contribute to our observations (Byrohl et al. 2021).However, this also suggests that our CGM study is 'contaminated' by sources other than the radio galaxy.For this paper, we exclude obvious emission of clearly detected companions (e.g.'Arrow galaxy ' Vernet et al. 2017).A systematic census of the companion galaxies in the MUSE HzRG fields will be reported in a future paper.

Gas distribution on large scales
So far, we have focused our discussion on morphological measurements sensitive to the high-surface brightness part of the Lyα nebulae and plausible explanation of environmental effect on ∼ 100 pkpc scales.It is likely that the jet shapes the nebula in the vicinity of the quasar to align with the jet axis (skewed distribution of |θ weight − PA radio | towards values < 30 • , Fig. 13) through interaction with the gaseous medium.Beyond the extent of the jets, we use |θ unweight − PA radio | (which is equally sensitive to the whole nebula) to inspect the relation between the jet axis and halo asymmetry.
While it is unlikely that the jet plays a major role shaping the morphology of the large-scale halo, Eales (1992) suggested that the observed direction of the radio jet is the result of an inhomogeneous gas density.When a jet is launched along the direction of higher gas density (e.g.n H ∼ 10 −2 cm −3 ), the jet 'working surface' (hot spots) will leading to brighter radio emission.Such a scenario would introduce a bias in observing jets preferably along directions of higher densities in flux-limited samples and we might expect a correlation between the morphology of a large scale halo and radio jet axis.Our observations presented in Figures 11 and 13 is in agreement with such a scenario.
The detection of molecular gas in HzRGs along the jet axis but beyond hot spots provides further evidence (< 20 • , which suggests the jet runs into filament of cold gas Emonts et al. 2014).If this direction indeed traces a filament of the cosmic web (e.g.West 1991;Umehata et al. 2019;Bacon et al. 2021) with a higher density of companion galaxies, then this is also in agreement with the discussion presented in Sect.5.4.Humphrey et al. (2007) found evidence of infalling CGM gas on to the HzRGs which bridges the radio galaxy and to the large-scale (IGM) gas.If the CGM gas is being accreted onto the central HzRGs through these higher density distributions, it would reflect the scenario proposed by Humphrey et al. (2007).

Inflow from the cosmic web?
The relatively large |θ − PA radio | (81.4 • and 45.3 • for weighted and un-weighted, respectively) detected in 4C+04.11 is inconsistent with other targets (see Figure 13).This suggests that the nebula of 4C+04.11 is elongated in the East-West direction on both small (< 20 pkpc, scope of jet) and large scales (≳ 20 pkpc,Sect. 4.3.1 and 5.4).The tile with the largest distance from AGN position (# 74,Fig. A.9) has a filamentary like structure.There is no known evidence that 4C+04.11resides at the center of a dense cluster-like environment (e.g.Kikuta et al. 2017) which makes it unlikely that the asymmetry is caused by contamination of nearby halos as discussed in Sect.5.4.
High velocity shocks (≳ 100 km s −1 , Allen et al. 2008) can heat the gas which will then cool by radiating UV photons.Shocks could be caused by inflows and power the observed Lyα emission.We estimated the dark matter halo of 4C+04.11 is on the order of M DM ∼ 10 13 M ⊙ from it stellar mass, M ⋆ ∼ 10 11 M ⊙ with an empirical M ⋆ − M DM relation (Wang et al. 2021b).The surface brightness measured in the farthest tile #74 is ∼ 1.2 × 10 −17 erg s −1 cm −2 arcsec −2 which is similar to (or slightly higher than) the simulated value from Rosdahl & Blaizot (2012).We converted our measurement to z = 3 which is redshift reported in the simulation.Our measurement and the simulation are consistent given that dark matter halo of 4C+04.11 is likely more massive than the one in the simulation (M DM ∼ 10 12 M ⊙ ).
Given the estimated dark matter halo mass, the virial radius and virial velocity can then be calculated as r vir ∼ 100 kpc and v vir ≃ 580 km s −1 .From the centre of the tile #74, we can derive its distance to the central AGN as ∼ 60 kpc.We note that this is the projected distance averaged for spaxels in the tile.The actual distance of the filament could be farther.Even though it is closer, Nelson et al. (2016) simulated gas accretion into 10 12 M ⊙ halos at z = 2 and found that the accretion flow structure can remain filamentary within r vir (∼ 0.5r vir ).The projected v 50 (velocity offset) of tile #74 is −536 km s −1 which is consistent with v vir .Thus, our measurements for the projected distance of the tile and its velocity offset are consistent with a scenario where the Lyα emission in tile #74 may be tracing shock driven by inflows.If the distance we observe is indeed < r vir , then the post-virial ac-cretion could be more complicated with multi gas phase components and fragmentation involved (e.g.Cornuault et al. 2018).The broad line width of 4C+04.11 at large radii (Fig. 7) may indicate the disturbed nature of the presumed accretion flow.
The discussion in this section based on the |θ − PA radio | and morphology of bin #74 of 4C+04.11.While the angle difference is a clear outlier in the sample, the tile shape may depend on the implemented method for nebula extent (Section 3.1.1).However, we checked the Lyα narrow band image (6690 Å to 6707 Å) directly collapsed from the data cube and confirm the indication of emission from this region.As for the nearby potential emission (north to bin #74, Fig. A.9), we confirm by spectral extraction that there is no S/N> 3 detection.Even if the Lyα emission is detected in this additional region, it would still be consistent with the accretion scenario.Complex filamentary structures are seen in simulations (see Rosdahl & Blaizot 2012, and their Fig.7) for accretion in higher mass dark matter halo.Nevertheless, we caution that this analysis only considers one possibility.As the data is near the detection limit, deeper observation are needed before more conclusive results can be derived.

Conclusions
In this paper, we study the intrinsic Lyα nebulae around a sample of eight high-redshift radio galaxies using optical IFU observations obtained with MUSE.We link our observations to results from the literature for other quasars (mainly type-1 quasars) at similar redshifts.
We have developed a new method to measure the maximum extent of the nebulae with improved sensitivity to low surface brightness emission.We also have developed a new method to tessellate the Lyα maps (Sect.3.1).We have detected the Lyα emission at scales ≳ 100 pkpc from the central AGN, down to an observed surface brightness of ∼ 2 − 20 × 10 −19 erg s −1 cm −2 arcsec −2 .
We summarise our results as follows: The Lyα emission line profiles of all sources are affected by multiple and deep absorption troughs out to the edge of the nebulae.We have corrected for this H i absorption (Sect.3.2) and constructed maps of the intrinsic Lyα (Sect.4.1) emission and also measured the kinematic properties of the Lyα emission spatially resolved.
We first investigated radial dependencies of the surface brightness, velocity shift and velocity width of our sample.We found that circularly averaged profiles of the intrinsic (absorption-corrected) Lyα surface brightness are brighter and more extended than type-1 quasar samples (Sect.4.2.1).We did not find major differences when we investigate the radial profiles along the approaching and receding direction of the radio jets, respectively (Sect.4.2.2).The surface brightness radial profiles can generally be described by an exponential drop-off for the inner high surface brightness part and a power law for the more extended part (Sect.4.2.3).We did not find evidence of ordered gas bulk motion from the v 50 radial profile (Sect.4.2.4).For four targets, the v 50 profiles at radii within the jet hot-spots range show a similar gradient as the jet-driven outflow (Fig. 7).The W 80 profiles show relatively large values (≳ 1500 km s −1 , Fig. 7, Sect.4.2.4) at all radii and three targets show a decrease beyond the distance of jet hot spots is indicative of jet-gas interactions.
We quantitatively studied the morphology of the HzRG nebulae (for both observed and intrinsic emission) and compared our results to other quasar samples (uncorrected for absorption, Sect.4.3).We found that our nebulae are, in general, more asymmetric than nebulae of type-1 quasars and are more similar to type-2 quasars (Sect.4.3.2),although, our sampled sources differ in their measure of asymmetry as a function of radius (Sect. 4.3.3).Furthermore, we found that the more asymmetric and larger the nebulae are, the greater the offset between the centroid of the nebulae and AGN positions (d max −d AGN−neb and e weight −d AGN−neb in Sect. 4.3.4,Fig. 11).
Lyα is a complicated emission line that can be heavily affected by absorption and resonant scattering which, as we demonstrated, needs to be accounted for.Assuming type-1 and type-2 quasars have similar intrinsic shapes of their nebulae, the difference of the nebulae asymmetry morphology between our sample and other quasars can be explained using the AGN unification model where the orientation of the ionisation cone is the fundamental parameter (Sect.5.2).Resonant scattering of the Lyα emission can result in the observed d max − d AGN−neb and e weight − d AGN−neb relation in our sample and partially explain the shape of the radial profile and the kinematic profiles (Sect.5.3).We also found evidence in our HzRGs that the extended nebulae are affected by Lyα emissions from nearby companions (Sect.5.4) and that CGM gas has higher density distribution along a specific direction which is coincident with the direction of the radio jet (Sect.5.5).There is also a possibility that, in our observations, we are capturing inflows from cosmic web (Sect.5.6).
In this paper we have shown that measurements of Lyα nebulae around high-redshift AGN can act as a probe of the environment of AGN and their host galaxies.We have found fundamental differences between the extended ionised nebulae of different types of quasars at Cosmic Noon and beyond.Due to its resonant nature, it is a challenge to use Lyα as a tracer of gas kinematics especially out to the ∼ 100 pkpc scale.Nevertheless, Lyα line observations offer some insight into the state of CGM gas at a time when it is simultaneously being impacted by more than one physical mechanism (e.g.quasar outflow, jet-gas interaction and cosmic inflow).These kind of observations will be particularly useful for future simulation works.The MUSE data only reveals the tip of the iceberg.In upcoming papers, the H i absorbers will be reported together with the high-resolution spectroscopy data.In addition, all the UV emission lines covered by MUSE will be studied in detail and this will provide crucial information on the ionisation, metallicity and AGN feedback in the quasar nebulae.Furthermore, scheduled JWST observations will be available for four HzRGs in our sample and this will provide unprecedented details of the AGN and host galaxy connection on scales of ≲ 1 kpc (∼ 0.05 ′′ ).The basic idea for detection map is to use the 3D information of the MUSE data cube to determine the maximum spatial extent of the Lyα emission.Smoothing technique is the key process in this procedure to guarantee the faint structures are captured.To assist with conveying the procedure, the process is shown schematically in Fig. A.1.In this procedure, we can select the faint end of the emission nebula to the surface brightness detection limit.
Before the generation of detection map, we do a step of continuum subtraction for each spectra spatially centred around the AGN and spectrally around their Lyα wavelength range (−5000 − 5000 km s −1 or −6000 − 6000 km s −1 for 4C+03.24which has broader line width) with the emission line masked.If there are sky-lines, emission or absorption features from foreground targets in the unmasked range, we do further masking.Continuum from the host galaxy, central AGN and foreground objects need to be subtracted to minimise their contamination.In this way, we exclusively focus on the line-emission nebula.The choice of flat or linear continuum do not have significant impact on the Lyα fitting (Wang et al. 2021b).Hence to better account for the slope from foreground stars, we subtract a firstorder polynomial from the cube.
The 3D adaptive smoothing follows the Martin et al. ( 2014) procedure with adaptation (e.g.Vernet et al. 2017).As described in Sect.3.1.1,it smooths each of the image planes of the continuum-subtracted cube adaptively with a Gaussian kernel over a wavelength range of ∼ 15 Å around the Lyα emission.For the image slice at each wavelength, the algorithm first takes average of adjacent n λ number of slices around this image slice along the wavelength dimension (Fig. A.1a).Then it adaptively smooths the averaged image with a Gaussian kernel starting from the smallest kernel size, σ 0 = 3/2.355pix, until the maximum, σ max , is reached (Fig. A.1b).Specifically at each step, the algorithm smooths the spaxels that are not passing the pre-set S/N threshold, T S/N .In the end, the algorithm set the voxels, after being maximally smoothed, that not pass the T S/N to 0 as no-signal containing (i.e. the voxels that potentially contain Lyα signals have positive value).In this way, we preliminarily detect where we have Lyα in the cube (Fig. A.1c).
To generate a detection map working on the spatial dimensions and guarantee line fitting, we perform a further check along the wavelength dimension for the smoothed cube.Specifically, if the smoothed spectrum at one spatial location (x, y) has n c numbers of consecutive wavelength pixels preserved (i.e. have pos-itive values), then we flag this spatial location (x, y) as signalcontained (see Fig. A.1d as an example).The others that do not pass this check are discarded.In this way, we can construct the detection map (Fig. A.1e).
The question left is to find a combination of the four parameters (T S/N , n λ , σ max and n c ) which returns an optimised detection map.For this, we generate a series of detection maps for each targets with different combination of the four parameters.Then we choose the one that optimises the spatial selection (i.e.captures large scale extent while avoid island-like structures).The best combination of the four parameters for each targets are presented in Table A.1.We check the detection maps constructed using different sets of parameters around our optimal combination (Table A.1) and find that the bulk (∼ 80%) of the selected spaxels remains the same.If using a stricter set of parameters (e.g. higher T S/N and/or n λ ), we will miss the low S/N part of the nebula by comparing with previous individual target studies (Swinbank et al. 2015;Vernet et al. 2017).If using a more relaxed set, the peripheral regions, mostly disconnected to the bulk of the nebula, with pure noise (after checking spectrum) will likely be selected.Thus, we are assured that the constructed detection map represented the observation to the detection limit and the method is robust against objectiveness.After constructing the map with the best parameter combination, we perform a final manual selection to exclude island-like regions (≳ 1 arcsec 2 in size) which are detached from the the main nebulae (> 1") and could be due to noise or companion.A further check of the spectra extracted from these regions is also conducted.Around 50% of the island regions only contain noise.The ones with potential Lyα signal detected are presented in Appendix A.3 (Fig. A.3 and A.4 ).Since these are detached from the extended nebula around the quasar, we do not include them in this analysis but point out they may trace companion emissions.We note that we refer to the smoothed cube obtained in this way (using T S /N , n λ and σ max in Table A ing a twofold Voronoi binning by using a S/N cut for doing separate binning for high and low S/N parts (M.Cappellari private communication).Solution (i) will return wedge shape tiles which is a known result (Cappellari & Copin 2003).It is, however, not desirable for our case since it smears out the radial structures which is one of the key interest of this work.Therefore, we adapt solution (ii); the twofold Voronoi binning for our sample.Firstly, we apply the detection map (Sect.3.1.1and Appendix A.1) to the continuum-subtracted un-smoothed cube to preserve only the Lyα signal detected spaxels for tessellation.To avoid complication and keep consistency for the whole sample, we perform the tessellation on a single narrow band image for each target.Since the Lyα spectra may have different observed peaks at different spatial locations, we choose the wavelength range, over which the narrow band image for tessellation will be produced (averaged by the number of wavelength pixels), to enclose as much of line emissions as possible and avoid adding pure noise.This range is ∼ 10Å wide.We note that for 4C+19.71,we select two wavelength ranges at both the blue and red sides of 5577 Å(sky-line), which is at roughly the systemic redshift of Lyα.In this way, we construct the S/N map from the narrow band image.
Secondly, a Gaussian smooth is performed to the S/N map to minimise the impact of randomly located spaxels dominated by noise (i.e.further avoid detached region problem).We use a Gaussian kernel with σ s = 3 in pixel unit.Then, we apply a S/N cut of 3 to the S/N map to select the high S/N regions for the first Voronoi binning.
Thirdly, for the selected S/N > 3 part, we reassign spaxels with S/N>6 to S/N=6 and perform the Voronoi binning with a target (S/N) vorbin,1st = 15 (12 for 4C+03.24,see later).The reason for the S/N reassignment is to control the minimum size of the tiles to avoid single spaxel tile and account for the seeing effect.Because the Lyα nebulae studied here have high S/N gradient from center to outskirt, performing Voronoi binning with a high S/N threshold (avoid single spaxel bin in the high S/N region) will result in tiles with large size for the low S/N part which is an overkill for the fitting and smears out resolution information.After reassigning S/N>6 spaxel and using (S/N) vorbin,1st =15, the minimum number of spaxels for one tile is 6 (= 3×2 pix 2 or 0.6× 0.4 arcsec 2 ) which is roughly half of the seeing disk.This is a compromise for being consistency for the whole sample and also physically connects the neighbour tiles which is useful for the implemented fitting procedure (Sec.3).As for 4C+03.24which was observed in the AO mode, we lower its (S/N) vorbin,1st =12 to have a minimum number of spaxels for one tile is 4 (= 2 × 2 pix 2 or 0.4 × 0.4 arcsec 2 ).stay with the 1-Gaussian fit for all tiles due to the complicated spatial variance of spectral shapes.We note that this selection is not crucial in this work since we do not distinguish and separately interpret different velocity components in the analysis.The purpose for this step is to consistently ensure that the nonsingle-Gaussian line shape is considered without missing flux.
Fifthly, we fit the spatial spectra with 2-Gaussian models for the tiles determined in the last step and 1-Gaussian for the others.The results from the previous steps are passed to the next as initial guesses.To keep consistence, we choose to use the same number of absorbers for the whole map due to the difficulties in determining where one absorber disappears.This is a good assumption to first order given that we observe most absorbers on  except for 4C+04.115 .For spectra extracted from inner tiles (≲ 10 ′′ , with high S/N), we constrain their column densities to vary within a ∼ 2 dex range from the initial input.We ease the column density constrain for the absorbers from outer tiles with distance to AGN ≳ 10 ′′ such that they can be given low column density (∼ 10 13 cm −2 ).In this way, the absorbers at low S/N regions can have negligible impact (< 0.1 %) on the reconstructed flux.We use the same number of absorbers reported in previous works for some of our targets (MRC 0943-242, TN J0121+1320, TN J1338-1942 Wilman et al. 2004;Swinbank et al. 2015).Otherwise, we use the number determined in the first step of 1D spectra fitting.For TN J0121+1320, we use 3 absorbers instead of 2 as in Wilman et al. (2004).For 4C+04.11, we use 9 absorbers (instead of 7 as in Wang et al. (2021b)) with 7 of them fixed to the redshifts determined in 1D spectra fit.
Finally, we run MCMC sampling (Foreman-Mackey et al. 2013) using the results from last step as initials and with larger boundaries to probe the probability distribution of the fitted parameters and uncertainties.We point out one caveat that the low S/N (especially for the tiles at the edge) will affect the reconstructed intrinsic Lyα flux.We implemented a test where we artificially decrease the S/N of the master spectrum and do the absorption correction fit.It shows that the reconstructed intrinsic flux can vary by a factor of ∼ 2 when the S/N decreases by one dex (e.g.∼ 100 to ∼ 10).This result is based on the fact that we have a relatively good initial guess for the fit.We note that in the aforementioned fitting procedure the fitting parameters are constrained by the neighbouring tiles to have them physically linked.There may still be the chance that the low S/N will affect the fitted flux leading to over-correction.This may be the case for 4C+04.11(Fig. 12).We present the flux ratio maps between intrinsic and observed in Appendix C which indicate this possibility (∼ 50 for the tiles at the edge).
Fig. 1.Mapping results of our MUSE HzRGs sample.(a) Master Lyα spectrum (blue shaded histogram) extracted from a r = 0.5 arcsec aperture at the AGN position with best fit (solid dark magenta line).Red dashed curve shows the intrinsic Lyα from fitting, i.e. corrected for absorption.The vertical black bars above the emission line mark the positions of the H i absorbers.The yellow shaded region (if any) indicates the 5 wavelength pixel range excluded in the fitting due to the contamination from the 5577Å sky-line.The flux density unit, F λ , is 10 −20 erg s −1 cm −2 Å −1 .We also show the scaled He iiλ1640Å spectrum extracted from the same position in green histogram.We scale the peak flux density of He ii to 0.3 − 0.7 (varied for different targets) of the maximum peak flux density of observed Lyα spectrum in −1000 to 1000 km s −1 .The ∆v = 0 km s −1 is the systemic redshift based on He ii or [C i] (Table1,Kolwa et al. 2023).(b) Intrinsic Lyα surface brightness map.The flux in each tile is the integrated flux of the line emission corrected for absorption, i.e. total flux of the one or two Gaussians, see Sect.3.2.The light blue circle shows the aperture where the Master spectrum is extracted from.Green triangles mark the positions of the radio lobes.We place a green bar linking the triangles on TN J0121+1320 to indicate the unresolved state of its radio emission.The length of the bar represents the linear size of the 3σ contour along the east-west direction.The white hatched regions are the ones where the flux uncertainty is higher than 50% of the fitted intrinsic flux.The white bar indicates the 50 pkpc at the redshift of the radio galaxy.The unit of the surface brightness is 10 −16 erg s −1 cm −2 arcsec −2 .We apply the same color scale for all targets.(c) v 50 map of the intrinsic Lyα nebula.The zero velocity used for each target is determined by the systemic redshift (Table1).Green contours show the morphology of the radio jet in arbitrary values.The green cross mark the AGN position (Table2).(d) W 80 map of the intrinsic Lyα nebula.The black hatched regions on (c, d) are the same as (b).The purple hatched regions (in 4C+03.24and TN J1338-1942) are manually excluded due to contamination from either foreground star or known companion (Arrow galaxy in the filed of MRC0316-257, see Vernet et al. 2017).We note that the color scales for panel (c) and (d) are customised.The purple hatched area (if any) indicates the manually excluded region affected by foreground star or known Lyα emitter.
Fig. 1.Mapping results of our MUSE HzRGs sample.(a) Master Lyα spectrum (blue shaded histogram) extracted from a r = 0.5 arcsec aperture at the AGN position with best fit (solid dark magenta line).Red dashed curve shows the intrinsic Lyα from fitting, i.e. corrected for absorption.The vertical black bars above the emission line mark the positions of the H i absorbers.The yellow shaded region (if any) indicates the 5 wavelength pixel range excluded in the fitting due to the contamination from the 5577Å sky-line.The flux density unit, F λ , is 10 −20 erg s −1 cm −2 Å −1 .We also show the scaled He iiλ1640Å spectrum extracted from the same position in green histogram.We scale the peak flux density of He ii to 0.3 − 0.7 (varied for different targets) of the maximum peak flux density of observed Lyα spectrum in −1000 to 1000 km s −1 .The ∆v = 0 km s −1 is the systemic redshift based on He ii or [C i] (Table1,Kolwa et al. 2023).(b) Intrinsic Lyα surface brightness map.The flux in each tile is the integrated flux of the line emission corrected for absorption, i.e. total flux of the one or two Gaussians, see Sect.3.2.The light blue circle shows the aperture where the Master spectrum is extracted from.Green triangles mark the positions of the radio lobes.We place a green bar linking the triangles on TN J0121+1320 to indicate the unresolved state of its radio emission.The length of the bar represents the linear size of the 3σ contour along the east-west direction.The white hatched regions are the ones where the flux uncertainty is higher than 50% of the fitted intrinsic flux.The white bar indicates the 50 pkpc at the redshift of the radio galaxy.The unit of the surface brightness is 10 −16 erg s −1 cm −2 arcsec −2 .We apply the same color scale for all targets.(c) v 50 map of the intrinsic Lyα nebula.The zero velocity used for each target is determined by the systemic redshift (Table1).Green contours show the morphology of the radio jet in arbitrary values.The green cross mark the AGN position (Table2).(d) W 80 map of the intrinsic Lyα nebula.The black hatched regions on (c, d) are the same as (b).The purple hatched regions (in 4C+03.24and TN J1338-1942) are manually excluded due to contamination from either foreground star or known companion (Arrow galaxy in the filed of MRC0316-257, see Vernet et al. 2017).We note that the color scales for panel (c) and (d) are customised.The purple hatched area (if any) indicates the manually excluded region affected by foreground star or known Lyα emitter.

Fig. 2 .
Fig. 2. Example for the intrinsic mapping of the Lyα nebula of TNJ0205+2242.The central panel shows the intrinsic surface brightness map of TNJ0205+2242 which is the same as Fig. 1b.The green cross and triangles mark the position of the AGN and jet lobes, respectively.In each of the side panel, we show the spectrum (blue shade histogram in normalised flux unit) extracted from the individual spatial bin whose number is labeled at the top left, and the best fit (dark magenta curve) and recovered intrinsic Lyα (dashed red line).The black vertical bars indicate the positions of the H i absorbers.
. The selected quasars are all observed by advanced IFU instruments (MUSE or KCWI) and cover a large range of redshift and physical properties.They are luminous radio-quiet quasars at z ∼ 3.2 quasar from Borisova et al. (2016) (profiles from Marino et al. 2019), luminous type-1 quasars at z ∼ 3.17 from Arrigoni Battaia et al. (2019), luminous quasars at z ∼ 2.3 from Cai et al. (2019), high redshift quasars at z ∼ 6.28 from Farina et al. (2019) and luminous quasars at z ∼ 3.8 from Fossati et al. (2021).The quasar nebulae do not show so many absorption features as in our HzRGs and the studies were preformed without absorption corrections (Sect.5.1).Nevertheless, since the comparison quasar samples are not corrected for absorption, we do not show them all again in our intrinsic profile (lower) panel.The two exceptions are Farina et al. (2019) (F19) and Vayner et al. (2023) (7C 1354+2552, V22 7C

Fig. 5 .
Fig. 5. Radial profiles of the Lyα nebulae extracted in circular annuli (Fig.D.1).For better comparison, we show the radial profile in comoving kpc (ckpc) and take the cosmological dimming into account by a factor of (1 + z) 4 , where z is the redshift of the target.The black dot-dashed curve and grey dotted line in both panels are the best fitted exponential and power law profiles of the Arrigoni Battaia et al. (2019) radio loud sample, respectively.The two vertical dotted lines mark the 50 and 300 ckpc, respectively.Upper panel: Radial profile of observed surface brightness map in thicker solid lines.In this panel, we also include the radial profiles of other quasar samples (dashed lines) for comparison: B16 -Borisova et al. (2016), AB19 -Arrigoni Battaia et al. (2019), C19 -Cai et al. (2019), Farina19 -Farina et al. (2019), Fossati21 -Fossati et al. (2021), L22 -Lau et al. (2022) and V22 -Vayner et al. (2023) (two sources, 4C09.17 and 7C 1354+2552).When it is available, we show the range spanned by the 25th and 75th of the comparison sample radial profile as the shaded region around median profile in the same color.The horizontal bars at the right-most indicate the observed surface brightness limits (scaled by area from Table2) for each target in the same color.Lower panel: Intrinsic radial profile in thicker solid lines.The shaded regions around each profiles indicates the uncertainty range of the surface brightness from fitting.In this panel, we show again the same profiles of F19 and V22 7C as in the upper panel for comparison.Our HzRGs are extended further with higher surface brightnesses (or flattening in some sources) at larger radii (∼ 300 ckpc) compared to other samples.

Fig. 8 .
Fig.8.Zoom-in intrinsic surface brightness maps (Fig.1b) of the HzRGs sample to 15 × 15 arcsec 2 (or ∼ 110 × 110 pkpc 2 ) around the central AGN (blue cross).In each panel, the blue+red and green solid lines indicate the direction of, and perpendicular to, the radio jet.The blue (red) color represents the direction of the approaching (receding) jet.The white dashed line shows the flux-weighted position angle of the nebula (θ weight ).The white dotted line shows the unweighted position angle of the nebula (θ unweight ).The white star indicates the intrinsic flux-weighted centroid of the Lyα nebula.The flux weighted measurement is sensitive to the morphology of the high surface brightness part of the nebula.The unweighted measurement quantifies the morphology of the whole nebula.We find that nebula is elongated along the jet axis for most of HzRGs.

Fig. 9 .
Fig. 9. (a) Flux-weighted Lyα nebula elliptical asymmetry measurement versus nebula luminosity, L Lyα .We show the intrinsic fluxweighted ellipticity (e weight ) in larger open symbols for our targets versus their intrinsic Lyα luminosity.We also show the observed flux-weighted ellipticity e weight,obs for our targets versus their observed Lyα luminosity in smaller filler symbols.The small grey symbols are data of comparison targets (AB19 -Arrigoni Battaia et al. (2019), C19 -Cai et al. (2019) and L22 -Lau et al. (2022)).We mark the median fluxweighted (not corrected for absorption) ellipticity, 0.72, of type-1s with the horizontal dashed line.We also show the median e weight , 0.69, of radio-loud type-1 quasars from Arrigoni Battaia et al. (2019) in red horizontal dash-dotted line.(b) Flux-unweighted Lyα nebula elliptical asymmetry measurement versus L Lyα .The larger symbol are measurements for our HzRGs while the grey symbols are comparison targets (type-1s: B16 -Borisova et al. (2016), L22 -Lau et al. (2022) and M21 -Mackenzie et al. (2021); type-2s: dB20 -den Brok et al. (2020) and S21 -Sanderson et al. (2021)).We mark the median e unweight for type-1s (0.69) and type-2s (0.80) in solid and dashed horizontal lines, separately.The e weight is sensitive to the morphology of the high surface brightness part of the nebula while the e unweight quantifies the morphology of the whole nebula.We note that for e → 0, the nebula is closer to round shape and vice versa.At the bottom right, we show the median uncertainty of the intrinsic L Lyα for our sample in logarithmic scale, 0.04.The ellipticity for our sample are higher compared to the other quasars for both high surface brightness part and whole nebula.There is no clear evidence that the nebula ellipticity correlates with luminosity.
ports the morphological parameters of our sample.Since the absolute flux-weighted centroid position and θ weight (and θ unweight ) are irrelevant, we report the projected distance between the nebula centroid and the AGN position (d AGN−neb ) and the difference in angles between θ weight (and θ unweight ) and the jet position angle (|θ weight − PA radio | and |θ unweight − PA radio |), respectively.The jet position angle (PA radio ) is shown in Fig. 8 and listed in Appendix D.

Fig. 9
Fig.9presents the ellipticity measurements as a function of their nebula Lyα luminosity for our targets and other quasars(Borisova et al. 2016;Arrigoni Battaia et al. 2019;Cai et al. 2019;Mackenzie et al. 2021;den Brok et al. 2020;Sanderson et al. 2021;Lau et al. 2022).We note that the L Lyα for comparison samples are not corrected for absorption.Part of the comparison samples are also used in Sect.4.2.1 for surface brightness radial profile analysis.We point out the newly included ones here: faint z ∼ 3.0 type-1 fromMackenzie et al. (2021) and type-2 AGN at z ∼ 3.4 from den Brok et al. (2020) and z ∼ 3.2 fromSanderson et al. (2021).The reason why they are not included in radial profile analysis is that they do not add new information.The HzRGs from our sample are measured to be asymmetric for their high surface brightness emission region (median e weight ≈ 0.78).Compared to the Arrigoni Battaia et al. (2019) andCai et al. (2019) samples, our HzRGs are consistent in asymmetry measurements and on the higher end of their distribution (type-1 median e weight ≈ 0.72, dashed horizontal in Fig.9a).The e weight of radio-loud type-1s from (ArrigoniBattaia et al. 2019) have a median of 0.69 which is even lower than the value of all type-1 targets (ArrigoniBattaia et al. 2019;Cai et al. 2019).This indicates that the radio emission in type-1 does not disturb the gaseous nebula as in our HzRGs at least along the plane of the sky.This further suggests that orientation is a critical factor (Sect. 5.3).By comparing the intrinsic flux-weighted and observed flux-weighted elliptical asymmetry measurements, we find the e weight can vary significantly (e.g.MRC 0316-257 and 4C+04.11).For MRC0316-257, we can already identify

Fig. 10 .
Fig. 10.Radial profiles of the surface brightness Fourier decomposition (asymmetry measurement).The c 0 , c 1 and c 2 are the 0th, 1st and 2nd modes Fourier decomposition coefficients of the surface brightness radial profile, respectively (see den Brok et al. 2020, for definition).The (c 1 /c 0 ) 2 + (c 2 /c 0 ) 2 is a measurement of nebula asymmetry along the radial distance from the AGN.Our HzRGs are shown in solid color lines.r AGN is the radial distance measured from the central AGN.For comparison, we include the same measurements for the 4 nebulae of type-2 quasars from den Brok et al. (2020) (grey dashed lines) and the type-2 from Sanderson et al. (2021) (grey dot-dashed line).We also include the type-1 measurements from Borisova et al. (2016) (dotted line represents the median and shaded region marks the 25th and 75th percentile, quantified by den Brok et al. 2020).The vertical shaded region is the 0.75 arcsec (∼ 25 ckpc) range affected by median seeing of our sample (the radial distance where the type-1 PSF is affected is ∼ 50 ckpc, see den Brok et al. 2020).The morphologies for most of the HzRGs nebulae are round (symmetric) ≲ 100 ckpc and become asymmetric at larger radial distances ∼ 100 ckpc (see text).

Fig. 11 .
Fig. 11.(a) Intrinsic flux-weighted ellipticity (sensitive to the morphology of the high surface brightness part of the nebula), e weight , versus distance offset between nebular centroid and AGN position, d AGN−neb .The typical uncertainty of d AGN−neb (σ d AGN−neb = 0.4 pkpc, Table 3) is shown at bottom left.(b) Maximum Lyα nebula linear extent, d max , versus d AGN−neb .In both panels, we give the Spearman correlation measurements for our sample at the top left (the star superscript indicates the correlation is calculated excluding MRC 0316-257 in panel (a)).We also include the data from Arrigoni Battaia et al. (2019) in both panels in grey (radio quiet) and pink (radio loud) dots.We note that the Arrigoni Battaia et al. (2019) sample shown in this figure is incomplete to concentrate on the relation for our targets.There are positive correlations detected for e weight -d AGN−neb and d max -d AGN−neb .

Fig. 12 .
Fig. 12. Distribution of f Lyα,full FOV /Σ f Lyα,tile fit for our HzRGs.The f Lyα,full FOV is the intrinsic Lyα flux resulted from fitting the spectrum summed over the entire FOV.The Σ f Lyα,tile fit is the summation of the intrinsic Lyα flux in each tessellation bin (Fig.1b).The smaller the value, the more likely the Lyα photon is being double-counted when correcting for absorption.4C+04.11 is the one with the smallest ratio which may also indicate over-correction (Appendix B).

Fig. 13 .
Fig. 13.Distribution of the angle difference between nebular position angle and radio jet position angle, |θ − PA radio |.The blue histogram shows the number distribution of flux-weighted angle difference, |θ weight − PA radio | (sensitive to high surface brightness nebula morphology).The magenta histogram represents the unweighted measurements, |θ unweight −PA radio | (quantifying the morphology of the whole nebula).The values are presented in Table3.We mark the obvious outliers in the two distributions in corresponding colors.This shows that most of our HzRGs have their Lyα nebula elongated along the jet axis.
sample for comparison.The (r-values, p-values) of the e weight − d AGN−neb relation for the Arrigoni Battaia et al. (2019) whole sample and radio-loud targets are (0.3, 0.02) and (0.3, 0.4), respectively, implying no strong correlations.The (r-values, p-values) of the d max − d AGN−neb relation for the Arrigoni Battaia et al. ( Fig. C.4 and e.g.Wang et al. 2021b).
(2020) based on Lyα nebula studies of type-2 quasars (also included as comparison sample in this paper).den Brok et al. (2020) suggested that the orientation difference based on unification model can explain the nebula morphology at radial distance ≳ 40 pkpc.Using the same radial asymmetry measurement in Fig.10 (Sect.4.3.3),we also find large nebula asymmetries at ≳ 40 pkpc (∼ 200 ckpc) in our HzRGs following den Brok et al. (2020) type-2s.den Brok et al. ( Notes.( †) Maximum Gaussian kernel scale in FWHM.The unit is in pixel unit.
Fig. A.1.Schematic presentation of the detection map construction.(a) Average n λ number of images around each of the image slices (cyan) from the continuum-subtracted cube (dark green cube).(b) Spatially smooth each of the averaged image with Gaussian kernel.The algorithm will increase the kernel size (σ 0 < σ 1 ) to smooth the spaxels that not passing the T S/N , S/N threshold, at each steps until the maximum, σ max , has been reached.(c) Combine the adaptively smoothed images into the smoothed cube (lighter green cube).(d) Check the smoothed spectrum (long black rectangular box) at location (x, y).The position is preserved when there are n c consecutive number of pixels selected with signal detection (red 'S') in previous steps.(f) Construct the detection map.

Fig
Fig. A.2. Smoothed nebula image (a) and tessellation map (b) of MRC0943-242.(a) The purple, yellow and blue colours represent median pseudonarrow band Lyα images collapsed arbitrarily from red part, middle and blue part, respectively, of the smoothed cubes.
Fig. A.3.Smoothed nebula image (a) and tessellation map (b) of MRC0316-257.The inset spectrum is extracted from the detached regions (hatched) which is selected having Lyα emissions but left out from the analysis following Sect.3.1.1.(a) The color-composed image is created in a manner similar to that of Figure A.2a.

Fig
Fig. A.4. Smoothed nebula image (a) and tessellation map (b) of TN J0205+2242.The inset spectrum is extracted from the detached region (hatched) which is selected having Lyα emissions but left out from the analysis following Sect.3.1.1.(a) The color-composed image is created in a manner similar to that of Figure A.2a.
Fig. A.5. Smoothed nebula image (a) and tessellation map (b) of TN J0121+1320.(a) The color-composed image is created in a manner similar to that of Figure A.2a.

Fig
Fig. A.6.Smoothed nebula image (a) and tessellation map (b) of 4C+03.24.(a) The color-composed image is created in a manner similar to that of Figure A.2a.

Fig
Fig. A.8. Smoothed nebula image (a) and tessellation map (b) of TN J1338-1942.(a) The color-composed image is created in a manner similar to that of Figure A.2a.

Fig
Fig. E.1.UV continuum map around MRC0316-257 (a) and He ii spectra from the X-ray position (central AGN, b) and UV continuum peak position (c).The UV continuum map is collapsed between the observed wavelength of N vλ1240 and O i+Si iiλ1305.The green contours show the radio jet in the same format as Fig. 1cd.The black contours in the step of [3σ, 5σ, 7σ, ...] trace the UV continuum emission, where σ is the background standard deviation.Blue and red circular regions indicate the r = 0.5 arcsec apertures where the systemic and redshifted He ii spectra are extract, respectively.The right panels (b, c) show the He ii spectra (histogram) along with their best Gaussian fitting (dark magenta line) results.The velocity zero (vertical black dotted line) in both panels is the systemic redshift.In the panel(c), we also mark the velocity shift (vertical red dotted line) of the redshifted He ii emission with respect to the systemic one.

Table 1 .
Details of the MUSE observation of the HzRG sample.