# 3D tomography of the giant Lyα nebulae of $z \approx 3-5$
radio-loud AGN

Wuji Wang, Dominika Wylezalek, Joël Vernet, Carlos de Breuck, Bitten Gullberg, Mark Swinbank, Montserrat Villar Martín, Matthew D Lehnert, Guillaume Drouart, Fabrizio Arrigoni Battaia, et al.

<figure><img src="image_1.png" /></figure>

► To cite this version:

Wuji Wang, Dominika Wylezalek, Joël Vernet, Carlos de Breuck, Bitten Gullberg, et al., 3D tomography of the giant Lyα nebulae of $z \approx 3-5$ radio-loud AGN. Astronomy & Astrophysics - A&A, 2023, 680, pp.A70. 10.1051/0004-6361/202346415. hal-05018796

HAL Id: hal-05018796

https://hal.science/hal-05018796v1

Submitted on 3 Apr 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# 3D tomography of the giant Lyα nebulae of $z \approx 3-5$ radio-loud AGN

Wuji Wang (王无忌)<fnref n="1" />, Dominika Wylezalek<fnref n="1" />, Joël Vernet<fnref n="2" />, Carlos De Breuck<fnref n="2" />, Bitten Gullberg<fnref n="3" />,<fnref n="4" />,
Mark Swinbank<fnref n="5" />, Montserrat Villar Martín<fnref n="6" />, Matthew D. Lehnert<fnref n="7" />, Guillaume Drouart<fnref n="8" />, Fabrizio Arrigoni Battaia<fnref n="9" />,
Andrew Humphrey<fnref n="10" />, Gaël Noirot<fnref n="11" />, Sthabile Kolwa<fnref n="12" />, Nick Seymour<fnref n="8" />, and Patricio Lagos<fnref n="10" />

Received 14 March 2023 / Accepted 26 September 2023

## ABSTRACT

