Radiative and mechanical energies in galaxies

J. A Villa-Vélez, Benjamin Godard, Pierre Guillard, Guillaume Pineau Des Forêts

► To cite this version:

J. A Villa-Vélez, Benjamin Godard, Pierre Guillard, Guillaume Pineau Des Forêts. Radiative and mechanical energies in galaxies. Astronomy & Astrophysics - A&A, 2024, 688, pp.A96. 10.1051/0004-6361/202449212. hal-04670415

HAL Id: hal-04670415

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

Submitted on 12 Aug 2024

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

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

<figure><img src="image_1.png" /></figure>
# Radiative and mechanical energies in galaxies

## I. Contributions of molecular shocks and PDRs in 3C 326 N

J. A. Villa-Vélez<fnref n="1" />, B. Godard<fnref n="2" />,<fnref n="1" />, B. P. Guillard<fnref n="3" />, and G. Pineau des Forêts<fnref n="4" />,<fnref n="2" />

Received 11 January 2024 / Accepted 2 May 2024

## ABSTRACT

**Context.** Atomic and molecular lines emitted from galaxies are fundamental tracers of the medium responsible for the emission and contain valuable information regarding the energy budget and the strength of the different feedback mechanisms.

**Aims.** The goal of this work is to provide a new framework for the interpretation of atomic and molecular lines originating from extragalactic sources and a robust method to deduce the mechanical and radiative energy budget from a set of observations.

**Methods.** Atomic and molecular lines detected in a given object are assumed to result from the combination of distributions of shocks and photo-dissociation regions (PDRs) within the observational beam. The emission of individual structures is computed using the Paris-Durham shock code and the Meudon PDR code over a wide range of parameters. The total emission is then calculated assuming probability distribution functions for shocks and PDRs. A distance between the observational dataset and the model is finally defined based on the ratios of the observed to the predicted intensities.

**Results.** As a test case scenario, we consider the radio galaxy 3C 326 N. The dataset is composed of 12 rotational and ro-vibrational lines of H₂ and the fine structure lines of C⁺ and O. Our interpretative framework shows that both shocks and PDRs are required to explain the line fluxes. Surprisingly, viable solutions are obtained at low density only ($n_H$ < 100 cm⁻<fnref n="3" />), indicating that most of the emission originates from diffuse interstellar matter. The optimal solution, obtained for $n_H$ = 10 cm⁻<fnref n="3" />, corresponds to a distribution of low-velocity shocks (between 5 and 20 km s⁻<fnref n="1" />) propagating in PDR environments illuminated by a UV radiation field ten times larger than that in the solar neighborhood. This solution implies that at least 4% of the total mass carried by the PDRs is shocked. The H₂ 0-0 S(0) 28 µm, [CII] 158 µm, and [OI] 63 µm lines originate from the PDR components, while all the other H₂ lines are mostly emitted by shocks. The total solid angles sustained by PDRs and shocks imply that the radiative and mechanical energies reprocessed by these structures are $L_{UV} = 6.3 \times 10^9 L_⊙$ and $L_K = 3.9 \times 10^8 L_⊙$, respectively, in remarkable agreement with the values of the IR luminosity deduced from the fit of the spectral energy distribution (SED) of 3C 326 N, and consistent with a small fraction of the active galactic nuclei (AGN) jet kinetic power dissipated in the interstellar medium (≈1%).

**Conclusions.** This work shows that the radiative and mechanical energy budget of galaxies can be derived from the sole observations of atomic and molecular lines. It reveals the unexpected importance of the diffuse medium for 3C 326 N, in contrast to previous studies. A last-minute comparison of the model to new JWST data obtained in 3C 326 N shows a striking agreement and demonstrates the ability of the model to make accurate predictions. This framework opens new prospects for the prediction and interpretation of extragalactic observations, in particular in the context of JWST observations.

**Key words.** shock waves – ISM: kinematics and dynamics – ISM: molecules – photon-dominated region (PDR) – galaxies: ISM – galaxies: jets

## 1. Introduction

Galaxies are dynamic systems where diverse energy sources and sinks shape their evolution over cosmic time. The intricate interplay between dynamics, energy dissipation, and star formation constitutes the essence of galaxy formation and evolution (e.g., Pfenniger 2010). For instance, in the hypothetical absence of radiative losses, the kinetic energy released by gravitational accretion onto halos would keep the gas unbound and galaxies would not form at all (Benson 2010). It is therefore of paramount importance to understand and quantify how and where radiative and mechanical energies are reprocessed and dissipated in the multiphase gas of galaxies.

The reprocessing of radiative and mechanical energy in the interstellar medium (ISM) regulates the distribution

of gas phases (e.g., Vazquez-Semadeni 2009), the formation of molecular clouds and stars (Klessen & Hennebelle 2010), and the feedback mechanisms (Cielo et al. 2018; Ostriker & Kim 2022). Mechanical processes encompass the momentum input from stellar winds (Krumholz et al. 2014), supernovae (Fierlinger et al. 2016), and active galactic nuclei (AGN) (Ciotti et al. 2010). Simultaneously, radiative processes involve cooling through thermal and nonthermal emissions across the electromagnetic spectrum, from radio waves to gamma rays. The dissipation of radiative and mechanical energy in the ISM is intrinsically linked to the gas heating and cooling mechanisms, governed by the microphysics. The thermodynamic and chemical states of the gas significantly influence ISM observables, including ionic, atomic, and molecular line emissions. Those gas properties depend on factors such as density,

<footnote n="1">¹ Laboratoire de Physique de l'École Normale Supérieure, ENS, Université PSL, CNRS, Sorbonne Université, Université Paris Cité, 75005 Paris, France
e-mail: jorge.villa@phys.ens.fr</footnote>

<footnote n="2">² Observatoire de Paris, PSL University, Sorbonne Université, LERMA, 75014 Paris, France</footnote>

<footnote n="3">³ Sorbonne Université, CNRS, UMR 7095, Institut d' Astrophysique de Paris, 98bis bd Arago, 75014 Paris, France</footnote>

<footnote n="4">⁴ Université Paris-Saclay, CNRS, Institut d'Astrophysique Spatiale, 91405 Orsay, France</footnote>
metal and dust abundance, and the strengths of heating and/or ionization sources such as the ultraviolet (UV) radiation field and cosmic rays.

Efforts to model line and continuum emission in galaxies provide a means to constrain dominant gas excitation mechanisms and to quantify the reprocessed energy required to reproduce observations. Nonetheless, the complexity of the ISM poses significant challenges in solving radiative transfer within galaxies (Kim et al. 2022). As a result, thermal equilibrium and energy balance arguments are frequently employed to simplify calculations. For instance, contemporary codes used to replicate the panchromatic spectral energy distributions (SEDs) of galaxies often assume that stellar UV and near-infrared (NIR) light is absorbed by interstellar dust and subsequently re-emitted in the far-infrared (FIR) domain (da Cunha et al. 2008; Carnall et al. 2018; Boquien et al. 2019), primarily in photo-dissociation regions (PDRs) situated at the interfaces between molecular clouds and diffuse regions of atomic and ionized gas. This assumption is commonly employed to constrain the amount of reprocessed radiative energy from FIR observations and deduce the properties of infrared (IR) galaxies (e.g., Casey et al. 2014).

The physical modeling of spectral line cooling from galaxies, arising from the de-excitation of ions, atoms, and molecules, typically focuses on specific chemical states of the gas (ionized, atomic, or molecular), ISM phases (particular temperature or density ranges), or energy types (electromagnetic radiation, cosmic rays, or mechanical dissipation). For instance, well-established photo-dissociation codes like CLOUDY (Ferland et al. 2017) or the Meudon PDR code (Le Petit et al. 2006) are commonly employed to model the reprocessing of radiative energy, thereby constraining excitation mechanisms and physical conditions of X-ray or UV-irradiated gas in galaxies (e.g., Polles et al. 2021). Conversely, the reprocessing of mechanical energy in ionized, atomic, or molecular media has been incorporated into codes such as MAPPINGS (Allen et al. 2008) or the Paris-Durham shock code (Flower 2003; Godard 2019). The majority of investigations employing the aforementioned models to characterize molecular line emissions from galaxies conclude that both radiative and mechanical processes contribute substantively to the observed gas excitation. As an illustration, the spectral line energy distribution (SLED) of the carbon monoxide (CO) molecule is frequently modeled through a mix of PDRs, X-ray dominated regions (XDRs), and shocks (Meijerink et al. 2013; Mingozzi et al. 2018; Esposito et al. 2022).

This paper addresses the reprocessing of energy in the molecular phase of the ISM. We present a method to estimate the relative contributions of shocks and PDRs to the observed line emission from galaxies, with a focus on H₂ lines. Given that H₂ can form on interstellar dust at elevated temperatures (≈50 to 250 K, Grieco et al. 2023) and that the H₂ rotational levels are separated by a few hundred Kelvin, it emerges as a significant coolant for warm gas and, consequently, a key participant in the reprocessing of radiation within PDRs (Habart et al. 2005; Shaw et al. 2009) and the dissipation of kinetic energy within molecular shocks (Flower 1999; Kristensen et al. 2023). A previous study by Lesaffre (2013) used probability distribution functions of shock velocities to model the rotational and ro-vibrational lines of H₂ in molecular gas heated by shock waves. These emission lines serve as indicators of mechanical energy dissipation (Lesaffre et al. 2020), particularly within galaxies where the total power of IR H₂ line emissions exceeds what can be accounted for by the

exclusive reprocessing of the available UV radiative luminosity (Guillard & Boulanger 2009; Ogle et al. 2010; Nesvadba et al. 2010; Herrera et al. 2012; Guillard et al. 2015a,b).

This paper presents an extension of the framework presented in Lehmann & Godard (2022), where we now combine distributions of shocks and PDR models. The state-of-the-art codes used to model the line emission are presented in Sect. 2. The toy model we use to describe the line fluxes from an ensemble of shocks and PDRs is described in Sect. 3, as is the methodology used to fit a suite of emission lines with grids of models. In Sect. 4 we apply this modeling to the 3C 326 N radio-galaxy, constraining UV- and shock-processed gas mass and energy budgets. In Sect. 5 and conclude in Sect. 6, we provide predictions of the suite of ro-vibrational H₂ emission lines that will be observed with the JWST. Throughout this work, we assume a modern flat cosmology (Planck Collaboration VI 2020) with $\Omega_m = 0.315$, $\Omega_\Lambda = 0.685$, $\Omega_b = 0.0493$, and $H_0 = 67.4 \text{ km s}^{-1} \text{ Mpc}^{-1}$.

## 2. Shock and PDR models

For this paper, two publicly available codes were used to model PDR and shock emissions: the Meudon PDR code<fnref n="1" /> (Le Petit et al. 2006) and the Paris-Durham shock code<fnref n="2" /> (Flower 2003; Godard 2019).

### 2.1. Grids of PDR models

The Meudon PDR code is a model designed to study the thermal and chemical structures of photo-dissociation regions with a constant proton density, $n_H$, or a constant thermal pressure. The geometry, the main assumptions, and the physical ingredients included in the code are described in Le Petit et al. (2006). In a nutshell, the model considers a 1D plane-parallel and static slab of gas and dust illuminated on both sides by semi-isotropic radiation fields, each of which is set to the standard interstellar radiation field (ISRF, Mathis et al. 1983) and scaled with a parameter $G_0$. The size of the cloud is set by its total extinction $A_{v \max}$. Within this geometry, the model solves the radiative transfer, taking into account the detailed absorption and emission processes induced by dust and gas compounds, and the thermal and chemical states of matter at equilibrium, through the inclusion of an extensive chemical network.

In addition to the chemical composition of the gas, the model also solves the excitation of atomic and molecular species including the populations of 250 ro-vibrational levels of H₂. In particular, the excitation of the fine structure levels of C⁺ and O is computed including excitation by collisions with H, H₂, and electrons (for C⁺) and with H, H⁺, H₂, He, and electrons (for O). The populations of the ro-vibrational levels of H₂ are computed taking into account excitation by collisions with H, He, H₂, and H⁺, the probability of exciting H₂ at formation on grain surfaces, and the excitation by radiative pumping of the electronic lines of molecular hydrogen followed by fluorescence. The emerging line fluxes of all these species are finally calculated by integrating the emissivities over the entire slab up to $A_{v \max}$.

The grid of models used in this paper is a subsample of the precomputed models obtained with the 1.5.4 version of the Meudon PDR code and publicly available on the Interstellar Medium DataBase service (ISMDB<fnref n="3" />). We consider PDR models at constant proton density only, with a left-side illumination

