Cosmic-ray induced ionization rates and non-thermal emissions from nuclei of starburst galaxies

Cosmic rays are the only agent capable of ionizing the interior of dense molecular clouds and, thus, they are believed to play an essential role in determining the physical and chemical evolution of star-forming regions. In this work, we aim to study cosmic-ray induced ionization rates in starburst environments using non-thermal emissions of cosmic rays from starburst nuclei. To this end, we first revisit cosmic-ray models which could explain data of non-thermal emissions from radio to X-ray and gamma-ray from nuclei of three prototypical starburst galaxies NGC 253, M82, and Arp 220. These models are then applied to predict ionization rates in starburst environments which gives values around $10^{-14}$ s$^{-1}$. Such a high value of the ionization rate, which is 2 to 3 orders of magnitude higher than the typical values found in the Milky Way, is probably due to relatively high rates of supernova explosions occurring within the nuclei of these starburst galaxies. We also discuss in more details the case of NGC 253 where our predicted ionization rate is found to be, in most cases, a few times smaller than the values inferred from molecular line observations of clouds in the starburst nucleus. The general framework provided in this work illustrates how the use of non-thermal emission data could help to provide more insights into ionization rates or, more generally, cosmic-ray impact in starburst environments.


INTRODUCTION
Starburst galaxies are galaxies with high star formation rates typically ranging from 10 to 10 3 times higher than that of the Milky Way (Gao & Solomon 2004).The intense star-forming activity also results in a high rate of supernova explosions and, as massive stars and supernova remnants (SNRs) are commonly believed to be cosmicray sources (Strong et al. 2007;Grenier et al. 2015;Gabici et al. 2019;Cristofari 2021), starburst galaxies are expected to be filled with cosmic rays (CRs).The connection between CRs and the star forming activity in starburst galaxies is further supported by tight correlations between the inferred star formation rate and the nonthermal luminosity, ascribed to CRs (see e.g.Lacki et al. 2010;Kornecki et al. 2020;Ajello et al. 2020;Kornecki et al. 2022, and references therein).Of particular interest are the starburst nuclei (SBNi, see e.g.Westmoquette et al. 2009) which have sizes of about a few hundred parsecs but within these small region exhibit rates of supernova explosions comparable to or even higher than that of the entire Milky Way.It is for this reason that SBNi are considered ideal laboratories to study CR impact on star-forming regions.In these environments, the density of gas, radiation and magnetic field are ★ vhmphan@obspm.fr† enrico.peretti.science@gmail.cominferred to be at least 10 2 times larger than the average interstellar medium (ISM) of the Milky Way.This implies that most of the injected non-thermal particles are expected to lose their energy before being able to escape (Yoast-Hull et al. 2013;Peretti et al. 2019).Such a calorimetric behavior for CRs, together with the enhanced star formation rate inferred at redshift 1-4 (Madau & Dickinson 2014) suggested that starbursts could be an ideal source class to substantially contribute to the diffuse flux of high-energy gammarays and neutrinos (Tamborra et al. 2014;Bechtol et al. 2017;Peretti et al. 2020;Ambrosone et al. 2021;Roth et al. 2021;Owen et al. 2022;Peretti et al. 2022;Condorelli et al. 2023).
It has long been suggested that CRs can play an essential role in setting the chemistry and even dynamics of star-forming regions (Padovani et al. 2020;Gabici 2022).This is because these particles could penetrate deep inside dense molecular clouds, where X-rays and UV photons cannot penetrate, to ionize the interior of these objects (e.g.Padovani et al. 2009;Ivlev et al. 2018;Phan et al. 2018;Owen et al. 2021).In other words, CRs control the ionization rate which is a key parameter in regulating the abundances of different chemical species in molecular clouds.The ionization rate also determines the coupling between gas and magnetic fields which is of critical importance for the process of star formation (Krumholz & Federrath 2019;Girichidis et al. 2020).Thus, the impact of CRs on star-forming regions can be partially quantified using the cosmic-ray ionization rate.
In fact, there are a few different variants for the definition of the ionization rate (Neufeld & Wolfire 2017).Here, we refer to the ionization rate as the production rate of H + 2 per hydrogen molecule which will be denoted as  (H 2 ).The observational determination of  (H 2 ) in the interstellar medium is a longstanding problem.It is usually solved by identifying a set of species whose abundance or abundance ratio is sensitive to it.As such, the method is intrinsically subject to various limitations.First, from the observational point of view, a sufficient number of lines and an appropriate radiative transfer treatment must be implemented in order to infer reliable column densities.Second, the model that is used to predict abundances or abundance ratios must contain all the chemical pathways relevant to the destruction and formation of the observed species, and all the excitation mechanisms (shocks, energetic photons in particular ionizing UV ones, and CRs) relevant to the observed regions.Last limitation, identifying which species is tracing which mechanism requires the computing of large grids of models purposefully covering all input parameters.Within this framework, various species have been found to allow for the determination of the cosmic ray ionization state in various environments within the limits of carefully spelt-out assumptions.Several reviews focus on the determination of the ionization rate in the dense and diffuse medium, like Indriolo & McCall (2013), Neufeld & Wolfire (2017), and Barger & Garrod (2020).In quiescent and dense regions completely shielded from dissociating UV radiation, where all hydrogen is in molecules, Guelin et al. (1977) and Wootten et al. (1979) analytically showed that the DCO + to HCO + abundance ratio can be used to measure the ionization rate in steadystate conditions.Their result was somehow confirmed by the use of a more complex chemical code with the same assumptions and applied to then-state-of-the-art observations by Caselli et al. (1998).In the diffuse molecular medium with electron abundance between 10 −7 and 10 −2 , H + 3 has been extensively used to measure the ionization rate, based on a simplistic description of its chemistry (e.g.McCall et al. 2003, Indriolo et al. 2007, Indriolo & McCall 2012).In fact, the modelling works of Le Petit et al. (2004) mitigated the conclusions by highlighting that constraining the ionization rate can not be done independently from other species.In other words, they showed that self-consistent models should be used to reproduce the abundances of both H + 3 and a number of atomic and molecular species in order to effectively provide a measure of the ionization rate.Le Petit et al. (2016) then detailed all the dependences of the H + 3 abundance to parameters of such self-consistent models and used them to measure the ionization rate in the central molecular zone, adding complementary constraints on hydride abundances.Comprehensive modelling was also the path chosen by Neufeld & Wolfire (2017) to provide constraints on the ionization rate in diffuse (both atomic and molecular) clouds of the Galactic disk.
In the NGC253 starburst galaxy, Holdship et al. (2022) (hereafter H22) used ALMA observations to constrain the ionization rate to between 10 −14 s −1 and 8 × 10 −13 s −1 , using rather strong assumptions.Their observations were performed at 1. s 6 angular resolution, that is about 30 pc spatial resolution, for which they used a 'single point model' with no photo-processes included.They ruled out the possibility that shocks or UV radiation could play a role in the heating or in the modelling of the column densities of SO and H 3 O + .Their justification for the shock process omission was that parametric, non-self-consistent shock models with high pre-shock density predicted values of the abundance ratio that they did not find in the observations.Their observed abundance ratio was obtained from the non-local thermodynamic equilibrium code RADEX modelling of 3 lines of H 3 O + and 37 lines of SO (see van der Tak et al. 2007, for more details on RADEX).For both species, collisions with other particles than H 2 were not considered, and for H 3 O + , the three transitions they used have the upper energy level between 79.5 and 169.1 K only, possibly biasing the results towards a hot component.Additionnally, Behrens et al. (2022) (from here on referred to as B22) used HNC and HCN observations from the same dataset as H22 to constrain the ionization rate to within the range from 10 −13 s −1 to 10 −12 s −1 .They used the same kind of chemical and radiative transfer modelling as H22, also neglecting UV radiation and shock processes.Given all these assumptions, related to the observations, radiative transfer and physico-chemical modelling, plus additional ones discussed in these two articles, we consider these ionization rate data as indications rather than definitive values.These results are, however, encouraging as they point to high values of  (H 2 ) which are qualitatively expected given the high supernova rate in the SBN of NGC 253.The relatively large difference in results of these similar analyses calls for a different approach to determine the ionization rate in these complex environments.In this paper, we will discuss a rather different approach to predict ionization rates in SBNi which relies on non-thermal emissions from these objects.
Many SBNi are bright gamma-ray sources both in the GeV and TeV energy range and some of them are also detected with X-ray telescopes (VERITAS Collaboration et al. 2009;Acero et al. 2009;Ackermann et al. 2012;Fleischhack & VERITAS Collaboration 2015).Several prototypical starburst galaxies, for example NGC 253, M82, or Arp 220, have been observed by both satellite and ground-based gamma-ray telescopes revealing their gamma-ray spectrum extending from a few hundred MeVs to about 10 TeV (see e.g.Ajello et al. 2020;Tibaldo et al. 2021, and references therein).The gamma-ray emission is likely dominated by the decay of neutral pions produced in interactions between CR protons and interstellar gas in SBNi independent on the transport conditions (Peretti et al. 2019;Krumholz et al. 2020).This means that the gamma-ray spectrum could be employed to extract the CR proton spectrum within these SBNi.In addition, CR electrons in these systems can also induce detectable emissions in the range extending from radio (with frequency around 1 GHz) to X-ray (around 1 keV) via synchrotron radiation.Observations in the X-ray domain could, however, be contaminated by unresolved sources such as X-ray binaries (Wik et al. 2014;Strickland & Heckman 2007;Paggi et al. 2017) and part of the radio emissions can also come CR electrons confined in a larger halo surrounding the SBN (Yoast-Hull et al. 2013).Nevertheless, the radio and X-ray spectrum can be used as an upper limit to constrain the CR electron spectrum.The combination of these non-thermal emissions from radio to X-ray and gamma-ray provides a powerful tool to study CRs in starburst galaxies and, ultimately, allows us to quantify the impact of these particles in the star forming activity of these systems.The paper will be structured as follows.In Section 2, we will introduce the transport model for CRs in SBNi which essentially provides the CR spectra used to study non-thermal emissions and also evaluate the ionization rates.The relevant radiative and ionization processes are introduced in Section 3. We then perform a fit of the transport model to non-thermal emission data from the nuclei of NGC 253, M82, and Arp 220 to derive the CR spectra in these systems which can be employed to predict ionization rates  (H 2 ) for molecular clouds of different column densities.Our predicted values of  (H 2 ) for NGC 253 seem to be, in most cases, a few times smaller than the values recently inferred by H22 and B22 using molecular line observations (see Section 4).We also discuss in Section 5 the potential implications of this discrepancy.

