11institutetext:  Astronomisches Rechen-Institut, Zentrum für Astronomie der
Universität Heidelberg, Mönchhofstr. 12-14, D-69120 Heidelberg, Germany
11email: wuji.wang@uni-heidelberg.de$;\,\,$11email:
<EMAIL_ADDRESS>22institutetext: European Southern Observatory,
Karl-Schwarzchild-Str. 2, D-85748 Garching, Germany 33institutetext: Cosmic
Dawn Center 44institutetext: DTU Space, Technical University of Denmark,
Elektronvej 327, DK-2800 Lyngby, Denmark 55institutetext: Centre for
Extragalactic Astronomy, Department of Physics, Durham University, South Road,
Durham DH1 3LE, UK 66institutetext: Centro de Astrobiología, CSIC-INTA, Ctra.
de Torrejón a Ajalvir, km 4, 28850 Torrejón de Ardoz, Madrid, Spain
77institutetext: Univ. Lyon, Univ. Lyon1, ENS de Lyon, CNRS, Centre de
Recherche Astrophysique de Lyon UMR5574, 69230 Saint-Genis-Laval, France
88institutetext: International Centre for Radio Astronomy Research, Curtin
University, 1 Turner Avenue, Bentley, WA 6102, Australia 99institutetext: Max-
Planck-Institut fur Astrophysik, Karl-Schwarzschild-Str 1, D-85748 Garching
bei München, Germany 1010institutetext: Instituto de Astrofísica e Ciências do
Espaço, Universidade do Porto, CAUP, Rua das Estrelas, PT4150-762 Porto,
Portugal 1111institutetext: Institute for Computational Astrophysics and
Department of Astronomy & Physics, Saint Mary’s University, 923 Robie Street,
Halifax, NS B3H 3C3, Canada 1212institutetext: Physics Department, University
of Johannesburg, 5 Kingsway Ave, Rossmore, Johannesburg 2092, South Africa

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

Wuji Wang 3D tomography of the giant Ly$\alpha$ nebulae of $z$$\approx$3–5
radio-loud AGN3D tomography of the giant Ly$\alpha$ nebulae of $z$$\approx$3–5
radio-loud AGN Dominika Wylezalek 3D tomography of the giant Ly$\alpha$
nebulae of $z$$\approx$3–5 radio-loud AGN3D tomography of the giant Ly$\alpha$
nebulae of $z$$\approx$3–5 radio-loud AGN Joël Vernet 3D tomography of the
giant Ly$\alpha$ nebulae of $z$$\approx$3–5 radio-loud AGN3D tomography of the
giant Ly$\alpha$ nebulae of $z$$\approx$3–5 radio-loud AGN Carlos De Breuck 3D
tomography of the giant Ly$\alpha$ nebulae of $z$$\approx$3–5 radio-loud AGN3D
tomography of the giant Ly$\alpha$ nebulae of $z$$\approx$3–5 radio-loud AGN
Bitten Gullberg, 3D tomography of the giant Ly$\alpha$ nebulae of
$z$$\approx$3–5 radio-loud AGN3D tomography of the giant Ly$\alpha$ nebulae of
$z$$\approx$3–5 radio-loud AGN3D tomography of the giant Ly$\alpha$ nebulae of
$z$$\approx$3–5 radio-loud AGN3D tomography of the giant Ly$\alpha$ nebulae of
$z$$\approx$3–5 radio-loud AGN Mark Swinbank 3D tomography of the giant
Ly$\alpha$ nebulae of $z$$\approx$3–5 radio-loud AGN3D tomography of the giant
Ly$\alpha$ nebulae of $z$$\approx$3–5 radio-loud AGN Montserrat Villar Martín
3D tomography of the giant Ly$\alpha$ nebulae of $z$$\approx$3–5 radio-loud
AGN3D tomography of the giant Ly$\alpha$ nebulae of $z$$\approx$3–5 radio-loud
AGN Matthew D. Lehnert 3D tomography of the giant Ly$\alpha$ nebulae of
$z$$\approx$3–5 radio-loud AGN3D tomography of the giant Ly$\alpha$ nebulae of
$z$$\approx$3–5 radio-loud AGN Guillaume Drouart 3D tomography of the giant
Ly$\alpha$ nebulae of $z$$\approx$3–5 radio-loud AGN3D tomography of the giant
Ly$\alpha$ nebulae of $z$$\approx$3–5 radio-loud AGN Fabrizio Arrigoni Battaia
3D tomography of the giant Ly$\alpha$ nebulae of $z$$\approx$3–5 radio-loud
AGN3D tomography of the giant Ly$\alpha$ nebulae of $z$$\approx$3–5 radio-loud
AGN Andrew Humphrey 3D tomography of the giant Ly$\alpha$ nebulae of
$z$$\approx$3–5 radio-loud AGN3D tomography of the giant Ly$\alpha$ nebulae of
$z$$\approx$3–5 radio-loud AGN Gaël Noirot 3D tomography of the giant
Ly$\alpha$ nebulae of $z$$\approx$3–5 radio-loud AGN3D tomography of the giant
Ly$\alpha$ nebulae of $z$$\approx$3–5 radio-loud AGN Sthabile Kolwa 3D
tomography of the giant Ly$\alpha$ nebulae of $z$$\approx$3–5 radio-loud AGN3D
tomography of the giant Ly$\alpha$ nebulae of $z$$\approx$3–5 radio-loud AGN
Nick Seymour 3D tomography of the giant Ly$\alpha$ nebulae of $z$$\approx$3–5
radio-loud AGN3D tomography of the giant Ly$\alpha$ nebulae of $z$$\approx$3–5
radio-loud AGN Patricio Lagos 3D tomography of the giant Ly$\alpha$ nebulae of
$z$$\approx$3–5 radio-loud AGN3D tomography of the giant Ly$\alpha$ nebulae of
$z$$\approx$3–5 radio-loud AGN

(Received xxxx; xxxx)

Ly$\alpha$ emission nebulae are ubiquitous around high-redshift galaxies and
are tracers of the gaseous environment on scales out to $\gtrsim 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$\alpha$ nebulae around eight HzRGs ($2.92<z<4.51$).
All the HzRGs have large scale Ly$\alpha$ emission nebulae with seven of them
extended over $100$ pkpc at the observed surface brightness limit ($\sim
2-20\times 10^{-19}\,\rm{erg\,s^{-1}\,cm^{-2}\,arcsec^{-2}}$). Because the
emission line profiles are significantly affected by neutral hydrogen
absorbers across the entire nebulae extent, we perform an absorption
correction to infer maps of the intrinsic Ly$\alpha$ surface brightness,
central velocity and velocity width, all at the last scattering surface of the
observed Ly$\alpha$ photons. We find: (i) The intrinsic surface brightness
radial profiles of our sample can be described by an inner exponential profile
and a power law in the low luminosity extended part; (ii) our HzRGs have
higher surface brightness and more asymmetric nebulae than both radio-loud and
radio-quiet type-1 quasars; (iii) intrinsic nebula kinematics of four HzRGs
show evidence of jet-driven outflows but we find no general trends for the
whole sample; (iv) a relation between the maximum spatial extent of the
Ly$\alpha$ nebula and the projected distance between the AGN and the centroids
of the Ly$\alpha$ nebula; (v) an alignment between radio jet position angles
and the Ly$\alpha$ 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$\alpha$ morphologies. Furthermore, we find evidence showing that the
outskirts of the ionised gas nebulae may be ‘contaminated’ by Ly$\alpha$
photons from nearby emission halos and that the radio jet affects the
morphology and kinematics of the nebulae. Overall, this work provides results
which allow us to compare Ly$\alpha$ 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 atom and molecular hydrogen, $\rm H_{2}$) is the
building block of the baryonic Universe. Studying $\rm H_{2}$ directly is
difficult due to lack of prominent transition lines. It is often probed using
low-$J$ CO transitions as a proxy which unfortunately results in added
uncertainties, for example in the conversion factor (e.g. Bolatto et al.,
2013). In contrast, neutral atomic hydrogen can be easily ionised ($E_{\rm
H^{0}}=13.6\,\rm{eV}$) and cascade with line emissions being produced. The H i
Ly$\alpha\lambda 1216$ (Ly$\alpha$ hereafter) line is the most prominent one
among them. For high-redshift galaxies, it is a commonly targeted emission
line which can easily be observed in the optical to near-infrared bands (e.g.
Hu & McMahon, 1996; Cowie & Hu, 1998; Shimasaku et al., 2006; Dawson et al.,
2007; Leclercq et al., 2017; Wisotzki et al., 2018; Umehata et al., 2019; Ono
et al., 2021; Ouchi et al., 2020, and reference therein). Ly$\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 (e.g. due to resonant nature of Ly$\alpha$ emission and various
potential ionising source acting at once), which is essential to understanding
the physics of the emitting gas observed on different scales and around
various types of objects (Dijkstra, 2019; Ouchi et al., 2020). This is further
complicated when AGN (active galactic nuclei) are present.

Active galaxies hosting AGN, especially the ones with quasar level activities
(bolometric luminosity, $L_{\rm bol}\gtrsim 10^{45}\rm{\,erg\,s^{-1}}$), at
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; Borisova et al., 2016; Cai et al., 2019;
Arrigoni Battaia et al., 2019). The central powerful AGN act as a main
ionising mechanism for the surrounding gas which is responsible for the
detection of these extended Ly$\alpha$ nebulae (as predicted by theoretical
works, e.g. Costa et al., 2022). In addition, the diffuse emission from
galaxies nearby to the AGN host can also contribute to the overall profile
observed of the central target (e.g. Byrohl et al., 2021). In some of the
giant nebulae, it is natural to find various mechanisms functioning at
different scales and positions (e.g. Vernet et al., 2017). Therefore, despite
leaving internal physics entangled, Ly$\alpha$ acts as a simpler tool for
detecting gaseous environment through out cosmic time.

