JWST’s PEARLS: TN J1338–1942 – I. Extreme jet-triggered star formation in a z  = 4.11 luminous radio galaxy Open Access

ABSTRACT

We present the first JWST observations of the z  = 4.11 luminous radio galaxy TN J1338–1942, obtained as part of the ‘Prime Extragalactic Areas for Reionization and Lensing Science’ (‘PEARLS’) project. Our NIRCam observations, designed to probe the key rest-frame optical continuum and emission line features at this redshift, enable resolved spectral energy distribution modelling that incorporates both a range of stellar population assumptions and radiative shock models. With an estimated stellar mass of log₁₀(M/M_⊙) ∼ 10.9, TN J1338–1942 is confirmed to be one of the most massive galaxies known at this epoch. Our observations also reveal extremely high equivalent-width nebular emission coincident with the luminous AGN jets that is best fit by radiative shocks surrounded by extensive recent star formation. We estimate the total star-formation rate (SFR) could be as high as |$\sim 1600\, \text{M}_{\odot }\, \text{yr}^{-1}$| , with the SFR that we attribute to the jet induced burst conservatively |$\gtrsim 500\, \text{M}_{\odot }\, \text{yr}^{-1}$| . The mass-weighted age of the star-formation, t_mass < 4 Myr, is consistent with the likely age of the jets responsible for the triggered activity and significantly younger than that measured in the core of the host galaxy. The extreme scale of the potential jet-triggered star-formation activity indicates the potential importance of positive AGN feedback in the earliest stages of massive galaxy formation, with our observations also illustrating the extraordinary prospects for detailed studies of high-redshift galaxies with JWST.

1 INTRODUCTION

The luminous high redshift radio galaxy (HzRG), TN J1338–1942 (TNJ1338 hereafter; De Breuck et al. 1999, z = 4.11), is one of the most powerful radio sources known in the early Universe. Initially selected based on its bright radio continuum emission and ultra-steep radio spectral index (De Breuck et al. 2000, S_ν ∝ ν^−1.31), TNJ1338 was subsequently discovered to reside at the heart of a significant overdensity of galaxies extending to Mpc scales (Venemans et al. 2002; Miley et al. 2004; Zirm et al. 2005; Intema et al. 2006; Overzier et al. 2008; Saito et al. 2015), making it one of the earliest known proto-cluster structures, and a potential progenitor of the most massive clusters (and brightest cluster galaxies) in the local Universe.

As host to extremely luminous extended Lyα emission (⁠|$L_{\text{Ly}\alpha } \sim 4\times 10^{44}\, \text{erg s}^{-1}$|⁠; De Breuck et al. 2000; Venemans et al. 2002; Swinbank et al. 2015), potential jet-induced star-formation activity (Miley et al. 2004; Zirm et al. 2005) and evidence of large active galactic nuclei (AGN)-driven outflows (Swinbank et al. 2015), TNJ1338 and its wider environment therefore represent a unique laboratory for studying both AGN feedback in action and the earliest stages of proto-cluster formation. Given these extreme properties and the wealth of multiwavelength observations available, TNJ1338 therefore represents a prime target for JWST observations, with Cycle 1 programmes in both imaging (GTO 1176, PI: Windhorst; Windhorst et al. 2022) and resolved spectroscopy (GO 1964, PIs: Overzier and Saxena) aiming to address a broad range of scientific questions.

One of the key outstanding challenges for radio-loud AGN such as TNJ1338 is understanding how their radio jets interact with their surrounding environments (see e.g. Fabian 2012; King & Pounds 2015, and references therein). For example, how, and to what degree, do the jets trigger or suppress star formation? One of the standout observational features of TNJ1338 is the presence of extended optical emission aligned with the most luminous radio lobe (Miley et al. 2004; Zirm et al. 2005). Based on its rest-frame UV emission, Zirm et al. (2005) estimated that non-stellar processes (scattered light, synchrotron, and inverse-Compton emission) can contribute only a few per cent of the rest-frame UV emission aligned with the radio jet; recent star-formation activity is therefore the dominant emission mechanism.

A primary objective of this study is therefore to resolve the nature of the extended emission aligned with the radio jet. Is star formation the dominant emission mechanism? If so, how does the age of this star-formation activity compare to the jet activity (and the accretion that powered this)? How significant is the triggered activity in the evolution of the host galaxy?

HzRGs are thought to host some of the most massive supermassive black holes (SMBHs) at a given epoch (Miley & De Breuck 2008). Crucially, unlike luminous quasars, the AGN in HzRG do not outshine their hosts, allowing detailed studies of the galaxies in which the most massive SMBH reside without being dazzled by the quasar light (e.g. Marshall et al. 2021). A secondary of this study is to exploit the high spatial resolution and sensitivity at rest-frame optical wavelengths now available from JWST to perform the first resolved study of the physical properties of a HzRG host galaxy, and begin placing high-redshift accretion activity into its wider context. Specifically, we want to advance beyond estimations of just the total stellar mass, establishing instead when and how the host galaxy formed.

Given these objectives, the filters for our JWST/NIRCam observations were selected to maximize the colour information available for sources at the systemic redshift of TNJ1338 (z = 4.11, derived from the observed He ii λ1640 emission line; Swinbank et al. 2015), with medium-band filters constraining key rest-frame optical features. Fortuitously, at z ∼ 4.1, the NIRCam SW filters at 1.5–2.1 μm bracket the stellar age-sensitive Balmer/D_4000Å; break, while the LW medium bands offer a clean measure of the rest-frame optical continuum and the strength of the Hα + N ii nebular emission line complex. With resolved constraints on both the older stellar population and recent star-forming activity respectively, NIRCam will enable us to probe the nature of TNJ1338 through detailed spectral energy distribution (SED) modelling in unprecedented detail for a galaxy at this redshift.

The rest of this paper is set out as follows: Section 2 outlines the NIRCam and accompanying ancillary photometric data used in this study, as well as the methodologies for data reduction and homogenization. In Section 3, we present a qualitative analysis of the rest-frame optical morphology of TNJ1338 revealed by NIRCam. Next, Section 4 outlines the methodology for our quantitative analysis of the resolved properties of TNJ1338 through SED modelling. Section 5 then presents the results of SED modelling analysis and our discussion. Finally, in Section 6, we summarize the results and our conclusions. Throughout this paper, all magnitudes are quoted in the AB system (Oke & Gunn 1983) unless otherwise stated. We also assume a Λ cold dark matter cosmology with H₀ = 70 km s⁻¹ Mpc⁻¹, Ω_m = 0.3 and Ω_Λ = 0.7 (corresponding to a scale of |$6.87\, \text{kpc}$| per″).

2 DATA

2.1 Optical imaging

