© The Authors 2025

1. Introduction

Emission-line diagnostics have been a very useful tool in investigating the different and complex mechanisms of gas excitation within galaxies. The most commonly used diagnostics for disentangling gas ionized by radiation of active galactic nuclei (AGN) from stellar-driven ionization by O and B stars are based on optical lines such as [O III]λ5007 Å, [O I]λ6300 Å, Hα, Hβ, [N II]λ6584 Å, and [S II]λλ 6716,6731 Å; they are known as the Baldwin, Phillips, and Terlevich (BPT) diagrams (Baldwin et al. 1981). Since their introduction to the astronomical community, BPT diagrams, have been enhanced by theoretical limits and empirical divisions to distinguish between AGN-dominated and star-forming sources, as well as between Seyfert and low-ionization nuclear emission-line region (LINER) galaxies (Kauffmann et al. 2003; Kewley et al. 2006; Schawinski et al. 2007).

Modeling the emission-line fluxes obtained from spectra is typically carried out using photoionization codes such as CLOUDY (Ferland et al. 2013, 2017) and MAPPINGS (Binette 1985; Sutherland & Dopita 1993; Dopita & Sutherland 1995), which provide detailed constraints on the physical properties of ionized gas. These codes incorporate various excitation mechanisms, including ultraviolet (UV) and X-ray radiation, and in the case of MAPPINGS, shocks are also included (e.g. Allen et al. 2008; Ferland et al. 2009; Feltre et al. 2016; Ferland et al. 2017; Feltre et al. 2023; Chatzikos et al. 2023; Zhu et al. 2023). Such models are used to interpret observed emission lines across BPT and other diagnostic diagrams. For example, photoionization models can reproduce line fluxes in individual objects, providing detailed physical constraints, as shown by Pérez-Montero & Díaz (2007) for Mrk 209 and by Oliveira et al. (2024) for nitrogen and oxygen abundances in LINER galaxies. The modeling of the specific position of the observations, is affected by the physical characteristics of the emitting sourc–such as stellar mass, star formation rate (SFR), metallicity, and dust-to-metal mass ratio (Hirschmann et al. 2017; Feltre et al. 2016, 2023), as well as density, shape and intensity of the radiation field (Feltre et al. 2016; Zhu et al. 2023), making the modeling process quite challenging, especially when examining large samples of galaxies. A common approach is to compare observational data with precomputed grids of models covering a range of physical parameters. This method forms the basis of widely used tools such as HII-CHI-MISTRY (Pérez-Montero 2014), IZI (Blanc et al. 2015), NEBULABAYES (Thomas et al. 2018), and HOMERUN (Marconi et al. 2024).

Moreover, cosmic rays (CRs)–highly energetic particles primarily generated by supernova remnants (SNRs) and black hole jets– represent another potential but less explored excitation source (Blasi 2013; Padovani et al. 2017; Veilleux et al. 2020; Wolfire et al. 2022; Kantzas et al. 2023). Their ability to induce excitation and ionization via secondary electrons deep within molecular clouds significantly impacts nebular emission lines (Spitzer & Tomasko 1968; McKee 1989; Padovani et al. 2018; Gabici 2022), and has similar effects to supersolar metallicity on the [N II]λ6584 Å BPT diagram (Koutsoumpou et al. 2025), creating a degeneracy that complicates the interpretation of metallicity and the ionization source.

Infrared emission-line ratio diagnostics have become essential for examining ionization mechanisms in dusty or obscured regions of galaxies, and have been used to disentangle photoionization due to AGN activity from young stars (see Spinoglio & Malkan 1992). Diagrams based on high- to low-ionization line ratios–such as [O IV]λ25.9 μ_m/[Ne II]λ12.8 μ_m and [Ne V]λ14.3 μ_m/[Ne II]λ12.8 μ_m– are typically used to estimate AGN activity, as they trace ionizing photons with energies beyond those produced by standard stellar populations (Genzel et al. 1998; Sturm et al. 2002; Groves et al. 2006). Yet notably, [Ne V] emission has also been observed in galaxies lacking clear AGN signatures (e.g. Thuan & Izotov 2005; Hatano et al. 2023, 2024; Mingozzi et al. 2025). Moreover, ratios such as [Ne III]/[Ne II], involving lower ionization lines serve as indicators of the ionizing spectrum hardness and have been used to characterize stellar populations in star-forming regions (Thornley et al. 2000; Pérez-Montero et al. 2024). With the development of new state-of-the-art facilities, these diagnostics can now be spatially resolved, even in fainter or more embedded sources, and extended with lines such as [Ne VI] and [Ar V] (Richardson et al. 2022; Feltre et al. 2023). These line ratios complement optical diagnostics by providing a dust-insensitive view of excitation conditions, revealing obscured or composite sources, and enabling a more complete interpretation of the ionization mechanisms in galaxies that could also be used to identify ionization due to CRs.

In this study, we extended our previous work Koutsoumpou et al. (2025), hereafter referred to as KFDS25, in order to examine whether CRs are capable of producing the observed emission in both the optical and mid-infrared (MIR) wavelength range. Specifically, we integrated MIR data from the Mid-Infrared Instrument (MIRI) aboard the James Webb Space Telescope (JWST), with optical data from the Multi Unit Spectroscopic Explorer (MUSE) on the Very Large Telescope (VLT) to explore the dual impact of photoionization and CR-induced ionization, in a single source, NGC 5728. We tested the effectiveness of the proposed diagnostics on NGC 5728, as it is a jetted and extensively studied AGN, with both optical and MIR observations available for comparison. Section 2 describes the main features and the observational data for NGC 5728. Section 3 outlines our modeling procedure and the steps taken to treat data from certain areas of interest within NGC 5728. Section 4 outlines the findings of the present analysis, and Section 5 explores possibilities for further research. We conclude the present study with the summary in Section 6.

2. Observational data

We selected NGC 5728, a Seyfert 2 galaxy, which has been observed thoroughly and exhibits clearly apparent jets (Durré & Mould 2018; Davies et al. 2024). Additionally, this source features a star formation ring with a significant number of young stars (Schommer et al. 1988; Shin et al. 2019; Shimizu et al. 2019). This renders it a perfect case to study the effects of both jet and supernova-related CRs in AGN and starburst-dominated regions. NGC 5728 proves to be an appropriate target, since it has publicly available VLT/MUSE observations in the European Southern Observatory (ESO) Science Archive Facility (see Section 2.1), as well as data from MIRI/JWST (see Section 2.2). Lying in close proximity, at a redshift of z = 0.00932¹, NGC 5728 provides a great opportunity to explore effect of CRs in different galactic environments through a variety of high-quality observational data.

2.1. MUSE data

VLT-MUSE raw data are handled by the ESO MUSE pipeline, which includes bias correction, dark subtraction, flat fielding, and background subtraction to improve image quality and reduce noise. Furthermore, the pipeline conducts flux and wavelength calibration for spectral accuracy, and glitch removal from CRs. Subsequently, spectral extraction yields the object’s spectra from spatially resolved data cubes. Weilbacher et al. (2014, 2020) provide a complete overview of these processes in the various operational stages. The MUSE data cube of NGC 5728 is publicly accessible in the ESO Archive², under program ID 097.B-0640 (PI: D. Gadotti). For further details on the data reduction performed on this data cube, see Mingozzi et al. (2019).

2.2. JWST data

The data processing pipeline for MIRI Medium Resolution Spectroscopy (MRS) on JWST meticulously processes raw observational data through several sophisticated stages to produce scientifically exploitable outputs. Initially, the Detector1 pipeline addresses raw, uncalibrated data from MIRI, correcting instrument-specific imprints and detector anomalies such as dark currents and CRs. This results in slope images that are corrected yet remain uncalibrated, representing the initial rate of signal change. Following this, the Spec2 pipeline performs calibration corrections including flat fielding, wavelength calibration, and spatial adjustments, alongside background subtraction to eliminate sky noise and other irrelevant nonastronomical signals, producing precisely calibrated 2D spectral images. Conclusively, the Spec3 pipeline combines data from multiple exposures and dither configurations, thereby enhancing the signal-to-noise ratio (S/N) and expanding both the spatial and spectral coverage. This process produces fully calibrated 3D data cubes for each observed target, meticulously prepared for comprehensive scientific exploration.

