MNRAS 000, 1–20 (2024) Preprint 9 August 2024 Compiled using MNRAS LATEX style file v3.

JWST/NIRSpec WIDE survey: a z=4.6 low-mass star-forming galaxy

hosting a jet-driven shock with low ionisation and solar metallicity
Francesco D’Eugenio 1,2⋆ , Roberto Maiolino 1,2,3 , Vijay H. Mahatma 2 , Giovanni Mazzolari 1,4,5 ,
Stefano Carniani 6 , Anna de Graaff 7 , Michael V. Maseda 8 , Eleonora Parlanti 6 , Andrew J.
Bunker 9 , Xihan Ji 1,2 , Gareth C. Jones 9 , Raffaella Morganti 10,11 , Jan Scholtz 1,2 , Sandro
Tacchella 1,2 , Clive Tadhunter 12 , Hannah Übler 1,2 and Giacomo Venturi 6
arXiv:2408.03982v1 [astro-ph.GA] 7 Aug 2024

1 KavliInstitute for Cosmology, University of Cambridge, Madingley Road, Cambridge, CB3 0HA, United Kingdom
2 Cavendish Laboratory - Astrophysics Group, University of Cambridge, 19 JJ Thomson Avenue, Cambridge, CB3 0HE, United Kingdom
3 Department of Physics and Astronomy, University College London, Gower Street, London WC1E 6BT, UK
4 Dipartimento di Fisica e Astronomia, Università di Bologna, Via Gobetti 93/2, I-40129 Bologna, Italy
5 INAF – Osservatorio di Astrofisica e Scienza dello Spazio di Bologna, Via Gobetti 93/3, I-40129 Bologna, Italy
6 Scuola Normale Superiore, Piazza dei Cavalieri 7, I-56126 Pisa, Italy
7 Max-Planck-Institut für Astronomie, Königstuhl 17, D-69117, Heidelberg, Germany
8 Department of Astronomy, University of Wisconsin-Madison, 475 N. Charter St., Madison, WI 53706 USA
9 Department of Physics, University of Oxford, Denys Wilkinson Building, Keble Road, Oxford OX1 3RH, UK
10 ASTRON, the Netherlands Institute for Radio Astronomy, Oude Hoogeveensedijk 4, 7991 PD, Dwingeloo, The Netherlands
11 Kapteyn Astronomical Institute, University of Groningen, Postbus 800, 9700 AV Groningen, The Netherlands
12 Department of Physics and Astronomy, University of Sheffield, Sheffield, S7 3RH, UK

Accepted XXX. Received YYY; in original form ZZZ

ABSTRACT
We present NIRSpec/MSA observations from the JWST large-area survey WIDE, targeting the rest-frame UV–optical spectrum
of Ulema, a radio-AGN host at redshift z = 4.6348. The low-resolution prism spectrum displays high equivalent width nebular
emission, with remarkably high ratios of low-ionisation species of oxygen, nitrogen and sulphur, relative to hydrogen; auroral O+
emission is clearly detected, possibly also C+ . From the high-resolution grating spectrum, we measure a gas velocity dispersion
of σ ∼ 400 km s−1 , broad enough to rule out star-forming gas in equilibrium in the gravitational potential of the galaxy.
Emission-line ratio diagnostics suggest that the nebular emission is due to a shock which ran out of pre-shock gas. To infer
the physical properties of the system, we model simultaneously the galaxy spectral energy distribution (SED) and shock-driven
line emission under a Bayesian framework. We find a relatively low-mass, star-forming system (M⋆ = 1.4 × 1010 M⊙ , SFR =
70 M⊙ yr−1 ), where shock-driven emission contributes 50 per cent to the total Hβ luminosity. The nebular metallicity is near
solar – three times higher than that predicted by the mass-metallicity relation at z = 4.6, possibly related to fast-paced chemical
evolution near the galaxy nucleus. We find no evidence for a recent decline in the SFR of the galaxy, meaning that, already at
this early epoch, fast radio-mode AGN feedback was poorly coupled with the bulk of the star-forming gas; therefore, most of
the feedback energy must end up in the galaxy halo, setting the stage for future quenching.
Key words: galaxies: active – galaxies: evolution – galaxies: formation – galaxies: high-redshift – galaxies: jets

1 INTRODUCTION has accumulated both from theory (Croton et al. 2006; Bower et al.
2012; Cresci & Maiolino 2018; Harrison et al. 2018; Piotrowska
There are many physical processes that can drastically reduce (or
et al. 2022) and observations (Bluck et al. 2022; Belli et al. 2023;
‘quench’) star formation in galaxies, making them (and keeping
D’Eugenio et al. 2023; Davies et al. 2024).
them) quiescent (or ‘passive’; e.g., Man & Belli 2018). For high-
mass galaxies, quenching and quiescence are likely caused by ac- However, even if there is growing evidence for AGN-driven
tive galactic nuclei (AGN), through feedback from accreting super- quenching, there is still the open problem of how exactly AGN
massive black holes (SMBH). This hypothesis was first proposed quench star formation. Statistical studies of observations and sim-
based on simple energy arguments (Silk & Rees 1998; Haehnelt ulations show that quiescence is most closely linked to SMBH mass
et al. 1998; Binney 2004). In the meantime, supporting evidence (a tracer of the time-integrated SMBH accretion; Bluck et al. 2022;
Piotrowska et al. 2022; Brownson et al. 2022), rather than AGN lu-
minosity or Eddington ratio (Stanley et al. 2015; Scholtz et al. 2018;
⋆ E-mail: [email protected] Ward et al. 2022). Evidence of a systematic metallicity difference

© 2024 The Authors

2 F. D’Eugenio et al.
between star-forming and quiescent galaxies implies that the last pe- gramme WIDE aims to observe a large sample of bright galaxies
riod of star formation proceeds with little or no influx of pristine gas at redshifts z = 1–7; the survey goals include both statistical stud-
(Peng et al. 2015; Trussler et al. 2020). This has been interpreted as ies of the bright-end of the luminosity function, as well as the study
evidence for ‘preventive’ feedback (Bower et al. 2006, 2008), which of rare or short-lived phenomena (Maseda et al. 2024). At the same
accumulates feedback energy not in the galaxy, but in its surround- time, modelling of the spectral energy distribution (SED) of galaxies
ing halo. This form of feedback is not effective at stopping ongoing has seen tremendous progress, with the development of many flexi-
star formation, but – by heating the galaxy halo – prevents cold gas ble tools (e.g., beagle, Chevallard & Charlot 2016; bagpipes, Carnall
from accumulating and fuelling future star formation. et al. 2019; and cigale, Boquien et al. 2019; to name a few). Mod-
The alternative possibility of ‘ejective’ AGN feedback, in contrast elling the emission of shocked gas is only a matter of correctly im-
to preventive feedback, removes the fuel of star formation through plementing existing shock models (e.g., Allen et al. 2008; Alarie &
fast, multi-phase outflows (e.g., Morganti et al. 2005; Hopkins et al. Morisset 2019) into these existing SED-fitting frameworks.
2008; Nesvadba et al. 2008; Feruglio et al. 2010; Harrison et al. In this work, we study the emission-line and stellar-population
2012; Maiolino et al. 2012; Cicone et al. 2015; Cresci et al. 2015). properties of a radio-loud AGN and its host at z = 4.6.
While there is no consensus on whether or not these outflows quench Thanks to WIDE’s combination of low- and high-spectral-resolution
star formation (see, e.g., Förster Schreiber et al. 2019; Carniani et al. JWST/NIRSpec spectroscopy (§ 2), we are able to study the ionisa-
2016; Maiolino et al. 2017; Perna et al. 2020; Scholtz et al. 2020, tion, metallicity and kinematic properties of the gas in detail. We
2021; Lamperti et al. 2021), recent results from JWST observations do so using a novel implementation that marries existing galaxy
of quiescent galaxies at redshifts z = 3–5 suggest rapid star forma- spectro-photometric modelling (prospector; Johnson et al. 2021) to
tion and fast AGN feedback at early cosmic epochs (i.e., within the shock-driven nebular emission (mappings v, Sutherland & Dopita
first few hundreds Myr after the Big Bang; e.g. Carnall et al. 2023; 2017; via the Mexican Million Models Database, 3MdB, Alarie &
Glazebrook et al. 2024; de Graaff et al. 2024; Carnall et al. 2024), Morisset 2019). With this new approach, we simultaneously esti-
which may need some form of ejective feedback. A few studies re- mate the physical properties of the shock and those of the galaxy
port observing ejective AGN feedback in action at z>2 in massive, host (§ 4), finding a puzzling combination of a strong shock but a
quenching galaxies (Belli et al. 2023; D’Eugenio et al. 2023; Davies relatively low-mass host. We discuss the implications of our find-
et al. 2024). ings in the context of jet-driven quenching (§ 5), and close with a
The suite of models that describe ejective feedback are gener- succinct summary in § 6.
ally broken into two main categories, where the feedback loop is Throughout this work, we assume the Planck Collaboration et al.
driven either by intense radiation from the central AGN, known as (2020) cosmology, which gives a physical scale of 6.66 kpc arcsec−1
radiative-mode feedback, or by jets, known as radio-mode feedback. at the redshift of the source, z = 4.6348. All masses are total stel-
Jet-driven outflows have been observed in all gas phases (e.g., Mor- lar mass formed. We use everywhere a Chabrier (2003) initial mass
ganti et al. 2005; Holt et al. 2006, 2008, 2011; Santoro et al. 2020; function, integrated between 0.1 and 120 M⊙ ; all wavelengths are
Murthy et al. 2022). However, radio-mode feedback is thought to vacuum wavelengths, but optical emission lines are named after their
be mainly responsible for preventive feedback, due to its inefficient air wavelengths. All equivalent widths are in the rest frame. We
coupling with the star-forming disc of the host galaxy (e.g., Croton adopt a mass-weighted solar metallicity Z⊙ = 0.0142, unless oth-
et al. 2006).Yet, more recently, both theoretical and observational erwise specified.
works have revealed that low-power radio jets at low inclination an-
gles relative to the galaxy disc can significantly affect the gas disc
of star-forming galaxies (e.g., Mukherjee et al. 2016, 2018; Venturi
et al. 2021), possibly suppressing star formation via increased tur- 2 OBSERVATIONS AND DATA ANALYSIS
bulence (e.g., Mandal et al. 2021; Roy et al. 2024), in addition to Our target WIDE-AEGIS-2018003848 (R.A.: 215.034109, Dec:
suppression via the combination of ejective and preventive feedback +52.984404; 3D-HST v4.1 ID: 26395, Brammer et al. 2012; Skel-
(Morganti et al. 2021). Nevertheless, connecting AGN and quench- ton et al. 2014; CANDELS ID: 20415, Grogin et al. 2011; Koeke-
ing directly has proven challenging, due in part to timing arguments moer et al. 2011) is drawn from the ‘AEGIS’ observations of
(e.g., Harrison et al. 2018), and in part to the co-fuelling of both the WIDE survey (Maseda et al. 2024), which leverages existing
SMBHs and star formation when cold gas is present (Vito et al. panchromatic observations of the ‘Extended Groth Strip’ field (EGS;
2014; Ward et al. 2022). Rhodes et al. 2000; Davis et al. 2007), obtained as part of the
Near-infrared spectroscopic studies of AGN-driven radio jets at All-Wavelength Extended Groth Strip International Survey (AEGIS;
redshifts z = 1–4 show that these systems are invariably found Davis et al. 2007). The source was originally detected by Spitzer as
in gas-rich galaxies, with complex, multi-component emission, i.e. EGSIRAC J142008.19+525903.7 (Barmby et al. 2008), and later as-
arising from both the host galaxy, the AGN, and shocked material sociated with 1.4-GHz emission using the VLA (Willner et al. 2012).
(Nesvadba et al. 2017, hereafter: N17). These remarkable works It was observed by WIDE as a priority-1 target, based on being a ra-
pushed the capabilities of ground-based spectrographs to their limits, dio AGN (alongside other ∼ 20 radio AGN; Maseda et al. 2024),
ultimately hitting two instrumental limitations; a limited wavelength and we characterise it here as an ‘ultra low-ionisation emission-line
coverage, and the inability to probe the strongest emission lines be- AGN’ (hereafter: Ulema).
yond redshift z = 2–3. At the same time, the lack of suitable mod-
els to capture the multi-component nature of these complex sources
hampered a clear view of their physical properties. In this work, we
2.1 WIDE spectroscopy data
show that all these limitations can be addressed to provide new in-
sights into the nature of AGN feedback at high redshift. WIDE uses the micro-shutter assembly (MSA) on JWST/NIRSpec
With the launch of JWST, we can finally study distant radio AGN (Ferruit et al. 2022; Jakobsen et al. 2022), and observes each galaxy
by exploring their full range of rest-frame optical emission lines using a ‘slitlet’ consisting of 3 shutters. The spectrograph was con-
(e.g., Saxena et al. 2024; Roy et al. 2024). The NIRSpec GTO pro- figured with three disperser–filter combinations, the low-resolution

MNRAS 000, 1–20 (2024)

 JWST shock emission at z=4.6 3
λ (obs.) [µm]
1 2 3 4 5

3
(b)
0500 52◦5900600

(a) N

σ
0.5

Fλ [10-20 erg s-1cm-2Å-1] y [arcsec]

2018003848

1 shutter
E 0.0
z = 4.64 3
σ
−0.5
Dec. (J2000)

5σ
8 (c )

C iii]

C ii]

[O ii]
[Ne iii]

Hγ

Hβ
[O iii]
[N i]

[O i]

Hα+[N ii]
0400

-5 0 5
S/N
5σ
6
B - F606W
0300

[O ii]
G - F814W 4
R - F160W

[S ii]
VLA L band 2
0200

LOFAR HBA 1 kpc

0
215◦0200600 0400 0200 0000 1500 2500 3500 4500 5500 6500 7500 8500
R.A. (J2000) λ (rest) [Å]

Figure 1. HST/ACS and WFC3 false-colour image (panel a), showing rest-frame UV photometry of Ulema. The NIRSpec/MSA shutters are overlaid in white,
the cyan dashed lines are the 3- and 5-σ contours from VLA 1.4 GHz (rms=30 µJy; the beam size is 3.8 arcsec), and the dotted line is the 5-σ contour from
LOFAR 144 GHz (rms=60 µJy; beam size 6 arcsec). The source is not centred on the shutters; also note the bright UV region south-east of the galaxy. It
is unclear whether this emission is associated with the target or if it is from an interloper, but its absence in F435W is consistent with the Lyα drop at the
redshift of Ulema. Panel b shows the 2-d signal-to-noise map of the NIRSpec/MSA prism spectrum, highlighting the extended nature of this galaxy. The
1-d spectrum (panel c) shows several low-ionisation and neutral lines ([O ii]λλ3726,3729, [O i]λλ6300,6364, [S ii]λλ6716,6731, [O ii]λλ7319–7331, possibly
C ii]λλ2324–2329) and weak high-ionisation lines ([O iii]λλ4959,5007), revealing the nebular emission is powered by shocks (cf. Fig. 4)

N
0500 52◦5900600

whereas the grating spectra used only two nods and integration times
of 27 and 30 min each for G235H and G395H, respectively. We used
E the NRSIRS2 readout mode (Rauscher et al. 2012, 2017). Our target
belongs to the WIDE field number 18 (Programme ID 1213), and
was observed on 23rd March 2023. The observations were inspected
3

for signatures of short circuits (Rawle et al. 2022) and none were
σ

found.
Dec. (J2000)

The data reduction follows the procedures established in JADES

(JWST Advanced Deep Extragalactic Survey; Eisenstein et al. 2023;
0400

Curtis-Lake et al. 2023; Cameron et al. 2023; Carniani et al. 2023,

and Bunker et al. 2023). We used a combination of the STScI and
5σ

ESA NIRSpec pipelines (Alves de Oliveira et al. 2018; Ferruit et al.

2022), as described in Carniani et al. (2023). In addition, the WIDE
0300

data reduction adds some survey-specific optimisations, described

in Maseda et al. (2024). Crucially, the slit-loss corrections assume
B - F150W a point-source morphology (Bunker et al. 2023; Curti et al. 2023),
G - F277W therefore we rely on photometry to obtain unbiased colour estimates.
R - F356W
0200