The Hubble Space Telescope (HST) Advanced Camera for Surveys (ACS) observations of the TNJ1338 field used in our analysis comprise 18 orbits split across four broad-band filters: 9400 s in F475W and F625W, 11 700 s in F775W, and 11 800 s in F850LP (Proposal 9291; Miley et al. 2004; Zirm et al. 2005). These four filters span the Lyman break and rest-frame UV at z ∼ 4.1, with the F625W-band probing the Lyα emission line. The reduced ACS observations were retrieved from the Hubble Legacy Archive (HLA),¹ with processing using the standard HLA reduction pipelines, including appropriate bias and dark current subtraction, flat-fielding, and geometric distortion corrections. In the HLA processing, the exposures for each band have been drizzled on to a consistent 0.04″ pixel scale.

In addition to the HST/ACS observations, TNJ1338 has also been observed with ground-based multiband optical imaging extending over a wider field of view (FOV) using the Subaru/Suprime-Cam imager (B, R, i′, and IA624 narrow-band; Intema et al. 2006; Saito et al. 2015). We retrieve the raw Suprime images from the SMOKA archive (Baba et al. 2002) and produce a stacked mosaic with a custom pipeline that includes astrometric corrections performed using Gaia EDR3 objects (Gaia Collaboration et al. 2021), taking into account their proper motion to the epoch of Suprime observations. The i′-band image, reaching a 5σ limiting magnitude of 26.56 is used to provide a reference catalogue for astrometric alignment of the JWST/NIRCam imaging as outlined below. Due to the increased depth and resolution afforded by the HST/ACS imaging, we do not include the Subaru ground-based imaging in our analysis of TNJ1338 directly.

2.2 NIRCam imaging

The JWST/NIRCam observations of the TNJ1338 field were taken as part of the ‘Prime Extragalactic Areas for Reionization and Lensing Science’ (PEARLS; Windhorst et al. 2022) Guaranteed Time Observations (GTO) programme. These observations, taken on 2022 July 1, consist of imaging in the F150W, F182M, F210M, F300M, F335M, and F360M filters with a uniform integration time of 1030.73 s in each filter (SHALLOW4 readout pattern and 4-point INTRAMODULEBOX dithers). For a full description of the data reduction steps, we refer the reader to Windhorst et al. (2022). Here we summarize the key steps, including those specific to the processing of the TNJ1338 imaging.

Initial calibration and processing is done using the standard STScI CALWEBB.² All analyses and images presented in this paper are based on the latest reductions using Pipeline version 1.7.2 and the context file ‘jwst_0995.pmap_filters’ made available in October 2022, which incorporates the on-orbit flat-fielding observations and improved flux calibrations. Detector level offsets in the SW filters and 1/f-noise in both SW and LW images are removed in individual frames using the ProFound code (Robotham et al. 2017, 2018). We also remove stray light ‘wisps’ in the SW imaging by subtracting a wisp template that best matches the amplitude observed in that specific image.

Before mosaicing, it is essential that individual frames are placed on a consistent and reliable astrometric frame. However, due to the unexpected early scheduling of the TNJ1338 visit,³ the observed roll angle of the NIRCam imaging resulted in reduced overlap with the existing HST/ACS observations described above. Image registration for the NIRCam frames is therefore done using the wide area Subaru imaging and associated catalogues on the Gaia DR3 reference frame (Gaia Collaboration et al. 2021, 2022), specifically the i′-band, which has the greatest sensitivity for redder objects. Positional offsets (adjusted for proper motion) between the NIRCam images and the reference catalogue are used to adjust each frame’s centre and positional angle until the offsets are minimized. The aligned images are then drizzled on to mosaics with a 0.03″ pixel scale using the AstroDrizzle package (Koekemoer et al. 2013; Avila et al. 2015).

As measured by Windhorst et al. (2022), the resulting mosaics reach 5σ limiting magnitudes of 27.7, 27.4, and 27.2 in the SW bands (F150W, F182M, and F210M, respectively), and 28.35, 28.25, and 28.16 in the LW bands (F300M, F335M, and F360M, respectively). In line with other observations, the measured sensitivities represent a small reduction in the SW bands when compared to pre-launch Exposure Time Calculations (∼0.1 mag). However, the measured LW bands represent a substantial improvement on expectations, typically reaching ∼1 mag deeper.

Finally, using the F335M NIRCam mosaic as our reference, we use the reproject package to place the ACS observations on to the unified pixel grid defined by the NIRCam mosaics. Comparing the positional offsets of sources detected in both the ACS and NIRCam imaging, we confirm agreement between the reprojected images and the drizzled mosaics is at the sub-pixel level. The combined HST/ACS and JWST/NIRCam filter coverage and corresponding limiting sensitivity is illustrated in Fig. 1, as well as an illustrative galaxy SED at z = 4.11 – highlighting both the key spectral features probed by the NIRCam bands and the unprecedented sensitivity at 1–4 μm.

2.3 Additional ancillary data

Due to its extreme nature, TNJ1338 and its surrounding protocluster environment have also been observed across the full electromagnetic spectrum. Following its initial detection in the low-frequency radio (De Breuck et al. 2000), TNJ1338 has also been observed with higher frequency radio continuum observations. In this study, we primarily make use of the highest available resolution imaging provided by the Karl G. Jansky Very Large Array (VLA) X-band observations presented by Pentericci et al. (2000). The observations are centreed at a frequency of 8.46 GHz, reaching an RMS sensitivity of 25|$\mathrm{ \mu} \text{Jybeam}^{-1}$|⁠, and have a spatial resolution of 0.23″, sufficient to resolve the jet structure on the scale of the observed optical emission (Zirm et al. 2005).

Additional observations spanning the infrared include 24 μm Spitzer/MIPS (De Breuck et al. 2010), 100–500 μm Herschel/PACS and SPIRE (Drouart et al. 2016), 450/850 μm JCMT/SCUBA (De Breuck et al. 2004), 1.2 mm IRAM/MAMBO observations (De Breuck et al. 2004), and ALMA Band-3 (Falkendal et al. 2019). However, none of these observations have sufficient spatial resolution to incorporate within any resolved analysis on the properties of TNJ1338. Never the less, together they can still place valuable constraints on the integrated IR emission from TNJ1338.

2.4 Point-spread function homogenization

To ensure accurate pixel by pixel colours and enable studies of the resolved SED, we homogenize all of the available ACS and NIRCam imaging to a common point-spread function (PSF). Empirical PSFs are first generated from stacked stars within the field, with the initial selection of bright stars from Gaia (Gaia Collaboration et al. 2022) visually inspected to exclude saturated sources, stars with close neighbours and stars near the edge of the field with incomplete coverage within the cutouts. Before generating the final stacked PSF, we use the Photutils.ePSF package (Bradley et al. 2016, 2022) to iteratively re-centre each star using a pixel oversampling of four. Due to the small number of suitable stars available in the field, we found that the full Photutils.ePSF effective PSF modelling approach was not able to accurately model the full structure of the JWST/NIRCam PSF. We therefore generate the final PSFs from a simple median stack of the re-centred star images with each star first normalized by the flux contained within a 0.18″ diameter aperture.