In our study, we employed uncalibrated data of NGC 5728 obtained under the JWST Proposal ID 1670 (PI: T. Taro Shimizu), which are publicly available. These data were retrieved from the ESA JWST Science Archive³, where they were downloaded and processed through all three stages of the JWST pipeline using version 1.16.0 and the Calibration References Data System (CRDS) context “jwst_1295.pmap”. The resulting data, spanning the continuous wavelength range of 4.9–27.5 μm across all MRS channels and sub-bands show no substantial differences from the Level 3 products available in the JWST archive or from recent studies that initially used different versions of the pipeline and CRDS contexts (Pereira-Santaella et al. 2022; Davies et al. 2024). The corresponding datasets are also available through the Mikulski Archive for Space Telescopes (MAST)⁴.

2.3. Spectrum extraction

The primary objective of this study was to investigate the parameters influencing the excitation of gas across different regions of NGC 5728, and specifically the combined effects of CRs with photoionization. Our approach involved selecting specific regions within NGC 5728 and extracting spectra from the MUSE and JWST data cubes using circular apertures. Each aperture was carefully positioned to assess the impact of CRs in both star-forming and jet-influenced areas.

We selected apertures that encompass gas interacting directly with jets and/or outflows, as well as regions indicative of ongoing star formation in the vicinity of NGC 5728. This was achieved by identifying areas exhibiting strong emission across multiple emission lines, as illustrated in Fig. 1. To ensure a robust comparison between the optical and MIR data, we adopted circular apertures with a radius of 0.7″ (∼138 pc), which closely matches the full width at half maximum (FWHM) of the point spread function (PSF) in the reddest channels of the MIRI and provides a suitable compromise between spatial resolution and S/N. Moreover, in Fig. 1, the selected apertures are depicted in varying shades of purple to indicate their distance from the nucleus. The nucleus itself (“N”) is shown in deep purple, while apertures at larger distances are labeled with increasing numbers and displayed in progressively lighter shades, ending with pale lilac for the most distant regions.

Fig. 1. Apertures chosen to extract spectra from the spectroscopic images of NGC 5728. These are depicted in different shades of purple from deep purple being the nuclear aperture “N” to pale lilac in increasing distance, also noted with numbers. The apertures are drawn over (a) the continuum in Channel 1 of MIRI, extracted in the rest-frame range of 5.90–6.05μm and over [Ar II]λ7 μm, [S IV]λ10.5 μm, [Ne II]λ12.8 μm, [Ne V]λ14.9 μm, [S III]λ18.7 μm, Hα, and [O III]λ5007Å (b–h) line emission maps with the stellar continuum subtracted. Panel (i) presents a schematic illustration of the adopted jet axis together with the three color-shaded regions used to group the apertures (see Sect. 2.3). These regions are defined by the projected vertical distance from the jet axis, relative to the aperture radius: JR (innermost, cyan) includes apertures within one radius (nuclear and/or jet regions); IJR (intermediate, lime green) covers one to two radii (intermediate jet-affected regions); and NJR (outermost, magenta) includes apertures beyond two radii (non-jet or star-forming regions). The multiwavelength data are aligned based on the astrometry of the nucleus.

In Panel (i) of Fig. 1, aperture observations are grouped into three categories based on their projected vertical distance from the jet axis, which is defined by the dichotomy line of the ionization cones shown in the middle panel of Fig. 12 in Shimizu et al. (2019). Regions within one aperture radius from the axis are considered nuclear-jet regions (JR), shaded cyan in panel (i) and marked with squares in subsequent figures. Apertures located beyond three and a half radii are classified as non-jet regions (NJR), more likely dominated by star formation, and are shaded magenta in panel (i) and represented by stars in other plots. Those situated between these two borders are regarded as intermediately jet-affected regions (IJR), shaded lime green in panel (i) and shown as diamonds elsewhere. If an aperture overlaps multiple regions, it is classified according to where the majority of its area lies.

Furthermore, from each aperture, we extracted the convolved spectra using MPDAF (Bacon et al. 2016) for the optical and PHOTUTILS (Bradley et al. 2024) for the infrared range. Both the continuum subtraction and the modeling of emission and absorption lines in the optical convolved spectra were performed using PYPLATEFIT (Bacon et al. 2023; Tremonti et al. 2004; Brinchmann et al. 2004), which fits and subtracts simple stellar population models (Bruzual & Charlot 2003; Brinchmann et al. 2013) as part of its integrated fitting procedure. Infrared emission lines, extracted from the same regions, were fitted with a single Gaussian component using the LMFIT package (Newville et al. 2015), while the continuum subtraction involved fitting a polynomial. Additionally, the spectra were smoothed using a median filter prior to fitting to reduce noise and eliminate substantial stellar and AGN contributions. A summary of the aperture properties–including positions, minimum S/N values, projected distances to the nucleus, and classifications based on Fig. 1i–is provided in Appendix A. All emission lines were initially selected for fitting by our routines using a 3σ detection threshold, with a few visually confirmed faint detections included in Table A.1. Finally, all emission lines detected were fitted in NGC 5728’s rest-frame wavelengths.

3. Methods

3.1. Constraining the parameter space

We constrained the ionization parameter, O/H, and N/O ratios for regions within NGC 5728 using the code HII-CHI-MISTRY⁵, developed by Pérez-Montero (2014). This tool applies Bayesian-like statistics to compare observed emission-line ratios against extensive grids of photoionization models. We specifically used the optical and infrared subroutines optimized for AGN (Pérez-Montero et al. 2019; Pérez-Díaz et al. 2022).

In the optical regime, HII-CHI-MISTRY relies on emission lines such as [O II]λ3727 Å, [Ne III]λ3868 Å, [O III]λ4363 Å, [O III]λ4959 Å, [O III]λ5007 Å, [N II]λ6584 Å, and [S II]λλ6716, 6731 Å. However, the MUSE spectral coverage (4800–9300 Å) excludes [O II]λ3727 Å and [Ne III]λ3868 Å. Thus, given the spectral coverage of MUSE, our analysis is necessarily limited to all the optical emission lines that fall within its wavelength range–specifically, [O III]λ4959 Å, [O III]λ5007 Å, [N II]λ6584 Å, and [S II]λλ6716, 6731 Å.

In the infrared regime, HII-CHI-MISTRY (Fernández-Ontiveros et al. 2021; Pérez-Díaz et al. 2022) employs emission lines including HIλ4.05 μm, HIλ7.46 μm, [S IV]λ10.5 μm, HIλ12.4 μm, [Ne II]λ12.8 μm, [Ne V]λ14.9 μm, [Ne III]λ15.5 μm, [S III]λ18.7 μm, [Ne V]λ24.3 μm, [O IV]λ25.9 μm, [S III]λ33.7  μm, [O III]λ52 μm, [N III]λ57 μm, [O III]λ88 μm, [N II]λ122 μm, and [N II]λ205 μm. While MIRI does not extend to the far-infrared lines included in the full HII-CHI-MISTRY grid, such as [O III]λ88 μm and [N II]λ122 μm, it still provides access to all the MIR emission lines needed for our analysis.

We use HII-CHI-MISTRY (v5.22 optical, v3.01 IR) with a double composite AGN continuum (α_UV = −1, α_OX = −0.8) and no ionization parameter constraint. The resulting parameter space, summarized in Table 1, guides our simulations. In the IR, the N/O ratio is unconstrained due to the lack of far-IR lines. See Section 3.2 and Table 2 for details. Furthermore, these results support our initial assumption of solar metallicity (1 Z_⊙), ensuring consistency with our previous analysis in KFDS25, and align with independent estimates reported by Pérez-Díaz et al. (2021). Finally, in Table 2 we present the relevant elemental abundances adopted as well.