1 kpc
We note that in the G235H grating we have only one line detec-
VLA L band tion of [O ii]λλ3726,3729, and it has low signal-to-noise ratio (SNR).
Therefore, in the remainder of this article, we use only the G395H
215◦0200600 0400 0200 0000 grating (hereafter, the grating). We do not apply a Galactic fore-
R.A. (J2000) ground extinction correction, but note that at the location of this tar-
get, the Galaxy has AV = 0.03 mag, implying a correction of 2–1
Figure 2. JWST/NIRCam false-colour image of the target, from CEERS
(Bagley et al. 2023). In the rest-frame optical, the target tentatively displays
per cent between 0.6–1 µm, and even smaller at longer wavelengths.
a smoother structure than in the rest-frame UV (cf. Fig. 1a). Because the The galaxy morphology and the prism spectrum are displayed
target centre lies outside of the NIRCam footprint, we do not use NIRCam in Fig. 1. We obtained the reduced photometric images from the
photometry in the analysis. The VLA L-band contours are the same as in version 1 of the DAWN JWST Archive (DJA; Valentino et al.
Fig. 1a. 2023). Panel a shows rest-frame UV photometry from HST/ACS
and WFC3/IR, with overlaid the position of the NIRSpec shut-
ters. Ulema presents a complex morphology, consisting of several
PRISM/CLEAR (hereafter; prism), and the high-resolution config- clumps. The MSA shutters do not cover the brightest clump, which
urations G235H/F170LP and G395H/F290LP. The prism spectrum may indicate the centre of the galaxy. However, the precise centre
covers 0.7–5.3 µm with spectral resolution R = 30–300, while the of the source is hard to determine in the rest-UV photometry. One
two gratings cover together 1.66–5.3 µm with R = 2, 700 (Jakob- of these clumps shows a distinctively blue colour. JWST/NIRCam
sen et al. 2022). All spectra were observed using nodding; the prism imaging is available through the Cosmic Evolution Early Release
observations use three nods and a total integration time of 41 min, Science (CEERS; Bagley et al. 2023), but only for part of the galaxy.

MNRAS 000, 1–20 (2024)

4 F. D’Eugenio et al.
From these images (Fig. 2), the galaxy seems significantly smoother Summary of literature photometry for Ulema.
and more extended at rest-frame optical wavelengths than in the rest- Instrument Filter Flux Source
UV image from HST. All together, these properties suggest that the [µJy]
UV image is shaped primarily by the distribution of dust and/or by HST/ACS f606w† 0.05 ± 0.01 (Koekemoer et al. 2011)
clumpy star formation. Note that the CEERS observations do not HST/ACS f814w 0.15 ± 0.02 "
fully cover the galaxy; the target centre (as estimated from HST) HST/WFC3 IR f125w 0.18 ± 0.02 "
seems to lie outside of the NIRCam footprint, which makes measur- HST/WFC3 IR f140w 0.31 ± 0.05 (Skelton et al. 2014)
ing accurate total fluxes even more challenging. For this reason, in HST/WFC3 IR f160w 0.25 ± 0.03 (Koekemoer et al. 2011)
CFHT/WIRCAM J† < 0.38 (Bielby et al. 2012)
the following we do not use NIRCam photometry, and supplement
CFHT/WIRCAM H 0.40 ± 0.11 "
our spectroscopy observations with legacy data from Stefanon et al.
CFHT/WIRCAM Ks < 0.43 "
(2017, see our Table 1). Spitzer/IRAC 3.6 µm 2.12 ± 0.11 (Ashby et al. 2015)
The 1-d spectral extraction (Fig. 1c) shows bright emission by Spitzer/IRAC 4.5 µm 1.27 ± 0.11 "
low-ionisation species like [O ii]λλ3726,3729, [O i]λλ6300,6364, Spitzer/IRAC 5.8 µm < 3.11 (Bielby et al. 2012)
[S ii]λλ6716,6731, and [O ii]λλ7319–7331. To our knowledge, this is Spitzer/IRAC 8.0 µm < 2.64 "
the most distant detection of the [O i] and auroral [O ii] lines to date.
All these lines are bright compared to [O iii]λλ4959,5007, which Table 1. All measurements are from the publicly available table of Stefanon
et al. (2017); the Source column indicates the original source of the data.
is a combination usually seen in high-metallicity, low-ionisation
We do not use existing JWST photometry because Ulema lies partially or
galaxies and in shock-dominated systems. Another line group typ-
completely outside of the NIRCam detector. † Unused bands, see main text.
ically seen in shocks is the auroral [O ii]λλ7319–7331, which is
clearly detected; these lines are coupled to tentative detections of
C ii]λλ2324–2329 and [N i]λλ5198,5200. The spectrum is clearly ex- From this value of α, we infer k-corrected radio luminosity den-
tended in the slit, and displays a clear gradient in both emission- sities of P1.4 GHz = 9.1 ± 1.8 × 1024 W Hz−1 and P144 MHz = 6.0 ±
line ratios and equivalent width (EW; e.g., at y = −0.5 arcsec 0.6 × 1025 W Hz−1 (following the approach of Pracy et al. 2016).
[O ii]λλ3726,3729 is fainter relative to the Hα +[N ii]λλ6549,6584 These values imply that Ulema is a low-luminosity radio galaxy.
complex, while EW(Hα +[N ii]λλ6549,6584) clearly increases from Nevertheless, the 1.4-GHz radio luminosity is two orders of mag-
y = 0 to −0.5 arcsec, due to the continuum becoming fainter). These nitude above the star-formation dominated threshold (Pracy et al.
strong emission-line gradients mean that applying a photometry- 2016), implying that this source is dominated by AGN emission.
based aperture correction introduces a bias against the average prop- This is confirmed by the fact that the radio-inferred SFR clearly
erties of the target. exceeds the value derived from Hα. Indeed, if the radio emission
was dominated by star formation, then from P1.4 GHz we would esti-
mate a SFR1.4 GHz of 300 M⊙ yr−1 (Davies et al. 2017), while from
2.2 Ancillary data: photometry and radio observations P144 MHz we would infer SFR144 MHz = 300–1, 000 M⊙ yr−1 (Heesen
et al. 2022, or even higher, SFR144 MHz ≈ 2, 000 M⊙ yr−1 using either
We use photometry spanning from optical to near-infrared, obtained eq. 1 or 2 from Smith et al. 2021). Both values are in excess of our
from the publicly available catalogue of Stefanon et al. (2017). Our aperture-corrected estimate of SFR10 ≈ 70 M⊙ yr−1 (Table 3), vali-
motivation for using photometry is to infer accurate aperture cor- dating our assessment that the radio emission has an AGN origin.
rections by combining the total fluxes from photometry with the de- The spectral index α = −0.83 is typically characterised as ‘steep
tailed physical information provided by the spectra. This comparison spectrum’ at these frequencies for sources with an active jet (e.g.,
is performed within the Bayesian inference framework of prospec- Zajaček et al. 2019): the high-energy components of radio galaxies
tor (§ 4.1 and Appendix A). (i.e. the jets and hotspots) are thought to have flat spectra at low ra-
Because we do not model the resonant Lyα line, we exclude the dio frequencies (αlow ≈ −0.5; associated with particle acceleration
HST/ACS F606W filter. Additionally, we do not use the J-band from strong shocks, e.g., Meisenheimer et al. 1989). However, sta-
upper limit from CFHT/WIRCAM (we have more informative de- tistical studies of radio galaxies as a population have median spec-
tections from HST/WFC3 IR in the same wavelength region). Our tral indices α1400
150 between -0.85 and -0.75 (e.g., Mahony et al. 2016;
photometry selection is summarised in Table 1, and spans observed Kutkin et al. 2023), placing Ulema close to the population median.
wavelengths between 0.8 µm (with HST/ACS) and 8 µm (with Given the uncertainties in the redshift evolution of α measured in
Spitzer/IRAC; aside from VLA and LOFAR, the longest-wavelength radio surveys (due, e.g., to surface-brightness dimming at high red-
detection is at 4.5 µm). shift, redshift-dependent inverse-Compton losses, and intrinsic evo-
This source is also detected at radio wavelengths by the VLA at lution), we explore the possibility that for Ulema we are seeing the
1.4 GHz, with a flux density S 1.4 GHz = 0.056 ± 0.011 mJy (Ivi- old, leftover plasma from previous jet activity, or that the jet is recur-
son et al. 2007; Willner et al. 2012). The source is also detected rent (consistent with studies on the radio spectra of low ionization
in the LOFAR Two-Metre Sky Survey (LoTSS DR2; Shimwell AGN with star-forming activity, e.g. Zajaček et al. 2019). At the
et al. 2022) at 144 MHz at 6-arcsec resolution, with a flux density redshift of the source (z = 4.6), we expect strong radiative losses
of S 144 MHz = 0.37 ± 0.037 mJy1 . Based on the measurements at at radio frequencies due to inverse-Compton scattering (low energy
144 MHz and 1.4 GHz, the radio spectral index (slope of the SED electrons up-scattering CMB photons), meaning a natural depletion
between 144 and 1.4 GHz, assuming a power-law energy distribu- of low energy electrons so that we ought to be sensitive only to the
tion) is −0.83 ± 0.09 (using the convention S ∝ να where α is the high energy, flat-spectrum components related to a jet. That we still
spectral index). observe a relatively steep radio spectrum (and not α ≈ −0.5) might
suggest that the radio jet is no longer active, or is restarting, and
we are seeing the leftover plasma. This would further support the
1 We have assumed a flux density error of 10 percent in line with typical conclusion drawn by the high flux ratio between low-ionisation and
errors in the survey release (Shimwell et al. 2019). high-ionisation species, which suggests an old-shock origin (Sec-

MNRAS 000, 1–20 (2024)

 JWST shock emission at z=4.6 5
tion 3 and Fig. 4). However, given that we do not have robust mea- Redshifts and nebular emission line fluxes from the WIDE NIRSpec/MSA
data of Ulema.
surements to determine an accurate radio spectral index at multi-
Redshift — 4.644 ± 0.001
ple radio frequencies, our interpretation is only tentative. Resolved
Wavelength Flux
observations at higher angular resolution are needed to determine Line(s)
[Å vacuum] [10−18 erg s−1 cm−2 ]
whether an active jet exists in Ulema. C iii]λλ1907,1909 † 1907.71 < 3.3
The target is undetected in X-rays, despite being covered by the C ii]λλ2324–2329 † 2325.40 5.4 ± 0.9
deepest Chandra X-ray observations available in the field (AEGIS- Mg iiλλ2796,2803 † 2799.94 < 2.4
XD, 800 ks; Nandra et al. 2015), implying an upper limit on the [O ii]λλ3726,3729 † 3728.49 11.0 ± 0.5
X-ray luminosity LX ≲ 5 ×1043 erg s−1 cm−2 . [Ne iii]λλ3869,3967 ‡ 3870.86
< 1.5

ppxf § 2.3
(ratio 3.32) 3968.59

PRISM
Hγ 4341.65 < 0.9
2.3 Emission-line fitting – prism spectrum Hβ 4862.64 1.9 ± 0.3
[O iii]λλ4959,5007 ‡ 4960.30
We fit the emission lines using the χ2 -minimisation algorithm ppxf 3.2 ± 0.5
(ratio 0.335) 5008.24
(Cappellari 2023), which models simultaneously the emission lines He iλ5875 5877.25 < 0.6
and the stellar continuum. The emission lines are modelled as Gaus- [O i]λλ6300,6364 ‡ 6302.05
sians, and the stellar continuum is a non-negative superposition of 7.6 ± 0.3
(ratio 3.03) 6363.67
simple stellar population (SSP) spectra spanning a grid of ages and Hα +[N ii]λλ6549,6584 † 6564.52 24.0 ± 0.4
metallicities. The SSP spectra employ MIST isochrones (Choi et al. [S ii]λλ6716,6731 † 6725.00 7.6 ± 0.3
2016) and C3K model atmospheres (Conroy et al. 2019). Both the He iλ7065 7067.14 < 0.9
emission-line and SSP spectra have been pre-convolved to match [O ii]λλ7319–7331 ‡ 7321.94
1.5 ± 0.2
the instrumental resolution of the prism spectrum (Jakobsen et al. (ratio 1.92) 7332.21
2022). Because the size of our target is larger than the width of Redshift — 4.6348 ± 0.0003
the MSA shutters, we use the nominal spectral resolution of each Velocity dispersion σ [km s−1 ] 400 ± 35
disperser (Jakobsen et al. 2022, i.e., we do not calculate an ad-hoc Wavelength Flux
line spread function, which would be appropriate for sources that ad-hoc model § 2.4 Line(s)
[Å vacuum] [10−18 erg s−1 cm−2 ]
are smaller than the shutters, de Graaff et al. 2023). The full list [O i]λλ6300,6364 ‡ 6302.05 6.6 ± 0.6
G395H

of emission lines and doublets we measure is presented in Table 2. (ratio 3.03) 6363.67 —
Among these lines, we highlight the presence of C ii]λλ2324–2329 [N ii]λλ6549,6584 ‡ 6549.86 —
(ratio 0.34) 6585.27 8.2 ± 0.6
and [O ii]λλ7319–7331– two auroral line groups characteristically
Hα 6564.52 11.6 ± 0.8
associated with strong shocks (e.g., Allen et al. 2008, their figs. 16
[S ii]λ6716 6718.30 3.7 ± 0.6
and 18). [S ii]λ6731 6732.67 4.2 ± 0.6

Table 2. Spectral measurements from the low resolution prism spectrum

2.4 Emission-line fitting – grating spectrum (top), and from the high-resolution G395H grating (bottom). Non detections
In the G395H spectrum, we consider only the brightest lines, are given as 3-σ upper limits. The flux of lines observed in both configura-
tions is statistically consistent.
[O i]λλ6300,6364, [N ii]λλ6549,6584, Hα and [S ii]λλ6716,6731. All † Sets of spectrally unresolved lines are treated as a single Gaussian at the
these lines are modelled as Gaussians, with the same redshift and ve-
reported effective wavelength.
locity dispersion. The [O i]λλ6300,6364 and [N ii]λλ6549,6584 dou- ‡ Doublets that we model as two Gaussians with fixed flux ratio; the adopted
blet ratios are fixed, whereas the [S ii]λ6717/[S ii]λ6731 ratio is con- ratio (quoted in parentheses below the doublet name) is always the ratio be-
strained between 0.44 and 1.47. We model the background as a first- tween the blue and red line.
order polynomial. In total, the model has nine free parameters; the
fluxes of [O i]λ6300, Hα, [N ii]λ6584 and [S ii]λ6733; the [S ii] dou-
blet ratio; the common redshift and velocity dispersion; and two We note that a multi-component fit – though certainly appropriate
parameters for the linear background. We adopt a Bayesian frame- physically – is not warranted by the data, due to the insufficient SNR.
work and estimate the posterior distribution of the parameters using We reach this conclusion because, when attempting to fit a multi-
a Markov-Chain Monte-Carlo (MCMC) integrator (emcee; Foreman- component model consisting of a set of narrow and broad emission
Mackey et al. 2013, and references therein). The results are reported lines, we find implausible line ratios. Two components are clearly
in Table 2. The data, fiducial model, and posterior distribution are identified, with velocity dispersion 200 ± 30 and 600 ± 100 km s−1 ,
shown in Fig. 3. and with the broad component blue-shifted by 80+40 −1
−120 km s . How-

