Cosmic-ray ionization of low-excitation lines in active galactic nuclei and starburst galaxies

E. Koutsoumpou; J. A. Fernández-Ontiveros; K. M. Dasyra; L. Spinoglio

doi:10.1051/0004-6361/202452232

Home

All issues

Volume 693 (January 2025)

A&A, 693 (2025) A215

Full HTML

Open Access

Issue		A&A Volume 693, January 2025


Article Number		A215
Number of page(s)		23
Section		Extragalactic astronomy
DOI		https://doi.org/10.1051/0004-6361/202452232
Published online		20 January 2025

A&A, 693, A215 (2025)

Cosmic-ray ionization of low-excitation lines in active galactic nuclei and starburst galaxies

E. Koutsoumpou¹^⋆, J. A. Fernández-Ontiveros², K. M. Dasyra¹ and L. Spinoglio³

¹ Section of Astrophysics, Astronomy & Mechanics, Department of Physics, National and Kapodistrian University of Athens, University Campus Zografos, GR, 15784 Athens, Greece
² Centro de Estudios de Física del Cosmos de Aragón (CEFCA), Plaza San Juan 1, E–44001 Teruel, Spain
³ Istituto di Astrofisica e Planetologia Spaziali (INAF–IAPS), Via Fosso del Cavaliere 100, I–00133 Roma, Italy

^⋆ Corresponding author; evkoutso@phys.uoa.gr

Received: 13 September 2024
Accepted: 25 November 2024

Abstract

Cosmic rays (CRs) can significantly impact dense molecular clouds in galaxies, heating the interstellar medium (ISM) and altering its chemistry, ionization, and thermal properties. Their influence is particularly relevant in environments with high CR rates, such as starburst galaxies with supernova remnants or jets and outflows in active galactic nuclei (AGN). CRs also transfer substantial energy to the ionized phase of the ISM far from the ionization source, preventing gas cooling and driving large-scale winds. In this work, we use CLOUDY photoionization models to investigate the effect of CRs on nebular gas which is an area of study that remains relatively under-explored, mainly focusing on cold molecular gas. Our models cover a broad range of density (1 to 10⁴ cm⁻³), ionization parameter (−3.5 ≤ log U ≤ −1.5), and CR ionization rate (10⁻¹⁶ s⁻¹ to 10⁻¹² s⁻¹). These are compared to VLT/MUSE observations of two prototypical AGN, Centaurus A (radio-loud) and NGC 1068 (radio-quiet), and the starburst NGC 253. We find that high CR rates (≳10⁻¹³ s⁻¹) typical of AGN and strong starburst galaxies can significantly alter the thermal structure of the ionized gas by forming a deep secondary low-ionization layer beyond the photoionization-dominated region. This enhances emission from low-ionization transitions, such as [N II]λ6584 Å, [S II]λλ6716,6731 Å, and [O I]λ6300 Å, affecting classical line-ratio diagnostics, metallicity, and ionization estimates. Unlike pure photoionization models, AGN simulations with high CR ionization rates reproduce the Seyfert loci in Baldwin, Phillips, and Terlevich (BPT) diagrams without requiring supersolar metallicities for the narrow-line region. Additionally, star-formation simulations with high CR ionization rates can explain line ratios in the LINER domain. We propose new maximum starburst boundaries for BPT diagrams in order to distinguish regions dominated by AGN photoionization from those that could be explained by star formation in conjunction with high CR ionization rates.

Key words: ISM: clouds / cosmic rays / galaxies: active / galaxies: jets / galaxies: starburst

© The Authors 2025

Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

This article is published in open access under the Subscribe to Open model. Subscribe to A&A to support open access publication.

1. Introduction

Several mechanisms can contribute to gas ionization in galaxies. First and foremost, there is photoionization caused by ultraviolet (UV) and X-ray emission from star-forming regions (e.g. Evans & Dopita 1985; McKee 1989) and supermassive black holes in active galactic nuclei (AGN; e.g. Groves et al. 2004a,b; Gallimore et al. 2004; Wolfire et al. 2022). Additionally, AGN- and starburst-driven outflows produce shock waves through jets and winds that heat and compress the gas (Bicknell et al. 2000; Fragile et al. 2004; Sutherland & Dopita 2017), resulting in excitation, ionization, and widening of the spectral line profiles. Diagnostic tools, such as the Baldwin, Phillips, and Terlevich diagrams (BPT; Baldwin et al. 1981), have been developed to analyze the gas excitation mechanisms in emission-line galaxies. The BPT diagram is built using the following optical line ratios: [O III]λ5007 Å/Hβ, [N II]λ6584 Å/Hα, [S II]λλ6716,6731 Å/Hα, and [O I]λ6300 Å/Hα. These diagrams have since been complemented with the theoretical “maximum starburst limit” according to photoionization models (Kewley et al. 2001) in order to identify AGN-dominated galaxies, an empirical separation between star-forming galaxies and Seyfert-H II composite nuclei (Kauffmann et al. 2003), and an empirical division between Seyfert and low-ionization nuclear emission-line region (LINER) galaxies (Kewley et al. 2006; Schawinski et al. 2007).

BPT diagrams are crucial in deciphering the complex physical processes of ionized nebulae and have been used to make detailed comparisons between observations and predictions from shock and photoionization models (Evans & Dopita 1985; Dopita & Evans 1986; Dopita & Sutherland 1995). In spite of that, to reconcile the distribution of optical line ratios observed for AGN galaxies with photoionization models in BPT diagrams, supersolar metallicities (2 − 4 Z_⊙) have been usually invoked in the past (e.g. Groves et al. 2004a,b). Other studies have shown that the observed line ratios could, at least partially, be explained by slightly supersolar metallicities (1.6 − 2.6 Z_⊙), along with steeper UV slopes in the ionizing continuum (Feltre et al. 2016), and a hotter accretion disk continuum, or a higher N/O abundance (Zhu et al. 2023). However, subsolar to solar abundances are obtained for the narrow-line region (NLR) in nearby AGN using photoionization models when additional optical line ratios, beyond those involved in the BPT, are considered (Pérez-Montero et al. 2019; Pérez-Díaz et al. 2021), or when infrared lines are included (Pérez-Díaz et al. 2022). These estimates also agree with semi-empirical calibrations of UV lines (Dors et al. 2019).

This discrepancy between models and observations can be caused by uncertainties in the ionizing continuum (Zhu et al. 2023) and/or additional sources of excitation, such as shocks (Allen et al. 2008). In this regard, a more elusive source of gas excitation, whose influence on the BPT line ratios has not yet been thoroughly studied, are cosmic rays (CRs). These are highly energetic particles mostly produced by supernova remnants (SNRs) and accreting black holes (Blasi 2013; Padovani et al. 2017; Veilleux et al. 2020; Wolfire et al. 2022; Kantzas et al. 2023) able to penetrate deeply into massive molecular clouds and deposit their energy at much greater depth than photons (McKee 1989; Padovani et al. 2018), determining the chemistry and physics of the cold interstellar medium (ISM; Dalgarno 2006). As they traverse through the ISM, the infiltration of low-energy CRs (≲1 GeV) is capable of producing excitation and ionization from secondary electrons due to collisions with atoms and molecules deep within the cloud (Spitzer & Tomasko 1968; Gredel et al. 1989; Gabici 2022). Thus, CRs are expected to have a direct influence on nebular emission lines (Ferland & Mushotzky 1984; Walker 2016), especially in galaxies with strong star formation activity (Phan et al. 2024) or close to AGN jets where a dense flow of these particles is carried along (Pacholczyk 1970; Guo & Mathews 2011).

While other excitation mechanisms like UV and X-ray photoionization and shocks have attracted most of the attention in previous studies of ionized gas (e.g. Allen et al. 2008; Ferland et al. 2017; Chatzikos et al. 2023), the more difficult detection and characterization of CRs in galaxies has prevented a more detailed analysis of the impact that these particles have on the physical conditions and the stratification of the nebular gas. Nevertheless, understanding the effects that CRs have on the ISM of galaxies is essential in order to characterize their net contribution to feedback processes and therefore to galaxy evolution (Ruszkowski et al. 2017; Buck et al. 2020; Peschken et al. 2023; Lin et al. 2023).

In the present work, we use photoionization simulations to investigate the impact of CRs on the ionization of the nebular gas in different AGN and star-formation scenarios, covering a large parameter space in terms of CR ionization rate, gas density, and intensity of the ionizing radiation. For this purpose, three main prototypical galaxies were selected in order to compare the effects of CRs in different environments where they are expected to have a relevant contribution: an AGN with prominent radio jets, Centaurus A, an AGN with star-forming activity, NGC 1068, and the nearby starburst galaxy, NGC 253. Additionally, the photoionization-dominated NLR in NGC 1320 is included for comparison. The main properties and related data of the three selected galaxies are presented in Section 2. We detail our models and the procedure for examining specific regions of interest on each galaxy in Section 3. In Section 4, we present the main results of this study, and in Section 5, we discuss the implications and possible future directions. Finally, our conclusions are summarized in Section 6.

2. Observational data

2.1. Galaxy sample

The purpose of this study is to examine whether CRs are capable of producing the observed emission line ratios. For this analysis, we have carefully selected galaxies with publicly available VLT/MUSE observations in the European Southern Observatory (ESO) Science Archive Facility (see Section 2.2), with distinct characteristics related to CRs, originating either from strong radio-loud jets or SNRs associated with intense star formation. To effectively examine the influence of CRs, it is crucial to choose galaxies that are representative of these phenomena.

To strengthen our case, we examine Centaurus A, a well-studied jetted radio-loud AGN with a prominent jet (Hardcastle et al. 2003; Kraft et al. 2008; Tingay & Lenc 2009). This is also the nearest jetted AGN, making it an ideal subject to study jet-induced CRs. We include in our study NGC 253, the nearest starburst galaxy which as of recently has direct CR rate measurements (Behrens et al. 2022; Holdship et al. 2022; Beck et al. 2023), most likely associated with the CR component from SNRs. NGC 1068 has both an AGN and a starburst component whose composite photoionized spectrum from the UV to the far-IR has been successfully modeled (Spinoglio et al. 2005). It serves as an intermediate example hosting both jets (Bland-Hawthorn et al. 1997; Mutie et al. 2024) and star-forming regions (García-Burillo et al. 2016), which makes it a particularly informative target for observing the effects of both jet and supernova related CRs. The effects of CRs in the far-IR molecular lines of H₂O⁺ and OH⁺ have been reported in the analysis of its submillimeter line spectrum (Spinoglio et al. 2012).

Furthermore, as a control source, NGC 1320, is used to counterbalance these cases. It stands out from the other galaxies, which are mostly influenced by CRs, because of its low CR-related behavior and characterized by a highly photoionized gas. For each of the aforementioned galaxies additional details are given in Table 1 and additional relevant literature is provided in the Appendix A.

Table 1.

Observational parameters for the selected galaxies, including coordinates, redshift values, z, and redshift-independent distances, D, provided by the NASA/IPAC Extragalactic Database.

These galaxies offer an excellent basis to assess the effects of CRs, considering that the expected CR presence is significant in these systems. All of these sources, lying in close proximity (see Table 1), have been observed by the Very Large Telescope (VLT) Multi Unit Spectroscopic Explorer (MUSE) providing high spatial resolution optical line imaging, offering high-quality data necessary for this study (see Table 2 for observational details). All in all, the selected galaxy sample guarantees a thorough examination in a variety of galactic environments.

Table 2.

Description of MUSE datacubes used.

2.2. MUSE data

Raw astronomical data have to be carefully processed to obtain verified and reliable products suitable for the final analysis. This is done via specific software, the pipeline, which involves several steps. VLT-MUSE data are processed by the ESO MUSE pipeline which performs a bias correction to remove camera electronic offsets, dark subtraction to compensate for the thermal noise of the detector, flat fielding to balance variability in detector sensitivity and illumination and background subtraction to remove ambient light and improve object visibility. The pipeline also includes flux calibration, which converts counts to standard flux units for comparability across different observations and wavelength calibration that ensures wavelength precision vital for spectral analysis. So as to get rid of artifacts that could impair image clarity, the removal of glitches due to the impact of CRs in the detector systems is also conducted. Finally, spectral extraction isolates the object’s spectra from 2D images, allowing detailed scientific analysis. The MUSE data processing procedures in different operational stages are detailed in Weilbacher et al. (2014, 2020).

In our analysis we make use of MUSE datacubes of Centaurus A, NGC 1068, NGC 253, and NGC 1320 which are publicly available in the ESO archive. For more specific information on the data reduction done on the galaxies of this work, see Mingozzi et al. (2019), Venturi et al. (2021) and Humire et al. (2022) for Centaurus A, NGC 1068 and NGC 253, respectively. An overview of the datasets utilized in this paper is provided in Table 2.

3. Methods

3.1. Spectrum extraction

The primary goal of this work is to explore the parameters that influence the excitation of gas across the various regions of galaxies, with a particular focus on the effects of CRs along with photoionization. Our method relies on selecting specific regions within the aforementioned galaxies and extracting the spectra of the MUSE datacubes from circular apertures. This technique allows an in depth examination of gas excitation within different galactic environments by use of the spatially resolved spectral information provided by the datacube of each galaxy. Each and every aperture is positioned strategically in order to understand the influence of CRs in star-forming regions and areas affected by the jets in the AGN cases.

We predict that areas in proximity of the AGN jets are mainly impacted by the CRs related to these phenomena. As opposite, regions situated near star-forming regions and farther from the galactic center are presumed to be dominated by photoionization from UV radiation produced by hot young stars. We intend to distinguish the corresponding contributions of photoionization and CRs in various regions of each galaxy by collecting and analyzing the spectra associated within the different apertures selected.

In the cases of Centaurus A and NGC 1068, we select apertures containing gas that interacts with the jets/outflows, as well as potential sites for star formation. In NGC 253 we inspect regions close to the central area of the galaxy, some of which are affected by strong stellar winds and others are potential sites of star formation. To achieve this, we identify regions that exhibit strong emission in both [O III]λ5007 Å and Hα along with the other emission lines required to create the BPT diagrams, namely [N II]λ6584 Å, [S II]λλ6716, 6731 Å and [O I]λ6300 Å. The BPT diagrams and their applications are introduced in Section 1 and fully detailed in Section 4.

In particular, for our analysis we implemented 18 apertures of a ∼50 pc radius in the vicinity of Centaurus A, 15 of the same radius in the vicinity of NGC 253, and 20 apertures of a ∼80 pc radius in the vicinity of NGC 1068, as the MUSE datacube field of view of NGC 1068 in kpc is the largest among the three galaxies considered (see Table 2). From all the apertures of each galaxy we extracted the convolved spectrum using MPDAF (Bacon et al. 2016). We also chose one aperture of ∼200 pc radius from the central area of NGC 1320, to establish the position of highly photoionized gas on the BPT diagrams to compare with the rest of our models and observations.