Before wide field integral field spectrographs (IFS) became available, narrow-
band imaging and long slit spectroscopy provided effective methods to detect
diffuse Ly$\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 observations (e.g. Multi-Unit
Spectroscopic Explorer or MUSE/VLT and Keck Cosmic Web Imager or KCWI/Keck)
allow us to measure the extent of the nebulae together with the information of
their dynamics. Numerous works of Ly$\alpha$ nebulae around quasars report
(10s of kpc 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 quasar, ERQ,
Christensen et al., 2006; Borisova et al., 2016; Arrigoni Battaia et al.,
2019; Cai et al., 2019; Farina et al., 2019; den Brok et al., 2020; Fossati et
al., 2021; Mackenzie et al., 2021; Lau et al., 2022; Vayner et al., 2023;
Zhang et al., 2023). This diversity in nebula properties suggest a range of
driving mechanisms, dependencies on orientation, and demonstrate that well-
selected samples are needed. Despite the effort has been done on this topic, a
link between the aforementioned types and type-2 radio-loud quasars on CGM
scale is missing.

Among the high-redshift quasar population, high-redshift radio galaxies
(HzRGs) are a unique sample despite smaller in number (see Miley & De Breuck,
2008, as a review). They host type-2 quasars and have powerful radio jets.
They have been shown to reside in dense protocluster environments (Venemans et
al., 2007; Wylezalek et al., 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 ($\sim 10^{44}\,\rm{erg\,s^{-1}}$,
$\gtrsim 100\,\rm{kpc}$, e.g. Hippelein & Meisenheimer, 1993; van Ojik et al.,
1996, 1997; Reuland et al., 2003; Villar-Martín et al., 2006; Villar-Martín et
al., 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; Villar-Martín et al., 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 (H i absorbers, e.g. Rottgering et al., 1995; van Ojik et al., 1997;
Jarvis et al., 2003; Wilman et al., 2004; Humphrey et al., 2008; Kolwa et al.,
2019).The absorbing gas is found to be extended on galaxy-wides scales and
likely related to outflowing gas from the host galaxy (e.g. Binette et al.,
2000; Swinbank et al., 2015; Silva et al., 2018a; Wang et al., 2021b). The
correction of these absorption is only possible through spectral observation.
Without careful treatment, a considerable amount (a factor of $\gtrsim 5$) of
flux will be missed, and inaccurate conclusions will be derived.
Alternatively, some absorption trough features might potentially be explained
by radiative transfer effects (Dijkstra, 2014; Gronke et al., 2015; Gronke &
Dijkstra, 2016; Gronke et al., 2016). Although it is interesting to compare
the different treatments of the observed Ly$\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 jet direction
(e.g. Chambers et al., 1987; Pentericci et al., 1999; Miley et al., 2004; Zirm
et al., 2005; Duncan et al., 2022) on relatively smaller scale (several kpc to
10s of kpc). Molecular gas detected around HzRGs was reported to be
distributed along the jet within and outside the hot spot which may suggest
several scenarii (e.g. jet-driven outflow, jet-induced gas cooling and jet
propagating into dense molecular gas medium, Emonts et al., 2014; Gullberg et
al., 2016; Falkendal et al., 2021). On Mpc scale, West (1991) found that the
radio jet often points towards nearby galaxies. Eales (1992) proposed a model
explaining the alignment effect, suggesting that the high-redshift radio
emission is often detected when the jet travels close to the major axis of
surrounding asymmetrically distributed gas. With the advanced IFS observation
and 100s kpc gas tracer of Ly$\alpha$, we are able to probe the intrinsic
(i.e. corrected for absorption) gaseous nebula around HzRGs in this work, test
its distribution with respect to the radio jets and seek evidence following
these pioneering works.

In this paper, we utilise the power of MUSE IFU to fully map the Ly$\alpha$
emission nebulae of a HzRGs sample over a redshift range of $2.92-4.51$ and
initiate a comparison with type-1 quasars and study of CGM-scale environments.
We introduce our sample of HzRGs, the MUSE observations and data reduction in
Sect. 2. We present how we measure the maximum extent of the nebulae in Sect.
3.1 and summarise the spectral fitting procedure in Sect. 3.2. We then present
the results of surface brightness, kinematics and morphology in Sect. 4
followed by a discussion in Sect. 5. Finally, we conclude in Sect. 6. In this
paper, we assume a flat $\Lambda$CDM cosmology with
$H_{0}=70\,\rm{km\,s^{-1}\,Mpc^{-1}}$ and $\Omega_{m}=0.3$. Following this
cosmology, $\rm{1\,arcsec\simeq 6.6-7.7\,pkpc}$ for our sample redshifts.
Throughout the paper, pkpc stands for proper kiloparsec and ckpc represents
comoving kiloparsec, ckpc$=(1+z)$pkpc. In this paper, we use ‘intrinsic’ to
refer to the absorption corrected Ly$\alpha$ emission.

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

HzRG | Redshift | UT date | Program ID | Mode | Exp. Time total | seeing
---|---|---|---|---|---|---
|  $z$\rm{{\dagger}}\rm{{\dagger}}$\rm{{\dagger}}$The redshifts are determined from the He ii$\lambda 1640\AA$ 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_{\rm sys}$ in this paper based on He ii fit from our MUSE data (Appendix E). | (dd/mm/yyyy) |  |  | hours | arcsec\rm{*}\rm{*}$\rm{*}$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\,\AA$) produced from the combined cube. We note that the stars used are red in color, i.e. the image quality in the Ly$\alpha$ wavelength should in general be worse than the reported seeing (e.g. larger by 10 to 20%).
MRC 0943-242 | 2.9230 | 21/02/2014 | 60.A-9100(A) | WFM-NOAO-E | 5.21 | 0.65aa$a$The seeing is determined from a star in the overlapping region of the two pointings.
- | - | 15/12/2015-18/01/2016 | 096.B-0752(A) | - | - | -
MRC 0316-257 | 3.1238 | 15/01-17/01/2015 | 094.B-0699(A) | WFM-NOAO-N | 4.24 | 0.61bb$b$There is no available star in the FoV. The seeing is determined from the fit of the most point-like source.
TN J0205+2242 | 3.5060 | 03/12-08/12/2015 | 096.B-0752(B) | WFM-NOAO-N | 4.24 | 0.73
TN J0121+1320 | 3.5190 | 06/10/2015 | 096.B-0752(C) | WFM-NOAO-N | 5.30 | 0.83
- | - | 08/08-28/08/2016 | 097.B-0323(C) | - | - | -
4C+03.24 | 3.5828 | 17/06-18/06/2017 | 60.A-9100(G) | WFM-AO-N | 1.25 | 0.63
4C+19.71 | 3.5892 | 08/06-02/09/2016 | 097.B-0323(B) | WFM-NOAO-N | 5.83 | 1.03
TN J1338-1942 | 4.0959 | 30/04-06/05/2014 | 60.A-9100(B) | WFM-NOAO-N | 8.93 | 0.77aa$a$The seeing is determined from a star in the overlapping region of the two pointings.
- | - | 30/06/2014 | 60.A-9318(A) | - | - | -
4C+04.11 | 4.5077 | 03-15/12/2015 | 096.B-0752(F) | WFM-NOAO-N | 4.24 | 0.88

111

##  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$\alpha$ to be covered by MUSE ; (ii) have
a known extended bright Ly$\alpha$ ($>10\arcsec$) emission nebula; and (iii)
be at DEC $<$ 25∘ to be observable by ground-based telescopes in the southern
hemisphere. This sample also has a wealth of high quality supporting data
obtained by our team, including deep Spitzer/IRAC and Spitzer/MIPS 24 $\mu$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}\rm M_{\odot}$ for
all targets, De Breuck et al., 2010) and the star formation rates ranging from
uppers limit of $<84\,\rm{M_{\odot}\,yr^{-1}}$ to constraints of 626
$\rm{M_{\odot}\,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$\alpha$ nebulae to literature works based on AGN bolometric luminosity.
There are different methods for estimating the bolometric luminosity of AGN,
$L_{\rm{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\,\rm{\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_{\rm{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_{\rm{bol,\,AGN}}$
for our sample. We report that the $L_{\rm{bol,\,AGN}}$ estimated for our
sample using those different methods varies from $10^{45.9}$ to
$10^{48.5}\,\mathrm{erg\,s^{-1}}$. Given this large uncertainty, we find it is
unreasonable to draw further conclusions from the comparison of
$L_{\rm{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_{\rm{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$\lambda 1640\AA$ (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 arcmin2 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\,\rm{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_{\rm{Ly\alpha,\,obs}}=4769\,\AA$). We use
the MUSE commissioning and science verification data of TN J1338-1942 under
the program IDs 60.A9100(B) and 60.A-9318(A) (e.g. Swinbank et al., 2015). For
4C+03.24, we adopt the data released from the MUSE WFM-AO commissioning
observations under the program ID 60.A-9100(G). The information of the
observations of our sample, in the order of redshift, is summarised in Table
1. For each object, observations consist of 1 (4C+03.24) to 6 (TN J 1338-1942)
observing blocks (OBs). Within each OB, the 2 or 3 exposures of 20$-$30
minutes were slightly dithered (with a $<1$″ amplitude pattern) and rotated by
90 degrees from each other.

###  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 et al., 2021). Two sources had no GAIA star
within the MUSE field-of-view (FoV). For TN J0121+1320 we use the SDSS DR16
(16th Data Release, Ahumada et al., 2020) catalogue instead. For MRC 0316-257,
we use Gaia EDR3 to first correct the astrometry of the HST/ACS F814W image
and then matched the MUSE cube to the HST image.

Using ZAP directly for sky-line residual removal without applying masks may
remove faint narrow Ly$\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 region222We 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.; (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\simeq 10\arcsec$
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}\,\rm{erg\,s^{-1}\,cm^{-2}\,\AA^{-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.

Table 2: HzRGs MUSE sample properties.

HzRG | RA (J2000)\rm{{\dagger}}\rm{{\dagger}}$\rm{{\dagger}}$Position of the AGN and/or host galaxy. | DEC (J2000)\rm{{\dagger}}\rm{{\dagger}}$\rm{{\dagger}}$Position of the AGN and/or host galaxy. | SB limit\rm{*}\rm{*}$\rm{*}$The surface brightness limit is the 2$\sigma$ limit extracted from continuum-source- and Ly$\alpha$-free regions in a narrow band image. The narrow band image is collapsed from $v_{05}$ to $v_{95}$ (see Section 4.1). The cgs unit is $\rm{erg\,s^{-1}\,cm^{-2}\,arcsec^{-2}}$. |  $L_{\rm{Ly\alpha,\,int}}$\rm{{\ddagger}}\rm{{\ddagger}}$\rm{{\ddagger}}$Intrinsic Ly$\alpha$ 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). |  $L_{\rm{Ly\alpha,\,obs}}$\rm{$\mathsection$}\rm{$\mathsection$}$\rm{$\mathsection$}$Observed Ly$\alpha$ luminosity of the nebula integrated over the entire area from $v_{05}$ to $v_{95}$ (see text).
---|---|---|---|---|---
| hh:mm:ss | dd:mm:ss |  $10^{-19}$ cgs | $10^{44}\,\rm{erg\,s^{-1}}$ | $10^{44}\,\rm{erg\,s^{-1}}$
MRC 0943-242 | 09:45:32.73 |  $-$24:28:49.65 aa$a$Spitzer (Seymour et al., 2007; De Breuck et al., 2010). | 17.0 | 4.7$\pm 0.1$ | 2.39$\pm$0.03
MRC 0316-257 | 03:18:12.07 |  $-$25:35:10.22 bb$b$Chandra Obs. ID 5734 (PI:Pentericci). Chandra Source Catalog (Evans et al., 2010). | 2.62 | 7.3$\pm 1.2$ | 1.04$\pm$0.04
TN J0205+2242 | 02:05:10.69 |  $+$22:42:50.4 aa$a$Spitzer (Seymour et al., 2007; De Breuck et al., 2010). cc$c$Radio (De Breuck et al., 2002). | 13.3 | 8.3$\pm 0.4$ | 3.82$\pm$0.06
TN J0121+1320 | 01:21:42.73 |  $+$13:20:58.0 aa$a$Spitzer (Seymour et al., 2007; De Breuck et al., 2010). | 5.64 | 2.0$\pm 0.1$ | 0.55$\pm$0.01
4C+03.24 | 12:45:38.37 |  $+$03:23:21.0 dd$d$Radio (van Ojik et al., 1996). | 11.5 | 28.8$\pm 0.8$ | 5.43$\pm$0.20
4C+19.71 | 21:44:07.56 |  $+$19:29:14.6 aa$a$Spitzer (Seymour et al., 2007; De Breuck et al., 2010). | 4.77 | 4.3$\pm 0.3$ | 1.48$\pm$0.22
TN J1338-1942 | 13:38:26.10 |  $-$19:42:31.1 ee$e$Radio (De Breuck et al., 1999). | 4.34 | 12.3$\pm 0.3$ | 5.89$\pm$0.10
4C+04.11 | 03:11:47.97 |  $+$05:08:03.74 ff$f$Radio (Parijskij et al., 2014). | 9.84 | 20.0$\pm 0.7$ | 2.89$\pm$0.07
333
References of the positions:

Figure 1: Mapping results of our MUSE HzRGs sample. (a) Master Ly$\alpha$
spectrum (blue shaded histogram) extracted from a $r=0.5\,\mathrm{arcsec}$
aperture at the AGN position with best fit (solid dark magenta line). Red
dashed curve shows the intrinsic Ly$\alpha$ 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\AA$ sky-line. The flux density unit, $F_{\lambda}$, is
$10^{-20}\,\mathrm{erg\,s^{-1}\,cm^{-2}\,\AA^{-1}}$. We also show the scaled
He ii$\lambda 1640\AA$ 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$\alpha$
spectrum in $-1000$ to $1000$ $\mathrm{km\,s^{-1}}$. The $\Delta
v=0\,\mathrm{km\,s^{-1}}$ is the systemic redshift based on He ii or [C i]
(Table 1, Kolwa et al., 2023). (b) Intrinsic Ly$\alpha$ 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\sigma$ contour along the east-west direction. The white
hatched regions are the ones where the flux uncertainty is higher than 50% of
the fitted intrinsic flux. The white bar indicates the 50 pkpc at the redshift
of the radio galaxy. The unit of the surface brightness is
$10^{-16}\,\mathrm{erg\,s^{-1}\,cm^{-2}\,arcsec^{-2}}$. We apply the same
color scale for all targets. (c) $v_{50}$ map of the intrinsic Ly$\alpha$
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_{80}$ map of the intrinsic Ly$\alpha$ 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 filed of MRC0316-257, see Vernet
et al., 2017). We note that the color scales for panel (c) and (d) are
customised. The purple hatched area (if any) indicates the manually excluded
region affected by foreground star or known Ly$\alpha$ emitter. Figure 1:
Continued. Figure 2: Example for the intrinsic mapping of the Ly$\alpha$
nebula of TNJ0205+2242. The central panel shows the intrinsic surface
brightness map of TNJ0205+2242 which is the same as Fig. 1b. The green cross
and triangles mark the position of the AGN and jet lobes, respectively. In
each of the side panel, we show the spectrum (blue shade histogram in
normalised flux unit) extracted from the individual spatial bin whose number
is labeled at the top left, and the best fit (dark magenta curve) and
recovered intrinsic Ly$\alpha$ (dashed red line). The black vertical bars
indicate the positions of the H i absorbers. Figure 3: Similar figure as Fig.
2 for TNJ0121+1320. Figure 4: Similar figure as Fig. 2 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\,\AA$ sky-line. Figure 5: Radial profiles of
the Ly$\alpha$ nebulae extracted in circular annuli (Fig.29). For better
comparison, we show the radial profile in comoving kpc (ckpc) and take the
cosmological dimming into account by a factor of $(1+z)^{4}$, where $z$ is the
redshift of the target. The black dot-dashed curve and grey dotted line in
both panels are the best fitted exponential and power law profiles of the
Arrigoni Battaia et al. (2019) radio loud sample, respectively. The two
vertical dotted lines mark the 50 and 300 ckpc, respectively. Upper panel:
Radial profile of observed surface brightness map in thicker solid lines. In
this panel, we also include the radial profiles of other quasar samples
(dashed lines) for comparison: B16 — Borisova et al. (2016), AB19 — Arrigoni
Battaia et al. (2019), C19 — Cai et al. (2019), Farina19 — Farina et al.
(2019), Fossati21 — Fossati et al. (2021), L22 — Lau et al. (2022) and V22 —
Vayner et al. (2023) (two sources, 4C09.17 and 7C 1354+2552). When it is
available, we show the range spanned by the 25th and 75th of the comparison
sample radial profile as the shaded region around median profile in the same
color. The horizontal bars at the right-most indicate the observed surface
brightness limits (scaled by area from Table 2) for each target in the same
color. Lower panel: Intrinsic radial profile in thicker solid lines. The
shaded regions around each profiles indicates the uncertainty range of the
surface brightness from fitting. In this panel, we show again the same
profiles of F19 and V22 7C as in the upper panel for comparison. Our HzRGs are
extended further with higher surface brightnesses (or flattening in some
sources) at larger radii ($\sim 300$ ckpc) compared to other samples. Figure
6: Surface brightness radial profiles for approaching (blue squares) and
receding (red circles) directions along the jet axis. The dotted curves in
corresponding colors show the exponential$+$power law fits for the two
directional profiles. We also include the fits for the circularly averaged
profile in solid magenta lines. In each panel, the magenta shaded region mark
again the same uncertainty range for the intrinsic surface brightness profile
as Fig. 5. The solid green curve is the normalised radial profile of a star
extracted up to $2\arcsec$ (the one in the FoV of MRC0316-257 is extracted
from a round galaxy due to no available star) showing the PSF (Table 1). The
vertical dashed lines indicate the distances of the jet hot spots in
corresponding colors. The profile along the receding side of the jet is
brighter than along approaching side for most sources within the extent of the
jets except 4C+03.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
$\gtrsim 100$ ckpc for MRC0943-242, MRC0316-257 and TNJ1338-1942 which may
related to nearby companions (see Sect. 5.4). Figure 7: Directional $v_{50}$
and $W_{80}$ profiles for approaching (blue squares) and receding (red
circles) sides along the jet axis extracted from the intrinsic maps (i.e.
corrected for absorption, Fig. 1cd). The vertical dashed lines indicate the
distances of the jet hot spots (blue for approaching, red for receding,
Appendix D). We note that the radio emission of TN J0121+1320 is unresolved.
The grey shaded regions are $>$5 arcsec from the host galaxy. The data points
in the shaded regions should be treated with caution given the large tile size
may smear kinematic structures. The horizontal black dotted line in the
$v_{50}$ panel marks the $0\,\mathrm{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_{\rm H\beta}\simeq
3.566$, $-1100\,\mathrm{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.

##  3 Data analysis

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

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

####  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 SNR threshold ($T_{\mathrm{S/N}}$) and
leaving to the next iterations only spaxels below this S/N threshold, until a
maximum smoothing radius is reached ($\sigma_{\mathrm{max}}$). The spaxels not
selected by the end of the iterative process are masked out. To further clean
the smoothed data cube from spurious noise features and make sure that a
proper line fitting can be made, we mask spatial positions selected by the
adaptive smoothing algorithm in less than $n_{\rm c}$ consecutive wavelength
bin.

To determine the optimal combination of the 4 parameters ($n_{\lambda}$,
$T_{\mathrm{S/N}}$, $\sigma_{\mathrm{max}}$ and $n_{\rm 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 ($\sim 15\,\AA$ 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 modeling

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 ($\sim$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 $\sim 30$ higher resolution than MUSE display sharp edges which
is fully consistent with a well-defined absorption profile rather than
radiative transfer effects. We note that the the term ‘radiative transfer’
used in the paper refers to the process where Ly$\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 Fig. 2, 3 and 4
where similar absorption features are seen in the selected spectra at larger
distance (10s of kpc) away from AGN. Similar maps of the remaining sources are
shown in Appendix C which are the ones have been previously published
(Swinbank et al., 2015; Gullberg et al., 2016; Vernet et al., 2017; Falkendal
et al., 2021; Wang et 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 modeling works (e.g.
Arrigoni Battaia et al., 2019; Chang et al., 2022). This could be a result of
prior radiative transfer effects of Ly$\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
assumptions that absorption rather than radiative transfer is primarily
responsible for the line profiles. We will discuss the limitation of this
treatment in Sect. 5.1.

We note that the ${O\textsc{v}}]\lambda\lambda 1213.8,1218.2$ (O v]) line
underneath the Ly$\alpha$ can affect the obtained flux especially in the
nuclear region where the ionisation parameter (and metallicity) is higher
(Humphrey, 2019). In our pilot study (Wang et al., 2021b) of 4C+04.11, we
found the contribution from O v] is negligible. Hence, we do not further
include O v] in our line fitting. We leave the inspection to future work when
data of metal lines (e.g. N v$\lambda\lambda 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$\alpha$ emission across the nebula, we fit
each spectrum in each tessellation bin (see Sect. 3.1.2) following the
procedure described in Sect. 3.2.1. We take into account the physical
connection between neighbouring tiles by using the fit results of a previous
connected bin as the starting parameters for the next bin (see Appendix A for
details on the ordering). We determine the number of absorbers based on the
Master spectrum where the S/N is the highest (Fig. 1a). We then use that same
number of absorbers across the nebula, where the centroid, column density and
Doppler parameter of absorbers are fitted in a given range (Appendix B). This
assumption is supported by the profile shapes at the largest spatial extents
(see Fig. 2, 3, and 4 and also Appendix C). We note that the number of
absorbers selected here may be incomplete but this has minor effects on the
results of this paper: (i) the absorbers that impact most the intrinsic flux
(i.e. spatially extended $\sim 10\arcsec$ and having higher column density
and/or larger Doppler parameter) are included; (ii) absorbers that seem to be
‘superfluous’ at the wings have only minor effects on the reconstructed flux
where S/N is low (Fig. 2, 3 and 4). Future work using high spectral resolution
data will address these issues also taking into account that some of these
absorbers have counterparts in metal lines covered by the MUSE data (e.g. N
v$\lambda\lambda 1238,1243$ and C IV$\lambda\lambda 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 $\chi^{2}$. The fit is done using a Least-squares method
followed by a Markov chain Monte Carlo sampling (MCMC, emcee Foreman-Mackey et
al., 2013). The uncertainties we report are either the direct output of the
$1\sigma$ error by the MCMC or the propagated $1\sigma$ 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
sources we show the Master spectrum together with its best fit in Fig. 1a as
an example (Sect. 3.2.1). We also show the non-resonant He ii spectrum
extracted in the same aperture (green histogram, not continuum subtracted)
which is used for systemic redshift ($\Delta v=0\,\mathrm{km\,s^{-1}}$, Table
1) determination. We note that there is no He ii detected at the AGN position
for 4C+03.24 (Sect. 4.2.4). In addition, to illustrate how fitting procedure
works spatially (Sect. 3.2.2), the selected exemplar individual fits are shown
in Fig. 2, 3 and 4 (also see Appendix C).

The intrinsic Ly$\alpha$ surface brightness maps are shown in Fig. 1b on the
same flux scale. Regions with larger fitting uncertainties ($\gtrsim 50\%$)
that should be treated with caution are indicated by the overlaid hatched
tiles. We report the total intrinsic Ly$\alpha$ luminosities
($L_{\mathrm{Ly\alpha,int}}$) of the nebulae in Table 2 and their maximum
linear extent, $d_{\rm 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 $\sim 347$ pkpc (MRC0316-257). TNJ0121+1320 is the only target with
nebula $<100$ pkpc ($\sim 72$ pkpc). The total intrinsic surface brightness
($L_{\mathrm{Ly\alpha,int}}$) of the nebulae ranges from $2$ to $29\times
10^{44}\,\mathrm{erg\,s^{-1}}$.

To characterise the kinematic information of the intrinsic nebulae that are
fitted with one or two Gaussians, we use a set of non-parametric emission line
measurements (see e.g. Liu et al., 2013) derived from the cumulative line flux
as a function of velocity $\Phi(v)$ defined as:

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

where $f(v^{\prime})$ is the flux density at $v^{\prime}$. 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$\alpha$ 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, for example the
southwest tile of 4C+04.11, displaying larger $W_{80}$ values ($\gtrsim
2500\,\mathrm{km\,s^{-2}}$). Except 4C+19.71, all targets show a line width of
$\sim 800-2500\,\rm{km\,s^{-1}}$. For 4C+19.71, due to the strong $5577\,\AA$
sky-line located close with the observed Ly$\alpha$ peak wavelength, its line
width should be treated as lower limit ($\gtrsim 600\mathrm{km\,s^{-1}}$)
especially for the tiles in the outskirts of the nebula.

We note that the non-parametric measurements used in this mapping are based on
intrinsic (= absorption-corrected) line profiles which are determined through
model fitting, same as Wang et al. (2021b). In Appendix C, we present the maps
of observed surface brightness and flux ratio as supplementary material.

###  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$\alpha$ 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. 29. 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 6 lists the extracted intrinsic profile
values.

Figure 5 shows the radial profiles after correction for cosmological dimming
and in comoving units for observed (upper panel) and intrinsic (lower panel,
corrected for absorption) maps. The dashed lines in the upper panel represent
the comparison quasar samples or single targets (ERQ and 2 radio-loud quasars
from Lau et al., 2022; Vayner et al., 2023, respectively). 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
preformed without absorption corrections (Sect. 5.1). Nevertheless, since the
comparison quasar samples are not corrected for absorption, we do not show
them all again in our intrinsic profile (lower) panel. The two exceptions are
Farina et al. (2019) (F19) and Vayner et al. (2023) (7C 1354+2552, V22 7C).
Those two are on the higher surface brightness end of the comparison samples,
and we will examine them quantitatively along with both the intrinsic and
observed HzRGs profiles. We note that Vayner et al. (2023) fitted the
Ly$\alpha$ absorbers from the spatially integrated 1D spectrum and found $\sim
10^{13.5}\,\mathrm{cm^{-2}}$ for the column densities. For Farina et al.
(2019), 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 ERQ from Lau et al. (2022) (which is also highly obscured),
type-1 radial profiles are dominated by direct emission from bright AGN point
source in the inner regions ($\sim$50 ckpc or $\sim$10 pkpc). Hence, due to
point spread function (PSF) subtraction, the inner-most radius covered in the
comparison samples is limited to $\sim$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
$\sim$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 $\sim 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 $\sim$50 ckpc, the intrinsic HzRG profiles are still a factor of
$1.1-15$ (or $4-40$ at $\sim 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 $\gtrsim 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_{\rm 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 $\sim 1\times
10^{-15}\,\rm{erg\,s^{-1}\,cm^{-2}\,arcsec^{-2}}$ (MRC0316-257, corrected for
cosmological dimming) which is higher than the faint end of the quasar samples
by a factor of $5-40$ (not at similar comoving distances). These indicate that
we are observing some of the most extend Ly$\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 $\sim 2-4$
($\sim 2-6$) at radii of $\sim$400 ($\sim$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 ($\gtrsim 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 $\pm
2000\,\rm{km\,s^{-1}}$ in Cai et al. (2019), $30\,\AA$ in Arrigoni Battaia et
al. (2019) and $\pm 500\,\rm{km\,s^{-1}}$ in Farina et al. (2019). Considering
the redshift difference between these samples, the integration range adopted
are consistent. For our study, particularly for the observed radial profile,
our extraction is based on the $v_{05}$ and $v_{95}$ which are determined
based on intrinsic fitting (Sect. 4.1). In this way, we can minimise the
uncertainties coming from the observed line width difference, for example
between the emission lines in the vicinity of the host and outskirts of the
nebula. Our velocity range ($v_{05}$ and $v_{95}$) used is basically the value
of $W_{90}$ which has the range of $\sim 800-2700\,\rm{km\,s^{-1}}$ for all
tiles of all targets. Nevertheless, we conduct a check by extracting observed
circular radial profiles through integration of $30\,\AA$ around the systemic
redshift of our targets for comparison. The results vary by $\sim\pm 10\%$ in
each annulus to the profiles in Fig. 5 (from $v_{05}$ and $v_{95}$),
especially for emissions at $>50\,\rm{ckpc}$ where the line width is narrower
comparing to the centre. The $30\,\AA$ 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$\alpha$ 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$\alpha$ 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$\alpha$ 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$\alpha$ 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.444We 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.
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

$\text{SB}(r)=\left\\{\begin{array}[]{lr}C_{\rm e}\mathrm{exp}(-r/r_{\rm
h})&\,r\leq r_{\rm b}\\\ C_{\rm p}r^{\alpha}&\,r>r_{\rm
b}\end{array}\right.\text{,}$ (2)

where $r_{\rm h}$ is the scale length of the exponential profile, $r_{\rm b}$
is the distance at which the inner and outer profiles separate and $C_{\rm e}$
and $C_{\rm p}$ are normalisation parameters for exponential and power law
profiles, respectively ($C_{\rm p}=C_{\rm e}\mathrm{exp(-r_{\rm b}}/r_{\rm
h})/r^{\alpha}_{\rm b}$). The determination of the piecewise function is
motivated by previous studies of quasar Ly$\alpha$ 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 7 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 $\sim 1$ dex lower than the receding side. This
could be an extreme case of uneven Ly$\alpha$ 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_{\rm b}$ (nor $r_{\rm h}$). As Sect. 4.2.1
described, our HzRGs are high in surface brightness (large $C_{\rm e}$); the
reasons for this include (i) our sample is composed of HzRGs with bright
Ly$\alpha$ 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_{\rm h}$ values derived for our sample
are mostly $<20$ pkpc (Table 7) which is consistent with the quasars. This
suggests a similarity between the central (high surface brightness) part of
HzRGs to other quasars (type-1 radio-loud and radio-quiet, type-2 radio-
quiet), despite the high surface brightness in our sample. We note that the
PSF-subtraction of quasar samples and resolution effects will impact the inner
part to the profile. We will further discuss the power law declining
(flattening) part of our nebula in Sect. 5 combining the information from
nebular morphology (Sect. 4.3).

####  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$\alpha$ 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$\alpha$ is an unreliable tracer of kinematics at
least on 10s to $100$s pkpc scale in HzRGs. We note that the tessellation
implemented, especially for tiles with larger size ($\sim
5\,\mathrm{arcsec^{2}}$) which are usually located in the low S/N region away
from the host galaxy, may smear out potential kinematic features. One possible
consequence of combining different kinematic components is broadening of the
line width. This may be the case for the receding side of MRC0316-257 and both
sides of TNJ0121+1320 and 4C+04.11. In general, the $W_{80}$ does not show an
increasing trend toward larger $r_{\rm 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_{\rm AGN}>5\arcsec$ 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$\alpha$ emission gas centred around
systemic velocity, at least in the vicinity of the AGN. Offsets of the
$v_{50}$ levels at $r_{\rm AGN}\sim 0$ ckpc from $0\,\mathrm{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$\alpha$ photons. The most noticeable case is
4C+03.24 which has an offset of $\sim 900\,\mathrm{km\,s^{-1}}$. We note that
its systemic redshift (Table 1) is based on [C i](1-0) emission (Kolwa et al.,
2023) due to lack of He ii from the AGN position in the MUSE data (Fig. 1a).
This offset can be eased if we use the redshift of H$\beta$, $z_{\rm
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\,\mathrm{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_{\rm 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$\alpha$ 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 ($\sim
1500\,\mathrm{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 $\gtrsim 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$\alpha$ 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_{\rm AGN}>5\arcsec$) should be treated with caution. We note that the FWHM
in Villar-Martín et al. (2003) was derived without correction for absorption.
Since we see a high spatial coverage of absorbers (Fig. 2, 3 and 4), the
correction indeed helps with recovering the close-to-intrinsic gas kinematics
at large radii. In Fig. 30, we show the position-velocity diagram extracted
based on our observed and intrinsic surface brightness maps along the jet as a
direct comparison with the long slit spectroscopic study. Although it
resembles the detection of Villar-Martín et al. (2003) at first glance, we
note that this is due to the tessellation and the contrast between high
surface brightness and low surface brightness part.

By considering both the radial profiles from $v_{50}$ and $W_{80}$, we can
generally find that the profiles within the jet extents have behave
differently compared to the profiles outside the jet extent. This again
suggests the jet is disturbing or entraining the Ly$\alpha$ emitting gas.
There are unclear signs of kinematics other than outflows or inflows seen
mostly at larger radial distance $\sim 300$ ckpc. For example, MRC0943-242
stays relatively flat (for both $v_{50}$ and $W_{80}$), while MRC0316-257 has
a decrease in $W_{80}$ followed by an increase beyond the jet extent on the
receding side. We reiterate that in this analysis we do not distinguish
between (potential) different kinematics components by using $v_{50}$ and
$W_{80}$ to quantify the overall velocity shift and line width. This may bring
bias of the measured values. Additionally, the measured kinematics farther
from the AGN are averaged from larger annulus (e.g. projected area of $\sim
4\times 10^{4}$ pkpc2 at $\sim 60$ pkpc or $\sim 300$ ckpc) which will bring
another bias. We point out again that we use grey shade to mark the data
$>5\arcsec$ from the center which has larger tile size. Nevertheless, we note
that the detected $W_{80}$ of $\sim 10^{3}\,\mathrm{km\,s^{-1}}$ (and abrupt
velocity shift) at large radii ($\sim 300$ ckpc) in some of the profiles could
be caused by the fact that the detected Ly$\alpha$ emission is dominated by
emission halos of nearby companions (e.g. Byrohl et al., 2021).

Figure 8: Zoom-in intrinsic surface brightness maps (Fig. 1b) of the HzRGs
sample to $15\times 15\,\rm{arcsec^{2}}$ (or $\sim 110\times
110\,\rm{pkpc^{2}}$) around the central AGN (blue cross). In each panel, the
blue$+$red and green solid lines indicate the direction of, and perpendicular
to, the radio jet. The blue (red) color represents the direction of the
approaching (receding) jet. The white dashed line shows the flux-weighted
position angle of the nebula ($\theta_{\rm weight}$). The white dotted line
shows the unweighted position angle of the nebula ($\theta_{\rm unweight}$).
The white star indicates the intrinsic flux-weighted centroid of the
Ly$\alpha$ 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. Figure 9: (a) Flux-weighted
Ly$\alpha$ nebula elliptical asymmetry measurement versus nebula luminosity,
$L_{\rm Ly\alpha}$. We show the intrinsic flux-weighted ellipticity
($e_{\rm{weight}}$) in larger open symbols for our targets versus their
intrinsic Ly$\alpha$ luminosity. We also show the observed flux-weighted
ellipticity $e_{\rm{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_{\rm{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_{\rm Ly\alpha}$. The larger symbol are
measurements for our HzRGs while the grey symbols are comparison targets
(type-1s: B16 — Borisova et al. (2016), L22 — Lau et al. (2022) and M21 —
Mackenzie et al. (2021); type-2s: dB20 — den Brok et al. (2020) and S21
—Sanderson et al. (2021)). We mark the median $e_{\rm{unweight}}$ for type-1s
(0.69) and type-2s (0.80) in solid and dashed horizontal lines, separately.
The $e_{\rm{weight}}$ is sensitive to the morphology of the high surface
brightness part of the nebula while the $e_{\rm{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_{\rm 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.  Figure
10: Radial profiles of the surface brightness Fourier decomposition (asymmetry
measurement). The $c_{0}$, $c_{1}$ and $c_{2}$ are the 0th, 1st and 2nd modes
Fourier decomposition coefficients of the surface brightness radial profile,
respectively (see den Brok et al., 2020, for definition). The
$(c_{1}/c_{0})^{2}+(c_{2}/c_{0})^{2}$ is a measurement of nebula asymmetry
along the radial distance from the AGN. Our HzRGs are shown in solid color
lines. $r_{\rm 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 ($\sim 25$ ckpc)
range affected by median seeing of our sample (the radial distance where the
type-1 PSF is affected is $\sim 50$ ckpc, see den Brok et al., 2020). The
morphologies for most of the HzRGs nebulae are round (symmetric) $\lesssim
100$ ckpc and become asymmetric at larger radial distances $\sim 100$ ckpc
(see text).  Figure 11: (a) Intrinsic flux-weighted ellipticity (sensitive to
the morphology of the high surface brightness part of the nebula), $e_{\rm
weight}$, versus distance offset between nebular centroid and AGN position,
$d_{\rm AGN-neb}$. The typical uncertainty of $d_{\rm AGN-neb}$ ($\sigma_{\rm
d_{\rm AGN-neb}}=0.4$ pkpc, Table 3) is shown at bottom left. (b) Maximum
Ly$\alpha$ nebula linear extent, $d_{\rm max}$, versus $d_{\rm AGN-neb}$. In
both panels, we give the Spearman correlation measurements for our sample at
the top left (the star superscript indicates the correlation is calculated
excluding MRC 0316-257 in panel (a)). We also include the data from Arrigoni
Battaia et al. (2019) in both panels in grey (radio quiet) and pink (radio
loud) dots. We note that the Arrigoni Battaia et al. (2019) sample shown in
this figure is incomplete to concentrate on the relation for our targets.
There are positive correlations detected for $e_{\rm weight}$-$d_{\rm AGN-
neb}$ and $d_{\rm max}$-$d_{\rm AGN-neb}$. Figure 12: Distribution of
$f_{\mathrm{Ly\alpha,full\,FOV}}/\Sigma f_{\mathrm{Ly\alpha,tile\,fit}}$ for
our HzRGs. The $f_{\mathrm{Ly\alpha,full\,FOV}}$ is the intrinsic Ly$\alpha$
flux resulted from fitting the spectrum summed over the entire FOV. The
$\Sigma f_{\mathrm{Ly\alpha,tile\,fit}}$ is the summation of the intrinsic
Ly$\alpha$ flux in each tessellation bin (Fig. 1b). The smaller the value, the
more likely the 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).  Figure 13: Distribution of the angle
difference between nebular position angle and radio jet position angle,
$|\theta-\rm{PA_{radio}}|$. The blue histogram shows the number distribution
of flux-weighted angle difference, $|\theta_{\rm{weight}}-\rm{PA_{radio}}|$
(sensitive to high surface brightness nebula morphology). The magenta
histogram represents the unweighted measurements,
$|\theta_{\rm{unweight}}-\rm{PA_{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 colors. This shows that
most of our HzRGs have their Ly$\alpha$ nebula elongated along the jet axis.
Table 3: HzRGs nebulae properties.

HzRG |  $e_{\rm{weight}}$\rm{{\dagger}}\rm{{\dagger}}$\rm{{\dagger}}$Intrinsic flux-weighted ($e_{\rm{weight}}$), observed flux-weighted ($e_{\rm{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.  |  $e_{\rm{weight,\,obs}}$\rm{{\dagger}}\rm{{\dagger}}$\rm{{\dagger}}$Intrinsic flux-weighted ($e_{\rm{weight}}$), observed flux-weighted ($e_{\rm{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.  |  $e_{\rm{unweight}}$\rm{{\dagger}}\rm{{\dagger}}$\rm{{\dagger}}$Intrinsic flux-weighted ($e_{\rm{weight}}$), observed flux-weighted ($e_{\rm{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.  |  $|\theta_{\rm{w.}}-\rm{PA_{radio}}|$\rm{*}\rm{*}$\rm{*}$Absolute difference between flux-weighted (unweighted) nebula position angle, $\theta_{\rm{weight}}$ ($\theta_{\rm{unweight}}$), and radio axis position angle, $\rm{PA_{radio}}$. |  $|\theta_{\rm{unw.}}-\rm{PA_{radio}}|$\rm{*}\rm{*}$\rm{*}$Absolute difference between flux-weighted (unweighted) nebula position angle, $\theta_{\rm{weight}}$ ($\theta_{\rm{unweight}}$), and radio axis position angle, $\rm{PA_{radio}}$. |  $d_{\rm AGN-neb}$\rm{{\ddagger}}\rm{{\ddagger}}$\rm{{\ddagger}}$Offset between the intrinsic flux-weighted centroid of Ly$\alpha$ nebula and AGN position. |  $d_{\rm max}$\rm{$\mathsection$}\rm{$\mathsection$}$\rm{$\mathsection$}$Maximum linear extent spanned by the Ly$\alpha$ nebula.
---|---|---|---|---|---|---|---
|  |  |  | deg | deg | pkpc | pkpc
MRC 0943-242 | 0.73$\pm$0.01 | 0.76 | 0.62 | 4.7$\pm$1.8 | 7.2 | 3.4$\pm$0.1 | 164
MRC 0316-257 | 0.28$\pm$0.38 | 0.96 | 0.87 | 28$\pm$54 | 20.8 | 13.5$\pm$0.2 | 347
TN J0205+2242 | 0.91$\pm$0.01 | 0.79 | 0.75 | 2.1$\pm$2.3 | 8.6 | 9.0$\pm$0.1 | 170
TN J0121+1320 | 0.66$\pm$0.05 | 0.74 | 0.48 | 10.5$\pm$8.7 | 72.1 | 3.1$\pm$0.1 | 72
4C+03.24 | 0.79$\pm$0.01 | 0.72 | 0.82 | 8.1$\pm$1.2 | 2.0 | 9.9$\pm$0.1 | 248
4C+19.71 | 0.85$\pm$0.02 | 0.83 | 0.51 | 20.1$\pm$9.1 | 29.3 | 6.7$\pm$0.1 | 143
TN J1338-1942 | 0.94$\pm$0.01 | 0.87 | 0.90 | 1.3$\pm$1.4 | 1.2 | 15.3$\pm$0.1 | 265
4C+04.11 | 0.76$\pm$0.01 | 0.20 | 0.59 | 81.4$\pm$1.3 | 45.3 | 5.7$\pm$0.1 | 163

555

###  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$\alpha$ nebulae (in our sample) can probe the
CGM gas beyond 100 pkpc. By visual inspection, we observe that the shape of
our Ly$\alpha$ 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$\alpha$ 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_{\rm{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_{\rm{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_{\rm{unweight}}$ following den
Brok et al. (2020). In this context, we also derive $\theta_{\rm{unweight}}$
(flux-unweighted position angle) to examine the jet-nebula relation with
respect to the entire nebula. Fig. 8 visualises the weighted (nebula centroids
and $\theta_{\rm{weight}}$) and unweighted ($\theta_{\rm{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_{\rm{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_{\rm weight}$ (and $\theta_{\rm
unweight}$) are irrelevant, we report the projected distance between the
nebula centroid and the AGN position ($d_{\rm AGN-neb}$) and the difference in
angles between $\theta_{\rm weight}$ (and $\theta_{\rm unweight}$) and the jet
position angle ($|\theta_{\rm weight}-\rm{PA_{radio}}|$ and $|\theta_{\rm
unweight}-\rm{PA_{radio}}|$), respectively. The jet position angle
($\rm{PA_{radio}}$) is shown in Fig. 8 and listed in Appendix D.

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

Fig. 9 presents the ellipticity measurements as a function of their nebula
Ly$\alpha$ 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_{\rm 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_{\rm{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_{\rm{weight}}\approx 0.72$, dashed horizontal
in Fig. 9a). The $e_{\rm{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_{\rm{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-257 is more towards the asymmetric end. The large
difference between its intrinsic and observed $e_{\rm{weight}}$ could be due
to the absorption correction elevates the flux difference between the high and
low S/N regions thus gives more weight to the central nebula. As for 4C+04.11,
this could be due to the potential over-correction of the absorption in the
low S/N regions given we use nine absorbers across the nebula (Sect. 3.2,
5.1). However, we point again that the absorption is necessary to reconstruct
the intrinsic flux given that the absorption features across the nebula were
observed (e.g. Swinbank et al., 2015; Wang et al., 2021b).

Our HzRGs have a median $e_{\rm{unweight}}\approx 0.70$ which is lower than
the measurement of type-2s and relatively consistent with the type-1s (median
$e_{\rm{unweight}}$ $\approx 0.80$ and $\approx 0.69$, respectively, Borisova
et al., 2016; den Brok et al., 2020; Sanderson et al., 2021; Mackenzie et al.,
2021). The ellipicity of the comparison type-2s is calculated from five
sources which may not be representative. We note that there is a large scatter
in the $e_{\rm{unweight}}$ measured for our sample with four clustered around
the type-2s and other four below the median $e_{\rm{unweight}}$ of type-1s.
This result could be biased by the implemented analysis methods (i.e.
detection map and tessellation, Sect. 3.1): (i) the construction of the
detection map may smooth out the nebula asymmetry at larger radii (lower S/N);
(ii) the tessellation results in lager bin sizes along the direction
perpendicular to the radio jet (lower surface brightness and S/N, Sect. 5.5)
which shapes the nebula morphology to be more round. This brings more effects
of the quantification of the entire nebula. Hence, the HzRGs nebulae will
likely have larger $e_{\rm{unweight}}$ (i.e. more asymmetric) than the
quantified value. In spite of that, we can conclude that at least half of our
sample, together with the type-2s, distribute at the most asymmetric end in
terms of the whole nebula. The rest has the possibility to be skewed to higher
$e_{\rm{unweight}}$. Our HzRGs have diverse properties and are not necessarily
representative for the entire HzRGs population. Inspection for individual
targets is required.

As for the $L_{\rm Ly\alpha}$, our HzRGs host the most luminous Ly$\alpha$
nebula compared with other quasar types. This is due to: (i) sample
selection666We can see that even the absorption uncorrected observed $L_{\rm
Ly\alpha}$, smaller symbols in Fig. 9, is on the higher end (Sect. 2.1.1).;
(ii) absorption correction; (iii) quasar PSF subtraction of comparison samples
(Sect 4.2.3).

We further compare our HzRGs to the ERQ from Lau et al. (2022). In Fig. 9,
both the $e_{\rm{weight}}$ (0.44) and $e_{\rm{unweight}}$ (0.69) of L22 are
lower than the measurements from HzRGs. Its $L_{\rm Ly\alpha}$ is also lower
than our HzRGs by $\sim 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_{k}(r)$ and
$b_{k}(r)$:

$\mathrm{SB_{Ly\alpha}}(r,\theta)=\sum_{k=0}^{\infty}\left[a_{k}(r)\cdot\cos(k\theta)+b_{k}(r)\cdot\sin(k\theta)\right]\text{,}$
(3)

where $\mathrm{SB_{Ly\alpha}}(r,\theta)$ is the surface brightness at given
polar coordinate $(r,\theta)$. The coefficients $a_{k}(r)$ and $b_{k}(r)$ are
defined as:

$\begin{split}a_{k}(r)&=\frac{1}{2\pi}\int^{2\pi}_{0}\mathrm{SB_{Ly\alpha}}(r,\theta)\cdot\cos(k\theta)d\theta\\\
b_{k}(r)&=\frac{1}{2\pi}\int^{2\pi}_{0}\mathrm{SB_{Ly\alpha}}(r,\theta)\cdot\sin(k\theta)d\theta\text{.}\end{split}$
(4)

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}}.$ (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 $k$th and 0th
modes, $c_{k}(r)/c_{0}(r)$. In practice, only the first three modes are used
since the higher mode coefficients are much smaller (den Brok et al., 2020).
Fig. 10 shows the radial profiles of this asymmetry measurement (represented
by $(c_{1}/c_{0})^{2}+(c_{2}/c_{0})^{2}$ on the y-axis) for our HzRGs. The
larger the value, the more asymmetric the morphology shift from circular shape
at a given radial distance. We also include type-2 measurements from den Brok
et al. (2020) (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 HzRGs 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 HzRGs ($\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 HzRGs 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 HzRGs ellipticity is larger
than type-1s. Together with the radial asymmetry profile in this section, we
can conclude that the asymmetry of the nebulae associated with HzRGs and
type-1s are different due to structures at larger radial distance ($>250$
ckpc) from AGN. This is also suggested by den Brok et al. (2020). Although we
caution the large $W_{80}$ at large distance in Sect 4.2.4 and mark the region
with $r>5$ arcsec in grey in Fig. 7, we can still find that most targets have
at least $W_{80}\sim 10^{3}\,\mathrm{km\,s^{-1}}.$ This may indicate that
these structures are likely disturbed and not ‘quiescent’ gas in the Cosmic
Web (Sect. 5). From the comparison with type-2s (den Brok et al., 2020;
Sanderson et al., 2021, which also have a steep rise), we find that the radial
asymmetry of our HzRGs rises at larger radial distance ($\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_{\rm 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 HzRGs with deeper observations.

####  4.3.4 Morphological correlations

In Fig. 11, we compare the $e_{\rm weight}$ (ellipticity that is sensitive to
high surface brightness nebula) and $d_{\rm max}$ (maximum nebula extent)
against $d_{\rm AGN-neb}$ (offset between AGN and nebula flux-weighted
centroid). For a consistent comparison for the unweighted measurements, we
should use the unweighted offset (centroid) and the unweighted ellipticity. We
note that the calculation of unweighted ellipticities, $e_{\rm unweight}$,
assumes that the centre of nebula corresponds to the AGN position (den Brok et
al., 2020). It would involve large uncertainties if we calculated the flux-
unweighted nebula centroids (i.e. the geometric center of the nebula) which
are entirely depended on the spaxel distribution from our detection maps
(Sect. 3.1.1). Hence, we report only the correlation between flux-weighted
ellipticity and offset, and use it as a proxy of the nebula even though they
are more sensitive to the high surface brightness part.

The (r-value, p-value) of Spearman’s correlation coefficients are (0.89,
0.007) and (0.88, 0.004) for the $e_{\rm weight}\,vs.\,d_{\rm AGN-neb}$ and
$d_{\rm max}\,vs.\,d_{\rm AGN-neb}$ relations. We note that due to the large
uncertainties of $e_{\rm weight}$ for MRC0316-257, we exclude it from the
correlation measurement. From the r-values and p-values, we can conclude that
$e_{\rm weight}$ and $d_{\rm max}$ are both correlated with $d_{\rm AGN-neb}$.
This suggests that the more asymmetric and more extended nebulae have larger
offsets between AGN position and nebular centroid. We also include the
Arrigoni Battaia et al. (2019) sample for comparison. The (r-values, p-values)
of the $e_{\rm weight}-d_{\rm AGN-neb}$ relation for the Arrigoni Battaia et
al. (2019) whole sample and radio-loud targets are (0.3, 0.02) and (0.3, 0.4),
respectively, implying no strong correlations. The (r-values, p-values) of the
$d_{\rm max}-d_{\rm AGN-neb}$ 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 will revisit these
correlations in Sect. 5.

####  4.3.5 Jet-nebula position angle distribution

We present the distribution of the $|\theta-\rm{PA_{radio}}|$ (Sect. 4.3.1,
Table 3) in Fig. 13. Both the flux-weighted and unweighted measurements are
shown which are sensitive to high-surface brightness parts and the entire
nebulae, respectively. Fig. 13 shows that most HzRGs have
$|\theta-\rm{PA_{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 will discuss the
indications behind the alignments further in 5.5. The exception is 4C+04.11
which has both large $|\theta_{\rm weight}-\rm{PA_{radio}}|$ ($81.4^{\circ}$)
and $|\theta_{\rm{unweight}}-\rm{PA_{radio}}|$ ($45.3^{\circ}$), 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^{\circ}$.
Given its round 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^{\circ}$, the probability for us to find 7 (or 6)
targets in a sample of 8 with $|\theta-\rm{PA_{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 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 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 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 Ly$\alpha$ photons along the line of sight.
The intrinsic Ly$\alpha$ flux reconstructed under our treatment is thus
assumed to be the one after radiative transfer processes have shaped the
Ly$\alpha$ nebula, and is approximated by a Gaussian profile. We also assume
the absorbed 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$\alpha$ 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. 26 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$\alpha$ photons (i.e. central AGN in
our case). This leaves the trough (or ‘main absorber’) at around the systemic
redshift of the AGN as predicted by radiative transfer studies (e.g. Verhamme
et al., 2006). Even though the scales probed in those simulations are
different (CGM scale versus sub-kpc), this might be the case in some of our
targets, for example MRC0943-242 and TNJ0121+1320 (Fig. 3 and 23). Therefore,
one consequence of our absorption treatment (with ‘lost’ photons) is that we
may double-count Ly$\alpha$ 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_{\mathrm{Ly\alpha,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_{\mathrm{Ly\alpha,tile\,fit}}$). Fig. 12 shows the result, where a value of
1.0 is expected if our treatment is fully accurate. Instead, we observe a
median value of 0.71 for our sample, indicating that there may be double-
counting of photons. However, several effects may also contribute to this.
First, our assumption of a Gaussian shape for the intrinsic profile may not be
accurate, as prior radiative transfer effects may have created troughs and
boosted the line wings (Verhamme et al., 2006). A future paper presenting high
spectral resolution UVES observations of our sample will help to better model
the profiles (Ritter et al. in prep.). Second, the assumption that all
absorbers extend across the entire nebula may also be oversimplified. While
many absorbers are indeed seen on 10s (or 100s) of kpc scales for our targets
(Fig. 2, 3 and 4 and also Appendix C), there may be exceptions.

The Ly$\alpha$ nebula properties of the type-1 (and radio-quiet type-2)
sources in the comparison samples are derived without correcting for
absorption. They are still good comparison samples given that not many
absorption features are detected in them (e.g. Arrigoni Battaia et al., 2019).
This fact alone is already an indication that there may be an environmental
differences between the nebulae of HzRGs and other quasars. Alternatively,
this difference may related to intrinsic differences between different AGN
types.

Overall, it is likely that both the absorption and radiative transfer effect
are working together shaping the the Ly$\alpha$ profile. Our analysis assumes
that the absorption correction is the dominant factor.

###  5.2 Ly$\alpha$ nebula and AGN unification model

The unification model of AGN (e.g. Antonucci, 1993) proposed that the
fundamental difference between type-1 and type-2 is due to the different
orientation of the obscuring dusty torus (ionisation cone). Specifically, the
ionisation cone of type-1s is more aligned with the observed line of sight
than in type-2s.

Within this unification model, we assume that the power of AGN is on a similar
level for type-1s and type-2s and that their gaseous nebulae therefore have
similar distributions and masses. If we further assume that the ionising
photons are primarily produced by the AGN, then the nebulae should have
similar morphologies and luminosities. In this picture, the Ly$\alpha$ nebulae
are elongated along the direction of the ionisation cone and have a ‘rugby-
ball’ shape. For type-1s, the ionisation cone would be pointing towards the
observer resulting in a rounder nebula morphology. For type-2s, the ionisation
cone would aligned along the plane of the sky resulting in the observed
elliptical morphology. Evidence for such a scenario was indeed reported in He
et al. (2018) using ionized gas nebulae but for local AGN and on small scales
(sub-kpc to kpc).

Using both the $e_{\rm weight}$ (sensitive to the high surface brightness
nebula) and $e_{\rm 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$\alpha$ nebula
(at least in the high surface brightness part where the photons of AGN are
presumably dominating the ionisation, Sect. 4.3.5 and Fig.13), we argue that
the AGN orientation between type-1s and type-2s (HzRGs) can explain the
observed morphological difference. The evidence for this claim comes mostly
from the observed ellipticity distribution (Fig. 9): the $e_{\rm 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_{\rm weight}$ (0.66) and $e_{\rm 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$\alpha$ 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 $\gtrsim 40$ pkpc. Using the same radial asymmetry measurement
in Fig.10 (Sect. 4.3.3), we also find large nebula asymmetries at $\gtrsim 40$
pkpc ($\sim 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 $\lesssim 20-30$ pkpc ($\sim 100$ ckpc)
could be a reason for the observed high $e_{\rm 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 $\sim
60$ ckpc in Fig. 10.

###  5.3 The role of Ly$\alpha$ resonant scattering

As presented in Sect. 4.1 and 4.3, Ly$\alpha$ nebulae around our HzRGs are
extended in size ($\gtrsim 150$ pkpc) and are asymmetric in shape.
Interestingly, there is a correlation between the maximum nebula extent
$d_{\rm max}$ ($e_{\rm weight}$: ellipticity measurement that is sensitive to
high surface brightness nebula) and the offset between AGN position and
nebulae’s flux-weighted centroid $d_{\rm 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$\alpha$ 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$\alpha$
photons, regardless of the powering source, could result in an offset between
the luminosity centroid of the nebula and the quasar position ($d_{\rm AGN-
neb}$). This offset can vary depending on the line of sight and can reach
$\sim 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_{\rm 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$\alpha$ 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_{\rm weight}$ and $d_{\rm AGN-neb}$. Specifically, when the
scattering is more efficient, we may observe the nebula to be more asymmetric
and more extended, at least in the high surface brightness regime.

The inner part (luminous) of the surface brightness radial profiles have an
exponential shape (Sect. 4.2.1). This gradual change in surface brightness and
profile slope at smaller radii is suggested to be due to scattering (Costa et
al., 2022). The high $W_{80}$ values out to radius of $\sim 50$ pkpc (or $\sim
230$ ckpc, Sect. 4.2.4) could also be related with efficient scattering (Fig.
B1 in Costa et al., 2022). The velocity shift of the nebulae can also be
complicated due to scattering at the similar radial distance (Fig. B1 in Costa
et al., 2022) which is the case of our $v_{50}$ (Fig. 1c, Sect. 4.1). As shown
by the stellar radial profiles (Fig. 6), the impact of seeing cannot be fully
responsible for the exponential shape at smaller radii. For the observed
kinematics, the jet (which is not included in the simulations) may also play a
role here. We overlay the jet hot spot distances on the radial profiles (Fig.
6 and 7). Qualitatively, the $v_{50}$ and $W_{80}$ show different behaviours
within and outside the extent of the jet hot spots at least for some targets
(e.g. MRC0316-257). We note that the distances marked by the vertical lines
are measured from the brightest radio emission positions, i.e. the full extent
of jet is even larger (Appendix D). The decrease of $W_{80}$ at large radii
($\gtrsim 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$\alpha$ properties, Costa et al. (2022) furthermore predicts that
Ly$\alpha$ nebulae of edge-on AGN should have lower surface brightness, larger
$d_{\rm 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_{\rm
max}-d_{\rm 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_{\rm 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_{\rm max}$ and $d_{\rm AGN-neb}$ despite the expectation that most
HzRG are observed under similar, large viewing angles, (i.e. close to edge-
on). All of our targets have clear two-sided radio lobe morphology (except TN
J0121+1320) and none of them shows clear signs of broad-line emission or
significant contamination by a bright point-like source in the centre (see
also Drouart et al., 2012, and Sect. 2.1.3).

Resonant scattering has the potential to shape the observed inner parts of the
radial profiles (surface brightness and kinematics) and the morphology
correlations.

###  5.4 Large scale environment: nearby Ly$\alpha$ emission halos

Byrohl et al. (2021) studied Ly$\alpha$ emission halos and their relation with
environmental factors using TNG50 simulation. They found a flattening in the
Ly$\alpha$ halo surface brightness radial profile at large radial distances
($\gtrsim 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$\alpha$ 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$\alpha$
photons. The filamentary Ly$\alpha$ 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_{\rm max}-d_{\rm AGN-neb}$ in Sect. 4.3.4. When a quasar resides in a
dense environment, its apparent size will likely be impacted by neighboring
Ly$\alpha$ 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_{\rm max}-d_{\rm AGN-neb}$ ($e_{\rm weight}-d_{\rm AGN-neb}$, Sect. 4.3.4,
Fig. 11) relations may be related to the distribution of the companion
emission halos. Contamination from neighboring halos will result in the
centroid of the nebula being offset from the AGN position and cause a more
asymmetric nebula morphology. The nearby companions can be the disturbing
sources resulting in the $\sim 1000\leavevmode\nobreak\ \mathrm{km,s^{-1}}$
line width seen at $\sim 100$ pkpc ($W_{80}$, Fig. 7).

The radial asymmetry profiles of our targets have a larger value than type-1s
and type-2s (Fig. 10) at larger radii, which is unlikely to be entirely due to
the AGN orientation (Sect. 4.3.3). This can now be explained by the impact and
contamination of nearby companions at 100s of pkcp. We point out that the
type-2 from Sanderson et al. (2021) has a projected size of $173$ pkpc, a
nebula centre offsets from the AGN of $\sim 17$ pkpc and high ellipticity
(0.8). This source is also found to have nearby companions ($\sim 60$ pkpc).
The difference of the radial asymmetry profiles (Fig. 10) and surface
brightness profiles between this source and our HzRGs may be related to the
jet (or large scale gas distribution, Sect. 5.5).

The stellar masses of the galaxies studied in Byrohl et al. (2021) are in the
range of $8.5<\log_{10}(M_{\star}/\mathrm{M_{\odot}})<10$ which is
approximately $1-2$ orders of magnitude lower than the stellar masses of our
HzRGs, implying lower dark matter halo masses, as well. Such lower masses may
explain the difference in the level of Ly$\alpha$ 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 $\sim
10^{-20}\,\mathrm{erg\,s^{-1}\,cm^{-2}\,arcsec^{-2}}$ (e.g. their Fig. 7)
while we observe it to happen at a level of $\sim
10^{-18}\,\mathrm{erg\,s^{-1}\,cm^{-2}\,arcsec^{-2}}$). The dark matter halo
mass difference can also explain the different radial distance at which the
companion emission dominates ($\sim 30$ pkpc in Byrohl et al., 2021, versus
$\sim 40$ pkpc, Table 7). We note that we use the $r_{\rm 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 $\gtrsim 100$ pkpc extents of the Ly$\alpha$ and the dense environments of
HzRGs (Wylezalek et al., 2013) suggest that nearby halos may contribute to our
observations (Byrohl et al., 2021). However, this also suggests that our CGM
study is ‘contaminated’ by sources other than the radio galaxy. For this
paper, we exclude obvious emission of clearly detected companions (e.g. ‘Arrow
galaxy’ Vernet et al., 2017). A systematic census of the companion galaxies in
the MUSE HzRG fields will be reported in a future paper.

###  5.5 Gas distribution on large scales

So far, we have focused our discussion on morphological measurements sensitive
to the high-surface brightness part of the Ly$\alpha$ nebulae and plausible
explanation of environmental effect on $\sim 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_{\rm weight}-\rm{PA_{radio}}|$
towards values $<30^{\circ}$, Fig.13) through interaction with the gaseous
medium. Beyond the extent of the jets, we use $|\theta_{\rm
unweight}-\rm{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_{\rm H}\sim
10^{-2}\,\mathrm{cm^{-3}}$), the jet ‘working surface’ (hot spots) will
leading to brighter radio emission. Such a scenario would introduce a bias in
observing jets preferably along directions of higher densities in flux-limited
samples and we might expect a correlation between the morphology of a large
scale halo and radio jet axis. Our observations presented in Figures 11 and 13
is in agreement with such a scenario.

The detection of molecular gas in HzRGs along the jet axis but beyond hot
spots provides further evidence ($<20^{\circ}$, which suggests the jet runs
into filament of cold gas Emonts et al., 2014). If this direction indeed
traces a filament of the cosmic web (e.g. West, 1991; Umehata et al., 2019;
Bacon et al., 2021) with a higher density of companion galaxies, then this is
also in agreement with the discussion presented in Sect. 5.4. Humphrey et al.
(2007) found evidence of infalling CGM gas on to the HzRGs which bridges the
radio galaxy and to the large-scale (IGM) gas. If the CGM gas is being
accreted onto the central HzRGs through these higher density distributions, it
would reflect the scenario proposed by Humphrey et al. (2007).

###  5.6 Inflow from the cosmic web?

The relatively large $|\theta-\rm{PA_{radio}}|$ ($81.4^{\circ}$ and
$45.3^{\circ}$ for weighted and un-weighted, respectively) detected in
4C+04.11 is inconsistent with other targets (see Figure 13). This suggests
that the nebula of 4C+04.11 is elongated in the East-West direction on both
small ($<20$ pkpc, scope of jet) and large scales ($\gtrsim 20$ pkpc, Sect.
4.3.1 and 5.4). The tile with the largest distance from AGN position (# 74,
Fig. 22) has a filamentary like structure. There is no known evidence that
4C+04.11 resides at the center of a dense cluster-like environment (e.g.
Kikuta et al., 2017) which makes it unlikely that the asymmetry is caused by
contamination of nearby halos as discussed in Sect. 5.4.

High velocity shocks ($\gtrsim 100\,\rm{km\,s^{-1}}$, Allen et al., 2008) can
heat the gas which will then cool by radiating UV photons. Shocks could be
caused by inflows and power the observed Ly$\alpha$ emission. We estimated the
dark matter halo of 4C+04.11 is on the order of $M_{\rm DM}\sim