To derive the convolution kernels required to transform the measured PSFs to a given target PSF we employ the pypher package (Boucaud et al. 2016). Nominally, the F360M filter exhibits the broadest PSF full width half-maximum (FWHM) within our filter set and would be a suitable target PSF for JWST/NIRCam-only studies. However, we find that that difference in structure between the wings of the HST/ACS and JWST/NIRCam PSFs can result in excess noise within the convolution kernel. Our target PSF is therefore defined to be a Moffat profile (f(r) ∝ (1 + (r/α)²)^−β) with β = 2 and α = 3.63 px, resulting in a FWHM = 0.14″ fractionally larger than that measured for our largest empirical PSF (F360M; 0.13″).

Measuring the curves of growth in circular apertures for all HST/ACS and JWST/NIRCam bands after PSF homogenization, we find that the PSFs of the convolved images agree to within |$\lt 5~{{\ \rm per\ cent}}$| at all radii (see Fig. 2).

2.5 Correcting uncertainties for correlated noise

While both the HST and JWST imaging allow calculation of uncertainties within an aperture or region based on the noise maps, it is well-understood that uncertainties derived by common photometry tools can be underestimated in the presence of correlated noise (Leauthaud et al. 2007). Since the convolution of the drizzled images with a PSF kernel will further increase any existed correlation between pixels, uncertainties derived from the HST/ACS or JWST/NIRCam error maps can therefore underestimate the total uncertainty for larger regions (i.e. apertures, isophotes, or source segments).

We estimate the corrections required to account for correlated noise in our subsequent photometry steps following an approach similar to other works in the literature (Bielby et al. 2012; Laigle et al. 2016; Mehta et al. 2018, see also Trenti et al. 2011 for an alternative approach). First, for each band we derive an initial source catalogue from the PSF-homogenized image, requiring a 5σ detection threshold with a minimum of 10 connected pixels. We then define 10⁴ random empty positions within the image, where the segmentation map derived from the detected sources is first dilated to ensure that the selected empty sky regions do not fall within extended sources. The correlated noise correction for a given aperture size is then estimated by comparing the 3σ-clipped scatter measured in empty sky apertures with the median uncertainty measured for sources using the same size aperture. Since our resolved analysis of TNJ1338 will employ flux measurements over a range of different photometric aperture or segment sizes, we follow the approach of Mehta et al. (2018) and derive the correlated noise uncertainty corrections as a function of aperture radius (area). The resulting measurements are shown in Fig. 3. For an aperture with radius r = 0.14″ (∼5 pixels; i.e. the FWHM of the PSF homogenized images), we find that the uncertainties are underestimated by a factor of 2−5, with a filter dependence in line with naive expectations (i.e. that bands with greater oversampling with respect to their native resolution require greater corrections).

In our subsequent analysis, we calculate the required uncertainty correction for a given segment, isophotal or Kron flux measurement by interpolating the correlated uncertainty measurements based on the segment area. We note that due to the extreme rest-frame optical luminosity of TNJ1338 and the sensitivity of the JWST observations, the limiting uncertainty on fluxes and derived colours is the remaining systematic uncertainty in the zeropoints between bands (see Appendix B1 of Windhorst et al. 2022). Never the less, since these zero-point uncertainties are systematic within a given detector and will apply uniformly to TNJ1338, accounting for the correlated noise during SED fitting analysis is necessary to ensure that resolved studies correctly account for signal to noise variations across the host galaxy.

3 JWST’S VIEW OF TN J1338-1942

The combination of the physical extent of TNJ1338 with the high spatial resolution and sensitivity of JWST enables us to conduct the first spatially resolved study of the rest-frame optical properties of an AGN host galaxy in the early Universe. Fig. 4 shows the combined PSF-homogenized imaging of TNJ1338, from the rest-frame UV probed by HST/ACS to the rest-frame optical probed by the NIRCam/LW filters.

A number of key features are immediately revealed by the NIRCam imaging. First, the JWST observations confirm the clear presence of a bright galaxy core located at the Southern side of the source. Although present in VLT/ISAAC K_S broad-band imaging (Zirm et al. 2005), the ground-based imaging was not able to fully resolve the morphology of the rest-frame optical emission, or sensitive enough to rule out more extensive obscured emission.

In the filters probing rest-frame wavelengths above the Balmer break (F210M to F360M; see Fig. 1) there are significant differences in colour for the galaxy core and jet-aligned emission (where the jet emission is illustrated in the Fig. 4). The extended rest-frame optical emission is visibly clumpy, with some of the distinct regions previously identified in Lyα visible more clearly in the F182M and F335M, filters that are dominated by the bright [O ii] and Hα emission lines, respectively (see annotations in Fig. 4).

To illustrate the variation and amplitude of the optical emission line strengths across TNJ1338, in Fig. 5 we show the observed F300M−F335M colour for the individual pixels with significant detections in both filters (>10 × the local RMS background in both images). To convert the measured colours into a physically meaningfully quantity independent of assumptions on the associated emission mechanism, we also estimate the |$\text{H}\alpha +\text{N}\,\small {II}$| emission line equivalent widths (EW) implied by the observed colours. For the simplifying assumption of a flat intrinsic optical continuum, where the emission line-free colour would result in F300M−F335M = 0, the EW of the combined emission line complex can be estimated as:

where Δ_mag is the change in magnitude in the F335M filter attributed to the |$\text{H}\alpha +\text{N}\,\small {II}$| emission line complex (see e.g. Mármol-Queraltó et al. 2016, with Δ_mag = F335M− F300M for our assumption), and the effective width, W_eff, for the F335M filter is taken to be 3389.42 Å. We caution that the resulting |$\text{EW}_{\text{H}\alpha +\text{N}\,\small {II}}$| estimates are only approximate, however, as is illustrated later in the paper (see Figs 9 and A1), the assumption of a flat intrinsic continuum is well-motivated for much of the extended emission.

There is clear evidence for a strong variation in |$\text{EW}_{\text{H}\alpha +\text{N}\,\small {II}}$| along the axis between the centre of the host galaxy (shown by the red plus) and the northern radio lobe, with the peak |$\text{EW}_{\text{H}\alpha +\text{N}\,\small {II}}$| (blue cross) situated extremely close to the peak of the 8.46-GHz continuum emission. Fig. 5 also reveals that there is also an increase in |$\text{EW}_{\text{H}\alpha +\text{N}\,\small {II}}$| in the direction of the knots of high-frequency radio emission in the south east, albeit with substantially lower amplitude than the northern lobe. The region of higher |$\text{EW}_{\text{H}\alpha +\text{N}\,\small {II}}$| emission south-east of the galaxy core is visible as a knot of bluer rest-UV emission in the original HST/ACS imaging and NIRCam F150W. The alignment of this higher |$\text{EW}_{\text{H}\alpha +\text{N}\,\small {II}}$| region relative to the southern radio lobe and the centre of the rest-frame optical continuum add additional evidence for jet activity driving the nebular line emission. However, given the systematic offset in astrometry for the VLA observations could be as large as 0.3″(Zirm et al. 2005), the original interpretation of the peak in 8.46 GHz emission as a radio core remains plausible (we note that such an offset would also place the peak of the northern lobe closer to the peak in |$\text{EW}_{\text{H}\alpha +\text{N}\,\small {II}}$|⁠).