We present the selected apertures overplotted on the Hα and [O III]λ5007 Å emission line maps in Fig. 1. Figs. 1b and 1e show the apertures in NGC 1068, Figs. 1a and 1d correspond to Centaurus A, while Figs. 1c and 1f correspond to NGC 253. We note that NGC 253 lacks of the “N” aperture, since this is not an AGN. In this case, the apertures are numbered from 1 to 15, with 1 being the closest to the central region aperture, as shown in Figs. 1c and 1f. In Centaurus A and NGC 1068, the jet direction is very well defined by radio continuum maps (Dufour & van den Bergh 1978; McKinley et al. 2018; Joseph et al. 2022; Gallimore et al. 1996, 2004). In Figs. 1a and 1d we show the VLA map at 8.4 GHz provided by Tingay & Lenc (2009). In the case of NGC 1068, Figs. 1b and 1e show the combined 5 GHz continuum map from VLA and e-Merlin, provided by Mutie et al. (2024).

Fig. 1.

Apertures chosen to extract spectra with MPDAF. The different shades of purple going from deep purple to pale lilac represent the ascending distance, and are also noted with numbers. The apertures are drawn over Hα and [O III]λ5007 Å, emission with the stellar continuum subtracted. In Centaurus A, the jets are depicted by use of VLA data at 8.4 GHz (Hardcastle et al. 2003; Tingay & Lenc 2009), in white contours at levels 7 × 10⁻⁴ Jy/beam, 10⁻² Jy/beam, 0.1 Jy/beam, and 6 Jy/beam. Similarly, for NGC 1068, the jets are illustrated using combined data from e-MERLIN and the VLA at 5 GHz (Gallimore et al. 1996; Muxlow et al. 1996; Mutie et al. 2024), in white contours at levels 10⁻⁴ Jy/beam, 5 × 10⁻² Jy/beam, and 10⁻² Jy/beam. All radio data were aligned with optical data based on astrometry. (a) Centaurus A, H_α line map. (b) NGC 1068, H_α line map. (c) NGC 253, H_α line map. (d) Centaurus A, [O III]λ5007 Å line map. (e) NGC 1068, [O III]λ5007 Å line map. (f) NGC 253, [O III]λ5007 Å line map.

The emission lines detected in the convolved spectra are fitted with PYPLATEFIT (Bacon et al. 2023), a Python package that fits the emission and absorption lines of MUSE spectra, inspired by the IDL routine PLATEFIT (Tremonti et al. 2004; Brinchmann et al. 2004). In this package, continuum subtraction is performed by fitting and subtracting simple stellar population models (Bruzual & Charlot 2003; Brinchmann et al. 2013). Subsequently, the key emission lines utilized in the BPT diagrams are fitted in each galaxy’s rest frame wavelengths. We present the respective fit of the BPT lines coming from Aperture 2 in Centaurus A as shown in Fig. 2. The fit of Hβ for the 3 galaxies is also presented in Fig. 3. As Hβ is a highly absorbed line and from the fits presented, it is clear that there is attenuation that needs to be addressed (see the discussion in Section 3.3).

Fig. 2.

BPT emission-line fits, using Pyplatefit (Bacon et al. 2023), in the rest frame of Centaurus A. (a) Hα. (b) [NII]λ6584 Å. (c) [S II]λλ6716,6731 Å. (d) [O I]λ6300 Å. (e) [O III]λ5007 Å.

Fig. 3.

Hβ emission line in each galaxy’s rest frame, fitted with Pyplatefit (Bacon et al. 2023). (a) Centaurus A. (b) NGC 1068. (c) NGC 253.

In particular, when examining the central region of NGC 1068, we found broad-line components associated with ∼100–200 km s⁻¹ outflows known in the vicinity of the nucleus (García-Burillo et al. 2016). In order to address this issue in the nuclear aperture N, we applied a multicomponent fitting using the LMFIT package (Newville et al. 2015). In this case, the continuum was subtracted after fitting a Chebyshev polynomial with SPECUTILS (Earl et al. 2023), which additionally smooths the spectrum using a median filter prior to the fit, to mitigate the effects of noise and remove significant stellar and AGN contributions. We fitted both narrow and wide emission line components with Gaussian functions. This procedure enabled us to acquire reliable measurements of the flux for every emission-line component, and specifically of the emission-line fluxes used in BPT diagrams that are necessary for the study of the central region in NGC 1068.

3.2. Constraining the parameter space

Given that the parameter space is fairly large, we have to restrict the main parameters before simulating the cloud physical conditions. One of the main parameters under study in the present work is the ionization parameter, U defined as the ratio of ionizing photons to the gas density and expressed as (Osterbrock & Ferland 2006):

$\begin{matrix} U = \frac{Q}{4 π r^{2} n_{H} c}, \end{matrix}$ $\begin{aligned} U = \frac{Q}{4 \pi r^2 n_{\rm H} c}, \end{aligned}$ (1)

where Q is the photon emission rate, r is the distance to the source, n_H is the gas density, and c is the speed of light. The other main parameters relevant to this study are the initial hydrogen density n_H, the metallicity Z, and the CR ionization rate ζ_CR.

As mentioned in Section 1, studying the excitation mechanisms is a strenuous task especially when having in mind all these factors that can affect the emission lines. Constraining the ionization parameter is necessary for grasping the roles that photoionization and its competitive excitation gas process, that of CRs, play in our models. With respect to n_H, our models cover a wide range of values in order to cover a variety of regions with different densities. It is, moreover, important to constrain the metallicity and have realistic values of it, as it has been a usual practice to use unrealistic supersolar metallicities to effectively model observations of Seyfert galaxies or NLR (Groves et al. 2004b).

Metallicity is an essential parameter in understanding galaxy evolution through the enrichment of the ISM. Metals are created in the cores of stars, transported by convection, and dispersed by outflows and supernovae (SNe). The abundance ratios of oxygen or nitrogen relative to hydrogen are considered indicators of metallicity. Since oxygen is abundant, as can be detected by prominent optical emission lines, the most popular choice to measure metallicity is using: 12 + log(O/H) (Asplund et al. 2009; Maiolino & Mannucci 2019). The N/O ratio also serves to differentiate between metals coming from helium burning in stellar cores and already existent in the ISM components (van Zee et al. 1998).

We are able to constrain the ionization parameter, O/H, and N/O ratios for each of the regions measured in the sample of galaxies, by employing HII-CHI-MISTRY¹, developed by Pérez-Montero (2014), uses Bayesian-like statistics to estimate these parameters by comparing emission-line ratios with a large grid of photoionization models. The code has been adapted to use metallicity tracers in different wavelength ranges; in particular, we employed the optical versions for star-forming galaxies (Pérez-Montero 2014) and Seyfert nuclei (Pérez-Montero et al. 2019).

In the optical wavelengths, HII-CHI-MISTRY employs the following emission lines, [O II]λ3727 Å, [Ne III]λ3868 Å, [O III]λ4363 Å, [O III]λ4959 Å, [O III]λ5007 Å, [N II]λ6584 Å, and [S II]λλ6716, 6731 Åto determine the chemical composition and physical conditions of ionized gas. However, [O II]λ3727 Å and [Ne III]λ3868 Å are not within the MUSE spectral range, which covers 4800 to 9300 Å. Therefore, the MUSE datacubes used (Section 2), limit the lines we use in the HII-CHI-MISTRY routines to [O III]λ4959 Å, [O III]λ5007 Å, [N II]λ6584 Å, and [S II]λλ6716, 6731 Å. These estimates result in metallicities from slightly subsolar to ∼1.2 Z_⊙ in accordance to metallicities found in already published literature (Pérez-Díaz et al. 2021) making the hypothesis of 1 Z_⊙ realistic. The N/O abundance ratio is also a valuable diagnostic tool to analyze galactic processes. In our dataset, we derive relative abundances in the range of −1.5 < log(N/O) < − 0.8 (Table 3), which are very close to the solar value of log(N/O)_⊙ ∼ −0.86 dex (Asplund et al. 2009).

Table 3.

Range of HII-CHI-MISTRY estimated values for the ionization parameter, the oxygen abundance, and the N/O relative abundance of the extracted regions, including the physical size of the aperture.

3.3. CLOUDY models

We model the line fluxes with the photoionization and radiative transfer code CLOUDY (Ferland et al. 2013, 2017) and make also use of the pyCLOUDY (Morisset 2013). We examine as possible excitation mechanisms photons, in a wide spectral range, depending on the ionizing spectrum, and CRs. With CLOUDY we examine CRs in the sense of the CR ionization rate ζ_CR integrating an electron CR component in the energy range 5 MeV − 10 GeV. CLOUDY does not include extinction and to compare the results of our simulations with observational data we transform the CLOUDY modeled fluxes to attenuated fluxes using Calzetti’s law (Calzetti 2001).

Initially, we develop a grid of CLOUDY models for five discrete CR ionization rates 10⁻¹⁶ s⁻¹, 10⁻¹⁵ s⁻¹, 10⁻¹⁴ s⁻¹, 10⁻¹³ s⁻¹, and 10⁻¹² s⁻¹ for hydrogen densities in the 1 ≤ n_H ≤ 10⁴ cm⁻³ range, with a step of 10^0.1 cm⁻³ and an ionization parameter −3.5 ≤ log U ≤ −1.5 with a step of 0.1 dex. Regarding the ionization parameter U values, we chose this range motivated by the estimations derived with HII-CHI-MISTRY (see Table 3), in agreement with values previously reported in the literature for AGN and star-forming regions (Pérez-Montero & Vílchez 2009; Pérez-Montero et al. 2019; Carvalho et al. 2020). The model grids have a solar chemical composition (Asplund et al. 2009). However, the models with CR ionization rates 10⁻¹⁵ s⁻¹ and 10⁻¹⁴ s⁻¹ do not show substantial differences; therefore, we present only the 10⁻¹⁴ s⁻¹ case. Moreover, the effect of CRs in the case of 10⁻¹⁶ s⁻¹ is small and will be presented for comparison in Sect. 5.3. The CR ionization rates values used in our models is supported by indirect measurements for both AGN and star-forming galaxies reported in the literature (González-Alfonso et al. 2013, 2018; Holdship et al. 2022; Behrens et al. 2022). A further justification of these CR values can be found in Section 5.5.

We employ two different models as input for the CLOUDY simulations. AGN models are used for the cases of Centaurus A and NGC 1068, while star-forming models are used for NGC 253. The star-forming models, generated with the BLACKBODY command, provide less energetic photons, as expected for a thermal continuum from stellar populations (Ferland et al. 2013). For simplicity, and to have better control over the shape of the incident continuum through the temperature, we adopted a black body (e.g. Pérez-Montero et al. 2024) instead of more sophisticated stellar population models. Nevertheless, we tested our predictions using Starburst99 models (Leitherer et al. 1999), and found that the differences in the BPT line ratios have a median relative difference of ∼7%, indicating a marginal deviation. The AGN models, on the other hand, are simulated using CLOUDY’s AGN command and assume a ionizing continuum typical of AGN (Ferland et al. 2013). AGN models can effectively reproduce the intense radiation and ionizing processes in the two jetted galaxies. The shapes of both the AGN and SF ionizing continua are shown in Fig. 4, with arbitrary scaling for νF_ν.

Fig. 4.

AGN and BB SEDs simulated with CLOUDY. The blue solid line represents the AGN SED used in the AGN models, while the green solid line represents the BB SED used in the SF models. The y-axis is νF_ν in arbitrary scaling.

More precisely, in the grid of models used for Centaurus A and NGC 1068, characterized as AGN models for the rest of the paper, the spectral energy distribution (SED) of the central photoionizing continuum is assumed to be consisting 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 may be expressed mathematically as the following expression:

$\begin{matrix} F_{ν} = ν^{α_{uv}} exp (- \frac{h ν}{k T_{bb}}) exp (- \frac{k T_{IR}}{h ν}) + α_{ox} ν^{α_{x}}, \end{matrix}$ $\begin{aligned} F_\nu = \nu ^{\alpha _{\text{ uv}}} \exp \left(-\frac{h\nu }{kT_{\text{ bb}}}\right)\exp \left(-\frac{kT_{\text{ IR}}}{h\nu }\right) + \alpha _{\text{ ox}}\nu ^{\alpha _{\text{ x}}}, \end{aligned}$ (2)

which is typical for AGN. In this formula F_ν is the flux density as a function of frequency ν, α_uv and α_x are the UV and X-ray spectral indexes respectively, T_bb is the representative temperature of the big blue bump, T_IR the temperature related to the infrared cutoff for the big blue bump, and α_ox describes the optical to X-ray spectral index. As a SED for our AGN models, we employed the following values, as they are both typical for AGN and are in accordance with the observational data of the two cases we study. We used α_uv = −0.5, α_x = −1.0, α_ox = −1.4, (T_BB = 10⁵) K, and (T_IR = 1.6 × 10³) K.

In the case of NGC 253, we simulate star-forming (hereafter SF) models using a black body (BB) with a characteristic temperature of T_BB = 5 × 10⁴ K as ionizing continuum. Such a high temperature is representative of the hot, massive O-type stars that dominate the emission in massive young star clusters (∼10⁵ M_⊙, ≲6 Myr) as those found in the nuclear starburst of this galaxy (Watson et al. 1996; Fernández-Ontiveros et al. 2009; Mills et al. 2021).

In all cases, we adopted solar elemental abundances, including 12 + log(O/H)_⊙ = 8.69, log(N/O)_⊙ = −0.86, log(C/O)_⊙ = −0.3, and log(He/H)_⊙ = 0.085 (Asplund et al. 2009). These values are consistent with HII-CHI-MISTRY abundances obtained for the different regions in our sample (Table 3) and with estimates obtained for the NLR of AGN in the local Universe from previous studies (e.g. Pérez-Díaz et al. 2021, 2022, 2024; Dors 2021; Dors et al. 2022). On the other hand, the presence of dust in NLR clouds has been a topic of debate, with studies both supporting (e.g. Groves et al. 2004a,b; Zhu et al. 2023) and opposing it (e.g. Ferguson et al. 1997; Nagao et al. 2006; Richardson et al. 2014). Nevertheless, Feltre et al. (2016) show that dust depletion in models with 1.5 Z_⊙ has a small impact on BPT line ratios, while differences between dusty and dust-free models in Zhu et al. (2023) are significant primarily at supersolar metallicities. We opted for dust-free models in our simulations, which are suggested to be more appropriate for radio galaxies (Matsuoka et al. 2009), where the effects of CR ionisation are expected to be more important. The presence of dust in NLR clouds is expected to increase the electron temperature through photoelectric heating and a lower cooling efficiency caused by the depletion of metals from the gas phase. Predictions for dusty models will be explored in a future study discussing the effect of CRs on mid-infrared transitions.

4. Results

4.1. BPT diagrams