Lyα emission nebulae are ubiquitous around high-redshift galaxies and are tracers of the gaseous environment on scales out to $\ge 100$ pkpc (proper kiloparsec). High-redshift radio galaxies (HzRGs, type-2 radio-loud quasars) host large-scale nebulae observed in the ionised gas differ 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 of the HzRGs have large-scale Lyα emission nebulae with seven of them extended over 100 pkpc at the observed surface brightness limit ($\sim 2-20 \times 10^{-19} \text{ erg s}^{-1} \text{ cm}^{-2} \text{ arcsec}^{-2}$). Because the emission line profiles are significantly affected by neutral hydrogen absorbers across the entire nebulae extent, we performed 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 the following: (i) that 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 a 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 active galactic nuclei (AGN) and the centroids of the Lyα nebula; and (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 affect 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 affects the morphology and kinematics of the nebulae. Overall, this work provides results that allow us to compare Lyα nebulae around various classes of quasars at and beyond cosmic noon ($z \sim 3$).

**Key words.** galaxies: active – galaxies: evolution – galaxies: high-redshift – galaxies: halos – galaxies: jets

## 1. Introduction

Being the most abundant element in the Universe, hydrogen (especially the cold gas, i.e. neutral hydrogen atoms and hydrogen molecules, $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 that 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_{H0} = 13.6$ eV) and cascade with line emissions being produced. The H I Ly$\alpha\lambda1216$ (Ly$\alpha$ hereafter) line is the most prominent

one among them. For high-redshift galaxies, it is a commonly targeted emission line that 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$\alpha$ 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$\alpha$ emission

<footnote n="1">¹ Astronomisches Rechen-Institut, Zentrum für Astronomie der Universität Heidelberg, Mönchhofstr. 12-14, 69120 Heidelberg, Germany
e-mail: wuji.wang@uni-heidelberg.de; wuji.wang_astro@outlook.com</footnote>

<footnote n="2">² European Southern Observatory, Karl-Schwarzschild-Str. 2, 85748 Garching, Germany</footnote>

<footnote n="3">³ Cosmic Dawn Center, Copenhagen, Denmark</footnote>

<footnote n="4">⁴ DTU Space, Technical University of Denmark, Elektronvej 327, 2800 Lyngby, Denmark</footnote>

<footnote n="5">⁵ Centre for Extragalactic Astronomy, Department of Physics, Durham University, South Road, Durham DH1 3LE, UK</footnote>

<footnote n="6">⁶ Centro de Astrobiología, CSIC-INTA, Ctra. de Torrejón a Ajalvir, km 4, 28850 Torrejón de Ardoz, Madrid, Spain</footnote>

<footnote n="7">⁷ Univ. Lyon, Univ. Lyon1, ENS de Lyon, CNRS, Centre de Recherche Astrophysique de Lyon UMR5574, 69230 Saint-Genis-Laval, France</footnote>

<footnote n="8">⁸ International Centre for Radio Astronomy Research, Curtin University, 1 Turner Avenue, Bentley, WA 6102, Australia</footnote>

<footnote n="9">⁹ Max-Planck-Institut fur Astrophysik, Karl-Schwarzschild-Str 1, 85748 Garching bei München, Germany</footnote>

<footnote n="10">¹⁰ Instituto de Astrofísica e Ciências do Espaço, Universidade do Porto, CAUP, Rua das Estrelas, 4150-762 Porto, Portugal</footnote>

<footnote n="11">¹¹ Institute for Computational Astrophysics and Department of Astronomy & Physics, Saint Mary's University, 923 Robie Street, Halifax, NS B3H 3C3, Canada</footnote>

<footnote n="12">¹² Physics Department, University of Johannesburg, 5 Kingsway Ave, Rossmore, Johannesburg 2092, South Africa</footnote>
(e.g. due to the resonant nature of $Ly\alpha$ emission and various potential ionising sources 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 active galactic nuclei (AGN) are present.

Active galaxies hosting AGN, especially the ones with quasar level activities (bolometric luminosity, $L_{bol} \ge 10^{45}$ erg s$^{-1}$), at a high redshift are known to host $Ly\alpha$ nebulae on scales of a few 100 kpc (e.g. Heckman et al. 1991a; Basu-Zych & Scharf 2004; Weidinger et al. 2004, 2005; Dey et al. 2005; Prescott et al. 2015; Cantalupo et al. 2014; Arrigoni Battaia et al. 2016, 2019; Borisova et al. 2016; Cai 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\alpha$ nebulae (as predicted by theoretical works, e.g. Costa et al. 2022). In addition, the diffuse emission from galaxies near 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\alpha$ acts as a simpler tool for detecting a gaseous environment throughout 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\alpha$ nebulae (e.g. Steidel et al. 2000; Francis et al. 2001; Matsuda et al. 2004; Saito et al. 2006; Yang et al. 2009, 2010; Cantalupo et al. 2012, 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 (IFU) observations – for example with the Multi-Unit Spectroscopic Explorer (MUSE/VLT) and Keck Cosmic Web Imager (KCWI/Keck) – allow us to measure the extent of the nebulae together with the information of their dynamics. Numerous works of $Ly\alpha$ nebulae around quasars report (tens of kiloparsecs to over 100 kpc) extended emission across a large range of redshifts ($z \sim 2$ to $z \sim 6.3$) and quasar types (e.g. radio-quiet and radio-loud type-1, radio-quiet type-2, and extremely red quasars, 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 that has been made regarding this topic, a link between the aforementioned types and type-2 radio-loud quasars on a CGM scale is missing.

Among the high-redshift quasar population, high-redshift radio galaxies (HzRGs) are a unique sample despite being 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. 2013, 2014; Noirot et al. 2016, 2018), which may evolve to modern galaxy clusters. HzRGs were among the first sources where giant $Ly\alpha$ nebulae were discovered (~$10^{44}$ erg s$^{-1}$, $\ge 100$ kpc, e.g. Hippelein & Meisenheimer 1993; van Ojik et al. 1996, 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\alpha$ 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\alpha$ 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\alpha$ 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\alpha$ nebulae of HzRGs are known to be partially absorbed by neutral hydrogen (H1 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 this absorption is only possible through spectral observation. Without careful treatment, a considerable amount (a factor of $\ge 5$) of flux would be missed, and inaccurate conclusions would be drawn. Alternatively, some absorption trough features might potentially be explained by radiative transfer effects (Dijkstra 2014; Gronke et al. 2015, 2016; Gronke & Dijkstra 2016). Although it is interesting to compare the different treatments of the observed $Ly\alpha$ 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 the jet direction (e.g. Chambers et al. 1987; Pentericci et al. 1999; Miley et al. 2004; Zirm et al. 2005; Duncan et al. 2023) on a relatively smaller scale (several kiloparsecs to tens of kiloparsecs). Molecular gas detected around HzRGs was reported to be distributed along the jet within and outside the hot spot, which may suggest several scenarios (e.g. jet-driven outflow, jet-induced gas cooling, and a jet propagating into a dense molecular gas medium, Emonts et al. 2014; Gullberg et al. 2016; Falkendal et al. 2021). On a megaparsec 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 hundreds of kiloparsec gas tracers of $Ly\alpha$, we were able to probe the intrinsic (i.e. corrected for absorption) gaseous nebula around HzRGs for this work, test its distribution with respect to the radio jets, and seek evidence following these pioneering works.

For this paper, we utilised the power of MUSE IFU to fully map the $Ly\alpha$ emission nebulae of a sample of HzRGs over a redshift range of 2.92–4.51 and initiated 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 measured 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
<table><caption>**Table 1.** Details of the MUSE observation of the HzRG sample.</caption><thead><tr><th>HzRG</th><th>Redshift z<sup>(†)</sup></th><th>UT date (dd/mm/yyyy)</th><th>Program ID</th><th>Mode</th><th>Exp. time total hours</th><th>Seeing arcsec<sup>(*)</sup></th></tr></thead><tbody><tr><td>MRC 0943-242</td><td>2.9230</td><td>21/02/2014</td><td>60.A-9100(A)</td><td>WFM-NOAO-E</td><td>5.21</td><td>0.65<sup>(a)</sup></td></tr><tr><td>—</td><td>—</td><td>15/12/2015-18/01/2016</td><td>096.B-0752(A)</td><td>—</td><td>—</td><td>—</td></tr><tr><td>MRC 0316-257</td><td>3.1238</td><td>15/01-17/01/2015</td><td>094.B-0699(A)</td><td>WFM-NOAO-N</td><td>4.24</td><td>0.61<sup>(b)</sup></td></tr><tr><td>TN J0205+2242</td><td>3.5060</td><td>03/12-08/12/2015</td><td>096.B-0752(B)</td><td>WFM-NOAO-N</td><td>4.24</td><td>0.73</td></tr><tr><td>TN J0121+1320</td><td>3.5190</td><td>06/10/2015</td><td>096.B-0752(C)</td><td>WFM-NOAO-N</td><td>5.30</td><td>0.83</td></tr><tr><td>—</td><td>—</td><td>08/08-28/08/2016</td><td>097.B-0323(C)</td><td>—</td><td>—</td><td>—</td></tr><tr><td>4C+03.24</td><td>3.5828</td><td>17/06-18/06/2017</td><td>60.A-9100(G)</td><td>WFM-AO-N</td><td>1.25</td><td>0.63</td></tr><tr><td>4C+19.71</td><td>3.5892</td><td>08/06-02/09/2016</td><td>097.B-0323(B)</td><td>WFM-NOAO-N</td><td>5.83</td><td>1.03</td></tr><tr><td>TN J1338-1942</td><td>4.0959</td><td>30/04-06/05/2014</td><td>60.A-9100(B)</td><td>WFM-NOAO-N</td><td>8.93</td><td>0.77<sup>(a)</sup></td></tr><tr><td>—</td><td>—</td><td>30/06/2014</td><td>60.A-9318(A)</td><td>—</td><td>—</td><td>—</td></tr><tr><td>4C+04.11</td><td>4.5077</td><td>03-15/12/2015</td><td>096.B-0752(F)</td><td>WFM-NOAO-N</td><td>4.24</td><td>0.88</td></tr></tbody></table>

**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 colour, 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.

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 \text{ km s}^{-1} \text{ Mpc}^{-1}$ and $\Omega_m = 0.3$. Following this cosmology, $1 \text{ arcsec} \approx 6.6-7.7 \text{ pkpc}$ for our sample redshifts. Throughout the paper, pkpc stands for proper kiloparsec and ckpc represents comoving kiloparsec, $\text{ckpc} = (1+z)\text{pkpc}$. In this paper, we use ‘intrinsic’ to refer to the absorption-corrected Lyα emission.

## 2. HzRGs sample, observations, and data processing

### 2.1. MUSE HzRGs sample

#### 2.1.1. 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 $\text{Dec} < 25^{\circ}$ 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_{\odot}$ for all targets, De Breuck et al. 2010) and the star formation rates ranging from uppers limit of $<84 M_{\odot} \text{yr}^{-1}$ to constraints of $626 M_{\odot} \text{yr}^{-1}$ (Falkendal et al. 2019).

#### 2.1.2. 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 estimat-

ing the bolometric luminosity of AGN, $L_{\text{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 \mu 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_{\text{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_{\text{bol,AGN}}$ for our sample. We report that the $L_{\text{bol,AGN}}$ estimated for our sample using those different methods varies from $10^{45.9}$ to $10^{48.5} \text{ erg s}^{-1}$. Given this large uncertainty, we find it is unreasonable to draw further conclusions from the comparison of $L_{\text{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_{\text{bol,AGN}}$ is beyond the scope of this work and may be done more thoroughly through multi-wavelength approach.

#### 2.1.3. 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, 2008, 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).

## 2.2. MUSE observations

In this work, we analyse data from MUSE integral field spectrograph (Bacon et al. 2010, 2014) mounted on the ESO Very Large Telescopes (VLT) Yepun (UT4). All observations were carried out in Wide-Field Mode (WFM) offering a $1 \times 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 $\Delta v \sim 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$\alpha$ emission ($L_{Ly\alpha, 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 min were slightly dithered (with a $<1''$ amplitude pattern) and rotated by 90 degrees from each other.

## 2.3. 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$\alpha$ 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 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 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$\alpha$ line emission at the outskirt of our sample. Since the Ly$\alpha$ 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$\alpha$ mask that covers most line-emission region<fnref n="1" />; (iii) re-run ZAP using this Ly$\alpha$ 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 \approx 10''$ circular aperture around the radio galaxy masking all continuum sources falling in the aperture; (ii) mask the Ly$\alpha$ 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.

# 3. Data analysis

## 3.1. Ly$\alpha$ nebulae extent and tessellation

To systematically study the Ly$\alpha$ nebulae of our HRGs 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 signal-to-noise ratio (S/N) using a tessellation technique (Sect. 3.1.2) before fitting the emission feature described in Sect. 3.2.

### 3.1.1. Maximum extent of the nebulae

To select the Ly$\alpha$ 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_\lambda$ neighbouring pixels. Then for each wavelength plane, the algorithm iteratively smoothes spatially with a growing gaussian kernel selecting pixels passing a given S/N threshold ($T_{S/N}$) and leaving to the next iterations only spaxels below this S/N threshold, until a maximum smoothing radius is reached ($\sigma_{max}$). The spaxels not selected by the end of the iterative process are

<footnote n="1">¹ We note that this mask is only used in this process to eliminate the impact of the Ly$\alpha$ signal on ZAP. The detection map for determining the maximum extent of the Ly$\alpha$ nebula is described in Sect. 3.1.1.</footnote>
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.

To determine the optimal combination of the four parameters ($n_\lambda$, $T_{S/N}$, $\sigma_{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).

### 3.1.2. Tessellation procedure

In order to increase the S/N to a level that allows fitting of the Ly$\alpha$ line, especially close to the detection limit at the periphery of the nebulae, we tessellate the Ly$\alpha$ detection map. To construct the tessellation map, we firstly use a S/N map based on a narrow-band image (~15 Å wide) extracted around the Ly$\alpha$ 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 two-step 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$\alpha$ 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.

## 3.2. Spectral fitting

### 3.2.1. Ly$\alpha$ absorption modelling

In this work we treat the Ly$\alpha$ 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$\alpha$ profiles. Indeed, such absorption troughs are observed in our Ly$\alpha$ 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 term ‘radiative transfer’ used in the paper refers to the process where Ly$\alpha$ 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$\alpha$ 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$\alpha$ 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).

We emphasise that the term ‘intrinsic Ly$\alpha$ emission’ throughout the work refers to the nebula Ly$\alpha$ emission corrected for intervening absorbers. The absorption troughs seen on the spectra (Fig. 1a) are due to the Ly$\alpha$ 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 observe these features continuously across the extent of the nebulae in most HzRGs indeed indicates they are coherent intergalactic-scale structures. This can be found in Figs. 2–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 al. 2021b). Our approach is a common interpretation in studies of HzRGs. Conversely, such absorbers are not often seen in the Ly$\alpha$ 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 modelling works (e.g. Arrigoni Battaia et al. 2019; Chang et al. 2023). This could be a result of prior radiative transfer effects of Ly$\alpha$ (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$\alpha$ (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$\alpha$ 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
<figure><img src="image_3.png" /></figure>

**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*<sub>λ</sub>, is 10<sup>-20</sup> erg s<sup>-1</sup> cm<sup>-2</sup> Å<sup>-1</sup>. 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<sup>-1</sup>. The Δ*v* = 0 km s<sup>-1</sup> is the systemic redshift based on He II or [C I] (Table 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<sup>-16</sup> erg s<sup>-1</sup> cm<sup>-2</sup> arcsec<sup>-2</sup>. We apply the same colour scale for all targets. (c) *v*<sub>50</sub> 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). (d) *W*<sub>80</sub> 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.24 and TN J1338-1942) are manually excluded due to contamination from either foreground star or known companion (Arrow galaxy in the field of MRC0316-257, see Vernet et al. 2017). We note that the colour scales for panels c and d are customised. The purple hatched area (if any) indicates the manually excluded region affected by foreground star or known Lyα emitter.

assumptions that absorption rather than radiative transfer is pri-
marily responsible for the line profiles. We discuss the limitation
of this treatment in Sect. 5.1.

We note that the O v]λλ1213.8, 1218.2 (O v]) line under-
neath the Lyα can affect the obtained flux especially in the
nuclear region where the ionisation parameter (and metallicic-
ity) 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 fit-
ting. 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.

3.2.2. Fitting procedure

To reconstruct the intrinsic Lyα emission across the nebula, we
fit each spectrum in each tessellation bin (see Sect. 3.1.2) follow-
ing 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
<figure><img src="image_4.png" /><figcaption>Fig. 1. continued.</figcaption></figure>

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 Figs. 2–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–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. Nvλλ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 χ². The fit is done using a least-squares method followed by a Markov chain Monte Carlo (MCMC, using the python package `emcee`, Foreman-Mackey et al. 2013) sampling. 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 Appendix 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.

# 4. Results

## 4.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 source 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 ($\Delta v = 0$ km s⁻¹, 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 selected exemplar individual fits are shown in Figs. 2–4 (also see Appendix C).

The intrinsic Lyα surface brightness maps are shown in Fig. 1b on the same flux scale. Regions with larger fitting
<figure><img src="image_5.png" /><figcaption>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 labelled 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.</figcaption></figure>

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 \times 10^{44}$ erg s⁻¹.

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 $\Phi(v)$ defined as:

$$ \Phi(v) = \int_{-\infty}^{v} f(v')dv' \qquad (1) $$

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, $\Phi(v_{50}) = 0.5\Phi(\infty)$. 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.71 and 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,
<figure><img src="image_6.png" /><figcaption>Fig. 3. Similar to Fig. 2, but for TNJ0121+1320.</figcaption></figure>

for example the south-west tile of 4C+04.11, displaying larger $W_{80}$ values (≥2500 km s⁻²). Except 4C+19.71, all targets show a line width of ~800–2500 km s⁻¹. 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 limit (≥600km s⁻¹) 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.

## 4.2. Radial profiles

### 4.2.1. 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 (an extremely red quasar and 2 radio-loud quasars from Lau et al. 2022; Vayner et al. 2023, respectively). 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 \sim 3.2$ quasar from Borisova et al. (2016; profiles from Marino et al. 2019), luminous type-1 quasars at $z \sim 3.17$ from Arrigoni Battaia et al. (2019), luminous quasars at $z \sim 2.3$ from Cai et al. (2019), high redshift quasars at $z \sim 6.28$ from Farina et al. (2019) and luminous quasars at $z \sim 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 performed 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, hereafter F19) and Vayner et al. (2023; 7C 1354+2552, V22 7C). Those two are on the higher surface brightness end of the comparison samples, and we 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⁻² for
<figure><img src="image_7.png" /><figcaption>**Fig. 4.** Similar to Fig. 2, but for 4C+03.24. The red box marks the secondary southern K-band continuum emission peak (van Breugel et al. 1998, Sect. 4.2.4). The yellow shaded regions show the wavelength range excluded due to the 5577 Å sky-line.</figcaption></figure>

<table><caption>**Table 2.** HzRGs MUSE sample properties.</caption><thead><tr><th>HzRG</th><th>RA (J2000) <sup>(†)</sup><br>hh:mm:ss</th><th>Dec (J2000) <sup>(†)</sup><br>dd:mm:ss</th><th>SB limit <sup>(*)</sup><br>10<sup>-19</sup> cgs</th><th>L<sub>Lyα, int</sub> <sup>(‡)</sup><br>10<sup>44</sup> erg s<sup>-1</sup></th><th>L<sub>Lyα, obs</sub> <sup>(§)</sup><br>10<sup>44</sup> erg s<sup>-1</sup></th></tr></thead><tbody><tr><td>MRC 0943-242</td><td>09:45:32.73</td><td>-24:28:49.65 <sup>(a)</sup></td><td>17.0</td><td>4.7±0.1</td><td>2.39±0.03</td></tr><tr><td>MRC 0316-257</td><td>03:18:12.07</td><td>-25:35:10.22 <sup>(b)</sup></td><td>2.62</td><td>7.3±1.2</td><td>1.04±0.04</td></tr><tr><td>TN J0205+2242</td><td>02:05:10.69</td><td>+22:42:50.4 <sup>(a), (c)</sup></td><td>13.3</td><td>8.3±0.4</td><td>3.82±0.06</td></tr><tr><td>TN J0121+1320</td><td>01:21:42.73</td><td>+13:20:58.0 <sup>(a)</sup></td><td>5.64</td><td>2.0±0.1</td><td>0.55±0.01</td></tr><tr><td>4C+03.24</td><td>12:45:38.37</td><td>+03:23:21.0 <sup>(d)</sup></td><td>11.5</td><td>28.8±0.8</td><td>5.43±0.20</td></tr><tr><td>4C+19.71</td><td>21:44:07.56</td><td>+19:29:14.6 <sup>(a)</sup></td><td>4.77</td><td>4.3±0.3</td><td>1.48±0.22</td></tr><tr><td>TN J1338-1942</td><td>13:38:26.10</td><td>-19:42:31.1 <sup>(e)</sup></td><td>4.34</td><td>12.3±0.3</td><td>5.89±0.10</td></tr><tr><td>4C+04.11</td><td>03:11:47.97</td><td>+05:08:03.74 <sup>(f)</sup></td><td>9.84</td><td>20.0±0.7</td><td>2.89±0.07</td></tr></tbody></table>

**Notes.** <sup>(†)</sup>Position of the AGN and/or host galaxy. <sup>(*)</sup>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 Sect. 4.1). The cgs unit is erg s$^{-1}$ cm$^{-2}$ arcsec$^{-2}$. <sup>(‡)</sup>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). <sup>(§)</sup>Observed Lyα luminosity of the nebula integrated over the entire area from $v_{05}$ to $v_{95}$ (see text).

**References.** <sup>(a)</sup>Spitzer (Seymour et al. 2007; De Breuck et al. 2010). <sup>(b)</sup>Chandra Obs. ID 5734 (PI:Pentericci). Chandra Source Catalogue (Evans et al. 2010). <sup>(c)</sup>Radio (De Breuck et al. 2002). <sup>(d)</sup>Radio (van Ojik et al. 1996). <sup>(e)</sup>Radio (De Breuck et al. 1999). <sup>(f)</sup>Radio (Parijskij et al. 2014).
<table><caption>**Table 3.** HzRGs nebulae properties.</caption><thead><tr><th>HzRG</th><th>e<sub>weight</sub><sup>(†)</sup></th><th>e<sub>weight, obs</sub><sup>(†)</sup></th><th>e<sub>unweight</sub><sup>(†)</sup></th><th>|θ<sub>w.</sub> - PA<sub>radio</sub>|<sup>(*)</sup><br>deg</th><th>|θ<sub>unw.</sub> - PA<sub>radio</sub>|<sup>(*)</sup><br>deg</th><th>d<sub>AGN-neb</sub><sup>(‡)</sup><br>pkpc</th><th>d<sub>max</sub><sup>(§)</sup><br>pkpc</th></tr></thead><tbody><tr><td>MRC 0943-242</td><td>0.73±0.01</td><td>0.76</td><td>0.62</td><td>4.7±1.8</td><td>7.2</td><td>3.4±0.1</td><td>164</td></tr><tr><td>MRC 0316-257</td><td>0.28±0.38</td><td>0.96</td><td>0.87</td><td>28±54</td><td>20.8</td><td>13.5±0.2</td><td>347</td></tr><tr><td>TN J0205+2242</td><td>0.91±0.01</td><td>0.79</td><td>0.75</td><td>2.1±2.3</td><td>8.6</td><td>9.0±0.1</td><td>170</td></tr><tr><td>TN J0121+1320</td><td>0.66±0.05</td><td>0.74</td><td>0.48</td><td>10.5±8.7</td><td>72.1</td><td>3.1±0.1</td><td>72</td></tr><tr><td>4C+03.24</td><td>0.79±0.01</td><td>0.72</td><td>0.82</td><td>8.1±1.2</td><td>2.0</td><td>9.9±0.1</td><td>248</td></tr><tr><td>4C+19.71</td><td>0.85±0.02</td><td>0.83</td><td>0.51</td><td>20.1±9.1</td><td>29.3</td><td>6.7±0.1</td><td>143</td></tr><tr><td>TN J1338-1942</td><td>0.94±0.01</td><td>0.87</td><td>0.90</td><td>1.3±1.4</td><td>1.2</td><td>15.3±0.1</td><td>265</td></tr><tr><td>4C+04.11</td><td>0.76±0.01</td><td>0.20</td><td>0.59</td><td>81.4±1.3</td><td>45.3</td><td>5.7±0.1</td><td>163</td></tr></tbody></table>

**Notes.** (†)Intrinsic flux-weighted ($e_{weight}$), observed flux-weighted ($e_{weight, obs}$) and unweighted elliptical asymmetric measurements. The flux-weighted 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 \to 0$, the nebula is closer to round shape and vice versa. (*)Absolute difference between flux-weighted (unweighted) nebula position angle, $\theta_{weight}$ ($\theta_{unweight}$), and radio axis position angle, PA<sub>radio</sub>. (‡)Offset between the intrinsic flux-weighted centroid of Ly$\alpha$ nebula and AGN position. (§)Maximum linear extent spanned by the Ly$\alpha$ nebula.

the column densities. For F19, there is not much evidence of absorption. Therefore, the comparison is legitimate. The best fit profiles to the observed Ly$\alpha$ nebulae of radio loud quasars in Arrigoni Battaia 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 extremely red quasar 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$\alpha$ nebulae (Sect. 2.1). We note that the jet compression is also known to result in high Ly$\alpha$ 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 similar extend Ly$\alpha$ 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$\alpha$ 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<sup>-15</sup> erg s<sup>-1</sup> cm<sup>-2</sup> arcsec<sup>-2</sup> (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 extended Ly$\alpha$ 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$\alpha$ 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 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<sup>-1</sup> in Cai et al. (2019), 30 Å in Arrigoni Battaia et al. (2019) and ±500 km s<sup>-1</sup> 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<sup>-1</sup> 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
<figure><img src="image_8.png" /></figure>

**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 colour. The horizontal bars at the right-most indicate the observed surface brightness limits (scaled by area from Table 2) for each target in the same colour. 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.
<figure><img src="image_9.png" /><figcaption>Fig. 6. Surface brightness radial profiles for approaching (blue squares) and receding (red circles) directions along the jet axis. The dotted curves in corresponding colours 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 colours. 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.24 and 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).</figcaption></figure>

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.

### 4.2.2. 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 direction-dependent 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<fnref n="2" />. Qualitatively 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.24 and 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.

### 4.2.3. 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

$$ \mathrm{SB}(r) = \begin{cases} C_e \exp(-r/r_h) & r \le r_b \\ C_p r^{\alpha} & r > r_b \end{cases}, \qquad (2) $$

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^\alpha$). The determination of the piecewise function is motivated by previous studies

<footnote n="2">² We note that the radio emission of TN J0121+1320 is unresolved. The 'jet size' represented by the vertical lines are linear size of the $3\sigma$ contour along the east-west elongation of the radio map.</footnote>
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 $\chi^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_h$ (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 further discuss the power law declining (flattening) part of our nebula in Sect. 5 combining the information from nebular morphology (Sect. 4.3).

### 4.2.4. 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 jet-driven outflows (e.g. Nesvadba et al. 2008, 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 towards 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} \sim 0$ ckpc from $0 \text{ 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.24 which has an offset of ~900 km s$^{-1}$. We note that its systemic redshift (Table 1) is based on [CII](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\beta$, $z_{H\beta} \simeq 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 \text{ km s}^{-1}$ with respective to the systemic velocity shift used in this paper. We caution that, however, the $H\beta$ 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\beta} \simeq 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.71 whose 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 $\ge 22$ pkpc) where most of our fit are done with a single Gaussian (Sect. 3.2). In 4C+03.24, 4C+19.71 and 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
<figure><img src="image_10.png" /><figcaption>**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. 1c). 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.24 indicates the velocity shift of H$\beta$ redshift ($z_{H\beta} \approx 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.71 and TN J1338-1942), $W_{80}$ decreases with increasing radial distance which may indicate that the jet is disturbing the gas.</figcaption></figure>

derived without correction for absorption. Since we see a high
spatial coverage of absorbers (Figs. 2–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 bright-
ness 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 bright-
ness 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<sub>50</sub> and W<sub>80</sub>), while MRC0316-257 has a decrease in W<sub>80</sub> 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<sub>50</sub> and W<sub>80</sub> 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<sup>4</sup> pkpc<sup>2</sup> 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 centre which has larger tile size. Nevertheless, we note that the detected W<sub>80</sub> of ~10<sup>5</sup> km s<sup>-1</sup> (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).

## 4.3. 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 asymmetry 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.

### 4.3.1. 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, $\alpha_{weight}$ (see Arrigoni Battaia et al. 2019, for the definition). This quantifies the asymmetry of the nebula in two perpendicular directions. Together with the asymmetry measurement, we also obtain the flux-weighted position angle $\theta_{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 $\alpha_{unweight}$ following den Brok et al. (2020). In this context, we also derive $\theta_{unweight}$ (flux-unweighted position angle) to examine the jet-nebula relation with respect to the entire nebula. Figure 8 visualises the weighted (nebula centroids and $\theta_{weight}$) and unweighted ($\theta_{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 flux-weighted and unweighted), $\alpha$, to an intuitive elliptical asymmetry measurement (or ellipticity) $e = \sqrt{1 - \alpha^2}$. For $e_{weight} \to 0$, the nebula is closer to round shape and vice versa. Table 3 reports the morphological parameters of our sample. Since the absolute flux-weighted centroid position and $\theta_{weight}$ (and $\theta_{unweight}$) are irrelevant, we report the projected distance between the nebula centroid and the AGN position ($d_{AGN-nev}$) and the difference in angles between $\theta_{weight}$ (and $\theta_{unweight}$) and the jet position angle ($|\theta_{weight} - \text{PA}_{\text{radio}}|$ and $|\theta_{unweight} - \text{PA}_{\text{radio}}|$), respectively. The jet position angle ($\text{PA}_{\text{radio}}$) is shown in Fig. 8 and listed in Appendix D.

### 4.3.2. Comparison of nebula asymmetry with other quasar samples

Figure 9 presents 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\alpha}$ 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 \sim 3.0$ type-1 from Mackenzie et al. (2021) and type-2 AGN at $z \sim 3.4$ from den Brok et al. (2020) and $z \sim 3.2$ from Sanderson 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} \approx 0.78$). Compared to the Arrigoni Battaia et al. (2019) and Cai et al. (2019) samples, our HzRGs are consistent in asymmetry measurements and on the higher end of their distribution (type-1 median $e_{weight} \approx 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 flux-weighted ellipticity. Hence, the morphology of MRC0316-
<figure><img src="image_11.png" /><figcaption>Fig. 8. Zoom-in intrinsic surface brightness maps (Fig. 1b) of the HzRGs sample to 15 × 15 arcsec² (or ~110 × 110 pkpc²) 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) colour represents the direction of the approaching (receding) jet. The white dashed line shows the flux-weighted position angle of the nebula (θ<sub>weight</sub>). The white dotted line shows the unweighted position angle of the nebula (θ<sub>unweight</sub>). 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.</figcaption></figure>

257 is more towards the asymmetric end. The large difference between its intrinsic and observed e<sub>weight</sub> 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).

Our HzRGs have a median e<sub>unweight</sub> ≈ 0.70 which is lower than the measurement of type-2s and relatively consistent with the type-1s (median e<sub>unweight</sub> ≈ 0.80 and ≈0.69, respectively, Borisova et al. 2016; den Brok et al. 2020; Sanderson et al. 2021; Mackenzie et al. 2021). The ellipticity 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<sub>unweight</sub> measured for our sample with four clustered around the type-2s and other four below the median e<sub>unweight</sub> 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<sub>unweight</sub> (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<sub>unweight</sub>. 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<sub>Lyα</sub>, our HzRGs host the most luminous Lyα nebula compared with other quasar types. This is due to the following: (i) sample selection<fnref n="3" />; (ii) absorption correction; and (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<sub>weight</sub> (0.44) and e<sub>unweight</sub> (0.69) of L22 are lower than the measurements from HzRGs. Its L<sub>Lyα</sub> 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.

### 4.3.3. 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<sub>k</sub>(r) and b<sub>k</sub>(r):

$$ \mathrm{SB}_{\mathrm{Ly}\alpha}(r, \theta) = \sum_{k=0}^{\infty} [a_k(r) \cdot \cos(k\theta) + b_k(r) \cdot \sin(k\theta)], \quad (3) $$

where SB<sub>Lyα</sub>(r, θ) is the surface brightness at given polar coordinate (r, θ). The coefficients a<sub>k</sub>(r) and b<sub>k</sub>(r) are defined as:

$$ \begin{align} a_k(r) &= \frac{1}{2\pi} \int_0^{2\pi} \mathrm{SB}_{\mathrm{Ly}\alpha}(r, \theta) \cdot \cos(k\theta) \mathrm{d}\theta \nonumber \\ b_k(r) &= \frac{1}{2\pi} \int_0^{2\pi} \mathrm{SB}_{\mathrm{Ly}\alpha}(r, \theta) \cdot \sin(k\theta) \mathrm{d}\theta. \tag{4} \end{align} $$

<footnote n="3">³ We can see that even the absorption uncorrected observed L<sub>Lyα</sub>, smaller symbols in Fig. 9, is on the higher end (Sect. 2.1.1).</footnote>
<figure><img src="image_12.png" /><figcaption>**Fig. 9.** Relation between $Ly\alpha$ nebula luminosity and asymmetry measurement. (a) Flux-weighted $Ly\alpha$ nebula elliptical asymmetry measurement versus nebula luminosity, $L_{Ly\alpha}$. We show the intrinsic flux-weighted ellipticity ($e_{weight}$) in larger open symbols for our targets versus their intrinsic $Ly\alpha$ luminosity. We also show the observed flux-weighted ellipticity $e_{weight,obs}$ for our targets versus their observed $Ly\alpha$ 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 flux-weighted (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\alpha$ nebula elliptical asymmetry measurement versus $L_{Ly\alpha}$. The larger symbol are measurements for our HRGs 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 \to 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\alpha}$ 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.</figcaption></figure>

The detailed description can be found in den Brok et al. (2020). Then we combine $a_k(r)$ and $b_k(r)$:

$$c_k(r) = \sqrt{a_k(r)^2 + b_k(r)^2}. \qquad (5)$$

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 practise, only the first three modes are used since the higher mode coefficients are much smaller (den Brok et al. 2020). Figure 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 HRGs. 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; 4 quasars) and Sanderson et al. (2021) as comparison. The type-1 sample from Borisova et al. (2016; which is quantified in den Brok et al. 2020) is also shown. Our HRGs generally show an increase of the asymmetry as a function of radial distance (some have a smaller secondary peak at $\lesssim 50$ ckpc) and a steep rise to $>1.5$ (7 out of 8) at $\sim 250$ ckpc. The most noticeable exception is MRC0316-257 which has a secondary peak ($\sim 0.8$) at $\sim 50$ ckpc ($\sim 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 HRGs ($\gtrsim 400$ ckpc) are larger than the comparison quasars ($\sim 300$ ckpc). If we limit the comparison to the largest extent ($\sim 250$ ckpc) reached by the type-1s from (Borisova et al. 2016), we find that our HRGs are similar in radial asymmetry measurements. Specifically, both have a relatively flat profile to $\sim 200$ ckpc which followed by a shallow rise to the value of $\sim 0.5$. However, in Sect. 4.3.2 (Fig. 9) we find that the HRGs 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 HRGs 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} \sim 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 HRGs rises at larger radial distance ($\gtrsim 250$ ckpc) and reaches higher asymmetry value ($>1.5$ compared with $\sim 1.4$ of type-2s). At least one of the radial asymmetry measurement (Cdfs 15, with the furthest extent $r_{AGN} \sim 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 HRGs with deeper observations.

### 4.3.4. 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-neh}$ (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 centre 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
<figure><img src="image_13.png" /><figcaption>**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 colour 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).</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**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-neh}$. The typical uncertainty of $d_{AGN-neh}$ ($\sigma_{d_{AGN-neh}}$ = 0.4 pkpc, Table 3) is shown at bottom left. (b) Maximum Ly$\alpha$ nebula linear extent, $d_{max}$, versus $d_{AGN-neh}$. 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-neh}$ and $d_{max}-d_{AGN-neh}$.</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>**Fig. 12.** Distribution of the angle difference between nebular position angle and radio jet position angle, $|\theta - \text{PA}_{\text{radio}}|$. The blue histogram shows the number distribution of flux-weighted angle difference, $|\theta_{\text{weight}} - \text{PA}_{\text{radio}}|$ (sensitive to high surface brightness nebula morphology). The magenta histogram represents the unweighted measurements, $|\theta_{\text{unweight}} - \text{PA}_{\text{radio}}|$ (quantifying the morphology of the whole nebula). The values are presented in Table 3. We mark the obvious outliers in the two distributions in corresponding colours. This shows that most of our HzRGs have their $\text{Ly}\alpha$ nebula elongated along the jet axis.</figcaption></figure>

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*<sub>weight</sub> versus *d*<sub>AGN-neb</sub> and *d*<sub>max</sub> versus *d*<sub>AGN-neb</sub> relations. We note that due to the large uncertainties of *e*<sub>weight</sub> for MRC0316-257, we exclude it from the correlation measurement. From the *r*-values and *p*-values, we can conclude that *e*<sub>weight</sub> and *d*<sub>max</sub> are both correlated with *d*<sub>AGN-neb</sub>. 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) sample for comparison. The (*r*-values, *p*-values) of the *e*<sub>weight</sub> − *d*<sub>AGN-neb</sub> 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*<sub>max</sub> − *d*<sub>AGN-neb</sub> relation for the Arrigoni Battaia et al. (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 revisit these correlations in Sect. 5.

### 4.3.5. Jet-nebula position angle distribution

We present the distribution of the $|\theta - \text{PA}_{\text{radio}}|$ (Sect. 4.3.1, Table 3) in Fig. 12. Both the flux-weighted and unweighted measurements are shown which are sensitive to high-surface brightness parts and the entire nebulae, respectively. Figure 12 shows that most HzRGs have $|\theta - \text{PA}_{\text{radio}}| < 30^\circ$. 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 discuss the indications behind the alignments further in Sect. 5.5. The exception is 4C+04.11 which has both large

<figure><img src="image_16.png" /><figcaption>**Fig. 13.** Distribution of $f_{\text{Ly}\alpha,\text{full FOV}}/\Sigma f_{\text{Ly}\alpha,\text{tile fit}}$ for our HzRGs. The $f_{\text{Ly}\alpha,\text{full FOV}}$ is the intrinsic $\text{Ly}\alpha$ flux resulted from fitting the spectrum summed over the entire FOV. The $\Sigma f_{\text{Ly}\alpha,\text{tile fit}}$ is the summation of the intrinsic $\text{Ly}\alpha$ flux in each tessellation bin (Fig. 1b). The smaller the value, the more likely the $\text{Ly}\alpha$ 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).</figcaption></figure>

$|\theta_{\text{weight}} - \text{PA}_{\text{radio}}|$ (81.4°) and $|\theta_{\text{unweight}} - \text{PA}_{\text{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 $\text{Ly}\alpha$ 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 $|\theta - \text{PA}_{\text{radio}}| < 30^\circ$ is only 0.2% (1.7%) using bimodal distribution.

## 5. Discussion

### 5.1. Absorption correction and radiative transfer

We stress that the analysis of our sampled $\text{Ly}\alpha$ nebulae in this paper is under the idealised assumptions stated in Sect. 3.2.1. Specifically, we interpret the trough features in the $\text{Ly}\alpha$ 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 $\text{Ly}\alpha$ 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 $\text{Ly}\alpha$ photons along the line of sight. The intrinsic $\text{Ly}\alpha$ flux reconstructed under our treatment is thus assumed to be the one after radiative transfer processes have shaped the $\text{Ly}\alpha$ nebula, and is approximated by a Gaussian profile. We also assume the absorbed $\text{Ly}\alpha$ 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 in Fig. C.4 and e.g. Wang et al. 2021b).

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 (Figs. 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, $\Sigma f_{Lyα,tile fit}$). Figure 13 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 (Figs. 2–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 Lyα profile. Our analysis assumes that the absorption correction is the dominant factor.

## 5.2. 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 (orientation 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 ionised gas nebulae but for local AGN and on small scales (sub-kiloparsec to kiloparsec).

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. 12), 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) 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. (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.
### 5.3. 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} \to 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.

Resonant scattering has the potential to shape the observed inner parts of the radial profiles (surface brightness and kinematics) and the morphology correlations. 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 (Figs. 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 positions, that is, 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).

### 5.4. Large-scale environment: Nearby Lyα emission halos

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\sigma$ 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.

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 neighbouring 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 neighbouring 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⁻¹ line width seen at ~100 pkpc (W₈₀, 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 HRGs 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_\star/M_\odot) < 10$ which is approximately 1 – 2 orders of magnitude lower than the stellar masses of our HRGs, 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⁻¹ cm⁻² arcsec⁻² (e.g. their Fig. 7) while we observe it to happen at a level of ~$10^{-18}$ erg s⁻¹ cm⁻² arcsec⁻²). 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 HRGs (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 HRG fields will be reported in a future paper.

## 5.5. Gas distribution on large scales

So far, we have focussed 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 $|\theta_{weight} - PA_{radio}|$ towards values <30°, Fig. 12) through interaction with the gaseous medium. Beyond the extent of the jets, we use $|\theta_{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 \sim 10^{-2}$ cm⁻³), 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 Figs. 11 and 12 is in agreement with such a scenario.

The detection of molecular gas in HRGs 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 HRGs which bridges the radio galaxy and to the large-scale (IGM) gas. If the CGM gas is being accreted onto the central HRGs through these higher density distributions, it would reflect the scenario proposed by Humphrey et al. (2007).

## 5.6. Evidence of capturing inflow from the cosmic web

The relatively large $|\theta - PA_{radio}|$ (81.4° and 45.3° for weighted and un-weighted, respectively) detected in 4C+04.11 is inconsistent with other targets (see Fig. 12). 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.11 resides at the centre 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⁻¹, 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 of the order of $M_{DM} \sim 10^{13} M_\odot$ from its stellar mass, $M_\star \sim 10^{11} M_\odot$ with an empirical $M_\star - M_{DM}$ relation (Wang et al. 2021b). The surface brightness measured in the farthest tile #74 is ~1.2 × $10^{-17}$ erg s⁻¹ cm⁻² arcsec⁻² 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} \sim 10^{12} M_\odot$).

Given the estimated dark matter halo mass, the virial radius and virial velocity can then be calculated as $r_{vir} \sim 100$ kpc and $v_{vir} \sim 580$ km s⁻¹. 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_\odot$ halos at z = 2 and found that the accretion flow structure can remain filamentary within $r_{vir}$ (~0.5 $r_{vir}$). The projected $v_{50}$ (velocity offset) of tile #74 is -536 km s⁻¹ 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 accretion 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 $|\theta - 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 (Sect. 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 are near the detection limit, deeper observation are needed before more conclusive results can be derived.

## 6. 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 an improved sensitivity to low surface brightness emission. We have also 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⁻¹⁹ erg s⁻¹ cm⁻² arcsec⁻².

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 spatially resolved $Lyα$ emission.

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 investigated 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-spot range show a similar gradient as the jet-driven outflow (Fig. 7). The $W_{80}$ profiles show relatively large values (≥1500 km s⁻¹, Fig. 7, Sect. 4.2.4) at all radii and three targets show that 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 that they 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 the 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 for 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 a 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 captured inflows from the 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 scales. 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 kinds of observations will be particularly useful for future simulation works. The MUSE data only reveal 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 of 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").

**Acknowledgements.** We would like to thank Bjorn Emonts and Zheng Cai for valuable discussion which improved this work. This work uses the NASA's Astrophysics Data System and a number of open source software other than the aforementioned ones such as Jupyter notebook (Kluyver et al. 2016); `matplotlib` (Hunter 2007); `SciPy` (Virtanen et al. 2020); `NumPy` (Harris et al. 2020); `Astropy` (Astropy Collaboration 2018); and `LMFIT` (Newville et al. 2016). B.G. acknowledges support from the Carlsberg Foundation Research Grant CF20-0644 'Physical Properties of the Interstellar Medium in Luminous Infrared Galaxies at High redshift: PRISM-LIGHT'. M.V.M. acknowledges support from grant PID2021-124665NB-I00 by the Spanish Ministry of Science and Innovation/State Agency of Research MCIN/AEI/ 10.13039/501100011033 and by "ERDF A way of making Europe". G.N. acknowledges funding support from the Natural Sciences and Engineering Research Council (NSERC) of Canada through a Discovery Grant and Discovery Accelerator Supplement, and from the Canadian Space Agency through grant 18JWST-GTO1. P.L. (contract DL57/2016/CP1364/CT0010) is supported by national funds through Fundação para a Ciência e Tecnologia (FCT) and the Centro de Astrofísica da Universidade do Porto (CAUP).

## References

* Adam, G., Rocca-Volmerange, B., Gerard, S., Ferruit, P., & Bacon, R. 1997, A&A, 326, 501

* Ahumada, R., Prieto, C. A., Almeida, A., et al. 2020, ApJS, 249, 3

* Allen, M. G., Groves, B. A., Dopita, M. A., Sutherland, R. S., & Kewley, L. J. 2008, ApJS, 178, 20

* Andrae, R., Schulze-Hartung, T., & Melchior, P. 2010, ArXiv e-prints [arXiv:1012.3754]

* Antonucci, R. 1993, ARA&A, 31, 473

* Arrigoni Battaia, F., Hennawi, J. F., Cantalupo, S., & Prochaska, J. X. 2016, ApJ, 829, 3

* Arrigoni Battaia, F., Hennawi, J. F., Prochaska, J. X., et al. 2019, MNRAS, 482, 3162
Astropy Collaboration (Price-Whelan, A. M., et al.) 2018, AJ, 156, 123

Bacon, R., Accardo, M., Adjali, L., et al. 2010, SPIE, 7735, 131

Bacon, R., Vernet, J., Borisova, E., et al. 2014, The Messenger, 157, 13

Bacon, R., Piqueras, L., Conseil, S., Richard, J., & Shepherd, M. 2016, Astrophysics Source Code Library [record ascl:1611.003]

Bacon, R., Mary, D., Garel, T., et al. 2021, A&A, 647, A107

Basu-Zych, A., & Scharf, C. 2004, ApJ, 615, L85

Best, P. N., Röttgering, H. J. A., & Longair, M. S. 2000, MNRAS, 311, 23

Bicknell, G. V., Sutherland, R. S., van Breugel, W. J. M., et al. 2000, ApJ, 540, 678

Binette, L., Kurk, J. D., Villar-Martín, M., & Röttgering, H. J. A. 2000, A&A, 356, 23

Bolatto, A. D., Wolfire, M., & Leroy, A. K. 2013, ARA&A, 51, 207

Borisova, E., Cantalupo, S., Lilly, S. J., et al. 2016, ApJ, 831, 39

Byrohl, C., Nelson, D., Behrens, C., et al. 2021, MNRAS, 506, 5129

Cai, Z., Cantalupo, S., Prochaska, J. X., et al. 2019, ApJS, 245, 23

Cantalupo, S., Lilly, S. J., & Haehnelt, M. G. 2012, MNRAS, 425, 1992

Cantalupo, S., Arrigoni-Battaia, F., Prochaska, J. X., Hennawi, J. F., & Madau, P. 2014, Nature, 506, 63

Cappellari, M., & Copin, Y. 2003, MNRAS, 342, 345

Carilli, C. L., Röttgering, H. J. A., van Ojik, R., Miley, G. K., & van Breugel, W. J. M. 1997, ApJS, 109, 1

Chambers, K. C., Miley, G. K., & van Breugel, W. 1987, Nature, 329, 604

Chang, S.-J., Yang, Y., Seon, K.-I., Zabludoff, A., & Lee, H.-W. 2023, ApJ, 945, 100

Christensen, L., Jahnke, K., Wisotzki, L., & Sánchez, S. F. 2006, A&A, 459, 717

Cigan, P. 2019, Astrophysics Source Code Library [record ascl:1909.002]

Cornuault, N., Lehnert, M. D., Boulanger, F., & Guillard, P. 2018, A&A, 610, A75

Costa, T., Arrigoni Battaia, F., Farina, E. P., et al. 2022, MNRAS, 517, 1767

Cowie, L. L., & Hu, E. M. 1998, AJ, 115, 1319

Dawson, S., Rhoads, J. E., Malhotra, S., et al. 2007, ApJ, 671, 1227

De Breuck, C., van Breugel, W., Minniti, D., et al. 1999, A&A, 352, L51

De Breuck, C., van Breugel, W., Stanford, S. A., et al. 2002, AJ, 123, 637

De Breuck, C., Seymour, N., Stern, D., et al. 2010, ApJ, 725, 36

den Brok, J. S., Cantalupo, S., Mackenzie, R., et al. 2020, MNRAS, 495, 1874

Dey, A., Bian, C., Soifer, B. T., et al. 2005, ApJ, 629, 654

Dijkstra, M. 2014, PASA, 31, e040

Dijkstra, M. 2019, in Correction to: Physics of Lyα Radiative Transfer, eds. A. Verhamme, P. North, S. Cantalupo, & H. Atek (Berlin, Heidelberg: Springer)

Drouart, G., De Breuck, C., Vernet, J., et al. 2012, A&A, 548, A45

Drouart, G., De Breuck, C., Vernet, J., et al. 2014, A&A, 566, A53

Duncan, K. J., Windhorst, R. A., Koekemoer, A. M., et al. 2023, MNRAS, 522, 4548

Eales, S. A. 1992, ApJ, 397, 49

Emonts, B. H. C., Norris, R. P., Feain, I., et al. 2014, MNRAS, 438, 2898

ESO CPL Development Team 2015, Astrophysics Source Code Library [record ascl:1504.003]

Evans, I. N., Primini, F. A., Glotfelty, K. J., et al. 2010, ApJS, 189, 37

Falkendal, T., De Breuck, C., Lehnert, M. D., et al. 2019, A&A, 621, A27

Falkendal, T., Lehnert, M. D., Vernet, J., De Breuck, C., & Wang, W. 2021, A&A, 645, A120

Farina, E. P., Arrigoni-Battaia, F., Costa, T., et al. 2019, ApJ, 887, 196

Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306

Fossati, M., Fumagalli, M., Lofthouse, E. K., et al. 2021, MNRAS, 503, 3044

Francis, P. J., Williger, G. M., Collins, N. R., et al. 2001, ApJ, 554, 1001

Gaia Collaboration (Brown, A. G. A., et al.) 2021, A&A, 649, A1

Gronke, M., & Dijkstra, M. 2016, ApJ, 826, 14

Gronke, M., Bull, P., & Dijkstra, M. 2015, ApJ, 812, 123

Gronke, M., Dijkstra, M., McCourt, M., & Oh, S. P. 2016, ApJ, 833, L26

Gullberg, B., De Breuck, C., Lehnert, M. D., et al. 2016, A&A, 586, A124

Hanuschik, R. W. 2003, A&A, 407, 1157

Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357

Harwood, J. J., Hardcastle, M. J., Croston, J. H., & Goodger, J. L. 2013, MNRAS, 435, 3353

He, Z., Sun, A.-L., Zakamska, N. L., et al. 2018, MNRAS, 478, 3614

Heckman, T. M., Lehnert, M. D., Miley, G. K., & van Breugel, W. 1991a, ApJ, 381, 373

Heckman, T. M., Lehnert, M. D., van Breugel, W., & Miley, G. K. 1991b, ApJ, 370, 78

Hennawi, J. F., & Prochaska, J. X. 2013, ApJ, 766, 58

Hippelein, H., & Meisenheimer, K. 1993, Nature, 362, 224

Hu, E. M., & McMahon, R. G. 1996, Nature, 382, 231

Humphrey, A. 2019, MNRAS, 486, 2102

Humphrey, A., Villar-Martín, M., Fosbury, R., Vernet, J., & di Serego Alighieri, S. 2006, MNRAS, 369, 1103

Humphrey, A., Villar-Martín, M., Fosbury, R., et al. 2007, MNRAS, 375, 705

Humphrey, A., Villar-Martín, M., Sánchez, S. F., et al. 2008, MNRAS, 390, 1505

Humphrey, A., Vernet, J., Villar-Martín, M., et al. 2013, ApJ, 768, L3

Hunter, J. D. 2007, Comput. Sci. Eng., 9, 90

Jarvis, M. J., Wilman, R. J., Röttgering, H. J. A., & Binette, L. 2003, MNRAS, 338, 263

Kikuta, S., Imanishi, M., Matsuoka, Y., et al. 2017, ApJ, 841, 128

Kluyver, T., Ragan-Kelley, B., Pérez, F., et al. 2016, in Positioning and Power in Academic Publishing: Players, Agents and Agendas, eds. F. Loizides, & B. Schmidt (IOS Press), 87

Kolwa, S., Vernet, J., De Breuck, C., et al. 2019, A&A, 625, A102

Kolwa, S., De Breuck, C., Vernet, J., et al. 2023, MNRAS, 525, 5831

Krogager, J.-K. 2018, ArXiv e-prints [arXiv:1803.01187]

Lau, M. W., Hamann, F., Gillette, J., et al. 2022, MNRAS, 515, 1624

Leclercq, F., Bacon, R., Wisotzki, L., et al. 2017, A&A, 608, A8

Li, Z., Steidel, C. C., Gronke, M., Chen, Y., & Matsuda, Y. 2022, MNRAS, 513, 3414

Liu, G., Zakamska, N. L., Greene, J. E., Nesvadba, N. P. H., & Liu, X. 2013, MNRAS, 436, 2576

Mackenzie, R., Pezzulli, G., Cantalupo, S., et al. 2021, MNRAS, 502, 494

Marino, R. A., Cantalupo, S., Pezzulli, G., et al. 2019, ApJ, 880, 47

Martin, D. C., Chang, D., Matuszewski, M., et al. 2014, ApJ, 786, 107

Matsuda, Y., Yamada, T., Hayashino, T., et al. 2004, AJ, 128, 569

McCarthy, P. J. 1993, ARA&A, 31, 639

McCarthy, P. J., van Breugel, W., & Kapahi, V. K. 1991, ApJ, 371, 478

Miley, G., & De Breuck, C. 2008, A&ARv, 15, 67

Miley, G. K., Overzier, R. A., Tsvetanov, Z. I., et al. 2004, Nature, 427, 47

Morais, S. G., Humphrey, A., Villar-Martín, M., et al. 2017, MNRAS, 465, 2698

Nelson, D., Genel, S., Pillepich, A., et al. 2016, MNRAS, 460, 2881

Nesvadba, N. P. H., Lehnert, M. D., De Breuck, C., Gilbert, A., & van Breugel, W. 2007, A&A, 475, 145

Nesvadba, N. P. H., Lehnert, M. D., De Breuck, C., Gilbert, A. M., & van Breugel, W. 2008, A&A, 491, 407

Nesvadba, N. P. H., De Breuck, C., Lehnert, M. D., Best, P. N., & Collet, C. 2017a, A&A, 599, A123

Nesvadba, N. P. H., Drouart, G., De Breuck, C., et al. 2017b, A&A, 600, A121

Newville, M., Stensitzki, T., Allen, D. B., et al. 2016, Astrophysics Source Code Library [record ascl:1606.014]

Noirot, G., Vernet, J., De Breuck, C., et al. 2016, ApJ, 830, 90

Noirot, G., Stern, D., Mei, S., et al. 2018, ApJ, 859, 38

Ono, Y., Itoh, R., Shibuya, T., et al. 2021, ApJ, 911, 78

Ouchi, M., Ono, Y., & Shibuya, T. 2020, ARA&A, 58, 617

Parijskij, Y. N., Thomasson, P., Kopylov, A. I., et al. 2014, MNRAS, 439, 2314

Pentericci, L., Röttgering, H. J. A., Miley, G. K., et al. 1999, A&A, 341, 329

Pentericci, L., Van Reeven, W., Carilli, C. L., Röttgering, H. J. A., & Miley, G. K. 2000, A&AS, 145, 121

Prescott, M. K. M., Martin, C. L., & Dey, A. 2015, ApJ, 799, 62

Rauch, M. 1998, ARA&A, 37, 267

Reuland, M., van Breugel, W., Röttgering, H., et al. 2003, ApJ, 592, 755

Reyes, R., Zakamska, N. L., Strauss, M. A., et al. 2008, AJ, 136, 2373

Rosdahl, J., & Blaizot, J. 2012, MNRAS, 423, 344

Röttgering, H. J. A., Hunstead, R. W., Miley, G. K., van Ojik, R., & Wieringa, M. H. 1995, MNRAS, 277, 389

Saito, T., Shimasaku, K., Okamura, S., et al. 2006, ApJ, 648, 54

Sanderson, K. N., Prescott, M. K. M., Christensen, L., Fynbo, J., & Møller, P. 2021, ApJ, 923, 252

Seymour, N., Stern, D., De Breuck, C., et al. 2007, ApJS, 171, 353

Shimasaku, K., Kashikawa, N., Doi, M., et al. 2006, PASJ, 58, 313

Silva, M., Humphrey, A., Lagos, P., et al. 2018a, MNRAS, 481, 1401

Silva, M., Humphrey, A., Lagos, P., et al. 2018b, MNRAS, 474, 3649

Soto, K. T., Lilly, S. J., Bacon, R., Richard, J., & Conseil, S. 2016, MNRAS, 458, 3210

Steidel, C. C., Adelberger, K. L., Shapley, A. E., et al. 2000, ApJ, 532, 170

Swinbank, A. M., Vernet, J. D. R., Smail, I., et al. 2015, MNRAS, 449, 1298

Tepper-García, T. 2006, MNRAS, 369, 2025

Tepper-García, T. 2007, MNRAS, 382, 1375

Tumlinson, J., Peeples, M. S., & Werk, J. K. 2017, ARA&A, 55, 389

Umehata, H., Fumagalli, M., Smail, I., et al. 2019, Science, 366, 97

van Breugel, W. J. M., Stanford, S. A., Spinrad, H., Stern, D., & Graham, J. R. 1998, ApJ, 502, 614

van Ojik, R., Roettgering, H. J. A., Carilli, C. L., et al. 1996, A&A, 313, 25

van Ojik, R., Roettgering, H. J. A., Miley, G. K., & Hunstead, R. W. 1997, A&A, 317, 358

Vayner, A., Zakamska, N. L., Riffel, R. A., et al. 2021, MNRAS, 504, 4445

Vayner, A., Zakamska, N. L., Sabhlok, S., et al. 2023, MNRAS, 519, 961

Veilleux, S., Cecil, G., & Bland-Hawthorn, J. 2005, ARA&A, 43, 769

Venemans, B. P., Röttgering, H. J. A., Miley, G. K., et al. 2007, A&A, 461, 823
Verhamme, A., Schaerer, D., & Maselli, A. 2006, A&A, 460, 397

Vernet, J., Lehnert, M. D., De Breuck, C., et al. 2017, A&A, 602, L6

Villar-Martín, M., Binette, L., & Fosbury, R. A. E. 1996, A&A, 312, 751

Villar-Martín, M., Vernet, J., di Serego Alighieri, S., et al. 2002, MNRAS, 336, 436

Villar-Martín, M., Vernet, J., di Serego Alighieri, S., et al. 2003, MNRAS, 346, 273

Villar-Martín, M., Sánchez, S. F., De Breuck, C., et al. 2006, MNRAS, 366, L1

Villar-Martín, M., Humphrey, A., De Breuck, C., et al. 2007a, MNRAS, 375, 1299

Villar-Martín, M., Sánchez, S. F., Humphrey, A., et al. 2007b, MNRAS, 378, 416

Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nat. Methods, 17, 261

Wang, W., De Breuck, C. A., Lehnert, M. D., et al. 2021a, JWST Proposal. Cycle 1, 1970

Wang, W., Wylezalek, D., De Breuck, C., et al. 2021b, A&A, 654, A88

Weidinger, M., Møller, P., & Fynbo, J. P. U. 2004, Nature, 430, 999

Weidinger, M., Møller, P., Fynbo, J. P. U., & Thomsen, B. 2005, A&A, 436, 825

Weilbacher, P. M., Palsa, R., Streicher, O., et al. 2020, A&A, 641, A28

West, M. J. 1991, ApJ, 379, 19

Wilman, R. J., Jarvis, M. J., Röttgering, H. J. A., & Binette, L. 2004, MNRAS, 351, 1109

Wisotzki, L., Bacon, R., Brinchmann, J., et al. 2018, Nature, 562, 229

Wylezalek, D., Galametz, A., Stern, D., et al. 2013, ApJ, 769, 79

Wylezalek, D., Vernet, J., De Breuck, C., et al. 2014, ApJ, 786, 17

Yang, Y., Zabludoff, A., Tremonti, C., Eisenstein, D., & Davé, R. 2009, ApJ, 693, 1579

Yang, Y., Zabludoff, A., Eisenstein, D., & Davé, R. 2010, ApJ, 719, 1654

Zakamska, N. L., Hamann, F., Pâris, I., et al. 2016, MNRAS, 459, 3144

Zhang, S., Cai, Z., Xu, D., et al. 2023, ApJ, 952, 124

Zirm, A. W., Overzier, R. A., Miley, G. K., et al. 2005, ApJ, 630, 68
# Appendix A: Detection map and tessellation procedures

## A.1. Procedure of detection map for nebulae extent

<table><caption>**Table A.1.** Optimal detection map parameters combination.</caption><thead><tr><th>HzRG</th><th>T<sub>S/N</sub></th><th>n<sub>λ</sub></th><th>2.355σ<sub>max</sub><sup>†</sup></th><th>n<sub>c</sub></th></tr></thead><tbody><tr><td>MRC 0943-242</td><td>3.5</td><td>4</td><td>25</td><td>2</td></tr><tr><td>MRC 0316-257</td><td>3.5</td><td>5</td><td>30</td><td>3</td></tr><tr><td>TN J0205+2242</td><td>3.5</td><td>5</td><td>25</td><td>2</td></tr><tr><td>TN J0121+1320</td><td>2.5</td><td>3</td><td>35</td><td>4</td></tr><tr><td>4C+03.24</td><td>3.5</td><td>5</td><td>20</td><td>2</td></tr><tr><td>4C+19.71</td><td>3.5</td><td>2</td><td>30</td><td>3</td></tr><tr><td>TN J1338-1942</td><td>3.5</td><td>2</td><td>30</td><td>2</td></tr><tr><td>4C+04.11</td><td>3.5</td><td>4</td><td>25</td><td>2</td></tr></tbody></table>

**Notes.** †Maximum Gaussian kernel scale in FWHM. The unit is in pixel unit.

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⁻¹ or −6000 − 6000 km s⁻¹ for 4C+03.24 which 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 first-order 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.355$ pix, 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 positive values), then we flag this spatial location (x, y) as signal-contained (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 would 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) would 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² in size) which are detached from 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.1) as 'optimally smoothed cube'. The smoothed cube are only used in constructing the detection map and not being further used in other analysis. In Appendix A.3 (Fig. A.2-A.9a), we present the smoothed nebula image extracted from the optimally smoothed cube.

## A.2. Tessellation procedure

The Lyα nebulae are irregular in morphology (e.g. Swinbank et al. 2015; Vernet et al. 2017), affected by absorption and having high S/N gradient from centre to outskirt and multiple dynamical components. Hence, neither the simple adoption of Voronoi tessellation (Cappellari & Copin 2003) nor the tessellation method from Swinbank et al. (2015) (which was invented for studying TN J1338-1942 and has also been applied to 4C+04.11 but only works for the central part, see Wang et al. 2021b) can be adopted without adaptation. We describe here the tessellation used in this work.

The Swinbank et al. (2015) tessellation is not optimal for some of our target with significantly irregular nebula shape (e.g. MRC0316-257) since it works with rectangular binning. Hence, we decide to adopt the Voronoi binning (Cappellari & Copin 2003), which is a sophisticated adaptive spatial binning algorithm implemented by various studies of IFU data, for this work. Directly performing Voronoi binning on the image will return detached regions (i.e. regions of spaxels belong to the same tile but are physically separated), as the images are S/N limited for the outer nebulae. The solution could be (i) turning off Weighted Voronoi Tessellation or Centroidal Voronoi Tessellation (see Cappellari & Copin 2003, for details); or (ii) perform-
<figure><img src="image_17.png" /><figcaption>**Fig. A.1.** Schematic presentation of the detection map construction. (a) Average $n_\lambda$ 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 ($\sigma_0 < \sigma_1$) to smooth the spaxels that not passing the $T_{S/N}$, S/N threshold, at each steps until the maximum, $\sigma_{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. (e) Construct the detection map.</figcaption></figure>

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.1 and Appendix A.1) to the continuum-subtracted un-smoothed cube to preserve only the $Ly\alpha$ 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\alpha$ 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\alpha$. 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 $\sigma_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\alpha$ nebulae studied here have high S/N gradient from centre 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² or 0.6×0.4 arcsec²) which is roughly half of the seeing disc. 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.24 which 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² or 0.4 × 0.4 arcsec²).