COSMIC-RAY TRANSPORT IN STARBURST NUCLEI
We will follow the approach presented in Peretti et al. (2019) and adopt the following transport equation to describe the differential number density   () for CRs of species  = p (protons) or  = e (electrons) with kinetic energy  in SBNi where   () is the injection spectrum of CRs from SNRs (or as secondary and tertiary products from CR interactions with the ISM of SBNi),   () is the energy loss rate due to interactions of CRs with SBNi's materials,  adv and  diff, () are respectively timescales for the escape of CRs from the SBNi due to advection and diffusion.The advection timescale is simply  adv = /  where  is the radius of the SBNi and   is the speed of galactic winds in SBNi.The diffusion timescale could be estimated as  diff, () =  2 /  () where   is the diffusion coefficient of CRs in SBNi.Here, we adopt Model A from Peretti et al. (2019) which assumes the diffusive motion of particles to be induced by magnetic turbulence following a Kolmogorov power spectrum and the diffusion coefficient scales as where  0 = 1 pc is the injection scale of turbulence,   =  2 / 2 = 1 is the ratio between the variance of turbulent magnetic fields  2 and the ordered field strength squared  2 ,   is the speed of CRs of species  with kinetic energy , and  , () is the Larmor radius of CRs of species  with kinetic energy .A complete list of energy loss processes for both protons and electrons can be found in Schlickeiser (2002) (see also Evoli et al. 2017, andPeretti et al. 2019).We note that the energy loss rates depend on many parameters characterizing the ISM of SBNi including magnetic field strength , ISM density  ISM , electron density  e , electron temperature  e , and interstellar radiation field (see Eq. 14 below).
Concerning the injection spectrum, we consider CR protons injected only from SNRs and their injection spectrum could be modelled as a power law in momentum where R SNR is the rate of supernova explosions within the SBNi,  CR,p is the fraction of supernova explosion kinetic energy converted into CR kinetic energy (also referred to as the acceleration efficiency),  SNR ≃ 10 51 erg is the typical kinetic energy of supernova explosions,  max,p ≃ 10 17 eV/c is the cut-off momentum for CR protons accelerated from SNRs,  p is proton mass,  is the speed of light,  = 4 3 /3 is the volume of the SBNi, and Note that the normalization factor Λ ensures that a fraction  CR,p of supernova explosion kinetic energy is converted into CR kinetic energy, meaning Throughout this work, we will assume  CR,p = 0.1 which is, roughly speaking, around the typical value adopted for Galactic SNRs in models attempting to fit CR spectra in the local ISM (see e.g.Evoli et al. 2019;Phan et al. 2021;Cristofari 2021;Mertsch et al. 2021).
For CR electrons, three different types of injection are taken into account where  SNR,e () is the injection spectrum of CR electrons from SNRs,  sec () is the injection spectrum of secondary electrons (and positrons) from the decay of  ± produced in proton-proton interactions, and  ter () is the injection spectrum of tertiary electrons (and positrons) created in interactions between CR-induced gamma-rays and low-energy photons in the ISM of SBNi.In the following, we assume also a power law injection spectrum for CR electrons from SNRs  SNR,e () ∝  CR,e  2−  exp(−/ max,e ) with  CR,e = 0.01 and  max,e ≃ 10 13 eV/c.For secondary electrons, we adopt the approach presented in Kelner et al. (2006) to model the injection spectrum of secondary electrons as follows where   ≃ 0.17 is the fraction of kinetic energy transferred from the parent proton to the single pion,  pp () is the total inelastic crosssection for interactions between CR protons with kinetic energy  and protons in the ISM of SBNi (Kafexhiu et al. 2014), and fe (/  ) is defined as in Kelner et al. (2006) (see also Appendix B of Peretti et al. 2019).Concerning tertiary electrons, the injection spectrum is where  nth () is the differential number density of non-thermal photons induced by CRs and   () is the opacity of gamma-rays due to interactions with low-energy photons in the ISM of SBNi which could be estimated as follows Here, we have introduced the volume emissivities for four main processes for non-thermal emissions induced by CRs  PPI ,  BRE ,  ICS , and  SYN which are respectively proton-proton interactions, bremsstrahlung radiation, inverse Compton scattering, and synchrotron radiation.These processes will be discussed in more detail in the next section.In principle, the source spectrum of tertiary electrons should depend on bremsstrahlung radiation and inverse Compton scattering induced by CR electrons themselves which makes the problem non-linear.We shall see later that the observed gamma-ray spectrum of SBNi might be dominated by the hadronic gammaray component and, thus, one can neglect the non-linearity in the CR transport equation for electrons by considering tertiary electrons coming from hadronic gamma-rays only.As for the opacity of gamma-rays in SBNi, we have introduced also the gamma-gamma interaction cross-section   (,  ph ) (Aharonian et al. 2013) and the differential number density of interstellar photons  ISRF ( ph ) which shall be introduced in more details in the next section (see Eq.

14).
Having discussed all the relevant ingredients, we could essentially proceed to the solution of the transport equation applicable for both CR protons and electrons .( 10) Note that the CR differential number density could be related to the CR spectra or CR flux as follows   () =   ()  /(4).