<footnote n="1">¹ Meudon PDR code: https://ism.obspm.fr/pdr.html</footnote>

<footnote n="2">² Paris-Durham code: https://ism.obspm.fr/shock.html</footnote>

<footnote n="3">³ Available at: http://ism.obspm.fr/index.html</footnote>
<table><caption>**Table 1.** PDR and shock model parameters.</caption><thead><tr><th>PDR parameter</th><th>Symbol</th><th>Value</th><th>Units</th></tr></thead><tbody><tr><td>Proton density</td><td>n<sub>H</sub></td><td>10, 10<sup>2</sup>, 10<sup>3</sup>, 10<sup>4</sup></td><td>cm<sup>-3</sup></td></tr><tr><td>Radiation field</td><td>G<sub>0</sub></td><td>1, 10</td><td>—</td></tr><tr><td>Total visual extinction</td><td>A<sub>v max</sub></td><td>10</td><td>mag</td></tr><tr><td>H<sub>2</sub> CR ionization rate</td><td>ζ<sub>H<sub>2</sub></sub></td><td>10<sup>-16</sup></td><td>s<sup>-1</sup></td></tr><tr><th>Shock parameter</th><th>Symbol</th><th>Value</th><th>Units</th></tr><tr><td>Proton density</td><td>n<sub>H</sub></td><td>10, 10<sup>2</sup>, 10<sup>3</sup>, 10<sup>4</sup></td><td>cm<sup>-3</sup></td></tr><tr><td>Radiation field <sup>(a)</sup></td><td>G<sub>0</sub></td><td>0, 1, 10</td><td>—</td></tr><tr><td>Initial visual extinction</td><td>A<sub>v</sub></td><td>1.0</td><td>mag</td></tr><tr><td>Shock velocity</td><td>v<sub>s</sub></td><td>1–40 in steps of 1</td><td>km s<sup>-1</sup></td></tr><tr><td>Magnetic parameter <sup>(b)</sup></td><td>b</td><td>0.1, 1.0</td><td>—</td></tr><tr><td>H<sub>2</sub> CR ionization rate</td><td>ζ<sub>H<sub>2</sub></sub></td><td>10<sup>-16</sup></td><td>s<sup>-1</sup></td></tr></tbody></table>

**Notes.** <sup>(a)</sup>Scaling factor of the standard ultraviolet radiation field of Mathis et al. (1983). <sup>(b)</sup>Transverse magnetic field $B_0 = b(n_H/cm^{-3})^{1/2}$ µG.

scaling factor $G_0 = 1$ or 10 and a right-side illumination scaling factor $G_0 = 1$. The range of parameters explored here is shown in Table 1.

## 2.2. Grids of shock models

The Paris-Durham shock code is a model designed to study the dynamical, thermal, and chemical evolution of molecular shocks propagating at a velocity $v_s$ in a magnetized and irradiated environment with an initial proton density $n_H$ and an initial magnetic field strength set by a magnetic parameter $b$. The geometry, the main assumptions, and the physical ingredients included in the code are described in Godard (2019). In a nutshell, the model considers a 1D plane-parallel shock wave at steady-state traveling in a photo-dissociation region. The shock is assumed to be irradiated upstream by a semi-isotropic radiation field set to the standard interstellar radiation field (ISRF, Mathis et al. 1983) and scaled with a parameter $G_0$, and to be initially located at a position within the PDR set by a visual extinction $A_v = 1.0$ (Kristensen et al. 2023). Within this geometry, the model follows in a Lagrangian frame, the out-of-equilibrium evolution of a fluid particle during its trajectory from the pre-shock medium to the post-shock gas, starting from initial conditions at thermal and chemical equilibrium.

In addition to the chemical composition of the gas, the model solves the excitation of atomic and molecular species and the out-of-equilibrium evolutions of the first 150 ro-vibrational levels of $H_2$. In particular, the excitation of the fine structure levels of $C^+$ and O and of the ro-vibrational levels of $H_2$ are computed taking into account all the processes included in the Meudon PDR code and described in the previous section. The emerging line fluxes of all these species are finally calculated by integrating the emissivities over the entire shock up to the point where 99.9% of the input mechanical energy is dissipated (see Eq. (1) of Kristensen et al. 2023). This cutting point is chosen to prevent the inclusion of large amounts of cold post-shocked gas in the computation of line intensities and to ensure that emerging intensities mostly come from the reprocessing of mechanical energy input and not from other sources of energy (UV and IR fields, cosmic rays).

Although large grids of models obtained with the Paris-Durham shock code have been recently computed by Kristensen et al. (2023) and are publicly available on the ISMDB

service, we decided here to adopt a more refined set of shock velocities. The values of the parameters covered by the grid are shown in Table 1. The shock velocity ranges between 1 and 40 km s⁻¹ in steps of 1 km s⁻¹.

## 3. Ensembles of shocks and PDRs

When an extragalactic object is observed, the instrumental beam captures the emission from all the phases of the ISM. As schema-tized in Fig. 1, these include the hot and diffuse ionized gas, the diffuse and irradiated neutral media, and the dense molecular gas. In principle, all these phases contribute to the reprocessing of the radiative energy input originating from the AGN and the stars, and the mechanical energy input originating, for instance from the interaction between the jet and the galactic disk. Since we focused on molecular lines and on typical tracers of the neutral interstellar medium, we neglected the contribution of the hot ionized phase. Under this approximation, we developed a simplified model in which the reprocessing of the input UV radiative energy occurs in idealized 1D PDRs, and the reprocessing of the input mechanical energy occurs in idealized 1D shock structures.

### 3.1. Toy model

As described in Appendix A, we consider that an ensemble of shock and PDR surfaces, subtending total solid angles $\Omega_S$ and $\Omega_P$, are captured within an observational beam of solid angle $\Omega_{obs}$ (see Fig. 1). The cumulative fluxes originating from shock and PDR structures collectively contribute to the observed emission lines. For simplicity, the total solid angles $\Omega_S$ and $\Omega_P$ are assumed to follow 1D probability distribution functions denoted by $f_S$ and $f_P$. The shock distribution is supposed to depend on the shock velocity $v_s$ and the PDR distribution is supposed to depend on the illumination factor $G_0$, yielding the specific functional forms $f_S(v_s) = d\Omega_S/dv_s$ and $f_P(G_0) = d\Omega_P/dG_0$, respectively.

Considering that the continuum absorption due to the shock and PDR surfaces is negligible and that no cross-absorption occurs between the surfaces<fnref n="4" /> (see Eqs. (A.8) and (A.9)), we derive a simple expression of the continuum subtracted line flux (see Eq. (A.11)). It follows that the integrated intensity of each line (in erg s⁻¹ cm⁻² sr⁻¹) is given by

$$I^m = \frac{1}{\Omega_{\text{obs}}} \left( \int_{v_s} f_S(v_s) I_S(v_s) dv_s + \int_{G_0} f_P(G_0) I_P(G_0) dG_0 \right), \quad (1)$$

where $I_S(v_s)$ and $I_P(G_0)$ are the line integrated intensities emitted by one shock of velocity $v_s$ and one PDR of illumination factor $G_0$, respectively. The superscript “m” stands for model.

It should be noted that the framework can be easily extended to other types of contributions to the line emission, for instance the XDR and HII region. In such cases, the right-hand side of Eq. (1) can be partitioned into multiple components.

### 3.2. Comparison between model and observations: Distance minimization

This toy model can be used to interpret an ensemble of integrated intensities of unresolved emission lines. The confrontation of the model to the observations only requires defining a

<footnote n="4">⁴ Despite this assumption, the radiative transfer inside each shock and PDR structure is computed taking into account absorption processes.</footnote>
<figure><img src="image_2.png" /><figcaption>Fig. 1. Schematic representation of the multiphase ISM in a galaxy. Left: radio emission observed in 3C 326 N with the Westerbork Synthesis Radio Telescope (WSRT) 21 cm (Willis & Strom 1978) and Spitzer/IRAC composite image (Ogle et al. 2007) at 3.6, 4.5, and 8.0 µm, showing the 3C 326 N galaxy (and its southern companion). The 3.7" × 3.7" (6.4 × 6.4 kpc) square (not-to-scale) represents the solid angle Ωobs over which line emissions are observed. Center: cutout showing the projected proton density from a numerical simulation of the interaction between a radio jet and the central part of a multiphase galactic disk (from Murthy et al. 2022). The different colors show the parsec-scale multiphase structure, with the dense irradiated medium in red (nH ≥ 10³ cm⁻³), the shock diffuse gas in yellow (i.e., nH ~ 10⁻³ cm⁻³), and the hot X-ray plasma in green (nH ~ 1 cm⁻³). Right: simplified model used in this paper where the observational beam is supposed to encompass distributions of shock and PDR surfaces (see Fig. A.1).</figcaption></figure>

distance between the modeled and the observed intensities. Minimizing this distance leads to the characterization of the distributions of shocks and PDRs within the beam, in particular the number of shocks and PDRs and the distribution of shock velocities required to reproduce the observations.

We consider an ensemble of N emission lines. According to Eq. (1), the modeled intensity of each line j writes

$$I_j^m = \frac{1}{\Omega_{\text{obs}}} \left( \int_{v_s} f_S(v_s) I_{j,s}(v_s) dv_s + \int_{G_0} f_P(G_0) I_{j,p}(G_0) dG_0 \right). \quad (2)$$

We define the distance between the modeled and the observed intensities of each line j as the difference between the model prediction and the observation in log space

$$d_j = \begin{cases} 0 & \text{if } I_j^o - \frac{W_j}{2} \le I_j^m \le I_j^o + \frac{W_j}{2} \\ \log_{10}(I_j^m) - \log_{10}\left(I_j^o - \frac{W_j}{2}\right) & \text{if } I_j^m < I_j^o - \frac{W_j}{2} \\ \log_{10}(I_j^m) - \log_{10}\left(I_j^o + \frac{W_j}{2}\right) & \text{if } I_j^m > I_j^o + \frac{W_j}{2} \end{cases}, \quad (3)$$

where I_j^o is the observed intensity of line j and W_j is the associated 1σ uncertainty. The total quadratic distance between the model and the ensemble of observations is calculated as the sum of the quadratic distances of all emission lines

$$d = \sqrt{\sum_{j=1}^{N} d_j^2}. \quad (4)$$

Equation (4) is nothing more than the distance, in log scale, between the model and a hypercube in N dimensions defined by the observed intensities and their associated uncertainties. This distance also bears a natural physical meaning. For instance, a distance of 1 dex implies that at most one observation is underestimated or overestimated by a factor of 10.

With this definition of the distance, the interpretation of an ensemble of atomic and molecular lines is performed using the following three-step procedure:

First, the 1D distribution functions f_S(v_s) and f_P(G_0) are assumed to have known functional forms. Based on the results of numerical studies (e.g., Lehmann et al. 2016), we adopt by default here an exponential distribution of shock velocities: f_S(v_s) = (Ω_S/σv_s) exp(-v_s/σv_s), where σv_s is the dispersion of shock velocities. For simplicity, we consider a Dirac delta distribution of the PDR illumination factor: f_P(G_0) = Ω_Pδ(G_0).

Second, the modeled intensities are computed using Eq. (2) over a wide range of the parameters Ω_S, σv_s, and Ω_P describing the distributions f_S and f_P, assuming that all the other parameters (e.g., nH, b; see Table 1) are constants. The associated distances d(Ω_S, σv_s, Ω_P) between the model and the set of observations are computed with Eq. (4).

Finally, the global and local minima of the distance are searched for to identify the distributions f_S and f_P that fall the closest to the observational dataset.

### 3.3. Mechanical an UV-reprocessed luminosities

The probability distribution functions f_S and f_P directly provide the amount of mechanical and radiative energies reprocessed by shocks and PDRs. The rate of dissipation of mechanical energy of a steady-state shock propagating at a velocity v_s in a medium with pre-shock mass density ρ_0 is

$$\mathcal{L}_K(v_s) = \frac{1}{2} \rho_0 v_s^3 \mathcal{A}, \quad (5)$$

where A stands for the shock surface area (see Fig. A.1). The total mechanical luminosity dissipated (in erg s⁻¹ sr) by a distribution of shocks at different velocities therefore writes

$$\mathcal{L}_K = \frac{1}{2} \rho_0 \left( \int_{v_s} f_S(v_s) v_s^3 dv_s \right) D_A^2, \quad (6)$$

where D_A is the angular diameter distance.

As explained in Sect. 2.1, a PDR is modeled as a plane-parallel cloud illuminated by a semi-isotropic ISRF with a fixed scaling factor on the backside and a varying radiation factor G_0 on the front side. The integrated flux of the Mathis et al. (1983)
ISRF (in erg cm⁻² s⁻¹) from 912 to 2400 Å (Le Petit et al. 2006) implies that the UV-reprocessed power by one PDR is