Another noteworthy morphological feature for which the physical interpretation is informed by the confirmed galaxy position is the extended cone of Lyα emission, previously dubbed ‘the Wedge’ by Zirm et al. (2005). This feature, visible in the HST/ACS F625W image in Fig. 4 (with dashed white lines illustrating the approximate bounds), was hypothesized to be the result of starburst-driven outflows emanating from the galaxy. However, we see no equivalent feature in the F335M filter that directly probes Hα, either in the raw imaging or when smoothing to enhance low surface-brightness emission. Furthermore, the focal point of the wedge (see annotation in Fig. 4) does not coincide with the confirmed location of the host galaxy and so is unlikely to be directly driven by a central starburst. We therefore interpret this emission to be a feature of the extended Lyα halo surrounding TNJ1338. However, given the observed extent of the Lyα halo (>100 kpc; Venemans et al. 2002), its complex kinematics (Swinbank et al. 2015) and the lack of additional information for this feature provided by the NIRCam imaging, we do not include any further analysis of ‘the Wedge’ in this study.

4 RESOLVED SED MODELLING

While the variation in observed colours seen in Fig. 5 indicates the presence of significant nebular contribution to TNJ1338’s rest-frame optical emission, more quantitative analysis of the full SED is necessary to place constraints on the nature of the emission and how it relates to the properties of the host galaxy. We therefore perform detailed SED modelling for individual regions within the galaxy.

4.1 Segmentation photometry

Given the irregular, clumpy, and strongly wavelength dependent morphology seen in Fig. 4, reliably decomposing the observed emission with parametric morphological models is not optimal. We therefore use a multistep process to segment the TNJ1338 into relevant physical regions for SED analysis. Our photometry methodology is intended to balance maximizing the spatial information available with the optimal signal-to-noise ratio (SNR) necessary to derive detailed stellar population properties with multiband photometry (see e.g. Leja et al. 2019).

We use the F335M medium band as our detection and segmentation image due to the presence of the Hα emission line within the band-pass, allowing extremely high SNR identification of multiple distinct clumps within the host galaxy, as well as the luminous galaxy core. Segmentation of the galaxy into these regions of interest is done in multiple stages using the Photutils package (v.1.5; Bradley et al. 2016, 2022). First, the main galaxy core is identified using a detection threshold of >120 × the local background RMS. After masking the full core, the peak of the jet-aligned emission and the galaxy outer core are selected with a lower noise threshold of >43 × local RMS (where the threshold was iteratively selected to ensure that the two areas were maximized without becoming connected). We further split the outer core in two relative to the peak of the radio continuum emission to allow us to compare the properties of the northern and southern edges of the host galaxy and the variation in physical properties relative to the radio jet. Similarly, to allow a more detailed study of how the SED varies close to the northern radio lobe, we split the secondary peak into four segments by first de-blending into two components and then splitting each of those in two relative to the peak of the radio emission.

Finally, the remaining galaxy emission is split into bins following the optimized Voronoi binning procedure of Cappellari & Copin (2003). Starting at a high target SNR per bin of >250, we iteratively apply the Cappellari & Copin (2003) algorithm, lowering the target SNR until the smallest Voronoi bin falls below an area of 50 pixels. The starting target SNR is set to this high threshold based on two factors: first, the Voronoi binning optimization only makes use of the noise provided by the noise map, which does not account for the correlated uncertainties as outlined in Section 2.5, and secondly the fact that the emission in F335M is substantially brighter than the others and by setting a high threshold we aim to ensure that the final SNR accounting for correlated uncertainties remains significant in all bands.The final multiband photometry is performed by constructing a combined segmentation map containing all of the regions defined above; segmentation photometry is then performed on each of the PSF-homogenized HST/ACS and JWST/NIRCam images. Uncertainties measured for each segment are corrected for correlated noise based on the segment area and the corresponding calculations outlined in Section 2.5. In Fig. 6, we illustrate the resulting image segments compared to the galaxy emission. Also shown in Fig. 6 is the mean signal to noise in each segment across all bands used in the SED fitting analysis (⁠|$\overline{\text{SNR}}$|⁠). The final photometry catalogue consists of 19 segments, with the mean SNR across the observed filters ranging from ∼7.5 in F775W to ∼105 in F335M.

4.2 SED modelling

4.2.1 Stellar population models

Modelling of the resolved TNJ1338 SED is performed using the Prospector Bayesian SED modelling code (Johnson et al. 2021b), which incorporates a full range of stellar population models from the Flexible Stellar Population Synthesis (FSPS; Conroy, Gunn & White 2009; Conroy & Gunn 2010) package (and the python-FSPS bindings; Johnson et al. 2021a). Nebular line and continuum emission associated with the stellar population are implemented in Prospector following the method outlined in Byler et al. (2017), which uses the photoionization models of Cloudy (Ferland et al. 2013) to self-consistently model the total nebular emission for the ionizing spectra of the FSPS stellar populations.

One of the key advantages of Prospector is its implementation of non-parametric star-formation histories (SFH), which have been demonstrated to provide more accurate constraints on the SFH derived from both photometric and spectroscopic measurements of galaxy SEDs, producing less biased measurements and more accurate estimates of the corresponding uncertainties (Leja et al. 2019). However, as highlighted by Johnson et al. (2021b) and Leja et al. (2019), the assumed priors can have a strong impact on the inferred SFH, and hence the corresponding stellar mass estimate. The impact of this prior is maximized in the presence of very young stellar populations (<30 Myr) due to the outshining of older stars, and also at higher redshifts where constraints on the rest-frame near-infrared stellar emission are not available (Tacchella et al. 2022; Topping et al. 2022; Whitler et al. 2023a, b).

From Fig. 5 and additional inspection of the measured SEDs in a number of regions within TNJ1338, there is clear evidence for extremely high |$\text{EW}_{\text{H}\alpha +\text{N}\,\small {II}}$| emission, potentially indicative of a very young stellar population. For our analysis, we therefore follow the approach of recent studies (see e.g. Endsley et al. 2022) and assume two sets of SFH parameters designed to bracket the potential extremes in stellar mass. First, we assume a non-parametric SFH using the ‘continuity’ prior, with five age bins:

Secondly, we assume a constant SFH (CSFH) where the time since the onset of SF, t_age, is allowed to vary between 1 Myr and 1.29 Gyr (i.e. a formation redshift of z = 20) with a uniform prior in log₁₀(t_age/Gyr).