Cosmic-ray induced gamma-rays and X-rays
As mentioned in the previous section, we will focus mostly on CR induced gamma-rays from the decay of  0 in proton-proton interactions, bremsstrahlung radiation, and inverse Compton scattering.In the following, we shall briefly discuss the volume emissivities for these processes (interested readers could find some more details in Peretti et al. 2019).
Concerning  0 decay, we shall follow the approach as presented in Kelner et al. (2006) but adopt the differential cross-section for proton-proton interactions d pp /d  and the nuclear enhancement factor  n (to correct for gamma-rays induced by CR nuclei) from Kafexhiu et al. (2014).The volume emissivity could be modelled as follows Note that the differential cross-section will be non-zero only for CR protons with kinetic energy satisfying  ≤   +  2   4 /(4  ).As for bremsstrahlung radiation, the volume emissivity could be estimated following Schlickeiser (2002) where the cross-section for bremsstrahlung radiation  BRE (,   ) (see Chapter 4 of Schlickeiser 2002 or Baring et al. 1999 for more details).
Another important process for gamma-rays induced by CR electrons is inverse Compton scattering where low-energy photons in the ISM are scattered by CR electrons and are boosted to higher energy.Modelling this process, thus, requires some knowledge of the lowenergy interstellar photons especially in the far infrared to optical energy range.We shall follow Peretti et al. (2019) and consider the interstellar radiation field made up mostly of four components: farinfrared (FIR), mid-infrared (MIR), near-infrared (NIR), and optical (OPT).The differential number density of each component could be modelled as follows where rad = FIR, MIR, NIR, or OPT,  rad is the energy density of the respective component,  rad is the effective photon temperature of the respective component, Γ() and  () are respective the Gamma and Riemann zeta functions,  is a spectral index which shall be set to be  = 0 for OPT and  ≃ 1.3 for the rest.In fact, the differential number density would become that of black-body radiation for the case where  = 0. Note also that we have used the notation   to indicate the photon energy for both the thermal and non-thermal photons in all energy domains (from radio to gamma-ray).The differential number density of interstellar photons is then written as follows where the sum is performed over all the components of the interstellar radiation field as mentioned above.The volume emissivity of inverse Compton scattering is then (Blumenthal & Gould 1970) where  ICS (,   ) is the cross section for the inverse Compton scattering (see Chapter 4 of Schlickeiser 2002 for more details).It is worth mentioning also that this process also has a threshold energy meaning that the emissivity is non-zero only for CR electrons with kinetic energy  ≳   1 + √︃ 1 +  2 e  4 /(   B  rad ) .In starburst galaxies, non-thermal radio and X-ray emissions can also come from synchrotron radiation induced by CR electrons.The volume emissivity of synchrotron radiation could be written as (Zirakashvili & Aharonian 2007) where  is the field strength of the magnetic field in the starburst nucleus, ℎ is the Planck constant, 2    .The flux of photons (for both thermal and non-thermal emissions) from these galaxies is then given as where  nth (  ) and  ISRF (  ) are respectively defined in Eq. 8 and Eq.14,  gal is the distance between the galaxy of interest and the Milky Way,  ff (  ) is the opacity due to free-free absorption relevant in the radio domain (see Appendix A), and   (  ) and  EBL ( SBN ,   ) are the opacities due to interactions of high-energy photons respectively with ISRF of the SBN and with the extragalactic background light relevant in the gamma-ray domain.We have defined   (  ) in Eq. 9 and adopt the analytic form of  EBL (  ) as presented in Appendix C of Peretti et al. (2019).Note that the factor 1/3 is due to the effective spherical shape assumed for the SBNi.The thernal and non-thermal flux shall be fitted later with observational data to obtain the source and transport parameters of CRs in SBNi.