Finally, for the S/N≤3 spaxels left in the first step, we apply a Voronoi binning with (S/N)$_{vorbin,2nd}$ =8 (11 for 4C+19.71 due to the impact of sky-line).
<figure><img src="image_18.png" /><figcaption>**Fig. A.2.** Smoothed nebula image (a) and tessellation map (b) of MRC0943-242. (a) The purple, yellow and blue colours represent median pseudo-narrow band Lyα images collapsed arbitrarily from red part, middle and blue part, respectively, of the smoothed cubes.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>**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 colour-composed image is created in a manner similar to that of Figure A.2a.</figcaption></figure>

In this way, we obtain the tessellation maps where each tile has S/N> 5. We note that the S/N stated here is calculated based on the narrow band image (see above) collapsed ~ 14 wavelength elements (may vary for different targets) and averaged by the number of wavelength elements. This narrow band image is chosen to enclose as much of the Lyα emission as possible which is a compensation for the different emission peak due to kinematics. Hence, a S/N> 5 in each tile selected in this way is feasible for further fitting as there will be other signal contained wavelength pixels outside the range given the broadness of the emission line. For targets where the extent of the nebula is overlapped with bright foreground stars (TN J1338-1942 and 4C+03.24), we manually assign a *d* = 1 arcsec circular mask at