3.2. CLOUDY models

Initially, we employed the dust-free AGN models detailed in section 3.3 of KFDS25. However, these models struggled to accurately reproduce high-ionization MIR emission lines such as [Ar V]λ13.1 μm, [Ne V]λ14.3 μm, [Ne V]λ24.3 μm, [Ne VI]λ7.7 μm, and [O IV]λ25.9 μm. To resolve this issue, we created a new grid of AGN models using CLOUDY v23.01 (Ferland et al. 2013, 2017; Chatzikos et al. 2023) with the intermediate AGN spectral energy distribution (SED), “AGN_Jin12_Eddr_mid.sed”, described by Ferland et al. (2020), also shown in teal color in Fig. B.1. This intermediate SED is a new option introduced in CLOUDY 23.01, while previous versions of CLOUDY (version 22.01 was employed in KFDS25) did not include this SED. We selected this intermediate SED because it provides a balanced and observationally validated scenario between lower Eddington ratio SEDs, which exhibit weaker far-UV emission and harder X-ray spectra, and higher Eddington ratio SEDs, marked by enhanced far-UV emission and softer X-ray spectra.

As in KFDS25, our CLOUDY simulations use a fixed hydrogen column density of N_H = 10²⁴ cm⁻², with the default temperature-based stopping criterion disabled. This ensures each model extends into cooler, deeper layers–crucial for capturing emission in strongly CR-heated models, where higher CR ionization rates produce more extended ionization structures. Using a smaller column or the default stopping criterion would prematurely truncate the emission from cooler regions.

Since HII-CHI-MISTRY utilizes AGN photoionization templates generated with CLOUDY 17.01 (see Fernández-Ontiveros et al. 2021; Pérez-Díaz et al. 2022), this is not identical to the new AGN SED described here, which provides improved modeling of high-ionization MIR lines. We tested HII-CHI-MISTRY in the MIR, by explicitly excluding the high-ionization lines [Ne V], and [O IV], to assess the impact on the derived metallicity and ionization parameter ranges. The exclusion of high-ionization lines does not introduce significant systematic shifts in the derived abundances or ionization parameter beyond the initial uncertainties, confirming that the AGN templates used in the HII-CHI-MISTRY routines remain reliable for metallicity and ionization parameter estimates. This robustness is consistent with the findings of Fernández-Ontiveros et al. (2021) and Pérez-Díaz et al. (2022), who demonstrate that HII-CHI-MISTRY produces stable metallicity and ionization parameter estimates even when certain diagnostic lines are absent or when only a subset of the MIR emission lines is available. Finally, in Section B, we show that the intermediate AGN template yields results that are identical to those obtained with the previous template in KFDS25.

The AGN model grid used is detailed in Table 2 and covers the parameter range previously outlined in section 3.3 of KFDS25, while the adopted abundances correspond to solar values (Asplund et al. 2009). These elemental abundances align with abundances derived via HII-CHI-MISTRY for NGC 5728 and are consistent with abundance estimates for AGN narrow line regions (NLR) in the local Universe (e.g. Pérez-Díaz et al. 2021; Dors 2021; Dors et al. 2022; Pérez-Díaz et al. 2022, 2024) Finally, since CLOUDY simulations do not include extinction, we applied Calzetti’s law (Calzetti 2001) to attenuate the modeled fluxes of the optical emission lines, enabling a direct comparison with observations.

4. Results

4.1. Optical line diagnostic diagrams

With the adoption of the intermediate SED, our models consistently reproduce the low-ionization optical lines–[N II], [S II], and [O I]–as effectively as in our previous work (KFDS25), across the same parameter space as shown in Table 2.

Furthermore, the intermediate SED shows enhanced consistency with MIR diagnostics, discussed in detail in Section 4.2, and thus represents a more complete modeling approach.

In Fig. C.1, it is evident that all jet-affected regions–whether strongly or moderately influenced-fall above the Kewley line, mainly in the Seyfert domain. CR ionization rates around 10⁻¹³ s⁻¹ better represent observations from nuclear and jet affected areas in the [S II] and [O I] BPT diagrams, while a rate of approximately 10⁻¹² s⁻¹ best models the parameter space of the observed [N II] ratios. These findings are in agreement with the results we acquired in KFDS25.

In KFDS25, we demonstrate that regions traditionally characterized as Seyfert and/or LINER sources can also be explained using the star-forming models incorporating CR ionization which led to a new parameter space demarcation line, the SFζ line. This boundary captures the upper limit of line ratios expected from star-forming regions influenced by elevated CR ionization rates, offering a more physically motivated and complementary perspective to classical starburst boundaries. Applying this framework allowed us to reassess excitation mechanisms in galaxies without the need to immediately attribute elevated line ratios to AGN activity or extreme metallicities (Groves et al. 2004).

Region 17 of NGC 5728 is indicated as a white star in Figure C.1, with the white color denoting its location at the greatest distance from the nucleus and the star symbol highlighting its star-forming nature. In the BPT diagrams, this region is located well within the star-forming zone, below the Kewley line and even further below the SFζ line, suggesting minimal influence from high CR ionization. This interpretation is further supported by its position away from the jets (magenta-shaded area in Figure 1i). As such, region 17 is important, as it exhibits line ratios that are consistent with pure star formation and minimal CR impact and provides a baseline for contrasting with jet-impacted or CR-enhanced regions within NGC 5728. Subsequently, regions 10 and 13, marked in lilac and almost white stars in Figure C.1, occupy the in-between zone above the Kewley line, while still falling bellow the SFζ boundary shown in Figure 2. These regions originate from areas closer to the jets and are therefore possibly affected by CRs. This confirms that regions such as 10 and 13 are representative of how CR-enhanced star-formation models can explain elevated line ratios in certain parts of NGC 5728, bridging the gap between classical star-forming excitation and AGN-like signatures without requiring a dominant AGN contribution.

Finally, regions closer to the nucleus and the jets, displayed as cyan squares in Fig. 2, lie above the SFζ line, making these areas consistent with CR-driven and AGN ionization scenario. Moreover, thin lime green diamonds in Fig. 2, representing intermediate regions far from the jets, lie near the SFζ line and can be reproduced either by AGN or star-formation models with CRs. These boundaries are useful in characterizing the intermediate region that is consistent with both the AGN-dominated models and those involving star formation combined with CR ionization (ζ_CR ∼ 10⁻¹³, s⁻¹). Our analysis of NGC 5728 underscores the value of such revised boundaries for interpreting emission-line diagnostics, particularly in spatially resolved regions influenced by CRs.

Fig. 2. BPT diagrams depicting [N II]/Hα, [S II]/Hα, and [O I]/Hα ratios. The observations from NGC 5728, and NGC 1320 are marked with stars, and a red diamond, respectively. The cyan squares, lime green diamonds, and magenta stars represent jet-affected (JR), intermediate jet-affected (IJR), and non-jet-affected (NJR) regions, respectively, and correspond to the shaded areas in Fig. 1. The Kewley and Schawinski lines are indicated with solid and dashed black lines, respectively, while the Koutsoumpou (SFζ) line, is depicted by the solid red line. In the background we show the line ratios measured for nearby galaxies from the Sloan Digital Sky Survey (SDSS) data release 7 (Abazajian et al. 2009).

4.2. MIR and optical–MIR hybrid line diagnostic diagrams

The first diagnostic plot we selected is [Ne V]λ14.3 μm/ [Ne II]λ12.8 μm versus [Ne III]λ15.6 μm/[Ne II]λ12.8 μm, depicted in the first column of Fig. 3. This frequently used diagnostic tool has been employed to compare observations with various models, such as pure photoionization models (Feltre et al. 2023) or photoionization and shock models (Zhang et al. 2024). Given the value of this plot, we reproduced it to verify that our models sufficiently described the parameter space. Our analysis reveals that as ζ_CR increases from 10⁻¹⁶ s⁻¹ in the top row to 10⁻¹² s⁻¹ in the bottom row, the observational data are best represented by models with ζ_CR values of ∼10⁻¹⁴ s⁻¹. These model cases spatially coincide with those presented in figures such as fig. 5 in Feltre et al. (2023) and Fig. 15 in Zhang et al. (2024), rendering our models good for reproducing the same parameter space. Subsequently, the observations in this plot are effectively captured by models with ζ_CR ∼ 10⁻¹⁴ − 10⁻¹³ s⁻¹. Some regions more directly (square) or indirectly (thin diamonds) influenced by the jets, are adequately modeled at ζ_CR ∼ 10⁻¹³ s⁻¹ or even 10⁻¹² s⁻¹.