Cosmic-ray induced ionization rates
Once the parameters determining CR spectra in SBNi have been fitted using gamma-ray and X-ray data, we could proceed to predict the ionization rate in molecular clouds that are embedded within these systems.We note, however, that the CR-induced ionization rate should actually vary for clouds with different column densities.Selfconsistent predictions of this quantity require a model to describe the transport of CRs into molecular clouds.In fact, CR transport In SBNi, the value of   is, in fact, not very well known.However, if we assume that the large scale magnetic turbulence in these systems are also generated by supernova explosions similar to that in the Milky Way (Berezinskii et al. 1990;Blasi 2013;Evoli et al. 2018), then the value of   should be comparable to the typical sizes of SNRs.We could take as a rough estimate for   the size of an SNR at the end of the Sedov-Taylor phase which, for a typical density of around 200 cm −3 in SBNi, is about 6 pc (assuming a typical ejecta mass of around 1  ⊙ ).Such a coherence length should be comparable to size of clouds in the SBNi of interest.It is for this reason that we shall adopt the ballistic model where the transport of CRs inside clouds is one-dimensional and the average CR spectra in a cloud of size   could be estimated as follows In fact, it can be shown that the average differential number density as defined in Eq. 18 depends only on the total column density of the cloud which, for dense molecular clouds, is     ≃ 2 (H 2 ).Thus, given the CR flux in the ISM of SBNi, we could predict the ionization rate inside clouds as a function of the H 2 column density  (H 2 ).
The H 2 ionization rate induced by CR protons and electrons could be obtained as in Padovani et al. (2009) (see also Chabot 2016;Phan et al. 2018;Recchia et al. 2019): where   ion is the ionization cross-section of CR species ,   () are the average secondary ionization per primary ionization computed as in Krause et al. (2015) (see also Padovani et al. 2009;Ivlev et al. 2021), and  (H 2 ) = 15.603eV is the ionization potential of H 2 .It should be noticed that, following Krause et al. (2015), the two ionization cross-sections are considered in the computation with fully relativistic corrections.The total CR-induced ionization rate is then It is worth mentioning also that we will ignore the ionization rate induced by X-rays.In fact, X-rays with energy from ∼ 1 to 10 keV from synchrotron radiation and inverse Compton scattering of CR electrons in nuclei of starburst galaxies could also contribute to the ionization rate but we have checked that  X−ray (H 2 ) ≲ 10 −17 s −1 for clouds with  (H 2 ) ≳ 10 23 cm −2 .
Having discussed all the relevant radiative processes in the previous subsection, we could now apply them to derive the CR spectra in several prototypical nearby SBNi, namely NGC 253, M82, and Arp 220 by fitting non-thermal emission data from these objects.This should then allow us to predict ionization rates within these SBNi using the ballistic model.