We find different redshifts between the prism and grating (more ever, the emission-line ratios of the narrow component are clearly
than 8-σ significance; Table 2). Note that a systematic difference be- too high – especially if we interpret this component as a star-forming
tween the NIRSpec prism and gratings has been already reported in disc (we find both [O i]λλ6300,6364 and [S ii]λλ6716,6731 to be
the literature (Bunker et al. 2023). The impact of this redshift dif- brighter than Hα). Therefore, we adopt the single-component fit
ference on distant-dependent quantities like mass and SFR is neg- as fiducial; in the next sections, we present a physically motivated
ligible compared to the relevant uncertainties. We therefore adopt multi-component approach.
everywhere the redshift zspec = 4.6348 from the grating. This mea-
surement has formal uncertainties (inter-percentile range from the
marginalised posterior) of 0.0004, corresponding to 20 km s−1 or 3 SHOCK-DOMINATED NEBULAR EMISSION
half a spectral pixel, but the systematic uncertainties due to the line
3.1 Emission-line diagnostic diagrams
decomposition are probably larger. When studying simultaneously
the prism and grating spectra, we rescale the wavelength array of Combining information from the prism and G395H gratings, we can
the prism to match zspec . populate the classic emission-line diagnostic diagrams (Fig. 4; Bald-

MNRAS 000, 1–20 (2024)

 6 F. D’Eugenio et al.

20

Fλ [10−20 erg s−1 cm−2 Å−1 ]

15

10

5
[10−18 erg s−1 cm−2 ] [10−18 erg s−1 cm−2 ] [10−18 erg s−1 cm−2 ]

.5
10
F ([N ii]λ6583)

0
0
9.

−5
5
7.

3.45 3.50 3.55 3.60 3.65 3.70 3.75 3.80 3.85

λ [µm]
16 6.0

F ([O i]λ6300) = 6.58+0.62

−0.61 10
−18
erg s−1 cm−2
F ([N ii]λ6583) = 8.18+0.63
−0.63 10
−18
erg s−1 cm−2
14
F (Hα)

F (Hα) = 11.56+0.88
−0.83 10
−18
erg s−1 cm−2
12

F ([S ii]λ6731) = 4.15+0.66

−0.60 10
−18
erg s−1 cm−2
F ([S ii]λ6716)/ F ([S ii]λ6731) = 0.89+0.27
10

−0.22

σ = 403.57+34.99
−30.33 km s
−1
F ([S ii]λ6731)

6
5
4
3
/ F ([S ii]λ6731)
F ([S ii]λ6716)

25
1.
00
1.
75
0.
0
56 0.5
0
σ [km s−1 ]

0
48
0
40
0
32

10

12

14

16

6
.5

50

75

00

25

0
4.

6.

7.

9.

6.

7.

9.

32

40

48

56
10

0.

0.

1.

1.

F ([O i]λ6300) F ([N ii]λ6583) F (Hα) F ([S ii]λ6731) F ([S ii]λ6716) σ [km s−1 ]
[10−18 erg s−1 cm−2 ] [10−18 erg s−1 cm−2 ] [10−18 erg s−1 cm−2 ] [10−18 erg s−1 cm−2 ] / F ([S ii]λ6731)

Figure 3. High spectral resolution grating observations of Ulema (black lines with uncertainty in grey) and fiducial model (solid black; the dashed spectra are
models of the individual emission lines and doublets, from left to right [O i]λλ6300,6364, blue; [N ii]λλ6549,6584, orange; Hα, green; [S ii]λλ6716,6731,
red). The corner plot reports the best-fit parameters of the emission-line model we use in the emission-line diagnostic diagrams (Fig. 4). Note that the
[S ii]λλ6716,6731 ratio is constrained to physically permitted values, and is otherwise unconstrained by the G395H observations.

win et al. 1981; Veilleux & Osterbrock 1987). The [O iii]λ5007/Hβ radio-loud AGN from N17 (shown as grey errorbars). The dot-
ratio is measured from the prism spectrum, and its measurement ted and dashed lines in Fig. 5a–c separate three photoionisation
uncertainty is highlighted in red; all other ratios were measured regimes:, star-forming, AGN, and shocks/low-ionisation. Compared
from the grating (highlighted in blue). As a reference, the dia- to SDSS, Ulema occupies a sparsely populated region of the dia-
grams in Fig. 4 also illustrate two datasets; the first is the distri- grams, with a rare combination of low [O iii]λ5007/Hβ and high
bution of local galaxies and AGN from SDSS (z ≲ 0.1; Abaza- [O i]λ6300/Hα (panel c). However, compared to luminous radio
jian et al. 2009, the contours enclose the 68th and 99th percentiles AGN at z = 2–3, Ulema is also unique; its [S ii]λλ6716,6731/Hα
of the sample; dots represent individual objects outside of the 99th - and [O i]λ6300/Hα are a factor of 2–3 higher than the N17 sample,
percentile contours); the second dataset is the sample of bright,

MNRAS 000, 1–20 (2024)

 JWST shock emission at z=4.6 7
whereas its [O iii]λ5007/[O ii]λλ3726,3729 and [O iii]λ5007/Hβ are models combining emission from the shock and the precursor and
an order of magnitude lower. with photoionisation (e.g., Heckman & Best 2014).
The position of Ulema in Fig. 4 is consistent with local low- The mismatch between our measurements and the models with a
ionisation emission-line regions (LIER; Heckman 1980; Ho 2008; precursor may be evidence for a shock that has run out of pre-shock
Belfiore et al. 2016), which are thought to be powered by evolved gas. This lack of high-ionisation emission typical of a shock precur-
stellar populations or interstellar shocks (Yan & Blanton 2012; sor lends itself to the scenario of an old shock, because – almost by
Belfiore et al. 2016; Lee et al. 2024). However, LIERs driven by stel- definition – a young shock is less likely to have completely swept
lar emission and interstellar shocks have generally low equivalent- the galaxy ISM.
width Hα emission, of the order of a few Å (e.g., ‘retired galaxies’ Complementary evidence of a strong shock component comes
in Cid Fernandes et al. 2011). In contrast, Ulema has an Hα equiv- from the diagnostic diagram of Best et al. (2000), which shows
alent width EW(Hα)=400 Å (measured from the prism spectrum, C iii]λλ1907,1909/C ii]λλ2324–2329 vs [Ne iii]λ3869/[Ne v] λ3426
and correcting for [N ii]λλ6549,6584 emission), more in line with (their fig. 1). While we do not have constraints on the Ne ratio, the
other high-redshift radio AGN (N17). In the ‘WHAN’ diagnostic upper limit on the ratio C iii]λλ1907,1909/C ii]λλ2324–2329 places
diagram, showing rest-frame EW(Hα) vs [N ii]λ6584/Hα (Cid Fer- our galaxy in the lower portion of their diagram, together with other
nandes et al. 2011), Ulema occupies the top right quarter – a po- compact, radio-loud AGN (see e.g., Best et al. 2000, their fig. 1; or
sition associated with shock-dominated emission (e.g., Drevet Mu- Allen et al. 2008, their fig. 26).
lard et al. 2023), yet its high EW(Hα) makes it an outlier. How- In summary, there is evidence that a strong shock component is
ever, unlike the sources from N17 – which all have high line EWs – needed to reproduce the line ratios observed in Fig. 4 (in Appendix A
the emission-line ratios in Ulema are skewed toward low-ionisation we also illustrate that a simple star-forming model is insufficient).
species. Among nearby sources, the most likely spectral analogue is However, in the next sections we illustrate that even shocks alone
the IRAS F23389+0303 ULIRG (Ultra-Luminous Infra-Red Galaxy; cannot explain our observations.
Spence et al. 2018). This ULIRG (also known as 4C+03.60, Gower
et al. 1967; SDSS J234130.31+031726.7) is a relatively nearby ob-
ject (z = 0.14515 ± 0.00013; Spence et al. 2018), hosting a nearby 3.2 Shock interpolator model
companion (Kim et al. 2002) and neutral-gas outflows (Veilleux
et al. 2002). It has a SFR = 100 M⊙ yr−1 (Fluetsch et al. 2021, To efficiently generate shock models, we use the library of pre-
comparable to our estimate for Ulema) and is a radio AGN with computed shock models from the database 3MdB (Alarie & Moris-
relatively low luminosity P1.4 GHz = 3.5 ± 0.1 × 1025 W Hz−1 (up set 2019). These models are obtained using mappings v (Suther-
to about five times higher than Ulema; S 1.4 GHz = 862 ± 25 mJy, land & Dopita 2017), and consist of both shock only and shock
Condon et al. 1998; we used a spectral index α = −0.81, mea- plus precursor emission (Allen et al. 2008). The models are calcu-
sured using S 365 MHz = 2, 563 ± 63 mJy from Douglas et al. lated on a discrete grid spanning magnetic field strengths 0.0001 <
1996). Spectrally, F23389+0303 displays all the characteristics ob- B < 10 µG, shock velocities 100 < v < 1, 000 km s−1 ,
served in Ulema, including an extremely low [O iii]λ5007/Hβ ra- pre-shock densities 1 < n < 10, 000 cm−3 and metallicities
tio of 2.2 ± 0.2, and high emission-line ratios between the three of 0.0001 < Z < 0.040 (hereafter: BnZv space). The chemical abun-
[O i]λλ6300,6364, [N ii]λλ6549,6584 and [S ii]λλ6716,6731 and Hα. dances are from Gutkin et al. (2016, we assumed solar carbon-to-
However, the [O ii]λλ3726,3729 emission is weaker than in Ulema, oxygen abundance ratio). In particular, we use the grid Gutkin16
and the auroral ratio [O ii]λλ7319–7331/[O ii]λλ3726,3729 is con- and all the abundance sets from Gutkin16_ISM0d0001_C1d00
siderably higher. to Gutkin16_ISM0d040_C1d00. To calculate the value of the
emission-line ratios at arbitrary locations in BnZv space, we use a
From the emission-line ratios in Ulema, we can measure an ex- linear interpolator (spanning logarithmically in B, n, Z and v). Val-
citation index of 0.4 dex (Buttiglione et al. 2010). The low ra- ues outside the grid are not allowed. Note that the shock model does
dio luminosity and excitation index place Ulema in the realm of not include a nebular continuum, which we discuss in § 5.3.
‘low excitation radio galaxies’ (LERG, Laing et al. 1994). Still, this The shock interpolator in BnZv (hereafter, shinbnzv) is demon-
classification is not without problems, because local LERGs have strated in Appendix B1. In Appendix C, we use shinbnzv and a sim-
EW([O iii]) < 5 Å (and typically ∼ 1 Å; Best & Heckman 2012), ple dust-attenuation curve to reproduce the line ratios observed in
whereas for Ulema we find 77 ± 19 Å. This comparison is compli- Ulema. We conclude that shocks, alone, cannot simultaneously ex-
cated by the fact that NIRSpec MSA may be targeting an off-nuclear plain the high [O ii]λλ3726,3729/Hβ (which requires low attenua-
region in Ulema (Fig. 1), unlike the observations of Buttiglione et al. tion) and the high Hα/Hβ (which requires high attenuation). To re-
(2010) and Laing et al. (1994). produce these otherwise irreconcilable observations we need a mix-
In Fig. 4 we also overlay a set of shock models from (Morisset ture model of low- and high-attenuation spectra. Following the anal-
et al. 2015, see § 3.2), from which we selected pre-shock density ysis of Moy & Rocca-Volmerange (2002), and the observations of
nH = 1 cm−3 , solar metallicity and the abundances of Gutkin et al. Drevet Mulard et al. (2023), we assume that this mixture consists of
(2016). The two grids show shock-plus-precursor models (purple, a combination of shocks and star-formation.
highest [O iii]λ5007/Hβ grids) and shock-only models (orange, low- To generate simultaneously the shock emission and a galaxy
est [O iii]λ5007/Hβ grids), with the grid nodes spanning a range of SED, we incorporate shinbnzv in the SED modelling tool prospec-
magnetic field strengths 0.0001 < B < 10 µG and shock veloci- tor (see § 4). We use the serialisation library pickle, and the
ties 200 < v < 1, 000 km s−1 . Our measurements lie clearly outside new class ShockSpecModel, derived (literally) from prospector’s
the model predictions including a precursor; in contrast, the shock- SpecModel. The ShockSpecModel class has eight new parameters:
only models lie reasonably close to all the observed measurements. • shock_type, is either "shock" or "shock+precursor", deter-
At the same time, all other high-redshift radio AGN present high mines which interpolator is used
[O iii]λ5007/Hβ, and are thus more in line with high-excitation radio • shock_elum, the luminosity of the Hβ line, in units of L⊙ /M⋆
galaxies (HERGs). These other galaxies agree much better with the • shock_eline_sigma, the velocity dispersion of the shock emis-

MNRAS 000, 1–20 (2024)

8 F. D’Eugenio et al.

log [N ii]λ6583/Hα log [S ii]λλ6716,6731/Hα

−1.5 −1.0 −0.5 0.0 0.5 −1.5 −1.0 −0.5 0.0 0.5
1.5

101 1.0

log [O iii]λ5007/Hβ
[O iii]λ5007/Hβ

0.5

100 0.0
G395H
prism
rAGN −0.5
Shock only
Shock plus
10−1 precursor
(a ) (b) −1.0
10−1 100 10−1 100
[N ii]λ6583/Hα [S ii]λλ6716,6731/Hα
log [O i]λ6300/Hα log [O iii]λ5007/[O ii]λλ3726,3729
−2.0 −1.5 −1.0 −0.5 0.0 −1 0 1
1.5
AV = 1 mag
Ns
101 Swa
mp o
f AG (reddening)
1.0

log [O iii]λ5007/Hβ
[O iii]λ5007/Hβ

0.5

100
star formation
Duchy of

0.0

s
ER
o f LI −0.5
un
ty cks
ho

10−1
Co
He
re
b es
(c ) (d )
10−2 10−1 100 10−1 100 101
[O i]λ6300/Hα [O iii]λ5007/[O ii]λλ3726,3729
Figure 4. Ulema occupies a region of the BPT and VO diagrams (Baldwin et al. 1981; Veilleux & Osterbrock 1987) associated with strong shocks. Our
measurements are the thick red or red and blue errobars; we combine measurements from both the prism (red errorbars) and G395H disperser (blue errorbars).
In panel d, in addition to our measurement (transparent red), we also show the effect of a representative dust reddening of AV = 1 mag to guide the eye (solid
red). Crosses are luminous, radio-loud AGN at redshifts z = 1.5–4 (log P1.4 GHz [W Hz−1 ] ∼ 28; N17). The contours and dots are local measurements from SDSS,
with the contours enclosing the 67th and 99th percentiles of the data. The dashed and dotted demarcation lines separate regions dominated by star-formation
photoionisation, AGNs and LIERs (Kewley et al. 2001; Kauffmann et al. 2003; Kewley et al. 2006; Schawinski et al. 2007). The purple (orange) grid is a set of
shock-plus-precursor (shock-only) models from mappings v, assuming solar metallicity, density n = 1 cm−3 and varying the shock velocity and strength of the
magnetic field (Alarie & Morisset 2019). Our target lies closer to the shock-only models.

sion lines; this value is in general different from the dispersion of • shock_dust ∈ [0, ∞] [τV ]
emission lines due to star formation, and from the shock velocity.
Here Z⊙ = 0.01542, the value from Gutkin et al. (2016, see § 4.1
• shock_logB ∈ [−4, 1] [dex µG]
on how we enforce consistency with the prospector value of Z⊙ ).
• shock_logn ∈ [0, 4] [dex cm−3 ]
The dust attenuation parameter of the shock is parametrised by the
• shock_logZ ∈ [−2.18, 0.419] [dex Z⊙ ]
V−band optical depth, assuming the attenuation law of Cardelli
• shock_logv ∈ [2, 3] [dex km s−1 ]
et al. (1989). To take into account any wavelength-dependent mis-

MNRAS 000, 1–20 (2024)

 JWST shock emission at z=4.6 9