For both types of SFH, we use a Kroupa (2001) initial mass function (IMF). Stellar metallicity is allowed to vary in the range −2 < log₁₀(Z/Z_⊙) < 0.19 and the gas-phase metallicity (for nebular emission) is fixed to that of the stellar metallicity. Additionally, for the nebular emission component of the FSPS models (Byler et al. 2017), we include both line and continuum emission, allowing the ionization parameter, U, to vary with a log-uniform prior in the range −4 < log₁₀(U) < −1. We choose to fit log₁₀(U) as a free parameter after initial tests demonstrated that lower log₁₀(U) values than the assumed default (= −2.5) were necessary to reproduce the F182M observations in many regions. Dust attenuation is assumed to be simple foreground attenuation with varying optical depth, assuming either Calzetti et al. (2000) or Small Magellanic Cloud (SMC; Gordon et al. 2003) dust attenuation laws. We allow the optical depth at 5500 Å to vary in the range 0 < τ_5500Å; < 3 with a uniform prior. Limitations and caveats of our SED fitting model assumptions will be discussed further in Section 5 within the context of the inferred physical properties of TNJ1338.

When fitting the observed photometry, we include an additional 5  per cent fractional uncertainty (added in quadrature to the measured uncertainties) to account for the estimated remaining systematic uncertainty in the NIRCam photometric zeropoints (Boyer et al. 2022). Because of the significant potential for scattering of Lyα photons into and out of the large extended Lyα emission around TNJ1338 (Venemans et al. 2002; Swinbank et al. 2015), we do not include the ACS F625W filter in the SED fitting analysis (however, we show the observed and predicted magnitudes in subsequent figures for reference).

4.2.2 Radiative shock models

In a number of powerful HzRGs, the jet-aligned emission has been attributed to radiative shocks (Best, Röttgering & Longair 2000; Inskip et al. 2002; Moy & Rocca-Volmerange 2002, in particular, in AGN with jets <150 kpc such as TNJ1338). To explore the potential for this mechanism driving the high-EW emission in TNJ1338, we also compare the observed photometry to predictions from the MAPPINGS III Library of Fast Radiative Shock Models (Allen et al. 2008). MAPPINGS III contains model predictions for radiative shocks and associated photoionized precursor for five different abundance patterns: SMC, LMC, Dopita et al. (2005), Solar, and 2 × Solar, with shock velocities ranging from 100 to 1000 km s⁻¹.

To fit the MAPPINGS III library to the observed photometry, we convolve the shock and shock + precursor SEDs with the HST/ACS and JWST/NIRCam filter profiles at the fixed redshift of z = 4.11 to obtain the predicted colours. The optimal scaling for each model to match the observed photometry from F775W to F360M for a given segment of TNJ1338 is then calculated using the standard analytic optimization (see e.g. equation (4) of Duncan et al. 2019), with the best-fitting model identified through χ² minimization of the full model set. For consistency with stellar population modelling above, we include an additional 5  per cent fractional uncertainty added in quadrature to the observed photometric uncertainties. We note that following Allen et al. (2008), the radiative shock models are assumed to be effectively dust free and no additional attenuation is applied to the shock model predictions.

5 RESOLVED PHYSICAL PROPERTIES OF TNJ1338

In total, our resolved SED modelling incorporates five different sets of assumptions (two sets of SFH each with two types of dust attenuation, plus one set of radiative shock models). To summarize the results of the SED fitting analysis, Fig. 7 illustrates the best-fitting model type for each region within TNJ1338. Colour shading shows which SFH assumption or shock model produces the best-fit for that region when assuming an SMC dust law, based on the χ² of the maximum a posteriori (MAP) value from Prospector or the minimum-χ² for the radiative shock models. Since the number of free parameters varies significantly between the assumed models, we also calculate the corresponding Bayesian Information Criterion (BIC) value for each model fit. Different colours around the edges of regions illustrate where the best model either changes for a Calzetti dust law or there is no significant evidence for the best model over the second best, i.e. where the difference in BIC between the two models is ΔBIC < 4, with the edge colour corresponding to the alternative possible model.

Supporting the qualitative picture presented in Section 3, Fig. 7 demonstrates that there are distinct regions within TNJ1338. The source regions closest to the host galaxy in the southern half of TNJ1338 show evidence for a more extended and complex SFH, marginally favouring the non-parametric SFH assumption. The SEDs of the region closest to the jet are instead more consistent with the radiative shock models or an extremely recent burst of SF activity. We note, however, that only one region (in which the peak of the 8.46 GHz radio continuum emission is located) strongly favours the radiative shock model regardless of the SFH or dust models assumed for the stellar population model (see Section 5.2 below). Regardless of the source of ionizing photons, the SEDs in this region are clearly dominated by strong nebular emission lines.

Fig. 8 presents two sets of resolved maps of key physical properties inferred from the stellar population modelling. The maps further highlight the distinction between the northern and southern halves of TNJ1338, with the stellar mass surface density, |$\Sigma _{\mbox{M}_{\star }}$|⁠, for the southern half of TNJ1338 substantially greater (by up to two dex) than the jet region, with a clear concentration around the galaxy core. A similar picture is seen in the inferred mass-weighted ages, |$t_{mass}$|⁠, where there is a clear dichotomy across the galaxy. Crucially, we observe the same spatial variation across the galaxy regardless of the assumed SFH or dust models. However, the exact normalization of some properties (e.g. Σ_SFR or |$t_{mass}$|⁠) do systematically vary with the assumed dust attenuation law. In the following subsections, we explore the details of Figs 7 and 8 and the corresponding properties of the jet-aligned emission and host galaxy in greater detail.

5.1 Nature of the jet-aligned emission

In their previous study of TNJ1338, Zirm et al. (2005) attribute the extended optical emission (rest-frame UV) to SF activity triggered by the radio jet. As illustrated in Fig. 7, the NIRCam observations reveal a more complex picture with multiple potential emission mechanisms. Resolved spectroscopic observations are required to fully disentangle the emission-line contributions from star formation, shocks, and AGN photoionization. Never the less, the resolved photometric SEDs are consistent with TNJ1338 being a massive galaxy with radiatively driven shocks and extensive recent star formation spatially correlated with the luminous radio jets.

In Fig. 9, we present the observed SED of the two regions closest to the peak of radio continuum emission and the adjacent region beyond the peak. For each region, we show the corresponding best-fitting stellar population (green lines and triangles), with the best-fitting radiative shock models (yellow dash–dotted line and diamonds) shown if consistent with the observations.

Because of the wide dynamic range present within the observed SEDs, visually identifying the colours that are driving the statistical preference for radiative shocks within these regions is non-trivial. Fig. 10 therefore illustrates the observed distribution of F210M− F300M and F300M−F335M colours across TNJ1338 compared to illustrations of the parameter space covered by the radiative shocks and stellar population models. Symbols for each region are coloured following the same convention as in Fig. 7. We can see that F210M−F300M is the colour driving the preference for shocks over SF models in the most extreme regions. Although a young stellar population can reproduce the most extreme observed F300M−F335M colour, the corresponding SEDs cannot produce sufficiently blue F210M−F300M colours even in the absence of dust. Without detailed measurements of the emission line ratios, we cannot make robust constraints on the properties of the shock itself. However, we note that the best-fitting shock velocities for the regions shown in Fig. 9, 450 and 375 km s⁻¹, are consistent with range of shock velocities measured in the best local analogues of TNJ1338 (Cresci et al. 2015; Capetti et al. 2022).