$$ \mathcal{L}_{UV}(G_0) = 1.92 \times 10^{-3} (G_0 + 1.0) \mathcal{A}. \qquad (7) $$

The total reprocessed UV-luminosity for an ensemble of PDRs that follow a 1D probability distribution $f_P(G_0)$, therefore writes

$$ \mathcal{L}_{UV} = 1.92 \times 10^{-3} (G_0 + 1.0) \left( \int_{G_0} f_P(G_0) dG_0 \right) D_A^2. \qquad (8) $$

### 3.4. Shock and PDR masses

The probability distribution functions $f_S$ and $f_P$ also directly provide the total masses of shocks and PDRs within the observational beam. The total mass of shocked material writes

$$ M_S = 1.4 m_H \left( \int_{v_S} f_S(v_S) N_H(v_S) dv_S \right) D_A^2, \qquad (9) $$

where $m_H$ is the hydrogen mass and $N_H(v_S)$ is the proton column density across a shock of velocity $v_S$. The total mass of PDR material depends on the typical PDR size or equivalently their typical visual extinction $A_v$. This mass writes

$$ M_P(A_v) = 1.4 m_H \left( \int_{G_0} f_P(G_0) dG_0 \right) N_H(A_v) D_A^2, \qquad (10) $$

where $N_H(A_v)$ is the proton column density evaluated at a given visual extinction $A_v$.

## 4. Application to the radio galaxy 3C 326 N

Active galaxies are natural sources to investigate the importance of mechanical heating (Begelman & Ruszkowski 2005). As a test-case scenario, we therefore apply the framework developed in Sect. 3 to the radio galaxy 3C 326 N. This particular choice is motivated by the unique features of 3C 326 N, which include a notably low star formation rate (SFR) (Ogle et al. 2010, ~0.07 $M_\odot$ yr⁻¹), a substantial reservoir of molecular content (Nesvadba et al. 2011, ~2 × 10⁹ $M_\odot$), and the conspicuous presence of H₂ emission lines (see Appendix B). The low SFR in this galaxy implies that the gas is weakly irradiated by the stellar UV radiation field, and therefore favors the contributions of shocks in line emission (Guillard et al. 2012, 2015b). As described in Appendix B, we consider a dataset composed of 14 emission lines, including H₂ pure-rotational and ro-vibrational transitions, alongside [CII] 158 µm and [OI] 63 µm emission lines. The observed fluxes and uncertainties used for the distance minimization (see Eq. (3)) are listed in Table B.2.

The grids of PDR and shock models are described in Table 1. Any observed emission line is assumed to result from the combination of PDRs and shocks within the observational beam.<fnref n="5" /> For simplicity, we consider that all PDRs and shocks have a unique and common density and pre-shock density $n_H$, and are illuminated by a unique and common UV radiation field characterized by a scaling factor $G_0$. As explained in Sect. 2, we also consider that shocks propagate in PDRs at a typical visual extinction<fnref n="6" /> $A_v = 1.0$. Finally, we assume that all shocks propagate in

a medium with a constant transverse magnetic field strength $B_0$, scaled by the magnetization parameter $b$ that follows the relation $B_0 = b(n_H/cm^{-3})^{1/2}$ µG.

As a fiducial model, we adopt $n_H = 10$ cm⁻³, $b = 0.1$, and $G_0 = 10$. The choice of density comes from the exploration of the parameters and the comparison with the observations shown in the following sections. The choice of the illumination factor comes from the measurements of the UV emission (Ogle et al. 2010) and of the spectral SED of dust in the IR (Guillard et al. 2015b), which suggest $G_0 = 6_{-2}^{+4}$ and $G_0 = 9 \pm 1$, respectively.

In this section we first present the minimization process and the optimal shock and PDR distributions obtained for the fiducial model, and then show how this optimal distribution depends on the fixed parameters. In principle, the effects of all fixed parameters should be discussed. However, we find that both $G_0$ and $b$ have a minor impact on the distribution of shocks and PDRs required to explain the observations as long as $1 \le G_0 \le 10$ and $b \le 1$. We therefore limit our discussion to the effect of the gas density only.

### 4.1. Likelihood analysis of the fiducial model

The results of the minimization procedure applied to the fiducial model are displayed in Fig. 2, which shows the distance computed with Eqs. (3) and (4) as a function of the parameters of the shock and PDR distributions, $\Omega_S$, $\sigma_{v_S}$, and $\Omega_P$. This analysis highlights several key features.

First and foremost, we find that the fiducial model leads to distances as small as 0.33 dex. As we show in Fig. 3, the value of the global minimum implies that the 14 observational constraints, which include atomic and molecular emissions, are all reproduced within a factor smaller than 2, a surprising result given the simplicity of the adopted distributions and the 1D plane-parallel geometry of the models. The global minimum is obtained for distributions of shocks and PDRs with total solid angles $\Omega_S = 1.9 \times 10^{-8}$ sr and $\Omega_P = 1.2 \times 10^{-9}$ sr, and a shock velocity distribution characterized by $\sigma_{v_S} \approx 4$ km s⁻¹. This result implies that there are about ten shock surfaces per PDR surface in the observational beam. It also indicates that the number of shocks rapidly decreases with the shock velocity in agreement with the distribution of shock velocities derived in magnetohydrodynamic (MHD) simulations of interstellar turbulence (see Fig. 9 of Lehmann et al. 2016). Interestingly, the optimal value of $\sigma_{v_S} \approx 4$ km s⁻¹ implies that the peak of dissipation of mechanical energy occurs for $v_S = 3\sigma_{v_S} \sim 12$ km s⁻¹ and that most of the dissipation is mediated by shocks with velocities between ~5 and ~20 km s⁻¹.

Another main feature of this analysis is the uniqueness of the optimal solution. As suggested by Fig. 2, we find that the entire parameter space is characterized by only one global minimum, which occupies a finite volume in the parameter space, and no local minima. This absence of local minima indicates that there are no complex degeneracies between the parameters of the shock and PDR distributions. In particular, the large distances obtained for small values of $\Omega_S$ or $\Omega_P$ imply that both shocks and PDRs are required to explain the observational dataset and that no acceptable solution (with a distance smaller than one dex) can be obtained with a sole distribution of shocks or PDRs. Nevertheless, Fig. 2 also reveals some degeneracy between $\Omega_S$ and $\sigma_{v_S}$: solutions with distances smaller than 0.5 dex can be obtained for $\Omega_S$ varying over one order of magnitude and $\sigma_{v_S}$ over a factor of 2. This degeneracy comes from the fact that the emission lines of H₂ mostly originate from shocks with a specific range of velocities (see Sect. 5) whose numbers are set by both $\Omega_S$ and $\sigma_{v_S}$.

<footnote n="5">⁵ We do not consider the heating of the molecular gas by the X-ray radiation from the AGN. This is not a strong limitation, as Ogle et al. (2010) showed that the H₂-to-X-ray luminosity ratio (>0.6) is so high that X-ray heating is not a primary driver of H₂ line emission.</footnote>

<footnote n="6">⁶ This choice of typical visual extinction has little impact on the results presented here because the predictions of the shock models regarding H₂ emission are insensitive to $A_v$ as long as the shocks run in a medium that is mostly molecular (Kristensen et al. 2023).</footnote>
<figure><img src="image_3.png" /><figcaption>**Fig. 2.** 2D cutoff of the likelihood distributions of the main parameters describing the shock and the PDR distributions obtained for the fiducial model. The parameters are the total shock solid angle $\Omega_S$, the dispersion of the shock velocities $\sigma_{v_s}$, and the PDR total solid angle $\Omega_p$. The color-coding represents the distance in dex units. The 0.5, 0.7, and 1.0 dex limits are highlighted with dashed black lines. The straight black dashed lines and the black points indicate the position of the global minimum. The gray areas represent the regions where the assumption on the radiative transfer regarding cross-absorption between surfaces breaks down. The criterion is more stringent for the [OI] 63 µm line (gray area above the white solid line) than for the [CII] 158 µm line (gray area above the white dashed line), and is always fulfilled for H$_2$ lines (see Appendix A).</figcaption></figure>

Figure 2 finally reveals the domain of validity of the radiative transfer algorithm. As detailed in Appendix A, the radiative transfer used in this work is valid only if Eqs. (A.8) and (A.9) are satisfied. The gray area displayed in Fig. 2 shows the regions where these conditions break down for different emission lines. Remarkably, we find that the global minimum obtained for the fiducial model falls within a valid zone. This result therefore justifies, a posteriori, the fact of using a simplified radiative transfer to study the emission of ensembles of shocks and PDRs in galaxies. The simplified radiative transfer, however, is not valid for all values of $\Omega_S$ and $\Omega_p$. The most stringent constraints on the validity of the radiative transfer are due to the 63 µm line of O and, to a lesser extent, to the 158 µm line of C$^+$. This feature simply comes from the high opacities of these two lines across shock and PDR surfaces. In practice, the gray areas should be analyzed

with a more sophisticated radiative transfer algorithm for [OI] emission (and if necessary [CII]) taking into account saturation effects induced by the cross-absorption between all the PDR and shock surfaces. However, such a treatment would be superfluous. As implied by the results shown in Sect. 4.3, these ranges of parameters lead to solutions in contradiction with the radiative and mechanical energy budgets of the galaxy.

## 4.2. Impact of the medium's density

The impact of the density of the pre-shock gas and of the PDRs on the interpretation of observations is shown in Fig. 3a, which displays the minimum distance between the observations and the distribution of shocks and PDR models for a density ranging from 10 to $10^4$ cm$^{-3}$. The smallest value of the minimum
<figure><img src="image_4.png" /></figure>

**Fig. 3.** Solution for a combination of shocks and PDRs, assuming an exponential probability distribution function. Panel (a): distance vs. the medium's density. The dashed line represents a 0.5 dex limit (i.e., a factor of 3 difference between observed and modeled line intensities) to guide the eye. Panel (b): derived mechanical-reprocessed luminosity from Eq. (6) as well as the UV-reprocessed luminosity from Eq. (8) as a function of density are shown in green and orange, respectively. The green area represents <1% of the estimated jet kinetic luminosity of 3C 326 N, while the orange area represents the SED estimation of the reprocessed UV-luminosity (see Appendix C).

distance is obtained for $n_H = 10 \text{ cm}^{-3}$, which corresponds to the fiducial model. The corresponding intensities of the 14 emission lines used as constraints predicted by this fiducial model are presented in Fig. 4.

As shown in Fig. 4, the combination of shocks and PDRs reproduces all the observational lines by a factor smaller than 2. The H₂ lines from 0-0 S(1) to 1-0 S(4) are mostly emitted by shocks. Conversely, the 0-0 S(1) line of H₂, the 158 µm line of C⁺, and the 63 µm line of O are mostly emitted by PDRs.

Figure 3a shows that the minimum distance rapidly increases with the gas density. For instance, a minimum distance larger than ~1.4 dex is obtained for $n_H = 10^3 \text{ cm}^{-3}$, meaning that at most one of the observational lines is underestimated or overestimated by a factor of ~25. This result is in contradiction with

those currently presented in the literature. Fits of the excitation diagram of the H₂ pure-rotational transitions with shock models performed by Nesvadba et al. (2010) for 3C 326 N lead to solutions with a density ranging between $10^3$ and $10^4 \text{ cm}^{-3}$. These discrepancies come from the fact that our observational dataset contains additional constraints, in particular the [CII] 158 µm, and [OI] 63 µm lines. As shown in Fig. D.2, pure shock models cannot reproduce the [OI] line emission, even at low density. Moreover, in Nesvadba et al. (2010) a density of $10^4 \text{ cm}^{-3}$ was needed to reproduce high-excitation lines (e.g., 0-0 S(7)), but consequently, the 0-0 S(0) line was overpredicted by a factor of ≈2.5. The fact that the [CII] 158 µm line originates from the diffuse gas is in agreement with Guillard et al. (2015b) who found similar results based on the observation of the [CII]-to-FIR ratio in 3C 326 N. It is also in line with Pathak et al. 2024 who recently showed that the diffuse medium has a substantial contribution to the overall emission of the dust in main sequence galaxies.

### 4.3. Reprocessed radiative and mechanical energies

