10^{13}\,\mathrm{M_{\odot}}$ from it stellar mass, $M_{\star}\sim
10^{11}\,\mathrm{M_{\odot}}$ with an empirical $M_{\star}-M_{\rm DM}$ relation
(Wang et al., 2021b). The surface brightness measured in the farthest tile #74
is $\sim 1.2\times 10^{-17}\,\mathrm{erg\,s^{-1}\,cm^{-2}\,arcsec^{-2}}$ which
is similar to (or slightly higher than) the simulated value from Rosdahl &
Blaizot (2012). We converted our measurement to $z=3$ which is redshift
reported in the simulation. Our measurement and the simulation are consistent
given that dark matter halo of 4C+04.11 is likely more massive than the one in
the simulation ($M_{\rm DM}\sim 10^{12}\,\mathrm{M_{\odot}}$).

Given the estimated dark matter halo mass, the virial radius and virial
velocity can then be calculated as $r_{\rm vir}\sim 100$ kpc and $v_{\rm
vir}\simeq 580\,\mathrm{km\,s^{-1}}$. From the centre of the tile #74, we can
derive its distance to the central AGN as $\sim 60$ kpc. We note that this is
the projected distance averaged for spaxels in the tile. The actual distance
of the filament could be farther. Even though it is closer, Nelson et al.
(2016) simulated gas accretion into $10^{12}\,M_{\odot}$ halos at $z=2$ and
found that the accretion flow structure can remain filamentary within $r_{\rm
vir}$ ($\sim 0.5r_{\rm vir}$). The projected $v_{50}$ (velocity offset) of
tile #74 is $-536\,\mathrm{km\,s^{-1}}$ which is consistent with $v_{\rm
vir}$. Thus, our measurements for the projected distance of the tile and its
velocity offset are consistent with a scenario where the Ly$\alpha$ emission
in tile #74 may be tracing shock driven by inflows. If the distance we observe
is indeed $<r_{\rm vir}$, then the post-virial accretion could be more
complicated with multi gas phase components and fragmentation involved (e.g.
Cornuault et al., 2018). The broad line width of 4C+04.11 at large radii (Fig.
7) may indicate the disturbed nature of the presumed accretion flow.