the position of the star to minimise its impact. For MRC0316-257, we mask the region where a known Lyα emitting galaxies at the similar redshift overlaps in spatial location with the nebula (see Vernet et al. 2017, Arrow galaxy). These manually masked regions are marked by purple hatches in Fig. 1bcd. The tessellation maps are shown in Appendix A.3. Through the tessellation, we reach ~ 2σ surface brightness limit (as reported in Table 2) in the faintest tile.

For the convenience of the spatial fitting (Sect. 3.2.2 and Appendix B), we run an automatic re-numbering algorithm to make physically attached tiles to be as consecutive in number as possible. In each tile, the spatial spectra from every spaxel are then extracted and summed according to the tessellation from the
<figure><img src="image_20.png" /><figcaption>**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 colour-composed image is created in a manner similar to that of Figure A.2a.</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>**Fig. A.5.** Smoothed nebula image (a) and tessellation map (b) of TN J0121+1320. (a) The colour-composed image is created in a manner similar to that of Figure A.2a.</figcaption></figure>

original cube (i.e. the one produced by the reduction, Sect. 2.3, without further continuum subtraction).

### A.3. Smooth nebula images and tessellation maps

In this Appendix, we present the smoothed nebula images (produced based on the optimally smoothed cube, Fig. A.2–A.9a) and the tessellation maps (Fig. A.2–A.9b) constructed following Appendix A.1 (Sect. 3.1.1) and A.2 (Sect. 3.1.2), respectively. The false-colour smooth images are generated using multi-colourfits (Cigan 2019). Each colour represents a pseudo-narrow band Lyα image constructed from the optimally smoothed cube