This new interpretation is valid only if the structures required to explain the observed intensities are in agreement with the energy budget of the source. This sanity check, which is often neglected, is paramount. The method described in Sect. 3.3 allows the computation of the mechanical energy dissipated by shocks $\mathcal{L}_K$, as well as the UV-reprocessed luminosity in PDRs $\mathcal{L}_{UV}$. The resulting luminosities are shown in Fig. 3b as functions of the pre-shock and PDR densities, and compared with the values estimated for 3C 326 N. For the fiducial model ($n_H = 10 \text{ cm}^{-3}$, $b = 0.1$, and $G_0 = 10$), we derive a dissipated mechanical power of $\mathcal{L}_K = 3.9 \times 10^8 L_\odot$ and a UV-reprocessed luminosity of $\mathcal{L}_{UV} = 6.3 \times 10^9 L_\odot$. The mechanical energy dissipated by the shocks is therefore as high as ~6% of the UV-reprocessed luminosity.

The UV luminosity predicted by the model corresponds to the radiation emitted by OB stars, which is reprocessed by interstellar matter. This radiation is expected to be mostly absorbed by dust and re-emitted at IR wavelengths. To obtain an independent estimate of the amount of UV radiation reprocessed, we performed a fit of the total SED of 3C 326 N, including 14 bands covering the UV-to-FIR range (see Appendix C), using the Bayesian code CIGALE<fnref n="7" /> (Boquien et al. 2019; Yang et al. 2020, 2022). This treatment leads to a reprocessed UV luminosity $L_{UV} = (7.1 \pm 1.4) \times 10^9 L_\odot$. This estimate is in agreement with Lanz et al. (2016), and is shown as an orange shaded region in Fig. 3b. This figure displays a remarkable agreement between the UV luminosity derived from the SED and that required to explain the atomic and molecular lines for gas densities ranging between $10^{-3} \text{ cm}^{-3}$. Interestingly, at $n_H \gtrsim 10^4 \text{ cm}^{-3}$ the UV-reprocessed luminosity obtained with the distribution of shocks and PDRs drops. This is due to the fact that at high densities the 0-0 S(0) line is mainly emitted by low-velocity shocks instead of PDRs.

Estimation of the available mechanical energy is less straightforward. In radio-galaxies, such as 3C 326 N, the interplay between the jet and the ISM results in a transfer of energy and momentum from large scales (the kiloparsec-sized galaxies inflated by the AGN-driven outflow, Hardcastle & Croston 2020) to smaller scales through a turbulent cascade, leading to an isotropization of the energy deposited by the AGN feedback, with turbulent velocity dispersions on the order of several hundreds of

<footnote n="7">⁷ Available at: https://cigale.lam.fr/</footnote>
<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Solution for a combination of shocks and PDRs, assuming an exponential probability distribution function. Solution for the fiducial model $n_H = 10 \text{ cm}^{-3}$, $b = 0.1$, $G_0 = 10$ shown in Fig. 3a. Observations are presented as blue stars (see Table B.2) and modeled intensities are shown as open black circles. Individual shock and PDR contributions are shown as green and orange solid lines, respectively.</figcaption></figure>

km s⁻¹ in the diffuse ionized phase (Wittor & Gaspari 2020). This cascade leads to the formation of low-velocity shocks which may account for the dissipation of a substantial fraction of the turbulent mechanical energy in the molecular gas (e.g., Lehmann et al. 2016; Park & Ryu 2019; Richard et al. 2022). According to simulations, approximately 20–30% of the kinetic energy of the jet is deposited into the ISM, depending on the AGN power (e.g., Mukherjee et al. 2016). Most of this energy is used to expand the gas along the jet and to drive gas outflows on subkiloparsec scales (Mukherjee et al. 2018; Morganti et al. 2023; Krause 2023). Recent observations obtained with the JWST suggest that only a minor fraction of this energy (less than 1%) feeds interstellar turbulence (Pereira-Santaella et al. 2022). Nesvadba et al. (2010) estimated a power of the jet of ~$10^{11}$ L⊙ from radio observations. The upper limit computed as 1% of this luminosity is shown as a green shaded area in Fig. 3b. For the fiducial model, the kinetic energy reprocessed in low-velocity shocks is found to be in agreement with the small fraction of the mechanical energy of the jet deposited in the ISM in 3C 326 N. Our results confirm that the dissipation of only a small fraction of the kinetic energy of the jet can account for the line emission of molecular hydrogen.

## 4.4. Total masses and volume of shocks and PDRs

The masses of shocks and PDRs in the fiducial model can be estimated using the prescription described in Sect. 3.4. The typical depth of PDRs is not known, but since the UV-impinging power needs to be entirely reprocessed in IR radiation, the associated visual extinction must be larger than 1. We assume two values of the visual extinction of PDRs, $A_v = 1$ and $A_v = 10$. Using Eqs. (9) and (10), we derive a total mass of shocks, $M_S = 7.8 \times 10^7 M_\odot$, and a total mass of PDR, $M_P = 1.9 \times 10^8 M_\odot$ (for $A_v = 1$) and $M_P = 1.9 \times 10^9 M_\odot$ (for $A_v = 10$). This result implies that between 4 and 40% of the total mass carried by the PDRs is shocked.

This estimation is in agreement with previous estimations of the molecular content of 3C 326 N. Ogle et al. (2007) and Nesvadba et al. (2010) derive a total mass of H₂ of $1.1 \times 10^9 M_\odot$.

and $1.3–2.7 \times 10^9 M_\odot$, respectively, from the excitation diagram of H₂ observed in 3C 326 N. Similarly, Guillard et al. (2015b) estimate a total molecular mass of ~$2 \times 10^9 M_\odot$ from the modeling of the FIR emission of dust grains and the emission of the CO(1-0) line.

The volume occupied by the PDRs can be derived from their total surface area and compared with the volume of the galaxy. For the fiducial model ($n_H = 10 \text{ cm}^{-3}$), we find that the volume of PDRs is comprised between $0.55 \text{ kpc}^3$ (for $A_v = 1$) and $5.5 \text{ kpc}^3$ (for $A_v = 10$). The volume of the galaxy is uncertain. NIR IFU spectroscopy with the VLT/SINFONI instrument show that the H₂ ro-vibrational emission arises from a molecular disk of about 2 kpc in diameter (see Nesvadba et al. 2011, footnote 1 page 3). The IR continuum and optical images show dust emission over 4 kpc (see Fig. 3 in Ogle et al. 2010). In addition, we estimate that the scale height of the galactic disk is between about 0.2 kpc (see Guillard et al. 2015b) and 0.5 kpc (see Fig. 3 in Ogle et al. 2010). This implies that the molecular gas is contained in a volume between 5 and $50 \text{ kpc}^3$. These estimations show that the models with densities $n_H \ge 10 \text{ cm}^{-3}$ are coherent because they lead to a volume filling factor smaller than 1. Given the results shown in Fig. 3, the question might arise of why we do not explore gas at densities lower than $10 \text{ cm}^{-3}$. However, this additional constraint shows that such models are ruled out because the molecular gas would occupy a volume larger than that of the galaxy.

## 4.5. Prediction of the CO(1-0) emission

The observed CO(1-0) emission line intensity given in Table B.1 is not included as a constraint in our model. Low-J CO emission is known to arise from the cold diffuse component of molecular gas (Mingozzi et al. 2018). Given that the total mass of shocks is much lower than that of PDRs, it is therefore expected that the low-J lines of CO are mostly emitted by PDRs.

The distribution of PDRs obtained in the fiducial model leads to a flux of the CO (1-0) line of $6.5 \times 10^{-18} \text{ erg s}^{-1} \text{ cm}^{-2}$, a value in fair agreement, yet slightly higher than that observed
<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Contributions of the different velocities of the shock distribution to the total intensities emitted by shocks in the fiducial model. The dominant shock velocity $v_{s,dom}$ that contributes the most to powering each line is highlighted in orange. Only the range 1–20 km s⁻¹ is shown for readability, even though the models include velocities up to 40 km s⁻¹.</figcaption></figure>

in 3C 326 N (3.5 × 10⁻¹⁸ erg s⁻¹ cm⁻²). This result reinforces our interpretation of this galaxy in which the 0-0 S(0), [CII] 158 µm, and [OI] 63 µm lines are mostly emitted by PDRs, while the remaining H₂ lines are mostly emitted by shocks.

# 5. Model predictions for JWST observations

## 5.1. Impact of shock velocities to the H₂ lines

From the fitting of observed line emission in 3C 326 N (Fig. 3), it is clear that H₂ lines are mainly powered by shocks, except the lowest excitation 0-0 S(0) line where the PDR contributions dominate. For the fiducial model ($n_H$ = 10 cm⁻³, $b$ = 0.1, and $G_0$ = 10, see Sect. 4), we now investigate, among the distribution of shocks, which shock velocities are the dominant contributors to the emission of H₂ lines.

Figure 5 displays the modeled shock integrated line intensities for the suite of 14 emission lines as a function of the shock velocity. A summary of the dominant velocity as a function of the upper-level energy of the H₂ lines is presented in Fig. 6. For the lines accessible with JWST observations (i.e., 0-0 S(1) to the 0-0 S(7), and the ro-vibrational 1-0 S(1) to 1-0 S(4)), we show

that the dominant shock velocities are remarkably low, in the $v_{s,dom}$ range of 6 km s⁻¹–15 km s⁻¹. The typical length scales over which those shocks dissipate energy is on the order of 0.001–0.01 pc (Kristensen et al. 2023). As previous studies pointed out, this suggests that a small fraction of the feedback kinetic energy cascades to very small scales in the molecular phase, powering the H₂ line emission (Nesvadba et al. 2010).

On the more technical side, the shocks responsible for the H₂ line emission are J-type shocks, where J stands for jump (see Kristensen et al. 2023, for a discussion of the shock model grid, and in particular their Fig. 6). These shocks are strong emitters of ro-vibrational emission lines, contrary to C-type (continuous) shocks, which were favored in Nesvadba et al. (2010) with fewer observational constraints. Therefore, JWST/NIRSpec observations of 3C 326 N will be a stringent test of our model predictions.

## 5.2. Impact of the distribution of shock velocities

The impact of the functional form of the distribution of shock velocities is presented in Appendix D.3 where we show the comparison between the model and the observations assuming
<figure><img src="image_7.png" /><figcaption>**Fig. 6.** Shock velocity that contributes the most to the H₂ line emissions as a function of the upper-level energy (orange points in Fig. 5). H₂ pure-rotational (green) and ro-vibrational (red) transitions are presented. $v_{s,dom}$ is derived for the fiducial model $n_H = 10 \text{ cm}^{-3}$, $b = 0.1$, and $G_0 = 10$.</figcaption></figure>

a Gaussian distribution for $f_S$. Appendix D.3, and in particular Figs. D.3 and D.4, show that a Gaussian distribution leads to the same conclusions regarding the media responsible for the different line emissions and the energy budget required to reproduce the observations. Figure D.3 reveals that assuming a more complex distribution of shock velocities only leads to additional degeneracies. For instance, we find that a narrow distribution of shock velocities centered at $\mu_{v_s} = 12 \text{ km s}^{-1}$ leads to the same distance as a broad distribution centered at $\mu_{v_s} = 3 \text{ km s}^{-1}$. This implies that the current ensemble of observations is insufficient to discriminate between different distributions of shock velocities as long as the shocks that dominate the ro-vibrational emission of H₂ (between 6 and 15 km s⁻¹, see Fig. 5) are included.

### 5.3. H₂ line predictions for JWST

The Paris-Durham shock and the Meudon PDR codes solve for the populations of a large number of H₂ levels (see Sect. 2), and therefore provide a predicted spectrum over the entire wavelength range observed by the JWST instruments. Figure 7 displays the spectrum predicted by the fiducial model for the strongest 1000 lines of H₂, accounting for the different spectral resolutions of the NIRSpec and MRS spectrographs. We fit the resolving power of the G235H NIRSpec disperser<fnref n="8" /> over the wavelength range 1.7–3.1 µm, finding the resolving power relation $R(\lambda) = 1225\lambda - 185$. For the MIRI instrument, we use the relation $R(\lambda) = 4603 - 128\lambda$ reported in Argyriou et al. (2023) valid over the wavelength range 5–28 µm. Outside the validity range in wavelength of these relations, we use a uniform resolving power of $\lambda/\Delta\lambda = 2500$. We note that, although displayed on the synthetic spectrum, the ground-state 0-0 S(0) emission line

is not observable with MIRI/MRS because its sensitivity drops above 26.5 µm.

We predict that the 3C 326 N spectrum should be dominated by ro-vibrational lines with $v_{up} = 1$ around 2–3 µm, as well as the strong 0-0 S(3) and 0-0 S(5) mid-infrared lines. The forest of NIR $v_{up} \ge 1$ ro-vibrational lines is thus predicted to be a major cooling channel for the dissipation of kinetic energy.