Our study intends to bridge the knowledge gap and offer an alternative method to model the observational data shown in the upper AGN area on the BPT diagrams. Our primary objective is to better understand the physical mechanisms that underlie ionized nebulae and their implications on BPT plots. We compare the location of the extracted spectra for the different regions on the BPT diagrams with the predictions obtained from the photoionization along with CRs models. We present the BPT diagrams for all three galaxies in the local universe that are examined in this study: Centaurus A, NGC 1068, and NGC 253.

The observational data of Centaurus A and NGC 1068 are compared to the AGN models, while those of NGC 253 are compared to the SF models; all models are presented in Section 3.3. The emission line ratios of all galaxies create data points on the BPT diagrams depicted as deep purple to white crosses that are located both in star-forming and composite regions, under the Kewley division line (Kewley et al. 2006), as well as above it, in the Seyfert/LINER area of the diagrams. The color of the crosses indicates their physical distance from either the nucleus of their host AGN or from the most central area in the starburst case, going from deep purple being the closest one to white being the most distant one and accompanied as well by the ascending aperture number. The apertures closest to the nucleus, depicted in darker shades of purple, tend to be in the Seyfert/LINER area of the BPT above the Kewley line, whereas the areas further away from the nuclei of both AGN, have emission line ratios that locate crosses in the star-forming area of the BPT diagrams (Figs. 5, 6). In NGC 253, we find that even though this is not an AGN, there are apertures with line ratios similar to those of AGN (Fig. 7), which means that the gas in the vicinity of NGC 253 has similar properties to gas in AGN. Additionally, there are several areas that exhibit star-forming characteristics in the two jetted AGN, Centaurus A and NGC 1068 (Figs. 5, 6).

Fig. 5.

BPT diagrams with the AGN photoionization models compared with the observations from the selected apertures in Centaurus A (Fig. 1a). 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 going from deep purple to pale lilac/white represent the ascending distance from the nucleus, as also noted with numbers, with “N” being the closest aperture. Also, from white to deep red, the different shades of red represent the range of ionization parameter values, −3.5 ≤ log U ≤ −1.5. All the models shown have solar abundances. The red diamonds represent the measured line ratios for the photoionization-dominated Seyfert 2 nucleus in NGC 1320. The Kewley, Kauffmann, and Schawinski lines correspond to the black, blue, and red solid lines, respectively. (a) ζ_CR = 10⁻¹⁴s⁻¹. (b) ζ_CR = 10⁻¹³s⁻¹. (c) ζ_CR = 10⁻¹²s⁻¹.

Fig. 6.

BPT diagrams with the AGN photoionization models compared with the observations from the selected apertures in NGC 1068 (Fig. 1b). 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 going from deep purple to pale lilac/white represent the ascending distance from the nucleus, as also noted with numbers, with “N” being the closest aperture. Also, from white to deep red, the different shades of red represent the range of ionization parameter values, −3.5 ≤ log U ≤ −1.5. All the models shown have solar abundances. The red diamonds represent the measured line ratios for the photoionization-dominated Seyfert 2 nucleus in NGC 1320. The Kewley, Kauffmann, and Schawinski lines correspond to the black, blue, and red solid lines, respectively. (a) ζ_CR = 10⁻¹⁴s⁻¹. (b) ζ_CR = 10⁻¹³s⁻¹. (c) ζ_CR = 10⁻¹²s⁻¹.

In Fig. 5 as well as in Fig. 6, we see the parameter space covered by our models. These models cover the ionization parameter range −3.5 ≤ log U ≤ −1.5 from white, being the lower ionization parameter, to dark red, being the highest (Figs. 5, 6, and 7). These models also span between 0 ≤ log n_H ≤ 4 from white, being the lower initial hydrogen density, to forest green, being the highest, as shown in Appendix B. These particular diagrams are produced for different CR ionization rates, from the top to the bottom in the figures: −14 ≤ log ζ_CR ≤ −12 for all galaxies; the instance of Centaurus A can be found in Figs. 5a, 5b, and 5c, respectively.

Fig. 7.

BPT diagrams with the SF photoionization models compared with the observed line ratios from the selected apertures in NGC 253 (Fig. 1c). 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 going from deep purple to pale lilac/white represent the ascending distance, as also noted with numbers, with “1” being the most central aperture. Also, from white to deep red, the different shades of red represent the range of ionization parameter values, −3.5 ≤ log U ≤ −1.5. All the models shown have solar abundances. The red diamonds represent the measured line ratios for the photoionization-dominated Seyfert 2 nucleus in NGC 1320. The Kewley, Kauffmann, and Schawinski lines correspond to the black, blue, and red solid lines, respectively. (a) ζ_CR = 10⁻¹⁴s⁻¹. (b) ζ_CR = 10⁻¹³s⁻¹. (c) ζ_CR = 10⁻¹²s⁻¹.

Not only the AGN models shown in Figs. 5, and 6 but the SF models as well, shown in Figs. 7, are relocated in the right area, the AGN area, of the BPT diagrams as the CR rate grows higher. These models mainly fall on the LINER part of the AGN region. As a general result (Figs. 5, 6, and 7), analyzing in detail the distribution of the models, it is clear that as the CR rate grows, there is more emission on [N II]λ6584 Å, [S II]λλ6716,6731 Å, and [O I]λ6300 Å. This extra emission is behind the relocation of models in the AGN locus of the BPT diagrams.

It is also evident that for each of the three aforementioned emission lines there is a density over which the emission of the line falls and then begins to increase again. This can be seen in Figs. B.1, B.2, and B.3 where the models seem to fold and change direction over a certain density. This is due to the effect of the CRs that can penetrate the clouds deep enough and ionize the gas in these regions but over a certain column density where the gas is also denser and more shielded they gradually lose their ability to excite.

The impact of CRs on ionized gas seems to be a plausible mechanism behind the emission line ratios in the AGN locus of the BPT diagrams in the absence of a strong radiation field. In the case of starburst nuclei high enough CR rates are able to model AGN-like emission in the vicinity of its host without assuming an AGN continuum. For CR ionization rate ∼10⁻¹³ s⁻¹ we acquire AGN and SF models capable of reproducing quite well the observations in the BPT with [S II]λλ6716,6731 Å, for all the cases studied (Figs. 5b, 6b, and 7b). Then, for 10⁻¹² s⁻¹ there are models in the far right area of the [N II]λ6584 Å BPT diagrams for both AGN and SF models (Figs. 5c, 6c, and 7c). Moreover, we find that the highest CR rates 10⁻¹² s⁻¹, fit the observed emission in the [N II]λ6584 Å BPT for both AGN and SF models in all galaxies (Figs. 5c, 6c, and 7c).

Finally, in the [O I]λ6300 Å BPT diagram, AGN models match the observations for 10⁻¹³ s⁻¹ in the case of Centaurus A (Fig. 5b), while the observed line ratios in NGC 1068 could be explained with 10⁻¹⁴ s⁻¹ (Fig. 6b). SF models also require 10⁻¹³ s⁻¹ to extend the simulated [O I]/Hα ratios beyond the star-forming limit and reproduce the observed ratios in NGC 253 (Fig. 7b). Most of the SF models with 10⁻¹⁴ s⁻¹ do not reach Seyfert/LINER loci in the BPT, which are occupied by some of the regions in NGC 253. At 10⁻¹² s⁻¹, both AGN and SF models overestimate the observed [O I]/Hα ratios. The ionization of low-excitation lines by CRs will be further discussed in Section 5.2.

It is also noteworthy that models with the highest CRs ionization rates, together with log U ≳ −2.0 and log n_H ≳ 1, drive both the AGN and SF models in the Seyfert/LINER area of the [N II]λ6584 Å BPT diagrams (Figs. 5c, 6c, 7c, 1c, 2c, and 3c). Overall, our photoionization models show that even the highest ionization parameter values, depicted with dark red circles, are not sufficient to reach the Seyfert/LINER domain of the BPT without the presence of a high enough CR ionization rate. This behavior is closely related to the sensitivity of low-ionization emission lines, such as [S II]λλ6716,6731 Å, [N II]λ6584 Å, and [O I]λ6300 Å to the low-energy CRs. The effects of CRs on these emission lines will be discussed in Section 4.2, while a comparison of these results with other studies using photoionization models is presented in Section 5.3.

Finally, in Centaurus A and NGC 1068 BPT diagrams (see Figs. 5, and 6), apertures located closer to the nucleus and the jets, in darker shades of purple, are mainly positioned in the LINER/Seyfert area and better reproduced by models of the higher CR ionization rates 10⁻¹³s⁻¹–10⁻¹²s⁻¹. This seems logical considering that due to their position, they are probably affected by the energetic particles originating from the jets, something that it is to be further discussed in Sect. 5.5.

Fig. 8.

Temperature vs. depth in the simulated cloud for AGN and SF models, for a given initial density (see subcaption), three ζ_CR values, 10⁻¹⁴ s⁻¹ (solid line), 10⁻¹³ s⁻¹ (dash-dotted), and 10⁻¹² s⁻¹ (dotted), and two log U values of −3.0 (green) and −2.0 (blue). The teal-shaded area indicates the approximate region where CR heating becomes dominant. (a) Temperature, AGN, n_H = 100 cm⁻³. (b) Temperature, AGN, n_H = 10³ cm⁻³. (c) Temperature, SF, n_H = 100 cm⁻³.

Fig. 9.

Line emissivities vs. depth for AGN and SF models, for three ζ_CR values of 10⁻¹⁴ s⁻¹ (solid line), 10⁻¹³ s⁻¹ (dash-dotted), and 10⁻¹² s⁻¹ (dotted), and two log U values of −3.0 (green) and −2.0 (blue). The teal-shaded area indicates the region where CRs heating becomes dominant. (a) [N II]λ6584 Å, AGN, n_H = 100 cm⁻³. (b) [N II]λ6584 Å, AGN, n_H = 10³ cm⁻³. (c) [N II]λ6584 Å, SF, n_H = 100 cm⁻³. (d) [S II]λλ6716,6731 Å, AGN, n_H = 100 cm⁻³. (e) [S II]λλ6716,6731 Å, AGN, n_H = 10³ cm⁻³. (f) [S II]λλ6716,6731 Å, SF, n_H = 100 cm⁻³. (g) [O I]λ6300 Å, AGN, n_H = 100 cm⁻³. (h) [O I]λ6300 Å, AGN, n_H = 10³ cm⁻³. (i) [O I]λ6300 Å, SF, n_H = 100 cm⁻³. (j) Hα, AGN, n_H = 100 cm⁻³. (k) Hα, AGN, n_H = 10³ cm⁻³. (l) Hα, SF, n_H = 100 cm⁻³.

Fig. 10.

Line emissivities vs. depth for AGN and SF models, for three ζ_CR values of 10⁻¹⁴ s⁻¹ (solid line), 10⁻¹³ s⁻¹ (dash-dotted), and 10⁻¹² s⁻¹ (dotted), and two log U values of −3.0 (green) and −2.0 (blue). The teal-shaded area indicates the region where CRs heating becomes dominant. (a) [O III]λ5007 Å, AGN, n_H = 100 cm⁻³. (b) [O III]λ5007 Å, AGN, n_H = 10³ cm⁻³. (c) [O III]λ5007 Å, SF, n_H = 100 cm⁻³. (d) Hβ, AGN, n_H = 100 cm⁻³. (e) Hβ, AGN, n_H = 10³ cm⁻³. (f) Hβ, SF, n_H = 100 cm⁻³.

4.2. Gas stratification diagrams

The ionization and excitation of the gas is most likely achieved with UV radiation coming from young hot stars (Kim et al. 2018). Nevertheless, it is interesting to note that UV photons of energies more than 13.6 eV, which is the ionization potential of hydrogen, are quickly absorbed and cannot pass through all the layers of the clouds (Gabici 2022). The UV radiation is significantly attenuated for column densities N_H ≥ 8 × 10²¹ cm⁻² and CRs become the prevailing ionizing mechanism (McKee 1989).

The effects of CRs are deeper in the clouds where examining the emissivity of emission lines explains why our models move to the right area on the BPT diagrams. In order to analyze this, we produce the following gas stratification diagrams (hereafter referred to as structure plots) for two ionization parameter values, that is, log U = { − 3.0, − 2.0} and for densities n_H = {100, 10³} cm⁻³ in the case of AGN models and the same for n_H = 100 cm⁻³ in the case of SF models. Fig. 8 makes quite evident the fact that the electron temperature at the illuminated face of the cloud is mainly affected by photoionization and the deeper layers show increase in temperature in the cases of the higher CR ionization rates 10⁻¹³ s⁻¹ and 10⁻¹² s⁻¹. The higher the ionization parameter is, the bigger is the depth up to which photons dominate. However, after a certain layer CRs are the main excitation mechanism.

Similarly when investigating simulations with higher initial hydrogen density n_H in the AGN models in Figs. 8a, and 8b the photoionization stops in the shallow parts of the cloud and the shielded layers, due to the higher densities, are again reached and excited only by high enough CR rates ∼10⁻¹³ − 10⁻¹² s⁻¹. This warm secondary ionized layer reaches electron temperature of T_e ∼ 8000 K, and enhances the emissivity of the low ionization lines. Similar work by Beck et al. (2023) on ionized, atomic and molecular emission lines is consistent with this result. In particular, in the structure plots (i.e., Figs. 8a, and 8b, 9a, and 9b, 9d, and 9e) it is also highlighted that when the gas density increases, CRs release their energy via collisions and thus heat the gas and boost emission lines in a less deep layer. For the lower ionization parameter, logU = −3.0, photoionization dominates the heating of the gas up to ∼10¹⁷ cm for n_H = 10³ cm⁻³ and then the CRs govern also shallower layers as shown in the teal shadowed area (Fig. 8b). On the other hand, the same ionization parameter value combined with CRs, for less dense gas of n_H = 100 cm⁻³ penetrates and excites the gas up to ∼10¹⁸ cm (Fig. 8a) and then the CRs dictate the ionization in deeper cloud areas. This behavior is consistent with the fact that the mean free path of a CR is inversely proportional to the number density of the target particles in the gas clouds.

The behavior of the electron temperature is mirrored in lower ionization lines such as [N II]λ6584 Å, [S II]λλ6716,6731 Å, and [O I]λ6300 Å, with ionization potentials of 14.53 eV, 10.36 eV, and 13.62 eV, respectively. The sensitivity of the aforementioned emission lines is also reflected on the behavior of CLOUDY models in the BPT diagrams, due to the increase of the CR ionization rate, as presented in Section 4.1, and it is evident that the lines with the lower ionization potentials are immediately affected by the increase in the CR ionization rate (Figs. 9a–9i).

Furthermore, we notice that photoionization prevails in the illuminated layers, while CRs become more significant in the deeper layers of the cloud, where the photoionization does not reach. All of the aforementioned emission lines’ structure diagrams, (Figs. 9, 10) show this deeper additional emission originating from inside of the cloud. The results of the presence of low-energy CRs are not the same on all emission lines. [O III]λ5007 Å, Hα, and Hβ are mildly affected by the higher CR rates (Figs. 9k–10f), confirming the relocation of the models towards the upper right area of the BPT diagrams and not in other direction.