Beyond the shock-dominated area, the galaxy regions outside of the galaxy core have SEDs consistent with extremely young stellar populations, with CSFH mass-weighted ages |$t_{mass} \lt 4$| Myr (log₁₀(t_mass/yr) < 6.6; Fig. 8). Although we cannot definitively demonstrate a direct causal link between the radio jet and the recent starburst activity, the measurements presented in this study do provide additional support for the scenario of jet-triggered activity. First, our photometric observations of TNJ1338 provide a qualitatively similar picture to that observed in local known examples of jet-triggered SF. In the luminous radio galaxy Coma A, resolved spectroscopic observations (Capetti et al. 2022) show a region of highly ionized gas extending along the jet axis with emission lines indicative of radiative shocks surrounded by young stellar populations. Similarly, in simulations of the jet-triggered SF observed in Centaurus A (Crockett et al. 2012), Gardner et al. (2017) demonstrate that the induced star-formation activity occurs in the transverse bow-shock of the radio jet. Such a scenario could explain the picture seen in Fig. 7, where the region beyond the shock has an SED consistent only with a young recent burst (see also Bicknell et al. 2000).

Secondly, the tight constraints on the mass-weighted age of the recent starburst spatially correlated with the jet (⁠|$t_{mass} \lt 4$| Myr) remain fully consistent with the expected time-scale of that jet activity (and the local analogues such as Centaurus A; Crockett et al. 2012). Comprehensively modelling the jet lifetime given the observed radio spectral information and losses from inverse-Compton scattering is beyond the scope of this work. However, for radio jets with projected sizes similar to that seen in TNJ1338 (∼7 kpc from galaxy core to the peak of the radio emission) and comparable luminosities, inferred jet ages range from ∼1 to >10 Myr (see e.g. fig. 14 of Hardcastle 2018), consistent with the observed burst of SF.

Where TNJ1338 differs compared to examples of jet-induced SF activity at lower redshift (Best, Longair & Rottgering 1997; Croft et al. 2006; Crockett et al. 2012) is the scale of the star formation potentially triggered by the jet. In Fig. 11, we present the combined SED fitting posteriors on the total stellar mass and star formation measured in TNJ1338 (see also Table 1). To produce the combined posteriors we make 3000 random draws from the Prospector nested sampling chain for each region, summing the total stellar mass and star-formation rate (SFR) across the subset of regions of interest. For regions associated with the host galaxy (blue shading in Fig. 7) we use the non-parametric SFH assumption, while for regions in the jet-aligned emission (green or yellow shaded regions) we use the fits from the CSFH assumption.

If we assume all regions are dominated by stellar emission (dark blue and red posteriors and histograms), the total estimated SFR amounts to |$\gtrsim 1000\, \text{M}_{\odot }\, \text{yr}^{-1}$|⁠, specifically |$\log _{10}(\text{SFR}/\text{M}_{\odot }\, \text{yr}^{-1}) = 3.21^{+0.06}_{-0.04}$| for Calzetti dust and |$\log _{10}(\text{SFR}/\text{M}_{\odot }\, \text{yr}^{-1}) = 3.02^{+0.07}_{-0.06}$| for SMC (see Table 1).⁴ Even excluding any region that could be dominated by radiative shock emission instead of SF (yellow shaded or bordered areas in Fig. 7) the total estimated SFR is reduced only by 0.18 dex, resulting in rates of |$\sim 700-1000\, \text{M}_{\odot }\, \text{yr}^{-1}$|⁠.The SFR estimated from the host galaxy regions alone amounts to only 30  per cent of the total SFR (42  per cent excluding shock regions) for the Calzetti et al. (2000) dust assumption (18 and 27  per cent for SMC), meaning that the on-going jet-triggered star-formation could be at a rate of at least |$\sim 500\, \text{M}_{\odot }\, \text{yr}^{-1}$|⁠, independent of the assumed dust attenuation law.

As noted in Section 2.3, a number of FIR to sub-mm observations have been taken of TNJ1338 that provide additional constraints on the total ongoing SFR. To explore whether our inferred physical properties are consistent with the existing long wavelength observations, in Fig. 12 we show the mid-IR to sub-mm SED predicted by the model assumptions that yield the most extreme estimated SFR. For the Draine & Li (2007) IR model implemented in Prospector, we can broadly reproduce the observed IR SED when assuming a minimum interstellar radiation field, |$\mathcal {U}_{\mbox{min}} = 10$| (in units of the Milky Way radiation field).⁵ The photometric observation that cannot be reproduced by the ongoing SF activity is the 24-μm measurement, with the observed flux for TNJ1338 over 1 magnitude greater than predicted by the SF emission. Fitting only the integrated IR SED, Falkendal et al. (2019) attribute the 24-μm emission entirely to hot dust associated with AGN activity (not included in our predicted SED). As our NIRCam observations do not probe the rest-frame near to mid-IR wavelengths where emission from an AGN IR power law can dominate the observed galaxy SED (Donley et al. 2012), additional resolved mid-IR observations (i.e. JWST/MIRI) are required before placing further constraints on the obscured accretion activity.

Given the large systematic uncertainties in robustly constraining dust properties at high redshift (see e.g. Hodge & da Cunha 2020, for a recent review), we cannot claim that the FIR observations provide additional evidence for the extreme SFRs estimated from our resolved SED modelling. Fig. 12 does, however, illustrate that a jet-triggered starburst with |$\text{SFR}\gtrsim 1000\, \text{M}_{\odot }\, \text{yr}^{-1}$| can be fully consistent with the existing FIR emission.

Although extreme relative to low redshift examples, our estimates for the induced SF in TNJ1338 are comparable to those measured in the other known high-redshift examples of jet-triggered star-formation, 4C 41.17 (z = 3.8; Dey et al. 1997; Steinbring 2014), which has a total SFR of |$\sim 650\, \text{M}_{\odot }\, \text{yr}^{-1}$| based on resolved mm observations (Nesvadba et al. 2020). Furthermore, as highlighted in Fig. 11 (and implicit in the stellar mass estimates shown in Table 1), despite the extremely high ongoing star formation observed, the very short lifetime of the recent burst (⁠|$t_{mass} \lt 4$| Myr) means that the total stellar mass added to the galaxy by the jet-triggered activity is negligible compared to the mass in situ. It therefore remains unclear whether the potential positive feedback from luminous radio jets such as those in TNJ1338 play a significant role in shaping the evolution of their host galaxies.

5.1.1 Limitations and systematic uncertainties on stellar population properties