Fig. 3. Diagrams with the AGN photoionization models compared with the observations from the selected apertures in NGC 5728 (Fig. 1). The different shades of purple ranging from deep purple to pale lilac and/or white, represent increasing distance from the nucleus, as also denoted by numbers, with “N” being the nuclear aperture. The different shapes-square, thin diamond, and star-represent the nucleus and/or jet impacted, intermediate and distant areas, respectively. The different shades from white to deep red represent the range of ionization parameter values, −3.5 ≤ log U ≤ −1.5. All the models have solar abundances. The panels from top to bottom correspond to ζCR = 10−16 s−1,  10−15 s−1,  10−14 s−1,  10−13 s−1, and 10−12 s−1, respectively.

The diagnostic plot of [Ne III]λ15.6 μm/([Ne II]λ12.8 μm+[Ar II]λ7.0 μm) against [Ne V]λ14.3 μm/ ([Ne III]λ15.6 μm+[Ar III]λ8.9 μm) displayed in the second column of Fig. 3, was chosen to highlight the substantial shift in the positions of the models within the parameter space as the CR ionization rate increases from 10⁻¹⁶ s⁻¹ in the top row to 10⁻¹² s⁻¹ in the bottom row. This pronounced change makes the selected ratios particularly effective for observing the maximal influence of CRs. The shift of the models toward the lower left corner of the plot leads us to conclude that the low ionization lines, namely [Ne II], [Ar II], [Ne III], and [Ar III], are significantly impacted by CRs. For this diagnostic, we note that models deviating from the folding trend (rightmost, bottom row in Fig. 3) correspond to the highest CR ionization rate (10⁻¹² s⁻¹) and the lowest densities (log n_H ≤ 0.1). At such low densities, strong CR ionization can suppress cooling and low-ionization lines (e.g., [Ar II], [Ne II]), disrupting the trend observed at higher densities in this diagram.

Combining optical and MIR emission-line ratios could represent a promising approach for disentangling the contributions of AGN, star formation, CRs, and shocks in composite galaxies. Optical diagnostics such as [O I]/Hα are sensitive to the hardness of the ionizing radiation and can trace shocks, while MIR ratios such as [O IV]/[Ne III] probe deeply embedded, high-ionization regions inaccessible in the optical due to dust. As shown by Feltre et al. (2023), diagnostic diagrams that combine these regimes–notably [O I]/Hα versus [O IV]/[Ne III]–offer enhanced diagnostic power, allowing degeneracies present in single-wavelength analyses to be resolved. In the third column of Fig. 3, we examine the impact of the CR ionization rate on this combination of line ratios. We find that increased CR ionization rates enhance [Ne III] and [O I] emission, resulting in trends that differ from those produced by shocks (see Fig. 7 of Feltre et al. (2023)). While both shocks and CRs can elevate [O I]/Hα, CR-driven models yield a broader distribution that better matches the observed parameter space (see Fig. 3). In contrast, AGN combined with star-formation models tend to underpredict, and AGN combined with shock models mostly overpredict, the observed ratios, providing less overlap with the data (Fig. 7 of Feltre et al. 2023). Nonetheless, this approach has limitations, as optical and MIR lines may trace different regions within the ionized clouds, sampling different temperatures and densities. Additionally, uncertainties in extinction corrections can affect the reliability of a combined diagnostic. The [O I]/Hα vs [O IV]/[Ne III] diagnostic is promising for distinguishing CR from shock excitation although it is strongly metallicity-dependent, as [O I]/Hα directly traces O/H. Since both our models and those of Feltre et al. (2023) use fixed metallicity, we likely underestimate the true spread in this ratio and may overlook CR-shock degeneracies. Exploring additional mixed optical-MIR diagnostics in the future will provide powerful tools for investigating the interplay between star formation, AGN activity, shocks, and CRs in complex galactic environments where these mechanisms coexist.

Additionally, the majority of observed emission-line ratios are accurately represented by models at ζ_CR ∼ 10⁻¹⁴ s⁻¹. Notably, star-forming regions, depicted with white and/or lilac stars based on their proximity to the nucleus, are well modeled without the influence of CRs, a behavior that is also present in the position of these stars in the BPT diagrams (Fig. C.1). However, data represented by squares and thin diamonds, which are either directly or indirectly affected by jets, are accurately described by even the highest CR ionization models, as illustrated in the second column and last row of Fig. 3. Finally, in Appendix D, we present additional MIR emission lines, including mid-ionization species such as [S IV]λ10.5 μm and [S III]λ18.7 μm, as well as high-ionization lines such as [Ar V]λ7.9 μm, [Ar V]λ13.1 μm, [Ne V]λ14.3 μm, [Ne V]λ24.3 μm, [Ne VI]λ7.7 μm, and [O IV]λ25.9 μm. The diagnostic plots presented across Fig. 3 and Fig. D.2 effectively lay out the varying impacts of CRs on emission lines with different ionization potentials. We confirm that while high-ionization lines are mainly influenced by photoionization and are thus not good indicators of CR presence, low-ionization lines show significant sensitivity to CRs. We further study the link between CRs and MIR emission by examining the temperature of nebular gas and the emissivity of different lines within the depth of an emission region in Section 4.3.

4.3. Gas stratification diagrams

In our study in KFDS25, we found that photoionization governs the illuminated face of the cloud (i.e. the nebular layers are defined with respect to the ionizing source), whereas CRs play a more critical role in the shielded regions where photoionization is minimal. Specifically, CR rates ∼10⁻¹³ s⁻¹ significantly impact the thermal dynamics of the ionized gas, leading to the formation of a secondary, low-ionization layer of warm ionized gas (∼8000 K) well beyond the region dominated by photoionization. This inner layer contributes to the emission of low-ionization lines such as [N II]λ6584 Å, [S II]λλ6716,6731 Å, and [O I]λ6300 Å. In contrast, the emissivities of higher ionization lines such as [O III]λ5007 Å, Hα, and Hβ, are only slightly increased (see Figs. 8, 9, and 10 in KFDS25). This analysis was replicated for different emission lines at MIR wavelengths to determine if low-ionization MIR emission lines were similarly affected as presented in Section 4.2.

As illustrated in Fig. 4, the interaction of CRs with gas significantly influences the MIR emission lines. For an initial hydrogen density of n_H = 100 cm⁻³ and an ionization parameter log U = −3.0, we observe varied sensitivities across different MIR emission lines under increasing CR ionization rates from 10⁻¹⁶ s⁻¹ to 10⁻¹² s⁻¹. We find that while photoionization influences the electron temperature at the cloud’s illuminated face, high CR rates ( ≳ 10⁻¹⁴ s⁻¹) dictate the excitation in deeper layers. Specifically, in the second and third row of the first column of Fig. 4, [Ar II]λ7 μm and [Ne II]λ 12.8 μm, with ionization potentials of 15.76 eV and 21.56 eV respectively, exhibit the most pronounced changes in emissivity due to higher ζ_CR, highlighting their strong sensitivity to CRs. The subsequent panels of the second column within Fig. 4 reveal that [S III]λ 18.7 μm, [Ar III]λ 8.9 μm, and [Ne III]λ 15.5 μm, having ionization potentials of 23.33eV, 27.63 eV, and 40.96 eV respectively, also experience considerable emissivity enhancement, though to a lesser degree than [Ar II] and [Ne II]. Conversely, displaying a gradient of decreasing responsiveness to CRs, [S IV]λ 10.5 μm with an ionization potential of 34.79 eV, shows some effects; while [Ar V], with an ionization potential of 59.81 eV, exhibits only minimal changes; and [Ne V], with an ionization potential of 97.12 eV, shows no sensitivity to CRs. These results indicate that the highest-ionization emission lines are largely unaffected by CR-driven ionization within the cloud.