Luckily, it came to our attention, while this paper was under peer review, that new data obtained with the NIRSpec and MIRI instruments aboard the JWST had just been published (Leftley et al. 2024). These data include 19 H₂ lines encompassing S, O, and Q series, as well as the pure-rotational lines 0-0 S(3), 0-0 S(5), and 0-0 S(6). It offered a unique and timely opportunity to test the predictive power of the model. The predictions of the fiducial model were taken as they were and directly compared to these new JWST data. Figure 8 displays the remarkable results that the fiducial model reproduces all the lines observed with the JWST within a factor of about 2. Also, as predicted by Fig. 7, the number of ro-vibrational lines present in the NIRSpec spectrum is abundant and is consistent with an emission powered by the reprocessing of mechanical energy from shocks.

## 6. Summary and conclusions

This paper presents a new framework for interpreting atomic and molecular emission lines observed in extragalactic sources. Emission lines collected within an observational beam are assumed to arise from distributions of shocks and PDRs. The comparison of the observations with the predictions of the model leads to the physical properties (e.g., density, UV illumination) of the emitting gas and provides the mechanical and radiative energy budget of the source. As a test case, we selected the radio galaxy 3C 326 N, a unique object with low SFR, large molecular content, and strong IR H₂ line emission. A total of 14 emission lines including eight H₂ rotational lines from 0-0 S(0) to 0-0 S(7) and four ro-vibrational lines from 1-0 S(1) to 1-0 S(4), as well as [CII] 158 µm and [OI] 63 µm emissions are used as observational constraints.

The framework developed in this paper is an oversimplification of the gas distribution in galaxies. Because the ISM is highly turbulent and contains a distribution of irradiation sources, both the density and $G_0$ should spread over a wide range of values. This diversity of ISM environments is obtained in numerical simulations (e.g., Mukherjee et al. 2018; Murthy et al. 2022) and is derived from the analysis of continuum emission (e.g., Pathak et al. 2024). We designed the model and the underlying distributions to reduce the number of free parameters. This model is not a representation of the actual state of the gas, but rather describes the mean physical conditions of the emitting medium. The remarkable result is that this oversimplified picture already leads to predictions that are in excellent agreement with the observations, and fulfills the constraints on the mechanical and radiative energy budgets. A unique minimum is found with a small degree of degeneracies between the parameters (namely $\Omega_S$ and $\sigma_{v_s}$). This implies that adopting a distribution of density or $G_0$, while more realistic, would only lead to additional degeneracies. The key results obtained within this framework are summarized as follows:

1. The analysis applied to 3C 326 N shows that the ensemble of observational constraints requires the combined contributions of both a distribution of shocks and a distribution of PDRs within the observational beam.

<footnote n="8">⁸ We use the wavelength calibration data from https://jwst-docs.stsci.edu/jwst-near-infrared-spectrograph/nirspec-instrumentation/nirspec-dispersers-and-filters</footnote>
<figure><img src="image_8.png" /><figcaption>**Fig. 7.** Predicted synthetic spectrum of H₂ in 3C 326 N obtained for the fiducial model ($n_H$ = 10 cm⁻³, *b* = 0.1, and $G_0$ = 10). The specific intensities of the strongest 1000 emission lines are shown as functions of the redshifted wavelength for *z* = 0.089. We applied the relations described in Sect. 5.3 for the resolving power of the NIRSpec and MIRI instruments. Outside the validity wavelength range, we fixed the resolving power to a value of $\lambda/\Delta\lambda$ = 2500. The NIRSpec and MIRI bandwidths correspond to the black solid lines. The wide-, medium-, and narrowband filters for NIRCam and MIRI are shown as gray and colored bars. Pure rotation lines are indicated in blue, while ro-vibrational transitions are shown in red and green.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Fig. 8.** Comparison of the intensities of the H₂ lines predicted by the fiducial model with the new JWST data observed by Leftley et al. 2024, which are corrected here by a factor (1 + z)³ (see Eq. (B.1)). Observations are shown as blue stars and the model predictions as open black circles. We emphasize that the model is obtained by fitting the *Spitzer* and VLT/SINFONI data, and not the JWST data.</figcaption></figure>

2. The predictions of the model are found to be highly dependent on the density of the gas and weakly dependent on the other parameters (the ambient radiation field and the magnetic field). Surprisingly, we find that only shocks and PDRs at low density ($n_H$ = 10–100 cm⁻³) can reproduce all the observational constraints. This is probably the most salient result of this study. It implies that both the mechanical and UV energy inputs are preferentially reprocessed in the low-density medium of this galaxy.

3. The optimal solution is obtained for $n_H$ = 10 cm⁻³. The likelihood, estimated from the ratios of the observed to the predicted intensities, is characterized by only one global minimum and no local minima. This global minimum provides

constraints on the total solid angle of PDRs, the total solid angle of shocks, and the range of shock velocities required to explain the observations. However, the suite of spectral lines used in this work is found to be insufficient to set apart the predictions obtained for an exponential or a Gaussian probability distributions of shock velocities.

4. The H₂ 0-0 S(0) 28 µm, [CII] 158 µm, and [OII] 63 µm lines are mostly emitted by PDRs. The resulting UV-reprocessed luminosity, $\mathcal{L}_{\text{UV}}$ = 6.3 × 10⁹ $L_\odot$, is found to be in remarkable agreement with the observed IR-luminosity.

5. Conversely, the rest of the H₂ lines from the 0-0 S(1) to the 1-0 S(4) are found to be emitted by low-velocity shocks
(5 < v_s < 20 km s⁻¹). The resulting mechanical luminosity dissipated by shocks, $\mathcal{L}_K = 3.9 \times 10^8 L_\odot$, is less than 1% of the estimated jet kinetic power, in agreement with the predictions of MHD simulations or recent JWST observations.

6. The total mass of shocks is $M_S = 7.8 \times 10^7 M_\odot$ and the total mass of PDRs is $M_P = 1.9 \times 10^9 M_\odot$. This budget is in line with previous estimations of the molecular content in 3C 326 N and implies that about 4% of the total mass carried by PDRs is shocked.

7. Finally, the predicted synthetic H₂ spectrum shows prominent ro-vibrational emission lines, originating from all the vibrational levels of H₂. This synthetic spectrum is found to be in excellent agreement with new JWST observations of 3C 326 N, underlining the ability of the model to make accurate predictions.

The new interpretative framework is applied here to a specific source to test the validity and relevance of the underlying method. In a companion paper, we will extend this methodology to a larger sample of objects (see, e.g., Fig. 7 in Guillard et al. 2012 and Fig. 12 in Ogle et al. 2024), such as star-forming galaxies, AGN, radio, or elliptical galaxies, and extract the radiative and mechanical energy budget of these sources from atomic and molecular line emissions.

The fact that most of the reprocessing of the input mechanical and radiative powers occurs in the diffuse gas in 3C 326 N is of particular importance. It remains to be seen whether this result holds for a larger sample of galaxies, which we will explore in a companion paper. It suggests that this methodology should be applied to ancillary interpretations of extragalactic sources. Previous studies (Meijerink et al. 2013; Mingozzi et al. 2018; Esposito et al. 2022) show that PDR and XDR models at high densities (e.g., $n_H \sim 10^4 \text{ cm}^{-3}$) could explain the CO, H₂, and [CII] 158 µm emission lines observed toward a few extragalactic sources. It would therefore be interesting to apply our methodology to those sources and estimate whether a combined distribution of shocks, PDRs, and XDRs at low densities could account for the observations with an acceptable energy budget. It is indeed timely because JWST is already observing H₂ and fine-structure emission lines in extragalactic sources, in particular in AGN (Pereira-Santaella et al. 2022; Buiten et al. 2024).

Acknowledgements. The research leading to these results has received funding from the European Research Council, under the European Community's Seventh framework Programme, through the Advanced Grant MIST (FP7/2017-2022, No 742719). The grid of simulations used in this work has been run on the computing cluster Totoro of the ERC MIST, administered by MesoPSL. We would also like to acknowledge the support from the Programme National "Physique et Chimie du Milieu Interstellaire" (PCMI) of CNRS/INSU with INC/INP co-funded by CEA and CNES. PG would like to thank the Institut Universitaire de France, the Centre National d'Etudes Spatiales (CNES), and the "Programme National de Cosmologie and Galaxies" (PNCG) of CNRS/INSU, with INC/INP co-funded by CEA and CNES, for their financial supports.

## References

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

Almeida, A., Anderson, S. F., Argudo-Fernández, M., et al. 2023, ApJS, 267, 44

Argyriou, I., Glasse, A., Law, D. R., et al. 2023, A&A, 675, A111

Begelman, M. C., & Ruszkowski, M. 2005, Phil. Trans. R. Soc. London Ser. A, 363, 655

Benson, A. J. 2010, Phys. Rep., 495, 33

Boquien, M., Burgarella, D., Roehlly, Y., et al. 2019, A&A, 622, A103

Bruzual, G., & Charlot, S. 2003, MNRAS, 344, 1000

Buiten, V. A., van der Werf, P. P., Viti, S., et al. 2024, ApJ, 966, 166

Carnall, A. C., McLure, R. J., Dunlop, J. S., & Davé, R. 2018, MNRAS, 480, 4379

Casey, C. M., Narayanan, D., & Cooray, A. 2014, Phys. Rep., 541, 45

Chabrier, G. 2003, PASP, 115, 763

Charlot, S., & Fall, S. M. 2000, ApJ, 539, 718

Cielo, S., Bieri, R., Volonteri, M., Wagner, A. Y., & Dubois, Y. 2018, MNRAS, 477, 1336

Ciotti, L., Ostriker, J. P., & Proga, D. 2010, ApJ, 717, 708

da Cunha, E., Charlot, S., & Elbaz, D. 2008, MNRAS, 388, 1595

Draine, B. T., Aniano, G., Krause, O., et al. 2014, ApJ, 780, 172

Esposito, F., Vallini, L., Pozzi, F., et al. 2022, MNRAS, 512, 686

Ferland, G. J., Chatzikos, M., Guzmán, F., et al. 2017, Rev. Mex. A&A, 53, 385

Fierlinger, K. M., Burkert, A., Normousi, E., et al. 2016, MNRAS, 456, 710

Flower, D. R., & Pineau des Forêts, G., 1999, MNRAS, 308, 271

Flower, D. R., & Pineau des Forêts, G., 2003, MNRAS, 343, 390

Godard, B., Pineau des Forêts, G., Lesaffre, P., et al. 2019, A&A, 622, A100

Grieco, F., Theulé, P., De Looze, I., & Dulieu, F. 2023, Nat. Astron., 7, 541

Guillard, P., Boulanger, F., & Pineau des Forêts, G., & Appleton, P. N., 2009, A&A, 502, 515

Guillard, P., Ogle, P. M., Emonts, B. H. C., et al. 2012, ApJ, 747, 95

Guillard, P., Boulanger, F., Lehnert, M. D., Appleton, P. N., & Pineau des Forêts, G. 2015a, in SF2A-2015:Proceedings of the Annual meeting of the French Society of Astronomyand Astrophysics, 81

Guillard, P., Boulanger, F., Lehnert, M. D., et al. 2015b, A&A, 574, A32

Habart, E., Walmsley, M., Verstraete, L., et al. 2005, Space Sci. Rev., 119, 71

Hardcastle, M. J., & Croston, J. H. 2020, New A Rev., 88, 101539

Herrera, C. N., Boulanger, F., Nesvadba, N. P. H., & Falgarone, E. 2012, A&A, 538, L9

Kim, J.-G., Gong, M., Kim, C.-G., & Ostriker, E. C. 2022, ApJS, 264, 10

Klessen, R. S., & Hennebelle, P. 2010, A&A, 520, A17

Krause, M. G. H. 2023, Galaxies, 11, 29

Kristensen, L. E., Godard, B., Guillard, P., Gusdorf, A., & Pineau des Forêts, G., 2023, A&A, 675, A86

Krumholz, M. R., Bate, M. R., Arce, H. G., et al. 2014, in Protostars and Planets VI, eds. H. Beuther, R. S. Klessen, C. P. Dullemond, & T. Henning, 243

Lanz, L., Ogle, P. M., Alatalo, K., & Appleton, P. N. 2016, ApJ, 826, 29

Leftley, J. H., Nesvadba, N. P. H., Bicknell, G., et al. 2024, A&A, in press, https://doi.org/10.1051/0004-6361/202449848

Lehmann, A., Federrath, C., & Wardle, M. 2016, MNRAS, 463, 1026