When comparing the effect of different ionizing sources on the emission lines, including high or low CR ionization rates, we find that AGN and SF models exhibit both distinct differences and similarities. The ionizing radiation in AGN originates in the accretion disk around supermassive black holes, emitting a broad spectrum from X-rays to UV light, often including a significant nonthermal contribution at high energies. In contrast, the UV radiation in star-forming galaxies comes from hot and young massive stars, characterized by a thermal continuum spectrum that drops sharply above the helium double ionization edge at 54 eV (e.g. Fig. 4 in Tumlinson & Shull 2000). The different ionizing continuum in both cases is probably the cause of the smoother temperature profile vs. depth seen in AGN (Fig. 8a), as compared with the sharp profile in SF models (Fig. 8c). This is likely due to the rapid recombination in the SF nebula, while the harder continuum in AGN produces a thicker layer.

Overall, CRs have a clear impact in both AGN and SF models by increasing the gas ionization, especially in regions where the ionized radiation cannot reach (Figs. 8, 9, 10). However, the change produced in the temperature profile vs. depth is stronger in SF models, as shown by Fig. 8, and therefore a more prominent increase in the line emissivity profiles is seen in Figs. 9 and 10.

5. Discussion

5.1. Thermal stability at high CR ionization rates

In Section 4.2 we found that high CR ionization rates create a high temperature layer deep in the simulated cloud (teal-shaded area in Fig. 8), enhancing the emission from lower ionization lines (Fig. 9). To further study the effects of CRs in the thermal stability of the ionized cloud, it is convenient to adopt the definition of the ionization parameter Ξ from Krolik et al. (1981), which probes the radiation to gas pressure ratio:

$\begin{matrix} Ξ = \frac{F_{ion}}{n_{H} k_{B} T_{e} c} = \frac{L_{ion}}{n_{H} r^{2} 4 π k_{B} T_{e} c} = \frac{ξ}{4 π k_{B} T_{e} c} \end{matrix}$ $\begin{aligned} \Xi = \frac{F_{\rm ion}}{n_{\rm H} k_B T_{\text{ e}} c} = \frac{L_{\rm ion}}{n_{\rm H} r^2 4 \pi k_B T_{\text{ e}} c} = \frac{\xi }{4 \pi k_B T_{\text{ e}} c} \end{aligned}$ (3)

where F_ion and L_ion are the ionizing flux between 1 and 10³ Ryd and its corresponding luminosity, r is the distance from the cloud to the ionizing source, T_e is the electron temperature, k_B is the Boltzmann constant, c is the speed of light, and ξ is the ionization parameter definition from Tarter et al. (1969).

Thus, following the approach in Krolik et al. (1981), we analyze the thermal stability of AGN and SF models including CRs by deriving the dependence of T_e on Ξ for different values of the CR ionization rate. For this purpose, we ran additional simulations with the same characteristics as those detailed in Section 3.3, but for a wider range of ionization parameter values, for three different CR rates 10⁻¹⁴ s⁻¹, 10⁻¹³ s⁻¹, 10⁻¹² s⁻¹. The results for AGN and SF models are shown in Figs. 11a and 11b, respectively.

Fig. 11.

Thermal stability plots, depicting electron temperature as a function of the ionization parameter Ξ for different values CR ionization rate, 10⁻¹⁴ s⁻¹, 10⁻¹³ s⁻¹, 10⁻¹² s⁻¹ illustrated with purple, blue and green lines respectively. The upper x-axis shows the corresponding values for the ionization parameter U. (a) T_e in AGN models for n_H = 100 cm⁻³. (b) T_e in SF models for n_H = 100 cm⁻³.

The thermally stable solutions correspond to positive slopes in the Ξ–T_e curve, since isobaric perturbations increasing the temperature in these regions correspond also to an increase in the gas cooling efficiency. On the other hand, vertical or negative slopes correspond to thermally unstable solutions, because an increase in temperature will produce a larger contribution of heating processes, thus causing a further increase in the temperature and the amplification of any perturbation. Above Ξ ≳ −4 (log U ≳ −5), the models corresponding to the three different CR rates behave similarly, with negligible differences among the stability curves in Figs. 11a and 11b. AGN models in Fig. 11a show, as expected, two different stable gas phases (Krolik et al. 1981): ionized gas at T_e ∼ 10⁴ K for −4 ≲ Ξ ≲ 0.5, and coronal gas at ∼10⁷ K for Ξ ≳ 1.5. Similarly, Fig. 11b shows a thermally stable solution for ionized gas in star-forming regions in the −4 ≲ Ξ ≲ 2 range. Higher temperatures (≳2 × 10⁴) can be reached at extreme ionization parameter values of Ξ ≳ 3, however the lack of a hard X-ray continuum prevents the coronal stability branch as shown in the case of AGN models.

However, below Ξ ≲ −4, the thermal stability of the gas shows a strong dependence on the CR ionization rate for both AGN and SF models. Models with the highest CR ionization rate value of 10⁻¹² s⁻¹ are able to keep the ionized gas at temperatures of about 10⁴ K even at very low Ξ values. Thus, their contribution to the gas heating extends the thermal stability of the simulated cloud to regions where the ionization processes are negligible. For an ionization rate of 10⁻¹³ s⁻¹, CRs can still keep the gas in a stable phase at 10³ K, while values below ≲10⁻¹⁴ s⁻¹ are no longer able to sustain a stable solution as depicted in Figs. 11a and 11b. These results are in agreement with the theoretical analysis in Walker (2016), who demonstrates that high fluxes of low-energy CRs (≲1 GHz) are able to efficiently heat the warm ionized medium and maintain a high gas temperature against radiative cooling.

The stability diagrams demonstrate that the effect of CRs is significant at high CR ionization rates and low ionization parameter values. This is relevant for ionization-bounded clouds when the flux of UV photons does not reach the deepest layers in the gas clouds, which are shielded from the radiation and thus are heated mainly via CRs. This is consistent with the results shown in the structure plots in Figs. 9 and 10, in Section 4.2, where we observe two separate domains in the simulated cloud, one dominated by photoionization, next to the illuminated face of the cloud, and a deeper one dominated by CRs where the electron temperature remains high. The strong dependence of the temperature on the CR ionization rate observed at large cloud depths in Figs. 9 and 10, and at low Ξ in Figs. 11a and 11b, is likely associated with the dissociation of H₂ gas by the CRs.

Models developed by Ferland et al. (2009) show that, as the fraction of free electrons rises due to H₂ dissociation and ionization, the gas heating becomes increasingly more efficient, due to the larger probability of CR secondary electrons to collide with the thermal electrons. This indicates that the greatest portion of the CR energy, beyond 10⁻¹³ s⁻¹, is directed towards heating rather than H₂ dissociation leading to thermal instabilities (Ferland et al. 2009). This effect, prominent in the ionized gas, produces a steep change in the gas temperature predicted above ζ_CR ≳ 3 × 10⁻¹³ s⁻¹ (see Fig. 8 in Ferland et al. 2009).

5.2. CR heating and low ionization nebular lines

The temperature boost caused by CR heating deep in the nebula at high ζ_CR values (Fig. 8) produces an extra contribution to the low-ionization nebular lines. As shown in Section 4.2, the electron heating by CRs increases significantly the emissivity of the low-ionization collisional lines, that is, [N II]λ6584 Å, [S II]λλ6716,6731 Å, and [O I]λ6300 Å (Figs. 9a–9c and 9d–9f). In contrast, the increase in emissivity of high-ionization and recombination lines such as [O III]λ5007 Å, Hα, and Hβ is much lower (Figs. 9j–10f). Thus, the different response of the BPT lines to CR heating drives both AGN and SF models to the right side along the horizontal axes in Figs. 5, 6, and 7.

Nevertheless, the three line ratios including low-ionization lines discussed here seem to be best reproduced by different ζ_CR values. The measured [N II]/Hα ratios in Centaurus A, NGC 1068 and NGC 253 (Figs. 5c, 6c and 7c, respectively) are reproduced by models with ζ_CR = 10⁻¹² s⁻¹, while the measured [S II]/Hα, and [O I]/Hα ratios for the same regions tend to be overestimated. These are best reproduced by ζ_CR = 10⁻¹³ s⁻¹ in Figs. 5b, 6b and 7b, and possibly by ζ_CR ≲ 10⁻¹⁴ s⁻¹ for [O I]/Hα in NGC 1068 (Fig. 6a). These transitions have different ionization potentials, with [N II]λ6584 Å showing the largest value (14.5 eV), followed by [S II]λλ6716,6731 Å (10.4 eV) and [O I]λ6300 Å (neutral gas). This implies that [S II]λλ6716,6731 Å and [O I]λ6300 Å emission can be produced also in neutral gas regions located farther from the ionizing source than the N⁺ gas, as has been reported in resolved H II regions (e.g. Scowen et al. 1998; Barman et al. 2022). Additionally, the [O I]λ6300 Å line is very sensitive to shocks (Dopita 1976) and X-ray radiation (Polles et al. 2021). Thus, a possible explanation for the lower ζ_CR values required by [S II]λλ6716,6731 Å and [O I]λ6300 Å could be related to inhomogeneities in the nebular structure and the different physical conditions in the regions where these lines form, which may receive on average a lower CR flux. Photoionization modeling of several emission-line ratios with a simple geometric configuration is a challenging task, as shown by the different coverage of the observed distribution in [N II]/Hα, [S II]/Hα and [O I]/Hα, for the same grid of models (e.g. Groves et al. 2004b; Feltre et al. 2016; Zhu et al. 2023). For instance, the models that best match the [O I]/Hα ratio in the BPT diagram typically have lower ionization parameters and lower metallicities when compared to the models that best match the [N II]/Hα distribution. This suggests that an accurate description of lines formed in the neutral gas region requires a more sophisticated treatment and possibly also a more complex geometric model.

5.3. Reconciling models with the observed metallicities

One of the main points of tension between models that aim to reproduce the Seyfert/LINER locus in the BPT diagram is the supersolar metallicities (∼2–3 Z_⊙) that most models require to match the observed distributions in the emission-line ratios (e.g. Groves et al. 2004a,b; Thomas et al. 2016; Feltre et al. 2016; Zhu et al. 2023). Such high metallicities are, however, in contrast with recent measurements of the chemical abundances in the NLR of hundreds of Seyfert and LINER galaxies in the local Universe, suggesting that the nebular gas in nearby AGN typically has slightly subsolar to solar-like abundances (Dors et al. 2019, 2020, 2022; Pérez-Montero et al. 2019, 2023; Pérez-Díaz et al. 2021, 2022, 2024). These results are supported by different methods applied to emission lines in the UV, optical, and IR ranges.

Photoionization models with solar abundances tend to underestimate the [N II]/Hα ratio, and to a lesser extent, the [S II]/Hα ratio (e.g. Fig. 3 in Zhu et al. 2023). Including CR ionization allows us to reproduce the LINER/Seyfert domain in BPT diagrams using solar abundances, while previous studies proposed higher than solar metallicities, supersolar N/O ratios, or harder ionizing SEDs with steeper power-law spectral indices (Feltre et al. 2016) or hotter accretion disks (Zhu et al. 2023). In this context, the main outcome of the present study is the significant role that CR heating may play in diagnostics based on line ratios, such as the BPT diagram. As discussed in Sections 4.2 and 5.2, the heating of electrons caused by CRs deep in the nebula produces an additional contribution to the emissivity of the low-ionization collisional lines, that is, [N II]λ6584 Å, [S II]λλ6716,6731 Å and [O I]λ6300 Å (Figs. 9a–9c and 9d–9f). In contrast, the emissivity increase of high-ionization and recombination lines such as [O III]λ5007 Å, Hα, and Hβ is much lower (Figs. 9f–10f). Consequently, the different response of the BPT diagnostic lines to CR heating drives AGN models towards the right along the horizontal axes in Figs. 5, 6, and 7, shifting solar-metallicity models into the Seyfert/LINER domain. This result is summarized in Fig. 12, which shows the three BPT diagrams including the measured line ratios for all the selected regions in the three sample galaxies presented in Fig. 1. The different contours represent the 10th, 50th, and 90th percentiles of the 2D distribution of AGN and SF models. The probability distribution was obtained from the discrete model grids with 1 ≤ log n_H ≤ 3.5 and −3.5 ≤ log U ≤ −1.5, using a Gaussian kernel density estimation (Silverman 1986). The gray-hatched area represents the region occupied by photoionization models including an average Galactic CR ionization rate of 10⁻¹⁶ s⁻¹ (Indriolo et al. 2007).

Fig. 12.

BPT diagrams depicting the area covered by AGN and SF models with solar abundances for −3.5 ≤ log U ≤ −1.5, shown with purple and red contours, with 1 ≤ log n_H ≤ 3.5 and, 1 ≤ log n_H ≤ 3 respectively. The solid contour lines map regions containing 10%, 50%, and 90% of the models with ζ_CR = 10⁻¹² s⁻¹ in the [N II] BPT diagram, and 10⁻¹³ s⁻¹ in both the [S II] and the [O I] panels. The gray hatched area stands for AGN models with ζ_CR = 10⁻¹⁶ s⁻¹, of the order of the average Galactic CR background (Indriolo et al. 2007). Cyan squares represent line ratios measured for the regions selected in Centaurus A, green circles for regions in NGC 1068, magenta thin diamonds for regions in NGC 253, and the red diamond correspond to the photoionization-dominated Seyfert 2 nucleus in NGC 1320. The Kewley, Kauffmann, and Schawinski lines correspond to the solid dark gray, the solid medium gray, and dashed light gray lines, respectively.

Models with solar abundances and a CR ionization rate of ζ_CR = 10⁻¹² s⁻¹ in Fig. 12a (purple and red contours for AGN and SF models, respectively) are in excellent agreement with the observed distribution in the [O III]/Hβ and [N II]/Hα ratios. The latter are stretched along a relatively narrow stripe from the star formation domain to the Seyfert/LINER division line, which is consistent with an increasing log U in the models (see Figs. 5c, 6c, and 7c). In contrast, solar metallicity models with the average Galactic ζ_CR = 10⁻¹⁶ s⁻¹ (gray-hatched area) underestimate the [N II]/Hα ratio by an order of magnitude. Alternatively, high N/O relative abundances could also explain the enhanced [N II]/Hα ratios observed (Pérez-Montero et al. 2019; Ji et al. 2020). However, the N/O ratios derived with HII-CHI-MISTRY for the regions in our sample of galaxies suggest instead values close to the solar abundance ratio of log(N/O)_⊙ ∼ −0.86 dex (see Table 3; Asplund et al. 2009). Although CR heating can explain an increase in the emissivity of low-ionization transitions for both AGN and SF models, as shown in Fig. 12, other mechanisms such as X-ray heating or shocks may also contribute and cannot be discarded. For instance, X-ray excitation is able to produce bright [N II]λ6584 Å, [S II]λλ6716, 6731 Å, and [O I]λ6300 Å emission (e.g. Polles et al. 2021; Wolfire et al. 2022; Lagos et al. 2022). Photoionization models with shocks can also explain the high [N II]/Hα, [S II]/Hα, and [O I]/Hα ratios in AGN and SF regions (Contini & Viegas 2001a,b; Allen et al. 2008; Dors et al. 2021).