Fig. 4. Temperature and line emissivity vs. depth in the simulated cloud for AGN models, for an initial density of nH = 100 cm−3, and for ζCR = 10−16 s−1, 10−15 s−1,  10−14 s−1,  10−13 s−1, and 10−12 s−1, and log U = −3.0. The different panels a-f correspond to kinetic temperature and the emissivity of [S III]18.7μm, [S IV]10.5μm,[Ar II]7μm,[Ar III]8.9 μm,[Ar V]7.9μm,[Ne II]12.8μm,[Ne III]15.5μm, and [Ne V]14.3μm, respectively. The shaded area indicates the approximate region in which CR heating becomes dominant. The CR-dominated area is depicted in lilac in the temperature plot and in teal in the emissivity plots.

This differentiation in sensitivity to CR rates among different emission lines and different ionization states of the same line is deeply rooted in the ionization energy required for each particular emission. The behavior of low-ionization lines in the MIR spectrum exhibits similar trends found in the low-ionization lines in the optical spectrum, as documented in our findings for KFDS25. Specifically, [N II]λ6584 Å, [S II]λλ6716,6731Å, and [O I]λ6300 Å in the optical spectrum, with ionization potentials of 14.53 eV, 10.36 eV, and 13.62 eV, respectively, exhibit increased emissivity under higher CR ionization rates, paralleling the effects seen in MIR lines such as [Ne II]λ12.8 μm and [S III]λ18.7 μm. This finding reinforces the interconnected nature of optical and MIR spectral behavior, demonstrating a unified response to CR ionization, which could be exploited in astrophysical environments not optically thin in the optical, such as dust-obscured nuclei of galaxies, star-forming regions, dense molecular clouds, or circumstellar disks.

5. Discussion

5.1. CR effect on MIR low ionization lines

CRs significantly influence the ionization and the heating of the ISM in environments such as AGN and starburst galaxies (McKee 1989; Padovani et al. 2009, 2018; Gabici 2022). Our analysis in KFDS25, demonstrates that the impact of CRs on low-ionization optical lines such as [N II], [S II], and [O I], is more evident than on higher-ionization lines such as Hα, Hβ, and [O III], due to a secondary low-ionization CR-induced layer beyond photoionization-dominated regions. We analyze that the effect of CRs on optical emission lines introduces new limits in the classical line-ratio BPT diagnostic diagrams and, as showcased in section 5.4 of KFDS25, can also affect metallicity and ionization parameter estimates.

Based on Section 4.2, the MIR emission lines most sensitive to CR ionization are [Ar II] and [Ne II]. This conclusion is supported by both their relatively low ionization potentials, 15.76eV, and 21.56eV, respectively, and their emissivity profiles in Fig. 4, which exhibit a pronounced enhancement with increasing ζ_CR. In contrast, lines from higher-ionization species such as [Ne III], [Ar III], [S IV], and even more in [Ar V], [Ne V] display minimal changes, remaining largely governed by photoionization (see Fig. 4). This differential response arises from the physical stratification within the cloud. In the layers closer to the ionizing source, where photoionization from the AGN dominates, the intense radiation field enhances the emission of high-ionization lines such as [Ar V] and [Ne V], while the ionic fractions of low-ionization species such as [Ar II] and [Ne II] remain low due to their efficient ionization to higher states. Deeper into the cloud, however, the radiation field becomes increasingly attenuated, allowing high-ionization ions to recombine and the ionic fractions of low-ionization species, like [Ar II] and [Ne II], to increase. In these shielded regions, CRs begin to dominate the ionization balance, sustaining a partially ionized phase and efficiently boosting emission from these low-ionization lines. In Fig. 5, Ar and Ne ionic fractions versus depth, clearly show that [Ar II] and [Ne II] increase in abundance in regions where CR heating is dominant and the higher ionization stages are no longer sustained. Hence, the combined ratio [Ar II] + [Ne II] offers a powerful diagnostic, effectively isolating the CR-influenced component of the ionized gas and distinguishing it from photoionized regions. This makes it a particularly suitable tool for probing environments where CRs are expected to be present.

Fig. 5. Ionic fraction vs. depth in the simulated cloud for AGN models, for an initial density of nH = 100 cm−3, and for ζCR = 10−16 s−1, 10−15 s−1,  10−14 s−1,  10−13 s−1, and 10−12 s−1, from left to right. Top row: Ar ionic fractions. Bottom row: Ne ionic fractions. The blue-shaded area indicates the approximate region where CR heating becomes dominant.

Additionally, we tested the effect of dust by running extra CLOUDY simulations using standard prescriptions that include either ISM grains or polycyclic aromatic hydrocarbons (PAHs) (grains ISM and grains PAH, respectively). These dusty models produced MIR line ratios nearly identical to dust-free models. Low-ionization lines ([Ne II]λ12.8 μm, [Ar II]λ7.0 μm, [S III]λ18.7 μm) showed comparable emissivities overall, but their emission decreased faster with depth due to more efficient cooling when dust is present. Additionally, dust slightly increase the electron temperatures near the illuminated face of the cloud, where hydrogen recombination dominates, particularly at lower ζ_CR. In contrast, the impact of dust in highly-excited lines ([Ne V]λ14.3 μm, [Ar V]λ7.9 μm) is negligible, since these lines originate predominantly in regions close to the ionizing source where AGN photoionization strongly dominates.

Following the discussion in KFDS25, we estimated the expected CR ionization rate for the central region of NGC 5728 using the measured radio continuum fluxes at 1.4 GHz and 5 GHz (Singh et al. 2013) and assuming equipartition between the magnetic field and relativistic electron energy densities. For a spectral index of α ∼ 1.09, and a ∼20″ beam at 5 GHz, corresponding to a size of ∼1.85 kpc, we obtain an equipartition magnetic field of ∼4 μG. Adopting p ∼ 3.2 for the electron energy distribution spectrum, derived from the synchrotron slope via p = 2α + 1 (Rybicki & Lightman 1979), and following the methodology described in Gabici (2022) and KFDS25, we obtain a CR ionization rate of ∼10⁻¹¹ s⁻¹ for the nuclear region of NGC 5728. A lower value of p = 2.4, corresponding to the standard slope in optically thin synchrotron sources, would result in ∼10⁻¹⁴ s⁻¹. These values are consistent with the CR rates inferred from our optical and MIR models, which reproduce the observed ratios in those regions of NGC 5728 most affected by the jet. Deeper radio observations, at higher angular resolution and particularly targeting optically thin emission, would improve CR ionization rate estimates.

5.2. Diagnostic power of MIR versus optical

MIR emission lines offer a significant observational advantage over their optical counterparts due to their much lower sensitivity to dust extinction and temperature variations. As shown in table 3 of Wang & Chen (2019), the relative extinction decreases from A_λ/A_V = 1.0 at 0.5525 μm (Johnson V) to ∼0.025 at 8 μm (Spitzer [8.0]). Thus, an optical extinction of A_V = 10 mag results in ∼0.25 mag of extinction at 8 μm, representing a reduction by a factor of 40. Conversely, producing just 0.1 mag of extinction in the MIR requires an optical extinction of about A_V ∼ 4 mag. This steep attenuation gradient highlights the power of MIR diagnostics for investigating heavily obscured environments, where optical line emission is significantly diminished.