Lehmann, A., Godard, B., & Pineau des Forêts, G., Vidal-García, A., & Falgarone, E., 2022, A&A, 658, A165

Le Petit, F., Nehmé, C., Le Bourlot, J., & Roueff, E. 2006, ApJS, 164, 506

Lesaffre, P., Pineau des Forêts, G., Godard, B., et al. 2013, A&A, 550, A106

Lesaffre, P., Todorov, P., Levrier, F., et al. 2020, MNRAS, 495, 816

Lilly, S. J., & Longair, M. S. 1984, MNRAS, 211, 833

Mathis, J. S., Mezger, P. G., & Panagia, N. 1983, A&A, 128, 212

Meijerink, R., Kristensen, L. E., Weiß, A., et al. 2013, ApJ, 762, L16

Mingozzi, M., Vallini, L., Pozzi, F., et al. 2018, MNRAS, 474, 3640

Morganti, R., Murthy, S., Guillard, P., Oosterloo, T., & Garcia-Burillo, S. 2023, Galaxies, 11, 24

Mukherjee, D., Bicknell, G. V., Sutherland, R., & Wagner, A. 2016, MNRAS, 461, 967

Mukherjee, D., Bicknell, G. V., Wagner, A. Y., Sutherland, R. S., & Silk, J. 2018, MNRAS, 479, 5544

Murthy, S., Morganti, R., Wagner, A. Y., et al. 2022, Nat. Astron., 6, 488

Nesvadba, N. P. H., Boulanger, F., Salomé, P., et al. 2010, A&A, 521, A65

Nesvadba, N. P. H., Boulanger, F., Lehnert, M. D., Guillard, P., & Salome, P. 2011, A&A, 536, L5

Ogle, P., Antonucci, R., Appleton, P. N., & Whysong, D. 2007, ApJ, 668, 699

Ogle, P., Boulanger, F., Guillard, P., et al. 2010, ApJ, 724, 1193

Ogle, P. M., López, I. E., Reynaldi, V., et al. 2024, ApJ, 962, 196

Ostriker, E. C., & Kim, C.-G. 2022, ApJ, 936, 137

Park, J., & Ryu, D. 2019, ApJ, 875, 2

Pathak, D., Leroy, A. K., Thompson, T. A., et al. 2024, AJ, 167, 39

Pereira-Santaella, M., Álvarez-Márquez, J., García-Bernete, I., et al. 2022, A&A, 665, L11

Pfenniger, D. 2010, in Galaxies in Isolation: Exploring Nature Versus Nurture, eds. L. Verdes-Montenegro, A. Del Olmo, & J. Sulentic, ASP Conf. Ser., 421, 183

Planck Collaboration VI. 2020, A&A, 641, A6

Polles, F. L., Salomé, P., Guillard, P., et al. 2021, A&A, 651, A13

Richard, T., Lesaffre, P., Falgarone, E., & Lehmann, A. 2022, A&A, 664, A193

Shaw, G., Ferland, G. J., Henney, W. J., et al. 2009, ApJ, 701, 677

Stalevski, M., Ricci, C., Ueda, Y., et al. 2016, MNRAS, 458, 2288

Vazquez-Semadeni, E. 2009, ArXiv e-prints [arXiv:0902.0820]

Willis, A. G., & Strom, R. G. 1978, A&A, 62, 375

Wittor, D., & Gaspari, M. 2020, MNRAS, 498, 4983

Yang, G., Boquien, M., Buat, V., et al. 2020, MNRAS, 491, 740

Yang, G., Boquien, M., Brandt, W. N., et al. 2022, ApJ, 927, 192
# Appendix A: Conditions of validity of the toy model

The toy model presented in Sect. 3 and used throughout this paper results from a simplification of the radiative transfer of an ensemble of shocks and PDRs within an observational beam. As illustrated in Fig. A.1, we consider that an observational beam of solid angle $\Omega_{\text{obs}}$ contains shock and PDR surfaces with a distribution of systemic velocities and with total solid angles $\Omega_S$ and $\Omega_P$. For simplicity, we set that the emission from one shock or one PDR in a line with a rest frequency $v_0$ follows a tophat profile over a width $\Delta v = v_0\Delta V/c$, with a velocity width $\Delta V = 5$ km s$^{-1}$ dominated by turbulent motions (Lehmann & Godard 2022). Both shock and PDR surfaces are assumed to have systemic velocities projected along the observed line of sight that follow a Gaussian distribution in frequency space

$$ \phi_v = \frac{1}{\sqrt{2\pi}\sigma_v^{\text{obs}}} \exp\left(-\frac{(v-v_0)^2}{2\sigma_v^{\text{obs}} 2}\right), \qquad (A.1) $$

where $\sigma_v^{\text{obs}}$ is the observed frequency dispersion of the line. Independently of this distribution, we also assume that $\Omega_S$ and $\Omega_P$ follow 1D probability distribution functions, $f_S(\Theta_S)$ and $f_P(\Theta_P)$, which both depend on a unique parameter, $\Theta_S$ for the shocks and $\Theta_P$ for the PDRs.<fnref n="9" /> The solid angles of shock and PDR surfaces that emit within a frequency bin of size $\Delta v$ centered at a frequency $v$ therefore write

$$ \Delta\Omega_S(v) = \Omega_S\phi_v\Delta v = \int_{\Theta_S} f_S(\Theta_S)d\Theta_S\phi_v\Delta v \qquad (A.2) $$

and

$$ \Delta\Omega_P(v) = \Omega_P\phi_v\Delta v = \int_{\Theta_P} f_P(\Theta_P)d\Theta_P\phi_v\Delta v. \qquad (A.3) $$

The number of shock and PDR layers within a frequency bin or, equivalently, the beam-filling factors per frequency bin of shocks and PDRs are

$$ \frac{\Delta\Omega_S}{\Omega_{\text{obs}}}(v) \qquad (A.4) $$

and

$$ \frac{\Delta\Omega_P}{\Omega_{\text{obs}}}(v). \qquad (A.5) $$

Their average total opacities are

$$ \bar{\tau}_S(v) = \int_{\Theta_S} f_S(\Theta_S) \tau_S(\Theta_S) d\Theta_S \phi_v \Delta v / \Omega_{\text{obs}} \qquad (A.6) $$

and

$$ \bar{\tau}_P(v) = \int_{\Theta_P} f_P(\Theta_P) \tau_P(\Theta_P) d\Theta_P \phi_v \Delta v / \Omega_{\text{obs}}, \qquad (A.7) $$

where $\tau_S(\Theta_S)$ and $\tau_P(\Theta_P)$ are the opacities across one shock layer of parameter $\Theta_S$ and one PDR layer of parameter $\Theta_P$ respectively (see Fig. A.2).

If there are no spatial overlaps of shocks and PDRs within any frequency bin, that is if

$$ \frac{\Delta\Omega_S}{\Omega_{\text{obs}}}(v_0) \ll 1 \quad \text{and} \quad \frac{\Delta\Omega_P}{\Omega_{\text{obs}}}(v_0) \ll 1, \qquad (A.8) $$

or if the averaged line opacities across the shock and PDR layers in any frequency bin are small,

$$ \bar{\tau}_S(v_0) \ll 1 \quad \text{and} \quad \bar{\tau}_P(v_0) \ll 1, \qquad (A.9) $$

absorption by shocks and PDRs can be neglected. The density flux (in erg cm$^{-2}$ s$^{-1}$ Hz$^{-1}$) collected in the observational beam therefore simply writes

$$ F_v = I_v^c \Omega_{\text{obs}} + \int_{\Theta_S} f_S(\Theta_S) \phi_v I_S(\Theta_S) d\Theta_S + \int_{\Theta_P} f_P(\Theta_P) \phi_v I_P(\Theta_P) d\Theta_P, \qquad (A.10) $$

where $I_v^c$ is the specific intensity of the background continuum (in erg s$^{-1}$ cm$^{-2}$ Hz$^{-1}$ sr$^{-1}$) and $I_S(\Theta_S)$ and $I_P(\Theta_P)$ are the line integrated intensities emitted by one shock of parameter $\Theta_S$ and one PDR of parameter $\Theta_P$, respectively (in erg s$^{-1}$ cm$^{-2}$ sr$^{-1}$). Integrating this equation over frequency finally gives the continuum subtracted line flux (in erg s$^{-1}$ cm$^{-2}$)

$$ F = \int_{\Theta_S} f_S(\Theta_S) I_S(\Theta_S) d\Theta_S + \int_{\Theta_P} f_P(\Theta_P) I_P(\Theta_P) d\Theta_P. \qquad (A.11) $$

This last equation is at the root of the toy model presented in Sect. 3 (Eq. 1). It results from the assumption that the absorption induced by shock and PDR surfaces within the observational beam is negligible. Throughout the paper, we consider that Eq. A.11 is always valid and can be used to find the optimal combination of models that reproduces a given set of observations (see Sect. 4). The validity of this assumption is checked a posteriori by verifying the conditions written in Eqs. A.8 and A.9. The domain of parameters where these conditions break down are highlighted in Figs. 2 and D.3.

<footnote n="9">⁹ This approach can easily be expanded to multidimensional distributions. In this paper we assume a distribution of shock velocities, that is $\Theta_S = v_S$, and a distribution of illumination factors for the PDRs, that is $\Theta_P = G_0$.</footnote>
<figure><img src="image_10.png" /></figure>

**Fig. A.1.** Schematic representation of ensembles of shocks and PDRs inside an observational beam with a solid angle $\Omega_{\text{obs}}$ (black and white circle). Shock or PDR surfaces are represented by the filled circles and filled squares, respectively, and are randomly placed within the beam. The color represents the systemic velocities of shock and PDR surfaces projected along the line of sight. Blue (resp. red) circles and squares correspond to blueshifted (resp. redshifted) shocks and PDRs with respect to the recession velocity of the galaxy. The *left* and *right* panels represent two equivalent scenarios, where the numbers of shock (resp. PDR) surfaces are different, but subtend the same total solid $\Omega_S$ (resp. $\Omega_P$). This figure is adapted from Lehmann & Godard (2022).

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

**Fig. A.2.** Opacities of the emission lines in shocks as functions of the pre-shock density and the shock velocity (top row) and in PDRs as a function of the gas density and the illumination factor (bottom row). All other parameters are set to their standard values given in Table 1. The left column is associated with the H$_2$ 1-0 S(1) line, the middle column with the [CII] 158 µm line, and the right column with [OI] 63 µm line. The color-coding indicates the central opacity of the lines $\tau_S(v_0)$ and $\tau_P(v_0)$, with a different color scale for shocks and PDRs.
## Appendix B: The radio galaxy 3C 326 N

3C 326 is a well-known radio source located at $z \sim 0.089 \pm 0.001$ (Ogle et al. 2007). Optical and IR imaging shows that it is a pair of interacting galaxies, 3C 326 N and 3C 326 S (Fig. 1). We focus on 3C 326 N, centered at $\alpha(J2000) = 15^h52^m9^s.178$, $\delta(J2000) = +20^\mathrm{deg}05^\mathrm{m}48^\mathrm{s}.23$, which hosts the AGN powering the radio jet observed on large scales, with a kinetic luminosity $L_{\mathrm{jet}} \sim 10^{11} L_{\odot}$ (Nesvadba et al. 2010). This galaxy has a remarkably low star formation rate, SFR = $(7 \pm 3) \times 10^{-2} M_{\odot} \mathrm{yr}^{-1}$ (Ogle et al. 2007), yet a large molecular gas content $M_{\mathrm{gas}} \sim 2 \times 10^9 M_{\odot}$ (Guillard et al. 2015b). Strong IR H$_2$ line emission has been detected, with a total luminosity $L(H_2) \approx 2.4 \times 10^8 L_{\odot}$ (Ogle et al. 2007). 3C 326 N represents a perfect target to test which type of physical mechanisms (e.g., shocks and/or PDRs) are responsible for powering the ensemble of H$_2$ lines, while putting constraints on the energy budget of the galaxy. Based on UV GALEX measurements, Ogle et al. (2010) derived a value of $G_0 = 6 \times 10^{-4}$ for the radiation field in 3C 326 N. Similarly, Guillard et al. (2015b) reported a value of $G_0 = 9 \pm 1$ from a fit of the SED in the IR range.

Table B.1 provides a summary of the relevant physical parameters of 3C 326 N.<fnref n="10" /> The observational solid angle is estimated from the source diameter in optical and IRAC images (Nesvadba et al. 2010). The bulk of the H$_2$ emission is contained inside this beam (Ogle et al. 2010; Nesvadba et al. 2011) as well as the [CII] 158 µm and [OI] 63 µm emissions (Guillard et al. 2015b). The observed fluxes given in Table B.2 are corrected for redshift as