(Appendix A.1) in arbitrary wavelength range with the goal to show different kinematic structures. Blue, yellow and red are relatively from blue wing, middle and red wing of the Lyα emission, respectively. We note that the smoothed nebula images are only used as representation and demonstration of how the algorithm in Appendix A.1 works. They are not included in the analysis of this work. For MRC0316-257 and TN J0205+2242, we detected line emissions around their Lyα wavelength at isolated regions from the main nebula. We did not include these regions in to account due to its possible origin from companion galaxies, but show the spectra and region where they are extracted in Figs. A.3 and A.4.
<figure><img src="image_22.png" /><figcaption>**Fig. A.6.** Smoothed nebula image (a) and tessellation map (b) of 4C+03.24. (a) The colour-composed image is created in a manner similar to that of Figure A.2a.</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>**Fig. A.7.** Smoothed nebula image (a) and tessellation map (b) of 4C+19.71. (a) The colour-composed image is created in a manner similar to that of Figure A.2a.</figcaption></figure>
<figure><img src="image_24.png" /><figcaption>**Fig. A.8.** Smoothed nebula image (a) and tessellation map (b) of TN J1338-1942. (a) The colour-composed image is created in a manner similar to that of Figure A.2a.</figcaption></figure>

<figure><img src="image_25.png" /><figcaption>**Fig. A.9.** Smoothed nebula image (a) and tessellation map (b) of 4C+04.11. (a) The colour-composed image is created in a manner similar to that of Figure A.2a.</figcaption></figure>
## Appendix B: Fitting procedure