Starburst nucleus of NGC 253
NGC 253 is a spiral galaxy located at a distance of  gal ≃ 3.8 Mpc (Rekola et al. 2005;Dalcanton et al. 2009) making it one of the closest objects to be classified as a starburst galaxy in the southern sky.This starburst galaxy is believed to have a star formation rate (SFR) of about 5  ⊙ yr −1 which is a few times higher than that of the Milky Way.About 70% of the star-forming activity occurs, however, in the starburst nucleus region (Melo et al. 2002).As a consequence, the SBN of NGC 253, which is of size  ≃ 100 pc, has a relatively high supernova rate with R SNR ≃ 0.03 yr −1 comparable to that of the entire Milky Way (Engelbracht et al. 1998).
Interestingly, the nucleus of NGC 253 has been observed in gamma-rays both in the GeV and TeV energy ranges with Fermi-LAT and H.E.S.S. telescopes (Abdo et al. 2010; H. E. S. S. Collaboration et al. 2018, see also Abramowski et al. 2012).There exist also upper limits in the energy range around 10 keV by the X-ray telescope NuSTAR (Wik et al. 2014).Recently, several molecular clouds in the central region of this SBN have also been targeted to study ionization rates by the ALCHEMI Collaboration using chemical surveys in the SBN of NGC 253 obtained using the source and transport parameters adopted for the fit of non-thermal emissions.Data of CR spectra observed in the local ISM is also presented for comparison.For CR electrons, we present different components including primary electrons from SNRs (orange line), secondary electrons from decays of  ± (green line), and tertiary electrons from interactions between CR-induced gamma-rays and interstellar radiation (blue line).
performed with ALMA data (H22; B22).We shall now estimate the ionization rates in the SBN of NGC 253 using non-thermal emissions and compare them to the ones derived from molecular line observations.
All the parameters relevant for modelling both thermal and nonthermal emissions from the SBN are presented in Table 1.We have fixed the distance, the SBN size, and the galactic wind speed as in Peretti et al. (2019) which are also motivated from several independent observations.The remaining parameters are fitted to data and upper limits from radio to gamma-ray and the results are shown in Fig. 1.The parameter choice leads to the expected flux being dominated by gamma-rays from  0 decay above about 100 MeV.Gamma-rays and X-rays from bremsstrahlung, inverse Compton scattering, and synchrotron radiation induced by CR electrons, on the other hand, contribute significantly to the flux below about 100 MeV with the expected flux consistent with the upper limit set by NuSTAR.Interestingly, the choice of the galactic wind speed and the diffusion coefficient as in Eq. 2 mean that the transport of CRs in this system is likely dominated by energy loss (calorimetric system) given a typical ISM density  ISM ≳ 100 cm −3 .In this case, the solution of the transport equation (see Eq. 10) simplifies to  () ∼ ()/() ∼  CR,p R SNR 2−  /( ISM  3 ) for  ≳ 300 GeV.If we notice also that the hadronic gamma-ray flux at energy   is also, roughly speaking proportional to the CR proton spectrum at  ≃ 10  , the SBN gas mass and the inverse of distance squared, we can show that  CR,p R SNR  2−   2 ∼ (  = /10) (see e.g.Abramowski et al. 2012 for similar discussions).Such a relation means that, given fixed values of  CR,p and , the index of the CR injection spectrum  and the supernova rate R SNR can be constrained respectively by the spectral index and the normalization of the high-energy gamma-ray spectrum.More importantly, the exact shape of the predicted gamma-ray spectrum around GeV energy is determined by the values of the ISM density  ISM .For NGC 253, we have found that  = 4.3, R SNR = 0.03 yr −1 , and  ISM = 170 cm −3 .The value of R SNR = 0.03 yr −1 is, in fact, consistent with values derived from Engelbracht et al. (1998) using spectroscopic data.We note, however, that uncertainties on gamma-ray data mean that this value can be uncertain within a factor of two.Also, the gas mass of the SBN derived from  ISM and  is about 7 × 10 7  ⊙ which is within the uncertainty range indicated by other estimates using molecular line observations (Bradford et al. 2003).
We present also the fit results in the frequency range from radio to optical in Fig. A1 of Appendix A with data from various observations retrieved from the NASA/IPAC Extragalactic Database 1 .It is clear that synchrotron radiation becomes relevant in the frequency range below a few tens GHz.Here, the magnetic field has been fitted to  = 120 G which is comparable to the magnetic field obtained by assuming equipartition of energy density between CRs and magnetic field and, thus, quite conservative as a lower limit for the values of .The fitted magnetic field strength also ensure that the synchrotron radiation from secondary CR electrons (fixed by the CR protons and the ISM density) do not surpass upper limits derived by Williams & Bower (2010) using data from the Allen Telescope Array around 1 GHz.In addition, we have chosen the electron acceleration efficiency  CR,e = 0.01 as commonly adopted in studying Galactic CRs and this lead to a subdominant contribution of primary CR electrons in the radio domain compatible with the available upper limits.The tight constraints in the GHz domain, in fact, leave little room for increasing the value of  CR,e much above the Galactic value.
The corresponding CR spectra derived from these fit parameters are also shown in Fig. 2 with CR data from the local ISM overlaid for comparison.It is clear the CR spectra in the SBN of NGC 253 are many orders of magnitude larger than that in the local ISM, especially in the MeV to GeV energy range.Thus, we expect the CR-induced ionization rate to be also much larger than typical Galactic values.
Using parameters obtained in the fit of non-thermal emissions, we could now predict the CR-induced ionization rate in the SBN of NGC 253.In this case, the ionization rate could reach about  (H 2 ) ≃ 1.5 × 10 −14 s −1 for clouds with column density  (H 2 ) = 10 23 cm −3 .The ionization rate for clouds of different column densities is shown in Fig. 3 together with the inferred values of ionization rates in several giant molecular clouds from H22 and B22.Separate contributions of CR protons and electrons to the total ionization rate in the SBN of NGC 253 are also shown in the lower panel of Fig. 3. Interestingly, the contributions from CR protons are always dominant over that of electrons.This result seems rather conservative given the tight constraints of CR electrons in the radio domain.It is clear the predicted ionization rate is lower than the values inferred from observations of molecular lines.The differences in most cases are, however, only a factor of a few to roughly one order of magnitude below the data points, except for two extreme cases from B22 where the ionization rates reach a value of a few 10 −12 s −1 .There might be many potential explanations for such a discrepancy which will be elaborated in Section 5.At this point, we would like however to provide a short discussion on ionization rate data derived using chemical modeling for molecular line observations in comparison to our predicted ionization rate.
Indeed, the cosmic ray ionization rates obtained by these studies are actually derived by fitting, or more precisely performing Bayesian inference for, a chemical model with a small number of parameters, including also  (H 2 ) and  (H 2 ), to molecular line observations.This gives, in the end, a posterior distribution of all the parameters which, ultimately, allow us to quantity the values of  (H 2 ) and  (H 2 ) together with their uncertainties.This procedure are adopted for both H22 and B22 and, in fact, the two analyses are performed for the same set of molecular clouds but study emissions from different molecules: H22 focus on H 3 O + and SO and B22 examine HCN and HNC.In other words, each green data point in the upper panel of Fig. 3 has a corresponding yellow data point and yet their values are different by a factor of a few to roughly one order of magnitude in most cases.This could mean that the quoted error bars from these data points might not fully reflect the uncertainties in chemical modelling adopted to derive the ionization rates.This is likely due to uncertainties intrinsic to: the observations (performed over clouds of almost 30 pc typical size, and sometimes possibly biased towards warm regions which is the case for H 3 O + , for which the lowerlying transitions are basically excluded from the analysis because satisfying fits cannnot be found), the radiative transfer of the species (with partial collisional processes implemented, neglecting collisions with electrons for instance), and decisions on chemical modelling (using a single point model for each cloud, excluding some reactions with unknown rates, approximating the sulfur depletion, and using equilibrium values) as well as on physical modelling (excluding the treatment of shocks or UV photons).In order to better illustrate these uncertainties, we provide in the lower panel of Fig. 3 an example where our predictions of the ionization rates are overlaid with the 1 and 2 contours for a particular molecular cloud, referred to as GMC 6 by both H22 and B22 (see Table A1 of B22 for the sky coordinate of this cloud).For the case of H22, our predicted ionization rates are actually within the 2 contour for this clouds if  (H 2 ) ∼ ×10 23 cm −2 .We have also checked that, roughly speaking, this is also true for all the data points from H22. Regarding the case of B22, the comparison is more complicated as there are, in fact, multiple regions of the plane  (H 2 ) - (H 2 ) where the chemical model gives good fit to molecular line data 2 .There are also regions of the posterior giving lower values for  (H 2 ) and  (H 2 ) which are less likely but compatible with our predictions for clouds with low column densities.
Nevertheless, we have illustrated that the ionization rate estimated from non-thermal emissions can be in most cases within a factor of a few different than that derived by molecular line observations.Thus, it can be used in complementarity to molecular line observations to provide more precise values for ionization rates which might be useful for chemical modeling of complex star-forming regions.