$$F^{\mathrm{corr}} = F^{\mathrm{obs}} (1 + z)^3 \quad (B.1)$$

and transformed into integrated intensities using the observational solid angle $\Omega_{\mathrm{obs}}$. We note that the total observed luminosity of the H$_2$, [CII] 158 µm, and [OI] 63 µm lines is $2.5 \times 10^8 L_{\odot}$, which corresponds to about 3% of the IR luminosity and 0.25% of the jet's luminosity.

### B.1. H$_2$ pure rotational and ro-vibrational lines

The observed fluxes and relevant information on the emission lines of molecular hydrogen are given in Table B.2. A total of 12 H$_2$ emission lines have been measured in 3C 326 N. Eight of these are the pure-rotational lines from H$_2$ 0-0 S(0) to H$_2$ 0-0 S(7), obtained with the mid-infrared *Spitzer* IRS low-resolution spectroscopy (R ~ 60 - 120), which did not resolve the lines (Ogle et al. 2010). Four ro-vibrational emission lines, from H$_2$ 1-0 S(1) to H$_2$ 1-0 S(4), were spatially and spectrally resolved by the VLT/SINFONI spectrograph (Nesvadba et al. 2011). However, since the main goal is to fit the observed integrated intensity of each line, the spectral resolution does not have any impact on our analysis.

### B.2. [CII] 158 µm, [OI] 63 µm, and CO(1-0)

The observed fluxes and relevant information on the emission lines of C$^+$, O, and CO are given in Table B.2. The [CII] 158 µm and the [OI] 63 µm emission lines are taken from observations obtained with *Herschel/PACS* and reported in Guillard et al. (2015b). The CO(1-0) line intensity is taken from observations performed with the IRAM PdBI and provided by Nesvadba et al.

<table><caption>Table B.1. Observational parameters of 3C 326 N.</caption><thead><tr><th>Parameters</th><th>Symbol</th><th>Value</th><th>Units</th></tr></thead><tbody><tr><td><sup>a</sup>Observational solid angle</td><td>Ω<sub>obs</sub></td><td>13.69</td><td>arcsec<sup>2</sup></td></tr><tr><td>Angular diameter distance</td><td>D<sub>A</sub></td><td>356</td><td>Mpc</td></tr><tr><td><sup>b</sup>Bolometric luminosity</td><td>L<sub>bol</sub></td><td>~ 2.6 × 10<sup>10</sup></td><td>L<sub>☉</sub></td></tr><tr><td><sup>c</sup>IR luminosity</td><td>L<sub>IR</sub></td><td>(7.1 ± 1.4) × 10<sup>9</sup></td><td>L<sub>☉</sub></td></tr><tr><td><sup>d</sup>Kinetic power available</td><td>L<sub>K</sub></td><td>~ 10<sup>9</sup></td><td>L<sub>☉</sub></td></tr><tr><td><sup>e</sup>Jet kinetic power</td><td>L<sub>jet</sub></td><td>~ 10<sup>11</sup></td><td>L<sub>☉</sub></td></tr><tr><td><sup>f</sup>Redshift</td><td>z</td><td>0.089</td><td>—</td></tr></tbody></table>

**Notes.** <sup>a</sup> Solid angle derived from the source diameter in optical and IRAC images (Nesvadba et al. 2010). <sup>b</sup> Bolometric luminosity reported in Nesvadba et al. (2010). <sup>c</sup> SED estimated value reported in Table C.3. <sup>d</sup> Observational derived values (Ogle et al. 2007; Guillard et al. 2015b). <sup>e</sup> Estimated value in Nesvadba et al. (2010). <sup>f</sup> Value reported in Ogle et al. (2007).

<table><caption>Table B.2. Emission line properties of 3C 326 N.</caption><thead><tr><th>Line</th><th>λ<sub>rest</sub><br>(μm)</th><th>Flux<br>(10<sup>-15</sup> erg s<sup>-1</sup> cm<sup>-2</sup>)</th><th>FWHM<br>(km s<sup>-1</sup>)</th><th>SNR</th></tr></thead><tbody><tr><td>0-0 S(0)</td><td>28.22</td><td>3.0 ± 0.6 <sup>a</sup></td><td>—</td><td>5.0</td></tr><tr><td>0-0 S(1)</td><td>17.03</td><td>6.9 ± 0.6 <sup>a</sup></td><td>—</td><td>11.5</td></tr><tr><td>0-0 S(2)</td><td>12.28</td><td>4.1 ± 0.4 <sup>a</sup></td><td>—</td><td>10.2</td></tr><tr><td>0-0 S(3)</td><td>9.66</td><td>12.6 ± 0.5 <sup>a</sup></td><td>—</td><td>25.2</td></tr><tr><td>0-0 S(4)</td><td>8.03</td><td>3.1 ± 0.6 <sup>a</sup></td><td>—</td><td>5.2</td></tr><tr><td>0-0 S(5)</td><td>6.91</td><td>4.6 ± 2.2 <sup>a</sup></td><td>—</td><td>2.1</td></tr><tr><td>0-0 S(6)</td><td>6.11</td><td>2.5 ± 0.9 <sup>a</sup></td><td>—</td><td>2.8</td></tr><tr><td>0-0 S(7)</td><td>5.55</td><td>2.9 ± 0.9 <sup>a</sup></td><td>—</td><td>3.2</td></tr><tr><td>1-0 S(1)</td><td>2.122</td><td>1.8 ± 0.1 <sup>b</sup></td><td>534 ± 32</td><td>13.8</td></tr><tr><td>1-0 S(2)</td><td>2.034</td><td>0.7 ± 0.1 <sup>b</sup></td><td>563 ± 45</td><td>14.0</td></tr><tr><td>1-0 S(3)</td><td>1.958</td><td>2.2 ± 0.2 <sup>b</sup></td><td>594 ± 41</td><td>11.0</td></tr><tr><td>1-0 S(4)</td><td>1.892</td><td>0.8 ± 0.1 <sup>b</sup></td><td>475 ± 45</td><td>16.0</td></tr><tr><td>[CII]</td><td>157.74</td><td>11.8 ± 1.0 <sup>c</sup></td><td>372 ± 12</td><td>11.8</td></tr><tr><td>[OI]</td><td>63.18</td><td>6.2 ± 1.4 <sup>c</sup></td><td>352 ± 31</td><td>4.7</td></tr><tr><td>CO(1-0)</td><td>2600</td><td>0.0035 ± 0.0007 <sup>d</sup></td><td>352 ± 58</td><td>5.0</td></tr></tbody></table>

**Notes.** <sup>a</sup> Ogle et al. (2010) from observations using *Spitzer* IRS. <sup>b</sup> Nesvadba et al. (2011) from observations using VLT/SINFONI. <sup>c</sup> Guillard et al. (2015b) from observations using *Herschel/PACS*. <sup>d</sup> Nesvadba et al. (2010) from observations using IRAM PdBI.

(2010). As explained in the main text, this CO(1-0) emission is not used to constrain the parameters of the model, but is used as a sanity check to test the validity of the solutions obtained.

## Appendix C: Spectral energy distribution

The SED of 3C 326 N is constrained using the photometric data spanning from the X-ray to the radio bands and available on the NASA/IPAC Extragalactic Database (NED). In this paper we use the measurements of 14 photometric bands obtained, in particular, by Lilly & Longair (1984), Ogle et al. (2010), and Guillard et al. (2015b), and given in Table C.1. These bands include the near-ultraviolet (NUV) GALEX band; the Sloan Digital Sky Surveys u, g, r, i, and z bands (Almeida et al. 2023); and the UKIRT J, UKIRT H, and UKIRT Ks NIR density fluxes. We finally include the *Spitzer* MIPS 24 µm and MIPS 70 µm bands, the *Herschel* PACS 100 µm and 160 µm bands, and the *Herschel* SPIRE 250 µm band.

<footnote n="10">¹⁰ With the adopted cosmology, the 3C 326 redshift of $z = 0.089$ corresponds to a luminosity distance $D_L = 422$ Mpc, an angular diameter distance $D_A = 356$ Mpc, and a scale of 1.725 kpc/arcsec.</footnote>
<table><caption>**Table C.1.** Photometry of 3C 326 N. The photometric bands derived from the Sloan Digital Sky Surveys (SDSS) were obtained using a circular aperture with a diameter of 8 arcsec on the SDSS DR18 (Almeida et al. 2023).</caption><thead><tr><th>Bands</th><th>Flux (mJy)</th></tr></thead><tbody><tr><td>NUV GALEX*</td><td>0.011 ± 0.002 <sup>a</sup></td></tr><tr><td>SDSS u</td><td>0.081 ± 0.017 <sup>b</sup></td></tr><tr><td>SDSS r</td><td>0.430 ± 0.040 <sup>b</sup></td></tr><tr><td>SDSS g</td><td>1.082 ± 0.063 <sup>b</sup></td></tr><tr><td>SDSS i</td><td>1.749 ± 0.080 <sup>b</sup></td></tr><tr><td>SDSS z</td><td>2.371 ± 0.093 <sup>b</sup></td></tr><tr><td>UKIRT J</td><td>2.36 ± 0.06 <sup>c</sup></td></tr><tr><td>UKIRT H</td><td>2.66 ± 0.07 <sup>c</sup></td></tr><tr><td>UKIRT K<sub>s</sub></td><td>2.39 ± 0.06 <sup>c</sup></td></tr><tr><td>Spitzer MIPS 24 µm</td><td>0.52 ± 0.07 <sup>d</sup></td></tr><tr><td>Spitzer MIPS 70 µm</td><td>6.12 ± 0.65 <sup>d</sup></td></tr><tr><td>Herschel PACS 100 µm</td><td>16.7 ± 0.6 <sup>d</sup></td></tr><tr><td>Herschel PACS 160 µm</td><td>15.9 ± 0.6 <sup>d</sup></td></tr><tr><td>Herschel SPIRE 250 µm</td><td>9.0 ± 3.5 <sup>d</sup></td></tr></tbody></table>

Notes. * Ogle et al. (2010) report a NUV AB magnitude of 21.47 (i.e., surface brightness of $1.59 \times 10^{-6}$ W m$^{-2}$ sr$^{-1}$) assuming a diameter area of 5 arcsec. The value is converted into flux density with an uncertainty equal to 20% of the observation. (a) Ogle et al. (2010); (b) This work; (c) Lilly & Longair (1984); (d) Guillard et al. (2015b)

### C.1. SED Fitting

The multiwavelength SED of 3C 326 N covering the UV to FIR range is fitted using the Code Investigating GALaxy Emission (CIGALE<fnref n="11" />; Boquien et al. 2019; Yang et al. 2020, 2022). CIGALE is based on an energy balance principle where ultraviolet-optical emission is partially absorbed by the dust and subsequently re-emitted in the IR. The code is based on a Bayesian statistical analysis approach that compares millions of models of the panchromatic emission of galaxies with the observations and provides estimates of the mean values and the standard deviations of the physical parameters of the galaxy.

A delayed star formation history without a burst is adopted to model this galaxy based on its low SFR with the simple stellar population models of Bruzual & Charlot (2003) and the initial mass function of Chabrier (2003). We adopt the attenuation law recipe of Charlot & Fall (2000) and the dust IR spectral templates from Draine et al. (2014). In contrast with the results presented in Lanz et al. (2016), we include the IR contribution of the AGN using the SKIRTOR clumpy torus models of

<figure><img src="image_12.png" /><figcaption>**Fig. C.1.** Spectral energy distribution of 3C 326 N. The MIR-to-FIR photometry is fitted using CIGALE. The stellar attenuated and unattenuated contributions are shown in blue and orange, respectively. The reprocessed IR emission is presented in purple and the AGN emission is in green. The best model is shown as a solid black line. The observed uncertainties are added in quadrature with a 10% floor value to account for possible instrumental errors.</figcaption></figure>

Stalevski et al. (2016) and following the configuration reported in Yang et al. (2022) with a fixed opening angle of 40° and varying type-2 viewing angles ($\theta$). All the parameters used for the fit are summarized in Table C.2.

### C.2. Physical parameters: 3C 326 N

Figure C.1 shows the results of the fit of the SED. The black solid line corresponds to the best SED and the different colored lines show the contribution of the stellar, dust, and AGN emission to the overall SED of the galaxy. We obtain a reduced $\chi^2$ of 1.2. The Bayesian estimated values of the stellar mass and SFR are presented in Table C.3 with new estimations of the UV and IR luminosities integrated between 0.0912 and 0.24 µm and between 1 and 1000 µm, respectively. The visual attenuation $A_V$ and the relative contribution of the AGN to the IR emission $f_{AGN}$ are also reported. The derived stellar mass and SFR are in good agreement (within 30%) with the values reported in the literature (Ogle et al. 2010; Nesvadba et al. 2010; Guillard et al. 2015b). The IR luminosity is also consistent with the fits performed by Guillard et al. (2015b).