In addition to differences in dust extinction, another key distinction between optical and MIR diagnostics lies in their sensitivity to the electron temperature, T_e. As shaded in pale pink in the first panel of Fig. 4, our models exhibit a pronounced drop in T_e at a depth of approximately 10¹⁸ cm into the cloud. This temperature decline coincides with the onset of recombination, where the decreasing number of free electrons reduces the overall thermal energy, leading to substantial cooling of the gas. Optical lines, such as [S II] and [O I], are particularly sensitive to T_e, as their emissivities rely heavily on collisional excitation, as argued in KFDS25. Consequently, this sharp drop in T_e leads to a corresponding decline in their emissivity at greater depths. This behavior is evident in the predictions of our AGN models with ζ_CR ≳ 10⁻¹³ s⁻¹ for log U = −3, where the optical line intensities diminish steeply beyond the photoionized surface layers (see Figs. 9, and 10 in KFDS25). In contrast, the MIR lines such as [Ar II]λ7.0 μm and [Ne II]λ12.8 μm (see Fig. 4) exhibit a much flatter emissivity profile across the same depth range. This difference arises because IR fine-structure transitions have much lower excitation energies and are less sensitive to the thermal electron population (see also Fig. 1 in Fernández-Ontiveros et al. 2021). Therefore, even as the gas cools and recombines, MIR lines can continue to be emitted efficiently in the deeper, CR-affected layers, making them reliable probes of low-ionization gas in regions where optical diagnostics are suppressed.

Moreover, recent studies of large, diffuse ionized nebulae with uncertain origin (Lumbreras-Calle et al. 2025), highlight the importance of resolving ionization mechanisms in faint emission-line structures. These findings support the need for diagnostics capable of identifying the excitation mechanisms in spatially extended regions. In this context, hybrid optical–MIR diagnostics offer valuable potential for resolving such degeneracies. In particular, the [O IV]/[Ne III] versus [O I]/Hα diagram emerges as a key tool.

5.3. CR excitation in AGN MIR diagnostics

CRs enhance the emission-line fluxes from low-ionization species relative to the pure photoionization case, leading to a systematic decrease in the high-to-low ionization line ratios. This trend arises from the unique ability of CRs to penetrate far deeper into cloud structures than ionizing photons, maintaining the ionization and excitation of low-ionization species in partially shielded regions where the radiation becomes ineffective.

Figure 6 illustrates this effect for the diagnostic diagram of [Ne V]/[Ne II] versus [Ne III]/[Ne II] through our models, which assume solar abundances and span a wide range in the ionization parameter (−3.5 ≤ log U ≤ −1.5) and hydrogen density (0 ≤ log n_H ≤ 4). The model predictions are presented as solid contours enclosing 10%, 50%, and 90% of all the predicted ratios for these input parameters, thus excluding outlier values. The contours colored in green, teal, blue, purple, and red, represent the increasing CR ionization rate, ζ_CR = 10⁻¹⁶, 10⁻¹⁵, 10⁻¹⁴, 10⁻¹³, and 10⁻¹² s⁻¹, respectively. These models are compared with the distribution of emission-line ratios measured for nearby AGN. The gray data points in Fig. 6 correspond to Spitzer/IRS spectroscopic observations for a sample of nearby Seyfert galaxies and LINERs, compiled by Fernández-Ontiveros et al. (2021) from the literature (Dudik et al. 2007, 2009; Goulding & Alexander 2009; Tommasin et al. 2008, 2010; Pereira-Santaella et al. 2010; Fernández-Ontiveros et al. 2016). This full dataset most likely captures the diversity of MIR line ratios across the local population of active nuclei, including AGN photoionization, shocks, AGN and/or starburst composites, and jet-dominated nuclei.

Fig. 6. Diagnostic diagrams depicting the area covered by AGN models with solar abundances for −3.5 ≤ log U ≤ −1.5 and for 0 ≤ log nH ≤ 4. The solid contour lines map regions containing 10%, 50%, and 90% of the models for ζCR = 10−16 s−1, 10−15 s−1,  10−14 s−1,  10−13 s−1, and 10−12 s−1, from left to right in green-, teal-, blue-, purple-, and red-colored contours, respectively. In the background are archival Spitzer/IRS measurements of Seyfert and LINER nuclei (Dudik et al. 2007, 2009; Goulding & Alexander 2009; Tommasin et al. 2008, 2010; Pereira-Santaella et al. 2010; Fernández-Ontiveros et al. 2016, 2021), along with recent nuclear line ratios from JWST/MIRI observations from the literature (Pereira-Santaella et al. 2022; Álvarez-Márquez et al. 2023; Armus et al. 2023; Alonso Herrero et al. 2024; Zhang et al. 2024; Hernandez et al. 2025; Goold et al. 2024; Hermosa Muñoz et al. 2025, Goold et al. in prep., Fernández-Ontiveros et al. in prep.)

Photoionization-dominated models with log U ≳ −2.5 are in agreement with the highest [Ne V]/[Ne II] and [Ne III]/[Ne II] ratios observed in the AGN population (Fig. 6a). As the CR ionization rate increases, the models extend toward progressively lower [Ne III]/[Ne II] values, highlighting the increasing contribution of CRs to low-ionization lines (Figs. 6c–6e). This contrasts with the behavior in pure AGN photoionization models at decreasing ionization parameter values (log U ≲ −3), which primarily affect the [Ne V]/[Ne II] ratio and deviate from the observed ratio (e.g. Fig. 6a). Thus, CR ionization may be a relevant ionization mechanism, together with shocks and star-formation contribution (e.g. Feltre et al. 2023), to explain the lower-excitation end of the observed distribution. The latter is mostly populated by LINER galaxies, which are expected to exhibit a higher incidence of compact jets compared to the Seyfert population (Baldi 2023). In particular, the distribution of observed ratios is well reproduced by models with CR ionization rates in the (ζ_CR ≳ 10⁻¹⁴–10⁻¹³ s⁻¹) range, in agreement with measurements obtained in the central region of active galaxies (e.g. González-Alfonso et al. 2013, 2018). This suggests that CRs may provide a dominant or complementary excitation mechanism in low-luminosity AGN where the radiation field alone fails to account for MIR diagnostics.

The nucleus of NGC 5728, consistent with AGN photoionization-dominated ratios (Figs. 2, 3, and C.1), is marked by a cyan square in Fig. 6. The locations of NGC 1068 and Centaurus A, whose optical line ratios were previously analyzed in KFDS25, are also highlighted by thin diamonds in magenta and lime green, respectively. Additional Seyfert and LINER galaxies with published line fluxes from nuclear measurements in recent JWST/MIRI observations (Pereira-Santaella et al. 2022; Álvarez-Márquez et al. 2023; Armus et al. 2023; Alonso Herrero et al. 2024; Zhang et al. 2024; Hernandez et al. 2025; Goold et al. 2024; Hermosa Muñoz et al. 2025, Goold et al. in prep., Fernández-Ontiveros et al. in prep.)

are indicated in blue and pink, respectively. The nucleus of Centaurus A, which exhibits prominent jet activity, is far from the AGN photoionization locus and can be described by models with CR ionization rates of ζ_CR ∼ 10⁻¹⁴–10⁻¹³ s⁻¹. Finally, it is noteworthy that Mrk231, a Seyfert galaxy exhibiting both jets and starburst activity (Aalto et al. 2015; Silpa et al. 2021), falls within the region of Fig. 6 that can be modeled with ζ_CR ≳ 10⁻¹⁴, s⁻¹. It is depicted as a blue thin diamond, and corresponds to the lowest [Ne III]/[Ne II] ratio among the JWST-observed Seyfert galaxies considered here, as reported by Alonso Herrero et al. (2024), supporting the idea that CRs play a role in shaping the emission-line properties of such galaxies.

6. Summary

We conducted a comprehensive study of NGC 5728 by integrating VLT/MUSE optical and JWST/MIRI MIR observations to investigate the impact of CRs on the ISM in AGN. Using updated CLOUDY models that include CR ionization and a more adequate AGN ionizing spectrum, we explored a broad parameter space in the ionization parameter, gas density, and CR ionization rate. Our models, assuming solar abundances, successfully reproduce both optical and MIR emission-line ratios without requiring supersolar metallicities.