The observed [S II]/Hα and [O I]/Hα ratios in Figs. 12b and 12c, respectively, are closer to the models with the average Galactic ζ_CR value, although an additional contribution to the [S II]λλ6716,6731 Å line is still required, as well as for [O I]λ6300 Å in some of the regions in Centaurus A. A significantly lower ζ_CR = 10⁻¹³ s⁻¹ value than in the [N II] case causes a substantial stretch in the model grids along the horizontal axes, while the line ratio distributions tend to favor models with moderate densities (n_H ∼ 100 cm⁻³; see Figs. B.1b, B.2b, and B.3b). Nevertheless, SF models with ζ_CR = 10⁻¹³ s⁻¹ can explain the [S II]/Hα and [O I]/Hα ratios measured for the regions in NGC 253 that are located on the right side of the Kewley division line, beyond the limit of pure SF photoionization models.

At very low densities (0 ≤ n_H ≤ 1) with ζ_CR ≳ 10⁻¹³ s⁻¹, the models predict that [S II]/Hα and [O I]/Hα ratios fall within the LINER domain (Figs. B.1b, B.2b, and B.3b) for most of the log U values (Figs. 5b, 6b, and 7b). This suggests that heating of low-density gas by CRs could also represent an important contribution to the Diffuse Ionized Gas emission in galaxies (Vandenbroucke et al. 2018), which is characterized by enhanced [N II]/Hα, [S II]/Hα, and [O I]/Hα ratios (e.g. Zhang et al. 2017).

5.4. CR effect on metallicity and ionization parameter

Line ratios involving [N II]λ6584 Å (Marino et al. 2013; Carvalho et al. 2020; Oliveira et al. 2024), and [S II]λλ6716,6731 Å (Christensen et al. 1997; Pérez-Montero et al. 2006) have been used as metallicity tracers for gaseous nebulae. In particular, the index N2 = [N II]/Hα is widely adopted as metallicity indicator (Marino et al. 2013; Carvalho et al. 2020), and also used to compute the ionization parameter (Baron & Netzer 2019). Given the CR influence on the emission lines from low-ionization species predicted by our models, a significant impact on the metallicity and ionization parameter is expected.

To quantify the effect of CRs on metallicity and ionization parameter estimations, we applied HII-CHI-MISTRY using as input the median values of different combinations of emission-line ratios as the CR ionization rate increases from 10⁻¹⁶ s⁻¹ to 10⁻¹² s⁻¹ for both AGN and SF models. When the auroral line [O III]λ4363 Å, is included, in the HII-CHI-MISTRY routines, the temperature is constrained (Pérez-Montero 2014) and thus the metallicity estimation does not show a clear trend with increasing ζ_CR, although the derived abundances are subsolar and the associated errors are quite large (∼0.7 dex). The N/O abundance is relatively robust, as CRs affect similarly both the [O II]λ3725 Å and [N II]λ6584 Å lines, whose ratio is the main tracer used by HII-CHI-MISTRY to determine N/O. Nevertheless, the ionization parameter tends to show lower values with increasing ζ_CR for AGN models, while SF models seem to be more stable.

Considering only the BPT lines, the enhancement caused by CRs on the [N II]/Hα line ratio (see Sect. 5.2) is not balanced now by the missing [O II]λ3725 Å, thus leading to underestimated N/O values and uncertain O/H and log U estimations. The increase in the derived metallicities between ζ_CR = 10⁻¹⁶ s⁻¹ and 10⁻¹² s⁻¹ is ∼25%, and ∼65% for AGN and SF models, respectively. This is likely caused by the lower [O III]/[N II] ratios as CR ionization rate grows larger.

On the other hand, using [N II]/Hα as the only input results in larger metallicities by approximately 40%, and 45% for AGN and SF models, respectively. This is due to the enhancement of the [N II]/Hα ratio with increasing ζ_CR, which is interpreted by HII-CHI-MISTRY as a decreasing log U. For instance, in our AGN and SF models, the log U varies by 1 dex and 0.2 dex, respectively. Thus, the resulting metallicities are overestimated in SF models due to the known empirical anti-correlation between log U and O/H, which HII-CHI-MISTRY exploits. Additionally, high [N II]/Hα ratios for AGN in HII-CHI-MISTRY are only explained by models with low log U and high O/H (see Fig. 4 in Pérez-Montero et al. 2019).

Finally, we applied the metallicity calibrations based on the N2 ratio from Carvalho et al. (2020), for AGN, and Marino et al. (2013), for star-forming galaxies, to the median N2 values in our AGN and SF models with CR ionization rates in the 10⁻¹⁶ s⁻¹ to 10⁻¹² s⁻¹ range. We find that the predicted metallicities are approximately 5 times larger for 10⁻¹² s⁻¹ models when compared to 10⁻¹⁶ s⁻¹ models in the AGN case, while in the star-forming case the difference is found to be approximately 1.5 times larger. Extreme models with densities outside the 1 ≤ n_H ≤ 10³ cm⁻³ range were excluded, in all aforementioned calculations, to avoid outliers.

As we incorporated different emission line ratios into HII-CHI-MISTRY we gain some insight of their diagnostic power in the metallicity when CRs are present. This drives us to the conclusion that [N II]/Hα alone, being strongly affected by CRs, can produce metallicity estimations biased in environments where CRs are prevalent. Therefore, having the temperature estimate through [O III]λ4363 Å makes the metallicity calibrations more reliable.

5.5. Comparison with CR rate estimates

The average CR ionization rate measured in the Milky Way is of the order of ζ_CR ∼ 10⁻¹⁶ s⁻¹ (Indriolo et al. 2007, 2009; Neufeld & Wolfire 2017), while the CR ionization rate at which photoionization models show a noticeable effect on the BPT diagnostics is considerably higher (≳10⁻¹³ s⁻¹). Although higher ζ_CR values are anticipated for active nuclei and strong SF environments, produced by nonthermal processes in AGN jets and SNRs (Meijerink et al. 2011; Guo & Mathews 2011; Phan et al. 2024), a difference of approximately three orders of magnitude raises the question of how realistic this scenario is for the regions under study in this work. Thus, in this Section we compare the values adopted in our models with observational measurements and independent estimates of CR ionization rates in AGN and SF galaxies found in the literature.

CR ionization rates cannot be directly measured, and therefore our comparison is based on available indirect estimates for nearby galaxies. The CR ionization and heating of H and H₂ initiates a series of chemical reactions in molecular gas clouds that lead to the formation of molecular ions such as OH⁺, H₂O⁺, and H₃O⁺ (Herbst & Klemperer 1973; González-Alfonso et al. 2013). Therefore, the observed transitions from these molecules can be compared with chemical models that predict their relative abundances as a function of ζ_CR. Following this approach, in Mrk 231, closest hosting quasar galaxy, CR ionization rates ∼5 × 10⁻¹³ in torus, and ∼10⁻¹² s⁻¹ in outflows were found (González-Alfonso et al. 2018). Likewise, Holdship et al. (2022) derives ζ_CR ∼ 10⁻¹³ s⁻¹ within the innermost few hundred parsecs in NGC 253. Slightly higher values (10⁻¹³–10⁻¹² s⁻¹) were obtained by Behrens et al. (2022) by modeling the HCN/HNC ratio for the same galaxy. These values are in agreement with the results found by González-Alfonso et al. (2013) in the starburst/AGN composite nuclei of NGC 4418 and Arp 220. Finally, values significantly higher than the Galactic average (10⁻¹⁴–10⁻¹³ s⁻¹) have also been inferred in diffuse molecular clouds near the Galactic center (Le Petit et al. 2016), where a high CR flux is expected.

An estimate of the CR ionization rate can also be derived from the nonthermal radio continuum emission (e.g. Yusef-Zadeh et al. 2013; Gabici 2022). In AGN, the synchrotron continuum measured along the jet can be used as a proxy for the number density of accelerated electrons. The energy distribution of relativistic electrons can be inferred from the slope of the synchrotron emission continuum (Rybicki & Lightman 1979),

$\begin{matrix} N_{e} d E \propto E^{- p} d E, \end{matrix}$ $\begin{aligned} N_{\text{ e}} \mathrm{d}E \propto E^{-p}\mathrm{d}E, \end{aligned}$ (4)

where N_e is the number density of electron with energies in the range E up to E + dE, with E being the energy of electron (Rybicki & Lightman 1979). The magnetic field strength is required to translate the radiated luminosity into a number density of particles (Padovani et al. 2018), which is done by assuming equipartition between the energy density of relativistic electrons and the magnetic field energy density of the synchrotron emitting area. Therefore, with this method, we infer the electron component of the CRs, by using the corresponding ionization rate from primary and secondary electrons – as displayed below – obtained for Rq. (40) in Gabici (2022), and we derive an estimate of the CR ionization rate of atomic hydrogen via the assumption that $1.5 ζ_{CR}^{H_{2}} ≃ 2.3 ζ_{CR}^{H}$ $1.5\,\zeta_{\text{ CR}}^{\text{ H}_2} \simeq 2.3\,\zeta_{\text{ CR}}^{\text{ H}}$ (Glassgold & Langer 1974).

$\begin{matrix} ζ_{CR}^{H} ≃ 1.3 π \int_{E_{min}}^{E_{max}} d E j_{e} (E) σ_{H_{2} - e}^{ion} (E) [1 + ϕ_{e, H_{2}} (E)], \end{matrix}$ $\begin{aligned} \zeta _{\text{ CR}}^{\text{ H}} \simeq 1.3\pi \int _{E_{\min }}^{E_{\max }} \mathrm{d}E \, j_{\text{ e}}(E) \, \sigma ^{\text{ ion}} _{\text{ H}_2-\text{ e}}(E) \left[1+\phi _{\text{ e}, \,\text{ H}_2}(E) \right], \end{aligned}$ (5)

where j_e(E) = cN_e/(4π), with c being the speed of light, is the differential energy spectrum of CR electrons, representing the number of electrons with energy E per unit energy, per unit area, per unit time, and per unit solid angle, ϕ_e, H₂(E) represents the secondary ionization by electrons, and we adopt the cross-section of H₂ − e⁻ interactions by Padovani et al. (2018), $σ_{H_{2} - e}^{ion} (E)$ $\sigma^{\text{ ion}}_{\text{ H}_2-{\text{ e}}}(E)$ . Furthermore, we assume E_min = 0.1 MeV as in Yusef-Zadeh et al. (2013) and E_max = 10 GeV as the upper energy limit which only weakly affects the outcome.

Specifically, in Centaurus A we calculate the equipartition magnetic field of the emitting radio source A1A (Goodger et al. 2010) and for an index p ∼ 2.7 coming from an indicative index α ∼ 0.85 based on the α indices of table 7 in (Goodger et al. 2010), we find it to be 300 μG. The SED and the emitting source radius used in this calculation are provided in Table 4 of Goodger et al. (2010), while the position of each source is depicted in Fig. 2 of the same work. We then acquire a CR ionization rate ∼2 × 10⁻¹⁰ s⁻¹ in this radio knot. We then repeat the calculation, assuming equipartition, and the same index p, in the northern lobe of Centaurus A in the area that covers A1A up to A2A radio knots (1 kpc × 100 pc) of Goodger et al. (2010), which corresponds to an effective radius of ∼300 pc. The flux used is the integrated flux of radio knots A1A up to A2A (see Table 4 in Goodger et al. 2010). Thus, we derive an equipartition magnetic field of ∼8 μG and a CR ionization rate ∼4.3 × 10⁻¹³ s⁻¹.

In an identical fashion, for the emitting source marked as C in in NGC 1068 (Mutie et al. 2024), we calculate an equipartition magnetic field of ∼700 μG. The SED, and emitting source radius used in this calculation are extracted from Table 2 in Mutie et al. (2024), while the position of the source can be seen in Fig. 2 of that work. Assuming an electron energy index p ∼ 2.8, derived from a spectral index α ∼ 0.9 (see Table 2 in Mutie et al. 2024), we determine a CR ionization rate of about ∼2 × 10⁻⁹ s⁻¹. For an effective radius of 160 pc, representative of the surface of the jet between the components NE up to S3 (500 pc × 50 pc), and using the integrated radio fluxes provided by Table 2 in Mutie et al. (2024), we estimate an equipartition magnetic field of ∼36 μG, and a CR ionization rate of ∼1.4 × 10⁻¹¹ s⁻¹.

The CR ionization rate values inferred for Centaurus A and NGC 1068 are indicative of the typical values that can be found along the trail of AGN jets. We expect these values to decrease as we move away from the jet trail, but the numbers obtained from the synchrotron continuum suggest that the order of magnitude of the CR rate in these galaxies is in agreement with the range of values explored in this work. It is also crucial to bear in mind that heavier particles such as protons are not traced by the synchrotron continuum, but can also contribute to gas heating and ionization, suggesting a higher value for the total CR rate (Padovani et al. 2018). Finally, the CR rates estimated from this method could be subject to future refinement and more detailed analysis.

Overall, both indirect measurements in SF galaxies based on the abundances of molecular transitions, and theoretical estimates based on the synchrotron continuum emission in AGN jets suggest that the CR ionization rates in starburst and active nuclei is expected to be about three to four orders of magnitude higher than the Galactic average, in agreement with the ζ_CR = 10⁻¹⁴–10⁻¹² s⁻¹ range probed by our models.

5.6. A new SFζ line for starburst galaxies with CRs

Our analysis shows that certain areas in the LINER or Seyfert domain of classical BPT diagrams, which are beyond the capabilities of pure photoionization star-forming models (Kewley et al. 2001), can instead be explained when typical CR ionization rate values, as those observed in starburst or expected in radio galaxies, are taken into account. For instance, several regions in NGC 253, Centaurus A, and NGC 1068 exhibit line ratios consistent with SF plus ζ_CR models (red contour in Fig. 12). Thus, we propose new maximum starburst lines for BPT diagrams to account for the potential contribution of CRs, which may alter the low-ionization line to Hα ratios beyond the boundaries defined by Kewley et al. (2001).

Fig. 13.

BPT diagrams depicting [N II]/Hα, [S II]/Hα, and [O I]/Hα ratios. Our proposed SFζ maximum starburst line is illustrated in red color, compared with the full observational dataset. Observations from Centaurus A, NGC 1068, NGC 253, and NGC 1320 are marked by cyan squares, green circles, magenta thin diamonds, and a red diamond, respectively. The Kewley and Schawinski lines are indicated with solid and dashed black, respectively, while the SFζ line is depicted by the red-solid 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).