match between photometry and spectroscopy data, we use a calibra- 4.2 prospector results
tion polynomial; this polynomial re-scales the flux of the spectrum
The fiducial prospector model is shown in Fig. 5, with the
to match the level of the photometry, because generally this kind
marginalised posterior distribution for a selected subset of the pa-
of calibration mismatch is due to aperture losses in the spectrum.
rameters (a more extensive list of parameters and their posterior
The implications of this choice are discussed later in the article.
probabilities is listed in Table 3). The target has a relatively low
These polynomials are implemented following the standard method
value of M⋆ (log M⋆ /M⊙ = 10.16+0.12 −0.09 Table 3), smaller than the
of prospector, i.e. through the new class ShockPolySpecModel,
value reported in Stefanon et al. (2017, log M⋆ /M⊙ = 10.4 dex)
which is derived from prospector’s PolySpecModel and from our
which was based on photometry alone. The redshift difference (Ste-
own ShockSpecModel, in this order.
fanon et al. 2017 used a photometric redshift z = 4.2) makes the
discrepancy even worse (by 25 per cent). The main driver of this dis-
crepancy is information from the spectrum, not including the emis-
sion from the shock; indeed, the star-forming model without shocks
finds very similar M⋆ to the fiducial run with shocks (Appendix A).
It is only by removing spectroscopy that the prospector model is
able to reproduce the higher M⋆ values reported in the literature. The
4 A DUSTY, INTERMEDIATE-MASS GALAXY HOST fiducial M⋆ value is remarkably small, considering that the knee of
the galaxy mass function at redshifts z = 4–5 is of order 3 × 1010 M⊙
4.1 The prospector fiducial model (for both star-forming and quiescent galaxies; Weaver et al. 2023).
The SFH of this fiducial model is traced by the solid red line in
To model the spectro-photometric data of Ulema, we use v2.0 of
Fig. 5f; the shaded pink region is the uncertainty, the blue line is the
prospector (Johnson et al. 2021). In addition to the photometry and
locus of the star-forming sequence for galaxies with mass equal to
prism spectrum, we also model simultaneously the high-resolution
the stellar mass formed at each epoch (Popesso et al. 2023, their
grating spectrum; fitting of multiple spectra is a feature implemented
equation 10), and the horizontal dashed red line is the quiescent
in prospector v2.0. We use the ShockPolySpecModel class, set-
threshold defined by sSFR < 0.2/tU (zspec ), where tU (z) is the age
ting shock_type="shock". This choice of shock-only interpola-
of the Universe at redshift z (e.g., Franx et al. 2008; Pacifici et al.
tor follows from our analysis of the emission-line ratios in § 3.1,
2016).
due to the low [O iii]λ5007/Hβ ratio and to the brightness of low-
Using more bins at recent times, or fewer bins at earlier times
ionisation and neutral species. We use flat, uninformative priors on
does not affect our results significantly. For instance, adding a time
all shock parameters, except for shock_dust, which we set to 0,
bin between look-back times of 0 and 10 Myr, most quantities re-
and for shock_logZ, which is tied to the nebular metallicity of the
main within 1 σ from their value in the fiducial run; the exceptions
star-forming prospector model. The metallicity of the shock com-
are log n and log Bsh , which both increase by 1.5 σ, and µ, which
ponent is scaled by the ratio between the solar value in prospector
decreases by 2 σ from 1.41 to 1.03 (i.e., it reverts to the mean of the
(Z⊙ = 0.0142) and the solar value in the models of Gutkin et al.
truncated-Gaussian prior, see Table 3). Similarly, we find little dif-
(2016, Z⊙ = 0.01542), such that the absolute metallicity is con-
ferences when changing the prior on the SFH to be more bursty; in-
sistent. We make no attempt to match the detailed chemical abun-
creasing the dispersion of the Student’s t priors from 0.3 to 1, we find
dances.
that the galaxy had higher SFRs at recent times (100–300 Myr), and
We use a non-parametric star-formation history (SFH) with con-
an overall shorter SFH, resulting in roughly the same M⋆ . We also
tinuity priors (Leja et al. 2019a). We rely on nine time bins, with
tested a model with a delayed exponential SFH plus a single, instan-
the two most recent bins having duration 30 and 100 Myr, followed
taneous burst (dotted grey line in Fig. 5f); the model finds a much
by seven logarithmically spaced bins between 130 Myr and redshift
higher value of f (Hβ)sh (0.9), but the attenuation is higher and over-
z = 20 (no stars are formed earlier). We use a flat prior on stellar
all M⋆ , the SFR and the shock parameters are almost unchanged.
mass and metallicity (in log space), the Kriek & Conroy (2013) at-
The relatively robust SFH derives from the shape of the continuum,
tenuation law, with extra attenuation towards birth clouds, following
including the weak Balmer break and mildly rising UV continuum.
Charlot & Fall (2000). Overall, the parameters of the host galaxy fol-
The presence of a Balmer break argues against an AGN-dominated
low the setup of Tacchella et al. (2022b, hereafter: T22), including
continuum. Overall, prospector finds a SFH that is strongly rising at
coupling ongoing star formation to nebular emission. A summary
recent times, placing the galaxy on or even above the star-forming
of the model parameters and their prior probability distributions is
sequence (Fig. 5f).
presented in Table 3.
To adjust the weight of photometry and spectroscopy, we use the
noise ‘jitter’ term (on the spectrum only), which can scale the in-
4.3 Modelling pitfalls
put noise vector by a uniform factor (with flat prior between 0.5
and 2). While prospector can model outliers with a mixture model, In Fig. 5b we compare the maximum a-posteriori prospector model
the emission-line ratios of Ulema are so extreme that we risk mask- with the photometry (diamonds) and prism spectroscopy (black
ing some lines. For this reason, we mask spectral pixels identified line). We find several discrepancies, which we discuss in the fol-
as outliers in the input spectrum, and do not use any outlier model lowing.
at runtime. The posterior probability is estimated using nested sam- There is a wavelength offset of ≈ 1–2 pixels at the location of
pling (Skilling 2004, using the library dynesty; Speagle 2020), but C ii]λλ2324–2329. A strong auroral C ii]λλ2324–2329 is a signature
we checked that our results do not change using the MCMC integra- of shocks (e.g., Best et al. 2000), where this line can be compara-
tor emcee. We provide an elementary test of this combination of data ble or even brighter than Hβ (Allen et al. 2008), so its detection in
and model in Appendix D, where we use a set of mock observations Ulema would not be surprising, and would indicate an overall low
to show that the model is capable of recovering the input parameters, dust attenuation. The wavelength mismatch could indicate a problem
albeit with some bias in Bsh and nsh . in the wavelength calibration of the prism spectrum, similar to what

MNRAS 000, 1–20 (2024)

10 F. D’Eugenio et al.

Table 3. Summary of the parameters, prior probabilities and posterior probabilities of the fiducial prospector model (see also Fig. 5).

Parameter Free Description Prior Posterior

(1) (2) (3) (4) (5)
zobs Y redshift N(zspec , 0.001) 4.6317+0.0002
−0.0002
log M⋆ [M⊙ ] Y total stellar mass formed U(7, 13) 10.16+0.12
Star-forming component

−0.09
log Z[Z⊙ ] Y stellar metallicity U(−2, 0.19) −0.33+0.12
−0.13
log SFR ratios Y ratio of the log SFR between adjacent bins of the non-parametric SFH T (0, 0.3, 2) —
σ⋆ [km s−1 ] Y stellar intrinsic velocity dispersion U(0, 450) 231+40
−37
n Y power-law modifier of the Kriek & Conroy (2013) dust curve (T22, their eq. 5) U(−1.0, 0.4) 0.03+0.06
−0.07
τV Y optical depth of the diffuse dust (T22, their eq. 5) G(0.3, 1; 0, 2) 0.74+0.09
−0.09
µ Y ratio between the optical depth of the birth clouds and τV (T22, their eq. 4) U(−1.0, 0.4) 1.08+0.18
−0.18
σgas [km s−1 ] Y intrinsic velocity dispersion of the star-forming gas U(30, 300) 231+32
−28
log Zgas [Z⊙ ] Y metallicity of the star-forming gas U(−2, 0.19) 0.03+0.02
−0.02
log U Y ionisation parameter of the star-forming gas U(−4, −1) −3.39+0.19
−0.11
L(Hβ)sh [L⊙ M⋆−1 ]‡
Shock component

Y luminosity of the Hβ line in the shock component U(0, 10) —

σsh [km s−1 ] Y intrinsic velocity dispersion for the shocked gas U(100, 1, 000) +15
478−18
log Bsh [µG] Y strength of the magnetic field U(−4, 1) 0.63+0.12
−0.21
log vsh [km s−1 ] Y shock velocity U(2, 3) 2.88+0.05
−0.09
log nsh [cm−3 ] Y pre-shock density U(0, 4) 1.19+0.10
−0.11
log Zsh [Z⊙ ] N shock metallicity tied to log Zgas —
τV,sh N shock optical depth, assuming the Cardelli et al. (1989) attenuation law 0 —
jspec Y multiplicative noise inflation term for spectrum U(0.5, 2) 1.0+0.3
−0.3
Other

f (Hβ)sh N flux ratio between the shock component and the total for Hβ — 0.54+0.02
−0.03
log S FR10 [M⊙ yr−1 ] N star-formation rate averaged over the last 10 Myr — 1.86+0.06
−0.08
log S FR100 [M⊙ yr−1 ] N star-formation rate averaged over the last 100 Myr — 1.83+0.08
−0.13

(1) Parameter name and units (where applicable). (2) Only parameters marked with ‘Y’ are optimised by prospector; parameters marked with ‘N’ are either
tied to other parameters (see Column 4), or are calculated after the fit from the posterior distribution (in this case, Column 4 is empty). (3) Parameter
description. (4) Parameter prior probability distribution; N(µ, σ) is the normal distribution with mean µ and dispersion σ; U(a, b) is the uniform distribution
between a and b; T (µ, σ, ν) is the Student’s t distribution with mean µ, dispersion σ and ν degrees of freedom; G(µ, σ, a, b) is the normal distribution with
mean µ and dispersion σ, truncated between a and b. (5) Median and 16th –84th percentile range of the marginalised posterior distribution; for some nuisance
parameters we do not present the posterior statistics (e.g., log SFR ratios). ‡ Units of solar luminosity per total stellar mass formed in M⊙ .

Emission-line fluxes from the fiducial prospector model (Table 3 and the [O ii]λλ3726,3729 is barely detected beyond the central shutter
Fig. 5).
(Fig. 5b). This indicates that the shock is spatially compact, rela-
Flux (SF) Flux (Shock) Ratio tive to the size of the galaxy – in agreement with the expectation for
Line(s)
[10−18 erg s−1 cm−2 ] Shock to total low-power radio AGN.
[O ii]λ3726 6.0 ± 0.8 32.1 ± 2.2 0.84
[O ii]λ3729 8.1 ± 1.1 17.3 ± 1.4 0.68 The model is biased to higher Hβ and lower [O ii]λλ7319–7331,
Hβ 6.4 ± 0.6 7.4 ± 0.4 0.54 which is explained in part by the low weight these spectral features
[O iii]λ5007 3.0 ± 1.3 6.7 ± 1.2 0.69 carry in the posterior probability (6 and 7-σ detection, respectively).
[O i]λ6300 1.5 ± 0.3 23.2 ± 1.1 0.94 Lowering artificially the noise around Hβ–[O iii]λλ4959,5007 and
Hα 30.3 ± 2.1 22.5 ± 1.3 0.43 [O ii]λλ7319–7331 (by a factor of ten), the model predicts these lines
[N ii]λ6583 16.5 ± 1.6 16.7 ± 2.0 0.50 correctly, but is unable to match the strong [O ii]λλ3726,3729 emis-
[S ii]λ6716 10.1 ± 1.8 6.0 ± 1.0 0.37 sion. To improve the fit around Hβ and [O ii]λλ7319–7331, we re-
[S ii]λ6731 7.8 ± 1.4 8.9 ± 1.7 0.53
quire introducing two more free parameters, the metallicity of the
[O ii]λλ7319–7331 0.4 ± 0.1 5.5 ± 0.5 0.92
shocked gas (independent from the metallicity of the star forming
Table 4. The model separates emission from shocks and from star-formation gas), and the shock dust attenuation. If we add only the shock atten-
photo-ionisation. We also report the fractional contribution of shocks. The uation, we find τV,sh = 3.4 (AV,sh = 3.7 mag), a value that is high,
shock fluxes assume no dust attenuation (Table 3). All fluxes are up-scaled but comparable to radio-loud AGN at high redshift (N17). Crucially,
to match the normalisation of the photometry; for this reason, all these fluxes the model predicts correctly [O ii]λλ7319–7331, but fails to predict
are ≈ 5 times brighter than the measurements (Table 2). both Hβ and [O ii]λλ3726,3729. The resulting stellar mass is 0.2 dex
lower than in the fiducial case. If we add only the shock metallic-
ity, the model correctly predicts [O ii]λλ3726,3729, but fails to pre-
reported by other works (D’Eugenio et al. 2024). Alternatively, the dict [O ii]λλ7319–7331; the resulting posterior has 0.2-dex higher
observed spectral feature is not C ii]λλ2324–2329, but an artefact; M⋆ , and different metallicities for the star-forming and shock gas
deeper observations would certainly clarify this point. (log Zgas /Z⊙ = −0.38 and log Zsh /Z⊙ = 0.30, respectively). These
Another mismatch between data and observations is the over- two alternative models give unrealistic values of either AV,sh or Zsh .
predicted K s -band flux, with the model value greatly exceeding As a third possibility, adding both τV,sh and a free Zsh , instead, pro-
the WIRCAM upper limit (panel 5b). This mismatch is certainly vides a physically plausible solution. The fiducial model parame-
due to the strong [O ii]λλ3726,3729 in the spectrum, and its much ters remain mostly unchanged, but the galaxy dust attenuation be-
lower equivalent width in the photometry; indeed, clear indications comes stronger – particularly for the birth-cloud component, where
of this spatial dependence are visible in the 2-d spectrum, where it reaches an optical depth of ≈ 2. The shock dust attenuation, in

MNRAS 000, 1–20 (2024)

 JWST shock emission at z=4.6 11
contrast, is relatively low, at τV,sh ≈ 0.9. This dust attenuation re- component should have even higher red-to-blue ratio than the em-
quires higher intrinsic [O ii]λλ3726,3729 flux, which the model is pirical measurement. However, as we have pointed out, the fiducial
able to reproduce with shock metallicity log Zsh /Z⊙ = 0.03; at the model under-predicts [O ii]λλ7319–7331, resulting in a lower ratio
same time, the requirements on the star-forming gas are lowered to than what is measured in the prism data (0.10 vs 0.14). The inter-
sub-solar values of log Zgas /Z⊙ = −1.4. This expanded model repro- section of the solid black line with the shaded region from [S ii]
duces the data slightly better than the fiducial model, but the reduced gives us log ne [cm−3 ] = 3.2 and T e = 20, 000 K, while the lim-
χ2 is indistinguishable from the fiducial model (to within less than its of the grey shaded region give 2.7 < log ne [cm−3 ] < 3.5 and
1 per cent), so we have no reason to prefer one model over the other. 44, 000 > T e > 13, 000 K.
Based on the equilibrium values of ne and T e , we obtain a metal-
licity 12 + log(O/H) = 7.59 dex, or 10 per cent solar. The large dis-
crepancy (one dex) between this value and the value inferred from
5 DISCUSSION the prospector model with shocks underscores the limitations of the
approach proposed in this section.
5.1 Comparison with photoionisation equilibrium models
Analysing the physical properties of line-emitting plasma is often
5.2 Dust attenuation and metallicity
done under the assumption that all nebular emission arises from a re-
gion with uniform density, temperature and chemical abundance, in While our models find solar metallicity, there are several caveats to
both thermal and ionisation equilibrium. In contrast, the shock mod- this result. First, the fiducial model does not reproduce perfectly the
els we use encompass a diverse range of densities, temperatures and observed Hβ flux, which is over-estimated by approximately 30 per
ionisation conditions (Allen et al. 2008; Alarie & Morisset 2019). cent; this suggests an even higher metallicity. Second, we have as-
It is therefore interesting to assess how our results compare to the sumed a uniform gas-phase metallicity for both star-forming regions
common approach of photoionisation equilibrium. and the shocked gas. In reality, the shocked gas has a different spa-
Using the emission-line measurements of Table 2, we would reach tial distribution compared to the star-forming gas (as evidenced by
implausibly high metallicity. This is driven by the high Balmer the spatial variation in the emission-line ratios, Fig. 1b). Therefore,
decrement of ≈ 6; interpreting this ratio under the usual assump- it is possible that the shock-driven emission arises from a different
tions of Case B recombination and standard density and temperature region of the galaxy, possibly a small, highly enriched central zone,
(intrinsic Hα/Hβ = 2.86, Osterbrock & Ferland 2006), we require near to the galaxy SMBH. It is possible – though still unclear – that
a dust attenuation of AV ≈ 2.4 mag. Even assuming higher values these regions undergo fast-track chemical enrichment (e.g., Hamann
of the intrinsic ratio (Hα/Hβ = 3; e.g., Tacchella et al. 2022a), the & Ferland 1992; Baldwin et al. 2003; Matsuoka et al. 2017; Ji et al.
dust attenuation required would still be high (AV = 2.2 mag). Such 2024). The alternative interpretation, i.e. that the entire ISM had so-
high attenuation increases the metallicity estimate in two ways (we lar metallicity, runs against the expectations of the mass–metallicity
assumed the Cardelli et al. 1989 attenuation law); on one hand, high relation at z ≈ 4 (e.g., Curti et al. 2023; Nakajima et al. 2023), which
attenuation lowers the ratio between auroral and strong [O ii] from for a galaxy with the mass of Ulema predicts a metallicity of 0.25–
0.14 to 0.02, which in turn lowers the temperature to T e ≈ 8, 000 K 0.3 Z⊙ ; admittedly, current sample sizes at z > 4 are still small, so the
– a low temperature increases the recombination rate of hydrogen scatter about the mass–metallicity relation is not known precisely.
and reduces the emissivity of collisionally excited metal lines. On Our data do not prefer the dual-metallicity over the single-
the other hand, high attenuation increases the [O ii]λλ3726,3729/Hβ metallicity model, but the dual-metallicity model gives a combina-
ratio from the observed value of ≈ 6 to 13. With these values and a tion of solar metallicity for the shocked gas, and an extremely low
mean density of 1,000 cm−3 (from [S ii]), we get a O+ /H+ abundance gas metallicity of log Zgas /Z⊙ = −1.4 for the star-forming gas. This
of 0.0016. Neglecting higher ionisation states, this corresponds to value is roughly ten times lower than the prediction of the mass-
12 + log(O/H) = 9.2 dex, or [O/H]=0.51 dex with the solar metal- metallicity relation, so we can conclude that the dual-metallicity
licity from Gutkin et al. (2016). This exercise illustrates the perils model does not alleviate the tension between our galaxy and the
of assuming the same physical properties and dust attenuation for metallicity scaling relations. The SNR of our data does not allow
all emission lines. For our case, part of the Balmer emission comes to constrain the metallicity of both regions.
from higher-ionisation regions than the [O ii] and [S ii] lines, so esti- Note that the dual-metallicity model also requires significant dust
mating the dust attenuation from the Balmer decrement and applying attenuation in the shock component – although this amount is sig-
the resulting correction to all emission lines results in implausible nificantly less than what found for the star-forming gas (τV,sh = 0.9
gas properties, like super-solar metallicity at z = 4.6. vs τV = 2 mag for the star-forming gas). Dust is not expected to sur-
Limiting ourselves to the shock emission, and still under the as- vive in shocks, due to a number of physical processes (e.g., Draine
sumptions of equilibrium, we use pyneb (Luridiana et al. 2015) to & Salpeter 1979; Dwek et al. 1996; Jones et al. 1996; Dwek et al.
calculate the emissivity of S+ and O+ on the density–temperature 1996; Pineau des Forêts & Flower 1997; Perna & Lazzati 2002).
plane, and to compare the predicted emission-line ratios to the val- This is particularly relevant for fast shocks like the case in hand,
ues measured in Ulema. In Fig. 6, the grid is rendered by thin dashed based on the inferred shock velocity of over 600 km s−1 (Table 3).
and solid lines, coloured respectively according to the ratios of the A modest amount of attenuation could come from a foreground dust
[S ii] and [O ii] lines; the observed ratios are the thick lines, either screen, although we note that there are local cases of spatially re-
dashed (for [S ii]) or solid ([O ii]). The measurement uncertainties solved dust attenuation in jet-driven shocks (Drevet Mulard et al.
on the [S ii] ratio are encompassed by the grey-shaded region; as we 2023), i.e., some dust could be still embedded within or very near to
pointed out, this measurement is not the most constraining, due to the shocked medium. The presence of blue emission in HST imaging
the combination of broad line profile and low SNR (Fig. 3). The around the Lyα drop (Fig. 1a) – if it arises from the galaxy and not
thick solid black line is the [O ii] ratio for the shock component of from an interloper – suggests a complex dust structure with incom-
the best-fit prospector model; due to the negligible contribution of plete covering. Dust destruction could also play a role in explaining
the star-forming component to the [O ii]λλ7319–7331 flux, the shock the high nebular metallicity. Typical dust-to-metal mass ratios are in