To minimise the uncertainties introduced by fitting the spectra in neighbour bin independently, and impact from foreground targets and the UV continuum, we have developed the following spatial fitting procedure after extracting spectra from the tessellation in Appendix A.2. In this way, we can link the fitting between adjacent bins to avoid un-physical jump in parameters:

To estimate the UV continuum, we firstly extract a narrow band image between the rest-frame Nνλ1240 and OI + Si IIλ1305 from the data cube leaving a rest-frame 10 Å buffer at each end of the central wavelength of the emission line to avoid line emission contamination. This step is implemented to select which tessellated bins are affected by the continuum such that we can model it in the fitting. Hence, this wavelength range is chosen because it is a line-free region closet to Lyα which could be a relatively accurate estimation of the UV continuum emission of the radio galaxy: since the blue wing of Lyα is suffered from absorption. We use this continuum image and the tessellation map (Sect. 3.1.2) to determine which tiles have continuum emissions.

To determine the initial value of the spatial fitting, we secondly extract the 1D Lyα master spectrum from a $d = 1$ arcsec (25 spaxels) at the position of the AGN centre (Table 2) and fit it with Gaussian + Voigt model (Sect. 3.2.1, Fig. 1a). In this way, the initial information of the absorbers (number, column density and positions) are determined. We note that for this fit a 0th-order polynomial for continuum and two Gaussian for both the systemic emission and broad (or blueshifted or redshifted) components are used.