We find that low-ionization lines such as [Ne II] and [Ar II] are particularly sensitive to CRs, while high-ionization lines such as [Ne V] and [Ar V] remain governed by photoionization. This allows high-ionization MIR diagnostics to serve as robust tracers of AGN activity, whereas low-ionization MIR lines provide valuable insight into the CR-dominated regime of the ISM. CR ionization rates ranging from ζ_CR ∼ 10⁻¹⁴–10⁻¹³ s⁻¹ are found to reproduce the observed MIR emission in regions affected by jets and outflows in NGC 5728, consistent with evidence of nonthermal activity from the synchrotron fitting of radio originating from the nucleus. Additionally, by combining optical and MIR diagnostics, we show that diagnostic diagrams involving [O I]/Hα and [O IV]/[Ne III] may offer a viable route toward disentangling the degeneracy between CR- and shock-induced excitation.

Finally, our models successfully reproduced the MIR line ratio parameter space of Centaurus A and NGC 1068, two benchmark systems with elevated CR activity previously examined in KFDS25. This demonstrates that photoionization along with CRs can account for the excitation conditions observed across both optical and MIR wavelengths. Centaurus A is consistent with ionization rates as high as ζ_CR ∼ 10⁻¹² s⁻¹, while NGC 5728 and NGC 1068 correspond to rates ranging from ζ_CR ∼ 10⁻¹⁴ to 10⁻¹³ s⁻¹. These results extend the conclusions of KFDS25 by demonstrating that AGN models with elevated CR ionization rates provide a plausible explanation for the excitation observed in low-luminosity AGN.

Notably, AGN-like emission-line ratios are consistent with high CR ionization rates (ζ_CR ≳ 10⁻¹⁴ s⁻¹), suggesting that CRs serve as a complementary excitation mechanism in such environments. Altogether, the consistency across wavelengths supports the scenario proposed in KFDS25, in which CRs, possibly originating from jets or SNRs, significantly contribute to both the ionization and heating of the gas in AGN.

Acknowledgments

We are grateful to the anonymous referee for their constructive comments for overall improving this work. We also sincerely thank Kameron Goold and the ReveaLLAGN collaboration for kindly providing their data for comparison with our models. EK gratefully acknowledges Prof. A. Mastichiadis for his early supervision and constructive support during the initial stages of this work, and Prof. D. Hatzidimitriou for her insightful feedback during the preparation of this manuscript. EK acknowledges full financial support by the State Scholarships Foundation (IKY) from the proceeds of the “N. D. Chrysovergis” bequest. JAFO acknowledges financial support by the Spanish Ministry of Science and Innovation (MCIN/AEI/10.13039/501100011033), by “ERDF A way of making Europe” and by “European Union NextGenerationEU/PRTR” through the grants PID2021-124918NB-C44 and CNS2023-145339; MCIN and the European Union – NextGenerationEU through the Recovery and Resilience Facility project ICTS-MRR-2021-03-CEFCA. This work is based on observations made with the NASA/ESA/CSA JWST, associated with program ID 1670 (PI: T. Taro Shimizu). The corresponding datasets are publicly available in MAST at the Space Telescope Science Institute, operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-03127.

References

Appendix A: Aperture Properties

Appendix B: Intermediate SED

Here we present the intermediate AGN spectral energy SED produced with CLOUDY 23.01 and introduced in Ferland et al. (2020) versus the SED used in KFDS25. The intermediate SED was adopted in our new model grid to better reproduce high-ionization MIR emission lines.

More precisely, in the grid of models used for the AGN cases in KFDS25, the SED of the central photoionizing continuum is assumed to consist of a large blue bump from the accretion disc, radio emission from synchrotron radiation in the jet, infrared emission from dust, a soft X-ray excess, and a power-law component at hard X-rays. This SED may be expressed mathematically by Equation (2) in KFDS25. In that formula, F_ν is the flux density as a function of frequency ν; α_uv and α_x are the UV and X-ray spectral indices, respectively; T_bb is the characteristic temperature of the big blue bump; T_IR is the temperature corresponding to the infrared cutoff of the big blue bump; and α_ox describes the optical-to-X-ray spectral index. In KFDS25 AGN models, we adopted the following values: α_uv = −0.5, α_x = −1.0, α_ox = −1.4, T_BB = 10⁵ K, and T_IR = 1.6 × 10³ K. In the frequency range between the SED peak and 10¹⁶eV, we also find a slope of approximately −1 for the blue curve in Fig. B.1.

For the intermediate SED, utilized in the present work, we adopt the similarly named in Ferland et al. (2020) SED, based on the analysis by Jin et al. (2012), energy-conserving SEDs for unobscured AGN, grouped according to their Eddington ratio L/L_Edd. The intermediate SED corresponds to a mean value of log(L/L_Edd)≃ − 0.55 and the template available in CLOUDY 23.01 displays a shallower slope of approximately −0.4 (i.e., α_uv ≃ −0.4) in the frequency range immediately after the SED peak up to 10¹⁷ eV.

Fig. B.1. AGN SEDs simulated with CLOUDY. The blue solid line represents the AGN SED used in KFDS25 while the teal solid line represents the AGN used in this paper. The y-axis is νFν in arbitrary scaling.

Appendix C: BPT diagnostics with the intermediate SED

Fig. C.1. BPT diagrams comparing AGN photoionization models with observations from selected apertures in NGC 5728 (Fig. 1). The BPT diagrams for [N II], [S II], and [O I] are shown on the left, middle, and right, respectively. The different shades of purple, ranging from deep purple to pale lilac and/or white, represent the increasing distance from the nucleus, as also indicated by numbers, with "N" denoting the nuclear aperture. The different shapes-square, thin diamond, and star-represent the nucleus and/or jet impacted, intermediate, and distant areas, respectively. The different shades from white to deep red represent the range of ionization parameter values, −3.5 ≤ log U ≤ −1.5. The red diamonds represent the measured line ratios for the photoionization-dominated Seyfert 2 nucleus in NGC 1320. The Kewley, Kauffmann, Schawinski, and Koutsoumpou (SFζ) lines correspond to the solid black, blue, red, and magenta lines, respectively. The panels from top to bottom correspond to ζCR = 10−14 s−1, 10−13 s−1, and 10−12 s−1, respectively.

Appendix D: Density estimates and extra MIR diagnostics

Both [Ne V] and [Ar V] ratios serve as density diagnostics; in the sampled parameter space, they exhibit ratios of 0.5 ≲ [Ne V]/[Ne V] ≲ 1.4 and 1 ≲ [Ar V]/[Ar V] ≲ 1.7, respectively (see Fig. D.1). These density values align well with the densities anticipated by the ratios obtained via our CLOUDY simulations, as it can be seen in the first and second column of Fig. D.2, and with the initial hydrogen densities of our models (see Table 2).

Fig. D.1. The [Ar V]13.1μm/[Ar V]7.9μm and [Ne V]24.3μm/[Ne V]14.3μm line ratios produced with PYNEB (Luridiana et al. 2014) as a function of electron density, ne, notated for different electron temperatures, Te.

Fig. D.2. Diagrams with the AGN photoionization models compared with the observations from the selected apertures in NGC 5728 (Fig.1). The different shades of purple, ranging from deep purple to pale lilac and/or white represent the increasing distance from the nucleus, as also denoted by numbers, with "N" corresponding to the nuclear aperture. The different shapes-square, thin diamond, and star, represent the nucleus and/or jet impacted, intermediate, and distant areas, respectively. The different shades from white to deep red represent the range of ionization parameter values, −3.5 ≤ log U ≤ −1.5. All the models have solar abundances. The panels from top to bottom correspond to ζCR = 10−16 s−1, 10−15, s−1, 10−14 s−1, 10−13 s−1, and 10−12 s−1, respectively.