MNRAS 000, 1–20 (2024)

 12 F. D’Eugenio et al.
(rest) [ m]
logM =10.16+0.12
0.09 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2 1.3 1.4
(a) photometry: spectroscopy: 0.8 (b) 1.5 G395H

F [10 19 erg s 1 cm 2 Å 1]
(d)

[N I]
C II]

H &[O III]

[O I]
C III]

H &[N II]
[S II]
[O II]
Data Data

19 erg s 1 cm 2 Å 1]
Model Model 1.0
Uncert. 0.6
logSFR100 =1.83+0.08 0.5
log SFR100[M yr 1]

0.13

0.4 0.0
0 5 0 5 2 4 6 8 0 5 0 5 0 5 0 5 0 .4 .8 .2 .6 5 0 5 0 5 .4 .6 .8 .0
0.5 0.7 1.0 1.2 1.5 2 2 2 2 0.0 0.2 0.5 0.7 2.6 2.6 2.6 2.6 2.7 0.4 0.5 0.5 0.6 0.0 0.0 0.0 0.1 0 0 1 1 0.4 0.6 0.7 0.9 1.0 1 1 1 2

3.5 3.6 3.7 3.8 3.9

(observed) [ m]
0.2

F [10
V =0.74+0.09
0.09

0.0
2.5 (c) 1 2 3 4 5 6 7 8
V

(observed) [5 m]
(e) G395H
=1.08+0.18
0.18
0.0 0
2.5 5 Residuals
Residuals 3.5 3.6 3.7 3.8 3.9
1 2 3 4 5 6 7 8
(observed) [ m]
logZgas =0.03+0.02
0.02
log Zgas[Z ]

102
(f)
f(H )sh =0.54+0.02
0.03 SFMS

SFR [M yr 1]
log nsh[cm 3] log vsh[km s 1] log Bsh[ G] log sh[km s 1] f(H )sh

101
log sh =2.68+0.01
0.02
Quiescence
100 101 102 103
logBsh =0.63+0.12
0.21
t [Myr]
SFH :
non param. (fiducial)
delayed + burst

logvsh =2.88+0.05
0.09
0 5 0 5 0 .6 .7 .8 .9

lognsh =1.19+0.10
0.11
10.05
10.20
10.35
10.50
.65
1.4
1.6
1.8
0.4.0
0.65
0.70
0.95
1.00
0.54
0.8
1.2
1.6

0.05
0.00
0.15
0

0.55
0.50
0.65
0
2.62
2.64
2.66
2.78
0.00
0.20
0.55
0.70
5
2.6
2.7
2.8
2.9

0. 0
1.075
1. 0
1.525
0
0.0

0.4

2.6

0.5
2
10

log M [M ] log SFR100[M yr 1] V log Zgas[Z ] f(H )sh log sh[km s 1] log Bsh[ G] log vsh[km s 1] log nsh[cm 3]

Figure 5. Posterior distribution for a selection of parameters of the spectro-photometric model of prospector, alongside a comparison between data and model
(panels b and d, for prism and G395H, respectively), the normalised residuals χ (panels c and e), and the reconstructed star-formation history (panel f). We find
that f (Hβ)sh , the Hβ flux fraction due to the shock, is equally split between shock emission and star-formation photoionisation, implying that a shock-only or
star-formation-only model may not capture the properties of this target. The gas metallicity is found to be approximately solar. The overall flux scaling is set by
the photometry (circles and errorbars). The presence of a Balmer/4000-Å break argues against an AGN-dominated continuum, and in favour of a stellar origin.
The SFH is dominated by a high SFR in the most recent 200 Myr, which places the galaxy on or above the star-forming sequence (dashed blue line Popesso et al.
2023) and clearly above the quiescence threshold (dashed red horizontal line). Using a non-parametric SFH (red solid line, fiducial) or a delayed-τ model plus
burst (grey dotted line) gives very similar results. We find no evidence of quenching, suggesting that radio-AGN feedback is not coupled with the star-forming
gas. See also Table 3.

MNRAS 000, 1–20 (2024)

 JWST shock emission at z=4.6 13
[S II]λ6716 [OII]λλ7319−7331 [OII]λλ7319−7331 The correlation between f (Hβ)sh and µ itself is harder to track
|
[S II]λ6731 obs [OII]λλ3726,3729 obs| |
[OII]λλ3726,3729 shock down, but a plausible explanation is that solutions with high f (Hβ)sh
and low µ cannot explain the high observed Hα flux. If this explana-

[OII]λλ7319−7331
105

[OII]λλ3726,3729
10-1
tion is correct, then galaxies with similar dust attenuation between
the shock- and star-forming emission may have strong degeneracies
10-2 between f (Hβ)sh and SFR10 and M⋆ .
A significant limitation of our approach is that we do not include
Te [K]

the shock nebular continuum. For low resolution spectroscopy, a

1.2
strong nebular continuum is significantly degenerate with stellar age,

[S II]λ6716
[S II]λ6731
1.0
104 and a higher nebular-continuum fraction in the UV can be hidden
0.8 by a strong Balmer break in stellar atmospheres. Indeed, there are
renowned examples of radio galaxies with strong nebular contin-
0.6
1.0 1.5 2.0 2.5 3.0 3.5 4.0 uum, with UV fractions as high as 50–80 per cent (e.g., 3C 368
−3
ne [cm ] Dickson et al. 1995; Stockton et al. 1996), and these have high
[O ii]λλ3726,3729/Hβ ratios, similar to Ulema (Inskip et al. 2002).
Figure 6. Predicted line ratios of [S ii] and [O ii] as a function of density
However, in these cases, the nebular continuum is readily seen in
and temperature, for an equilibrium model. The dashed (solid) coloured lines
the spectrum via the Balmer jump, which is not the case in Ulema
show the predicted ratios of [S ii] ([O ii]), according to the colourbar; the ob-
served [S ii] ratio is the thick dashed line (with the grey region encompassing (Fig. 1c; although, due to low SNR and spectral resolution, we can-
the uncertainties); the observed [O ii] ratio is the thick solid grey line, while not fully rule out a combination of shock-driven Balmer jump and
the thick solid black line is the [O ii] ratio from the prospector shock model. stellar Balmer break). Moreover, a strong nebular continuum is nor-
Emissivities were calculated using pyneb. mally associated with strong recombination lines like Hβ, while our
target has remarkably weak Hβ. Regardless, generalising the present
model to include the shock continuum would require a different ap-
the range 0.1–0.3 (De Vis et al. 2019); a significant suppression of proach than what used for shinbnzv.
this ratio, due to the shock destroying dust, could therefore transfer The stellar mass and SFH we infer are remarkably robust against
metals from the solid to the gas phase, hence increasing the gas- different modelling choices, from including/removing shock emis-
phase metallicity by up to an order of magnitude. sion, to changing the SFH priors and parametrisations (Fig. 5f). The
SFH parametrisation should carry a systematic error of 0.2–0.3 dex
(e.g., Leja et al. 2019b; Carnall et al. 2019; Lower et al. 2020), but
5.3 Host-galaxy properties
in our case the non-parametric and delayed-exponential models are
One of the goals of our simultaneous model of the host-galaxy SED very similar – with the difference being mostly at very early epochs.
and shocks-driven emission is to quantify the impact of neglecting Note that the delayed-exponential model also includes a burst, able
shock emission, and degeneracies between the model parameters de- to decouple old from recent SFH, but the optimal model has negligi-
scribing the host and the shocks. By comparing the fiducial model ble burst fraction. This too is in agreement with the non-parametric
to the no-shocks model (Figs. 5 and A1), we find that neglecting model, which finds a very similar solution when adding more time
the shock emission has a negligible impact on M⋆ and SFR. This is bins at recent look-back times. This spectacular agreement cannot
surprising, given the high equivalent width of the nebular emission be general (see again e.g., Leja et al. 2019b), and requires an ex-
lines. Inspecting the posterior probability distribution from Fig. 5a, planation. The most likely scenario is that the combination of weak
we find that most of the shock and galaxy parameters are fairly de- Balmer break and moderate rest-UV flux requires a fairly constant
coupled, while the usual degeneracies within each parameter subset SFR in the last few hundred Myr, which can be conveniently cap-
are readily recovered (e.g., M⋆ –SFR, SFR–τV ; or Bsh –nsh ). How- tured by both a delayed-exponential (with long e-folding time) and
ever, we also find significant degeneracy between f (Hβ)sh and SFR, a non-parametric SFH with continuity prior.
τV and µ (the extra attenuation toward stellar birth clouds). Using Finally, we turn our attention to the AGN. Outside the high radio
a partial-correlations approach (e.g., Bait et al. 2017; Vallat 2018; luminosity, we find no additional evidence of an AGN. HST imag-
Bluck et al. 2019), we identify the correlation between f (Hβ)sh and ing does not show a bright point source (Fig. 1a; JWST imaging is
µ as the most significant, with a Spearman rank correlation coeffi- severely affected by edge effects). The presence of a Balmer break
cient ρ = 0.7 after controlling for τV and SFR. The positive cor- argues against a dominant AGN continuum at UV–optical wave-
relations between f (Hβ)sh and SFR10 and SFR100 are less strong lengths.
(ρ = 0.39 and 0.44, respectively). Naively, one would expect an Similarly, we find no evidence for an upturn of the continuum at
anti-correlation between f (Hβ)sh and SFR – particularly with SFR10 ; red wavelengths, around 8 µm (rest-frame 1.5 µm; these wavelengths
this is because the Hβ flux is constrained by the data, so a higher can already show evidence of AGN dust emission; e.g., D’Eugenio
f (Hβ)sh value decreases the fraction of Hβ flux available for star et al. 2023, although, admittedly, the sensitivity of our Spitzer/IRAC
formation. What is happening, however, is that increasing f (Hβ)sh observations is low at 6–8 µm). Repeating the fiducial fit includ-
also increases the attenuation µ towards birth clouds; this means that ing an AGN dust torus in the model returns an AGN MIR lumi-
the observed fraction f (Hβ)sh increases, the total dust-corrected Hβ nosity consistent with 0. All together, these lines of evidence argue
flux also increases, thus increasing SFR10 . This is underscored by against the presence of a radiatively efficient AGN, in agreement
the fact that the correlations between f (Hβ)sh and SFR are signifi- with most local radio galaxies (e.g., Heckman & Best 2014), and
cantly reduced or even reversed after controlling for µ (partial corre- as expected from the low-ionisation nebular emission. However, in
lation coefficient ρ = −0.50 and 0.26 for SFR10 and SFR100 , respec- the local Universe, radiatively inefficient AGN are found mostly in
tively). Ultimately, it is this correlation between f (Hβ)sh and µ that high-mass radio-loud elliptical galaxies, while Ulema is clearly a
drives the robustness of M⋆ and SFR against including or excluding gas-rich, lower-mass system. This combination of LERG and star-
shocks. forming host however agrees with observations, which find the frac-

MNRAS 000, 1–20 (2024)

14 F. D’Eugenio et al.
tion of quiescent-host radio-AGN decreases with redshift (Maglioc- The lack of evidence for a strong AGN continuum (from X-rays,
chetti et al. 2016) and at z ≳ 1 LERG hosts are predominantly star- optical and MIR), coupled with the low-power radio luminosity
forming galaxies, in agreement with the evolving fraction of quies- (P1.4 GHz = 6.5 ± 1.3 × 1024 W Hz−1 ), suggests a stage of radiatively
cent galaxies (Kondapally et al. 2022). inefficient accretion onto the supermassive black hole. In the local
Universe, this accretion mode has been interpreted as due to ineffi-
cient Bondi-Hoyle accretion of hot gas, supported by morphological
5.4 A strong shock
evidence and low gas fractions of galaxies hosting low-power radio
Ulema is known to harbour an AGN, due to its radio luminosity ex- AGN (Hardcastle et al. 2006; Smolčić & Riechers 2011; Heckman
ceeding the limit of star formation (Pracy et al. 2016). The emission- & Best 2014). Ulema certainly struggles to fit this hot-gas picture;
line ratios can only be explained by a strong shock. The G395H its high SFR and low stellar mass imply a high gas fraction.
spectrum, alone, is insufficient to adequately measure the properties In terms of energy, it is difficult to accurately measure the total
of the host galaxy vs shock component. While we do find evidence energy of the ionised gas, because – under the assumption of an
of two components (with velocity dispersions 200 and 600 km s−1 ; old shock – an inordinate amount of gas could have already recom-
§ 2.4), these have unrealistic line ratios for the narrow component, bined, thus contributing little power to nebular emission. Using the
therefore we cannot trust the resulting kinematics either. Our joint standard approach, we can estimate the mass and kinetic energy of
modelling of the star-forming and shocked gas using both prism the ionised gas (e.g., Venturi et al. 2021) from the Hα luminosity,
and G395H data prefers a fast shock solution (vsh = 750 km s−1 ; from the gas density and from the gas velocity dispersion. We infer
Table 3 and Fig. 5a). Assuming the typical conversion between de- the Hα luminosity from the Hα flux of the shock component in the
projected velocity and line-of-sight velocity and velocity dispersion model (Table 4), and assume a density of 1, 000 cm−3 , giving a mass
that is used in outflows (i.e., vsh = |vlos | + 2 · σlos ; e.g., Rupke et al. of Mion ∼ 2 × 107 M⊙ . We further assume a velocity dispersion of
2005), our vsh is in loose agreement with the relatively high velocity 400–600 km s−1 , spanning the range between the single-component
dispersion of the gas, which we constrain independently from the fit and the broad component of the multiple-component fit (§ 2.4),
G395H spectrum (σ = 400 ± 35 km s−1 ; Table 2 and Fig. 3). which gives a kinetic energy of Ekin ∼ 3–8 × 1055 erg. If we as-
As an alternative explanation, we ask whether virial equilibrium sume that the radio emission is due to a relativistic jet, we estimate
can explain the observed line broadening. To address this question, the kinetic power of the jet to be Pjet = 4 × 1044 erg s−1 , by em-
we compare the velocity dispersion from the virial theorem to the ploying the local scaling relation between radio power and cavity
measured velocity dispersion (we use the calibration of Cappellari (jet) power observed in cluster galaxies showing cavities in the X-
et al. 2013). With the stellar mass estimate of 1010.2 M⊙ and a half- ray haloes filled by radio emission (Cavagnolo et al. 2010). This
light radius of 2 kpc (from the mass–size relation of Ward et al. estimate comes with large uncertainties, though we note that the cal-
2024), we estimate a virial dispersion σvir = 90 km s−1 – clearly in- ibrations based on synchrotron emission and X-ray cavities tend to
consistent with the observations. Overall, we can conclude that both agree (Cavagnolo et al. 2010; Willner et al. 2012; Tadhunter 2016).
the emission-line ratios and the gas kinematics argue for the pres- Due to the low spatial resolution of the radio observations (3.8 arc-
ence of a strong shock. sec; Fig. 1a; Ivison et al. 2007), we are unable to study the spatial
The weakness of [O iii]λ5007 relative to lower-ionisation oxygen structure of the warm ionised gas relative to the radio emission, un-
species tells us that this shock has run out of precursor gas, either by like in other high-redshift examples (Nesvadba et al. 2017; Saxena
running through the entire galaxy ISM (as seen in some local radio et al. 2024). At face value, the above estimates of the gas mass and
galaxies; e.g., Drevet Mulard et al. 2023), or by ‘breaking free’, i.e., kinetic energy and jet power align with estimates from local, well
crossing from the ISM to the low-density CGM. The first interpre- resolved galaxies (Venturi et al. 2021). This agreement suggests that
tation seems more problematic, vis-á-vis the model preference for the radio AGN is capable of powering the observed shock emission,
solar metallicity; in this case, in fact, we would have a highly en- though the uncertainties on our estimates are large. To estimate the
riched galaxy at z = 4.6 and at relatively low (sub-knee, e.g. Weaver total energy of the jet, we assume an age of tjet < 6.5 Myr, ob-
et al. 2023) stellar mass, possibly in tension with scaling relations, tained by dividing the LOFAR FWHM by a speed of 0.01 c, corre-
as we have seen. The other possibility – the shock breaking free sponding to the speed of the advancing radio hotspot (O’Dea et al.
from the galaxy ISM – seems more plausible; in this case, the shock 2009). This is again an order-of-magnitude estimate, due to the radio
metallicity would be tracing a relatively small central region of the emission being spatially unresolved. We lean toward an old shock
galaxy, likely close to the centre. Fast enrichment around AGN has age from the lack of nebular emission from high-ionisation species,
been suggested in the past, and there are many processes which can which is characteristic of shocks without precursor gas. With this
expedite a fast-track chemical evolution for these regions (Baldwin estimate, the ratio between the gas kinetic energy and the jet energy
et al. 2003; Cameron et al. 2023; Huang et al. 2023). is Ekin /(Pjet · tjet ) ∼ 0.001.
The large energy difference between the AGN and shock implies
a low-efficiency coupling between the AGN and the galaxy ISM – as
5.5 AGN feedback
seen in most local radio-loud AGN. Based on numerical simulations
As for the force driving the shock, the presence of a radio-loud AGN, and observations (Mukherjee et al. 2018; Venturi et al. 2021), this in-
coupled with the high shock velocity, strongly suggests an AGN ori- efficient coupling would require a near-polar jet alignment (relative
gin. The clumpy nature of the rest-UV image (Fig. 1a) may suggest to the galaxy disc). The inefficient coupling agrees with the galaxy
an ongoing merger, but UV morphology is often complex and is sen- SFH, which we see increasing in the last 100 Myr, suggesting that
sitive to the spatial distribution of even little foreground material, AGN feedback did not alter the evolutionary trajectory of the galaxy.
which masquerades the mass distribution. Rest-optical Spitzer imag- However, the large amount of jet energy freed from the black hole
ing has insufficient spatial resolution to investigate the galaxy mor- may couple with the galaxy CGM on halo-wide scales, seeding a
phology, but NIRCam imaging from CEERS tentatively suggests a quiescent future.
smooth rest-optical light distribution (Fig. 2), so a major merger is Alternatively, the SFR we infer may be connected to the AGN
disfavoured. A strong conclusion would require better image quality. by positive feedback (in agreement with many other cases, e.g.,

MNRAS 000, 1–20 (2024)

 JWST shock emission at z=4.6 15
Maiolino et al. 2017; Gallagher et al. 2019), for instance if the shock Our observations showcase the amazing capability of WIDE to re-
compressed the ISM and triggered increased star formation. This veal the detailed properties of bright, rare sources. Deeper observa-
scenario however seems unlikely, given that highly star-forming, tions may reveal additional emission lines, enabling a more detailed
dust-obscured emission is seen beyond the reach of the shock-driven characterisation of the physical properties of this system. Spatially
emission (as traced by [O ii]λλ3726,3729, Fig. 1b). resolved spectroscopy with NIRSpec/IFS would greatly improve our
In any case, radio-AGN feedback seems at best unable to interrupt understanding of the geometry of this system.
ongoing star formation, and at worst, even promoting it. This argues
against a major role of ejective feedback in Ulema. On the other
hand, radio AGN were likely already contributing to halo heating
one Gyr after the Big Bang, starting the work needed to achieve ACKNOWLEDGEMENTS
preventive feedback. The authors are grateful to Bernd Husemann for helping with the ini-
Larger samples and more complete spatial coverage of this and tial target selection of WIDE – particularly for AGN. FDE, RM, XJ
similar sources may help clarify the role of radio AGN in the for- and JS acknowledge support by the Science and Technology Facili-
mation and evolution of galaxies in the first two billion years of the ties Council (STFC), by the ERC through Advanced Grant 695671
Universe. “QUENCH”, and by the UKRI Frontier Research grant RISEand-
FALL. RM also acknowledges funding from a research professor-
ship from the Royal Society. GM acknowledges financial support
from the grant PRIN MIUR 2017PH3WAT (‘Black hole winds and
6 SUMMARY AND CONCLUSIONS
the baryon life cycle of galaxies’). SC and GV acknowledge sup-
In this work, we study the low-luminosity radio-loud AGN Ulema at port by European Union’s HE ERC Starting Grant No. 101040227 -
redshift z = 4.6, observed with JWST/NIRSpec as part of the large WINGS. MVM is supported by the National Science Foundation via
GTO programme WIDE (Maseda et al. 2024). AAG grant 2205519, the Wisconsin Alumni Research Foundation
The prism spectrum reveals a rich set of emission lines, includ- via grant MSN251397, and NASA via STScI grant JWST-GO-4426.
ing several metals; oxygen, nitrogen, sulphur and possibly carbon AJB and GCJ acknowledge funding from the “FirstGalaxies” Ad-
(Fig. 1a; Table 2). Emission-line ratio diagnostics show high val- vanced Grant from the European Research Council (ERC) under the
ues for low-ionisation and neutral species, indicating a large reser- European Union’s Horizon 2020 research and innovation program
voir of warm yet neutral gas, typical of shock-driven emission (Grant agreement No. 789056). ST acknowledges support by the
(Fig. 4). The low luminosity of higher ionisation species (most no- Royal Society Research Grant G125142. CT acknowledges support
tably [O iii]λλ4959,5007) favours an old shock. from STFC grants ST/R000964/1 and ST/V000853/1. HÜ gratefully
We study simultaneously the star-formation history and shock acknowledges support by the Isaac Newton Trust and by the Kavli
properties of Ulema under a Bayesian framework; we combine the Foundation through a Newton-Kavli Junior Fellowship. This study
grid of pre-computed shock models from the database 3MdB (Alarie makes use of data from AEGIS, a multiwavelength sky survey con-
& Morisset 2019) with the flexible galaxy inference tool prospector ducted with the Chandra, GALEX, Hubble, Keck, CFHT, MMT,
(Johnson et al. 2021), using a static linear interpolator. With this Subaru, Palomar, Spitzer, VLA, and other telescopes and supported
model, we are able to measure the properties of the AGN host (stel- in part by the NSF, NASA, and the STFC. The imaging data prod-
lar mass and star-formation history) while taking into account any ucts presented herein were retrieved from the Dawn JWST Archive
degeneracy with the shock parameters. The strongest degeneracies (DJA). DJA is an initiative of the Cosmic Dawn Center (DAWN),
are between the shock luminosity fraction and the dust attenuation which is funded by the Danish National Research Foundation under
of the star-forming gas (Fig. 5). grant DNRF140.
We infer a low stellar mass log M⋆ /M⊙ = 10.16 dex and an in- LOFAR is the Low Frequency Array designed and constructed
creasing SFH (Fig. 5; Table 3), the latter suggesting that AGN feed- by ASTRON. It has observing, data processing, and data storage
back did not reduce the host star-formation rate. These parameters facilities in several countries, which are owned by various par-
are robust against different modelling assumptions – even neglect- ties (each with their own funding sources), and which are collec-
ing shocks – supporting the conclusion that M⋆ and SFR of low- tively operated by the ILT foundation under a joint scientific policy.
luminosity radio-AGN hosts can be inferred successfully from SED The ILT resources have benefited from the following recent major
modelling (Kondapally et al. 2022). We find solar gas metallicity, funding sources: CNRS-INSU, Observatoire de Paris and Univer-
which we tentatively interpret as evidence for fast-track chemical sité d’Orléans, France; BMBF, MIWF-NRW, MPG, Germany; Sci-
enrichment in the central regions of the host galaxy, probably related ence Foundation Ireland (SFI), Department of Business, Enterprise
to the AGN. and Innovation (DBEI), Ireland; NWO, The Netherlands; The Sci-
The high shock velocity and the presence of a radio AGN suggest ence and Technology Facilities Council, UK; Ministry of Science
that the shock is AGN driven, similar to local radio-galaxies with ra- and Higher Education, Poland; The Istituto Nazionale di Astrofisica
dio jets. The increasing SFH suggests that outflows, if present, have (INAF), Italy.
low mass loading and are unable to quench the galaxy, contrary to This research made use of the Dutch national e-infrastructure with
what observed in quenching galaxies at z = 2–3 (e.g., Belli et al. support of the SURF Cooperative (e-infra 180169) and the LOFAR
2023; D’Eugenio et al. 2023; Davies et al. 2024). The low kinetic e-infra group. The Jülich LOFAR Long Term Archive and the Ger-
energy of the ionised gas (relative to the radio power), as well as man LOFAR network are both coordinated and operated by the
the increasing SFH, both suggest the coupling between the radio- Jülich Supercomputing Centre (JSC), and computing resources on
AGN and the galaxy ISM to be inefficient. This would leave a large the supercomputer JUWELS at JSC were provided by the Gauss
amount of energy to be ultimately dumped into the galaxy halo, a Centre for Supercomputing e.V. (grant CHTB00) through the John
necessary step towards achieving preventive feedback, which is sug- von Neumann Institute for Computing (NIC).
gested by both observations (Terrazas et al. 2017; Bluck et al. 2022, This research made use of the University of Hertfordshire high-
2024) and by theoretical models (e.g., Piotrowska et al. 2022). performance computing facility and the LOFAR-UK computing fa-

MNRAS 000, 1–20 (2024)

16 F. D’Eugenio et al.
cility located at the University of Hertfordshire and supported by Allen M. G., Groves B. A., Dopita M. A., Sutherland R. S., Kewley L. J.,
STFC [ST/P000096/1], and of the Italian LOFAR IT computing in- 2008, ApJS, 178, 20
frastructure supported and operated by INAF, and by the Physics Alves de Oliveira C., et al., 2018, in Observatory Operations: Strate-
Department of Turin university (under an agreement with Consorzio gies, Processes, and Systems VII. p. 107040Q (arXiv:1805.06922),
Interuniversitario per la Fisica Spaziale) at the C3S Supercomputing doi:10.1117/12.2313839
Ashby M. L. N., et al., 2015, ApJS, 218, 33
Centre, Italy.
Astropy Collaboration et al., 2013, A&A, 558, A33
This work made extensive use of the freely available Debian Bagley M. B., et al., 2023, ApJ, 946, L12
GNU/Linux operative system. We used the Python programming Bait O., Barway S., Wadadekar Y., 2017, MNRAS, 471, 2687
language (van Rossum 1995), maintained and distributed by the Baldwin J. A., Phillips M. M., Terlevich R., 1981, PASP, 93, 5
Python Software Foundation. We made direct use of Python pack- Baldwin J. A., Hamann F., Korista K. T., Ferland G. J., Dietrich M., Warner
ages astropy (Astropy Collaboration et al. 2013), corner (Foreman- C., 2003, ApJ, 583, 649
Mackey 2016), dynesty (Speagle 2020; Koposov et al. 2023), em- Barber C. B., Dobkin D. P., Huhdanpaa H., 1996, ACM Trans. Math. Softw.,
cee (Foreman-Mackey et al. 2013), jwst (Alves de Oliveira et al. 22, 469–483
2018), matplotlib (Hunter 2007), numpy (Harris et al. 2020), pin- Barmby P., Huang J. S., Ashby M. L. N., Eisenhardt P. R. M., Fazio G. G.,
gouin (Vallat 2018), ppxf (Cappellari & Emsellem 2004; Cappellari Willner S. P., Wright E. L., 2008, ApJS, 177, 431
Belfiore F., et al., 2016, MNRAS, 461, 3111
2017, 2022), prospector (Johnson et al. 2021) v2.0, pyneb (Luridi-
Belli S., et al., 2023, arXiv e-prints, p. arXiv:2308.05795
ana et al. 2015), python-fsps (Johnson et al. 2023), scipy (Jones et al.
Best P. N., Heckman T. M., 2012, MNRAS, 421, 1569
2001), and smplotlib (Li 2023). dynesty uses the nested sampling al- Best P. N., Röttgering H. J. A., Longair M. S., 2000, MNRAS, 311, 23
gorithms of Skilling (2004, 2006) and Higson et al. (2019). shinbnzv Bielby R., et al., 2012, A&A, 545, A23
uses the qhull triangulation algorithm (Barber et al. 1996). We also Binney J., 2004, MNRAS, 347, 1093
used the softwares fsps (Conroy et al. 2009; Conroy & Gunn 2010), Bluck A. F. L., et al., 2019, MNRAS, 485, 666
topcat, (Taylor 2005), fitsmap (known as ‘mapviewer’ in the DJA; Bluck A. F. L., Maiolino R., Brownson S., Conselice C. J., Ellison S. L.,
Hausen & Robertson 2022) and ds9 (Joye & Mandel 2003). The data Piotrowska J. M., Thorp M. D., 2022, A&A, 659, A160
reduction on the DJA makes use of grizli 2 . Bluck A. F. L., et al., 2024, ApJ, 961, 163
Boquien M., Burgarella D., Roehlly Y., Buat V., Ciesla L., Corre D., Inoue
A. K., Salas H., 2019, A&A, 622, A103
Bower R. G., Benson A. J., Malbon R., Helly J. C., Frenk C. S., Baugh C. M.,
DATA AVAILABILITY Cole S., Lacey C. G., 2006, MNRAS, 370, 645
Bower R. G., McCarthy I. G., Benson A. J., 2008, MNRAS, 390, 1399
The reduced and calibrated spectra are available as part of the first
Bower R. G., Benson A. J., Crain R. A., 2012, MNRAS, 422, 2816
public data release of WIDE (Maseda et al. 2024). Reduced im- Brammer G. B., et al., 2012, ApJS, 200, 13
ages presented herein were retrieved from the Dawn JWST Archive Brownson S., Bluck A. F. L., Maiolino R., Jones G. C., 2022, MNRAS, 511,
(DJA). DJA is an initiative of the Cosmic Dawn Center (DAWN), 1913
which is funded by the Danish National Research Foundation un- Bunker A. J., et al., 2023, arXiv e-prints, p. arXiv:2306.02467
der grant DNRF140. The photometric measurements used in the Buttiglione S., Capetti A., Celotti A., Axon D. J., Chiaberge M., Macchetto
prospector fits are available from Stefanon et al. (2017). F. D., Sparks W. B., 2010, A&A, 509, A6
This work is based on observations made with the Cameron A. J., Katz H., Rey M. P., Saxena A., 2023, MNRAS, 523, 3516
NASA/ESA/CSA James Webb Space Telescope. The data were ob- Cappellari M., 2017, MNRAS, 466, 798
tained from the Mikulski Archive for Space Telescopes at the Space Cappellari M., 2022, arXiv e-prints, p. arXiv:2208.14974
Cappellari M., 2023, MNRAS, 526, 3273
Telescope Science Institute, which is operated by the Association of
Cappellari M., Emsellem E., 2004, PASP, 116, 138
Universities for Research in Astronomy, Inc., under NASA contract Cappellari M., et al., 2013, MNRAS, 432, 1862
NAS 5-03127 for JWST. These observations are associated with Cardelli J. A., Clayton G. C., Mathis J. S., 1989, ApJ, 345, 245
programmes PID 1213 (NIRSpec WIDE MOS Survey – AEGIS; Carnall A. C., et al., 2019, MNRAS, 490, 417
PI N. Lützgendorf) and PID 1345 (The Cosmic Evolution Early Carnall A. C., et al., 2023, Nature, 619, 716
Release Science (CEERS) Survey; PI S. L. Finkelstein). Carnall A. C., et al., 2024, arXiv e-prints, p. arXiv:2405.02242
The software developed for this work will be made avail- Carniani S., et al., 2016, A&A, 591, A28
able after review. It includes the prospector derived classes Carniani S., et al., 2023, arXiv e-prints, p. arXiv:2306.11801
ShockSpecModel and ShockPolySpecModel, and a static imple- Cavagnolo K. W., McNamara B. R., Nulsen P. E. J., Carilli C. L., Jones C.,
mentation of shinbnzv. At the time of this writing, the documen- Bîrzan L., 2010, ApJ, 720, 1066
Chabrier G., 2003, PASP, 115, 763
tation is wanting, but we encourage the reader to address ques-
Charlot S., Fall S. M., 2000, ApJ, 539, 718
tions to the package maintainers. The package shows how to build
Chevallard J., Charlot S., 2016, MNRAS, 462, 1415
a custom interpolator model to use different input grids and differ- Choi J., Dotter A., Conroy C., Cantiello M., Paxton B., Johnson B. D., 2016,
ent emission-line selections; this requires access the library of pre- ApJ, 823, 102
computed shock models from the database 3MdB (Alarie & Moris- Cicone C., et al., 2015, A&A, 574, A14
set 2019), which can be requested to the database maintainers. Cid Fernandes R., Stasińska G., Mateus A., Vale Asari N., 2011, MNRAS,
413, 1687
Condon J. J., Cotton W. D., Greisen E. W., Yin Q. F., Perley R. A., Taylor
G. B., Broderick J. J., 1998, AJ, 115, 1693
REFERENCES
Conroy C., Gunn J. E., 2010, FSPS: Flexible Stellar Population Syn-
Abazajian K. N., et al., 2009, ApJS, 182, 543 thesis, Astrophysics Source Code Library, record ascl:1010.043
Alarie A., Morisset C., 2019, Rev. Mex. Astron. Astrofis., 55, 377 (ascl:1010.043)
Conroy C., Gunn J. E., White M., 2009, ApJ, 699, 486
Conroy C., Naidu R. P., Zaritsky D., Bonaca A., Cargile P., Johnson B. D.,
2 10.5281/zenodo.1146904 Caldwell N., 2019, ApJ, 887, 237

MNRAS 000, 1–20 (2024)

 JWST shock emission at z=4.6 17
Cresci G., Maiolino R., 2018, Nature Astronomy, 2, 179 Jones E., Oliphant T., Peterson P., et al., 2001, SciPy: Open source scientific
Cresci G., et al., 2015, ApJ, 799, 82 tools for Python, http://www.scipy.org/
Croton D. J., et al., 2006, MNRAS, 365, 11 Joye W. A., Mandel E., 2003, in Payne H. E., Jedrzejewski R. I., Hook R. N.,
Curti M., et al., 2023, arXiv e-prints, p. arXiv:2304.08516 eds, Astronomical Society of the Pacific Conference Series Vol. 295,
Curtis-Lake E., et al., 2023, Nature Astronomy, 7, 622 Astronomical Data Analysis Software and Systems XII. p. 489
D’Eugenio F., et al., 2023, arXiv e-prints, p. arXiv:2308.06317 Kauffmann G., et al., 2003, MNRAS, 346, 1055
D’Eugenio F., et al., 2024, arXiv e-prints, p. arXiv:2404.06531 Kewley L. J., Dopita M. A., Sutherland R. S., Heisler C. A., Trevena J., 2001,
Davies L. J. M., et al., 2017, MNRAS, 466, 2312 ApJ, 556, 121
Davies R. L., et al., 2024, MNRAS, 528, 4976 Kewley L. J., Groves B., Kauffmann G., Heckman T., 2006, MNRAS, 372,
Davis M., et al., 2007, ApJ, 660, L1 961
De Vis P., et al., 2019, A&A, 623, A5 Kim D. C., Veilleux S., Sanders D. B., 2002, ApJS, 143, 277
Dickson R., Tadhunter C., Shaw M., Clark N., Morganti R., 1995, MNRAS, Koekemoer A. M., et al., 2011, ApJS, 197, 36
273, L29 Kondapally R., et al., 2022, MNRAS, 513, 3742
Douglas J. N., Bash F. N., Bozyan F. A., Torrence G. W., Wolfe C., 1996, AJ, Koposov S., et al., 2023, joshspeagle/dynesty: v2.1.3,
111, 1945 doi:10.5281/zenodo.8408702
Draine B. T., Salpeter E. E., 1979, ApJ, 231, 77 Kriek M., Conroy C., 2013, ApJ, 775, L16
Drevet Mulard M., et al., 2023, A&A, 676, A35 Kutkin A. M., et al., 2023, A&A, 676, A37
Dwek E., Foster S. M., Vancura O., 1996, ApJ, 457, 244 Laing R. A., Jenkins C. R., Wall J. V., Unger S. W., 1994, in Bicknell G. V.,
Eisenstein D. J., et al., 2023, arXiv e-prints, p. arXiv:2306.02465 Dopita M. A., Quinn P. J., eds, Astronomical Society of the Pacific Con-
Ferruit P., et al., 2022, A&A, 661, A81 ference Series Vol. 54, The Physics of Active Galaxies. p. 201
Feruglio C., Maiolino R., Piconcelli E., Menci N., Aussel H., Lamastra A., Lamperti I., et al., 2021, A&A, 654, A90
Fiore F., 2010, A&A, 518, L155 Lee M. Y. L., Yan R., Ji X., Lin Z., 2024, in American Astronomical Society
Fluetsch A., et al., 2021, MNRAS, 505, 5753 Meeting Abstracts. p. 402.20
Foreman-Mackey D., 2016, The Journal of Open Source Software, 1, 24 Leja J., Carnall A. C., Johnson B. D., Conroy C., Speagle J. S., 2019a, ApJ,
876, 3
Foreman-Mackey D., Hogg D. W., Lang D., Goodman J., 2013, PASP, 125,
Leja J., Carnall A. C., Johnson B. D., Conroy C., Speagle J. S., 2019b, ApJ,
306
876, 3
Förster Schreiber N. M., et al., 2019, ApJ, 875, 21
Li J., 2023, AstroJacobLi/smplotlib: v0.0.8, doi:10.5281/zenodo.7966831,
Franx M., van Dokkum P. G., Förster Schreiber N. M., Wuyts S., Labbé I.,
https://doi.org/10.5281/zenodo.7966831
Toft S., 2008, ApJ, 688, 770
Lower S., Narayanan D., Leja J., Johnson B. D., Conroy C., Davé R., 2020,
Gallagher R., Maiolino R., Belfiore F., Drory N., Riffel R., Riffel R. A., 2019,
arXiv e-prints, p. arXiv:2006.03599
MNRAS, 485, 3409
Luridiana V., Morisset C., Shaw R. A., 2015, A&A, 573, A42
Glazebrook K., et al., 2024, Nature, 628, 277
Magliocchetti M., Lutz D., Santini P., Salvato M., Popesso P., Berta S., Pozzi
Gower J. F. R., Scott P. F., Wills D., 1967, Mem. RAS, 71, 49
F., 2016, MNRAS, 456, 431
Grogin N. A., et al., 2011, ApJS, 197, 35
Mahony E. K., et al., 2016, MNRAS, 463, 2997
Gutkin J., Charlot S., Bruzual G., 2016, MNRAS, 462, 1757
Maiolino R., et al., 2012, MNRAS, 425, L66
Haehnelt M. G., Natarajan P., Rees M. J., 1998, MNRAS, 300, 817
Maiolino R., et al., 2017, Nature, 544, 202
Hamann F., Ferland G., 1992, ApJ, 391, L53
Man A., Belli S., 2018, Nature Astronomy, 2, 695
Hardcastle M. J., Evans D. A., Croston J. H., 2006, MNRAS, 370, 1893
Mandal A., Mukherjee D., Federrath C., Nesvadba N. P. H., Bicknell G. V.,
Harris C. R., et al., 2020, Nature, 585, 357 Wagner A. Y., Meenakshi M., 2021, MNRAS, 508, 4738
Harrison C. M., et al., 2012, ApJ, 760, L15 Maseda M. V., et al., 2024, arXiv e-prints, p. arXiv:2403.05506
Harrison C. M., Costa T., Tadhunter C. N., Flütsch A., Kakkad D., Perna M., Matsuoka K., Nagao T., Maiolino R., Marconi A., Park D., Taniguchi Y.,
Vietri G., 2018, Nature Astronomy, 2, 198 2017, A&A, 608, A90
Hausen R., Robertson B. E., 2022, Astronomy and Computing, 39, 100586 Meisenheimer K., Roser H. J., Hiltner P. R., Yates M. G., Longair M. S.,
Heckman T. M., 1980, A&A, 87, 152 Chini R., Perley R. A., 1989, A&A, 219, 63
Heckman T. M., Best P. N., 2014, ARA&A, 52, 589 Morganti R., Tadhunter C. N., Oosterloo T. A., 2005, A&A, 444, L9
Heesen V., et al., 2022, A&A, 664, A83 Morganti R., Oosterloo T., Tadhunter C., Bernhard E. P., Raymond Oonk
Higson E., Handley W., Hobson M., Lasenby A., 2019, Statistics and Com- J. B., 2021, A&A, 656, A55
puting, 29, 891 Morisset C., Delgado-Inglada G., Flores-Fajardo N., 2015, Rev. Mex. Astron.
Ho L. C., 2008, ARA&A, 46, 475 Astrofis., 51, 103
Holt J., Tadhunter C., Morganti R., Bellamy M., González Delgado R. M., Moy E., Rocca-Volmerange B., 2002, A&A, 383, 46
Tzioumis A., Inskip K. J., 2006, MNRAS, 370, 1633 Mukherjee D., Bicknell G. V., Sutherland R., Wagner A., 2016, MNRAS,
Holt J., Tadhunter C. N., Morganti R., 2008, MNRAS, 387, 639 461, 967
Holt J., Tadhunter C. N., Morganti R., Emonts B. H. C., 2011, MNRAS, 410, Mukherjee D., Bicknell G. V., Wagner A. Y., Sutherland R. S., Silk J., 2018,
1527 MNRAS, 479, 5544
Hopkins P. F., Hernquist L., Cox T. J., Kereš D., 2008, ApJS, 175, 356 Murthy S., Morganti R., Wagner A. Y., Oosterloo T., Guillard P., Mukherjee
Huang J., Lin D. N. C., Shields G., 2023, MNRAS, 525, 5702 D., Bicknell G., 2022, Nature Astronomy, 6, 488
Hunter J. D., 2007, Computing in Science and Engineering, 9, 90 Nakajima K., Ouchi M., Isobe Y., Harikane Y., Zhang Y., Ono Y., Umeda H.,
Inskip K. J., Best P. N., Röttgering H. J. A., Rawlings S., Cotter G., Longair Oguri M., 2023, arXiv e-prints, p. arXiv:2301.12825
M. S., 2002, MNRAS, 337, 1407 Nandra K., et al., 2015, ApJS, 220, 10
Ivison R. J., et al., 2007, ApJ, 660, L77 Nesvadba N. P. H., Lehnert M. D., De Breuck C., Gilbert A. M., van Breugel
Jakobsen P., et al., 2022, A&A, 661, A80 W., 2008, A&A, 491, 407
Ji X., et al., 2024, arXiv e-prints, p. arXiv:2404.04148 Nesvadba N. P. H., De Breuck C., Lehnert M. D., Best P. N., Collet C., 2017,
Johnson B. D., Leja J., Conroy C., Speagle J. S., 2021, ApJS, 254, 22 A&A, 599, A123
Johnson B., et al., 2023, dfm/python-fsps: v0.4.4, O’Dea C. P., Daly R. A., Kharb P., Freeman K. A., Baum S. A., 2009, A&A,
doi:10.5281/zenodo.8230430, https://doi.org/10.5281/zenodo. 494, 471
8230430 Osterbrock D. E., Ferland G. J., 2006, Astrophysics of gaseous nebulae and
Jones A. P., Tielens A. G. G. M., Hollenbach D. J., 1996, ApJ, 469, 740 active galactic nuclei. University Science Books, Sausalito, California

MNRAS 000, 1–20 (2024)

18 F. D’Eugenio et al.
Pacifici C., et al., 2016, ApJ, 832, 79 (rest) [ m]
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2 1.3 1.4
Peng Y., Maiolino R., Cochrane R., 2015, Nature, 521, 192
0.8 (a) 1.5 G395H

F [10 19 erg s 1 cm 2 Å 1]
(c)

[N I]
C II]