To define the new maximum starburst lines, we used our SF models with densities in the 1 ≤ log n_H ≤ 3 range, which are typical values observed in star-forming regions (e.g. Pérez-Montero & Díaz 2003). Subsequently, we consider only the area including 90% of the 2D probability density kernel associated with the model distribution (red contour in Fig. 12), to exclude extreme values from the definition of the new boundaries. The new maximum line for starburst plus CR contribution, hereafter referred as SFζ is obtained by fitting the 90% contour boundary in each case, taking also into account the Kewley et al. (2001) model values at the lowest [N II]/Hα, [S II]/Hα, and [O I]/Hα ratios. For this purpose, we define the following X and Y variables:

$\begin{matrix} X = {\frac{[N II]}{H α}, \frac{[S II]}{H α}, \frac{[O I]}{H α}}, Y = \frac{[O III]}{H β} \end{matrix}$ $\begin{aligned} X = \left\{ \frac{[{\text{ N}}{\small { {\text{ II}}}}]}{\mathrm{H}\alpha }, \frac{[{\text{ S}}{\small { {\text{ II}}}}]}{\mathrm{H}\alpha },\frac{[{\text{ O}}{\small {{\text{ I}}}}]}{\mathrm{H}\alpha } \right\} , \quad Y = \frac{[{\text{ O}}{\small { {\text{ III}}}}]}{\mathrm{H}\beta } \end{aligned}$ (6)

while the relationship between them is characterized by the equation:

$\begin{matrix} X = (a Y^{3} + b Y^{2} + c Y + d) e^{fY} . \end{matrix}$ $\begin{aligned} X = \left(a Y^3 + b Y^2 + c Y + d\right) \, \mathrm{e} ^{f Y}. \end{aligned}$ (7)

The best-fit parameters obtained for the different BPT diagrams are listed in Table 4. The proposed SFζ boundaries, shown by the solid red lines in Fig. 13, can be used to distinguish regions that are unambiguously dominated by AGN photoionization from those that could be explained by starburst plus CR ionization, up to values of about ζ_CR ≲ 10⁻¹² s⁻¹. The latter scenario is perhaps less likely when the line fluxes are measured within large apertures in relatively quiescent galaxies, as suggested by the lack of SDSS galaxies near the upper side of the SFζ line in Fig. 13. However, the new separation line is more appropriate for line fluxes derived from resolved regions within galaxies, since the emission lines from low-excitation gas in these cases can be significantly affected when the regions are exposed to nearby sources with high CR rates. This is the case for some regions in the center of NGC 253, which are located above the Kewley line in Fig. 13, suggesting an additional source of excitation apart from star formation. According to the new SFζ line, these regions are consistent with starburst plus CR ionization.

Table 4.

Parameters defining the SFζ boundaries in BPT diagrams, i.e., the maximum contribution predicted for models including ionization from star-forming regions and CRs.

5.7. Further remarks and considerations

Our analysis shows that the addition of CRs to photoionization models can significantly increase the nebular emission from low-ionization species, due to the temperature increase in deep gas layers located beyond the Strömgren radius (see Section 4.2). The predicted temperatures drop rapidly as ζ_CR decreases, reaching ≲10³ K for ≲10⁻¹³ s⁻¹ (Fig. 11). At such low temperatures, the emissivities of optical nebular lines decrease by orders of magnitude, because the collisions do not have enough energy to populate the energy levels involved in these transitions. However, the infrared nebular lines are produced by transitions from levels located much closer to the ground state, which become populated at low temperatures (≲1000 K), and therefore their emissivities have a negligible dependency with the temperature (e.g. Fig. 1 in Fernández-Ontiveros et al. 2021). Thus, IR lines from low-ionization species should be more sensitive to the effect of CR heating. We will further study this aspect in a future work.

Beyond their role on nebular emission lines, CRs also play a significant role in the dynamical stability of gas clouds and the regulation of star formation (Naab & Ostriker 2017; Girichidis et al. 2024). Considerable progress has been achieved in simulating the influence of CRs as a driving mechanism for both warm and cool outflows (Breitschwerdt et al. 1991; Armillotta et al. 2024; Rodríguez Montero et al. 2024). Simulations of CR driven winds are nowadays 3D magnetohydrodynamical (MHD) and quite elaborate, taking into account magnetic fields, radiative cooling, self-gravity, star formation, the dynamical role of CRs injected by SNe, and CR transport processes–anisotropic diffusion and CR streaming along the magnetic fields (Uhlig et al. 2012; Rosdahl et al. 2015; Girichidis et al. 2018; Peschken et al. 2023; Armillotta et al. 2024). Simulations like this support that the presence of CRs quenches star formation due to them being able to move cooler gas away from the galactic disk (Jubelgas et al. 2008; Ruszkowski et al. 2017; Peschken et al. 2023). Our work supports the idea that CRs might be able to sustain ionization in outflowing material farther away, potentially extending the observable impact and physical reach of outflows. We found that CRs deposit energy and sustain higher temperatures in deeper parts of clouds, and by doing so, they possibly assist in their dispersion and therefore in the quenching of star formation processes.

6. Summary

In this work, we use photoionization simulations with CLOUDY to investigate the impact of CRs on the ionized gas phase of AGN and starburst galaxies. The models sample a wide range in density (1 to 10⁴ cm⁻³), ionization parameter (−3.5 ≤ log U ≤ −1.5), and CR ionization rate (10⁻¹⁵ to 10⁻¹² s⁻¹), and are compared with VLT/MUSE spectroscopic observations of two jet-dominated AGN (Centaurus A and NGC 1068) and a prototypical starburst galaxy (NGC 253).

We find that CR rates of the order of 10⁻¹³ s⁻¹ are able to significantly affect the thermal structure of the ionized gas by creating a secondary low-ionization layer of warm ionized gas deep within the cloud, and far beyond the photoionization-dominated region surrounding the central source. This secondary ionized gas layer produces bright emission from low-ionization transitions such as [N II]λ6584 Å, [S II]λλ6716,6731 Å, and [O I]λ6300 Å, whereas [O III]λ5007 Å, Hα, and Hβ are marginally enhanced. These findings are supported by analyzing the line emissivities of these transitions as a function of cloud depth. Thus, high CR ionization rates can significantly affect classical line-ratio-based diagnostics such as BPT diagrams. CR ionization rates of the order of 10⁻¹³ s⁻¹ have been measured in both jet-dominated AGN and strong starburst nuclei, which is in agreement with the values adopted in our simulations.

In contrast with pure photoionization models for NLR, which require between two and three times the solar metallicity to reproduce the Seyfert loci in BPT diagrams, simulations for AGN photoionization with high CR ionization rates are able to reproduce the Seyfert distribution in BPT diagrams with solar-like metallicities due to the CR boost of low-ionization transitions. Metallicity estimates derived for the regions extracted from the MUSE datacubes in our sample of galaxies suggest solar to subsolar values (0.4–1.2 Z_⊙; see Section 3.2 and Table 3). Additionally, our photoionization simulations of star-forming regions with high CR ionization rates can also explain the presence of non-AGN sources in the LINER domain. Consequently, we propose a new maximum starburst SFζ boundary (Section 5.6) that extends the previous limit defined by Kewley et al. (2001) to distinguish regions dominated by AGN photoionization from those that can be explained by star formation plus CR ionization.

The results of this study suggest that the gas ionization and heating of the ISM can be largely affected in environments with high CR rates. In particular, CRs can play a major role in the excitation of low-ionization nebular gas in environments close to jet-dominated AGN and starburst galaxies, where nonthermal processes are important. SNRs originated by strong star-forming activity and particle acceleration and shocks along AGN jets are the main sources of CRs in galaxies. CRs can also significantly impact the metallicity and physical parameters derived from line ratios involving low-ionization and neutral species. Therefore, future studies should consider the potential contribution of CR heating when high rates of these particles are expected.

¹

HII-CHI-MISTRY is publicly available at: https://home.iaa.csic.es/~epm/HII-CHI-mistry.html

Acknowledgments

The authors are grateful to the anonymous referee for their insightful and constructive comments, for suggesting the inclusion of the SFζ line, and for overall improving this work. EK extends heartfelt gratitude to Dr. Steven Tingay and Dr. Emil Lenc for generously providing VLA data of Centaurus A, and Mr. Isaac Mutie for providing the combined e-Merlin plus VLA data of NGC 1068. EK is thankful to Prof. Despina Hatzidimitriou and Prof. Apostolos Mastichiadis for insightful discussions and their guidance and to Dr. Christophe Morisset for all his help regarding CLOUDY and pyCLOUDY. 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. Based on data obtained from the ESO Science Archive Facility with DOI: https://doi.org/10.18727/archive/41. Funding for the Sloan Digital Sky Survey V has been provided by the Alfred P. Sloan Foundation, the Heising-Simons Foundation, the National Science Foundation, and the Participating Institutions. DSS acknowledges support and resources from the Center for High-Performance Computing at the University of Utah. SDSS telescopes are located at Apache Point Observatory, funded by the Astrophysical Research Consortium and operated by New Mexico State University, and at Las Campanas Observatory, operated by the Carnegie Institution for Science. The SDSS website is www.sdss.org.

References