Although the choice of models and associated priors used in our SED fitting analysis are well-motivated for the available photometric information and source redshift, they never the less represent an implicit prior on the resulting inferred physical properties. A number of known systematic uncertainties could potentially alter our interpretation of the mechanisms powering the jet-aligned emission.

First, the assumed Kroupa (2001) IMF implemented in FSPS is limited to |$\lt 120\, \mbox{M}_{\odot }$| mass stars, while recent evidence suggests that stars as massive as |$300\, \mbox{M}_{\odot }$| exist in local star-forming complexes (Crowther et al. 2010). Secondly, the stellar isochrones and spectral libraries we assume (MIST; Choi et al. 2016, and MILES; Sánchez-Blázquez et al. 2006, respectively) do not include the impact of interacting binaries on the resulting SED. Including both higher mass stars and the effects of binary interaction leads to both bluer intrinsic stellar continuum SEDs and an increase in the production of ionizing photons at a given population age (Stanway, Eldridge & Becker 2016).

The use of alternative stellar population models such as the Binary Population and Spectral Synthesis code (BPASS; Eldridge et al. 2017), which incorporates both of these changes, could therefore lead to older inferred mass-weighted ages, higher permitted |$\text{EW}_{\text{H}\alpha +\text{N}\,\small {II}}$| and allow for bluer rest-frame optical colours (potentially leading to SF models compatible with the SED currently only fit by radiative shock models). The amount of recent SF required to produce the observed |$\text{EW}_{\text{H}\alpha +\text{N}\,\small {II}}$| could also be significantly lower than our current estimates (Fig. 11). Never the less, given the uncertainties on jet lifetimes, such changes would not rule out that the SF activity is induced by the jet.

5.1.2 Alternative emission mechanisms

A key scenario not explored in the preceding analysis is the possibility that the jet-aligned emission is powered by AGN photoionization (e.g. Tsvetanov et al. 1996; Keel et al. 2015). Using the He ii λ1640 emission line observed in long-slit spectroscopy as a proxy for Hβ, Zirm et al. (2005) previously estimated that the total nebular continuum emission from photoionized gas could only produce a total emission of m_F850LP = 27.4 over the whole galaxy, significantly fainter than most individual components (see Figs 9, 13, and A1). Additionally, in Fig. 10, we illustrate the NIRCam colours predicted for a pure AGN photoionization model representative of emission from a potential AGN ionization cone. The AGN photoionization tracks show Ferland et al. (2013) predictions for the nebular continuum and emission line SED produced by a harder ionizing spectrum (Mathews & Ferland 1987) over a range of ionization parameters (−2.5 < log₁₀(U) < −1.5) and for the same range in metallicities as shown for the CSFH predictions. The AGN photoionization emission models are unable to reproduce the colours observed for the regions with the highest |$\text{EW}_{\text{H}\alpha +\text{N}\,\small {II}}$| , exhibiting redder F210M−F300M colours than the young stellar populations and radiative shock models.

Alternatively, the Hα nebulae associated with cooling flows in massive clusters at z < 1 have been observed to correlate spatially with radio jets within the cluster (McDonald et al. 2010; Tremblay et al. 2015). The exact nature of these Hα nebulae is however still an on-going debate, with evidence for both shock and star-formation-driven sources of nebular emission within the known population of cooling flow nebulae (McDonald, Veilleux & Rupke 2012), as well as additional ionization from X-rays (McDonald et al. 2010). This scenario offers a potential mechanism for the observed high |$\text{EW}_{\text{H}\alpha +\text{N}\,\small {II}}$| emission in TNJ1338 to be spatially associated with, but not necessarily directly triggered by, the jet activity. While it is unlikely that a cooling flow such as those seen in virialized clusters could be in place as early as z = 4.11 (Santos et al. 2008; McDonald et al. 2013), the same phenomenon may also arise in the cold accretion flows thought to fuel massive galaxy growth at high redshift (Kereš et al. 2005; Dekel et al. 2009).

Never the less, without spectroscopic observations to measure detailed line ratios and disentangle different kinematic components (e.g. narrow lines, out-flowing gas, and shock components), we cannot exclude an additional contribution from the AGN within either the host galaxy emission or the jet-aligned emission. Further JWST observations with JWST/NIRSpec IFU (GO Programme 1964) should enable conclusive and detailed measurements on the relative contributions of stellar, shock, and AGN-driven emission across TNJ1338, as well as resolved kinematic information (see for example, the complex scenario revealed by NIRSpec around a luminous z ∼ 3 red quasar; Vayner et al. 2023). Additional detailed analysis from our NIRCam photometry alone is therefore not warranted at this time.

5.2 TNJ1338 host galaxy properties

As shown in Fig. 11, the total stellar mass we infer for TNJ1338 (including all regions) is log₁₀(M/M_⊙) = 10.97 ± 0.06 assuming a Calzetti et al. (2000) dust attenuation law (or log₁₀(M/M_⊙) = 10.88 ± 0.04 assuming SMC; Gordon et al. 2003). Excluding any region with emission that could be dominated by a radiative shock only reduces the total mass by 0.02 dex. Despite the significant differences in input data and modelling assumptions, these values are in agreement with the previous mass estimates for TNJ1338 of log₁₀(M/M_⊙) ∼ 11 presented by Overzier et al. (2009). Our NIRCam observations therefore confirm that TNJ1338 is amongst the most massive galaxies at this epoch, with a mass significantly above the knee of the galaxy stellar mass function (log₁₀(M^⋆/M_⊙) ≈ 10.5; Duncan et al. 2014; Song et al. 2016). For a typical field galaxy, this stellar mass would place TNJ1338 in a dark matter halo with mass log₁₀(M_h/M_⊙) ∼ 12.5 (Legrand et al. 2019). However, with TNJ1338 also residing with a substantial galaxy overdensity (Venemans et al. 2002; Overzier et al. 2008; Saito et al. 2015), this may represent a conservative estimate. The full properties of the TNJ1338 proto-cluster revealed by NIRCam, the associated dark matter halo(s), and corresponding cosmological context will be explored in subsequent studies.