The discussion in this section based on the $|\theta-\rm{PA_{radio}}|$ and
morphology of bin #74 of 4C+04.11. While the angle difference is a clear
outlier in the sample, the tile shape may depend on the implemented method for
nebula extent (Section 3.1.1). However, we checked the Ly$\alpha$ narrow band
image ($6690\,\AA$ to $6707\,\AA$) directly collapsed from the data cube and
confirm the indication of emission from this region. As for the nearby
potential emission (north to bin #74, Fig. 22), we confirm by spectral
extraction that there is no S/N$>3$ detection. Even if the Ly$\alpha$ emission
is detected in this additional region, it would still be consistent with the
accretion scenario. Complex filamentary structures are seen in simulations
(see Rosdahl & Blaizot, 2012, and their Fig. 7) for accretion in higher mass
dark matter halo. Nevertheless, we caution that this analysis only considers
one possibility. As the data is near the detection limit, deeper observation
are needed before more conclusive results can be derived.

##  6 Conclusions

In this paper, we study the intrinsic Ly$\alpha$ nebulae around a sample of
eight high-redshift radio galaxies using optical IFU observations obtained
with MUSE. We link our observations to results from the literature for other
quasars (mainly type-1 quasars) at similar redshifts.

We have developed a new method to measure the maximum extent of the nebulae
with improved sensitivity to low surface brightness emission. We also have
developed a new method to tessellate the Ly$\alpha$ maps (Sect. 3.1). We have
detected the Ly$\alpha$ emission at scales $\gtrsim 100\,\rm{pkpc}$ from the
central AGN, down to an observed surface brightness of $\sim 2-20\times
10^{-19}\,\rm{erg\,s^{-1}\,cm^{-2}\,arcsec^{-2}}$.

We summarise our results as follows: The Ly$\alpha$ emission line profiles of
all sources are affected by multiple and deep absorption troughs out to the
edge of the nebulae. We have corrected for this H i absorption (Sect.3.2) and
constructed maps of the intrinsic Ly$\alpha$ (Sect. 4.1) emission and also
measured the kinematic properties of the Ly$\alpha$ emission spatially
resolved.

We first investigated radial dependencies of the surface brightness, velocity
shift and velocity width of our sample. We found that circularly averaged
profiles of the intrinsic (absorption-corrected) Ly$\alpha$ surface brightness
are brighter and more extended than type-1 quasar samples (Sect. 4.2.1). We
did not find major differences when we investigate the radial profiles along
the approaching and receding direction of the radio jets, respectively (Sect.
4.2.2). The surface brightness radial profiles can generally be described by
an exponential drop-off for the inner high surface brightness part and a power
law for the more extended part (Sect. 4.2.3). We did not find evidence of
ordered gas bulk motion from the $v_{50}$ radial profile (Sect. 4.2.4). For
four targets, the $v_{50}$ profiles at radii within the jet hot-spots range
show a similar gradient as the jet-driven outflow (Fig. 7). The $W_{80}$
profiles show relatively large values ($\gtrsim 1500\,\mathrm{km\,s^{-1}}$,
Fig. 7, Sect. 4.2.4) at all radii and three targets show a decrease beyond the
distance of jet hot spots is indicative of jet-gas interactions.

We quantitatively studied the morphology of the HzRG nebulae (for both
observed and intrinsic emission) and compared our results to other quasar
samples (uncorrected for absorption, Sect. 4.3). We found that our nebulae
are, in general, more asymmetric than nebulae of type-1 quasars and are more
similar to type-2 quasars (Sect. 4.3.2), although, our sampled sources differ
in their measure of asymmetry as a function of radius (Sect. 4.3.3).
Furthermore, we found that the more asymmetric and larger the nebulae are, the
greater the offset between the centroid of the nebulae and AGN positions
($d_{\rm max}-d_{\rm AGN-neb}$ and $e_{\rm weight}-d_{\rm AGN-neb}$ in Sect.
4.3.4, Fig. 11).

Ly$\alpha$ is a complicated emission line that can be heavily affected by
absorption and resonant scattering which, as we demonstrated, needs to be
accounted for. Assuming type-1 and type-2 quasars have similar intrinsic
shapes of their nebulae, the difference of the nebulae asymmetry morphology
between our sample and other quasars can be explained using the AGN
unification model where the orientation of the ionisation cone is the
fundamental parameter (Sect. 5.2). Resonant scattering of the Ly$\alpha$
emission can result in the observed $d_{\rm max}-d_{\rm AGN-neb}$ and $e_{\rm
weight}-d_{\rm AGN-neb}$ relation in our sample and partially explain the
shape of the radial profile and the kinematic profiles (Sect. 5.3). We also
found evidence in our HzRGs that the extended nebulae are affected by
Ly$\alpha$ emissions from nearby companions (Sect. 5.4) and that CGM gas has
higher density distribution along a specific direction which is coincident
with the direction of the radio jet (Sect. 5.5). There is also a possibility
that, in our observations, we are capturing inflows from cosmic web (Sect.
5.6).

In this paper we have shown that measurements of Ly$\alpha$ nebulae around
high-redshift AGN can act as a probe of the environment of AGN and their host
galaxies. We have found fundamental differences between the extended ionised
nebulae of different types of quasars at Cosmic Noon and beyond. Due to its
resonant nature, it is a challenge to use Ly$\alpha$ as a tracer of gas
kinematics especially out to the $\sim 100$ pkpc scale. Nevertheless,
Ly$\alpha$ line observations offer some insight into the state of CGM gas at a
time when it is simultaneously being impacted by more than one physical
mechanism (e.g. quasar outflow, jet-gas interaction and cosmic inflow). These
kind of observations will be particularly useful for future simulation works.
The MUSE data only reveals the tip of the iceberg. In upcoming papers, the H i
absorbers will be reported together with the high-resolution spectroscopy
data. In addition, all the UV emission lines covered by MUSE will be studied
in detail and this will provide crucial information on the ionisation,
metallicity and AGN feedback in the quasar nebulae. Furthermore, scheduled
JWST observations will be available for four HzRGs in our sample and this will
provide unprecedented details of the AGN and host galaxy connection on scales
of $\lesssim 1$ kpc ($\sim 0.05$″).

###### Acknowledgements.

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

## References

  * Adam et al. (1997) Adam, G., Rocca-Volmerange, B., Gerard, S., Ferruit, P., & Bacon, R. 1997, A&A, 326, 501
  * Ahumada et al. (2020) Ahumada, R., Prieto, C. A., Almeida, A., et al. 2020, ApJS, 249, 3
  * Allen et al. (2008) Allen, M. G., Groves, B. A., Dopita, M. A., Sutherland, R. S., & Kewley, L. J. 2008, ApJS, 178, 20
  * Andrae et al. (2010) Andrae, R., Schulze-Hartung, T., & Melchior, P. 2010, arXiv e-prints, arXiv:1012.3754
  * Antonucci (1993) Antonucci, R. 1993, ARA&A, 31, 473
  * Arrigoni Battaia et al. (2016) Arrigoni Battaia, F., Hennawi, J. F., Cantalupo, S., & Prochaska, J. X. 2016, ApJ, 829, 3
  * Arrigoni Battaia et al. (2019) Arrigoni Battaia, F., Hennawi, J. F., Prochaska, J. X., et al. 2019, MNRAS, 482, 3162
  * Astropy Collaboration et al. (2018) Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ, 156, 123
  * Bacon et al. (2010) Bacon, R., Accardo, M., Adjali, L., et al. 2010, in Ground-based and Airborne Instrumentation for Astronomy III, ed. I. S. McLean, S. K. Ramsay, & H. Takami, Vol. 7735, International Society for Optics and Photonics (SPIE), 131 – 139
  * Bacon et al. (2021) Bacon, R., Mary, D., Garel, T., et al. 2021, A&A, 647, A107
  * Bacon et al. (2016) Bacon, R., Piqueras, L., Conseil, S., Richard, J., & Shepherd, M. 2016, MPDAF: MUSE Python Data Analysis Framework, Astrophysics Source Code Library, record ascl:1611.003
  * Bacon et al. (2014) Bacon, R., Vernet, J., Borisova, E., et al. 2014, The Messenger, 157, 13
  * Basu-Zych & Scharf (2004) Basu-Zych, A. & Scharf, C. 2004, ApJ, 615, L85
  * Best et al. (2000) Best, P. N., Röttgering, H. J. A., & Longair, M. S. 2000, MNRAS, 311, 23
  * Bicknell et al. (2000) Bicknell, G. V., Sutherland, R. S., van Breugel, W. J. M., et al. 2000, ApJ, 540, 678
  * Binette et al. (2000) Binette, L., Kurk, J. D., Villar-Martín, M., & Röttgering, H. J. A. 2000, A&A, 356, 23
  * Bolatto et al. (2013) Bolatto, A. D., Wolfire, M., & Leroy, A. K. 2013, ARA&A, 51, 207
  * Borisova et al. (2016) Borisova, E., Cantalupo, S., Lilly, S. J., et al. 2016, ApJ, 831, 39
  * Byrohl et al. (2021) Byrohl, C., Nelson, D., Behrens, C., et al. 2021, MNRAS, 506, 5129
  * Cai et al. (2019) Cai, Z., Cantalupo, S., Prochaska, J. X., et al. 2019, ApJS, 245, 23
  * Cantalupo et al. (2014) Cantalupo, S., Arrigoni-Battaia, F., Prochaska, J. X., Hennawi, J. F., & Madau, P. 2014, Nature, 506, 63
  * Cantalupo et al. (2012) Cantalupo, S., Lilly, S. J., & Haehnelt, M. G. 2012, MNRAS, 425, 1992
  * Cappellari & Copin (2003) Cappellari, M. & Copin, Y. 2003, MNRAS, 342, 345
  * Carilli et al. (1997) Carilli, C. L., Röttgering, H. J. A., van Ojik, R., Miley, G. K., & van Breugel, W. J. M. 1997, ApJS, 109, 1
  * Chambers et al. (1987) Chambers, K. C., Miley, G. K., & van Breugel, W. 1987, Nature, 329, 604
  * Chang et al. (2022) Chang, S.-J., Yang, Y., Seon, K.-I., Zabludoff, A., & Lee, H.-W. 2022, arXiv e-prints, arXiv:2212.09630
  * Christensen et al. (2006) Christensen, L., Jahnke, K., Wisotzki, L., & Sánchez, S. F. 2006, A&A, 459, 717
  * Cigan (2019) Cigan, P. 2019, MultiColorFits: Colorize and combine multiple fits images for visually aesthetic scientific plots, Astrophysics Source Code Library, record ascl:1909.002
  * Cornuault et al. (2018) Cornuault, N., Lehnert, M. D., Boulanger, F., & Guillard, P. 2018, A&A, 610, A75
  * Costa et al. (2022) Costa, T., Arrigoni Battaia, F., Farina, E. P., et al. 2022, MNRAS, 517, 1767
  * Cowie & Hu (1998) Cowie, L. L. & Hu, E. M. 1998, AJ, 115, 1319
  * Dawson et al. (2007) Dawson, S., Rhoads, J. E., Malhotra, S., et al. 2007, ApJ, 671, 1227
  * De Breuck et al. (2010) De Breuck, C., Seymour, N., Stern, D., et al. 2010, ApJ, 725, 36
  * De Breuck et al. (1999) De Breuck, C., van Breugel, W., Minniti, D., et al. 1999, A&A, 352, L51
  * De Breuck et al. (2002) De Breuck, C., van Breugel, W., Stanford, S. A., et al. 2002, AJ, 123, 637
  * den Brok et al. (2020) den Brok, J. S., Cantalupo, S., Mackenzie, R., et al. 2020, MNRAS, 495, 1874
  * Dey et al. (2005) Dey, A., Bian, C., Soifer, B. T., et al. 2005, ApJ, 629, 654
  * Dijkstra (2014) Dijkstra, M. 2014, PASA, 31, e040
  * Dijkstra (2019) Dijkstra, M. 2019, Correction to: Physics of Ly$\alpha$ Radiative Transfer, ed. A. Verhamme, P. North, S. Cantalupo, & H. Atek (Berlin, Heidelberg: Springer Berlin Heidelberg), C1–C2
  * Drouart et al. (2012) Drouart, G., De Breuck, C., Vernet, J., et al. 2012, A&A, 548, A45
  * Drouart et al. (2014) Drouart, G., De Breuck, C., Vernet, J., et al. 2014, A&A, 566, A53
  * Duncan et al. (2022) Duncan, K. J., Windhorst, R. A., Koekemoer, A. M., et al. 2022, arXiv e-prints, arXiv:2212.09769
  * Eales (1992) Eales, S. A. 1992, ApJ, 397, 49
  * Emonts et al. (2014) Emonts, B. H. C., Norris, R. P., Feain, I., et al. 2014, MNRAS, 438, 2898
  * ESO CPL Development Team (2015) ESO CPL Development Team. 2015, EsoRex: ESO Recipe Execution Tool, Astrophysics Source Code Library, record ascl:1504.003
  * Evans et al. (2010) Evans, I. N., Primini, F. A., Glotfelty, K. J., et al. 2010, ApJS, 189, 37
  * Falkendal et al. (2019) Falkendal, T., De Breuck, C., Lehnert, M. D., et al. 2019, A&A, 621, A27
  * Falkendal et al. (2021) Falkendal, T., Lehnert, M. D., Vernet, J., De Breuck, C., & Wang, W. 2021, A&A, 645, A120
  * Farina et al. (2019) Farina, E. P., Arrigoni-Battaia, F., Costa, T., et al. 2019, ApJ, 887, 196
  * Foreman-Mackey et al. (2013) Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306
  * Fossati et al. (2021) Fossati, M., Fumagalli, M., Lofthouse, E. K., et al. 2021, MNRAS, 503, 3044
  * Francis et al. (2001) Francis, P. J., Williger, G. M., Collins, N. R., et al. 2001, ApJ, 554, 1001
  * Gaia Collaboration et al. (2021) Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2021, A&A, 649, A1
  * Gronke et al. (2015) Gronke, M., Bull, P., & Dijkstra, M. 2015, ApJ, 812, 123
  * Gronke & Dijkstra (2016) Gronke, M. & Dijkstra, M. 2016, ApJ, 826, 14
  * Gronke et al. (2016) Gronke, M., Dijkstra, M., McCourt, M., & Oh, S. P. 2016, ApJ, 833, L26
  * Gullberg et al. (2016) Gullberg, B., De Breuck, C., Lehnert, M. D., et al. 2016, A&A, 586, A124
  * Hanuschik (2003) Hanuschik, R. W. 2003, A&A, 407, 1157
  * Harris et al. (2020) Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357
  * Harwood et al. (2013) Harwood, J. J., Hardcastle, M. J., Croston, J. H., & Goodger, J. L. 2013, MNRAS, 435, 3353
  * He et al. (2018) He, Z., Sun, A.-L., Zakamska, N. L., et al. 2018, MNRAS, 478, 3614
  * Heckman et al. (1991a) Heckman, T. M., Lehnert, M. D., Miley, G. K., & van Breugel, W. 1991a, ApJ, 381, 373
  * Heckman et al. (1991b) Heckman, T. M., Lehnert, M. D., van Breugel, W., & Miley, G. K. 1991b, ApJ, 370, 78
  * Hennawi & Prochaska (2013) Hennawi, J. F. & Prochaska, J. X. 2013, ApJ, 766, 58
  * Hippelein & Meisenheimer (1993) Hippelein, H. & Meisenheimer, K. 1993, Nature, 362, 224
  * Hu & McMahon (1996) Hu, E. M. & McMahon, R. G. 1996, Nature, 382, 231
  * Humphrey (2019) Humphrey, A. 2019, MNRAS, 486, 2102
  * Humphrey et al. (2013) Humphrey, A., Vernet, J., Villar-Martín, M., et al. 2013, ApJ, 768, L3
  * Humphrey et al. (2007) Humphrey, A., Villar-Martín, M., Fosbury, R., et al. 2007, MNRAS, 375, 705
  * Humphrey et al. (2006) Humphrey, A., Villar-Martín, M., Fosbury, R., Vernet, J., & di Serego Alighieri, S. 2006, MNRAS, 369, 1103
  * Humphrey et al. (2008) Humphrey, A., Villar-Martín, M., Sánchez, S. F., et al. 2008, MNRAS, 390, 1505
  * Hunter (2007) Hunter, J. D. 2007, Computing in Science & Engineering, 9, 90
  * Jarvis et al. (2003) Jarvis, M. J., Wilman, R. J., Röttgering, H. J. A., & Binette, L. 2003, MNRAS, 338, 263
  * Kikuta et al. (2017) Kikuta, S., Imanishi, M., Matsuoka, Y., et al. 2017, ApJ, 841, 128
  * Kluyver et al. (2016) Kluyver, T., Ragan-Kelley, B., Pérez, F., et al. 2016, in Positioning and Power in Academic Publishing: Players, Agents and Agendas, ed. F. Loizides & B. Schmidt, IOS Press, 87 – 90
  * Kolwa et al. (2023) Kolwa, S., De Breuck, C., Vernet, J., et al. 2023, arXiv e-prints, arXiv:2309.00459
  * Kolwa et al. (2019) Kolwa, S., Vernet, J., De Breuck, C., et al. 2019, A&A, 625, A102
  * Krogager (2018) Krogager, J.-K. 2018, arXiv e-prints, arXiv:1803.01187
  * Lau et al. (2022) Lau, M. W., Hamann, F., Gillette, J., et al. 2022, MNRAS, 515, 1624
  * Leclercq et al. (2017) Leclercq, F., Bacon, R., Wisotzki, L., et al. 2017, A&A, 608, A8
  * Li et al. (2022) Li, Z., Steidel, C. C., Gronke, M., Chen, Y., & Matsuda, Y. 2022, MNRAS, 513, 3414
  * Liu et al. (2013) Liu, G., Zakamska, N. L., Greene, J. E., Nesvadba, N. P. H., & Liu, X. 2013, MNRAS, 436, 2576
  * Mackenzie et al. (2021) Mackenzie, R., Pezzulli, G., Cantalupo, S., et al. 2021, MNRAS, 502, 494
  * Marino et al. (2019) Marino, R. A., Cantalupo, S., Pezzulli, G., et al. 2019, ApJ, 880, 47
  * Martin et al. (2014) Martin, D. C., Chang, D., Matuszewski, M., et al. 2014, ApJ, 786, 107
  * Matsuda et al. (2004) Matsuda, Y., Yamada, T., Hayashino, T., et al. 2004, AJ, 128, 569
  * McCarthy (1993) McCarthy, P. J. 1993, ARA&A, 31, 639
  * McCarthy et al. (1991) McCarthy, P. J., van Breugel, W., & Kapahi, V. K. 1991, ApJ, 371, 478
  * Miley & De Breuck (2008) Miley, G. & De Breuck, C. 2008, A&A Rev., 15, 67
  * Miley et al. (2004) Miley, G. K., Overzier, R. A., Tsvetanov, Z. I., et al. 2004, Nature, 427, 47
  * Morais et al. (2017) Morais, S. G., Humphrey, A., Villar-Martín, M., et al. 2017, MNRAS, 465, 2698
  * Nelson et al. (2016) Nelson, D., Genel, S., Pillepich, A., et al. 2016, MNRAS, 460, 2881
  * Nesvadba et al. (2017a) Nesvadba, N. P. H., De Breuck, C., Lehnert, M. D., Best, P. N., & Collet, C. 2017a, A&A, 599, A123
  * Nesvadba et al. (2017b) Nesvadba, N. P. H., Drouart, G., De Breuck, C., et al. 2017b, A&A, 600, A121
  * Nesvadba et al. (2007) Nesvadba, N. P. H., Lehnert, M. D., De Breuck, C., Gilbert, A., & van Breugel, W. 2007, A&A, 475, 145
  * Nesvadba et al. (2008) Nesvadba, N. P. H., Lehnert, M. D., De Breuck, C., Gilbert, A. M., & van Breugel, W. 2008, A&A, 491, 407
  * Newville et al. (2016) Newville, M., Stensitzki, T., Allen, D. B., et al. 2016, Lmfit: Non-Linear Least-Square Minimization and Curve-Fitting for Python
  * Noirot et al. (2018) Noirot, G., Stern, D., Mei, S., et al. 2018, ApJ, 859, 38
  * Noirot et al. (2016) Noirot, G., Vernet, J., De Breuck, C., et al. 2016, ApJ, 830, 90
  * Ono et al. (2021) Ono, Y., Itoh, R., Shibuya, T., et al. 2021, ApJ, 911, 78
  * Ouchi et al. (2020) Ouchi, M., Ono, Y., & Shibuya, T. 2020, ARA&A, 58, 617
  * Parijskij et al. (2014) Parijskij, Y. N., Thomasson, P., Kopylov, A. I., et al. 2014, MNRAS, 439, 2314
  * Pentericci et al. (1999) Pentericci, L., Röttgering, H. J. A., Miley, G. K., et al. 1999, A&A, 341, 329
  * Pentericci et al. (2000) Pentericci, L., Van Reeven, W., Carilli, C. L., Röttgering, H. J. A., & Miley, G. K. 2000, A&AS, 145, 121
  * Prescott et al. (2015) Prescott, M. K. M., Martin, C. L., & Dey, A. 2015, ApJ, 799, 62
  * Rauch (1998) Rauch, M. 1998, ARA&A, 36, 267
  * Reuland et al. (2003) Reuland, M., van Breugel, W., Röttgering, H., et al. 2003, ApJ, 592, 755
  * Reyes et al. (2008) Reyes, R., Zakamska, N. L., Strauss, M. A., et al. 2008, AJ, 136, 2373
  * Rosdahl & Blaizot (2012) Rosdahl, J. & Blaizot, J. 2012, MNRAS, 423, 344
  * Rottgering et al. (1995) Rottgering, H. J. A., Hunstead, R. W., Miley, G. K., van Ojik, R., & Wieringa, M. H. 1995, MNRAS, 277, 389
  * Saito et al. (2006) Saito, T., Shimasaku, K., Okamura, S., et al. 2006, ApJ, 648, 54
  * Sanderson et al. (2021) Sanderson, K. N., Prescott, M. K. M., Christensen, L., Fynbo, J., & Møller, P. 2021, ApJ, 923, 252
  * Seymour et al. (2007) Seymour, N., Stern, D., De Breuck, C., et al. 2007, ApJS, 171, 353
  * Shimasaku et al. (2006) Shimasaku, K., Kashikawa, N., Doi, M., et al. 2006, PASJ, 58, 313
  * Silva et al. (2018a) Silva, M., Humphrey, A., Lagos, P., et al. 2018a, MNRAS, 481, 1401
  * Silva et al. (2018b) Silva, M., Humphrey, A., Lagos, P., et al. 2018b, MNRAS, 474, 3649
  * Soto et al. (2016) Soto, K. T., Lilly, S. J., Bacon, R., Richard, J., & Conseil, S. 2016, MNRAS, 458, 3210
  * Steidel et al. (2000) Steidel, C. C., Adelberger, K. L., Shapley, A. E., et al. 2000, ApJ, 532, 170
  * Swinbank et al. (2015) Swinbank, A. M., Vernet, J. D. R., Smail, I., et al. 2015, MNRAS, 449, 1298
  * Tepper-García (2006) Tepper-García, T. 2006, MNRAS, 369, 2025
  * Tepper-García (2007) Tepper-García, T. 2007, MNRAS, 382, 1375
  * Tumlinson et al. (2017) Tumlinson, J., Peeples, M. S., & Werk, J. K. 2017, ARA&A, 55, 389
  * Umehata et al. (2019) Umehata, H., Fumagalli, M., Smail, I., et al. 2019, Science, 366, 97
  * van Breugel et al. (1998) van Breugel, W. J. M., Stanford, S. A., Spinrad, H., Stern, D., & Graham, J. R. 1998, ApJ, 502, 614
  * van Ojik et al. (1996) van Ojik, R., Roettgering, H. J. A., Carilli, C. L., et al. 1996, A&A, 313, 25
  * van Ojik et al. (1997) van Ojik, R., Roettgering, H. J. A., Miley, G. K., & Hunstead, R. W. 1997, A&A, 317, 358
  * Vayner et al. (2021) Vayner, A., Zakamska, N. L., Riffel, R. A., et al. 2021, MNRAS, 504, 4445
  * Vayner et al. (2023) Vayner, A., Zakamska, N. L., Sabhlok, S., et al. 2023, MNRAS, 519, 961
  * Veilleux et al. (2005) Veilleux, S., Cecil, G., & Bland-Hawthorn, J. 2005, ARA&A, 43, 769
  * Venemans et al. (2007) Venemans, B. P., Röttgering, H. J. A., Miley, G. K., et al. 2007, A&A, 461, 823
  * Verhamme et al. (2006) Verhamme, A., Schaerer, D., & Maselli, A. 2006, A&A, 460, 397
  * Vernet et al. (2017) Vernet, J., Lehnert, M. D., De Breuck, C., et al. 2017, A&A, 602, L6
  * Villar-Martin et al. (1996) Villar-Martin, M., Binette, L., & Fosbury, R. A. E. 1996, A&A, 312, 751
  * Villar-Martín et al. (2007a) Villar-Martín, M., Humphrey, A., De Breuck, C., et al. 2007a, MNRAS, 375, 1299
  * Villar-Martín et al. (2006) Villar-Martín, M., Sánchez, S. F., De Breuck, C., et al. 2006, MNRAS, 366, L1
  * Villar-Martín et al. (2007b) Villar-Martín, M., Sánchez, S. F., Humphrey, A., et al. 2007b, MNRAS, 378, 416
  * Villar-Martín et al. (2003) Villar-Martín, M., Vernet, J., di Serego Alighieri, S., et al. 2003, MNRAS, 346, 273
  * Villar-Martín et al. (2002) Villar-Martín, M., Vernet, J., di Serego Alighieri, S., et al. 2002, MNRAS, 336, 436
  * Virtanen et al. (2020) Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nature Methods, 17, 261
  * Wang et al. (2021a) Wang, W., De Breuck, C. A., Lehnert, M. D., et al. 2021a, Zooming into the monster’s mouth: tracing feedback from their hosts to circumgalactic medium in z=3.5 radio-loud AGN, JWST Proposal. Cycle 1, ID. #1970
  * Wang et al. (2021b) Wang, W., Wylezalek, D., De Breuck, C., et al. 2021b, A&A, 654, A88
  * Weidinger et al. (2004) Weidinger, M., Møller, P., & Fynbo, J. P. U. 2004, Nature, 430, 999
  * Weidinger et al. (2005) Weidinger, M., Møller, P., Fynbo, J. P. U., & Thomsen, B. 2005, A&A, 436, 825
  * Weilbacher et al. (2020) Weilbacher, P. M., Palsa, R., Streicher, O., et al. 2020, A&A, 641, A28
  * West (1991) West, M. J. 1991, ApJ, 379, 19
  * Wilman et al. (2004) Wilman, R. J., Jarvis, M. J., Röttgering, H. J. A., & Binette, L. 2004, MNRAS, 351, 1109
  * Wisotzki et al. (2018) Wisotzki, L., Bacon, R., Brinchmann, J., et al. 2018, Nature, 562, 229
  * Wylezalek et al. (2013) Wylezalek, D., Galametz, A., Stern, D., et al. 2013, ApJ, 769, 79
  * Wylezalek et al. (2014) Wylezalek, D., Vernet, J., De Breuck, C., et al. 2014, ApJ, 786, 17
  * Yang et al. (2010) Yang, Y., Zabludoff, A., Eisenstein, D., & Davé, R. 2010, ApJ, 719, 1654
  * Yang et al. (2009) Yang, Y., Zabludoff, A., Tremonti, C., Eisenstein, D., & Davé, R. 2009, ApJ, 693, 1579
  * Zakamska et al. (2016) Zakamska, N. L., Hamann, F., Pâris, I., et al. 2016, MNRAS, 459, 3144
  * Zhang et al. (2023) Zhang, S., Cai, Z., Xu, D., et al. 2023, ApJ, 952, 124
  * Zirm et al. (2005) Zirm, A. W., Overzier, R. A., Miley, G. K., et al. 2005, ApJ, 630, 68

##  Appendix A Detection map and tessellation procedures

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

Table 4: Optimal detection map parameters combination.

HzRG | $T_{S/N}$ | $n_{\lambda}$ |  $2.355\sigma_{\rm max}$\rm{{\dagger}}\rm{{\dagger}}$\rm{{\dagger}}$Maximum Gaussian kernel scale in FWHM. The unit is in pixel unit. | $n_{\rm c}$
---|---|---|---|---
MRC 0943-242 | 3.5 | 4 | 25 | 2
MRC 0316-257 | 3.5 | 5 | 30 | 3
TN J0205+2242 | 3.5 | 5 | 25 | 2
TN J0121+1320 | 2.5 | 3 | 35 | 4
4C+03.24 | 3.5 | 5 | 20 | 2
4C+19.71 | 3.5 | 2 | 30 | 3
TN J1338-1942 | 3.5 | 2 | 30 | 2
4C+04.11 | 3.5 | 4 | 25 | 2
777

Figure 14: Schematic presentation of the detection map construction. (a)
Average $n_{\lambda}$ number of images around each of the image slices (cyan)
from the continuum-subtracted cube (dark green cube). (b) Spatially smooth
each of the averaged image with Gaussian kernel. The algorithm will increase
the kernel size ($\sigma_{0}<\sigma_{1}$) to smooth the spaxels that not
passing the $T_{\rm S/N}$, S/N threshold, at each steps until the maximum,
$\sigma_{\rm max}$, has been reached. (c) Combine the adaptively smoothed
images into the smoothed cube (lighter green cube). (d) Check the smoothed
spectrum (long black rectangular box) at location (x, y). The position is
preserved when there are $n_{\rm c}$ consecutive number of pixels selected
with signal detection (red ’S’) in previous steps. (f) Construct the detection
map.

The basic idea for detection map is to use the 3D information of the MUSE data
cube to determine the maximum spatial extent of the Ly$\alpha$ emission.
Smoothing technique is the key process in this procedure to guarantee the
faint structures are captured. To assist with conveying the procedure, the
process is shown schematically in Fig. 14. In this procedure, we can select
the faint end of the emission nebula to the surface brightness detection
limit.

Before the generation of detection map, we do a step of continuum subtraction
for each spectra spatially centred around the AGN and spectrally around their
Ly$\alpha$ wavelength range ($-5000-5000\,\rm{km\,s^{-1}}$ or
$-6000-6000\,\rm{km\,s^{-1}}$ for 4C+03.24 which has broader line width) with
the emission line masked. If there are sky-lines, emission or absorption
features from foreground targets in the unmasked range, we do further masking.
Continuum from the host galaxy, central AGN and foreground objects need to be
subtracted to minimise their contamination. In this way, we exclusively focus
on the line-emission nebula. The choice of flat or linear continuum do not
have significant impact on the Ly$\alpha$ fitting (Wang et al. 2021b). Hence
to better account for the slope from foreground stars, we subtract a first-
order polynomial from the cube.

The 3D adaptive smoothing follows the Martin et al. (2014) procedure with
adaptation (e.g. Vernet et al. 2017). As described in Sect. 3.1.1, it smooths
each of the image planes of the continuum-subtracted cube adaptively with a
Gaussian kernel over a wavelength range of $\sim 15\,\AA$ around the
Ly$\alpha$ emission. For the image slice at each wavelength, the algorithm
first takes average of adjacent $n_{\lambda}$ number of slices around this
image slice along the wavelength dimension (Fig. 14a). Then it adaptively
smooths the averaged image with a Gaussian kernel starting from the smallest
kernel size, $\sigma_{0}=3/2.355$ pix, until the maximum, $\sigma_{\rm max}$,
is reached (Fig. 14b). Specifically at each step, the algorithm smooths the
spaxels that are not passing the pre-set S/N threshold, $T_{\rm S/N}$. In the
end, the algorithm set the voxels, after being maximally smoothed, that not
pass the $T_{\rm S/N}$ to 0 as no-signal containing (i.e. the voxels that
potentially contain Ly$\alpha$ signals have positive value). In this way, we
preliminarily detect where we have Ly$\alpha$ in the cube (Fig. 14c).

To generate a detection map working on the spatial dimensions and guarantee
line fitting, we perform a further check along the wavelength dimension for
the smoothed cube. Specifically, if the smoothed spectrum at one spatial
location (x, y) has $n_{\rm c}$ numbers of consecutive wavelength pixels
preserved (i.e. have positive values), then we flag this spatial location (x,
y) as signal-contained (see Fig. 14d as an example). The others that do not
pass this check are discarded. In this way, we can construct the detection map
(Fig. 14e).

The question left is to find a combination of the four parameters ($T_{\rm
S/N}$, $n_{\lambda}$, $\sigma_{\rm max}$ and $n_{\rm c}$) which returns an
optimised detection map. For this, we generate a series of detection maps for
each targets with different combination of the four parameters. Then we choose
the one that optimises the spatial selection (i.e. captures large scale extent
while avoid island-like structures). The best combination of the four
parameters for each targets are presented in Table 4. We check the detection
maps constructed using different sets of parameters around our optimal
combination (Table 4) and find that the bulk ($\sim 80\%$) of the selected
spaxels remains the same. If using a stricter set of parameters (e.g. higher
$T_{\rm S/N}$ and/or $n_{\lambda}$), we will miss the low S/N part of the
nebula by comparing with previous individual target studies (Swinbank et al.
2015; Vernet et al. 2017). If using a more relaxed set, the peripheral
regions, mostly disconnected to the bulk of the nebula, with pure noise (after
checking spectrum) will likely be selected. Thus, we are assured that the
constructed detection map represented the observation to the detection limit
and the method is robust against objectiveness. After constructing the map
with the best parameter combination, we perform a final manual selection to
exclude island-like regions ($\gtrsim 1$ arcsec2 in size) which are detached
from the the main nebulae ($>1"$) and could be due to noise or companion. A
further check of the spectra extracted from these regions is also conducted.
Around 50% of the island regions only contain noise. The ones with potential
Ly$\alpha$ signal detected are presented in Appendix A.3 (Fig. 16 and 17 ).
Since these are detached from the extended nebula around the quasar, we do not
include them in this analysis but point out they may trace companion
emissions. We note that we refer to the smoothed cube obtained in this way
(using $T_{S/N}$, $n_{\lambda}$ and $\sigma_{\rm max}$ in Table 4) as
’optimally smoothed cube’. The smoothed cube are only used in constructing the
detection map and not being further used in other analysis. In Appendix A.3
(Fig. 15-22a), we present the smoothed nebula image extracted from the
optimally smoothed cube.

###  A.2 Tessellation procedure

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

The Swinbank et al. (2015) tessellation is not optimal for some of our target
with significantly irregular nebula shape (e.g. MRC0316-257) since it works
with rectangular binning. Hence, we decide to adopt the Voronoi binning
(Cappellari & Copin 2003), which is a sophisticate adaptive spatial binning
algorithm implemented by various studies of IFU data, for this work. Directly
performing Voronoi binning on the image will return detached regions (i.e.
regions of spaxels belong to the same tile but are physically separated), as
the images are S/N limited for the outer nebulae. The solution could be (i)
turning off Weighted Voronoi Tessellation or Centroidal Voronoi Tessellation
(see Cappellari & Copin 2003, for details); or (ii) performing a twofold
Voronoi binning by using a S/N cut for doing separate binning for high and low
S/N parts (M. Cappellari private communication). Solution (i) will return
wedge shape tiles which is a known result (Cappellari & Copin 2003). It is,
however, not desirable for our case since it smears out the radial structures
which is one of the key interest of this work. Therefore, we adapt solution
(ii); the twofold Voronoi binning for our sample.

Firstly, we apply the detection map (Sect. 3.1.1 and Appendix A.1) to the
continuum-subtracted un-smoothed cube to preserve only the Ly$\alpha$ signal
detected spaxels for tessellation. To avoid complication and keep consistency
for the whole sample, we perform the tessellation on a single narrow band
image for each target. Since the Ly$\alpha$ spectra may have different
observed peaks at different spatial locations, we choose the wavelength range,
over which the narrow band image for tessellation will be produced (averaged
by the number of wavelength pixels), to enclose as much of line emissions as
possible and avoid adding pure noise. This range is $\sim 10\AA$ wide. We note
that for 4C+19.71, we select two wavelength ranges at both the blue and red
sides of 5577 Å(sky-line), which is at roughly the systemic redshift of
Ly$\alpha$. In this way, we construct the S/N map from the narrow band image.

Secondly, a Gaussian smooth is performed to the S/N map to minimise the impact
of randomly located spaxels dominated by noise (i.e. further avoid detached
region problem). We use a Gaussian kernel with $\sigma_{s}=3$ in pixel unit.
Then, we apply a S/N cut of 3 to the S/N map to select the high S/N regions
for the first Voronoi binning.

Thirdly, for the selected S/N $>3$ part, we reassign spaxels with S/N¿6 to
S/N=6 and perform the Voronoi binning with a target (S/N)${}_{\rm
vorbin,1st}=$ 15 (12 for 4C+03.24, see later). The reason for the S/N
reassignment is to control the minimum size of the tiles to avoid single
spaxel tile and account for the seeing effect. Because the Ly$\alpha$ nebulae
studied here have high S/N gradient from center to outskirt, performing
Voronoi binning with a high S/N threshold (avoid single spaxel bin in the high
S/N region) will result in tiles with large size for the low S/N part which is
an overkill for the fitting and smears out resolution information. After
reassigning S/N¿6 spaxel and using (S/N)${}_{\rm vorbin,1st}=$15, the minimum
number of spaxels for one tile is 6 ($=3\times 2$ pix2 or $0.6\times 0.4$
arcsec2) which is roughly half of the seeing disk. This is a compromise for
being consistency for the whole sample and also physically connects the
neighbour tiles which is useful for the implemented fitting procedure (Sec.
3). As for 4C+03.24 which was observed in the AO mode, we lower its
(S/N)${}_{\rm vorbin,1st}=$12 to have a minimum number of spaxels for one tile
is 4 ($=2\times 2$ pix2 or $0.4\times 0.4$ arcsec2).

Finally, for the S/N$\leq$3 spaxels left in the first step, we apply a Voronoi
binning with (S/N)${}_{\rm vorbin,2nd}=$8 (11 for 4C+19.71 due to the impact
of sky-line).

In this way, we obtain the tessellation maps where each tile has S/N$>5$. We
note that the S/N stated here is calculated based on the narrow band image
(see above) collapsed $\sim 14$ wavelength elements (may vary for different
targets) and averaged by the number of wavelength elements. This narrow band
image is chosen to enclose as much of the Ly$\alpha$ emission as possible
which is a compensation for the different emission peak due to kinematics.
Hence, a S/N$>5$ in each tile selected in this way is feasible for further
fitting as there will be other signal contained wavelength pixels outside the
range given the broadness of the emission line. For targets where the extent
of the nebula is overlapped with bright foreground stars (TN J1338-1942 and
4C+03.24), we manually assign a $d=1\,\rm arcsec$ circular mask at the
position of the star to minimise its impact. For MRC0316-257, we mask the
region where a known Ly$\alpha$ emitting galaxies at the similar redshift
overlaps in spatial location with the nebula (see Vernet et al. 2017, Arrow
galaxy). These manually masked regions are marked by purple hatches in Fig.
1bcd. The tessellation maps are shown in Appendix A.3. Through the
tessellation, we reach $\sim 2\sigma$ surface brightness limit (as reported in
Table 2) in the faintest tile.

For the convenience of the spatial fitting (Sect. 3.2.2 and Appendix B), we
run an automatic re-numbering algorithm to make physically attached tiles to
be as consecutive in number as possible. In each tile, the spatial spectra
from every spaxel are then extracted and summed according to the tessellation
from the original cube (i.e. the one produced by the reduction, Sect. 2.3,
without further continuum subtraction).

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

In this Appendix, we present the smoothed nebula images (produced based on the
optimally smoothed cube, Fig. 15-22a) and the tessellation maps (Fig. 15-22b)
constructed following Appendix A.1 (Sect. 3.1.1) and A.2 (Sect. 3.1.2),
respectively. The false-colour smooth images are generated using
multicolorfits (Cigan 2019). Each colour represents a pseudo-narrow band
Ly$\alpha$ image constructed from the optimally smoothed cube (Appendix. A.1)
in arbitrary wavelength range with the goal to show different kinematic
structures. Blue, yellow and red are relatively from blue wing, middle and red
wing of the Ly$\alpha$ emission, respectively. We note that the smoothed
nebula images are only used as representation and demonstration of how the
algorithm in AppendixA.1 works. They are not included in the analysis of this
work. For MRC0316-257 and TN J0205+2242, we detected line emissions around
their Ly$\alpha$ wavelength at isolated regions from the main nebula. We did
not include these regions in to account due to its possible origin from
companion galaxies, but show the spectra and region where they are extracted
in Fig. 16 and Fig. 17.

Figure 15: Smoothed nebula image (a) and tessellation map (b) of MRC0943-242.
(a) The purple, yellow and blue colours represent median pseudo-narrow band
Ly$\alpha$ images collapsed arbitrarily from red part, middle and blue part,
respectively, of the smoothed cubes. Figure 16: Smoothed nebula image (a) and
tessellation map (b) of MRC0316-257. The inset spectrum is extracted from the
detached regions (hatched) which is selected having Ly$\alpha$ emissions but
left out from the analysis following Sect. 3.1.1. (a) The color-composed image
is created in a manner similar to that of Figure 15a. Figure 17: Smoothed
nebula image (a) and tessellation map (b) of TN J0205+2242. The inset spectrum
is extracted from the detached region (hatched) which is selected having
Ly$\alpha$ emissions but left out from the analysis following Sect. 3.1.1. (a)
The color-composed image is created in a manner similar to that of Figure 15a.
Figure 18: Smoothed nebula image (a) and tessellation map (b) of TN
J0121+1320. (a) The color-composed image is created in a manner similar to
that of Figure 15a.  Figure 19: Smoothed nebula image (a) and tessellation map
(b) of 4C+03.24. (a) The color-composed image is created in a manner similar
to that of Figure 15a. Figure 20: Smoothed nebula image (a) and tessellation
map (b) of 4C+19.71. (a) The color-composed image is created in a manner
similar to that of Figure 15a. Figure 21: Smoothed nebula image (a) and
tessellation map (b) of TN J1338-1942. (a) The color-composed image is created
in a manner similar to that of Figure 15a. Figure 22: Smoothed nebula image
(a) and tessellation map (b) of 4C+04.11. (a) The color-composed image is
created in a manner similar to that of Figure 15a.

##  Appendix B Fitting procedure

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

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

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

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

To get a the first spatial fit, we fourthly fit the spatial spectra from each
tiles in the order determined in Appendix A.2 with the Least-square method. In
the fitting, we pass the results from the previous tile fit to the next as
initials to minimize the potential randomness introduced by fitting each
spectra freely without any constrain. We note that (i) the constrains from the
previous fit may be relaxed more if the ratio of the integrated observed
fluxes between the tiles are significant ($\sim 50\%$) and/or the distance
ratio is large (¿1.1)888Distance ratio, $\frac{r_{1}+r_{2}}{d_{1,2}}$, is the
ratio between the sum of the radii, $r_{1}+r_{2}$, of the two tiles (estimated
by $r=\sqrt{A/\pi}$, where A is the area of the tile) and the distance between
the geometric center of the two tiles, $d_{1,2}$. and/or the area of the tile
is large ($>9\pi$ arcsec2); (ii) the Doppler parameter, $b$, (Tepper-García
2006) is set to a broad range (Kolwa et al. 2019; Wang et al. 2021b,
$40<b<400\,\rm{km\,s^{-1}}$) because of the $b-N$ degeneracy (e.g. Silva et
al. 2018b). We only constrain the range of column density,
$N_{\rm{{H\textsc{i}}}}$. The continuum model is included in the fitting with
a 0th-order-polynomial for tiles overlapped with continuum detected region
determined in first step. A fixed step function with the step at the
wavelength of the systemic Ly$\alpha$ is used for the highest redshift target,
4C+04.11, due to the heavily absorbed continuum on the blue-side which may be
due to Ly$\alpha$ forest, e.g. Rauch (1998).

Due to the scattering nature of Ly$\alpha$ and presence of broad (outflow)
components, it is important to select which region contains more than one
emission component to better study the broader wing. We fit the set of spatial
spectra two times: (i) first using 1-Gaussian only for the emission and (ii)
then using 2-Gaussian for all spectra. In this way, we can do a simple
$\chi^{2}$ ratio selection between the 2-Gaussian vs. the 1-Gaussian fitting
results and select the tiles favor 2-Gaussian fit with a threshold of $\sim
0.80-0.98$ (depends on different targets). We point that the $\chi^{2}$ value
is not a robust measurement of the fitting quality (e.g. Andrae et al. 2010).
For a quick test in our case, however, it is good enough to the first order.
We note that for 4C+03.24, we stay with the 1-Gaussian fit for all tiles due
to the complicated spatial variance of spectral shapes. We note that this
selection is not crucial in this work since we do not distinguish and
separately interpret different velocity components in the analysis. The
purpose for this step is to consistently ensure that the non-single-Gaussian
line shape is considered without missing flux.

Fifthly, we fit the spatial spectra with 2-Gaussian models for the tiles
determined in the last step and 1-Gaussian for the others. The results from
the previous steps are passed to the next as initial guesses. To keep
consistence, we choose to use the same number of absorbers for the whole map
due to the difficulties in determining where one absorber disappears. This is
a good assumption to first order given that we observe most absorbers on large
extent. For example in Fig. 4, we see the tiles at the nebula edge are also
affected by absorbers. The column density of the absorbers is a free parameter
during the fit. Given the degeneracy between the column density and Doppler
parameter, $b$ (Silva et al. 2018a), we leave the constraint of $b$ to a broad
range following Kolwa et al. (2019) ($40$ to $400\,\mathrm{km\,s^{-1}}$). The
centroid (redshift) of the absorbers are also allowed to vary ($\Delta
z\approx\pm 0.001$) except for 4C+04.11999For this case, we follow the
procedure in Wang et al. (2021b) and fix the positions of five absorbers on
the blue low S/N wing to the value determined from the Master spectrum.. For
spectra extracted from inner tiles ($\lesssim 10\arcsec$, with high S/N), we
constrain their column densities to vary within a $\sim 2$ dex range from the
initial input. We ease the column density constrain for the absorbers from
outer tiles with distance to AGN $\gtrsim 10\arcsec$ such that they can be
given low column density ($\sim 10^{13}\,\mathrm{cm^{-2}}$). In this way, the
absorbers at low S/N regions can have negligible impact ($<0.1$ %) on the
reconstructed flux. We use the same number of absorbers reported in previous
works for some of our targets (MRC 0943-242, TN J0121+1320, TN J1338-1942
Wilman et al. 2004; Swinbank et al. 2015). Otherwise, we use the number
determined in the first step of 1D spectra fitting. For TN J0121+1320, we use
3 absorbers instead of 2 as in Wilman et al. (2004). For 4C+04.11, we use 9
absorbers (instead of 7 as in Wang et al. (2021b)) with 7 of them fixed to the
redshifts determined in 1D spectra fit.

Finally, we run MCMC sampling (Foreman-Mackey et al. 2013) using the results
from last step as initials and with larger boundaries to probe the probability
distribution of the fitted parameters and uncertainties.

We point out one caveat that the low S/N (especially for the tiles at the
edge) will affect the reconstructed intrinsic Ly$\alpha$ flux. We implemented
a test where we artificially decrease the S/N of the master spectrum and do
the absorption correction fit. It shows that the reconstructed intrinsic flux
can vary by a factor of $\sim 2$ when the S/N decreases by one dex (e.g. $\sim
100$ to $\sim 10$). This result is based on the fact that we have a relatively
good initial guess for the fit. We note that in the aforementioned fitting
procedure the fitting parameters are constrained by the neighbouring tiles to
have them physically linked. There may still be the chance that the low S/N
will affect the fitted flux leading to over-correction. This may be the case
for 4C+04.11 (Fig. 12). We present the flux ratio maps between intrinsic and
observed in Appendix C which indicate this possibility ($\sim 50$ for the
tiles at the edge).

##  Appendix C Supplementary maps

Figure 23: Similar figure as Fig. 2$-$4 for MRC0943-242. Figure 24: Similar
figure as 2$-$4 for MRC0316-257. Figure 25: Similar figure as 2$-$4 for
4C+19.71. Figure 26: Similar figure as 2$-$4 for TNJ1338-1942. Figure 27:
Similar figure as 2$-$4 for 4C+04.11. Figure 28: (a) Observed surface
brightness maps in the unit of
$10^{-16}\,\mathrm{erg\,s^{-1}\,cm^{-2}\,arcsec^{-2}}$. (b) Flux ratio maps of
intrinsic and observed surface brightness. The green cross and contours are
the same as Fig. 1 for individual targets, respectively.

We show the intrinsic Ly$\alpha$ maps with selected spectra for MRC0943-242,
MRC0316-257, 4C+19.71, TNJ1338-1942 and 4C+04.11 in Fig. 23, 24, 25, 26 and
27, respectively. The spatial study of the Ly$\alpha$ nebulae for these
sources with MUSE were published previously (e.g. Gullberg et al. 2016; Vernet
et al. 2017; Falkendal et al. 2021; Swinbank et al. 2015; Wang et al. 2021b).
In this paper, we performed a consistent analysis for the full sample with the
new method of correcting for absorption. With the individual spectra mostly
showing the tiles at larger distance away from the AGN, we can find that the
absorption features are observed nearly across the entire nebula.

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

##  Appendix D Supplements of nebula radial profile and morphology analysis

Table 5: Position angle and hot spot distance of the radio jet.

HzRG | PAradio\rm{{\dagger}}\rm{{\dagger}}$\rm{{\dagger}}$Jet position angle. This is measured east from north. We note that we do not quantify the uncertainties which could be $\pm 1$ deg. |  $r_{\rm app.}$\rm{*}\rm{*}$\rm{*}$Distance between the hot spot and the AGN position. $r_{\rm app.}$ and $r_{\rm rec.}$ are distance for approaching and receding hot spot, respectively. |  $r_{\rm rec.}$\rm{*}\rm{*}$\rm{*}$Distance between the hot spot and the AGN position. $r_{\rm app.}$ and $r_{\rm rec.}$ are distance for approaching and receding hot spot, respectively.
---|---|---|---
| deg | arcsec | arcsec
MRC 0943-242 | 74 | 1.79 | 1.76
MRC 0316-257 | 53 | 5.56 | 1.15
TN J0205+2242 | 150 | 1.47 | 1.24
TN J0121+1320 | 90 | 0.82aafootnotemark: $a$ | 0.82aafootnotemark: $a$
4C+03.24 | 146 | 2.00bbfootnotemark: $b$ | 1.75bbfootnotemark: $b$
4C+19.71 | 176 | 3.82 | 4.58
TN J1338-1942 | 150 | 3.78 | 1.18
4C+04.11 | 158 | 1.15ccfootnotemark: $c$ | 1.70ccfootnotemark: $c$
101010 $\rm{a}$$\rm{a}$footnotetext: The radio emission of TN J0121+1320 is
spatially unresolved. The ‘hot spot’ distance is represented by the distance
from the host galaxy position to the either side of $3\sigma$ contour of the
radio image along the east-west direction.$\rm{b}$$\rm{b}$footnotetext:
Approaching hot spot is A1 and receding hot spot is B1 (named after van Ojik
et al. 1996). $\rm{c}$$\rm{c}$footnotetext: Approaching hot spot is knot8 and
receding hot spot is knot1 (named after Parijskij et al. 2014).

Figure 29: Circular apertures for radial profile (Sect. 4.2.1) extraction on
top of intrinsic surface brightness maps (Fig. 1a) of the HzRGs sample. In
each panel, the blue concentric annuli show the apertures where the radial
profile is extracted. The annuli are centered around the host galaxy (central
AGN) position (Table 2). The radii of the annuli in each panel are the same in
the unit of arcsec for consistence except the smallest radius which is set to
be $0.75\times$ seeing for each target. The black bar at the bottom left
corner of each panel indicate the 10 arcsec scale. Figure 30: Ly$\alpha$
position-velocity diagrams (i.e. 2D spectra) of our sample targets extracted
along the radio jet axis. For each target, the left panel shows the diagram
constructed from tessellated observed cube (not continuum-subtracted for the
host galaxy) and the right panel shows the diagram from absorption-corrected
intrinsic cube. The zero offset is set to the position of the central AGN
(Table 2). The direction of approaching and receding sides of the jet (see
Sect. 2.1.3) are marked in the left panel by the blue and red arrows,
respectively. The white horizontal dashed lines represent the furthest extent
of the jet. The dotted green contours are given in arbitrary steps which are
used to guide the eye for the high brightness part. The vertical black shaded
regions in the observed position-velocity diagrams of 4C+03.24 and 4C+19.71
indicate the wavelength ranges affected by the $5577\,\AA$ sky-line.

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

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

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

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

We show the exponential and power law fitted parameters (Sect. 4.2.3) of the
circularly averaged and directional surface brightness profiles in Table 7. In
Fig. 30, we show the position-velocity diagram of the observed and intrinsic
cubes along the jet direction. This can be used as a direct comparison with
Villar-Martín et al. (2003). On each target, we also mark the largest extent
of the radio lobes in both directions in horizontal dashed lines. This is
determined as the distance from the AGN position along the direction of the
jet position angle to the farthest position reached by the $3sigma$ radio flux
contour. We note that the broader line width observed visually is due to the
contrast between the high and low surface brightness parts. We find in some
cases (e.g. 4C+04.11 and TN J1338-1942) that there is a sharp surface
brightness drop at the distance of the radio jet boundary. A discontinuity in
the line width and surface brightness of the Ly$\alpha$ nebulae across the
extent of the radio source has been previously reported in Villar-Martín et
al. (2003). The changes is not as sharp as shown in 30. The sharpness seen in
the intrinsic panels is due to tessellation as it is the transition from high
to low surface brightness part. Hence, we assure that these are not all
artificial effect due to our analysis methods.

Table 6: Intrinsic surface brightness from circular aperture

R | MRC0943-242 | MRC0316-257 | TN J0205+2242 | TN J0121+1320 | 4C+03.24 | 4C+19.71 | TN J1338-1942 | 4C+04.11
---|---|---|---|---|---|---|---|---
$[\rm{arcsec}]$ | $(1+z)^{4}\rm{SB}\,[10^{-16}\,{\rm erg\,s^{-1}\,cm^{-2}\,arcsec^{-2}}]$
0.0-$r_{\rm s}$\rm{{\dagger}}\rm{{\dagger}}$\rm{{\dagger}}$$r_{\rm s}=0.75\times$ seeing which is specified for each target Table 1. The median seeing is 0.75 arcsec. | 2621.4 | 1534.1 | 1886.6 | 1452.9 | 2955.1 | 261.9 | 5177.9 | 6464.6
$r_{\rm s}$\rm{{\dagger}}\rm{{\dagger}}$\rm{{\dagger}}$$r_{\rm s}=0.75\times$ seeing which is specified for each target Table 1. The median seeing is 0.75 arcsec. -1.0 | 1792.0 | 833.8 | 1764.7 | 889.7 | 2102.8 | 224.9 | 3617.2 | 4686.2
1.0-1.8 | 623.3 | 390.1 | 953.4 | 250.2 | 1884.4 | 195.4 | 1829.5 | 2659.6
1.8-3.3 | 95.3 | 60.9 | 232.4 | 48.5 | 1020.3 | 130.1 | 385.0 | 950.4
3.3-4.5 | 24.6 | 36.4 | 120.3 | 38.6 | 465.3 | 105.4 | 106.4 | 367.7
4.5-6.0 | 15.2 | 31.5 | 97.2 | 25.4 | 227.4 | 68.3 | 58.4 | 229.6
6.0-7.3 | 14.7 | 32.2 | 87.1 | 0.0 | 194.8 | 51.3 | 29.2 | 168.9
7.3-8.9 | 13.1 | 35.9 | 45.0 | 0.0 | 142.9 | 43.0 | 25.0 | 97.6
8.9-10.8 | 11.9 | 38.5 | 23.3 | 0.0 | 117.3 | 25.1 | 26.3 | 65.5
10.8-13.2 | 13.0 | 37.2 | 23.3 | 0.0 | 82.8 | 20.7 | 33.9 | 58.3
13.2-16.0 | 0.0 | 38.3 | 23.3 | 0.0 | 65.0 | 0.0 | 33.9 | 31.2
16.0-19.5 | 0.0 | 33.0 | 0.0 | 0.0 | 45.5 | 0.0 | 17.0 | 0.0
19.5-23.8 | 0.0 | 14.5 | 0.0 | 0.0 | 0.0 | 0.0 | 29.2 | 0.0
23.8-28.9 | 0.0 | 12.3 | 0.0 | 0.0 | 0.0 | 0.0 | 36.6 | 0.0

111111
Intrinsic surface brightness values extracted from circular annuli following
Fig. 29. To be consistent, the first column gives the radii of each annulus in
the unit of arcsec instead of ckpc. The reported surface brightness values are
corrected for cosmological dimming.

Table 7: Fit parameters of surface brightness profiles.

HzRG |  $r_{\rm h}$\rm{{\dagger}}\rm{{\dagger}}$\rm{{\dagger}}$Scale length of the exponential profile. |  $C_{\rm e}$\rm{*}\rm{*}$\rm{*}$Normalisation parameter of the exponential profile. The cgs unit is $\mathrm{erg\,s^{-1}\,cm^{-2}\,arcsec^{-2}}$. |  $r_{\rm b}$\rm{{\ddagger}}\rm{{\ddagger}}$\rm{{\ddagger}}$Distance where inner exponential profile changes to power law profile. |  $\alpha$\rm{$\mathsection$}\rm{$\mathsection$}$\rm{$\mathsection$}$Power law index.
---|---|---|---|---
| pkpc |  $10^{-16}$ cgs | pkpc |
MRC 0943-242 | 6.9$\pm$0.2 | 20.7$\pm$1.0 | 39.8$\pm$1.5 | -0.24$\pm$0.18
” ” app. | 5.6$\pm$0.2 | 18.4$\pm$1.1 | 32.0$\pm$1.3 | -0.13$\pm$0.15
” ” rec. | 7.5$\pm$0.2 | 25.8$\pm$1.0 | 42.9$\pm$1.3 | -0.15$\pm$0.19
MRC 0316-257 | 7.0$\pm$0.6 | 8.9$\pm$1.2 | 27.3$\pm$3.2 | -0.62$\pm$0.24
” ” app. | 1.3$\pm$0.3 | 69.0$\pm$70.0 | 7.7$\pm$1.0 | -0.46$\pm$0.07
” ” rec. | 7.5$\pm$3.3 | 14.2$\pm$15.7 | 31.6$\pm$14.4 | -0.01$\pm$0.49
TN J0205+2242 | 48.9$\pm$30.4 | 5.7$\pm$0.4 | 10.9$\pm$21.4 | -1.90$\pm$0.58
” ” app. | 8.8$\pm$2.0 | 6.2$\pm$0.7 | 13.4$\pm$10.8 | -1.52$\pm$0.26
” ” rec. | 10.8$\pm$3.1 | 8.3$\pm$0.3 | 48.7$\pm$0.7 | -0.97$\pm$0.14
TN J0121+1320 | 9.8$\pm$1.9 | 5.6$\pm$1.2 | 6.2$\pm$10.0 | -2.05$\pm$0.71
” ” app. | 4.4$\pm$1.0 | 8.9$\pm$2.3 | 19.6$\pm$6.6 | -0.67$\pm$0.76
” ” rec. | 5.1$\pm$1.3 | 10.2$\pm$2.1 | 10.9$\pm$9.3 | -2.16$\pm$0.75
4C+03.24 | 17.6$\pm$1.8 | 7.9$\pm$0.6 | 20.0$\pm$14.9 | -1.64$\pm$0.21
” ” app. | 11.2$\pm$0.9 | 10.3$\pm$0.9 | 38.5$\pm$4.2 | -0.87$\pm$0.18
” ” rec. | 27.7$\pm$2.0 | 6.5$\pm$0.4 | 23.4$\pm$10.8 | -1.67$\pm$0.18
4C+19.71 | 29.5$\pm$7.2 | 0.8$\pm$0.1 | 23.1$\pm$9.7 | -1.20$\pm$0.16
” ” app. | 63.6$\pm$20.3 | 0.4$\pm$0.1 | 38.1$\pm$7.8 | -1.88$\pm$0.37
” ” rec. | 23.3$\pm$24.5 | 1.0$\pm$0.2 | 15.1$\pm$15.3 | -1.24$\pm$0.35
TN J1338-1942 | 7.3$\pm$0.3 | 13.6$\pm$0.7 | 41.0$\pm$2.7 | -0.12$\pm$0.26
” ” app. | 3.6$\pm$7.0 | 20.7$\pm$8.8 | 7.5$\pm$1.7 | -2.10$\pm$0.07
” ” rec. | 7.3$\pm$8.0 | 19.0$\pm$5.1 | 40.0$\pm$14.7 | -0.52$\pm$0.77
4C+04.11 | 8.2$\pm$1.5 | 11.8$\pm$1.3 | 15.2$\pm$11.7 | -2.04$\pm$0.24
” ” app. | 9.9$\pm$1.0 | 11.0$\pm$0.9 | 38.9$\pm$12.0 | -1.77$\pm$0.37
” ” rec. | 7.5$\pm$1.6 | 11.6$\pm$1.6 | 16.0$\pm$8.3 | -2.14$\pm$0.28
121212 See Section 4.2.3 for the fitting equations (Eq. 2):

##  Appendix E MRC0316-257 systemic redshift

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

For MRC0316-257, two velocity components of the He ii emission are detected in
our MUSE observation which are also separated spatially. The one detected at
the position of the X-ray emission peak from Chandra observation (Table 2) is
believed to be the systemic one while the one at north-west position that is
coincident with UV continuum emission peak may trace jet-gas interaction. We
show the UV continuum map of MRC0316-257 and the fits of the two He ii spectra
in Fig. 31. We note that the UV continuum map is constructed from the MUSE
cube using the wavelength in observed frame between N v$\lambda 1240$ and O
i+Si ii$\lambda 1305$ which is a emission-line-free region of HzRGs (McCarthy
1993). The bright continuum emission object east of the central AGN position
peak is a foreground galaxy (Vernet et al. 2017). We report here that the
systemic redshift detected is $z_{\rm sys}=3.1238\pm 0.0002$ and the redshift
of the component at the UV continuum peak is $z_{\rm red}=3.1323\pm 0.0002$
which is redshifted of $v=620\,\rm{km\,s^{-1}}$ from the systemic one. The
velocity gradient of the He ii agrees with the Ly$\alpha$ $v_{50}$ (Fig. 1c)
and [O iii]$\lambda 5007$ (Nesvadba et al. 2008) within the scope of the jet.
The UV continuum at this position may suggest the younger stellar population
distribution. Combine with the jet kinematics, we may seeing jet induced star
forming activities. However, there is also the possibility that the UV
continuum could be produced by the inverse Compton processes. The redshifted
He ii near the west radio lobe could then be due to the ionisation emission
from the shock region exerting on the un-shocked gas. This is supported by the
relatively narrow with of the redshift He ii (FWHM $\sim
600\,\mathrm{km\,s^{-1}}<1000\,\mathrm{km\,s^{-1}}$) which could indicate that
is has not been impacted by shocks (Best et al. 2000; Allen et al. 2008). We
note that a detailed verification for this scenario is beyond the scope of
this paper which involves spectral aging inspection of the radio jet hot-spot
(e.g. Harwood et al. 2013). Hence, we simply point these possibilities and
leave them to future study with multi-wavelength observations combined.