Abazajian, K. N., Adelman-McCarthy, J. K., Agüeros, M. A., et al. 2009, ApJS, 182, 543 [Google Scholar]
Allen, M. G., Groves, B. A., Dopita, M. A., Sutherland, R. S., & Kewley, L. J. 2008, ApJS, 178, 20 [Google Scholar]
Armillotta, L., Ostriker, E. C., Kim, C.-G., & Jiang, Y.-F. 2024, ApJ, 964, 99 [NASA ADS] [CrossRef] [Google Scholar]
Asplund, M., Grevesse, N., Sauval, A. J., & Scott, P. 2009, ARA&A, 47, 481 [NASA ADS] [CrossRef] [Google Scholar]
Bacon, R., Piqueras, L., Conseil, S., Richard, J., & Shepherd, M. 2016, Astrophysics Source Code Library [record ascl:1611.003] [Google Scholar]
Bacon, R., Brinchmann, J., Conseil, S., et al. 2023, A&A, 670, A4 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Baldwin, J. A., Phillips, M. M., & Terlevich, R. 1981, PASP, 93, 5 [Google Scholar]
Baloković, M., Comastri, A., Harrison, F. A., et al. 2014, ApJ, 794, 111 [CrossRef] [Google Scholar]
Barbosa, F. K. B., Storchi-Bergmann, T., McGregor, P., Vale, T. B., & Rogemar Riffel, A. 2014, MNRAS, 445, 2353 [Google Scholar]
Barman, S., Neelamkodan, N., Madden, S. C., et al. 2022, ApJ, 930, 100 [NASA ADS] [CrossRef] [Google Scholar]
Baron, D., & Netzer, H. 2019, MNRAS, 486, 4290 [Google Scholar]
Beck, A., Lebouteiller, V., Madden, S. C., et al. 2023, A&A, 680, A55 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Behrens, E., Mangum, J. G., Holdship, J., et al. 2022, ApJ, 939, 119 [NASA ADS] [CrossRef] [Google Scholar]
Bicknell, G. V., Sutherland, R. S., van Breugel, W. J. M., et al. 2000, ApJ, 540, 678 [NASA ADS] [CrossRef] [Google Scholar]
Bland-Hawthorn, J., Lumsden, S., Voit, G., Cecil, G., & Weisheit, J. 1997, Astrophys. Space Sci., 248, 177 [NASA ADS] [CrossRef] [Google Scholar]
Blasi, P. 2013, A&ARv, 21, 70 [Google Scholar]
Breitschwerdt, D., McKenzie, J. F., & Voelk, H. J. 1991, A&A, 245, 79 [NASA ADS] [Google Scholar]
Brinchmann, J., Charlot, S., White, S. D. M., et al. 2004, MNRAS, 351, 1151 [Google Scholar]
Brinchmann, J., Charlot, S., Kauffmann, G., et al. 2013, MNRAS, 432, 2112 [CrossRef] [Google Scholar]
Bruzual, G., & Charlot, S. 2003, MNRAS, 344, 1000 [NASA ADS] [CrossRef] [Google Scholar]
Buck, T., Pfrommer, C., Pakmor, R., Grand, R. J. J., & Springel, V. 2020, MNRAS, 497, 1712 [NASA ADS] [CrossRef] [Google Scholar]
Calzetti, D. 2001, PASP, 113, 1449 [Google Scholar]
Carilli, C. L., Holdaway, M. A., Ho, P. T. P., & de Pree, C. G. 1992, ApJ, 399, L59 [Google Scholar]
Carvalho, S. P., Dors, O. L., Cardaci, M. V., et al. 2020, MNRAS, 492, 5675 [NASA ADS] [CrossRef] [Google Scholar]
Cecil, G., Bland, J., & Tully, R. B. 1990, ApJ, 355, 70 [Google Scholar]
Chatzikos, M., Bianchi, S., Camilloni, F., et al. 2023, Rev. Mex. Astron. Astrofis., 59, 327 [Google Scholar]
Christensen, T., Petersen, L., & Gammelgaard, P. 1997, A&A, 322, 41 [NASA ADS] [Google Scholar]
Contini, M., & Viegas, S. M. 2001a, ApJS, 137, 75 [NASA ADS] [CrossRef] [Google Scholar]
Contini, M., & Viegas, S. M. 2001b, ApJS, 132, 211 [NASA ADS] [CrossRef] [Google Scholar]
Dalgarno, A. 2006, Proc. Nat. Academy Sci., 103, 12269 [NASA ADS] [CrossRef] [Google Scholar]
De Robertis, M. M., & Osterbrock, D. E. 1986, ApJ, 301, 98 [NASA ADS] [CrossRef] [Google Scholar]
Dopita, M. A. 1976, ApJ, 209, 395 [NASA ADS] [CrossRef] [Google Scholar]
Dopita, M. A., & Evans, I. N. 1986, ApJ, 307, 431 [NASA ADS] [CrossRef] [Google Scholar]
Dopita, M. A., & Sutherland, R. S. 1995, ApJ, 455, 468 [Google Scholar]
Dors, O. L. 2021, MNRAS, 507, 466 [NASA ADS] [CrossRef] [Google Scholar]
Dors, O. L., Monteiro, A. F., Cardaci, M. V., Hägele, G. F., & Krabbe, A. C. 2019, MNRAS, 486, 5853 [NASA ADS] [CrossRef] [Google Scholar]
Dors, O. L., Maiolino, R., Cardaci, M. V., et al. 2020, MNRAS, 496, 3209 [NASA ADS] [CrossRef] [Google Scholar]
Dors, O. L., Contini, M., Riffel, R. A., et al. 2021, MNRAS, 501, 1370 [Google Scholar]
Dors, O. L., Valerdi, M., Freitas-Lemes, P., et al. 2022, MNRAS, 514, 5506 [NASA ADS] [CrossRef] [Google Scholar]
Dufour, R. J., & van den Bergh, S. 1978, ApJ, 226, L73 [NASA ADS] [CrossRef] [Google Scholar]
Earl, N., Tollerud, E., O’Steen, R., et al. 2023, https://doi.org/10.5281/zenodo.10016569 [Google Scholar]
Eilek, J. A. 2014, New J. Phys., 16, 045001 [NASA ADS] [CrossRef] [Google Scholar]
Evans, I. N., & Dopita, M. A. 1985, ApJS, 58, 125 [NASA ADS] [CrossRef] [Google Scholar]
Feltre, A., Charlot, S., & Gutkin, J. 2016, MNRAS, 456, 3354 [Google Scholar]
Ferguson, J. W., Korista, K. T., Baldwin, J. A., & Ferland, G. J. 1997, ApJ, 487, 122 [NASA ADS] [CrossRef] [Google Scholar]
Ferland, G. J., & Mushotzky, R. F. 1984, ApJ, 286, 42 [NASA ADS] [CrossRef] [Google Scholar]
Ferland, G. J., Fabian, A. C., Hatch, N. A., et al. 2009, MNRAS, 392, 1475 [NASA ADS] [CrossRef] [Google Scholar]
Ferland, G. J., Porter, R. L., van Hoof, P. A. M., et al. 2013, Rev. Mex. Astron. Astrofis., 49, 137 [Google Scholar]
Ferland, G. J., Chatzikos, M., Guzmán, F., et al. 2017, Rev. Mex. Astron. Astrofis., 53, 385 [NASA ADS] [Google Scholar]
Fernández-Ontiveros, J. A., Prieto, M. A., & Acosta-Pulido, J. A. 2009, MNRAS, 392, L16 [CrossRef] [Google Scholar]
Fernández-Ontiveros, J. A., Pérez-Montero, E., Vílchez, J. M., Amorín, R., & Spinoglio, L. 2021, A&A, 652, A23 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Ferruit, P., Wilson, A. S., & Mulchaey, J. 2000, ApJS, 128, 139 [Google Scholar]
Fragile, P. C., Murray, S. D., Anninos, P., & van Breugel, W. 2004, ApJ, 604, 74 [NASA ADS] [CrossRef] [Google Scholar]
Gabici, S. 2022, A&ARv, 30, 4 [NASA ADS] [CrossRef] [Google Scholar]
Galliano, E., Alloin, D., Pantin, E., Lagage, P. O., & Marco, O. 2005, A&A, 438, 803 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Gallimore, J. F., Baum, S. A., O’Dea, C. P., & Pedlar, A. 1996, ApJ, 458, 136 [Google Scholar]
Gallimore, J. F., Baum, S. A., & O’Dea, C. P. 2004, ApJ, 613, 794 [Google Scholar]
Gallimore, J. F., Yzaguirre, A., Jakoboski, J., et al. 2010, ApJS, 187, 172 [NASA ADS] [CrossRef] [Google Scholar]
Gallimore, J. F., Elitzur, M., Maiolino, R., et al. 2016, ApJ, 829, L7 [Google Scholar]
García-Burillo, S., Combes, F., Ramos Almeida, C., et al. 2016, ApJ, 823, L12 [Google Scholar]
Girichidis, P., Naab, T., Hanasz, M., & Walch, S. 2018, MNRAS, 479, 3042 [NASA ADS] [CrossRef] [Google Scholar]
Girichidis, P., Werhahn, M., Pfrommer, C., Pakmor, R., & Springel, V. 2024, MNRAS, 527, 10897 [Google Scholar]
Glassgold, A. E., & Langer, W. D. 1974, ApJ, 193, 73 [NASA ADS] [CrossRef] [Google Scholar]
González-Alfonso, E., Fischer, J., Bruderer, S., et al. 2013, A&A, 550, A25 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
González-Alfonso, E., Fischer, J., Bruderer, S., et al. 2018, ApJ, 857, 66 [Google Scholar]
Goodger, J. L., Hardcastle, M. J., Croston, J. H., et al. 2010, ApJ, 708, 675 [NASA ADS] [CrossRef] [Google Scholar]
Gredel, R., Lepp, S., Dalgarno, A., & Herbst, E. 1989, ApJ, 347, 289 [Google Scholar]
Groves, B. A., Dopita, M. A., & Sutherland, R. S. 2004a, ApJS, 153, 9 [NASA ADS] [CrossRef] [Google Scholar]
Groves, B. A., Dopita, M. A., & Sutherland, R. S. 2004b, ApJS, 153, 75 [NASA ADS] [CrossRef] [Google Scholar]
Guo, F., & Mathews, W. G. 2011, ApJ, 728, 121 [NASA ADS] [CrossRef] [Google Scholar]
Hardcastle, M. J., Worrall, D. M., Kraft, R. P., et al. 2003, ApJ, 593, 169 [NASA ADS] [CrossRef] [Google Scholar]
Harris, G. L. H., Rejkuba, M., & Harris, W. E. 2010, PASA, 27, 457 [NASA ADS] [CrossRef] [Google Scholar]
Heesen, V., Beck, R., Krause, M., & Dettmar, R.-J. 2008, A&A, 494, 563 [Google Scholar]
Heesen, V., Beck, R., Krause, M., & Dettmar, R. J. 2011, A&A, 535, A79 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Herbst, E., & Klemperer, W. 1973, ApJ, 185, 505 [Google Scholar]
Holdship, J., Mangum, J. G., Viti, S., et al. 2022, ApJ, 931, 89 [NASA ADS] [CrossRef] [Google Scholar]
Humire, P. K., Henkel, C., Hernández-Gómez, A., et al. 2022, A&A, 663, A33 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Indriolo, N., Geballe, T. R., Oka, T., & McCall, B. J. 2007, ApJ, 671, 1736 [NASA ADS] [CrossRef] [Google Scholar]
Indriolo, N., Fields, B. D., & McCall, B. J. 2009, ApJ, 694, 257 [NASA ADS] [CrossRef] [Google Scholar]
Israel, F. 1998, A&ARv., 8, 237 [NASA ADS] [CrossRef] [Google Scholar]
Israel, F. P., Güsten, R., Meijerink, R., Requena-Torres, M. A., & Stutzki, J. 2017, A&A, 599, A53 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Janssen, M., Falcke, H., Kadler, M., et al. 2021, Nat. Astron., 5, 1017 [NASA ADS] [CrossRef] [Google Scholar]
Ji, X., Yan, R., Riffel, R., Drory, N., & Zhang, K. 2020, MNRAS, 496, 1262 [NASA ADS] [CrossRef] [Google Scholar]
Joseph, P., Sreekumar, P., Stalin, C. S., et al. 2022, MNRAS, 516, 2300 [NASA ADS] [CrossRef] [Google Scholar]
Jubelgas, M., Springel, V., Enßlin, T., & Pfrommer, C. 2008, A&A, 481, 33 [CrossRef] [EDP Sciences] [Google Scholar]
Kantzas, D., Markoff, S., Cooper, A. J., et al. 2023, MNRAS, 524, 1326 [NASA ADS] [CrossRef] [Google Scholar]
Kapinska, A. D., Staveley-Smith, L., Crocker, R., et al. 2017, ApJ, 838, 68 [CrossRef] [Google Scholar]
Kauffmann, G., Heckman, T. M., Tremonti, C., et al. 2003, MNRAS, 346, 1055 [Google Scholar]
Kewley, L. J., Dopita, M. A., Sutherland, R. S., Heisler, C. A., & Trevena, J. 2001, ApJ, 556, 121 [Google Scholar]
Kewley, L. J., Groves, B., Kauffmann, G., & Heckman, T. 2006, MNRAS, 372, 961 [Google Scholar]
Kim, J.-G., Kim, W.-T., & Ostriker, E. C. 2018, ApJ, 859, 68 [NASA ADS] [CrossRef] [Google Scholar]
Kraemer, S. B., Schmitt, H. R., Crenshaw, D. M., et al. 2011, ApJ, 727, 130 [NASA ADS] [CrossRef] [Google Scholar]
Kraft, R. P., Hardcastle, M. J., Sivakoff, G. R., et al. 2008, ApJ, 677, L97 [NASA ADS] [CrossRef] [Google Scholar]
Krolik, J. H., McKee, C. F., & Tarter, C. B. 1981, ApJ, 249, 422 [NASA ADS] [CrossRef] [Google Scholar]
Lagos, P., Loubser, S. I., Scott, T. C., et al. 2022, MNRAS, 516, 5487 [NASA ADS] [CrossRef] [Google Scholar]
Le Petit, F., Ruaud, M., Bron, E., et al. 2016, A&A, 585, A105 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Leitherer, C., Schaerer, D., Goldader, J. D., et al. 1999, ApJS, 123, 3 [Google Scholar]
Lenc, E., & Tingay, S. J. 2006, AJ, 132, 1333 [NASA ADS] [CrossRef] [Google Scholar]
Lin, Y.-H., Yang, H. Y. K., & Owen, E. R. 2023, MNRAS, 520, 963 [CrossRef] [Google Scholar]
Lucero, D. M., Carignan, C., Elson, E. C., et al. 2015, MNRAS, 450, 3935 [NASA ADS] [CrossRef] [Google Scholar]
Maiolino, R., & Mannucci, F. 2019, A&ARv, 27, 3 [Google Scholar]
Marino, R. A., Rosales-Ortega, F. F., Sánchez, S. F., et al. 2013, A&A, 559, A114 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Matsuoka, K., Nagao, T., Maiolino, R., Marconi, A., & Taniguchi, Y. 2009, A&A, 503, 721 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
McKee, C. F. 1989, ApJ, 345, 782 [NASA ADS] [CrossRef] [Google Scholar]
McKinley, B., Tingay, S. J., Carretti, E., et al. 2018, MNRAS, 474, 4056 [NASA ADS] [CrossRef] [Google Scholar]
Meijerink, R., Spaans, M., Loenen, A. F., & van der Werf, P. P. 2011, A&A, 525, A119 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Mills, E. A. C., Gorski, M., Emig, K. L., et al. 2021, ApJ, 919, 105 [NASA ADS] [CrossRef] [Google Scholar]
Mingozzi, M., Cresci, G., Venturi, G., et al. 2019, A&A, 622, A146 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Morganti, R. 2010, PASA, 27, 463 [NASA ADS] [CrossRef] [Google Scholar]
Morisset, C. 2013, Astrophysics Source Code Library [record ascl:1304.020] [Google Scholar]
Mutie, I. M., Williams-Baldwin, D., Beswick, R. J., et al. 2024, MNRAS, 527, 11756 [Google Scholar]
Muxlow, T. W. B., Pedlar, A., Holloway, A. J., Gallimore, J. F., & Antonucci, R. R. J. 1996, MNRAS, 278, 854 [NASA ADS] [CrossRef] [Google Scholar]
Naab, T., & Ostriker, J. P. 2017, ARA&A, 55, 59 [Google Scholar]
Nagao, T., Maiolino, R., & Marconi, A. 2006, A&A, 447, 863 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Neufeld, D. A., & Wolfire, M. G. 2017, ApJ, 845, 163 [Google Scholar]
Newville, M., Stensitzki, T., Allen, D. B., & Ingargiola, A. 2015, https://doi.org/10.5281/zenodo.11813 [Google Scholar]
Oliveira, C. B., Krabbe, A. C., Dors, O. L., et al. 2024, MNRAS, 531, 199 [NASA ADS] [CrossRef] [Google Scholar]
Osterbrock, D. E., & Ferland, G. J. 2006, Astrophysics of gaseous nebulae and active galactic nuclei (University Science Books) [Google Scholar]
Pacholczyk, A. G. 1970, Radio astrophysics. Nonthermal processes in galactic and extragalactic sources (W.H.Freeman& Co Ltd.) [Google Scholar]
Padovani, P., Alexander, D. M., Assef, R. J., et al. 2017, A&ARv, 25, 2 [Google Scholar]
Padovani, M., Ivlev, A. V., Galli, D., & Caselli, P. 2018, A&A, 614, A111 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Pérez-Díaz, B., Masegosa, J., Márquez, I., & Pérez-Montero, E. 2021, MNRAS, 505, 4289 [CrossRef] [Google Scholar]
Pérez-Díaz, B., Pérez-Montero, E., Fernández-Ontiveros, J. A., & Vílchez, J. M. 2022, A&A, 666, A115 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Pérez-Díaz, B., Pérez-Montero, E., Fernández-Ontiveros, J. A., et al. 2024, A&A, 685, A168 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Pérez-Montero, E. 2014, MNRAS, 441, 2663 [CrossRef] [Google Scholar]
Pérez-Montero, E., & Díaz, A. I. 2003, MNRAS, 346, 105 [CrossRef] [Google Scholar]
Pérez-Montero, E., & Vílchez, J. M. 2009, MNRAS, 400, 1721 [CrossRef] [Google Scholar]
Pérez-Montero, E., Díaz, A. I., Vílchez, J. M., & Kehrig, C. 2006, A&A, 449, 193 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Pérez-Montero, E., Dors, O. L., Vílchez, J. M., et al. 2019, MNRAS, 489, 2652 [CrossRef] [Google Scholar]
Pérez-Montero, E., Amorín, R., Pérez-Díaz, B., Vílchez, J. M., & García-Benito, R. 2023, MNRAS, 521, 1556 [CrossRef] [Google Scholar]
Pérez-Montero, E., Fernández-Ontiveros, J. A., Pérez-Díaz, B., et al. 2024, A&A, 684, A40 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Peschken, N., Hanasz, M., Naab, T., Wóltański, D., & Gawryszczak, A. 2023, MNRAS, 522, 5529 [NASA ADS] [CrossRef] [Google Scholar]
Phan, V. H. M., Peretti, E., Cristofari, P., Gusdorf, A., & Mertsch, P. 2024, MNRAS, 531, 2930 [NASA ADS] [CrossRef] [Google Scholar]
Polles, F. L., Salomé, P., Guillard, P., et al. 2021, A&A, 651, A13 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Rejkuba, M. 2004, A&A, 413, 903 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Rekola, R., Richer, M. G., McCall, M. L., et al. 2005, MNRAS, 361, 330 [NASA ADS] [CrossRef] [Google Scholar]
Richardson, C. T., Allen, J. T., Baldwin, J. A., Hewett, P. C., & Ferland, G. J. 2014, MNRAS, 437, 2376 [NASA ADS] [CrossRef] [Google Scholar]
Rodríguez Montero, F., Martin-Alvarez, S., Slyz, A., et al. 2024, MNRAS, 530, 3617 [CrossRef] [Google Scholar]
Rosdahl, J., Schaye, J., Teyssier, R., & Agertz, O. 2015, MNRAS, 451, 34 [CrossRef] [Google Scholar]
Ruszkowski, M., Yang, H.-Y. K., & Zweibel, E. 2017, ApJ, 834, 208 [NASA ADS] [CrossRef] [Google Scholar]
Rybicki, G. B., & Lightman, A. P. 1979, Radiative processes in astrophysics (John Wiley& Sons) [Google Scholar]
Saito, T., Takano, S., Harada, N., et al. 2022, ApJ, 927, L32 [NASA ADS] [CrossRef] [Google Scholar]
Salomé, Q., Salomé, P., Combes, F., & Hamer, S. 2016, A&A, 595, A65 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Santoro, F., Oonk, J. B. R., Morganti, R., & Oosterloo, T. 2015, A&A, 574, A89 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Schawinski, K., Thomas, D., Sarzi, M., et al. 2007, MNRAS, 382, 1415 [Google Scholar]
Scowen, P. A., Hester, J. J., Sankrit, R., et al. 1998, AJ, 116, 163 [NASA ADS] [CrossRef] [Google Scholar]
Sharp, R. G., & Bland-Hawthorn, J. 2010, ApJ, 711, 818 [Google Scholar]
Silverman, B. W. 1986, Density estimation for statistics and data analysis (London: Chapman and Hall/CRC) [Google Scholar]
Spinoglio, L., Malkan, M. A., Smith, H. A., González-Alfonso, E., & Fischer, J. 2005, ApJ, 623, 123 [NASA ADS] [CrossRef] [Google Scholar]
Spinoglio, L., Pereira-Santaella, M., Busquet, G., et al. 2012, ApJ, 758, 108 [Google Scholar]
Spinoglio, L., Fernández-Ontiveros, J. A., & Malkan, M. A. 2024, ApJ, 964, 117 [NASA ADS] [CrossRef] [Google Scholar]
Spitzer, L. J., & Tomasko, M. G. 1968, ApJ, 152, 971 [NASA ADS] [CrossRef] [Google Scholar]
Strickland, D. K., Heckman, T. M., Weaver, K. A., Hoopes, C. G., & Dahlem, M. 2002, ApJ, 568, 689 [NASA ADS] [CrossRef] [Google Scholar]
Sutherland, R. S., & Dopita, M. A. 2017, ApJS, 229, 34 [Google Scholar]
Tarter, C. B., Tucker, W. H., & Salpeter, E. E. 1969, ApJ, 156, 943 [Google Scholar]
Thomas, A. D., Groves, B. A., Sutherland, R. S., et al. 2016, ApJ, 833, 266 [NASA ADS] [CrossRef] [Google Scholar]
Tingay, S. J., & Lenc, E. 2009, AJ, 138, 808 [NASA ADS] [CrossRef] [Google Scholar]
Tingay, S. J., Jauncey, D. L., Reynolds, J. E., et al. 1998, AJ, 115, 960 [NASA ADS] [CrossRef] [Google Scholar]
Tremonti, C. A., Heckman, T. M., Kauffmann, G., et al. 2004, ApJ, 613, 898 [Google Scholar]
Tumlinson, J., & Shull, J. M. 2000, ApJ, 528, L65 [NASA ADS] [CrossRef] [Google Scholar]
Uhlig, M., Pfrommer, C., Sharma, M., et al. 2012, MNRAS, 423, 2374 [CrossRef] [Google Scholar]
Ulvestad, J. S., & Antonucci, R. R. J. 1997, ApJ, 488, 621 [NASA ADS] [CrossRef] [Google Scholar]
van Zee, L., Salzer, J. J., & Haynes, M. P. 1998, ApJ, 497, L1 [NASA ADS] [CrossRef] [Google Scholar]
Vandenbroucke, B., Wood, K., Girichidis, P., Hill, A. S., & Peters, T. 2018, MNRAS, 476, 4032 [NASA ADS] [CrossRef] [Google Scholar]
Veilleux, S., Maiolino, R., Bolatto, A. D., & Aalto, S. 2020, A&ARv., 28, 1 [CrossRef] [Google Scholar]
Venturi, G., Cresci, G., Marconi, A., et al. 2021, A&A, 648, A17 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Walker, M. A. 2016, ApJ, 818, 23 [NASA ADS] [CrossRef] [Google Scholar]
Watson, A. M., Gallagher, J. S. I., Holtzman, J. A., et al. 1996, AJ, 112, 534 [NASA ADS] [CrossRef] [Google Scholar]
Weilbacher, P. M., Streicher, O., Urrutia, T., et al. 2014, ASP Conf. Ser., 485, 451 [NASA ADS] [Google Scholar]
Weilbacher, P. M., Palsa, R., Streicher, O., et al. 2020, A&A, 641, A28 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Wolfire, M. G., Vallini, L., & Chevance, M. 2022, ARA&A, 60, 247 [NASA ADS] [CrossRef] [Google Scholar]
Yusef-Zadeh, F., Hewitt, J. W., Wardle, M., et al. 2013, ApJ, 762, 33 [Google Scholar]
Zhang, K., Yan, R., Bundy, K., et al. 2017, MNRAS, 466, 3217 [NASA ADS] [CrossRef] [Google Scholar]
Zhu, P., Kewley, L. J., & Sutherland, R. S. 2023, ApJ, 954, 175 [NASA ADS] [CrossRef] [Google Scholar]