To minimise the impact of undesired features, we thrively mask the strong 5577 Å sky-line, which affects 4C+03.24 and 4C+19.71 most, by 5 spectral resolution units. We also mask the foreground or background objects line emission wavelength range (none of them overlapped with the Lyα of the radio galaxy to our current knowledge) and replace with the noise from the cube variance extension (for example of the reported objects in the vincity of Lyα nebula see, e.g. Falkendal et al. 2021).

To get a first spatial fit, we fourthly fit the spatial spectra from each tiles in the order determined in Appendix A.2 with the Least-square method. In the fitting, we pass the results from the previous tile fit to the next as initials to minimise the potential randomness introduced by fitting each spectra freely without any constrain. We note that (i) the constrains from the previous fit may be relaxed more if the ratio of the integrated observed fluxes between the tiles are significant (~50%) and/or the distance ratio is large (>1.1)<fnref n="4" /> and/or the area of the tile is large (> 9π arcsec²); (ii) the Doppler parameter, *b*, (Tepper-García 2006) is set to a broad range (Kolwa et al. 2019; Wang et al. 2021b, 40 < *b* < 400 km s⁻¹) because of the *b* – *N* degeneracy (e.g. Silva et al. 2018b). We only constrain the range of column density, $N_{H_1}$. The continuum model is included in the fitting with a 0th-order-polynomial for tiles overlapped with continuum detected region determined in first step. A fixed step function with the step at the wavelength of the systemic Lyα is used for the highest redshift target, 4C+04.11, due to the heavily absorbed continuum on the blue-side which may be due to Lyα forest, e.g. Rauch (1998).

Due to the scattering nature of Lyα and presence of broad (outflow) components, it is important to select which region

