11institutetext: Laboratoire de Physique de l’École Normale Supérieure, ENS,
Université PSL, CNRS, Sorbonne Université, Université Paris Cité, F-75005
Paris, France 22institutetext: Observatoire de Paris, PSL University, Sorbonne
Université, LERMA, 75014 Paris, France 33institutetext: Sorbonne Université,
CNRS, UMR 7095, Institut d’Astrophysique de Paris, 98bis bd Arago, 75014
Paris, France 44institutetext: Université Paris-Saclay, CNRS, Institut
d’Astrophysique Spatiale, 91405 Orsay, France

# Radiative and mechanical energies in galaxies

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

J. A. Villa-Vélez , E-mail<EMAIL_ADDRESS>B. Godard 2211 P.
Guillard 33 G. Pineau des Forêts 4422

(Received ; accepted )

###### 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 (PDR) 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 and 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 H2 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 ($\rm n_{H}<100\leavevmode\nobreak\
cm^{-3}$), indicating that most of the emission originates from diffuse
interstellar matter. The optimal solution, obtained for $\rm
n_{H}=10\leavevmode\nobreak\ cm^{-3}$, corresponds to a distribution of low-
velocity shocks (between $5$ and $20$ km s-1) propagating in PDR environments
illuminated by a UV radiation field $10$ 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 H2 0-0 S(0) $28\leavevmode\nobreak\ \mu$m, [CII]
$158\leavevmode\nobreak\ \mu$m, and [OI] $63\leavevmode\nobreak\ \mu$m lines
originate from the PDR components, while all the other H2 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
$\mathcal{L}_{\rm UV}=6.3\times 10^{9}\leavevmode\nobreak\ {\rm L_{\odot}}$,
and $\mathcal{L}_{\rm K}=3.9\times 10^{8}\leavevmode\nobreak\ {\rm
L_{\odot}}$, respectively, in remarkable agreement with the values of the
infrared luminosity deduced from the fit of the spectral energy distribution
of 3C 326 N, and consistent with a small fraction of the AGN jet kinetic power
dissipated in the interstellar medium ($\approx 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 show a striking agreement and demonstrate
the ability of the model to make accurate predictions. This framework
therefore 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: photon-dominated region
(PDR) – ISM: molecules – 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 multi-phase 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 non-thermal emissions across the electromagnetic spectrum,
from radio waves to gamma-rays. The dissipation of radiative and mechanical
energy in the ISM is intricately 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, metal and dust abundance, and the strengths of
heating/ionization sources like the ultra-violet (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 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 & Pineau
des Forêts, 2003; Godard et al., 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. Illustratively, the Spectral
Line Energy Distribution (SLED) of the carbon monoxide (CO) molecule is
frequently modeled through a mix of Photo-Dissociation Regions (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 H2 lines.
Given that H2 can form on interstellar dust at elevated temperatures
($\approx$ 50 to 250 K, Grieco et al., 2023) and that the H2 rotational levels
are separated by a few hundreds of 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 & Pineau des
Forêts, 1999; Kristensen et al., 2023). A previous study by Lesaffre et al.
(2013) used probability distribution functions (PDFs) of shock velocities to
model the rotational and ro-vibrational lines of H2 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 H2 line emissions exceeds what can be accounted for by the
exclusive reprocessing of the available UV radiative luminosity (Guillard et
al., 2009; Ogle et al., 2010; Nesvadba et al., 2010; Herrera et al., 2012;
Guillard et al., 2015b, a).

This paper presents an extension of the framework presented in Lehmann et al.
(2022), where we now combine distributions of shocks and PDRs models. The
state-of-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 well as the methodology 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, and finally, we provide predictions of the suite of ro-
vibrational H2 emission lines that will be observed with the JWST. Throughout
this work, we assume a modern flat cosmology (Planck Collaboration et al.,
2020) with $\Omega_{\rm m}$ = 0.315, $\Omega_{\Lambda}$ = 0.685, $\Omega_{\rm
b}$ = 0.0493, and H0 = 67.4 km s-1 Mpc-1.

##  2 Shock and PDR models

Throughout this paper, two publicly available codes are used to model PDR and
shock emissions: the Meudon PDR code111Meudon PDR code:
https://ism.obspm.fr/pdr.html (Le Petit et al., 2006) and the Paris-Durham
shock code222Paris-Durham code: https://ism.obspm.fr/shock.html (Flower &
Pineau des Forêts, 2003; Godard et al., 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, $\rm
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 $\rm G_{0}$. The size
of the cloud is set by its total extinction $\rm A_{v\leavevmode\nobreak\
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
rovibrational levels of H2. In particular, the excitation of the fine
structure levels of C+ and O is computed including excitation by collisions
with H, H2, and electrons (for C+) and with H, H+, H2, He, and electrons (for
O). The populations of the ro-vibrational levels of H2 are computed taking
into account excitation by collisions with H, He, H2, and H+, the probability
of exciting H2 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 $\rm
A_{v\leavevmode\nobreak\ max}$.

Figure 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,\leavevmode\nobreak\ {\rm
and}\leavevmode\nobreak\ 8.0\leavevmode\nobreak\ \mu{\rm m}$ showing the 3C
326 N galaxy (and its South companion). The $3.7^{\prime\prime}\times
3.7^{\prime\prime}$ ($6.4\times 6.4$ kpc) square (not-to-scale) represents the
solid angle $\Omega_{\rm 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
($\rm n_{H}\geq 10^{3}\leavevmode\nobreak\ cm^{-3}$), the shock diffuse gas in
yellow (i.e., $\rm n_{H}\sim 10-10^{3}\leavevmode\nobreak\ cm^{-3}$), and the
hot X-ray plasma in green ($\rm n_{H}\sim 1\leavevmode\nobreak\ cm^{-3}$).
Right: simplified model used in this paper where the observational beam is
supposed to encompass distributions of shock and PDR surfaces (see Fig. 9).

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 (ISMDB333Available at:
http://ism.obspm.fr/index.html). We consider PDR models at constant proton
density only, with a left-side illumination scaling factor $\rm G_{0}=1$ or 10
and a right-side illumination scaling factor $\rm G_{0}=1$. The range of
parameters explored here is shown in Table 1.

Table 1: PDR and shock model parameters.

PDR parameter | Symbol | Value | Units
---|---|---|---
Proton density | $\rm n_{H}$ | 10,$10^{2}$,$10^{3}$,$10^{4}$ | $\rm cm^{-3}$
Radiation field | $\rm G_{0}$ | $1,10$ | —
Total visual extinction | $\rm A{{}_{\rm v\leavevmode\nobreak\ max}}$ | $10$ | $\rm mag$
H2 CR ionization rate | $\zeta_{\rm H_{2}}$ | $10^{-16}$ | $\rm s^{-1}$
Shock parameter | Symbol | Value | Units
Proton density | $\rm n_{H}$ | 10,$10^{2}$,$10^{3}$,$10^{4}$ | $\rm cm^{-3}$
Radiation field(a) | $\rm G_{0}$ | $0,1,10$ | —
Initial visual extinction | $\rm A{{}_{\rm v}}$ | $1.0$ | $\rm mag$
Shock velocity | $\rm v_{s}$ |  $1-40$ in steps of 1 | $\rm km\leavevmode\nobreak\ s^{-1}$
Magnetic parameter(b) | $\rm b$ | $0.1,1.0$ | —
H2 CR ionization rate | $\zeta_{\rm H_{2}}$ | $10^{-16}$ | $\rm s^{-1}$

Notes.

(a) Scaling factor of the standard ultraviolet radiation field of Mathis et
al. (1983).

(b) Transverse magnetic field $\rm
B_{0}=b(n_{H}/cm^{-3})^{1/2}\leavevmode\nobreak\ \mu G$.

###  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
$\rm v_{s}$ in a magnetized and irradiated environment with an initial proton
density, $\rm n_{H}$, and an initial magnetic field strength set by a magnetic
parameter $\rm b$. The geometry, the main assumptions, and the physical
ingredients included in the code are described in Godard et al. (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 $\rm G_{0}$, and to be initially located at a position within the
PDR set by a visual extinction $\rm 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 H2. In particular, the
excitation of the fine structure levels of C+ and O and of the ro-vibrational
levels of H2 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 decide 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-1 in steps of 1 km
s-1.

##  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 schematized in Fig. 1, those
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 focus on molecular
lines and on typical tracers of the neutral interstellar medium, we neglect
here the contribution of the hot ionized phase. Under this approximation, we
develop 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_{\rm S}$ and $\Omega_{\rm P}$,
are captured within an observational beam of solid angle $\Omega_{\rm 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_{\rm S}$ and $\Omega_{\rm P}$ are assumed to follow
1D probability distribution functions denoted by $f_{\rm S}$ and $f_{\rm P}$.
The shock distribution is supposed to depend on the shock velocity $\rm v_{s}$
and the PDR distribution is supposed to depend on the illumination factor $\rm
G_{0}$, yielding the specific functional forms $f_{\rm S}({\rm
v_{s}})=d\Omega_{\rm S}/d{\rm v_{s}}$ and $f_{\rm P}({\rm G_{0}})=d\Omega_{\rm
P}/d{\rm G_{0}}$, respectively.

Considering that the continuum absorption due to the shock and PDR surfaces is
negligible and that no cross-absorption occurs between the
surfaces444Regardless of this assumption, the radiative transfer inside each
shock and PDR structure is computed taking into account absorption processes.
(see Eqs. 18 and 19), we derive a simple expression of the continuum
subtracted line flux (see Eq. 21). It follows that the integrated intensity of
each line (in $\rm erg\leavevmode\nobreak\ s^{-1}\leavevmode\nobreak\
cm^{-2}\leavevmode\nobreak\ sr^{-1}$) is given by

$I^{\rm m}=\frac{1}{\Omega_{\rm obs}}\left(\int_{\rm v_{s}}f_{\rm S}({\rm
v_{s}})I_{\rm s}({\rm v_{s}})d{\rm v_{s}}+\int_{\rm G_{0}}f_{\rm P}({\rm
G_{0}})I_{\rm p}({\rm G_{0}})d{\rm G_{0}}\right),$ (1)

where $I_{\rm S}(\rm v_{s})$ and $I_{\rm P}(\rm G_{0})$ are the line
integrated intensities emitted by one shock of velocity $\rm v_{s}$ and one
PDR of illumination factor $\rm G_{0}$, respectively. The superscript ‘m’
stands for model.

Note that the framework can be easily extended to other types of contributions
to the line emission, such as XDRs and HII-region for instance. 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 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.

Let’s consider an ensemble of N emission lines. According to Eq. 1, the
modeled intensity of each line $j$ writes

$I_{j}^{\rm m}=\frac{1}{\Omega_{\rm obs}}\left(\int_{\rm v_{s}}f_{\rm S}({\rm
v_{s}})I_{j,\rm s}({\rm v_{s}})d{\rm v_{s}}+\int_{\rm G_{0}}f_{\rm P}({\rm
G_{0}})I_{j,\rm p}({\rm G_{0}})d{\rm G_{0}}\right).$ (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

${\rm d}_{j}=\left\\{\begin{array}[]{@{}l@{\thinspace}l}0&\leavevmode\nobreak\
\leavevmode\nobreak\ \text{if}\leavevmode\nobreak\ \leavevmode\nobreak\ I^{\rm
o}_{j}-\frac{W_{j}}{2}\leq I^{\rm m}_{j}\leq I^{\rm o}_{j}+\frac{W_{j}}{2}\\\
\log_{10}(I^{\rm m}_{j})-\log_{10}\left(I^{\rm
o}_{j}-\frac{W_{j}}{2}\right)&\leavevmode\nobreak\ \leavevmode\nobreak\
\text{if}\leavevmode\nobreak\ \leavevmode\nobreak\ I^{\rm m}_{j}<I^{\rm
o}_{j}-\frac{W_{j}}{2}\\\ \\\ log_{10}(I^{\rm m}_{j})-\log_{10}\left(I^{\rm
o}_{j}+\frac{W_{j}}{2}\right)&\leavevmode\nobreak\ \leavevmode\nobreak\
\text{if}\leavevmode\nobreak\ \leavevmode\nobreak\ I^{\rm m}_{j}>I^{\rm
o}_{j}+\frac{W_{j}}{2},\\\ \end{array}\right.$ (3)

where $I^{\rm o}_{j}$ is the observed intensity of line $j$ and $W_{j}$ is the
associated one $\sigma$ 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

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

Eq. 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. Such 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 procedure:

  1. 1.

The 1D distribution functions $f_{\rm S}({\rm v_{s}})$ and $f_{\rm P}({\rm
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_{\rm S}({\rm
v_{s}})=(\Omega_{\rm S}/\sigma_{\rm v_{s}})\exp(-{\rm v_{s}}/\sigma_{\rm
v_{s}})$, where $\sigma_{\rm v_{s}}$ is the dispersion of shock velocities.
For simplicity, we consider a Dirac delta distribution of PDRs illumination
factor: $f_{\rm P}({\rm G_{0}})=\Omega_{\rm P}\delta(\rm G_{0})$.

  2. 2.

The modeled intensities are computed using Eq. 2 over a wide range of the
parameters $\Omega_{\rm S}$, $\sigma_{\rm v_{s}}$, and $\Omega_{\rm P}$
describing the distributions $f_{\rm S}$ and $f_{\rm P}$, assuming that all
the other parameters (e.g., $\rm n_{H}$, b, see Table 1) are constants. The
associated distances $\rm d(\Omega_{\rm S},\sigma_{\rm v_{s}},\Omega_{\rm P})$
between the model and the set of observations are computed with Eq. 4.

  3. 3.

Global and local minima of the distance are finally searched for to identify
the distributions $f_{\rm S}$ and $f_{\rm P}$ which fall the closest to the
observational dataset.

###  3.3 Mechanical an UV-reprocessed luminosities

The probability distribution functions $f_{\rm S}$ and $f_{\rm 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 $\rm v_{s}$ in a medium with pre-shock mass density
$\rho_{0}$ is

$\mathcal{L}_{\rm K}({\rm v_{s}})=\frac{1}{2}\rho_{0}{\rm
v_{s}^{3}}\mathcal{A},$ (5)

where $\mathcal{A}$ stands for the shock surface area (see Fig. 9). The total
mechanical luminosity dissipated (in $\rm erg\leavevmode\nobreak\
s^{-1}\leavevmode\nobreak\ sr$) by a distribution of shocks at different
velocities therefore writes

$\mathcal{L}_{\rm K}=\frac{1}{2}\rho_{0}\left(\int_{\rm v_{s}}f_{\rm S}({\rm
v_{s}}){\rm v^{3}_{s}}d{\rm v_{s}}\right){\rm D^{2}_{A}},$ (6)

where ${\rm 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 $\rm G_{0}$ on the front side. The
integrated flux of the Mathis et al. (1983) ISRF (in $\rm
erg\leavevmode\nobreak\ cm^{-2}\leavevmode\nobreak\ s^{-1}$) from 912 to 2400
Å (Le Petit et al., 2006) implies that the UV-reprocessed power by one PDR is

$\mathcal{L}_{\rm UV}({\rm G_{0}})=1.92\times 10^{-3}\leavevmode\nobreak\
({\rm G_{0}+1.0})\mathcal{A}.$ (7)

The total reprocessed UV-luminosity for an ensemble of PDRs that follow a 1D
probability distribution $f_{\rm P}({\rm G_{0}})$, therefore writes

$\mathcal{L}_{\rm UV}=1.92\times 10^{-3}\leavevmode\nobreak\ ({\rm
G_{0}+1.0})\left(\int_{\rm G_{0}}f_{\rm P}({\rm G_{0}})d{\rm G_{0}}\right){\rm
D^{2}_{A}}.$ (8)

###  3.4 Shock and PDR masses

The probability distribution functions $f_{\rm S}$ and $f_{\rm P}$ also
directly provide the total masses of shocks and PDRs within the observational
beam. The total mass of shocked material writes

${\rm M}_{\rm S}=1.4{\rm m_{H}}\left(\int_{\rm v_{s}}f_{\rm S}({\rm
v_{s}}){\rm N_{H}(v_{s})}d{\rm v_{s}}\right){\rm D^{2}_{A}},$ (9)

where $\rm m_{H}$ is the hydrogen mass and $\rm N_{H}(v_{s})$ is the proton
column density across a shock of velocity $\rm v_{s}$. The total mass of PDR
material depends on the typical PDR size or equivalently their typical visual
extinction $\rm A_{v}$. This mass writes

${\rm M}_{\rm P}({\rm A_{v}})=1.4{\rm m_{H}}\left(\int_{\rm G_{0}}f_{\rm
P}({\rm G_{0}})d{\rm G_{0}}\right){\rm N_{H}(A_{v})}{\rm D^{2}_{A}},$ (10)

where $\rm N_{H}(A_{v})$ is the proton column density evaluated at a given
visual extinction $\rm 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 (Ogle et al., 2010, $\sim
0.07\leavevmode\nobreak\ {\rm M_{\odot}\leavevmode\nobreak\ yr^{-1}}$), a
substantial reservoir of molecular content (Nesvadba et al., 2011, $\sim
2\times 10^{9}\leavevmode\nobreak\ \rm{M_{\odot}}$), and the conspicuous
presence of H2 emission lines (see Appendix B). The low star formation rate 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 which include H2 pure-
rotational and ro-vibrational transitions, alongside [CII]
$158\leavevmode\nobreak\ \mu$m and [OI] $63\leavevmode\nobreak\ \mu$m emission
lines. The observed fluxes and uncertainties used for the distance
minimization (see Eq. 3) are listed in Table 3.

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 beam555We 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 H2-to X-ray luminosity ratio
($>0.6$) is so high that X-ray heating is not a primary driver of H2 line
emission.. For simplicity, we consider that all PDRs and shocks have a unique
and common density and pre-shock density, $\rm n_{H}$, and are illuminated by
a unique and common UV radiation field characterized by a scaling factor $\rm
G_{0}$. As explained in Sect. 2, we also consider that shocks propagate in
PDRs at a typical visual extinction666Note that this choice of typical visual
extinction has little impact on the results presented here because the
predictions of the shock models regarding H2 emission are insensitive to $\rm
A_{v}$ as long as the shocks run in a medium mostly molecular (Kristensen et
al., 2023). $\rm A_{v}=1.0$. Finally, we assume that all shocks propagate in a
medium with a constant transverse magnetic field strength $\rm B_{0}$, scaled
by the magnetization parameter b that follows the relation $\rm
B_{0}=b(n_{H}/cm^{-3})^{1/2}\leavevmode\nobreak\ \mu G$.

As a fiducial model, we adopt $\rm n_{H}=10\leavevmode\nobreak\ cm^{-3}$, $\rm
b=0.1$, and $\rm 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 $\rm G_{0}=6^{+4}_{-2}$
and $\rm 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 $\rm G_{0}$ and $\rm b$ have a minor impact on the distribution of shocks
and PDRs required to explain the observations as long as $1\leq\rm G_{0}\leq
10$ and $\rm b\leq 1$. We will therefore limit our discussion to the effect of
the gas density only.

###  4.1 Likelihood analysis of the fiducial model

Figure 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_{\rm S}$, the
dispersion of the shock velocities $\sigma_{\rm v_{s}}$, and the PDR total
solid angle $\Omega_{\rm 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. Gray areas represent the regions where the
assumption on the radiative transfer regarding cross absorption between
surfaces breaks down. Note that the criterion is more stringent for the [OI]
$63\leavevmode\nobreak\ \mu$m line (gray area above the white solid line) than
for the [CII] $158\leavevmode\nobreak\ \mu$m line (gray area above the white
dashed line) and is always fulfilled for H2 lines (see Appendix A).

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_{\rm
S}$, $\sigma_{\rm v_{s}}$, and $\Omega_{\rm 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 1D,
plane-parallel geometry of the models. The global minimum is obtained for
distributions of shocks and PDRs with total solid angles $\Omega_{\rm
S}=1.9\times 10^{-8}\leavevmode\nobreak\ {\rm sr}$ and $\Omega_{\rm
P}=1.2\times 10^{-9}\leavevmode\nobreak\ {\rm sr}$, and a shock velocity
distribution characterized by $\sigma_{\rm v_{s}}\approx 4\leavevmode\nobreak\
{\rm km/s}$. This result implies that there are about $10$ 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_{\rm v_{s}}\approx
4\leavevmode\nobreak\ {\rm km/s}$ implies that the peak of dissipation of
mechanical energy occurs for $\rm v_{s}=3\sigma_{\rm v_{s}}\sim
12\leavevmode\nobreak\ {\rm km/s}$ and that most of the dissipation is
mediated by shocks with velocities between $\sim 5$ and $\sim
20\leavevmode\nobreak\ {\rm 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_{\rm S}$ or $\Omega_{\rm 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_{\rm S}$ and $\sigma_{\rm v_{s}}$: solutions
with distances smaller than $0.5$ dex can be obtained for $\Omega_{\rm S}$
varying over one order of magnitude and $\sigma_{\rm v_{s}}$ over a factor of
$2$. This degeneracy comes from the fact that the emission lines of H2 mostly
originate from shocks with a specific range of velocities (see Sect. 5) whose
numbers are set by both $\Omega_{\rm S}$ and $\sigma_{\rm v_{s}}$.

Fig. 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. 18 and 19 are satisfied. The grey 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 is, however, not valid for all values of $\Omega_{\rm S}$ and
$\Omega_{\rm P}$. The most stringent constraints on the validity of the
radiative transfer are due to the $63\leavevmode\nobreak\ \mu$m line of O and,
to a lesser extent, to the $158\leavevmode\nobreak\ \mu$m line of C+. This
feature simply comes from the high opacities of these two lines across shock
and PDR surfaces. In practice, the grey 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 treatment would be
superfluous. Indeed, 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

Figure 3: Solution for a combination of shocks and PDRs, assuming an
exponential PDF. Panel (a): distance versus the medium’s density. The dash-
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): 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
colors, 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).

The impact of the density of the pre-shock gas and of the PDRs on the
interpretation of observations is shown in Fig. 3 (a) which displays the
minimum distance between the observations and the distribution of shocks and
PDRs models for a density ranging from $10$ to $\rm 10^{4}\leavevmode\nobreak\
cm^{-3}$. The smallest value of the minimum distance is obtained for $\rm
n_{H}=10\leavevmode\nobreak\ 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 shown 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 H2 lines from 0-0 S(1) to
1-0 S(4) are mostly emitted by shocks. Conversely, the 0-0 S(1) line of H2,
the $158\leavevmode\nobreak\ \mu$m line of C+, and the $63\leavevmode\nobreak\
\mu$m line of O are mostly emitted by PDRs.

Fig. 3 (a) shows that the minimum distance rapidly increases with the gas
density. For instance, a minimum distance larger than $\sim 1.4$ dex is
obtained for $\rm n_{H}=10^{3}\leavevmode\nobreak\ cm^{-3}$, meaning that at
most one of the observational lines is underestimated or overestimated by a
factor $\sim 25$. This result is in contradiction with those currently
presented in the literature. Fits of the excitation diagram of the H2 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
$\rm 10^{4}\leavevmode\nobreak\ cm^{-3}$. These discrepancies come from the
fact that our observational dataset contains additional constraints, in
particular the [CII] $158\leavevmode\nobreak\ \mu$m, and [OI]
$63\leavevmode\nobreak\ \mu$m lines. As shown in Fig. 13, pure shock models
cannot reproduce the [OI] line emission, even at low density. Moreover, in
Nesvadba et al. (2010), a density of $\rm 10^{4}\leavevmode\nobreak\ cm^{-3}$
was needed to reproduce high-excitation lines, e.g. 0-0 S(7), but
consequently, the 0-0 S(0) line was over-predicted by a factor of $\approx
2.5$. The fact that the [CII] $158\leavevmode\nobreak\ \mu$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 to compute the mechanical energy dissipated by
shocks, $\mathcal{L}_{\rm K}$, as well as the UV-reprocessed luminosity in
PDRs, $\mathcal{L}_{\rm UV}$. The resulting luminosities are shown in Fig. 3
(b) as functions of the pre-shock and PDR densities and compared with the
values estimated for 3C 326 N. For the fiducial model ($\rm
n_{H}=10\leavevmode\nobreak\ {\rm cm}^{-3}$, $\rm b=0.1$, and $\rm G_{0}=10$),
we derive a dissipated mechanical power of $\mathcal{L}_{\rm K}=3.9\times
10^{8}\leavevmode\nobreak\ {\rm L_{\odot}}$ and a UV-reprocessed luminosity of
$\mathcal{L}_{\rm UV}=6.3\times 10^{9}\leavevmode\nobreak\ {\rm L_{\odot}}$.
The mechanical energy dissipated by the shocks is therefore as high as $\sim
6$% of the UV-reprocessed luminosity.

Figure 4: Solution for a combination of shocks and PDRs, assuming an
exponential PDF. Panel (c): solution for the fiducial model $\rm
n_{H}=10\leavevmode\nobreak\ cm^{-3}$, $\rm b=0.1$, $\rm G_{0}=10$ shown in
Fig. 3 (a). Observations are presented as blue stars (see Table 3) and modeled
intensities are shown as open black circles. Individual shock and PDR
contributions are shown as green and orange solid lines, respectively. In the
lower panel, the residuals are shown as black squares. The dash and solid
lines give the $\rm\pm 0.5,1\leavevmode\nobreak\ dex$ limits, respectively.

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 infrared wavelengths.
To obtain an independent estimate of the amount of UV radiation reprocessed,
we perform 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 CIGALE777Available
at: https://cigale.lam.fr/ (Boquien et al., 2019; Yang et al., 2020, 2022).
This treatment leads to a reprocessed UV luminosity $\rm L_{UV}=(7.1\pm
1.4)\times 10^{9}\leavevmode\nobreak\ L_{\odot}$. This estimate is in
agreement with Lanz et al. (2016), and is shown as an orange shaded region in
Fig. 3 (b). 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$ to $\rm
10^{3}\leavevmode\nobreak\ {cm}^{-3}$. Interestingly, at $\rm n_{H}\gtrsim
10^{4}\leavevmode\nobreak\ 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 kpc-sized
cavities 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
of the order of several hundreds of $\rm km\leavevmode\nobreak\ s^{-1}$ 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 sub-kpc 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 $\rm\sim 10^{11}\leavevmode\nobreak\ L_{\odot}$ from
radio observations. The upper limit computed as 1% of this luminosity is shown
as a green shaded area in Fig. 3 (b). 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 infrared
radiation, the associated visual extinction must be larger than 1. We assume
two values of the visual extinction of PDRs, $\rm A_{v}=1$ and $\rm A_{v}=10$.
Using Eqs. 9 and 10, we derive a total mass of shocks, ${\rm M}_{\rm
S}=7.8\times 10^{7}\leavevmode\nobreak\ {\rm M}_{\odot}$, and a total mass of
PDR, ${\rm M}_{\rm P}=1.9\times 10^{8}\leavevmode\nobreak\ {\rm M}_{\odot}$
(for $\rm A_{v}=1$) and ${\rm M}_{\rm P}=1.9\times 10^{9}\leavevmode\nobreak\
{\rm M}_{\odot}$ (for $\rm 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 H2 of $\rm 1.1\times 10^{9}\leavevmode\nobreak\ M_{\odot}$ and
$\rm 1.3-2.7\times 10^{9}\leavevmode\nobreak\ M_{\odot}$, respectively, from
the excitation diagram of H2 observed in 3C 326 N. Similarly, Guillard et al.
(2015b) estimate a total molecular mass of $\sim\rm 2\times
10^{9}\leavevmode\nobreak\ M_{\odot}$ from the modeling of the FIR emission of
dust grains and the emission of the CO(1-0) line.

Figure 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, $\rm v_{s\leavevmode\nobreak\ dom}$, that contributes the most
to powering each line is highlighted in orange. We display only a range
between 1-20 km/s for readability, even though the models include velocities
up to 40 km/s.

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 ($\rm
n_{H}=10\leavevmode\nobreak\ cm^{-3}$), we find that the volume of PDRs is
comprised between $\rm 0.55\leavevmode\nobreak\ kpc^{3}$ (for $\rm A_{v}=1$)
and $\rm 5.5\leavevmode\nobreak\ kpc^{3}$ (for $\rm A_{v}=10$). The volume of
the galaxy is uncertain. Near-infrared IFU spectroscopy with the VLT/SINFONI
instrument show that the H2 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 $\rm 50\leavevmode\nobreak\
kpc^{3}$. These estimations show that the models with densities $\rm
n_{H}\geqslant 10\leavevmode\nobreak\ cm^{-3}$ are coherent because they lead
to a volume filling factor smaller than $1$. Given the results shown in Fig. 3
one could wonder why we do not explore gas at densities lower than $\rm
10\leavevmode\nobreak\ 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 2 is not included
as a constraint in our model. Low-J CO emission are known to arise from the
cold diffuse component of molecular gas (Mingozzi et al., 2018). Given that
the total mass of shocks is much smaller 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}$ erg s-1 cm-2, a value in fair agreement,
yet slightly larger than that observed in 3C 326 N ($3.5\times 10^{-18}$ erg
s-1 cm-2). This result reinforces our interpretation of this galaxy in which
the 0-0 S(0), [CII] $158\leavevmode\nobreak\ \mu$m, and [OI]
$63\leavevmode\nobreak\ \mu$m are mostly emitted by PDRs, while the remaining
H2 lines are mostly emitted by shocks.

##  5 Model predictions for JWST observations

###  5.1 Impact of shock velocities to the H2 lines

From the fitting of observed line emission in 3C 326 N (Fig. 3), it is clear
that H2 lines are mainly powered by shocks, except the lowest excitation 0-0
S(0) line where the PDRs contribution dominate. For the fiducial model ($\rm
n_{H}=10\leavevmode\nobreak\ {\rm cm}^{-3}$, $\rm b=0.1$, and $\rm G_{0}=10$,
see Sect. 4), we now investigate, amongst the distribution of shocks, which
shock velocities are the dominant contributors to the emission of H2 lines.

Fig. 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 H2 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 range $\rm
v_{s,dom}=6\leavevmode\nobreak\ km/s$ to $\rm 15\leavevmode\nobreak\ km/s$.
The typical length scales over which those shocks dissipate energy is of 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 H2 line
emission (Nesvadba et al., 2010).

On a fairly more technical side, the shocks responsible for the H2 line
emission are J-type (J for Jump) shocks (see Kristensen et al., 2023, for a
discussion of the shock model grid, and in particular their Fig. 6). Those
shocks are strong emitters of ro-vibrational emission lines, contrary to
C-type (C for 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.

Figure 6: Shock velocity that contributes the most to the H2 line emissions as
a function of the upper-level energy (orange points in Fig. 5). H2 pure-
rotational (green) and ro-vibrational (red) transitions are presented. $\rm
v_{s,dom}$ is derived for the fiducial model $\rm n_{H}=10\leavevmode\nobreak\
cm^{-3}$, $\rm b=0.1$, and $\rm G_{0}=10$.

###  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 a Gaussian distribution for $f_{\rm S}$. Appendix
D.3, and in particular Figs. 14 and 15, 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. In fact, Fig. 14 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_{\rm v_{s}}=12$ km s-1 leads to the same distance as a broad distribution
centered at $\mu_{\rm v_{s}}=3$ km s-1. This implies that the current ensemble
of observations is insufficient to discriminate different distributions of
shock velocities as long as the shocks that dominate the rovibrational
emission of H2 (between 6 and 15 km s-1, see Fig. 5) are included.

###  5.3 H2 line predictions for JWST

Figure 7: Predicted synthetic spectrum of H2 in 3C 326 N obtained for the
fiducial model ($\rm n_{H}=10\leavevmode\nobreak\ {cm}^{-3}$, $\rm b=0.1$, and
$\rm G_{0}=10$). The specific intensities of the strongest $1000$ emission
lines are shown as functions of the redshifted wavelength for $\rm z=0.089$.
We apply the relations described in Sect. 5.3 for the resolving power of the
NIRSpec and MIRI instruments. Outside the validity wavelength range, we fix
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 narrow-band filters for NIRCam and MIRI are shown as gray and color bars.
Pure rotation lines are indicated in blue while ro-vibrational transitions are
shown in red and green.  Figure 8: Comparison of the intensities of the H2
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)^{3}$ (see Eq.
22). Observations are shown as blue stars and the model predictions as opened
black circles. We stress that the model is obtained by fitting the Spitzer and
VLT/SINFONI data, and not the JWST data.

The Paris-Durham shock and the Meudon PDR codes solve for the populations of a
large number of H2 levels (see Sect. 2) and therefore provide a predicted
spectrum over the entire wavelength range observed by the JWST instruments.
Fig. 7 displays the spectrum predicted by the fiducial model for the strongest
1000 lines of H2, accounting for the different spectral resolutions of the
NIRSpec and MRS spectrographs. We fit the resolving power of the G235H NIRSpec
disperser888We use wavelength calibration data from https://jwst-
docs.stsci.edu/jwst-near-infrared-spectrograph/nirspec-
instrumentation/nirspec-dispersers-and-filters over the wavelength range
$1.7-3.1\leavevmode\nobreak\ \mu{\rm 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\leavevmode\nobreak\ \mu{\rm m}$. Outside the validity
range in wavelength of these relations, we use a uniform resolving power of
$\lambda/\Delta\lambda=2500$. 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 $\mu$m.

We predict that the 3C 326 N spectrum should be dominated by ro-vibrational
lines with $\upsilon_{\rm up}=1$ around $2-3\leavevmode\nobreak\ \mu{\rm m}$,
as well as the strong 0-0 S(3) and 0-0 S(5) mid-infrared lines. The forest of
near-infrared $\upsilon_{\rm up}\geq 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 H2 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. Fig. 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 to interpret 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 (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 H2 line emission. A total of 14
emission lines including eight H2 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\leavevmode\nobreak\ \mu$m and [OI] $63\leavevmode\nobreak\ \mu$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 $\rm G_{0}$ should
spread over a wide range of values. Indeed, this diversity of ISM environments
is obtained in numerical simulations (e.g., Mukherjee et al. 2018; Murthy et
al. 2022) and 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 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_{\rm S}$ and $\sigma_{\rm v_{s}}$). This implies
that adopting a distribution of density or $\rm G_{0}$, while more realistic,
would only lead to additional degeneracies. The key results obtained within
this framework are summarized as follows.

  1. 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.

  2. 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 ($\rm n_{H}=10-100\leavevmode\nobreak\ cm^{-3}$) 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. 3.

The optimal solution is obtained for $\rm n_{H}=10\leavevmode\nobreak\
cm^{-3}$. The likelihood, estimated from the ratios between the observed and
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. 4.

The H2 0-0 S(0) $28\leavevmode\nobreak\ \mu$m, [CII] $158\leavevmode\nobreak\
\mu$m, and [OI] $63\leavevmode\nobreak\ \mu$m lines are mostly emitted by
PDRs. The resulting UV-reprocessed luminosity, $\mathcal{L}_{\rm UV}=6.3\times
10^{9}\leavevmode\nobreak\ {\rm L_{\odot}}$, is found to be in remarkable
agreement with the observed IR-luminosity.

  5. 5.

Conversely, the rest of the H2 lines from the 0-0 S(1) to the 1-0 S(4) are
found to be emitted by low-velocity shocks ($\rm 5<v_{s}<20$ km s-1). The
resulting mechanical luminosity dissipated by shocks, $\mathcal{L}_{\rm
K}=3.9\times 10^{8}\leavevmode\nobreak\ {\rm 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. 6.

The total mass of shocks is ${\rm M}_{\rm S}=7.8\times
10^{7}\leavevmode\nobreak\ {\rm M}_{\odot}$ and the total mass of PDRs is
${\rm M}_{\rm P}=1.9\times 10^{9}\leavevmode\nobreak\ {\rm 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. 7.

Finally, the predicted synthetic H2 spectrum shows prominent ro-vibrational
emission lines, originating from all the vibrational levels of H2. 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., $\rm
n_{H}\sim 10^{4}\leavevmode\nobreak\ cm^{-3}$) could explain the CO, H2, and
[CII] $158\leavevmode\nobreak\ \mu$m emission lines observed towards 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 H2 and fine-structure emission lines in extragalactic sources, in
particular in AGN (Pereira-Santaella et al., 2022; Buiten et al., 2023).

###### Acknowledgements.

We are grateful to the referee for their thorough reading and their insightful
comments. 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 et al. (2008) Allen, M. G., Groves, B. A., Dopita, M. A., Sutherland, R. S., & Kewley, L. J. 2008, ApJS, 178, 20
  * Almeida et al. (2023) Almeida, A., Anderson, S. F., Argudo-Fernández, M., et al. 2023, ApJS, 267, 44
  * Argyriou et al. (2023) Argyriou, I., Glasse, A., Law, D. R., et al. 2023, A&A, 675, A111
  * Begelman & Ruszkowski (2005) Begelman, M. C. & Ruszkowski, M. 2005, Philosophical Transactions of the Royal Society of London Series A, 363, 655
  * Benson (2010) Benson, A. J. 2010, Phys. Rep, 495, 33
  * Boquien et al. (2019) Boquien, M., Burgarella, D., Roehlly, Y., et al. 2019, A&A, 622, A103
  * Bruzual & Charlot (2003) Bruzual, G. & Charlot, S. 2003, MNRAS, 344, 1000
  * Buiten et al. (2023) Buiten, V. A., van der Werf, P. P., Viti, S., et al. 2023, arXiv e-prints, arXiv:2312.01945
  * Carnall et al. (2018) Carnall, A. C., McLure, R. J., Dunlop, J. S., & Davé, R. 2018, MNRAS, 480, 4379
  * Casey et al. (2014) Casey, C. M., Narayanan, D., & Cooray, A. 2014, Phys. Rep, 541, 45
  * Chabrier (2003) Chabrier, G. 2003, PASP, 115, 763
  * Charlot & Fall (2000) Charlot, S. & Fall, S. M. 2000, ApJ, 539, 718
  * Cielo et al. (2018) Cielo, S., Bieri, R., Volonteri, M., Wagner, A. Y., & Dubois, Y. 2018, MNRAS, 477, 1336
  * Ciotti et al. (2010) Ciotti, L., Ostriker, J. P., & Proga, D. 2010, ApJ, 717, 708
  * da Cunha et al. (2008) da Cunha, E., Charlot, S., & Elbaz, D. 2008, MNRAS, 388, 1595
  * Draine et al. (2014) Draine, B. T., Aniano, G., Krause, O., et al. 2014, ApJ, 780, 172
  * Esposito et al. (2022) Esposito, F., Vallini, L., Pozzi, F., et al. 2022, MNRAS, 512, 686
  * Ferland et al. (2017) Ferland, G. J., Chatzikos, M., Guzmán, F., et al. 2017, Rev. Mexicana Astron. Astrofis., 53, 385
  * Fierlinger et al. (2016) Fierlinger, K. M., Burkert, A., Ntormousi, E., et al. 2016, MNRAS, 456, 710
  * Flower & Pineau des Forêts (1999) Flower, D. R. & Pineau des Forêts, G. 1999, MNRAS, 308, 271
  * Flower & Pineau des Forêts (2003) Flower, D. R. & Pineau des Forêts, G. 2003, MNRAS, 343, 390
  * Godard et al. (2019) Godard, B., Pineau des Forêts, G., Lesaffre, P., et al. 2019, A&A, 622, A100
  * Grieco et al. (2023) Grieco, F., Theulé, P., De Looze, I., & Dulieu, F. 2023, Nature Astronomy, 7, 541
  * Guillard et al. (2015a) 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 Astronomy and Astrophysics, 81–85
  * Guillard et al. (2015b) Guillard, P., Boulanger, F., Lehnert, M. D., et al. 2015b, A&A, 574, A32
  * Guillard et al. (2009) Guillard, P., Boulanger, F., Pineau des Forêts, G., & Appleton, P. N. 2009, A&A, 502, 515
  * Guillard et al. (2012) Guillard, P., Ogle, P. M., Emonts, B. H. C., et al. 2012, ApJ, 747, 95
  * Habart et al. (2005) Habart, E., Walmsley, M., Verstraete, L., et al. 2005, Space Sci. Rev., 119, 71
  * Hardcastle & Croston (2020) Hardcastle, M. J. & Croston, J. H. 2020, New A Rev., 88, 101539
  * Herrera et al. (2012) Herrera, C. N., Boulanger, F., Nesvadba, N. P. H., & Falgarone, E. 2012, A&A, 538, L9
  * Kim et al. (2022) Kim, J.-G., Gong, M., Kim, C.-G., & Ostriker, E. C. 2022, The Astrophysical Journal Supplement Series, 264, 10
  * Klessen & Hennebelle (2010) Klessen, R. S. & Hennebelle, P. 2010, A&A, 520, A17
  * Krause (2023) Krause, M. G. H. 2023, Galaxies, 11, 29
  * Kristensen et al. (2023) Kristensen, L. E., Godard, B., Guillard, P., Gusdorf, A., & Pineau des Forêts, G. 2023, A&A, 675, A86
  * Krumholz et al. (2014) Krumholz, M. R., Bate, M. R., Arce, H. G., et al. 2014, in Protostars and Planets VI, ed. H. Beuther, R. S. Klessen, C. P. Dullemond, & T. Henning, 243–266
  * Lanz et al. (2016) Lanz, L., Ogle, P. M., Alatalo, K., & Appleton, P. N. 2016, ApJ, 826, 29
  * Le Petit et al. (2006) Le Petit, F., Nehmé, C., Le Bourlot, J., & Roueff, E. 2006, ApJS, 164, 506
  * Leftley et al. (2024) Leftley, J. H., Nesvadba, N. P. H., Bicknell, G., et al. 2024, arXiv e-prints, arXiv:2404.04341
  * Lehmann et al. (2016) Lehmann, A., Federrath, C., & Wardle, M. 2016, MNRAS, 463, 1026
  * Lehmann et al. (2022) Lehmann, A., Godard, B., Pineau des Forêts, G., Vidal-García, A., & Falgarone, E. 2022, A&A, 658, A165
  * Lesaffre et al. (2013) Lesaffre, P., Pineau des Forêts, G., Godard, B., et al. 2013, A&A, 550, A106
  * Lesaffre et al. (2020) Lesaffre, P., Todorov, P., Levrier, F., et al. 2020, MNRAS, 495, 816
  * Lilly & Longair (1984) Lilly, S. J. & Longair, M. S. 1984, MNRAS, 211, 833
  * Mathis et al. (1983) Mathis, J. S., Mezger, P. G., & Panagia, N. 1983, A&A, 128, 212
  * Meijerink et al. (2013) Meijerink, R., Kristensen, L. E., Weiß, A., et al. 2013, ApJ, 762, L16
  * Mingozzi et al. (2018) Mingozzi, M., Vallini, L., Pozzi, F., et al. 2018, MNRAS, 474, 3640
  * Morganti et al. (2023) Morganti, R., Murthy, S., Guillard, P., Oosterloo, T., & Garcia-Burillo, S. 2023, Galaxies, 11, 24
  * Mukherjee et al. (2016) Mukherjee, D., Bicknell, G. V., Sutherland, R., & Wagner, A. 2016, MNRAS, 461, 967
  * Mukherjee et al. (2018) Mukherjee, D., Bicknell, G. V., Wagner, A. Y., Sutherland, R. S., & Silk, J. 2018, MNRAS, 479, 5544
  * Murthy et al. (2022) Murthy, S., Morganti, R., Wagner, A. Y., et al. 2022, Nature Astronomy, 6, 488
  * Nesvadba et al. (2011) Nesvadba, N. P. H., Boulanger, F., Lehnert, M. D., Guillard, P., & Salome, P. 2011, A&A, 536, L5
  * Nesvadba et al. (2010) Nesvadba, N. P. H., Boulanger, F., Salomé, P., et al. 2010, A&A, 521, A65
  * Ogle et al. (2007) Ogle, P., Antonucci, R., Appleton, P. N., & Whysong, D. 2007, ApJ, 668, 699
  * Ogle et al. (2010) Ogle, P., Boulanger, F., Guillard, P., et al. 2010, ApJ, 724, 1193
  * Ogle et al. (2024) Ogle, P. M., López, I. E., Reynaldi, V., et al. 2024, ApJ, 962, 196
  * Ostriker & Kim (2022) Ostriker, E. C. & Kim, C.-G. 2022, ApJ, 936, 137
  * Park & Ryu (2019) Park, J. & Ryu, D. 2019, ApJ, 875, 2
  * Pathak et al. (2024) Pathak, D., Leroy, A. K., Thompson, T. A., et al. 2024, AJ, 167, 39
  * Pereira-Santaella et al. (2022) Pereira-Santaella, M., Álvarez-Márquez, J., García-Bernete, I., et al. 2022, A&A, 665, L11
  * Pfenniger (2010) Pfenniger, D. 2010, in Astronomical Society of the Pacific Conference Series, Vol. 421, Galaxies in Isolation: Exploring Nature Versus Nurture, ed. L. Verdes-Montenegro, A. Del Olmo, & J. Sulentic, 183
  * Planck Collaboration et al. (2020) Planck Collaboration, Aghanim, N., Akrami, Y., et al. 2020, A&A, 641, A6
  * Polles et al. (2021) Polles, F. L., Salomé, P., Guillard, P., et al. 2021, A&A, 651, A13
  * Richard et al. (2022) Richard, T., Lesaffre, P., Falgarone, E., & Lehmann, A. 2022, A&A, 664, A193
  * Shaw et al. (2009) Shaw, G., Ferland, G. J., Henney, W. J., et al. 2009, ApJ, 701, 677
  * Stalevski et al. (2016) Stalevski, M., Ricci, C., Ueda, Y., et al. 2016, MNRAS, 458, 2288
  * Vazquez-Semadeni (2009) Vazquez-Semadeni, E. 2009, arXiv e-prints, arXiv:0902.0820
  * Willis & Strom (1978) Willis, A. G. & Strom, R. G. 1978, A&A, 62, 375
  * Wittor & Gaspari (2020) Wittor, D. & Gaspari, M. 2020, MNRAS, 498, 4983
  * Yang et al. (2022) Yang, G., Boquien, M., Brandt, W. N., et al. 2022, ApJ, 927, 192
  * Yang et al. (2020) Yang, G., Boquien, M., Buat, V., et al. 2020, MNRAS, 491, 740

##  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. 9, we consider that an
observational beam of solid angle $\Omega_{\rm obs}$ contains shock and PDR
surfaces with a distribution of systemic velocities and with total solid
angles $\Omega_{\rm S}$ and $\Omega_{\rm P}$. For simplicity, we set that the
emission from one shock or one PDR in a line with a rest frequency $\nu_{0}$
follows a tophat profile over a width $\Delta\nu=\nu_{0}\Delta V/c$, with a
velocity width $\Delta V=5$ km s-1 dominated by turbulent motions (Lehmann et
al. 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_{\nu}=\frac{1}{\sqrt{2\pi}\sigma_{\nu}^{\rm obs}}{\rm
exp}\left(-\frac{(\nu-\nu_{0})^{2}}{2\sigma_{\nu}^{{\rm obs}\,2}}\right),$
(11)

where $\sigma_{\nu}^{\rm obs}$ is the observed frequency dispersion of the
line. Independently of this distribution, we also assume that $\Omega_{\rm S}$
and $\Omega_{\rm P}$ follow 1D probability distribution functions, $f_{\rm
S}(\Theta_{\rm S})$ and $f_{\rm P}(\Theta_{\rm P})$, which both depend on a
unique parameter, $\Theta_{\rm S}$ for the shocks and $\Theta_{\rm P}$ for the
PDRs999This approach can easily be expanded to multidimensional distributions.
In this paper, we assume a distribution of shock velocities, that is
$\Theta_{\rm S}={\rm v_{S}}$, and a distribution of illumination factors for
the PDRs, that is $\Theta_{\rm P}=G_{0}$.. The solid angles of shock and PDR
surfaces that emit within a frequency bin of size $\Delta\nu$ centered at a
frequency $\nu$ therefore write

$\Delta\Omega_{\rm S}(\nu)=\Omega_{\rm S}\phi_{\nu}\Delta\nu=\int_{\Theta_{\rm
S}}f_{\rm S}(\Theta_{\rm S})d\Theta_{\rm S}\phi_{\nu}\Delta\nu$ (12)

and

$\Delta\Omega_{\rm P}(\nu)=\Omega_{\rm P}\phi_{\nu}\Delta\nu=\int_{\Theta_{\rm
P}}f_{\rm P}(\Theta_{\rm P})d\Theta_{\rm P}\phi_{\nu}\Delta\nu.$ (13)

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_{\rm S}}{\Omega_{\rm obs}}(\nu)$ (14)

and

$\frac{\Delta\Omega_{\rm P}}{\Omega_{\rm obs}}(\nu).$ (15)

Their average total opacities are

$\overline{\tau}_{\rm S}(\nu)=\int_{\Theta_{\rm S}}f_{\rm S}(\Theta_{\rm
S})\tau_{\rm S}(\Theta_{\rm S})d\Theta_{\rm S}\phi_{\nu}\Delta\nu/\Omega_{\rm
obs}$ (16)

and

$\overline{\tau}_{\rm P}(\nu)=\int_{\Theta_{\rm P}}f_{\rm P}(\Theta_{\rm
P})\tau_{\rm P}(\Theta_{\rm P})d\Theta_{\rm P}\phi_{\nu}\Delta\nu/\Omega_{\rm
obs},$ (17)

where $\tau_{\rm S}(\Theta_{\rm S})$ and $\tau_{\rm P}(\Theta_{\rm P})$ are
the opacities across one shock layer of parameter $\Theta_{\rm S}$ and one PDR
layer of parameter $\Theta_{\rm P}$ respectively (see Fig. 10).

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

$\frac{\Delta\Omega_{\rm S}}{\Omega_{\rm obs}}(\nu_{0})\ll 1\quad{\rm
and}\quad\frac{\Delta\Omega_{\rm P}}{\Omega_{\rm obs}}(\nu_{0})\ll 1,$ (18)

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

$\overline{\tau}_{\rm S}(\nu_{0})\ll 1\quad{\rm and}\quad\overline{\tau}_{\rm
P}(\nu_{0})\ll 1,$ (19)

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_{\nu}=I_{\nu}^{c}\Omega_{\rm obs}+\int_{\Theta_{\rm S}}f_{\rm
S}(\Theta_{\rm S})\phi_{\nu}I_{\rm S}(\Theta_{\rm S})d\Theta_{\rm
S}+\int_{\Theta_{\rm P}}f_{\rm P}(\Theta_{\rm P})\phi_{\nu}I_{\rm
P}(\Theta_{\rm P})d\Theta_{\rm S}$ (20)

where $I_{\nu}^{c}$ is the specific intensity of the background continuum (in
erg s-1 cm-2 Hz-1 sr-1) and $I_{\rm S}(\Theta_{\rm S})$ and $I_{\rm
P}(\Theta_{\rm P})$ are the line integrated intensities emitted by one shock
of parameter $\Theta_{\rm S}$ and one PDR of parameter $\Theta_{\rm 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_{\rm S}}f_{\rm S}(\Theta_{\rm S})I_{\rm S}(\Theta_{\rm
S})d\Theta_{\rm S}+\int_{\Theta_{\rm P}}f_{\rm P}(\Theta_{\rm P})I_{\rm
P}(\Theta_{\rm P})d\Theta_{\rm P}.$ (21)

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. 21 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. 18 and 19. The domain of parameters where these
conditions break down are highlighted in Figs. 2 and 14.

Figure 9: Schematic representation of ensembles of shocks and PDRs inside an
observational beam with a solid angle $\Omega_{\rm 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_{\rm S}$ (resp. $\Omega_{\rm P}$). This
figure is adapted from Lehmann et al. (2022).
Figure 10: Opacities of the emission lines in shocks as functions of the pre-
shock density and the shock velocity (top row) and in PDRs as 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 H2 1-0 S(1) line, the middle column with the [CII]
$158\leavevmode\nobreak\ \mu$m line, and the right column to [OI]
$63\leavevmode\nobreak\ \mu$m line. The color code indicates the central
opacity of the lines $\tau_{\rm S}(\nu_{0})$ and $\tau_{\rm P}(\nu_{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 $\rm 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) = $\rm 15^{h}52^{m}9^{s}.178$, $\delta$(J2000) =
$\rm+20^{deg}05^{m}48^{s}.23$, which hosts the AGN powering the radio jet
observed on large scales, with a kinetic luminosity $\rm L_{jet}\sim
10^{11}\leavevmode\nobreak\ L_{\odot}$ (Nesvadba et al. 2010). This galaxy has
a remarkably low star-formation rate $\rm SFR=(7\pm 3)\times
10^{-2}\leavevmode\nobreak\ M_{\odot}\leavevmode\nobreak\ yr^{-1}$ (Ogle et
al. 2007), yet a large molecular gas content $\rm M_{gas}\sim 2\times 10^{9}$
M⊙ (Guillard et al. 2015b). Strong IR H2 line emission has been detected, with
a total luminosity ${\rm L(H_{2})}\approx 2.4\times 10^{8}$ L⊙ (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 H2 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^{+4}_{-2}$ 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 2 provides a summary of the relevant physical parameters of 3C 326
N101010With the adopted cosmology, the 3C 326 redshift of $z=0.089$
corresponds to a luminosity distance DL = 422 Mpc, an angular diameter
distance DA = 356 Mpc, and a scale of 1.725 kpc/arcsec.. The observational
solid angle is estimated from the source diameter in optical and IRAC images
(Nesvadba et al. 2010). The bulk of the H2 emission is contained inside this
beam (Ogle et al. 2010; Nesvadba et al. 2011) as well as the [CII]
$158\leavevmode\nobreak\ \mu$m and [OI] $63\leavevmode\nobreak\ \mu$m
emissions (Guillard et al. 2015b). The observed fluxes given in Table 3 are
corrected for redshift as

${F}^{\rm corr}={F}^{\rm obs}\,\,(1+\rm{z})^{3}$ (22)

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

###  B.1 H2 pure rotational and ro-vibrational lines

The observed fluxes and relevant information on the emission lines of
molecular hydrogen are given in Table 3. A total of twelve H2 emission lines
have been measured in 3C 326 N. Eight of these are the pure-rotational lines
from H2 0-0 S(0) to H2 0-0 S(7), obtained with the mid-infrared Spitzer IRS
low-resolution spectroscopy ($\rm R\sim 60-120$), which did not resolve the
lines (Ogle et al. 2010). Four ro-vibrational emission lines, from H2 1-0 S(1)
to H2 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\leavevmode\nobreak\ \mu{\rm m}$, [OI]
$63\leavevmode\nobreak\ \mu{\rm m}$ and CO(1-0)

The observed fluxes and relevant information on the emission lines of C+, O,
and CO are given in Table 3. The [CII] $158\leavevmode\nobreak\ \mu$m and the
[OI] $63\leavevmode\nobreak\ \mu$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. (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.

Table 2: Observational parameters of 3C 326 N.

Parameters | Symbol | Value | Units
---|---|---|---
aObservational solid angle | $\Omega_{\rm obs}$ | $13.69$ | arcsec2
Angular diameter distance | DA | $356$ | Mpc
bBolometric luminosity | Lbol | $\sim 2.6\times 10^{10}$ | L⊙
cIR luminosity | LIR | $(7.1\pm 1.4)\times 10^{9}$ | L⊙
dKinetic power available | LK | $\sim 10^{9}$ | L⊙
eJet kinetic power | Ljet | $\sim 10^{11}$ | L⊙
fRedshift | z | $0.089$ | —

Notes.

a Solid angle derived from the source diameter in optical and IRAC images
(Nesvadba et al. 2010).

b Bolometric luminosity reported in Nesvadba et al. (2010).

c SED estimated value reported in Table 6.

d Observational derived values (Ogle et al. 2007; Guillard et al. 2015b).

e Estimated value in Nesvadba et al. (2010).

f Value reported in Ogle et al. (2007).

Table 3: Emission line properties of 3C 326 N.

Line | $\lambda_{\rm rest}$ | Flux | FWHM | SNR
---|---|---|---|---
| ($\mu{\rm m}$) | ($\rm 10^{-15}\leavevmode\nobreak\ erg\leavevmode\nobreak\ s^{-1}\leavevmode\nobreak\ cm^{-2}$) | (km s-1) |
0-0 S(0) | 28.22 |  $3.0\pm 0.6$ a | — | 5.0
0-0 S(1) | 17.03 |  $6.9\pm 0.6$ a | — | 11.5
0-0 S(2) | 12.28 |  $4.1\pm 0.4$ a | — | 10.2
0-0 S(3) | 9.66 |  $12.6\pm 0.5$ a | — | 25.2
0-0 S(4) | 8.03 |  $3.1\pm 0.6$ a | — | 5.2
0-0 S(5) | 6.91 |  $4.6\pm 2.2$ a | — | 2.1
0-0 S(6) | 6.11 |  $2.5\pm 0.9$ a | — | 2.8
0-0 S(7) | 5.55 |  $2.9\pm 0.9$ a | — | 3.2
1-0 S(1) | 2.122 |  $1.8\pm 0.1$ b | $534\pm 32$ | 13.8
1-0 S(2) | 2.034 |  $0.7\pm 0.1$ b | $563\pm 45$ | 14.0
1-0 S(3) | 1.958 |  $2.2\pm 0.2$ b | $594\pm 41$ | 11.0
1-0 S(4) | 1.892 |  $0.8\pm 0.1$ b | $475\pm 45$ | 16.0
[CII] | 157.74 |  $11.8\pm 1.0$ c | $372\pm 12$ | 11.8
$[\rm OI]$ | 63.18 |  $6.2\pm 1.4$ c | $352\pm 31$ | 4.7
CO(1-0) | 2600 |  $0.0035\pm 0.0007$ d | $352\pm 58$ | 5.0

Notes.

a Ogle et al. (2010) from observations using Spitzer IRS.

b Nesvadba et al. (2011) from observations using VLT/SINFONI.

c Guillard et al. (2015b) from observations using Herschel/PACS.

d Nesvadba et al. (2010) from observations using IRAM PdBI.

##  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 4. These bands include the
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 near-
infrared density fluxes. We finally include the Spitzer MIPS 24 $\mu$m and
MIPS 70 $\mu$m bands, the Herschel PACS 100 $\mu$m and 160 $\mu$m bands, and
the Herschel SPIRE 250 $\mu$m band.

Table 4: 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).

Bands | Flux (mJy)
---|---
NUV GALEX⋆ |  $0.011\pm 0.002$ a
SDSS u |  $0.081\pm 0.017$ b
SDSS r |  $0.430\pm 0.040$ b
SDSS g |  $1.082\pm 0.063$ b
SDSS i |  $1.749\pm 0.080$ b
SDSS z |  $2.371\pm 0.093$ b
UKIRT J |  $2.36\pm 0.06$ c
UKIRT H |  $2.66\pm 0.07$ c
UKIRT Ks |  $2.39\pm 0.06$ c
Spitzer MIPS 24 $\mu$m |  $0.52\pm 0.07$ d
Spitzer MIPS 70 $\mu$m |  $6.12\pm 0.65$ d
Herschel PACS 100 $\mu$m |  $16.7\pm 0.6$ d
Herschel PACS 160 $\mu$m |  $15.9\pm 0.6$ d
Herschel SPIRE 250 $\mu$m |  $9.0\pm 3.5$ d

Notes.

⋆ Ogle et al. (2010) report a NUV AB magnitude of $21.47$ (i.e., surface
brightness of $\rm 1.59\times 10^{-6}\leavevmode\nobreak\
W\leavevmode\nobreak\ m^{-2}\leavevmode\nobreak\ 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.

References ${\rm(a)}$ Ogle et al. (2010); ${\rm(b)}$ This work; ${\rm(c)}$
Lilly & Longair (1984); ${\rm(d)}$ Guillard et al. (2015b)

###  C.1 SED Fitting

Table 5: CIGALE modules and input parameters

Parameter | Symbol | Value
---|---|---
Delayed Star Formation History |  |
e-folding time of main stellar population (Myr) | $\tau_{\mathrm{main}}$ | 250, 500, 1000, 2000, 4000, 6000, 8000,10000,12000,15000
age of main stellar population (Myr) | agemain | 250, 500, 1000, 2000, 4000, 8000, 10000, 12000
Stellar Populations; Bruzual & Charlot (2003) |  |
initial mass function | IMF | Chabrier (2003)
metallicity | Z | 0.008
Dust attenuation. Values adopted by Charlot & Fall (2000) |  |
ISM attenuation in the V-band (mag) | A${}_{\mathrm{V}}^{\mathrm{ISM}}$ | 0.05, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6
A${}^{\mathrm{ISM}}_{\mathrm{V}}$/(A${}^{\mathrm{ISM}}_{\mathrm{V}}+$A${}^{\mathrm{BC}}_{\mathrm{V}}$) | $\mu$ | 0.44
Dust Emission; Draine et al. (2014) |  |
mass fraction of PAH | qPAH | 1.12, 2.5, 3.19, 4.58, 5.95
minimum radiation field | Umin | 5.0, 10.0, 25.0
power-law slope $\mathrm{dU/dM}\propto\mathrm{U}^{\alpha}$ | $\alpha$ | 2.0
fraction illuminated from $\mathrm{U_{min}}$ to $\mathrm{U_{max}}$ | $\gamma$ | 0, 0.02
AGN (UV-to-IR) SKIRTOR; Stalevski et al. (2016) |  |
AGN contribution to IR luminosity | fAGN | 0.05, 0.1, 0.2, 0.3
Viewing angle (deg) | $\theta$ | 60, 70, 80, 90
Polar-dust color excess (mag) | $\rm E(B-V)_{PD}$ | 0, 0.2, 0.4

The multi-wavelength SED of 3C 326 N covering the UV-to-FIR range is fitted
using the Code Investigating GALaxy Emission (CIGALE111111Code is available
here: https://cigale.lam.fr/; 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
infrared. 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.

Figure 11: 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.

A delayed star formation history without a burst is adopted to model this
galaxy based on its low star-formation rate 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 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 5.

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

Fig. 11 shows the results of the fit of the SED. The black solid line
corresponds to the best SED and the different color 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 6 with new estimations of
the UV and IR luminosities integrated between $0.0912$ and
$0.24\leavevmode\nobreak\ \mu{\rm m}$ and between $1$ and
$1000\leavevmode\nobreak\ \mu{\rm m}$, respectively. The visual attenuation,
$\rm A_{v}$, and the relative contribution of the AGN to the IR emission, $\rm
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).

Table 6: Physical parameters derived from the fit of the SED. Parameter | Value
---|---
$\rm SFR\leavevmode\nobreak\ (M_{\odot}\leavevmode\nobreak\ yr^{-1})$ | $(5\pm 3)\times 10^{-2}$
$\rm M_{*}\leavevmode\nobreak\ (M_{\odot})$ | $(0.7\pm 0.1)\times 10^{11}$
$\rm L_{IR}\leavevmode\nobreak\ (L_{\odot})$ | $(7.1\pm 1.4)\times 10^{9}$
$\rm L_{UV}\leavevmode\nobreak\ (L_{\odot})$ | $(8.5\pm 0.4)\times 10^{8}$
$\rm A_{V}\leavevmode\nobreak\ (mag)$ | $0.1$
$\rm f_{AGN}$ | $0.1$

##  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 on panel (a) of Fig. 12 which displays the minimum
distance between the observations and the distribution of PDR models for a
density ranging from $10$ to $\rm 10^{4}\leavevmode\nobreak\ cm^{-3}$. The
panel (b) shows the UV-reprocessed luminosity as a function of the density and
the panel (c) displays the comparisons between the observations and the
intensities of the lines predicted by the model at $\rm
n_{H}=10^{3}\leavevmode\nobreak\ cm^{-3}$ (which corresponds to the smallest
distance). Fig. 12 shows that sole PDR models cannot explain the observational
dataset and raise two main issues. First, the minimum value of the distance is
large ($\sim 2.3$). As shown in panel (c), the need to reproduce the observed
intensities of all the rovibrational lines of H2 leads to overestimations of
the intensities of the [CII] $158\leavevmode\nobreak\ \mu$m and [OI]
$63\leavevmode\nobreak\ \mu$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 larger 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.

Figure 12: PDR only scenario: summary of the minimum distance, the UV-
reprocessed luminosity, and the intensities predicted with a sole distribution
of PDR models. Panel (a): minimum distance as a function of the gas density.
The dash-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 $\rm n_{H}=10^{3}\leavevmode\nobreak\
cm^{-3}$, $\rm b=0.1$, $\rm G_{0}=10$. Observations are presented as blue
stars (see Table 3) and modeled intensities are shown as open black circles.
In the lower panel, the residuals are shown as orange squares. The dash and
solid lines give the $\rm\pm 0.5,1\leavevmode\nobreak\ dex$ limits,
respectively.

###  D.2 Shock-only scenario

Figure 13: Shock only scenario: summary of the minimum distance, the
mechanical luminosity, and the intensities predicted with a sole distribution
of shock models. Panel (a): minimum distance as a function of the pre-shock
density. The dash-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 $\rm n_{H}=10\leavevmode\nobreak\
cm^{-3}$, $\rm b=0.1$, $\rm G_{0}=10$. Observations are presented as blue
stars (see Table 3) and modeled intensities are shown as open black circles.
In the lower panel, the residuals are shown as green squares. The dash and
solid lines give the $\rm\pm 0.5,1\leavevmode\nobreak\ dex$ limits,
respectively.

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. 13 (panel a) which displays
the minimum distance between the observations and the distribution of shocks
models for a density ranging from $10$ to $\rm 10^{4}\leavevmode\nobreak\
cm^{-3}$. The panel (b) shows the mechanical luminosity as a function of the
density and the panel (c) displays the comparisons between the observations
and the intensities of the lines predicted by the model at $\rm
n_{H}=10\leavevmode\nobreak\ 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 $\rm 10^{3}\leavevmode\nobreak\ cm^{-3}$ (see panel b).
Fig. 12 shows, however, that a sole distribution of shocks cannot explain the
entire observational dataset. As shown in panel (c), the need to reproduce the
observed intensities of the rovibrational lines of H2 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_{\rm S}$ assuming that the shock velocities follow a Gaussian
distribution

$f_{\rm S}=\frac{1}{\sqrt{2\pi}\sigma_{\rm v_{s}}}{\rm exp}\left(-\frac{(\rm
v_{s}-\mu_{\rm v_{s}})^{2}}{2\sigma_{\rm v_{s}}^{2}}\right)$ (23)

rather than an exponential distribution. Fig. 14 shows the distance computed
with Eqs. 3 and 4 as a function of the parameters of the Gaussian distribution
of shocks, $\rm\Omega_{S}$, $\mu_{\rm v_{s}}$, and $\sigma_{\rm v_{s}}$, and
the parameter of the Dirac distribution of PDRs, $\Omega_{\rm P}$.

The most striking feature of Fig. 14 is the strong degeneracy between
$\Omega_{S}$, $\mu_{\rm v_{s}}$, and $\sigma_{\rm v_{s}}$. A narrow
distribution of shock velocities ($\sigma_{\rm v_{s}}\sim 6$ km s-1) centered
at $\mu_{\rm v_{s}}=12$ km s-1 leads to the same distance than a broad
distribution ($\sigma_{\rm v_{s}}\sim 15$ km s-1) centered at $\mu_{\rm
v_{s}}=3$ km s-1. This strong degeneracy is not surprising and is simply due
to the fact that the intensities of the rovibrational lines of H2 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. 15 (panel a) which displays the minimum distance between the
observations and the distribution of shocks and PDRs for a density ranging
from $10$ to $\rm 10^{4}\leavevmode\nobreak\ cm^{-3}$. The panel (b) shows the
UV-reprocessed luminosity and the mechanical luminosity as functions of the
density and the panel (c) displays the comparisons between the observations
and the intensities of the lines predicted by the model at $\rm
n_{H}=10\leavevmode\nobreak\ 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
($\rm n_{H}\leqslant 100$ cm-3) and that the 0-0 S(0) line of H2 and the fine
structure lines of C+ and O are mostly emitted by PDRs while all the other H2
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. Fig. 14 reveals that assuming a more complex distribution
only leads to additional degeneracies. This means that the current ensemble of
observations is insufficient to discriminate different distributions of shock
velocities.

Figure 14: 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_{\rm S}$, the central shock velocity $\mu_{\rm v_{s}}$ and the
dispersion of the shock velocities $\sigma_{\rm v_{s}}$, and the PDR total
solid angle $\Omega_{\rm 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. Gray areas represent the regions where the
assumption on the radiative transfer regarding cross absorption between
surfaces breaks down. Note that the criterion is more stringent for the [OI]
$63\leavevmode\nobreak\ \mu$m line (gray area above the white solid line) than
for the [CII] $158\leavevmode\nobreak\ \mu$m line (gray area above the white
dashed line) and is always fulfilled for H2 lines (see Appendix A).

Figure 15:  PDR and shock scenario with a Gaussian distribution of shock
velocities: summary of the minimum distance, the UV-reprocessed and mechanical
luminosities, and the intensities predicted. Panel (a): minimum distance as a
function of the gas density. The dash-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 $\rm
n_{H}=10\leavevmode\nobreak\ cm^{-3}$, $\rm b=0.1$, $\rm G_{0}=10$.
Observations are presented as blue stars (see Table 3) and modeled intensities
are shown as open black circles. In the lower panel, the residuals are shown
as green squares. The dash and solid lines give the $\rm\pm
0.5,1\leavevmode\nobreak\ dex$ limits, respectively.