Solving the multi-messenger puzzle of the AGN-starburst composite galaxy NGC 1068

Multi-wavelength observations indicate that some starburst galaxies show a dominant non-thermal contribution from their central region. These active galactic nuclei (AGN)-starburst composites are of special interest, as both phenomena on their own are potential sources of highly-energetic cosmic rays and associated gamma-ray and neutrino emission. In this work, a homogeneous, steady-state two-zone multi-messenger model of the non-thermal emission from the AGN corona as well as the circumnuclear starburst region is developed and subsequently applied to the case of NGC 1068, which has recently shown some first indications of high-energy neutrino emission. Here, we show that the entire spectrum of multi-messenger data - from radio to gamma-rays including the neutrino constraint - can be described very well if both, starburst and AGN corona, are taken into account. Using only a single emission region is not sufficient.


INTRODUCTION
The high rate of star formation and supernova explosions together with their relatively large abundance make starburst galaxies some of the most promising sources of high-energy cosmic rays (CRs) in the nearby Universe. Multiwavelength observations reveal that in some starburst galaxies the dominant non-thermal emission component originates in their central regions, indicating the presence of an active super-massive black hole. These active galactic nuclei (AGN)-starburst composites are of special interest, as both phenomena on their own are potential sources of highly-energetic CRs which may contribute to the extragalactic CR component observed at Earth, whose origin remains unknown.
The Fermi -Large Area Telescope (LAT) has detected high-energy gamma-rays from a number of nearby starburst galaxies (Ackermann et al. 2012), whereas a handful of those (M82, NGC 253, NGC 1068) have been detected all the way up to TeV gamma-ray energies with imaging atmospheric Cherenkov telescopes (Acero et al. 2009;VERITAS Collaboration et al. 2009;Abdalla et al. 2018;Acciari et al. 2019), confirming that at least a fraction of the starburst galaxy population accelerate particles to very high energies. In addition, starburst galaxies have long been considered prime environments for high-energy neutrino production, e.g. Loeb & Waxman (2006); Becker et al. (2009), due to having higher magnetic fields and gas densities than Milky Way-like galaxies. A first observational hint of possible neutrino production in the starburst-AGN composite NGC 1068 has been seen in 10 years of data from the IceCube Neutrino Observatory: an excess of neutrinos has been found in a region centered 0.35 • away from the coordinates of NGC 1068 after a catalogue-based search using ten years of events from the point source analysis (Aartsen et al. 2020). The excess is inconsistent with background expectations at the 2.9σ level after accounting for trials. NGC 1068 is the brightest and one of the closest Seyfert type 2 galaxies. Due to its proximity, at a distance of 14.4 Mpc (Meyer et al. 2004), it forms a prototype of its class. It has also been predicted as one of the brightest neutrino sources in the northern hemisphere (Murase & Waxman 2016).
In order to understand the emission mechanisms active in these composite sources in general, and in NGC 1068 in particular, and the possible connection to neutrinos, it is necessary to distinguish the non-thermal emission from the different acceleration sites. In particular in the case of NGC 1068, there are strong observational hints for multiple emission sites: firstly, in the near-infrared to radio the ALMA experiment has observed a strong flux cutoff towards smaller frequencies emerging from the inner parsecs of that source (García-Burillo et al. 2016;García-Burillo et al. 2019). Further non-thermal radio emission is associated with the more extended region of a few hundreds of parsecs (Wilson & Ulvestad 1982;Sajina et al. 2011). Secondly, the gamma-ray flux extends out to 100 GeV (Acciari et al. 2019).
Indications that a one-zone model is not enough to explain the multi-messenger emission from NGC1068 have been present for several years. Using radio and gamma-ray observations, it has been shown in previous works, such as Eichmann & Becker Tjus (2016); Yoast-Hull et al. (2014), that the circumnuclear starburst environment alone is not powerful enough (by a factor 10) to account for the observed gamma-ray luminosity. Since NGC 1068 also possesses a large-scale jet as indicated by centimeter radio observations (e.g. Wilson & Ulvestad 1982;Gallimore et al. 2004Gallimore et al. , 2006 it has been suggested (Lenain et al. 2010) that the observed gamma-rays originate in these jets as a result of CR electrons inverse Compton scattering off the infrared (IR) radiation of the surrounding environment. Another possible origin of the observed gamma-ray emission that has been suggested by Lamastra et al. (2019) are CR particles accelerated by the AGN-driven wind that is observed in the circumnuclear molecular disk of NGC 1068. This model predicts a rather hard spectrum extending up to a several TeV which, however, can be excluded due to the upper limits by the MAGIC telescope (Acciari et al. 2019). The potential neutrino emission indicated from the IceCube excess at about 1 TeV poses another problem, as it is at least an order of magnitude stronger than the GeV photon flux. These different components of the high energy phenomena can hardly be explained by a single zone modela high energy neutrino signal with a missing gamma-ray counterpart indicates an optically thick source environment such as the AGN corona (Inoue et al. 2020a;Murase et al. 2020;Kheirandish et al. 2021), whereas the presence of a gamma-ray flux up to some tens of GeV that is only slightly steeper than ∝ E −2 can only be explained by a second source environment. In this work we model the multi-messenger emission of NGC 1068 in the context of a two-zone model which considers the AGN corona and circumnuclear starburst region. In Section 2 we present the details of the two-zone model and the formalism employed to model the multi-messenger emission of the source. In Section 3 we present our fitting approach as well as its results and conclude in Section 4.

THE TWO-ZONE AGN-STARBURST MODEL
The multi-messenger observations of NGC 1068 dictate the need for multiple emission zones. In this work, we will capture its inner µpcs referring to a spherically symmetric structure for the corona of the AGN as well as an outer starburst ring with a radius of ∼ 1 kpc (see Fig. 1). Note, that in between these two emission sites NGC 1068 shows strong indications of a jet structure on scales of up to about 1 kpc (e.g. Wilson & Ulvestad 1982;Gallimore et al. 2004Gallimore et al. , 2006, which however is not included in this work. Due to mathematical convenience we treat both spatial regions as homogeneous. For particle acceleration processes that take place on considerably shorter timescales than the energy loss in these zones, we can disentangle these processes and only describe the steady-state transport of non-thermal, accelerated electrons and protons. Hereby, we suppose that in both zones some acceleration mechanism yields a differential source rate q(T ) of relativistic protons and primary electrons that can be described by a powerlaw distribution in momentum space up to a certain maximal kinetic energyT , which depends on the competing energy loss timescales in these zones. In case of the starburst zone, we suppose that a certain fraction f SN of the total energy of the supernova -that release about 10 51 erg and occur with an approximate rate 1 (Veilleux et al. 2005) ν SN 0.34 [L IR /(10 11 L )]yr −1 dependent on the IR luminosity L IR -gets accelerated into CRs according to diffusive shock acceleration (DSA, e.g. Drury 1983;Protheroe 1999) by individual supernova remnants (e.g. Bell 2014, and references therein). In general, many starburst galaxies-NGC 253 is a prominent example-show a galactic superwind (e.g. Veilleux et al. 2005, and references therin) as a result of the large number of core-collapse supernovae. These winds introduce another source of acceleration 2 (e.g. Anchordoqui et al. 1999;Romero et al. 2018), however, we are not aware of any observational indications of such a superwind in the starburst ring of NGC 1068. For the AGN corona, we suppose that a fraction f inj 1 of the mass accretion rateṀ = L bol /(η rad c 2 ), with a radiation efficiency of η rad = 0.1 (Kato et al. 2008), goes into relativistic protons via stochastic diffuse acceleration (SDA, e.g. Lemoine & Malkov 2020, and references therein). For both zones, the non-thermal primary electrons are normalized by the non-thermal proton rates due to the requested quasi-neutral total charge number of the injection spectra of primary CRs above a characteristic kinetic energy ofŤ 10 keV (Schlickeiser 2002;Eichmann & Becker Tjus 2016;Merten et al. 2017). Note that this corresponds to a quasi-neutral acceleration site, however CR transport can subesequently Figure 1. Sketch of the two-zone AGN-starburst model. The sizes of the regions are not to scale as the (inner) coronal region typically extends out to about some hundreds of µpc, whereas the starburst ring is located at about 1 kpc distance from the nucleus.
remove CR electrons and protons in different amounts from the non-thermal energy regime, nevertheless their charge stays conserved. Transforming the source rates from momentum space into kinetic energy T , we obtain for the injected non-thermal protons, and q e ± (T ) = q e − ,0 for the non-thermal electrons (e − ) and positrons (e + ). Here, the latter term q 2nd e ± (T ) introduces the source rate of secondary electrons and positrons that are generated by hadronic interaction processes, as discussed in the following.
Thus, the steady-state behavior of the differential non-thermal electron and proton density n(T ) in the AGN corona and the starburst zone, respectively, can be approximated by Here, τ cool refers to the total continuous energy-loss timescale, which in case of the relativistic electrons is given by the inverse of the sum of the synchrotron (syn), inverse Compton (IC), non-thermal Bremsstrahlung (brems), and Coulomb (C) loss rates, according to and in case of the relativistic protons we use including the photopion (π), Bethe-Heitler pairs (BH), and hadronic pion (pp) production loss rates. Proton synchrotron losses-as well as the associated radiation-are negligible for the considered environments. Note that these processes require additional information on the associated interaction medium, which is one of the following targets: (i) A magnetic field, which is assumed to be uniform on small scales (with respect to the particles' gyro radius) and randomly orientated on significantly larger scales (due to isotropic Alfvénic turbulence).
(ii) A photon target, which is in the case of the starburst zone dominated by the thermal IR emission due to the re-scattered starlight by dust grains with a temperature θ dust and can be described by an isotropic, diluted modified blackbody radiation field where the dust clouds become optically thick above a critical energy E 0 = 8.2 meV (Yun & Carilli 2002). The constant dilution factor C dil is determined from the observed IR luminosity L IR according to the relation L IR /(πR 2 str c) = dE E n IR (E). 3 In case of the coronal region we used a parametrized model (Ho 2008) above 1 eV that accounts for the optical and UV emission by the disk as well as the Comptonized X-ray emission by hot thermal electrons in the corona. Hereby, the parametrization depends on the Eddington ratio (L bol /L Edd ), i.e. the ratio of the bolometric over the Eddington luminosity, and we adopt the relation of Hopkins et al. (2007) to determine L bol based on the intrinsic X-ray luminosity L X between (2 − 10) keV.
(iii) A thermal gas target with a given temperature θ, which is due to mathematical convenience assumed to be homogeneously distributed in both regions. For the starburst ring Spinoglio et al. (2012) determine θ = 127 K, a gas density of n(H 2 ) = 10 2.9 cm −3 , and a molecular hydrogen mass of M (H 2 ) ∼ 3.5 × 10 8 M More details on the individual energy-loss timescales can be found in the Appendix A.
In addition, we have to account for catastrophic particle losses according to the escape of particles from the considered zones in a total time τ esc . Here, the total escape rate via gyro-resonant scattering through turbulence with a power spectrum ∝ k −κ and a turbulence strength η −1 as well as a bulk stream flow can be approximated by (Murase et al. 2020) with respect to the characteristic size R and the magnetic field strength B of the zones. Here, e denotes the elementary charge and c refers to the speed of light. In the case of the starburst zone, the bulk motion is given by a galactic wind with a velocity v w ; in the case of the AGN corona, we account for the infall timescale, which is expected to be similar to the advection dominated accretion flow (Murase et al. 2020), with a viscosity parameter a vis and the Keplerian velocity v K = G M BH /R. For the spectral index of the turbulence spectrum it is common to assume a value of either κ = 5/3 referring to Kolmogorov turbulence, or κ = 3/2 (Kraichnan turbulence) which can be motivated from isotropic MHD turbulence in the magnetically dominated regime. In the following, we will adopt Kolmogorov turbulence for both regions unless stated otherwise. As discussed by several previous works, e.g. Kheirandish et al. (2021); Inoue et al. (2019), the stochastic diffuse acceleration (SDA) in the coronal region typically appears to be inefficient compared to the cooling rates. Therefore, it has been suggested that there needs to be some other acceleration mechanism such as magnetic reconnection at work. However, there is very little known about the actual acceleration efficiency of this process in the AGN corona as well as the resulting spectral shape of the CR energy distribution, so that we choose to stick to the SDA process at first. Hereby, we adopt the ansatz that the same scattering process that yields the diffusive escape is also responsible for the stochastic acceleration. In the case of the starburst region, the acceleration is expected to be introduced by a multitude of supernova remnants (SNRs) via diffusive shock acceleration (DSA).
Here, the CRs are kept within the accelerating shock region by the Bohm diffusion-which gets introduced by the CR self-generated Bell instability (Bell 2004)-where the wave turbulence typically inherits the flat energy spectrum of the generating cosmic rays, i.e. κ = 1. Further, we assume that the average shock speed v sh equals the wind speed v w . Thus, we use the acceleration timescale (e.g. Murase et al. 2020;Romero et al. 2018) for the starburst,    Table 1. The thin solid black line refers to the total energy loss timescale (τ cool ) and the thin dashed black line refers to the total escape timescale (τesc). The abbreviated individual process are as follows: pp=hadronic pion production, pγ=photopion production, BH=Bethe Heitler pair production, syn=synchrotron radiation, diff=diffusion, adv=advection, acc=acceleration (DSA for the starburst and SDA for the corona), Coul.=Coulomb losses, Brems.=Bremsstrahlung, IC=inverse Compton scattering. as the gas density or the magnetic field strength-the different timescales differ significantly between the starburst and the corona zone, as shown in Fig. 2. With respect to the maximal CR energyT it is shown that in case of the starburst DSA can typically provide a maximal primary proton (electron) energy of several tens of TeV (hundreds of GeV). In the AGN corona a high turbulence strength (η ∼ 1) or a rather flat turbulence spectrum (κ 3/2) is needed to obtain CR proton energies of about 100 TeV or more which would be necessary to stay within the limits of the IceCube observations at the indicated potential flux (Aartsen et al. 2020). Still, the coronal CR electrons typically suffer from significant synchrotron/ IC losses at energies as low as a few tens of MeV. In addition, these primary electrons also have to overcome the Coulomb losses at low energies, and hence, need to be injected into the SDA process at about few keV. 4

Spectral energy distribution of CR electrons and protons
Solving the transport equation (3) by fundamental methods (see e.g. Eichmann & Becker Tjus 2016) provides the differential non-thermal electron and proton density n(T ) in the starburst as well as the coronal zone, as shown in Fig. 3. Despite the assumed quasi-neutrality of the primary source rates, the resulting fraction of CR protons is significantly higher than the one of the CR electrons, especially in case of the coronal region. This is due to significantly smaller energy loss time scales of the relativistic electrons according to synchrotron and IC losses (see Fig. 2 ). Both regions are perfect calorimeters (except for CR protons of the starburst region with energies above about 1 PeV), so that the energy distributions steepen at certain characteristic energies either due to one of the energy loss processes or due to the exponential cut-off introduced by the acceleration. In Fig. 3 we consider a parameter scenario that results from a fit to the data (as introduced in Sect. 3) and it can be seen that, for this particular scenario, the secondary electrons have a major contribution to the total CR electron spectrum. The individual contributions to its source rate q 2nd e will be introduced in the following.

Non-thermal emission
Using the spectral energy distribution of CR protons and electrons we determine the non-thermal emission of those particles from radio to gamma-ray energies. Hereby, both emission regions are considered to be spherically symmetric, with the inner coronal region constrained to a radius R cor = r cor R s , where typically r cor ∈ [1, 100] and R s = 2GM/c 2 denotes the Schwarzschild radius, and the surrounding starburst ring extends between an inner radius R in str ∼ 1 kpc (Rico-Villas et al. 2021) and an outer radius R out str ∈ [1.3, 2] kpc. In the following we include the non-thermal photon emission by synchrotron radiation, inverse Compton scattering, non-thermal bremsstrahlung from the CR electrons according to a spectral emissivity syn , ic , and brems , respectively, as well as hadronic and photo pion production by the CR protons with a spectral emissivity pp and pγ , respectively. These hadronic processes also introduce high-energy neutrinos as well as secondary electrons and positrons with a source rate q 2nd e = (e) pp + (e) pγ + BH + γγ , where the last two terms introduce electrons/positrons by γγ and Bethe-Heitler pair production.
In addition to these non-thermal emission processes, we also account for the free-free emission ( ff ) by the considered thermal target gas. Note that we use the spatially averaged thermal gas distribution in both zones. The actual gas distribution might differ. This would change the corresponding free-free emission. Further there are hints (Inoue et al. 2020b) of a colder, less dense gas distribution in the extended coronal region which might contribute at radio/IR frequencies. But since there is no need for such an additional component in our model, we neglect the thermal emission by the extended coronal region to limit the number of free parameters. The total spectral emission rate of the considered messenger particle m = (γ, ν) is given by Moreover, it is necessary-in particular for the coronal region-to account for different absorption processes such as synchrotron-self absorption (α syn ) and free-free absorption (α ff ) at radio/IR energies as well as γγ pair production (α γγ ) at gamma-ray energies, yielding a total absorption coefficient Note that we consider both source regions as optically thin for the neutrinos, which does not necessarily need to be the case for the dense coronal region at the highest energies. More details on the individual emissivities as well as the absorption coefficients can be found in the Appendix B and C. For a homogeneous source for which m and α m are constants, the spectral energy flux at a source distance d is generally given by (e.g. Gould 1979) with the effective emission volume where r s represents an arbitrary position on the surface of the emission volume. Due to the different geometry of the considered emission sites also the effective emission volume differs. Considering the directional symmetries of both systems the previous equation can be simplified to with and a ± = (R out str ± R in str )/2. Note that for optically thin emission, i.e. in the case that α m is significantly smaller than the characteristic length scales of the system, we obtain as expected.
3. EXPLAINING THE MULTI-MESSENGER DATA

Data and Fit Procedure
In the following we use the previously introduced model to explain the multi-messenger data of NGC 1068. Hereby, we quantify the goodness of the fit by the chi-squared value ) denotes the model prediction (observation) of the flux of photons and neutrinos at different energies. As the model prediction depends in general on a large parameter set (see Table 1), we use the differential evolution algorithm 5 to find a global minimum within a dedicated subset of this parameter set. Hereby, we keep fixed those parameters that have either a minor impact on the resulting flux prediction-such as the wind speed v w and the characteristic temperature θ dust of the dust grains in the starburst region, the radiative efficiency of the disk, and the viscous parameter µ vis in the accretion flow (where we adopt the same value as suggested by Schartmann et al. 2010)-or are rather well defined from observations: we use the intrinsic 2−10 keV X-ray luminosity L X = 10 43 L 43 erg/s of the coronal region as derived from two different analyses, as this luminosity carries large uncertainties due to the high column density of the source. Hence, a detailed modelling of the impact of the torus is necessary to obtain the actual fraction of X-rays that is scattered into the line of sight. In the energy range between 2 and 10 keV NuSTAR and XMM-Newton monitoring campaigns (Marinucci et al. 2016) yield an intrinsic luminosity of L 43 = 7 +7 −4 , whereas a different analysis by Ricci et al. (2017) that combines 70-month averaged Swift/BAT data with different measurements in this soft X-ray band obtains L 43 = 0.9.
In general, it is expected (e.g. Mayers et al. 2018) that a higher X-ray luminosity corresponds to higher black hole mass M BH than a low X-ray luminosity. As shown by different observations (see GRAVITY Collaboration et al. 2020, and references therein) the expected range of the black hole mass yields M BH ∈ [0.8, 1.7] × 10 7 M and since a higher X-ray luminosity generally indicates a higher value of M BH (e.g. Mayers et al. 2018), we use M BH (L 43 = 0.9) = 0.8 × 10 7 M and M BH (L 43 = 7) = 1.7 × 10 7 M , respectively, yielding an Eddington luminosity of L Edd (L 43 = 0.9) = 1.0 × 10 45 erg s −1 and L Edd (L 43 = 7) = 2.1 × 10 45 erg s −1 , respectively. And based on the relation by Hopkins et al. (2007) we obtain a bolometric luminosity of L bol (L 43 = 0.9) = 0.2 × 10 45 erg s −1 and L bol (L 43 = 7) = 2.9 × 10 45 erg s −1 , which is nicely within the range that has been found by others (see GRAVITY Collaboration et al. 2020, and references therein). Hence, this also has a direct consequence on the adopted black hole properties, so that the Schwarzschild radius R s as well as the mass accretion rateṀ decrease for a decreasing L X value. The surrounding starburst region is dominant in the IR, due to scattering off of dust grains with a characteristic temperature θ dust . In the following we adopt θ dust = 127 K as observed by Spinoglio et al. (2012), although also higher temperatures of up to 250 K have been observed recently (Rico-Villas et al. 2021). But as previously mentioned the impact of θ dust on our results is negligible. Further, we use the observed bolometric mid-infrared luminosity L IR of the outer starburst ring of about 30 in diameter, which is found to account for almost half its total mid-infrared luminosity (Bock et al. 2000).
But also for the other, non-fixed parameters we account for physical constraints based on observations or numerical simulations. For example the optical depth of the coronal region is typically ω T ∼ 0.1 − 1 (Ricci et al. 2018;Merloni & Fabian 2001) constraining the gas density according to n gas = ω T /(σ T R) if the coronal plasma is not dominated by electron-positron pairs. And based on the virial gas temperature θ gas = GM BH m p /(3R cor k B ) = m p c 2 /(6r cor k B ) of the protons we are able to draw constraints on the coronal magnetic field strength according to B = 8πn gas k B θ gas /β, where a low plasma beta (β ∼ 0.1 − 3) is expected from numerical MHD simulations (Jiang et al. 2019(Jiang et al. , 2014Miller & Stone 2000). All of the subsequently used parameter values and the constraints, respectively, are summarized in Table 1. The observational data that were used in the fitting procedure are as follows: Table 1. Summary of fixed and free parameters that describe the radiation zones of NGC 1068 in our model. The top lines give the parameters that describe the starburst region, namely the spectral index of the electron and proton distributions, s, the fraction of supernova energy that gets converted to CRs, fSN, the gas density, ngas, the magnetic field strength, B, the radius of the starburst ring, Rstr, the turbulence strength parameter, η, the spectral index of the turbulence spectrum, κ, the minimal kinetic energy of the CR particles, Tinj, the average wind speed, vw, and the infrared luminosity, LIR. The bottom two lines give the parameters that describe the AGN environment which include some of the starburst related parameters and in addition, the fraction of X-ray emissivity that goes into relativistic protons, finj, the radius of the corona, Rcor, the viscous parameter, µvis, and the intrinsic coronal X-ray luminosity, LX. For the free parameters, the scanning range is reported in the (i) The radio data of the starburst region have first been observed by Wilson & Ulvestad (1982), who could, however, only disentangle the large-scale emission at 20 cm. Therefore, Wynn-Williams et al. (1985) re-observed the large scale emission of that source with VLA at 2, 6 and 20 cm using 6 .5 FWHM beams. Based on the given range of the flux density of (4 − 8) mJy per beam at 2 cm as well as the spectral behavior, we obtain an integral flux of 0.10 ± 0.03, 0.31 ± 0.10 and 0.81 ± 0.27 Jy at 2, 6 and 20 cm, respectively. Here, we consider the spatial region at a distance between 10 and 25 with respect to the nucleus due to the observed plateau in the flux density.
Using the latest distance measurements-so that 1 ∼ 70 pc-this large scale emission arises from a size between about 0.7 and 1.8 kpc. The resulting flux at 20 cm is only smaller by a few percent than what has been found earlier by Wilson & Ulvestad (1982), although these authors considered distances up to 2 . Further, both works conclude that the origin of this emission is almost certainly synchrotron emission. We are not aware of any additional recent data of this spatial region in the radio band, which could be used to constrain our starburst model at low energies.
(ii) The radio and IR data of the coronal region, require a high spatial resolution of the observational instrument to exclude any additional contribution, e.g. by the torus which emits predominantly in the IR. Two components have been identified within the central parsecs of that source. One of those is the compact (0.5 − 1.4) pc sized core of the AGN which we will associate with the coronal region in the following. Still these spatial scales are at least three order of magnitude larger than the inner corona, that we are modeling here. Thus, there is the chance that the so-called extended corona region at distance of about mpc to pc from the black hole, provides an additional contribution. Previous works (Roy et al. 1998;Gallimore et al. 2004;Inoue et al. 2020b) have shown that free-free emission by a gas with an electron temperature of a few ×10 6 K of this extended corona is actually able to explain the VLBA radio observation. Despite the recent objections (which we will discuss in more detail in Sect. 4) by Baskin & Laor (2021), we adopt this simple approach using a gas with an (electron/proton) density of 2.5 × 10 5 cm −3 and an electron temperature of 10 6 K to explain the resolution matched flux densities (Gallimore et al. 2004) of < 0.7, 5.9 ± 0.5 and 5.4 ± 0.5 mJy at 1.4, 5 and 8.4 GHz, respectively (see Fig. 5). Note, that the assumed gas becomes optically thick at about 5 GHz which introduces a slight tension with the observed flux at that frequency, but yields a flux at 1.4 GHz that is in agreement with the upper limit. To explain the steep flux increase in the IR, however, an additional contribution is needed which we suggest to be given by the inner corona. At these frequencies the ALMA observatory has determined a flux density of 6.6 ± 0.3 and 13.8 ± 1.0 mJy at 256 and 694 GHz, respectively, using a beam sizes of 20 and 60 mas, respectively (Impellizzeri et al. 2019;García-Burillo et al. 2016). In addition, we also include the recent results of the data analysis of the continuum fluxes at 224, 345, and 356 GHz with a beam sizes of 30 mas (Inoue et al. 2020b). The different resolutions of these observations introduce a mismatch in particular if an additional, strong flux contribution by the torus emerges at about a few tens of pc. Hence, the flux prediction at 224 GHz is expected to become smaller then the one from Impellizzeri et al. (2019) at 256 GHz, if the mismatched coverage would be taken into account properly. So, we account heuristically for this effect by introducing an additional, lower uncertainty ∆F add = (1 − 20 mas/δ bs ) F obs of the observed flux F obs dependent on the beam size δ bs . 6 (iii) The gamma-ray data are taken up to 100 GeV from the fourth Fermi-LAT catalog of gamma-ray sources (Fermi-LAT collaboration et al. 2022), and at higher energies we include the upper limits from the MAGIC telescope (Acciari et al. 2019). At these frequencies one cannot resolve individual spatial regions and the data need to be explained by the total flux of both regions.
(iv) The high energy neutrino flux that corresponds to the 2.9 σ excess observed by the IceCube Neutrino Observatory is taken from Figure 7 of Aartsen et al. (2020). In terms of the chi-squared calculation we only account for the most well constrained flux value at about 1 TeV as well as flux at 28 TeV to account for the steep spectral behavior.
Note that we do not account for the flux attenuation from the coronal region by the torus, which is mostly relevant for 1 MeV. Here detailed modelling by Ricci et al. (2017) of the torus absorption 7 with respect to the broadband X-ray characteristics has shown that in the soft (hard) X-ray at 2 − 10 keV (14 − 195 keV) the coronal flux gets attenuated by a factor 0.018 (0.17). Since the corona is a perfect CR calorimeter, as indicated by Fig. 2, there is no additional (hadronic or leptonic) emission by CR interactions in the torus region, so that we can completely neglect this region in the following. In addition to the inner corona, we also account for the free-free emission by the outer corona, that extends up to about 1 pc. As previously mentioned in (ii), its parameters are not changed by the fit algorithm, but fixed to explain the compact radio data.

Fit Results
Using the (2×6)-dimensional parameter space of constrained fit parameters-as introduced by the first six columns in Table 1-we obtain a robust global chi-squared minimum of min(χ 2 ) 8 (10) for an intrinsic coronal X-ray luminosity of L X = 7 (0.9) × 10 43 erg/s. Even though the resulting best-fit spectra are almost equal for both of those two casesexcept for the resulting neutrino flux at 1 TeV-the resulting best-fit parameter space of the corona shows some differences: For L X = 0.9 × 10 43 erg/s the inner corona needs to extend about (150 − 200) R s , which however, is about the same absolute size as for the case of a high X-ray luminosity. To further obtain a sufficient amount of CRs, a higher value of f inj is needed due to the comparably small mass accretion rate. In addition the initial CR spectrum in the corona needs to be slightly softer (s ∼ 2) to explain the IR and neutrino data, due to the smaller loss rate by IC scattering as well as photopion and Bethe-Heitler pair production. But the rest of the resulting parameter space is similar to what is described in the following, in particular for the starburst ring.
In the case of L X = 7 × 10 43 erg/s, the Fig. 4 shows the goodness of the fit for a certain range of the parameter space, where the chosen evolution strategy has converged. Here, the so-called 'best1bin' strategy is used, where two members of the population are randomly chosen and the difference is used to mutate the best member. Hence, the algorithm tends to increase the number of chi-squared function evaluation if the algorithm converges towards its minimum. An extended minimum in the parameter space is, in general, still favored with respect to a narrow one, so that it cannot be excluded, especially in such a multi-dimensional parameter space, that the resulting minimum is actually not a global but a local one. Since these are two-dimensional representations of the (2 × 6)-dimensional parameter space, we have to marginalise over the other dimensions, which is done by using its minimal chi-squared value χ 2 min . A systematic scan of the whole parameter space would be needed to expose all of the details of the χ 2 -distribution. Still we can conclude that almost all of the 12 fit parameters have a significant impact on the goodness of the fit and the best-fit parameters (with χ 2 10) can be well constrained. However, this does not apply to the outer radius R out str of the starburst ring, to which the fit results are not sensitive. Further, this best-fit parameter range does not represent any extreme scenarios, however, the inner corona needs a rather high gas density and magnetic field strength which extends up to about (80 − 100) R s -note that the observed black hole mass of NGC 1068 is rather small compared to other AGN yielding a small Schwarzschild radius, so that 100 R s 1.5 × 10 −4 pc. Still, in combination with a high gas density of ∼ 10 10 cm −3 this yields a rather high value of the optical Thomson depth of ω T 1. In addition, these fit results suggest a strongly magnetized inner corona with a plasma beta β 1. In case of the starburst region a high gas density of n gas ∼ 800 cm −3 , as suggested by Spinoglio et al. (2012), yields good agreement with the data, even though the best-fit result is obtained for a gas density that is about a factor of three smaller.
In the high (low) coronal X-ray luminosity case, the best-fit parameters of the coronal gas density, radius and magnetic field strength correspond to an optical Thomson depth of ω T 2.6 (4.4) as well as a coronal plasma beta of β 0.1 (in both cases). Here, a CR luminosity of L CR = 9.55 (2.07) × 10 43 erg s −1 is needed which is about 4.5 % (2.1 %) of the Eddington luminosity. Further, we obtain a coronal CR-to-thermal gas density pressure of P CR /P gas = 0.30 (0.10). We checked that also for the case of Kraichnan turbulence (κ = 3/2) quite similar fit results can be obtained, however, some of the parameter values change-such as e.g. slightly larger plasma beta values (of up to 0.2) and smaller CRto-thermal gas density pressure ratios (of 12 and 6% in case of a high and low X-ray luminosity, respectively) for the corona. Fig. 5 shows that the minimal χ 2 fit to the data yields an almost perfect agreement with the data-except for the 4FGL data point at the highest energy. The large scale radio data are described by synchrotron radiation of predominantly secondary electrons from the starburst region, whereas the small scale radio/IR data result from the coronal synchrotron radiation of mostly primary electrons. Due to the optical thickness by synchrotron self-absorption, this spectrum cuts off sharply towards small frequencies in the FIR. At the high energy end of the spectrum, where the impact of the torus attenuation vanishes, the γ-ray contribution of the corona and the starburst are at about the same level at about 50 MeV, so that at low γ-ray energies both regions are needed to describe the multi-wavelength data. At about few ×100 MeV the coronal γ-ray emission becomes subdominant and the γ-rays of the starburst ring are sufficient to explain the observed γ-ray data above about 1 GeV. Its high γ-ray luminosity is mostly a consequence of the estimated SN rate of 0.5 yr −1 . However, we verified that even for a significantly lower rate the data can still be explained-although the chi-squared value increases due to increasing deviations in the radio and γ-ray band. In this case a higher coronal γ-ray flux is needed that compensates the lack of γ-rays from the starburst ring and explains the data up to about 1 GeV. The observed IceCube neutrinos can be explained by the corona, as the neutrinos-in contrast to the associated γ-rays-are able to leave this central region. However, in the the best-fit case of a low coronal X-ray luminosity (L X = 0.9 × 10 43 erg/s) the resulting neutrino flux at 1 TeV is about a factor of five smaller  than what is shown in Fig. 5 (but still matches the potential IceCube flux at 28 TeV). Independent of the adopted coronal X-ray field or the particular fit scenario, our model suggests a hardening of the neutrino flux below about 1 TeV. An alternative fit scenario with χ 2 ∼ 14, that enables a higher neutrino flux at 1 TeV for the low X-ray luminosity case is briefly summarized in the following: Using a higher injection fraction (f inj ∼ 0.06) and a smaller radius (∼ 95 R s ) for the corona, it becomes possible to match the potential neutrino flux also at around 1 TeV. However, in that case the corona yields a higher γ-ray flux at a few ×100 MeV, so that the starburst ring needs to be negligible at these energies. Hence a somewhat smaller gas density (n gas ∼ 100 cm −3 ) and harder initial CR spectrum (s ∼ 2.1) in the starburst ring is needed, but still the data at about 100 MeV is slightly overshot due to the additional minor contribution by the starburst ring. At low energies the data is still explained quite accurately, in which the coronal IR emission results from synchrotron radiation of secondary electrons, whereas primary synchrotron emission is no longer present. However, this scenario yields a much higher CR pressure (P CR /P gas 0.50) in the corona, so that altogether we consider this alternative scenario to be less likely than the best-fit scenarios that have been described previously.   Figure 5. The model predictions of the photon and neutrino SED of NGC 1068 with respect to the data-red markers refer to a beam size of ∼ (0.02 − 0.06) arcsec, and black or blue markers indicate a beam size of 10 arcsec. The light red area shows the free-free emission from the extended corona. The dark grey area indicates the internal flux of the background (target) photon fields (disk-and torus emission as well as Comptonized X-rays of the AGN corona) of the central AGN and the light grey area indicates the thermal IR emission by dust grains of the starburst region. Note that the torus attenuation is not taken into account here.

CONCLUSIONS AND DISCUSSION
In this work, we introduced a spatially homogeneous, spherically symmetric, steady state two-zone model for AGNstarburst composite galaxies. Using the multi-messenger data of NGC 1068 from the radio up the the γ-ray band as well as its recent indications of high-energy neutrino emission, we present a first application of this model. Hereby, we perform a global parameter optimization within the (2 × 6)-dimensional parameter space and manage to perfectly explain all data-except for some minor deviations of the γ-ray flux at about 10 GeV. So, the γ-ray emission above a few ×100 MeV results predominantly from the starburst region, whereas the high-energy neutrinos at TeV energies must originate from the coronal region. As already discussed in other works (e.g. Murase et al. 2020;Inoue et al. 2020b;Kheirandish et al. 2021) the corona is optically thick for the associated γ-rays, which introduced a cascade of secondary electrons that dominate the emission at 0.1 eV E γ 100 MeV via synchrotron and IC radiation-in addition to the strong free-free emission of the hot gas. In contrast to these other works we however manage to explain the high-energy neutrino emission by using an acceleration scenario where the CRs are scattered off stochastically by Alfvénic turbulence that shows either a small spectral index (κ 3/2) or a turbulence strength parameter η ∼ 1. Hence, there is no need for an alternative acceleration process such as magnetic reconnection, as studied extensively in Kheirandish et al. (2021), even though, such an alternative acceleration scenario in the AGN corona would relax the need for strong Alfvénic turbulence. Further, the resulting gas density, radius and magnetic field strength of the corona yield a rather high optical Thomson depth of ω T 1 as well as a low plasma beta of β ∼ 0.1, which however, is within the range of expectations (e.g. Ricci et al. 2018;Miller & Stone 2000). Some of our best-fit scenarios suggest a rather large CR pressure of about 30% of the thermal gas pressure, hence, a huge amount of the gravitational binding energy goes into CRs. But this becomes less extreme if we account for the additional energy that is supplied by the disk.
Finally, we manage to explain all data well in case of a strongly magnetized corona and a starburst ring with a high supernova rate (∼ 0.5 yr −1 ). Such a full multi-messenger fit from radio to TeV energies in photons plus the potential neutrino flux has not been attempted before. In particular, using the pure AGN core model has difficulties explaining the full high-energy signatures (Murase et al. 2020;Inoue et al. 2020b;Kheirandish et al. 2021). Including the additional contribution from the starburst ring obviously helps to explain the photon emission above about 100 MeV, but also with respect to the coronal high-energy neutrino emission the detailed fitting approach enables us to find a much better agreement to the potential neutrino flux.
In total we showed that the broadband multi-messenger data of NGC 1068 can only be explained if we account for the non-thermal emission by the outer starburst ring as well as the inner corona. However, we are not able to explain the VLBA radio data of the central region by the inner corona region, neither via free-free emission (due to the high electron temperature), nor via synchrotron radiation (due to the optically thickness at these frequencies for a magnetic field strength of > 10 G). Therefore, we followed the common assumption (see e.g. Roy et al. 1998;Gallimore et al. 2004;Inoue et al. 2020b) and introduced the extended coronal region (extending up to 0.7 pc) to explain these data.
Here we suppose that this extended corona is filled with a thermal gas with an electron temperature of 10 6 K that emits free-free radiation and becomes optically thick at about 5 GHz, which yields an appropriate agreement with the VLBA data. Based on the effect of radiation pressure compression on an ionized gas Baskin & Laor (2021) recently showed that the brightness temperature of a dusty gas is limited to 2 × 10 5 K at 5 GHz. In addition, they showed that a hot free-free emitting gas (with electron temperatures of 10 7 K) also over produces the observed X-ray luminosity of NGC 1068. Therefore, they exclude optically thin free-free emission in NGC 1068, on the sub pc scales. We noticed, that due to the given maximal extension of this region an electron temperature of at least 10 6 K is needed and the upper flux limit at 1.4 GHz can only be satisfied if the free-free emission already becomes optically thick at about 5 GHz. In that case the free-free emitting gas is subdominant at X-rays and necessary electron temperature could be realized for a dustless gas even under consideration of radiation pressure compression.
As previous models, our model is limited with respect to the spatial description of the two emission zones, so that inhomogeneities and magnetic field structures cannot be taken into account. Hence, a more accurate treatment of the spatial structures as well as the three-dimensional transport might change some of the details of these results and should be taken into account in future work.
In general, more data in particular in the range of about (1 − 100) MeV would be very useful to further constrain the model. At lower energies the coronal emission is expected to be attenuated by the torus, so that it would become necessary to account for the physical processes in the torus region, if the data cannot resolve the sub-torus structures. As the coronal region is a perfect CR calorimeter, we do not expect any additional non-thermal emission from that region. Hence, it is not expected that the model prediction benefits from the inclusion of the torus, as this involves another significant expansion of the parameter space. But in case CR protons get also accelerated up to TeV energies in the torus-such as by winds from the coronal region that impact the torus and trigger shocks as proposed recently by Inoue et al. (2022)-the observed GeV photons could also originate from hadronic pion production with the torus gas.
Another sub-structure of NGC 1068-which we do not take into account-is its jet that has been observed by centimeter radio observations on scales of a few ×100 pc (e.g. Wilson & Ulvestad 1982;Gallimore et al. 2004Gallimore et al. , 2006. Therefore, we cannot exclude that there is some minor contamination of the considered large scale radio data by the jet. In addition, this jet has previously been discussed (Lenain et al. 2010;Lamastra et al. 2019) as a possible origin of the γ-ray signal, as the necessary γ-ray luminosity of NGC 1068 is typically too high to be explained only by the nuclear starburst activity (see also Eichmann & Becker Tjus 2016;Yoast-Hull et al. 2014). Based on the bolometric IR luminosity of the starburst ring we estimate an average SN rate of 0.5 yr −1 , which is almost an order of magnitude higher than what has been supposed for the nuclear starburst activity and right within the range of (0.1 − 1) yr −1 that has been suggested by Mannucci et al. (2003); Wilson et al. (1991). Thus, if about (4 − 5)% of that SN energy is converted into CRs-which is similar to what has been found in numerical simulations (e.g. Haggerty & Caprioli 2020)-there is no need for the non-thermal jet emission to explain the γ-ray data. But since the jet as well as the torus region is not taken into account in this work, we cannot exclude that they provide some contribution to the observed γ-ray flux.
We would like to thank the anonymous referee for constructive comments that helped to improve the original version of this paper. Thanks also go to M. Kachelrieß, M. Zacharias, as well as E. Kun for fruitful discussions in terms of different aspects of this work. Further, we acknowledge funding from the German Science Foundation DFG, within the Collaborative Research Center SFB 1491 "Cosmic Interacting Matters -From Source to Signal". In addition, BE acknowledges support by the DFG grant EI 963/2-1. Software: Some of the results in this paper have been derived using the software packages Numpy (Harris et al. 2020), Scipy (Virtanen et al. 2020), Pandas (McKinney 2010, Matplotlib (Hunter 2007), Seaborn (Waskom 2021). APPENDIX A. ENERGY LOSSES CR electrons lose energy due to bremsstrahlung, ionization, synchrotron radiation and inverse Compton scattering, while for CR protons the processes to consider are synchrotron radiation, hadronic pion production, Bethe Heitler pair production and photopion processes. Before illustrating these energy loss rates, the Lorentz factor γ of CR electrons and protons, respectively, and the dimensionless velocity β are introduced: Here T stands for the relativistic kinetic energy of a CR particle with a rest energy E 0 . In the following, we add an index to those quantities to specify that they refer to the CR electron (e) or protons (p).
For a fully ionized medium with density n Z , the relativistic electron bremsstrahlung e + p −→ e + p + γ energy loss rate is given by the expression (Dermer & Menon 2009) where α F is the fine structure constant, σ T is the Thomson cross section and Z is the atomic charge number of the particles characterising the medium.
Another way high energy particles lose energy is by ionization and excitation to bound atomic levels of the matter they travel through. In a fully ionized plasma, the energy loss is due to scattering of individual plasma electrons and due to the excitations of large-scale systematic or collective motions of many plasma electrons. The energy loss rate due to this process is (Schlickeiser 2002 where n e denotes the thermal electron density, which is about equal to the gas density in a fully ionized plasma. The low energy part of the NGC 1068 emission spectrum is strongly affected by synchrotron radiation, created when charged particles move through a magnetic field. The electron synchrotron cooling rate is given by (Blumenthal & Gould 1970) The inverse Compton process involves the scattering of low energy photons to high energies by ultra-relativistic electrons so that the photons gain and the electrons lose energy. Integrating the inverse Compton power of photons over all scattered photon energies E γ , we obtain the energy loss of a single relativistic electron due to inverse Compton scattering (Schlickeiser 2002) where (Blumenthal & Gould 1970) F KN (q, Γ) = 2q ln(q) + 1 + q − 2q 2 + (Γq) 2 (1 − q) 2 (1 + Γq) (A6) with Γ = 4E γ e /(m e c 2 ) and q = E/[Γ γ e m e c 2 − E ]. The terms E and E are, respectively, the energies of the photon before and after the inverse Compton scattering.
When considering protons, the energy loss rate due to synchrotron radiation needs to be rescaled due to the decrease of the cross-section, so that Also the energy loss rate of a fast proton due to Coulomb interactions with the fully ionized plasma has a different cross-section yielding (Schlickeiser 2002) (τ where x m = 0.0286 (θ e /2 × 10 6 K) 1/2 . In addition, relativistic protons can interact with the gas protons and produce pions, in the so called hadronic pion production process: p + p −→ π + X (X is anything else created in the pp collision (e.g. Becker 2008; Dermer & Menon 2009)). The particles energy loss rate can be approximated in the range 1.2 GeV < E ≤ 10 8 GeV by (Krakau & Schlickeiser 2015) τ where n gas is the interstellar gas density and H() denotes the Heaviside function to account for these losses only above 1.2 GeV. Note that we included the β p dependence to enable an extrapolation towards mildly relativistic energies.
Proton interactions with the background (target) photons can produce e ± pairs (p + γ −→ p + e + + e − ), which can lead to electromagnetic cascades. This process is called Bethe-Heitler pair production and it is described by the characteristic particles energy loss rate of (Zheng et al. 2016) where E = γ p E (1 − β p cos λ) is the energy of the photon in the rest frame of the proton with the angle between the proton and photon directions λ and the proton velocity β p is in units of c. The terms E min and E max correspond respectively to 1 MeV/m e c 2 and 2γ p E, while E γ,min = 1 MeV /(2γ p ).
When the interaction between a proton and a target photon produces a pion (p + γ −→ π + p), the typical energy loss rate of the initial particle is given by (Dermer & Menon 2009) where γ p is the proton Lorentz factor and K pγ (E ) is the fraction of energy lost by the ultrarelativistic proton (γ p 1 and β p −→ 1) in the interaction, therefore the inelasticity of the collision.

B. EMISSIVITIES
The emissivities (in units of cm −3 s −1 eV −1 ) of secondary particles (photons, electrons/positrons and neutrinos) with an energy E that we introduce in the following account for (i) synchrotron radiation, inverse Compton scattering and bremsstrahlung of CR electrons; (ii) hadronic and photo-hadronic pion production as well as Bethe-Heitler pair production of CR protons; (iii) free-free emission of the thermal gas; and (iv) γγ pair production. Table 2. Multiplicities ζ and mean fractional energies χ of secondaries formed in photomeson production.
Charged pions formed by the photomeson process decay into leptons and neutrinos, and neutral pions decay into γ-rays. The emissivity of these secondaries, assuming isotropy of the photon and ultra high energetic cosmic ray proton spectra, is described by the following expression (Dermer & Menon 2009) Hereby, for secondaries generated by single pion production we use σ s πγ = 340 µb,Ẽ l = 390,Ẽ u = 980, whereas for those produced in multi pion production we adopt σ m πγ = 120 µb,Ẽ l = 980,Ẽ u −→ ∞. Further, ζ i denotes the multiplicity of secondary i, χ i is the mean fractional energy of the produced secondary compared to the incident primary proton (see Table 2), n p (γ p ) is the differential intensity of CR protons with a Lorentz factorγ p ≡ E/(χ i m p c 2 ). The dimensionless target photon energy E γ is described in units of E e,0 .
The free-free radiation emission of a thermal gas with about the same number density of electrons and ions is defined as (Padmanabhan 2000) where θ is the temperature and g f f is the so called Gaunt factor, which is given by (Padmanabhan 2000): π ln 2k B θ (1.78me) 3/2 2me 1.78q 2 e πE/h , for θ < 8.9 × 10 5 K , √ 3 π ln 4k B θ 1.78E , for θ ≥ 8.9 × 10 5 K .
For the hadronic pion production process the emissivity of secondary particles of type a-either photons (γ), neutrinos (ν) or electrons (e)-is given by (Koldobskiy et al. 2021) a pp (E) = c n gas in the case of an homogeneous source, where n gas denotes the number density of target protons. The differential cross-section dσ a /dE is adopted from the Kamae et al. parametrisation (Kamae et al. 2006) below a threshold energy of 4 GeV, and above that energy by the AAfrag parametrisation (Koldobskiy et al. 2021).
For inverse Compton radiation the γ-ray emissivity is determined by (Schlickeiser 2002) ic where γ min = E 2m e c 2 1 + 1 + m 2 e c 4 E E and F KN (q, Γ) is previously defined in the inverse Compton energy loss rate.
The Bethe-Heitler process emissivity is given by (Zheng et al. 2016) where n p is the differential CR proton number density and τ −1 BH is the energy loss rate as introduced in Eq. (A10). Another important process to take into account in terms of the production of secondary electrons and positrons is γγ pair production (γ + γ −→ e + + e − ) as the corona becomes optically thick at high energy. Here, γ-rays that are in our case predominantly generated by the CR protons via hadronic pion production (with an emissivity γ pp ) interact with the thermal background photons (which are Comptonized up to X-ray energies within the corona) and introduce another source of electrons/ positrons. According to Böttcher et al. (2013) the steady state production rate of these e + /e − pairs can be estimated by γγ (γ e ) = f abs Ẽ 1 · γ pp Ẽ 1 + f abs Ẽ 2 · γ pp Ẽ 2 (B21) whereẼ 1 = γ e /f γ andẼ 2 = γ e /(1 − f γ ) denote the characteristic dimensionless γ-ray energies (in units of m e c 2 ), and the absorption fraction f abs ≡ 1 − [1 − exp(−ω γγ (E))]/ω γγ (E) is determined from the given optical depth ω γγ (E) due to this process (see Eq. C24). We assume that one of the produced particles obtains the major fraction f γ of the photon energy, where f γ = 0.9 is found from Monte Carlo simulations. Hence, an electron/positron pair is produced with energies γ 1 = f γ E and γ 2 = (1 − f γ )E.

C. ABSORPTION COEFFICIENTS
The effective optical thickness for an homogeneous medium is given by ω = R (α syn + α ff + α γγ ), where, for the considered astrophysical environments, we need to account for the synchrotron-self absorption (α syn ), the free free absorption (α ff ) and the γγ pair production (α γγ ) process.
The synchrotron self-absorption can be described by (Schlickeiser 2002) α syn (ν, θ) = − c 2 8πm e c 2 ν 2 ∞ 1 dγ e γ 2 e P syn (ν, γ e ) d dγ e γ −2 e n e (γ e ) where n e (γ e ) denotes the differential number density of CR electrons, which are assumed to be isotropic, and P syn (ν, γ e ) is the total spontaneously emitted spectral synchrotron power of a single relativistic electron that has already been introduced in Eq. (B12).
The thermal bremsstrahlung (free free) radiation emitted by an electron moving in the field of an ion-such as given by Eq. (B16)-can subsequently also get absorbed by the associated absorption process. This so-called free free absorption coefficient is given by (Rybicki & Lightman 1991) where θ and n gas stand for the gas temperature and density, respectively. We assume a quasi-neutral gas with about the same number of ions and electrons.
In contrast to the previous processes the γγ pair production process attenuates the photon intensity typically at γ-ray energies. The corresponding absorption coefficient for the γγ pair can be approximated by (Dermer & Menon 2009) where r e is the classical electron radius,Ẽ 1 = E/(m e c 2 ) denotes the dimensionless photon energy passing through a background of photons with energyẼ, n γ Ẽ is the isotropic photon field andφ s 0 cm is defined as ϕ(s 0 cm ) = 2 Here s 0 cm ≡ẼẼ 1 and the γγ pair production cross section σ γγ dependent on the dimensionless interaction energy s cm is given by σ γγ (s cm ) = 1 2 πr 2 e 1 − β 2 cm 3 − β 4 cm ln where β cm = (1 − γ −2 cm ) 1/2 = 1 − s −1 cm and γ cm = √ s cm denotes the center-of-momentum frame Lorentz factor of the produced electron/ positron.