H &[O III]

[O I]
C III]

H &[N II]
[S II]
[O II]
19 erg s 1 cm 2 Å 1]
Perna R., Lazzati D., 2002, ApJ, 580, 261
Perna M., et al., 2020, A&A, 643, A139 1.0
Pineau des Forêts G., Flower D., 1997, IAU Symposium, 178, 113 0.6
0.5
Piotrowska J. M., Bluck A. F. L., Maiolino R., Peng Y., 2022, MNRAS, 512,
0.4 0.0
1052 3.5 3.6 3.7 3.8 3.9
Planck Collaboration et al., 2020, A&A, 641, A6 (observed) [ m]
Popesso P., et al., 2023, MNRAS, 519, 1526 0.2

F [10
Pracy M. B., et al., 2016, MNRAS, 460, 2
Rauscher B. J., et al., 2012, in Holland A. D., Beletic J. W., eds, Society of 0.0
1 2 3 4 5 6 7 8
Photo-Optical Instrumentation Engineers (SPIE) Conference Series Vol. 5 (b) (observed) [5 m]
(d) G395H
8453, High Energy, Optical, and Infrared Detectors for Astronomy V. p. 0 0
84531F, doi:10.1117/12.926089 5 5
Rauscher B. J., et al., 2017, PASP, 129, 105003 3.5 3.6 3.7 3.8 3.9
1 2 3 4 5 6 7 8
Rawle T. D., et al., 2022, in Coyle L. E., Matsuura S., Perrin M. D., eds, (observed) [ m]
Society of Photo-Optical Instrumentation Engineers (SPIE) Conference
Series Vol. 12180, Space Telescopes and Instrumentation 2022: Opti-
cal, Infrared, and Millimeter Wave. p. 121803R (arXiv:2208.04673), Figure A1. Modelling the data without including shocks we find very
doi:10.1117/12.2629231 similar M⋆ as the fiducial case, but this is by happenstance; the model
Rhodes J., Refregier A., Groth E. J., 2000, ApJ, 536, 79 clearly mismatches some of the observed emission-line ratios, particularly
Roy N., et al., 2024, arXiv e-prints, p. arXiv:2401.11612 [O ii]λλ3726,3729, and Hα +[N ii]λλ6549,6584, resulting in an overesti-
Rupke D. S., Veilleux S., Sanders D. B., 2005, ApJS, 160, 115 mated stellar continuum. The mismatch is partly due to the competing re-
Santoro F., Tadhunter C., Baron D., Morganti R., Holt J., 2020, A&A, 644, quirements to have high dust attenuation (to reproduce the observed Balmer
A54 decrement) and high [O ii]λλ3726,3729 flux, which requires low attenuation,
Saxena A., et al., 2024, arXiv e-prints, p. arXiv:2401.12199 for reasonable metallicity. [O i]λλ6300,6364 is also not reproduced; its ob-
Schawinski K., Thomas D., Sarzi M., Maraston C., Kaviraj S., Joo S.-J., Yi served ratio with Hα is beyond the capabilities of star-forming models. The
S. K., Silk J., 2007, MNRAS, 382, 1415 lines and symbols are the same as Fig. 5b–e.
Scholtz J., et al., 2018, MNRAS, 475, 1288
Scholtz J., et al., 2020, MNRAS, 492, 3194
Scholtz J., et al., 2021, MNRAS, 505, 5469 Weaver J. R., et al., 2023, A&A, 677, A184
Shimwell T. W., et al., 2019, A&A, 622, A1 Willner S. P., et al., 2012, ApJ, 756, 72
Shimwell T. W., et al., 2022, A&A, 659, A1 Yan R., Blanton M. R., 2012, ApJ, 747, 61
Silk J., Rees M. J., 1998, A&A, 331, L1 Zajaček M., et al., 2019, A&A, 630, A83
Skelton R. E., et al., 2014, ApJS, 214, 24 de Graaff A., et al., 2023, arXiv e-prints, p. arXiv:2308.09742
Skilling J., 2004, in Fischer R., Preuss R., Toussaint U. V., eds, American de Graaff A., et al., 2024, arXiv e-prints, p. arXiv:2404.05683
Institute of Physics Conference Series Vol. 735, Bayesian Inference and van Rossum G., 1995, CWI Technical Report, CS-R9526
Maximum Entropy Methods in Science and Engineering: 24th Interna-
tional Workshop on Bayesian Inference and Maximum Entropy Methods
in Science and Engineering. pp 395–405, doi:10.1063/1.1835238
Skilling J., 2006, Bayesian Anal., 1, 833 APPENDIX A: SED MODELLING WITHOUT SHOCKS
Smith D. J. B., et al., 2021, A&A, 648, A6
Smolčić V., Riechers D. A., 2011, ApJ, 730, 64 As a reference, we measure the physical properties of Ulema without
Speagle J. S., 2020, MNRAS, 493, 3132 modelling the shock-driven nebular emission. We model simultane-
Spence R. A. W., Tadhunter C. N., Rose M., Rodríguez Zaurín J., 2018, MN- ously the photometry and prism spectrum using the standard version
RAS, 478, 2438 of prospector (Johnson et al. 2021). As we have seen (§ 3–4), half
Stanley F., Harrison C. M., Alexander D. M., Swinbank A. M., Aird J. A., of the nebular emission in Ulema is due to shocks.
Del Moro A., Hickox R. C., Mullaney J. R., 2015, MNRAS, 453, 591 We use the prospector PolySpecModel, with the same model
Stefanon M., et al., 2017, ApJS, 229, 32 parameters and prior probability distributions of Table 3, except
Stockton A., Ridgway S. E., Kellogg M., 1996, AJ, 112, 902
for the shock-specific parameters, which are absent. The resulting
Sutherland R. S., Dopita M. A., 2017, ApJS, 229, 34
Tacchella S., et al., 2022a, MNRAS, 513, 2904
model is shown in Fig. A1; as expected, the model fails to repro-
Tacchella S., et al., 2022b, ApJ, 926, 134 duce the flux of high-temperature and/or low-ionisation emission,
Tadhunter C., 2016, A&ARv, 24, 10 namely C ii]λλ2324–2329, [O ii]λλ3726,3729, [O i]λλ6300,6364 and
Taylor M. B., 2005, in Shopbell P., Britton M., Ebert R., eds, Astronomical [O ii]λλ7319–7331. The observed line strengths are so far from what
Society of the Pacific Conference Series Vol. 347, Astronomical Data a star-forming model can generate, that the optimisation resorts
Analysis Software and Systems XIV. p. 29 to increasing the level of the continuum to reduce the residuals
Terrazas B. A., Bell E. F., Woo J., Henriques B. M. B., 2017, ApJ, 844, 170 (Fig. A1b, around wavelengths 1.5 and 3.2–3.3 µm). In addition, the
Trussler J., Maiolino R., Maraston C., Peng Y., Thomas D., Goddard D., Lian model cannot reproduce the observed Balmer decrement, tending to
J., 2020, MNRAS, 491, 5406 both underestimate Hα and overestimate Hβ. This is due to the fact
Valentino F., et al., 2023, ApJ, 947, 20
that higher dust attenuation would make it impossible to reproduce
Vallat R., 2018, Journal of Open Source Software, 3, 1026
the already under-estimated flux of [O ii]λλ3726,3729.
Veilleux S., Osterbrock D. E., 1987, ApJS, 63, 295
Veilleux S., Kim D. C., Sanders D. B., 2002, ApJS, 143, 315 Despite this shortcoming, the resulting model has log M⋆ [M⊙ ] =
Venturi G., et al., 2021, A&A, 648, A17 10.1 ± 0.1 dex and log SFR100 [M⊙ yr−1 ] = 1.6 ± 0.1 dex, statis-
Vito F., et al., 2014, MNRAS, 441, 1059 tically consistent with the fiducial value. The agreement in SFR is
Ward S. R., Harrison C. M., Costa T., Mainieri V., 2022, MNRAS, 514, 2936 remarkable, given that in this model the emission lines are entirely
Ward E., et al., 2024, ApJ, 962, 176 due to star formation, whereas in the fiducial model 50 per cent of