Starburst nuclei of M82 and Arp 220
We can also apply the framework presented above to study the ionization rate in SBNi of M82 and Arp 220.These two starburst galaxies are also relatively nearby with high SFRs (and correspondingly high supernova rates) such that they are also visible by gamma-ray telescopes in the GeV and TeV energy ranges and by the X-ray telescope Chandra in the keV energy range.We will again model the underlying CR spectra which could account for the non-thermal emissions and employ these spectra to predict ionization rates in the SBNi.
Concerning M82, this is a nearby starburst galaxy (about 3.9 Mpc from the Milky Way, Sakai & Madore 1999) and it is very well known for hosting a galactic superwind with a wind speed reaching several hundred kilometers per second (Strickland & Heckman 2009).This galaxy also has an active compact starburst nucleus (of size  ≃ 150 pc, Völk et al. 1996) which is believed to form due to its interactions with the nearby spiral galaxy M81 (Yun et al. 1994).The SBN has been inferred to have an SFR which is about 10 times higher than that of the Milky Way and the corresponding supernova rate in this compact starburst region is R SNR ≃ 0.1-0.3yr −1 (Kronberg et al. 1985;Fenech et al. 2008).
Since the distance, size, and supernova rate of this SBN are very similar to that of NGC 253, the nucleus of M82 is also expected to be a bright gamma-ray and X-ray source.In fact, it has been observed also in the GeV energy range by Fermi-LAT (Acero et al. 2015) and in the TeV energy range by VERITAS (VERITAS Collaboration et al. 2009).Observations with several telescopes including Chandra, XMM-Newton, and NuSTAR also reveal a very high flux of hard Xray around 10 keV from the M82 SBN (Strickland & Heckman 2007;Ranalli et al. 2008;Bachetti et al. 2014).However, it has been argued by Peretti et al. (2019) that there might exist unresolved X-ray sources in the SBN and, thus, the observed X-ray flux should be treated as an upper limit for the expected flux of non-thermal emissions induced by CRs.
The resulting fit of the non-thermal emissions is presented in Fig. 4 with the fit parameters as reported in Table 1.We could also see that gamma-rays of energy above about 100 MeV are mostly induced by CR protons.Below about 100 MeV, leptonic processes namely bremsstrahlung radiation, inverse Compton scattering, and synchrotron radiation start to dominate the non-thermal emissions.As before, we adopt these fit parameters to predict the ionization rate for the SBN of M82 (see Fig. 3) which could reach  (H 2 ) ≃ 6.9 × 10 −15 s −1 for clouds with  (H 2 ) = 10 23 cm −2 .The ionization rate in this case is about 2 times lower than that of NGC 253 even though the nucleus of M82 has a higher supernova rate.This is because the radius of the SBN of M82 is slightly larger and also it has a larger wind speed which results in particles escaping the nucleus more quickly via advection.
Another interesting system to be considered is the galaxy Arp 220 which is located at a distance of  gal ≃ 77 Mpc (Scoville et al. 1998).In fact, this galaxy is a merger of two galaxies and, thus, it has two dense nuclei which are about a few hundred parsecs apart from each other.We shall follow Peretti et al. (2019) and treat the Arp 220 SBN approximately as one nucleus with an effective size  ≃ 165 pc.Interestingly, the nucleus of this galaxy has been observed at a few different wavelengths in the radio domain and the spectral analysis of several detected sources allow us to infer the supernova rate of R SNR ≃ 2-6 yr −1 .We fit the non-thermal spectrum of the Arp 220 nucleus to X-ray data from Chandra (Paggi et al. 2017) and gamma-ray data from Fermi-LAT and VERITAS (Peng et al. 2016;Fleischhack & VERITAS Collaboration 2015).The results are shown in Fig. 4. The underlying CR spectra leads to the CR-induced ionization rate of about  (H 2 ) ≃ 2.5 × 10 −14 s −1 for clouds with  (H 2 ) = 10 23 cm −2 .The ionization rate versus column density is shown in Fig. 3.
We summarize the ionization rate for clouds of column density  (H 2 ) = 10 23 cm −2 for all these prototypical SBNi in Table 2.The ionization rates from CR protons and electrons are also shown separately.It is interesting to note also that the contribution from CR protons is always slightly more dominant than that of electrons, especially at large column densities (see also the lower panel of Fig. 3 for the case of NGC 253).Another point worth mentioning is that our predicted ionization rates for the three SBNi considered only differ slightly.These similarities, in fact, come from similarities in the fitted CR proton spectra.For the discussion of these similarities, let's parametrize for this discussion the CR proton spectrum as follows  p () =  p,300 GeV ()/( = 300 GeV) where () is a function describing the shape of the spectrum and  p,300 GeV is the normalization fixed at  = 300 GeV.We should first notice that the CR transport is mostly dominated by energy loss (calorimetric systems) and, thus, the form of () is determined mostly by the injection spectral index of SNRs  (see Eq. 3) which has rather similar values (between 4.2 and 4.4) as constrained by the gamma-ray spectral index.More importantly, since the gamma-ray fluxes are expected to be dominated by hadronic gamma rays, the normalization of the CR spectra should be, roughly speaking, proportional to the gamma-ray flux, distance squared, and the inverse of the total gas mass of the SBN (as also mentioned in the previous subsection).In fact, it can be shown that  p,300 GeV ∼   (  = 30 GeV) 2 / gas ∼   (  = 30 GeV) 2 /( ISM  3 ).As the values of   (  = 30 GeV) 2 /( ISM  3 ) are only different by a factor of a few for the three SBNi, the differences in the ionization rates are also expected to be of this order.