We do not observe a large variation in the overall SED across the three regions that contain the bulk of the host galaxy stellar mass (see Fig. 13). This consistency can also be seen in the non-parametric SFHs inferred for the three galaxy regions (Fig. 14, which all support a relatively flat or smoothly rising SFH until a recent burst of SF within the most recent age bin (consistent with the potentially jet-induced SF activity). However, the tight constraints on |$\text{EW}_{\text{H}\alpha +\text{N}\,\small {II}}$| provided by the NIRCam medium bands are able to reveal subtle differences in the young stellar populations present in these regions. The reddening in F300M−F335M along the southern radio jet axis illustrated in Fig. 5 can be seen in the integrated SEDs (Fig. 13) with a significant enhancement in |$\text{EW}_{\text{H}\alpha +\text{N}\,\small {II}}$| compared to the core or northern galaxy. This stronger emission line contribution results in a much greater upturn in recent SF activity for the 0 < t < 10 Myr bin, with the SFR in the south galaxy region increasing by a factor ∼5 compared to ∼2 in the core. We also note that the southern side of the galaxy core is the only region of the host galaxy that has strong evidence favouring the complex SFH regardless of dust assumption (Fig. 7).

Given the co-location of the bright radio continuum emission with the enhancement of recent SF in the host galaxy, we interpret these observations as further evidence supporting jet-induced SF as the dominant mechanism driving the extended high |$\text{EW}_{\text{H}\alpha +\text{N}\,\small {II}}$| emission. As above, we defer a more complete analysis of the formation history of TNJ1338 in the context of both its surrounding protocluster and the wider galaxy population at z ∼ 4 to subsequent studies (Duncan et al. in preparation).

6 SUMMARY AND FUTURE PROSPECTS

We have used JWST/NIRCam to study the resolved rest-frame optical properties of the z = 4.11 high-redshift radio galaxy TNJ1338. Our results are summarized as follows:

• Deep NIRCam imaging confirms the clear presence of a bright galaxy core in the South of the emission region. We observe significant variation in the F300M−F335M colour (a proxy for |$\text{EW}_{\text{H}\alpha +\text{N}\,\small {II}}$|⁠) across the associated jet-aligned emission, with enhanced emission line contribution associated with the Northern luminous radio lobe and the fainter Southern jet.

• JWST’s unique combination of resolution and sensitivity allow us to study the detailed rest-frame optical properties of the jet-aligned emission and host galaxy for the first time. Resolved SED modelling incorporating both a range of stellar population assumptions and shock models reveal that the SED of the emission closest to the jet is best-fit by radiative shock models, surrounded by regions more consistent with an extremely young burst of star formation.

• Thanks to the tight constraints on the rest-frame |$\text{EW}_{\text{H}\alpha +\text{N}\,\small {II}}$| provided by our NIRCam medium-bands, we are able to constrain the mass-weighted ages of the star formation associated with the jet activity to t_mass < 4 Myr, consistent with the expected lifetime of radio jets with the projected size observed in TNJ1338. Together, these new observations support the previously hypothesized physical picture of TNJ1338 being a massive galaxy undergoing significant obscured AGN activity that is powering the luminous radio jets, which are themselves driving radiatively driven shocks and extensive recent star formation in its wake. However, additional detailed spectroscopic observations are required to conclusively confirm this scenario.

• Under the assumption that all of the jet-aligned emission is dominated by stellar populations, the total current SFR estimated ranges from |$\log _{10}(\text{SFR}/\text{M}_{\odot }\, \text{yr}^{-1}) = 3.02^{+0.07}_{-0.06}$| to |$\log _{10}(\text{SFR}/\text{M}_{\odot }\, \text{yr}^{-1}) = 3.21^{+0.06}_{-0.04}$| depending on the assumed dust attenuation law. Even conservatively excluding any region of the jet-aligned optical emission that may be dominated by radiative shocks, the SFR attributed to the jet-triggered starburst amounts to at least |$\sim 500\, \text{M}_{\odot }\, \text{yr}^{-1}$| (∼5 × than previously estimated from rest-frame UV emission alone).

• Our NIRCam observations confirm previous estimates placing TNJ1338 as one of the most massive galaxies at this epoch, with total stellar mass estimates of log₁₀(M/M_⊙) = 10.88 ± 0.04 to log₁₀(M/M_⊙) = 10.97 ± 0.06 assuming SMC (Gordon et al. 2003) or Calzetti et al. (2000) dust attenuation, respectively.

The scale of the jet-triggered SF activity revealed by our NIRCam observations highlights the role that positive AGN feedback could play in the formation of massive galaxies. The observations presented here, and the complementary studies (Cheng et al. 2022; Windhorst et al. 2022), represent only a fraction of the full potential contained in the NIRCam observations. In the second paper in this series (Duncan et al. in preparation), we will exploit the extraordinary sensitivity of NIRCam to Hα emission at z ∼ 4.1 to produce a complete census of the proto-cluster environment surrounding TNJ1338. Additionally, Cycle 1 NIRSpec IFU observations (GO Programme 1964) will add a wealth of new information, including kinematics and resolved maps of emission line diagnostics, that will allow conclusive information on the precise physical mechanisms within the jet-aligned emission. Full spectro-photometric SED modelling (Johnson et al. 2021b) combining the NIRCam and NIRSpec observations should enable more detailed constraints on the formation history of the TNJ1338 host galaxy, alleviating existing degeneracies between the dust, stellar ages, and emission line contributions.

Finally, we highlight that these observations demonstrate the extraordinary potential for further studies of HzRGs extending to higher redshifts and into the epoch of reionization. Forthcoming massively multiplexed spectroscopic surveys such as WEAVE-LOFAR (Smith et al. 2016) will soon provide large homogeneous samples of luminous HzRG out to z > 6, that combined with the power of JWST will enable systematic and statistical studies of SMBH host galaxies and the role of AGN feedback in shaping the most massive galaxies in the Universe.

ACKNOWLEDGEMENTS

KJD thanks Philip Best, Dan Smith, Roderik Overzier, Aayush Saxena, Rebecca Larson, and Fergus Cullen for valuable discussions and input in the course of this analysis. KJD acknowledges funding from the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement number 892117 (HIZRAD) and support from the STFC through an Ernest Rutherford Fellowship (grant number ST/W003120/1). This work is based on observations made with the NASA/ESA/CSA James Webb Space Telescope. The data were obtained from the Mikulski Archive for Space Telescopes at the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5–03127 for JWST. These observations are associated with JWST programmes 1176 and 2738. RAW, SHC, and RAJ acknowledge support from NASA JWST Interdisciplinary Scientist grants NAG5-12460, NNX14AN10G, and 80NSSC18K0200 from GSFC. Work by CJC acknowledges support from the European Research Council (ERC) Advanced Investigator Grant EPOCHS (788113). BLF thanks the Berkeley Center for Theoretical Physics for their hospitality during the writing of this paper. MAM acknowledges the support of a National Research Council of Canada Plaskett Fellowship, and the Australian Research Council center of Excellence for All Sky Astrophysics in 3 Dimensions (ASTRO 3D), through project number CE17010001. CNAW acknowledges funding from the JWST/NIRCam contract NASS-0215 to the University of Arizona. TAH is supported by an appointment to the NASA Postdoctoral Program (NPP) at NASA Goddard Space Flight Center, administered by Oak Ridge Associated Universities under contract with NASA.

DATA AVAILABILITY

The raw data underlying this article will be freely available on the Mikulski Archive for Space Telescopes (https://mast.stsci.edu) following the 12-month exclusive access period. In the interim, raw, and processed imaging data as well as resolved photometry catalogues will be shared on request to the principal investigator of the PEARLS programme (R. A. Windhorst) and the corresponding author (KJD).

Footnotes

References

APPENDIX: RESOLVED SEDS FOR TNJ1338

Fig. A1 presents the observed SEDs and best-fitting stellar population models for the resolved segments not presented in Figs 9 or 13.