<footnote n="11">¹¹ The code is available at: https://cigale.lam.fr/</footnote>
<table><caption>**Table C.2.** CIGALE modules and input parameters</caption><thead><tr><th>Parameter</th><th>Symbol</th><th>Value</th></tr></thead><tbody><tr><td><strong>Delayed Star Formation History</strong></td><td></td><td></td></tr><tr><td>e-folding time of main stellar population (Myr)</td><td>τ<sub>main</sub></td><td>250, 500, 1000, 2000, 4000, 6000, 8000, 10000, 12000, 15000</td></tr><tr><td>age of main stellar population (Myr)</td><td>age<sub>main</sub></td><td>250, 500, 1000, 2000, 4000, 8000, 10000, 12000</td></tr><tr><td><strong>Stellar Populations</strong>; Bruzual &amp; Charlot (2003)</td><td></td><td></td></tr><tr><td>initial mass function</td><td>IMF</td><td>Chabrier (2003)</td></tr><tr><td>metallicity</td><td>Z</td><td>0.008</td></tr><tr><td><strong>Dust attenuation</strong>; Values adopted by Charlot &amp; Fall (2000)</td><td></td><td></td></tr><tr><td>ISM attenuation in the V-band (mag)</td><td>A<sup>ISM</sup><sub>V</sub></td><td>0.05, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6</td></tr><tr><td>A<sup>ISM</sup><sub>V</sub>/(A<sup>ISM</sup><sub>V</sub>+A<sup>BC</sup><sub>V</sub>)</td><td>μ</td><td>0.44</td></tr><tr><td><strong>Dust Emission</strong>; Draine et al. (2014)</td><td></td><td></td></tr><tr><td>mass fraction of PAH</td><td>qPAH</td><td>1.12, 2.5, 3.19, 4.58, 5.95</td></tr><tr><td>minimum radiation field</td><td>U<sub>min</sub></td><td>5.0, 10.0, 25.0</td></tr><tr><td>power-law slope dU/dM ∝ U<sup>α</sup></td><td>α</td><td>2.0</td></tr><tr><td>fraction illuminated from U<sub>min</sub> to U<sub>max</sub></td><td>γ</td><td>0, 0.02</td></tr><tr><td><strong>AGN (UV-to-IR)</strong> SKIRTOR; Stalevski et al. (2016)</td><td></td><td></td></tr><tr><td>AGN contribution to IR luminosity</td><td>f<sub>AGN</sub></td><td>0.05, 0.1, 0.2, 0.3</td></tr><tr><td>Viewing angle (deg)</td><td>θ</td><td>60, 70, 80, 90</td></tr><tr><td>Polar-dust color excess (mag)</td><td>E(B − V)<sub>PD</sub></td><td>0, 0.2, 0.4</td></tr></tbody></table>

<table><caption>**Table C.3.** Physical parameters derived from the fit of the SED.</caption><thead><tr><th>Parameter</th><th>Value</th></tr></thead><tbody><tr><td>SFR (M<sub>☉</sub> yr<sup>-1</sup>)</td><td>(5 ± 3) × 10<sup>-2</sup></td></tr><tr><td>M<sub>*</sub> (M<sub>☉</sub>)</td><td>(0.7 ± 0.1) × 10<sup>11</sup></td></tr><tr><td>L<sub>IR</sub> (L<sub>☉</sub>)</td><td>(7.1 ± 1.4) × 10<sup>9</sup></td></tr><tr><td>L<sub>UV</sub> (L<sub>☉</sub>)</td><td>(8.5 ± 0.4) × 10<sup>8</sup></td></tr><tr><td>A<sub>V</sub> (mag)</td><td>0.1</td></tr><tr><td>f<sub>AGN</sub></td><td>0.1</td></tr></tbody></table>
# Appendix D: Alternative scenarios

## D.1. PDR-only scenario

We test the possibility of reproducing the ensemble of observations in 3C 326 N using only PDR models. The impact of the PDR density on the interpretation of observations is shown in panel (a) of Fig. D.1, which displays the minimum distance between the observations and the distribution of PDR models for a density ranging from 10 to $10^4$ cm$^{-3}$. Panel (b) shows the UV-reprocessed luminosity as a function of the density and panel (c) displays the comparisons between the observations and the intensities of the lines predicted by the model at $n_H = 10^3$ cm$^{-3}$ (which corresponds to the smallest distance). Figure D.1 shows that PDR-only models cannot explain the observational dataset, and raise two main issues. First, the minimum value of the distance is large (~ 2.3). As shown in panel (c), the need to reproduce the observed intensities of all the ro-vibrational lines of H$_2$ leads to overestimations of the intensities of the [CII] 158 µm and [OI] 63 µm lines by one to two orders of magnitude. Second, the predicted UV-reprocessed luminosity is found to be more than one order of magnitude higher than that deduced from the IR luminosity of the galaxy (see panel b). This result demonstrates the importance of estimating the radiative energy budget of galaxies for the interpretation of atomic and molecular emission lines.

## D.2. Shock-only scenario

We also consider the possibility of reproducing all the observations using only shock models. The impact of the density of the pre-shock gas on the interpretation of observations is shown in Fig. D.2 (panel a) which displays the minimum distance between the observations and the distribution of shocks models for a density ranging from 10 to $10^4$ cm$^{-3}$. Panel (b) shows the mechanical luminosity as a function of the density and panel (c) displays the comparisons between the observations and the intensities of the lines predicted by the model at $n_H = 10^3$ cm$^{-3}$ (which corresponds to the smallest distance). Interestingly, the predicted mechanical luminosity is found to be in agreement with the available energy deduced from the jet's power as long as the density is below $10^3$ cm$^{-3}$ (see panel b). However, Fig. D.1 shows that a shock-only distribution cannot explain the entire observational dataset. As shown in panel (c), the need to reproduce the observed intensities of the ro-vibrational lines of H$_2$ and of the fine structure line of C$^+$ leads to an underestimation of the [OI] emission by about one order of magnitude. This result demonstrates the importance of including both

atomic and molecular lines for the interpretation of extragalactic observations.

## D.3. Gaussian distribution of the shock velocities

We finally explore the impact of the form of the probability distribution function $f_S$ assuming that the shock velocities follow a Gaussian distribution

$$f_S = \frac{1}{\sqrt{2\pi}\sigma_{v_s}}\exp\left(-\frac{(v_s - \mu_{v_s})^2}{2\sigma_{v_s}^2}\right) \qquad (D.1)$$

rather than an exponential distribution. Figure D.3 shows the distance computed with Eqs. 3 and 4 as a function of the parameters of the Gaussian distribution of shocks ($\Omega_S$, $\mu_{v_s}$, and $\sigma_{v_s}$) and the parameter of the Dirac distribution of PDRs ($\Omega_P$).

The most striking feature of Fig. D.3 is the strong degeneracy between $\Omega_S$, $\mu_{v_s}$, and $\sigma_{v_s}$. A narrow distribution of shock velocities ($\sigma_{v_s} \sim 6$ km s$^{-1}$) centered at $\mu_{v_s} = 12$ km s$^{-1}$ leads to the same distance as a broad distribution ($\sigma_{v_s} \sim 15$ km s$^{-1}$) centered at $\mu_{v_s} = 3$ km s$^{-1}$. This strong degeneracy is not surprising and is simply due to the fact that the intensities of the ro-vibrational lines of H$_2$ are mostly emitted by shocks with velocities between 6 and 15 km s$^{-1}$, regardless of the number of shocks below 6 km s$^{-1}$ (see Fig. 5).

The impact of the gas density on the interpretation of observations is shown in Fig. D.4 (panel a) which displays the minimum distance between the observations and the distribution of shocks and PDRs for a density ranging from 10 to $10^4$ cm$^{-3}$. Panel (b) shows the UV-reprocessed luminosity and the mechanical luminosity as functions of the density and panel (c) displays the comparisons between the observations and the intensities of the lines predicted by the model at $n_H = 10^3$ cm$^{-3}$ (which corresponds to the smallest distance).

This ensemble of figures is remarkably similar to those obtained for an exponential distribution of shock velocities (see Fig. 3). In particular, they indicate that the emission likely originates from the diffuse interstellar gas ($n_H \le 100$ cm$^{-3}$) and that the 0-0 S(0) line of H$_2$ and the fine structure lines of C$^+$ and O are mostly emitted by PDRs, while all the other H$_2$ lines are emitted by shocks, with radiative and mechanical energy budgets in agreement with observational constraints. All these results confirm the interpretation obtained in the main text with an exponential distribution of shock velocities. Figure D.3 reveals that assuming a more complex distribution only leads to additional degeneracies. This means that the current ensemble of observations is insufficient to discriminate between different distributions of shock velocities.
<figure><img src="image_13.png" /></figure>

**Fig. D.1.** Summary of the minimum distance, the UV-reprocessed luminosity, and the intensities predicted with a sole distribution of PDR models for the PDR-only scenario. **Panel (a)**: Minimum distance as a function of the gas density. The dashed line represents a 0.5 dex limit (i.e., a factor of three difference between observed and modeled line intensities) to guide the eye. **Panel (b)**: UV-reprocessed luminosity as a function of the gas density. The orange shaded area represents the SED estimation of the reprocessed UV-luminosity (see Appendix C). **Panel (c)**: Comparison between the observed line intensities and those obtained for the model at $n_H = 10^3$ cm$^{-3}$, $b = 0.1$, $G_0 = 10$. The observations are presented as blue stars (see Table B.2) and the modeled intensities are shown as open black circles. In the lower panel, the residuals are shown as orange squares. The dashed and solid lines give the ±0.5, 1 dex limits, respectively.
<figure><img src="image_14.png" /></figure>

**Fig. D.2.** Summary of the minimum distance, the mechanical luminosity, and the intensities predicted with a sole distribution of shock models for the shock-only scenario. **Panel (a)**: Minimum distance as a function of the pre-shock density. The dashed line represents a 0.5 dex limit (i.e., a factor of three difference between observed and modeled line intensities) to guide the eye. **Panel (b)**: Mechanical luminosity as a function of the pre-shock density. The green area represents < 1% of the estimated jet kinetic luminosity of 3C 326 N (see Sect. 4.3). **Panel (c)**: Comparison between the observed line intensities and those obtained for the model at $n_H = 10 \text{ cm}^{-3}$, $b = 0.1$, $G_0 = 10$. Observations are presented as blue stars (see Table B.2) and modeled intensities are shown as open black circles.
<figure><img src="image_15.png" /></figure>

**Fig. D.3.** 2D cutoff of the likelihood distributions of the parameters describing the shock and the PDR distributions obtained for a Gaussian distribution of shock velocities. The parameters are the total shock solid angle $\Omega_S$, the central shock velocity $\mu_{v_s}$ and the dispersion of the shock velocities $\sigma_{v_s}$, and the PDR total solid angle $\Omega_P$. The color code represents the distance in dex units. The 0.5, 0.7, and 1.0 dex limits are highlighted with dashed black lines. The straight black dashed lines and the black points indicate the position of the global minimum. The gray areas represent the regions where the assumption on the radiative transfer regarding cross-absorption between surfaces breaks down. The criterion is more stringent for the [OI] 63 µm line (gray area above the white solid line) than for the [CII] 158 µm line (gray area above the white dashed line), and is always fulfilled for the H$_2$ lines (see Appendix A).
<figure><img src="image_16.png" /></figure>

**Fig. D.4.** Summary of the minimum distance, the UV-reprocessed and mechanical luminosities, and the intensities predicted for the PDR and shock scenario with a Gaussian distribution of shock velocities. **Panel (a)**: Minimum distance as a function of the gas density. The dashed line represents a 0.5 dex limit (i.e., a factor of three difference between observed and modeled line intensities) to guide the eye. **Panel (b)**: UV-reprocessed luminosity (orange) and mechanical luminosity (green) as functions of the gas density. The orange shaded area represents the SED estimation of the reprocessed UV-luminosity (see Appendix C). The green area represents < 1% of the estimated jet kinetic luminosity of 3C 326 N (see Sect. 4.3). **Panel (c)**: Comparison between the observed line intensities and those obtained for the model at $n_H = 10 \text{ cm}^{-3}$, $b = 0.1$, $G_0 = 10$. The observations are presented as blue stars (see Table B.2) and the modeled intensities are shown as open black circles.