DISCUSSIONS AND CONCLUSIONS
We have performed fits for the non-thermal emissions to derive the CR spectra in the nuclei of some prototypical starburst galaxies namely NGC 253, M82, and Arp 220.These spectra are then implemented in the ballistic model (see Section 3.2 and also Padovani et al. 2009) to describe the of CRs into dense molecular clouds which allows us to predict the ionization rate for clouds of different column densities.Our predicted ionization rate varies around 10 −14 s −1 for all the prototypical SBNi considered and the values can decrease slightly with increasing column density.
In the case of NGC 253, our predicted ionization rate  (H 2 ) are compared to inferred values from molecular line observations by H22 and B22.In most cases, the inferred ionization rates are a few times to about an order of magnitude higher than our predicted values.Such a discrepancy can be due to (i) Difference in the regions probed by the observations: One important difference between our analysis and the ones from H22 and B22 is the size over which the modelling is performed.The inference of ionization rates by H22 and B22 focused in various molecular clouds of size comparable to the ALMA telescope synthetised beam of 1. s 6, that is roughly 30 pc.Our study, on the other hand, is based on observing constraints coming mostly from gammaray telescopes with large point-spread functions, typically of a few arcminutes (relatively low spatial resolution compared to radio observations of molecular line emissions).This requires us to assume a uniform CR density over the entire SBN in our modeling which should be a good approximation given the high supernova rate within the system.Variations of CR density on small scales, however, might exist on scales comparable to remnant size Phan et al. (2021Phan et al. ( , 2023) ) which should lead to corresponding variations on ionization rates.Modelling such variations might require a better description of not only CR transport but also of the large-scale ISM within these SBNi (relevant for energy loss processes of CRs) and will be examined in our future works.
(ii) Uncertainties in chemical modeling of line observations: The two analyses, H22 and B22, examine the same set of clouds using different line emissions; H22 study H 3 O + and SO and B22 focus on HCN and HNC.The results, however, contain values of ionization rates different by up to an order of magnitude which could mean that the uncertainties in these inferred values are not fully reflected by the errors on ionization rates.In this scenario, it would be interesting to perform a combined analysis taking into account not only line emission data but also non-thermal emission data.Such an analysis might help to reduce the uncertainties on the inferred ionization rates.
(iii) Uncertainties in gas mass and supernova rate: As mentioned in Section 4, the normalization CR proton spectrum as constrained by gamma-ray data should be, roughly speaking, proportional to the gamma-ray flux and the inverse of the SBN gas mass, i.e.
p,300 GeV ∼   (  = 30 GeV)/ gas .In this work, the SBN size and the fitted ISM density correspond to  gas ≃ 7 × 10 7  ⊙ which is actually within the range of the gas mass estimated independently from molecular line observations (between 2.5 × 10 7  ⊙ and 4 × 10 8  ⊙ as presented respectively by Harrison et al. 1999 andHoughton et al. 1997, see also the discussion by Bradford et al. 2003).However, the uncertainty on the gas mass estimate might also mean that its value is actually a few times smaller than expected (see e.g.Mauersberger et al. 1996) which should accordingly require a larger density of CR protons in order to fit the gamma-ray data and, as a result, lead to predicted ionization rates being higher.This can help to improve the agreement between our predictions and the measurements from H22 and B22.Similarly, the uncertainty on the gamma-ray flux, which can be translated into the uncertainty on the fitted value of the supernova rate (see discussions in Section 4), can also be a source of the discrepancy.Indeed, there exists also the possibility that the SBN gas mass is close to the higher end of its uncertainty range and the supernova rate is lower than expected which could further increase the difference between our predictions and measurements.More precise estimates of these quantities are, therefore, essential to improve our understanding of this discrepancy.
(iv) Local sources of MeV CRs: For the SBNi considered, the nonthermal emissions in the GeV and TeV energy, where data are most constraining, are contributed mostly from the decay of  0 created in proton-proton interactions.The production of  0 , however, has a threshold  th ≃ 280 MeV meaning that these gamma-ray data could not probe CR protons with  ≲  th .In other words, there might exist a class of sources accelerating mostly CRs in the energy range of around a few hundred MeVs, e.g.wind termination shocks of stars (Scherer et al. 2008), protostellar jets embedded within molecular clouds (Padovani et al. 2015(Padovani et al. , 2016;;Gaches & Offner 2018), or even HII regions (Padovani et al. 2019;Meng et al. 2019), that contribute to the ionization rate in these systems but could not be observed with GeV and TeV gamma-ray telescopes.We note also that if these MeV sources exist in SBNi and they are sufficiently abundant, they might contribute to the gamma-ray emissions in the MeV energy range, particularly relevant for future missions like eASTROGRAM or AMEGO (de Angelis et al. 2018;McEnery et al. 2019).
Further investigations are required to understand the discrepancies between our predicted ionization rates and the values inferred from molecular observations.The difference by a factor of a few in most cases, however, mean that our predictions for ionization rates in SBNi can be potentially very useful for future chemical and dynamical studies of these rather complex star-forming regions.
Table 1.Parameters for both non-thermal and thermal emissions from the SBNi of NGC 253, M82, and Arp 220.We have fixed the distance  gal , the SBN size , and the galactic wind speed   as in Peretti et al. (2019) as motivated by independent observations.The other parameters are fitted using data from radio to gamma-ray observations (see the comments for more details on parameters and their corresponding most constraining data).Total ionization rate 1.5 × 10 −14 6.9 × 10 −15 2.5 × 10 −14

Figure 1 .
Figure 1.Non-thermal emissions from the SBN of NGC 253 from hard X-ray to TeV gamma-ray domains from  0 decay (dashed red line), inverse Compton scattering (dotted magenta line), bremsstrahlung radiation (dash-dotted green line), and synchrotron radiation (solid blue line).The flux is fitted to gammaray data from Fermi-LAT and H.E.S.S. (H.E. S. S. Collaboration et al. 2018) and upper limits from NuSTAR (Wik et al. 2014).
) where   () is the differential number density for CRs of species  ( = p or  = e respectively for protons and electrons) in the ISM of SBNi as obtained from Eq. 1,  , () ≃     ()/ ISM is the energy loss rate of CRs inside clouds (  is the gas density inside the cloud),  01 =  0 (, ) and  02 =  0 (  − , ) where the function  0 (depending on the species of CRs considered) is the initial energy of CRs as they enter the cloud and is obtained by solving the following equation  = ∫  0    d   () .