MNRAS 000, 1–20 (2024)

 JWST shock emission at z=4.6 19
the Hβ luminosity is due to shocks (Table 3). However, for the star- log [N ii]λ6583/Hα
−1.0 −0.8 −0.6 −0.4 −0.2 0.0
forming model, the higher Hα and Hβ luminosities are compensated
by lower dust attenuation than in the fiducial model with shocks Original grid : Interpolated : 1.2
(τV = 0.58 ± 0.06 and µ = 0.5 ± 0.2), hence the total SFR is still B = const. B = const.
very close to the fiducial value. Similarly, the agreement in M⋆ is 101 vsh = const. vsh = const. 1.0
due to the combination of lower dust attenuation (which decreases
M⋆ ) and a systematically over-predicted continuum level (which in- 0.8

log [O iii]λ5007/Hβ
creases M⋆ ).

[O iii]λ5007/Hβ
Compared to previous studies, we find a lower M⋆ ; we trace this 0.6
difference to the inclusion of the spectrum in our analysis. Indeed,
using photometry only to constrain the above model, we too find a 0.4
higher value, log M⋆ [M⊙ ] = 10.4 ± 0.1, the same as the literature
value of 10.4 (Stefanon et al. 2017). 0.2

100 0.0
APPENDIX B: SHOCK INTERPOLATOR
−0.2
The interpolator we have built for this study predicts the flux of var- n = 1 cm−3 Z = 0.01542
ious emission lines. Specifically, we retrieve the line fluxes from the Gutkin et al. (2016) abundances
−0.4
3MdB database, then construct a grid in BnZv space. We transform 10−1 100
B, n and Z from linear to log-space (§ 3.2). To interpolate, we use [N ii]λ6583/Hα
the qhull algorithm (Barber et al. 1996) to triangulate the input grid,
Figure B1. Our linear interpolator shinbnzv correctly predicts the BPT line
then perform linear interpolation between the barycentre of each tri-
ratios from the input grid, down to the machine precision. We show the grid
angle. Gutkin16 from the database 3MdB (Alarie & Morisset 2019). We selected
The specific interpolator realisation we used in this work pre- the model with Z = 0.01542 and pre-shock density nsh = 1 cm−3 , with abun-
dicts the logarithm of the flux of 52 emission lines. These in- dances from Gutkin et al. (2016). The thin dashed lines are the original grid
clude the brightest recombination lines in the rest UV–NIR range from 3MdB, the thick solid lines are the grid interpolated with shinbnzv. Mis-
(i.e., hydrogen and helium), strong, non-resonant metal lines, matches between the two grids are due to the choice of interpolating linearly,
and the brightest metal auroral lines. Notable lines that we and are caused by the different rate of change of the numerator and denom-
do not model are Lyα, C ivλλ1548,1551 and Mg iiλλ2796,2803. inator. The grids show eight logarithmically spaced values of magnetic field
Aside from these three resonant transitions, all strong lines from strength B and shock velocity vsh . The top and bottom grids show respec-
He iiλ1640 to Paα are modelled. In addition, we model the low- tively the model including a precursor and the shock-only model.
ionisation [N i]λλ5198,5200 and the temperature-sensitive auro-
ral lines C ii]λλ2324–2329, [S ii]λ4070, [O iii]λ4363, [O i]λ5577,
[N ii]λ5755, and [O ii]λλ7319–7331. The interpolator itself does not et al. (1989) attenuation curve, parametrised by the V-band attenua-
apply any dust reddening, nor does it attempt to model the spectral tion AV . This model has five free parameters: AV and the four shock
profile of the lines; where relevant in this work, dust reddening and parameters already introduced in § 3.2 (like for the prospector
line broadening are delegated to the caller program. model, we rescale the metallicity values to Z⊙ = 0.0142). We use a
This specific interpolator can be easily generalised to use more or Bayesian approach, where the probability of observing each ratio is
less lines, limited only by their availability in 3MdB. Gaussian. The prior probability distribution of the shock parameters
To evaluate the accuracy of the interpolator, we compare the input are flat in log space, spanning the range of parameters from the
grids from 3MdB (dashed lines) to the grids predicted by shinbnzv grid. The prior of AV is flat between 0 ≤ AV ≤ 4 mag. The posterior
(solid lines), on the BPT diagram (Fig. B1). By construction, the probability is integrated using emcee (Foreman-Mackey 2016),
model predicts the emission-line fluxes exactly at every grid node; and at each step the likelihood is calculated using the shinbnzv
deviations between the input and interpolated ratios are due to the interpolator.
common practice of connecting linearly between grid nodes in the The results are shown in Fig. C1; the corner diagram displays the
2-d space of emission-line ratios. marginalised posterior distribution, exhibiting the known degener-
acy between the magnetic field strength B and the pre-shock density
n. Panel b shows the data (black circles) and the model prediction
(red diamonds); we also report the intrinsic, de-reddened line ratios
APPENDIX C: PURE SHOCK MODEL
as grey squares. The posterior parameter values are all reasonably
A shock-only model – without photoionisation due to star formation close to the prospector model from § 4.1. But despite this agree-
– cannot explain the line ratios observed in Ulema. Here we con- ment, it is clear that the shock model cannot satisfactorily reproduce
sider the shock models of § 3.2, and fit eight emission-line ratios, both the low-ionisation species and the strong Balmer decrement;
taking the measurements from Table 2: [O ii]λλ3726,3729/Hβ, the outcome is an under-predicted Hα/Hβ ratio. Note that this so-
[O iii]λλ4959,5007/Hβ, [O i]λλ6300,6364/Hβ, [N ii]λ6584/Hα, lution is driven primarily by the strong prior in AV ; allowing AV to
Hα/Hβ, [S ii]λλ6716,6731/Hα, [S ii]λ6716/[S ii]λ6731 and assume values higher than 2 mag, the model favours high dust atten-
[O ii]λλ7319–7331/[O ii]λλ3726,3729. We tested shock plus uation (thus matching the Balmer decrement) but at the expense of
precursor models, but found that shock-only models agreed much the [O ii]λλ3726,3729/Hβ ratio, which is then under-predicted. Be-
better with the data, as discussed in § 3.1. For this reason, here we tween these alternatives, it is clear that the low-attenuation solution
report just the shock-only models. Our setup consists in combining is the most credible, because the fraction of the Hα flux that is unex-
the emission-line ratios from the shock models with the Cardelli plained by the shock can easily be ascribed to dust-obscured star for-