contains more than one emission component to better study the broader wing. We fit the set of spatial spectra two times: (i) first using 1-Gaussian only for the emission and (ii) then using 2-Gaussian for all spectra. In this way, we can do a simple $\chi^2$ ratio selection between the 2-Gaussian versus the 1-Gaussian fitting results and select the tiles favour 2-Gaussian fit with a threshold of ~0.80–0.98 (depends on different targets). We point that the $\chi^2$ value is not a robust measurement of the fitting quality (e.g. Andrae et al. 2010). For a quick test in our case, however, it is good enough to the first order. We note that for 4C+03.24, we 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 non-single-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 large extent. For example in Fig. 4, we see the tiles at the nebula edge are also affected by absorbers. The column density of the absorbers is a free parameter during the fit. Given the degeneracy between the column density and Doppler parameter, *b* (Silva et al. 2018a), we leave the constraint of *b* to a broad range following Kolwa et al. (2019) (40 to 400 km s⁻¹). The centroid (redshift) of the absorbers are also allowed to vary ($\Delta z \approx \pm 0.001$) except for 4C+04.11⁵. For spectra extracted from inner tiles ($\le 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 $\ge 10''$ such that they can be given low column density (~$10^{13}$ cm⁻²). 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

<footnote n="4">⁴ Distance ratio, $r_1+r_2/d_{1,2}$, is the ratio between the sum of the radii, $r_1 + r_2$, of the two tiles (estimated by $r = \sqrt{A/\pi}$, where A is the area of the tile) and the distance between the geometric centre of the two tiles, $d_{1,2}$.</footnote>

<footnote n="5">⁵ For this case, we follow the procedure in Wang et al. (2021b) and fix the positions of five absorbers on the blue low S/N wing to the value determined from the Master spectrum.</footnote>
over-correction. This may be the case for 4C+04.11 (Fig. 13). We present the flux ratio maps between intrinsic and observed in Appendix C which indicate this possibility (~ 50 for the tiles at the edge).

## Appendix C: Supplementary maps

We show the intrinsic Lyα maps with selected spectra for MRC0943-242, MRC0316-257, 4C+19.71, TNJ1338-1942 and 4C+04.11 in Fig. C.1, C.2, C.3, C.4 and C.5, respectively. The spatial study of the Lyα nebulae for these sources with MUSE were published previously (e.g. Gullberg et al. 2016; Vernet et al. 2017; Falkendal et al. 2021; Swinbank et al. 2015; Wang et al. 2021b). In this paper, we performed a consistent analysis for the full sample with the new method of correcting for absorption. With the individual spectra mostly showing

the tiles at larger distance away from the AGN, we can find that the absorption features are observed nearly across the entire nebula.

To better show the difference before and after absorption correction, we present the observed surface brightness maps and flux ratio maps between the intrinsic flux and observed flux. For the observed flux in each tile, we use the flux integrated from $v_{05}$ to $v_{95}$ and show in Fig. C.6a for each target. For consistency, the $v_{05}$ (and $v_{95}$) is determined based on the intrinsic line since the $v_{05}$ (and $v_{95}$) of the observed lines in low S/N regions are affected by noise more. To ensure the intrinsic and observed surface brightness are comparable as in Fig. C.6b, we take the ratio between the intrinsic and observed flux all integrated between $v_{05}$ and $v_{95}$. Generally, the tiles at the outskirt of the nebulae having larger ratios (> 10). The observed nebular properties derived in the main text are based on the $v_{05}-v_{95}$ observed maps obtained here.

<figure><img src="image_26.png" /><figcaption>Fig. C.1. Similar to Fig. 2–4, but for MRC0943-242.</figcaption></figure>
<figure><img src="image_27.png" /><figcaption>Fig. C.2. Similar to Fig. 2–4, but for MRC0316-257.</figcaption></figure>
<figure><img src="image_28.png" /><figcaption>Fig. C.3. Similar to Fig. 2–4, but for 4C+19.71.</figcaption></figure>
<figure><img src="image_29.png" /><figcaption>Fig. C.4. Similar to Fig. 2–4, but for TNJ1338-1942.</figcaption></figure>
<figure><img src="image_30.png" /><figcaption>Fig. C.5. Similar to Fig. 2–4, but for 4C+04.11.</figcaption></figure>
<figure><img src="image_31.png" /><figcaption>**Fig. C.6.** (a) Observed surface brightness maps in the unit of $10^{-16}$ erg s$^{-1}$ cm$^{-2}$ arcsec$^{-2}$. (b) Flux ratio maps of intrinsic and observed surface brightness. The green cross and contours are the same as Fig. 1 for individual targets, respectively.</figcaption></figure>
## Appendix D: Supplements of nebula radial profile and morphology analysis

<table><caption>**Table D.1.** Position angle and hot spot distance of the radio jet.</caption><thead><tr><th>HzRG</th><th>PA<sub>radio</sub><sup>†</sup><br>deg</th><th>r<sub>app.</sub><sup>*</sup><br>arcsec</th><th>r<sub>rec.</sub><sup>*</sup><br>arcsec</th></tr></thead><tbody><tr><td>MRC 0943-242</td><td>74</td><td>1.79</td><td>1.76</td></tr><tr><td>MRC 0316-257</td><td>53</td><td>5.56</td><td>1.15</td></tr><tr><td>TN J0205+2242</td><td>150</td><td>1.47</td><td>1.24</td></tr><tr><td>TN J0121+1320</td><td>90</td><td>0.82<sup>a</sup></td><td>0.82<sup>a</sup></td></tr><tr><td>4C+03.24</td><td>146</td><td>2.00<sup>b</sup></td><td>1.75<sup>b</sup></td></tr><tr><td>4C+19.71</td><td>176</td><td>3.82</td><td>4.58</td></tr><tr><td>TN J1338-1942</td><td>150</td><td>3.78</td><td>1.18</td></tr><tr><td>4C+04.11</td><td>158</td><td>1.15<sup>c</sup></td><td>1.70<sup>c</sup></td></tr></tbody></table>

**Notes.** †Jet position angle. This is measured east from north. We note that we do not quantify the uncertainties which could be ±1 deg. *Distance between the hot spot and the AGN position. r<sub>app.</sub> and r<sub>rec.</sub> are distance for approaching and receding hot spot, respectively. <sup>a</sup>The radio emission of TN J0121+1320 is spatially unresolved. The ‘hot spot’ distance is represented by the distance from the host galaxy position to the either side of 3σ contour of the radio image along the east-west direction. <sup>b</sup>Approaching hot spot is A1 and receding hot spot is B1 (named after van Ojik et al. 1996). <sup>c</sup>Approaching hot spot is knot8 and receding hot spot is knot1 (named after Parijskij et al. 2014).

In this appendix, we show the supplementary information accompanied with studying nebulae radial profiles and morphology.

We present the annuli used in Sect. 4.2.1 in Fig. D.1. For each targets, the smallest one aperture has the radius equals to 0.75 of its seeing. The radii are in logarithmic steps from centre to outskirt, and the values extract from the annuli are presented in Table D.2.

For the directional radial surface brightness profile analysed in Sect. 4.2.2, it is extracted in half annuli along approaching and receding directions (Fig. 8) with the same step as Fig. D.1.

The position angle of the radio axis is obtained from the two jet hot-spot positions (e.g. van Ojik et al. 1996; Carilli et al. 1997; Pentericci et al. 2000; Parijskij et al. 2014, and unpublished radio maps) and presented in Table D.1. For TN J0121+1320 which has a compact radio emission, we assign an east-west jet position angle to it. We also present the distance of the radio hot spot from the AGN position in Table D.1. This is the value shown in Fig. 6 and 7 (after converted to ckpc). The hot spot is determined to be located at position of the brightest radio emission. 4C+03.24 and 4C+04.11 show multiple radio flux peaks in their radio data. We calculate the jet position angle and hot spot distance based on A1 and B1 knots for 4C+03.24 (named after van Ojik et al. 1996) and knots 1 and 8 for 4C+04.11 (named after Parijskij et al. 2014). We note that the radio jet of TN J0121+1320 is unresolved. Hence, we use the farthest distance reached by the jet (i.e. 3σ contour of the radio image) along the jet direction as a proxy.

We show the exponential and power law fitted parameters (Sect. 4.2.3) of the circularly averaged and directional surface brightness profiles in Table D.3. In Fig. D.2, we show the position-velocity diagram of the observed and intrinsic cubes along the jet direction. This can be used as a direct comparison with Villar-Martín et al. (2003). On each target, we also mark the largest extent of the radio lobes in both directions in horizontal dashed lines. This is determined as the distance from the AGN position along the direction of the jet position angle to the farthest position reached by the 3σ radio flux contour. We note that the broader line width observed visually is due to the contrast between the high and low surface brightness parts. We find in some cases (e.g. 4C+04.11 and TN J1338-1942) that there is a sharp surface brightness drop at the distance of the radio jet boundary. A discontinuity in the line width and surface brightness of the Lyα nebulae across the extent of the radio source has been previously reported in Villar-Martín et al. (2003). The changes is not as sharp as shown in D.2. The sharpness seen in the intrinsic panels is due to tessellation as it is the transition from high to low surface brightness part. Hence, we assure that these are not all artificial effect due to our analysis methods.
<figure><img src="image_32.png" /><figcaption>**Fig. D.1.** Circular apertures for radial profile (Sect. 4.2.1) extraction on top of intrinsic surface brightness maps (Fig. 1a) of the HzRGs sample. In each panel, the blue concentric annuli show the apertures where the radial profile is extracted. The annuli are centred around the host galaxy (central AGN) position (Table 2). The radii of the annuli in each panel are the same in the unit of arcsec for consistency except the smallest radius which is set to be 0.75× seeing for each target. The black bar at the bottom left corner of each panel indicate the 10 arcsec scale.</figcaption></figure>

<table><caption>**Table D.2.** Intrinsic surface brightness from circular aperture</caption><thead><tr><th>R [arcsec]</th><th>MRC0943-242</th><th>MRC0316-257</th><th>TN J0205+2242 (1+z)<sup>4</sup>SB [10<sup>-16</sup> erg s<sup>-1</sup> cm<sup>-2</sup> arcsec<sup>-2</sup>]</th><th>TN J0121+1320</th><th>4C+03.24</th><th>4C+19.71</th><th>TN J1338-1942</th><th>4C+04.11</th></tr></thead><tbody><tr><td>0.0-r<sub>s</sub><sup>†</sup></td><td>2621.4</td><td>1534.1</td><td>1886.6</td><td>1452.9</td><td>2955.1</td><td>261.9</td><td>5177.9</td><td>6464.6</td></tr><tr><td>r<sub>s</sub><sup>†</sup> - 1.0</td><td>1792.0</td><td>833.8</td><td>1764.7</td><td>889.7</td><td>2102.8</td><td>224.9</td><td>3617.2</td><td>4686.2</td></tr><tr><td>1.0-1.8</td><td>623.3</td><td>390.1</td><td>953.4</td><td>250.2</td><td>1884.4</td><td>195.4</td><td>1829.5</td><td>2659.6</td></tr><tr><td>1.8-3.3</td><td>95.3</td><td>60.9</td><td>232.4</td><td>48.5</td><td>1020.3</td><td>130.1</td><td>385.0</td><td>950.4</td></tr><tr><td>3.3-4.5</td><td>24.6</td><td>36.4</td><td>120.3</td><td>38.6</td><td>465.3</td><td>105.4</td><td>106.4</td><td>367.7</td></tr><tr><td>4.5-6.0</td><td>15.2</td><td>31.5</td><td>97.2</td><td>25.4</td><td>227.4</td><td>68.3</td><td>58.4</td><td>229.6</td></tr><tr><td>6.0-7.3</td><td>14.7</td><td>32.2</td><td>87.1</td><td>0.0</td><td>194.8</td><td>51.3</td><td>29.2</td><td>168.9</td></tr><tr><td>7.3-8.9</td><td>13.1</td><td>35.9</td><td>45.0</td><td>0.0</td><td>142.9</td><td>43.0</td><td>25.0</td><td>97.6</td></tr><tr><td>8.9-10.8</td><td>11.9</td><td>38.5</td><td>23.3</td><td>0.0</td><td>117.3</td><td>25.1</td><td>26.3</td><td>65.5</td></tr><tr><td>10.8-13.2</td><td>13.0</td><td>37.2</td><td>23.3</td><td>0.0</td><td>82.8</td><td>20.7</td><td>33.9</td><td>58.3</td></tr><tr><td>13.2-16.0</td><td>0.0</td><td>38.3</td><td>23.3</td><td>0.0</td><td>65.0</td><td>0.0</td><td>33.9</td><td>31.2</td></tr><tr><td>16.0-19.5</td><td>0.0</td><td>33.0</td><td>0.0</td><td>0.0</td><td>45.5</td><td>0.0</td><td>17.0</td><td>0.0</td></tr><tr><td>19.5-23.8</td><td>0.0</td><td>14.5</td><td>0.0</td><td>0.0</td><td>0.0</td><td>0.0</td><td>29.2</td><td>0.0</td></tr><tr><td>23.8-28.9</td><td>0.0</td><td>12.3</td><td>0.0</td><td>0.0</td><td>0.0</td><td>0.0</td><td>36.6</td><td>0.0</td></tr></tbody></table>

**Notes.** $r_s$ = 0.75× seeing which is specified for each target Table 1. The median seeing is 0.75 arcsec. Intrinsic surface brightness values extracted from circular annuli following Fig. D.1. To be consistent, the first column gives the radii of each annulus in the unit of arcsec instead of ckpc. The reported surface brightness values are corrected for cosmological dimming.
<figure><img src="image_33.png" /><figcaption>**Fig. D.2.** Lyα position-velocity diagrams (i.e. 2D spectra) of our sample targets extracted along the radio jet axis. For each target, the left panel shows the diagram constructed from tessellated observed cube (not continuum-subtracted for the host galaxy) and the right panel shows the diagram from absorption-corrected intrinsic cube. The zero offset is set to the position of the central AGN (Table 2). The direction of approaching and receding sides of the jet (see Sect. 2.1.3) are marked in the left panel by the blue and red arrows, respectively. The white horizontal dashed lines represent the furthest extent of the jet. The dotted green contours are given in arbitrary steps which are used to guide the eye for the high brightness part. The vertical black shaded regions in the observed position-velocity diagrams of 4C+03.24 and 4C+19.71 indicate the wavelength ranges affected by the 5577 Å sky-line.</figcaption></figure>
<table><caption>**Table D.3.** Fit parameters of surface brightness profiles.</caption><thead><tr><th>HzRG</th><th>r<sub>h</sub><sup>†</sup><br>pkpc</th><th>C<sub>e</sub><sup>*</sup><br>10<sup>-16</sup> cgs</th><th>r<sub>b</sub><sup>‡</sup><br>pkpc</th><th>α<sup>§</sup></th></tr></thead><tbody><tr><td>MRC 0943-242</td><td>6.9±0.2</td><td>20.7±1.0</td><td>39.8±1.5</td><td>-0.24±0.18</td></tr><tr><td>" " app.</td><td>5.6±0.2</td><td>18.4±1.1</td><td>32.0±1.3</td><td>-0.13±0.15</td></tr><tr><td>" " rec.</td><td>7.5±0.2</td><td>25.8±1.0</td><td>42.9±1.3</td><td>-0.15±0.19</td></tr><tr><td>MRC 0316-257</td><td>7.0±0.6</td><td>8.9±1.2</td><td>27.3±3.2</td><td>-0.62±0.24</td></tr><tr><td>" " app.</td><td>1.3±0.3</td><td>69.0±70.0</td><td>7.7±1.0</td><td>-0.46±0.07</td></tr><tr><td>" " rec.</td><td>7.5±3.3</td><td>14.2±15.7</td><td>31.6±14.4</td><td>-0.01±0.49</td></tr><tr><td>TN J0205+2242</td><td>48.9±30.4</td><td>5.7±0.4</td><td>10.9±21.4</td><td>-1.90±0.58</td></tr><tr><td>" " app.</td><td>8.8±2.0</td><td>6.2±0.7</td><td>13.4±10.8</td><td>-1.52±0.26</td></tr><tr><td>" " rec.</td><td>10.8±3.1</td><td>8.3±0.3</td><td>48.7±0.7</td><td>-0.97±0.14</td></tr><tr><td>TN J0121+1320</td><td>9.8±1.9</td><td>5.6±1.2</td><td>6.2±10.0</td><td>-2.05±0.71</td></tr><tr><td>" " app.</td><td>4.4±1.0</td><td>8.9±2.3</td><td>19.6±6.6</td><td>-0.67±0.76</td></tr><tr><td>" " rec.</td><td>5.1±1.3</td><td>10.2±2.1</td><td>10.9±9.3</td><td>-2.16±0.75</td></tr><tr><td>4C+03.24</td><td>17.6±1.8</td><td>7.9±0.6</td><td>20.0±14.9</td><td>-1.64±0.21</td></tr><tr><td>" " app.</td><td>11.2±0.9</td><td>10.3±0.9</td><td>38.5±4.2</td><td>-0.87±0.18</td></tr><tr><td>" " rec.</td><td>27.7±2.0</td><td>6.5±0.4</td><td>23.4±10.8</td><td>-1.67±0.18</td></tr><tr><td>4C+19.71</td><td>29.5±7.2</td><td>0.8±0.1</td><td>23.1±9.7</td><td>-1.20±0.16</td></tr><tr><td>" " app.</td><td>63.6±20.3</td><td>0.4±0.1</td><td>38.1±7.8</td><td>-1.88±0.37</td></tr><tr><td>" " rec.</td><td>23.3±24.5</td><td>1.0±0.2</td><td>15.1±15.3</td><td>-1.24±0.35</td></tr><tr><td>TN J1338-1942</td><td>7.3±0.3</td><td>13.6±0.7</td><td>41.0±2.7</td><td>-0.12±0.26</td></tr><tr><td>" " app.</td><td>3.6±7.0</td><td>20.7±8.8</td><td>7.5±1.7</td><td>-2.10±0.07</td></tr><tr><td>" " rec.</td><td>7.3±8.0</td><td>19.0±5.1</td><td>40.0±14.7</td><td>-0.52±0.77</td></tr><tr><td>4C+04.11</td><td>8.2±1.5</td><td>11.8±1.3</td><td>15.2±11.7</td><td>-2.04±0.24</td></tr><tr><td>" " app.</td><td>9.9±1.0</td><td>11.0±0.9</td><td>38.9±12.0</td><td>-1.77±0.37</td></tr><tr><td>" " rec.</td><td>7.5±1.6</td><td>11.6±1.6</td><td>16.0±8.3</td><td>-2.14±0.28</td></tr></tbody></table>

**Notes.** See Section 4.2.3 for the fitting equations (Eq. 2): †Scale length of the exponential profile. *Normalisation parameter of the exponential profile. The cgs unit is erg s⁻¹ cm⁻² arcsec⁻². ‡Distance where inner exponential profile changes to power law profile. §Power law index.
**Appendix E: MRC0316-257 systemic redshift**

For MRC0316-257, two velocity components of the He II emission are detected in our MUSE observation which are also separated spatially. The one detected at the position of the X-ray emission peak from *Chandra* observation (Table 2) is believed to be the systemic one while the one at north-west position that is coincident with UV continuum emission peak may trace jet-gas interaction. We show the UV continuum map of MRC0316-257 and the fits of the two He II spectra in Fig. E.1. We note that the UV continuum map is constructed from the MUSE cube using the wavelength in observed frame between N vλ1240 and O I+Si IIλ1305 which is an emission-line-free region of HzRGs (McCarthy 1993). The bright continuum emission object east of the central AGN position peak is a foreground galaxy (Vernet et al. 2017). We report here that the systemic redshift detected is $z_{\text{sys}} = 3.1238 \pm 0.0002$ and the redshift of the component at the UV continuum peak is $z_{\text{red}} = 3.1323 \pm 0.0002$ which is

redshifted of $v = 620 \text{ km s}^{-1}$ from the systemic one. The velocity gradient of the He II agrees with the $Ly\alpha v_{50}$ (Fig. 1c) and [O III]λ5007 (Nesvadba et al. 2008) within the scope of the jet. The UV continuum at this position may suggest the younger stellar population distribution. Combine with the jet kinematics, we may see jet induced star forming activities. However, there is also the possibility that the UV continuum could be produced by the inverse Compton processes. The redshifted He II near the west radio lobe could then be due to the ionisation emission from the shock region exerting on the un-shocked gas. This is supported by the relatively narrow with of the redshift He II (FWHM ~ 600 km s$^{-1}$ < 1000 km s$^{-1}$) which could indicate that it has not been impacted by shocks (Best et al. 2000; Allen et al. 2008). We note that a detailed verification for this scenario is beyond the scope of this paper which involves spectral ageing inspection of the radio jet hot-spot (e.g. Harwood et al. 2013). Hence, we simply point these possibilities and leave them to future study with multi-wavelength observations combined.

<figure><img src="image_34.png" /></figure>

**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 extracted, 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.