Figure 2 .
Figure2.Spectra of CR protons (upper panel) and electrons (lower panel) in the SBN of NGC 253 obtained using the source and transport parameters adopted for the fit of non-thermal emissions.Data of CR spectra observed in the local ISM is also presented for comparison.For CR electrons, we present different components including primary electrons from SNRs (orange line), secondary electrons from decays of  ± (green line), and tertiary electrons from interactions between CR-induced gamma-rays and interstellar radiation (blue line).

Figure 3 .
Figure 3. Upper: Cosmic-ray induced ionization rate expected in the SBN of NGC 253 (solid red line), M82 (dotted blue line), and Arp220 (dashed black line).Inferred values of the ionization rates for a few molecular clouds in the central region of NGC 253 from H22 (green squares) and B22 (orange squares) are also overlaid for comparison.Lower: Cosmic-ray induced ionization rate expected for NGC 253 is shown again together with the 1 and 2 contours of the ionization rate inferred for a particular molecular cloud in the SBN, referred to as GMC 6 in the analyses of both H22 and B22.Seperate contributions of CR protons (dotted red line) and electrons (dashed red line) to the total ionization rate in the SBN of NGC 253 are also presented.

Figure 4 .
Figure 4. Upper panel: Non-thermal emissions from the SBN of M82 from hard X-ray to TeV gamma-ray domains from  0 decay (dashed red line), inverse Compton scattering (dotted magenta line), bremsstrahlung radiation (dash-dotted green line), and synchrotron radiation (solid blue line).The flux is fitted to gamma-ray data from Fermi-LAT and VERITAS (Acero et al. 2015; VERITAS Collaboration et al. 2009) and upper limits from Chandra (Strickland & Heckman 2007).Lower panel: Same as the upper panel but for Arp 220.Gamma-ray data from Fermi-LAT and VERITAS (Peng et al. 2016; Fleischhack & VERITAS Collaboration 2015) and upper limits from Chandra (Paggi et al. 2017) are also overlaid.

Figure A1 .
Figure A1.Fitted spectra in the domain from radio to optical frequency for the three starburst nuclei of NGC 253, M82, and Arp 220 overlaid with data from various observations retrieved from the NASA/IPAC Extragalactic Database.