Appendix A: Galaxies

A.1. Centaurus A

NGC 5128, commonly known as Centaurus A is a nearby Seyfert galaxy and among the closest AGN, and it is situated at a distance of ∼3.84 Mpc (Rejkuba 2004; Harris et al. 2010). In this AGN ionized, atomic, neutral, and molecular gas is oriented along the dust-lane and is observed at many wavelengths (Morganti 2010).

Two subrelativistic jets that span a few kiloparsecs have been detected in radio and X-ray observations of Centaurus A (Hardcastle et al. 2003; Kraft et al. 2008). The jet structures have been resolved down to subparsec scale by utilizing Very Long Baseline Interferometry (VLBI) and Event Horizon Telescope (EHT) observations (Tingay et al. 1998; Janssen et al. 2021). The jet’s northern side is more prominent than its southern one. From the inner lobes ∼5 kpc up to the outer lobes at 250 kpc from the nucleus, the jet expands into radio lobes (Israel 1998; Eilek 2014).

The jet and the interstellar medium (ISM) interact in a complicated way at this source (Santoro et al. 2015; Salomé et al. 2016). Aligned with the jets’ direction, there is a ∼500 pc area of ionized gas with [O III]λ5007 Å, emission (Kraft et al. 2008; Sharp & Bland-Hawthorn 2010). From far-infrared (FIR) emission lines in ionized gas, as well as in neutral atomic and molecular gas phases, an outflow is observed, spanning over ∼200 pc (Israel et al. 2017).

A.2. NGC 1068

NGC 1068 is another one of the closest and brightest active galaxies, classified as a Seyfert 2, located at a distance of ∼14.4 Mpc. NGC 1068 played a key role in the establishment of the unified model for AGN (Bland-Hawthorn et al. 1997), which reconciles the two populations of Seyfert galaxies (types 1 and 2) to a simple orientation effect of an obscuring torus, as it shows broad optical permitted lines in its polarized spectrum. It has a relatively high bolometric luminosity of L_bol ∼ 2.6 × 10⁴⁴ ergs s⁻¹ (Spinoglio et al. 2024).

From X-ray to radio wavelengths, NGC 1068 is known to contain radio jets (Bland-Hawthorn et al. 1997; Mutie et al. 2024). Its activity appears to be limited to ionization cones centered on the AGN, which cross the disk plane at an inclination of ∼45° and are roughly aligned with the inner radio jet (Gallimore et al. 1996; Bland-Hawthorn et al. 1997). This means that the nuclear radiation field immediately affects the interstellar medium of the disk situated in the north-east and south-west direction (Cecil et al. 1990; Bland-Hawthorn et al. 1997; Gallimore et al. 2004).

There are evident outflows in the ionized (Cecil et al. 1990; Barbosa et al. 2014), atomic (Saito et al. 2022) and molecular gas (García-Burillo et al. 2016; Gallimore et al. 2016; Saito et al. 2022), possibly driven by the radio jet intercepting the disk (Cecil et al. 1990; Mutie et al. 2024). Strong starburst activity is seen in the galaxy, centered on a notable starburst ring that has a radius of around 1-1.5 kpc (Spinoglio et al. 2012; García-Burillo et al. 2016).

A.3. NGC 253

The nearby starburst galaxy NGC 253 has a distance of ∼ 3.94Mpc (Rekola et al. 2005). NGC 253 is notable for massive dust areas surrounding its core and by an apparent bipolar outflow cone in both X-ray and ionized gas emissions (Strickland et al. 2002). The outflows in this starburst galaxy are of two kinds, a multiphase outflow due to stellar feedback, also known as a superwind, and a disk wind (Strickland et al. 2002; Heesen et al. 2008).

Studies conducted on NGC 253 have been over a broad range of frequencies, including polarization, HI observations and broadband radio continuum (Ulvestad & Antonucci 1997; Carilli et al. 1992; Lenc & Tingay 2006; Heesen et al. 2011; Lucero et al. 2015). In the infrared (IR) spectrum hotspots have been identified in the nuclear starburst area (Galliano et al. 2005; Beck et al. 2023). Previous radio observations detected over 60 distinct radio sources (Ulvestad & Antonucci 1997), most of which were oriented north-east—south-west and most likely connected to a 300 pc nuclear ring (Fernández-Ontiveros et al. 2009).

Thermal emission from electrons colliding with ions in the ionized interstellar medium (ISM) surrounding hot stars and nonthermal synchrotron emission from relativistic electrons traveling around the interstellar magnetic field lines are the primary causes of radio emission in starburst galaxies (Kapinska et al. 2017). SNRs accelerate CRs and are the main origin of nonthermal emission, and they subsequently create a notable synchrotron radio halo in NGC 253 (Carilli et al. 1992; Kapinska et al. 2017). A CR rate of ζ_CR ∼ 10⁻¹² s⁻¹ (greater than 10⁴ times the average Galactic rate) was recently measured in the central giant molecular clouds (GMCs) of NGC 253, based on measurements of HCN and HNC emission of rotational transitions observed with the Atacama Large Millimeter/submillimeter Array via the ALCHEMI Large Program (Behrens et al. 2022).

A.4. NGC 1320

NGC 1320 is a barred spiral galaxy situated at around ∼37.7 Mpc (Kraemer et al. 2011). Lying in the constellation of Eridanus, it shows star formation, that is primarily restricted to the narrow, inner Northwest spiral arm (De Robertis & Osterbrock 1986). It is a high-inclination galaxy that has been observed with HST in [O III]λ5007 Å and [N II]6584 Å (Ferruit et al. 2000). It also contains a low luminosity, high-ionization, Seyfert type 2 active nucleus (L_bol = 6.3 × 10⁴³ ergs s⁻¹) (Gallimore et al. 2010; Spinoglio et al. 2024). It exhibits very narrow emission lines and a faint, featureless continuum (De Robertis & Osterbrock 1986; Baloković et al. 2014).

Appendix B: BPT diagrams with density notation

Our models span between 0 ≤ log n_H ≤ 4 from white being the smaller initial hydrogen density to forest green being the highest (Figs. B.1, B.2, and B.3), and we can also see the range of the same models over the ionization parameter −3.5 ≤ logU ≤ −1.5 from white being the smaller ionization parameter to deep red being the highest (Figs. 5, 6, and 7). The same diagrams are produced for different CR ionization rates −14 ≤ log ζ_CR ≤ −12 for all galaxies.

Fig. B.1.

BPT diagrams with the AGN photoionization models compared with the observations from the selected apertures in Centaurus A (Fig. 1a). 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 going from deep purple to pale lilac/white represent the ascending distance from the nucleus, as also noted with numbers, with "N" being the closest aperture. Also from white to deep green, the different shades of green represent the range of AGN models’ densities 1 ≤ n_H ≤ 10⁴ cm⁻³. All the models shown have solar abundances. The red diamonds represent the measured line ratios for the photoionization-dominated Seyfert 2 nucleus in NGC 1320. The Kewley, Kauffmann, and Schawinski lines correspond to the black, blue, and red solid lines, respectively.

Fig. B.2.

BPT diagrams with the AGN photoionization models compared with the observations from the selected apertures in NGC 1068 (Fig. 1b). 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 going from deep purple to pale lilac/white represent the ascending distance from the nucleus, as also noted with numbers, with "N" being the closest aperture. Also, from white to deep green, the different shades of green represent the range of AGN models’ densities 1 ≤ n_H ≤ 10⁴cm⁻³. All the models shown have solar abundances. The red diamonds represent the measured line ratios for the photoionization-dominated Seyfert 2 nucleus in NGC 1320. The Kewley, Kauffmann, and Schawinski lines correspond to the black, blue, and red solid lines, respectively.

Fig. B.3.

BPT diagrams with the SF photoionization models compared with the observed line ratios from the selected apertures in NGC 253 (Fig. 1c). 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 going from deep purple to pale lilac/white represent the ascending distance, as also noted with numbers, with "1" being the most central aperture. Also from white to deep green, the different shades of green represent the range of SF models’ densities 1 ≤ n_H ≤ 10⁴ cm⁻³. All the models shown have solar abundances. The red diamonds represent the measured line ratios for the photoionization-dominated Seyfert 2 nucleus in NGC 1320. The Kewley, Kauffmann, and Schawinski lines correspond to the black, blue, and red solid lines, respectively.

All Tables

Table 1.

Observational parameters for the selected galaxies, including coordinates, redshift values, z, and redshift-independent distances, D, provided by the NASA/IPAC Extragalactic Database.

	Fig. 2. BPT emission-line fits, using Pyplatefit (Bacon et al. 2023), in the rest frame of Centaurus A. (a) Hα. (b) [NII]λ6584 Å. (c) [S II]λλ6716,6731 Å. (d) [O I]λ6300 Å. (e) [O III]λ5007 Å.
In the text

	Fig. 3. Hβ emission line in each galaxy’s rest frame, fitted with Pyplatefit (Bacon et al. 2023). (a) Centaurus A. (b) NGC 1068. (c) NGC 253.
In the text

	Fig. 4. AGN and BB SEDs simulated with CLOUDY. The blue solid line represents the AGN SED used in the AGN models, while the green solid line represents the BB SED used in the SF models. The y-axis is νF_ν in arbitrary scaling.
In the text

Cosmic-ray ionization of low-excitation lines in active galactic nuclei and starburst galaxies

1. Introduction

2. Observational data

2.1. Galaxy sample

2.2. MUSE data

3. Methods

3.1. Spectrum extraction

3.2. Constraining the parameter space

3.3. CLOUDY models

4. Results

4.1. BPT diagrams

4.2. Gas stratification diagrams

5. Discussion

5.1. Thermal stability at high CR ionization rates

5.2. CR heating and low ionization nebular lines

5.3. Reconciling models with the observed metallicities

5.4. CR effect on metallicity and ionization parameter

5.5. Comparison with CR rate estimates

5.6. A new SFζ line for starburst galaxies with CRs

5.7. Further remarks and considerations

6. Summary

Acknowledgments

References

Appendix A: Galaxies

A.1. Centaurus A

A.2. NGC 1068

A.3. NGC 253

A.4. NGC 1320

Appendix B: BPT diagrams with density notation

All Tables

All Figures