MNRAS 000, 1–20 (2024)

 20 F. D’Eugenio et al.
(a) log B[µG] = 0.46 +0 .17
−0.12 [OI]λλ6300, 6364 (b) (rest) [ m]
log v [km s −1 ] = 2.84 +0 .09
−0.09
Hβ 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2 1.3 1.4
log ne [cm −3 ] = 0.37 +0 .28 [OII]λλ3726, 3729 1.0 (a) 1.5

19 erg s 1 cm 2 Å 1]
(c) G395H

[N I]
C II]

H &[O III]

[O I]
C III]

H &[N II]
[S II]
[O II]
−0.17
Hβ
log Z [Z ¯ ] = −0.03 +0

19 erg s 1 cm 2 Å 1]
.04
−0.04 [OIII]λλ4959, 5007 1.0
AV [mag] = 1.35 +0 .31
−0.23 Hβ 0.8
[OII]λλ7319 − 7331 0.5
] log v [km s −1 ]

[OII]λλ3726, 3729 0.6

.25 .50 .75 1.00 2.6 2.7 2.8 2.9

[S II]λλ6716, 6731 0.0

F [10
Hα 0.4 3.5 3.6 3.7 3.8 3.9
(observed) [ m]
[N II]λ6583
Hα 0.2

F [10
[S II]λ6716
−3

[S II]λ6731
0.0
log Z [Z ] 0log0 ne0 [cm

Hα/Hβ
2 0 2 2.5 (c) 1 2 3 4
(observed) [5 m]
(d)
5 6 7
G395H
8
Observed ratio 0.0 0
2.5 5
3.5 3.6 3.7 3.8 3.9
1.2 1.6 2.0 2.4 0.16 0.08 0.00 0.08
¯

1 2 3 4 5 6 7 8
Data (observed) [ m]
Model
No dust
Figure D1. Mock observation and re-fit of the fiducial model (Table 3 and
Fig. 5), with input parameters equal to fiducial model posterior, but setting
A [mag]

log L(Hβ)sh = 0.06. The symbols and panels are the same as Fig. A. Over-
all, the mock data (grey) is of similar quality as the WIDE observations
V

(cf. Fig. 5b–d). Notable differences are the absence of outliers in the G395H
spectrum (outliers are masked in the analysis), and the much higher K s -band
0.2
0.4
0.6
0.8

2.6
2.7
2.8
2.9

0.55
0.70
1.05
0
6
8
0
8
1.2
1.6
2.0
2.4
0.2

0.1
0.0
0.0
0.0

log B[µG] log v [km s −1 ] log ne [cm −3 ] log Z [Z ¯ ] AV [mag] flux.

1 σ, except for log Bsh and log nsh , which tend to have larger excur-
Figure C1. Marginalised posterior probability distribution for the shock-only
model to the observed emission-line ratios. Panel b shows the observed line
sions of 1.5 σ. Physically, the strongest constraint on these two pa-
ratios (black), predicted ratios (red diamonds) and de-reddened ratios (grey rameters comes from multi-level forbidden line sets. Unfortunately,
squares). Without star formation, shocks are unable to simultaneously re- none of these emission features are spectrally resolved in the prism
produce all the observed ratios (reduced χ2 = 3.9), due to the incompatible spectrum. In the G395H spectrum, the only detection useful for dis-
requirements of high observed [O ii]λλ3726,3729/Hβ and Hα/Hβ. Overall, entangling log Bsh and log nsh is the [S ii]λλ6716,6731 doublet, but
the best-fit shock parameters are very similar to the results of the prospector the pernicious combination of high dispersion and low signal-to-
composite model, which combines shocks and star formation. noise ratio prevents us from drawing strong constraints; indeed, the
posterior probability distribution of the [S ii]λ6716/[S ii]λ6731 ratio
in Fig. 3 spans the full allowed range.
mation. In contrast, the high-attenuation solution leaves unexplained Overall, the outcome of our tests confirm that the software imple-
the [O ii]λλ3726,3729 flux, which is hard or even impossible to pro- mentation is self consistent, and that the parameter recovery process
duce without also altering the other line ratios. does not introduce bias.
A final word of caution is that – in the presence of different phys-
ical mechanisms powering nebular emission – the Hα/Hβ ratio may
This paper has been typeset from a TEX/LATEX file prepared by the author.
not necessarily give an unbiased estimate of dust attenuation.

APPENDIX D: ELEMENTARY RECOVERY TESTS

To assess the ability of the software to infer the input parameters,
we conduct a set of simple recovery tests. Due to the computational
cost of these tests, we focus on a single model, the fiducial prospec-
tor model present in § 4.1. We create mock observations using the
fiducial model parameters listed in Table 3, and varying one pa-
rameter at a time, as follows: log L(Hβ)sh = 0.06, log vsh = 2.6,
log M⋆ = 11.3, log Zgas = −0.5, τV = 0.0, τV = 1.3, log ngas = 2,
and log Bgas = −1 (all units are as in Table 3). The mocks are gen-
erated by first creating the model spectrum, then matching the spec-
tral resolution and wavelength grid of the observations, and finally
adding random noise to mimic the noise of the input data. For pho-
tometry, we draw the mock data from a Gaussian distribution with
mean equal to the model flux and standard deviation equal to the
observed noise. For the spectra, noise is generated by adding the
residuals of the fiducial model to the mock spectrum; for each spec-
tral pixel, we randomly draw the residuals in a narrow wavelength
window around the pixel (but we do not use residuals identified as
outliers). An example mock model is shown in Fig. D1.
Overall, we find that the input parameters are recovered within

MNRAS 000, 1–20 (2024)