In the first and second column of Fig. D.2, we present the [S IV]λ10.5μm/[S III]λ18.7μm ratio plotted against the Ar V]λ13.1μm/[Ar V]λ7.9μm and the [Ne V]λ24.3μm/[Ne V]λ14.3μm ratios, respectively. The ratios on the horizontal axes, [Ar V]λ13.1μm/[Ar V]λ7.9μm and [Ne V]λ24.3μm/[Ne V]λ14.3μm, are density tracers if the gas is assumed to be in thermodynamical equilibrium and, appear relatively stable as the CR ionization rate increases from the top to the bottom row. The vertical axis on both first and second column of Fig. D.2, the [S IV]λ10.5μm/[S III]λ18.7μm ratio, is mildly affected by the increasing ζ_CR, due to the sensitivity of [S III] to CRs and to a lesser degree that of [S IV] (see Fig. 4b, c). Finally, in the third column of Fig. D.2, the [Ne V] ratio combined with the [Ne VI]λ7.7μm /[O IV]λ25.9μm ratio, is unaffected by CRs suggesting that high-ionization lines are primarily influenced by photoionization.

All Tables

All Figures

Fig. 1. Apertures chosen to extract spectra from the spectroscopic images of NGC 5728. These are depicted in different shades of purple from deep purple being the nuclear aperture “N” to pale lilac in increasing distance, also noted with numbers. The apertures are drawn over (a) the continuum in Channel 1 of MIRI, extracted in the rest-frame range of 5.90–6.05μm and over [Ar II]λ7 μm, [S IV]λ10.5 μm, [Ne II]λ12.8 μm, [Ne V]λ14.9 μm, [S III]λ18.7 μm, Hα, and [O III]λ5007Å (b–h) line emission maps with the stellar continuum subtracted. Panel (i) presents a schematic illustration of the adopted jet axis together with the three color-shaded regions used to group the apertures (see Sect. 2.3). These regions are defined by the projected vertical distance from the jet axis, relative to the aperture radius: JR (innermost, cyan) includes apertures within one radius (nuclear and/or jet regions); IJR (intermediate, lime green) covers one to two radii (intermediate jet-affected regions); and NJR (outermost, magenta) includes apertures beyond two radii (non-jet or star-forming regions). The multiwavelength data are aligned based on the astrometry of the nucleus.

In the text

Fig. 2. BPT diagrams depicting [N II]/Hα, [S II]/Hα, and [O I]/Hα ratios. The observations from NGC 5728, and NGC 1320 are marked with stars, and a red diamond, respectively. The cyan squares, lime green diamonds, and magenta stars represent jet-affected (JR), intermediate jet-affected (IJR), and non-jet-affected (NJR) regions, respectively, and correspond to the shaded areas in Fig. 1. The Kewley and Schawinski lines are indicated with solid and dashed black lines, respectively, while the Koutsoumpou (SFζ) line, is depicted by the solid red line. In the background we show the line ratios measured for nearby galaxies from the Sloan Digital Sky Survey (SDSS) data release 7 (Abazajian et al. 2009).

In the text

Fig. 3. Diagrams with the AGN photoionization models compared with the observations from the selected apertures in NGC 5728 (Fig. 1). The different shades of purple ranging from deep purple to pale lilac and/or white, represent increasing distance from the nucleus, as also denoted by numbers, with “N” being the nuclear aperture. The different shapes-square, thin diamond, and star-represent the nucleus and/or jet impacted, intermediate and distant areas, respectively. The different shades from white to deep red represent the range of ionization parameter values, −3.5 ≤ log U ≤ −1.5. All the models have solar abundances. The panels from top to bottom correspond to ζCR = 10−16 s−1,  10−15 s−1,  10−14 s−1,  10−13 s−1, and 10−12 s−1, respectively.

In the text

Fig. 4. Temperature and line emissivity vs. depth in the simulated cloud for AGN models, for an initial density of nH = 100 cm−3, and for ζCR = 10−16 s−1, 10−15 s−1,  10−14 s−1,  10−13 s−1, and 10−12 s−1, and log U = −3.0. The different panels a-f correspond to kinetic temperature and the emissivity of [S III]18.7μm, [S IV]10.5μm,[Ar II]7μm,[Ar III]8.9 μm,[Ar V]7.9μm,[Ne II]12.8μm,[Ne III]15.5μm, and [Ne V]14.3μm, respectively. The shaded area indicates the approximate region in which CR heating becomes dominant. The CR-dominated area is depicted in lilac in the temperature plot and in teal in the emissivity plots.

In the text

	Fig. 5. Ionic fraction vs. depth in the simulated cloud for AGN models, for an initial density of nH = 100 cm−3, and for ζCR = 10−16 s−1, 10−15 s−1,  10−14 s−1,  10−13 s−1, and 10−12 s−1, from left to right. Top row: Ar ionic fractions. Bottom row: Ne ionic fractions. The blue-shaded area indicates the approximate region where CR heating becomes dominant.
In the text

Fig. 6. Diagnostic diagrams depicting the area covered by AGN models with solar abundances for −3.5 ≤ log U ≤ −1.5 and for 0 ≤ log nH ≤ 4. The solid contour lines map regions containing 10%, 50%, and 90% of the models for ζCR = 10−16 s−1, 10−15 s−1,  10−14 s−1,  10−13 s−1, and 10−12 s−1, from left to right in green-, teal-, blue-, purple-, and red-colored contours, respectively. In the background are archival Spitzer/IRS measurements of Seyfert and LINER nuclei (Dudik et al. 2007, 2009; Goulding & Alexander 2009; Tommasin et al. 2008, 2010; Pereira-Santaella et al. 2010; Fernández-Ontiveros et al. 2016, 2021), along with recent nuclear line ratios from JWST/MIRI observations from the literature (Pereira-Santaella et al. 2022; Álvarez-Márquez et al. 2023; Armus et al. 2023; Alonso Herrero et al. 2024; Zhang et al. 2024; Hernandez et al. 2025; Goold et al. 2024; Hermosa Muñoz et al. 2025, Goold et al. in prep., Fernández-Ontiveros et al. in prep.)

In the text

	Fig. B.1. AGN SEDs simulated with CLOUDY. The blue solid line represents the AGN SED used in KFDS25 while the teal solid line represents the AGN used in this paper. The y-axis is νFν in arbitrary scaling.
In the text

Fig. C.1. BPT diagrams comparing AGN photoionization models with observations from selected apertures in NGC 5728 (Fig. 1). The BPT diagrams for [N II], [S II], and [O I] are shown on the left, middle, and right, respectively. The different shades of purple, ranging from deep purple to pale lilac and/or white, represent the increasing distance from the nucleus, as also indicated by numbers, with "N" denoting the nuclear aperture. The different shapes-square, thin diamond, and star-represent the nucleus and/or jet impacted, intermediate, and distant areas, respectively. The different shades from white to deep red represent the range of ionization parameter values, −3.5 ≤ log U ≤ −1.5. The red diamonds represent the measured line ratios for the photoionization-dominated Seyfert 2 nucleus in NGC 1320. The Kewley, Kauffmann, Schawinski, and Koutsoumpou (SFζ) lines correspond to the solid black, blue, red, and magenta lines, respectively. The panels from top to bottom correspond to ζCR = 10−14 s−1, 10−13 s−1, and 10−12 s−1, respectively.

In the text

	Fig. D.1. The [Ar V]13.1μm/[Ar V]7.9μm and [Ne V]24.3μm/[Ne V]14.3μm line ratios produced with PYNEB (Luridiana et al. 2014) as a function of electron density, ne, notated for different electron temperatures, Te.
In the text

Fig. D.2. Diagrams with the AGN photoionization models compared with the observations from the selected apertures in NGC 5728 (Fig.1). The different shades of purple, ranging from deep purple to pale lilac and/or white represent the increasing distance from the nucleus, as also denoted by numbers, with "N" corresponding to the nuclear aperture. The different shapes-square, thin diamond, and star, represent the nucleus and/or jet impacted, intermediate, and distant areas, respectively. The different shades from white to deep red represent the range of ionization parameter values, −3.5 ≤ log U ≤ −1.5. All the models have solar abundances. The panels from top to bottom correspond to ζCR = 10−16 s−1, 10−15, s−1, 10−14 s−1, 10−13 s−1, and 10−12 s−1, respectively.

In the text