A publishing partnership

The following article is Open access

Solving the Multimessenger Puzzle of the AGN-starburst Composite Galaxy NGC 1068

, , , , and

Published 2022 November 2 © 2022. The Author(s). Published by the American Astronomical Society.
, , Citation Björn Eichmann et al 2022 ApJ 939 43DOI 10.3847/1538-4357/ac9588

Download Article PDF
DownloadArticle ePub

You need an eReader or compatible software to experience the benefits of the ePub3 file format.

0004-637X/939/1/43

Abstract

Multiwavelength observations indicate that some starburst galaxies show a dominant nonthermal 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 γ-ray and neutrino emission. In this work, a homogeneous, steady-state two-zone multimessenger model of the nonthermal 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 multimessenger data—from radio to γ-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.

Export citation and abstractBibTeXRIS

Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

1. Introduction

The high rate of star formation and supernova (SN) 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 nonthermal emission component originates in their central regions, indicating the presence of an active supermassive 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 that may contribute to the extragalactic CR component observed at Earth, whose origin remains unknown.

The Fermi Large Area Telescope (LAT) has detected high-energy γ-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 γ-ray energies with imaging atmospheric Cherenkov telescopes (VERITAS Collaboration et al. 2009; Acero et al. 2009; Abdalla et al. 2018; Acciari et al. 2019), confirming that at least a fraction of the starburst galaxy population accelerates 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 AGN-starburst composite NGC 1068 has been seen in 10 yr of data from the IceCube Neutrino Observatory: an excess of neutrinos has been found in a region centered 0fdg35 away from the coordinates of NGC 1068 after a catalog-based search using 10 yr 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 nonthermal emission from the different acceleration sites. In the case of NGC 1068, there are strong observational hints for multiple emission sites: first, in the near-infrared to radio, the ALMA experiment has observed a strong flux cutoff toward smaller frequencies emerging from the inner parsecs of that source (Garcia-Burillo et al. 2016, 2019). Further nonthermal radio emission is associated with the more extended region of a few hundred parsecs (Wilson & Ulvestad 1982; Sajina et al. 2011). Second, the γ-ray flux extends out to 100 GeV (Acciari et al. 2019).

Indications that a one-zone model is not enough to explain the multimessenger emission from NGC1068 have been present for several years. Using radio and γ-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 γ-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. 2004, 2006) it has been suggested (Lenain et al. 2010) that the observed γ-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 γ-ray emission that has been suggested by Lamastra et al. (2019) is 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 several TeV which, however, can be excluded due to the upper limits of 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 model; a high-energy neutrino signal with a missing γ-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 γ-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 multimessenger 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 multimessenger emission of the source. In Section 3 we present our fitting approach as well as its results and conclude in Section 4.

2. The Two-zone AGN-starburst Model

The multimessenger observations of NGC 1068 dictate the need for multiple emission zones. In this work, we will capture its inner microparsecs 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 Figure 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. 2004, 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 nonthermal, 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 power-law distribution in momentum space up to a certain maximal kinetic energy , which depends on the competing energy-loss timescales in these zones. In the case of the starburst zone, we suppose that a certain fraction of the total energy of the SN that releases about 1051 erg and occurs with an approximate rate 5 (Veilleux et al. 2005) dependent on the IR luminosity LIR gets accelerated into CRs according to diffusive shock acceleration (DSA; e.g., Drury 1983; Protheroe 1999) by individual supernova remnants (SNR; 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 therein) as a result of a large number of core-collapse SNe. These winds introduce another source of acceleration 6 (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 finj ≪ 1 of the mass accretion rate , 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 nonthermal primary electrons are normalized by the nonthermal proton rates due to the requested quasi-neutral total charge number of the injection spectra of primary CRs above a characteristic kinetic energy of (Schlickeiser 2002; Eichmann & Becker Tjus2016; Merten et al. 2017). Note that this corresponds to a quasi-neutral acceleration site; however, CR transport can subsequently remove CR electrons and protons in different amounts from the nonthermal energy regime. Nevertheless their charge stays conserved. Transforming the source rates from momentum space into kinetic energy T, we obtain

for the injected nonthermal protons, and

for the nonthermal electrons (e) and positrons (e+). Here, the latter term 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 nonthermal 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 the case of the relativistic electrons is given by the inverse of the sum of the synchrotron (syn), inverse Compton (IC), nonthermal Bremsstrahlung (brems), and Coulomb (C) loss rates, according to

and in the 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 oriented 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 rescattered 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 E0 = 8.2 meV (Yun & Carilli 2002). The constant dilution factor Cdil is determined from the observed IR luminosity LIR according to the relation . 7 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 parameterization depends on the Eddington ratio (Lbol/LEdd), i.e., the ratio of the bolometric over the Eddington luminosity, and we adopt the relation of Hopkins et al. (2007) to determine Lbol based on the intrinsic X-ray luminosity LX between (2 and 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(H2) = 102.9 cm−3, and a molecular hydrogen mass of M(H2) ∼ 3.5 × 108 M.

More details on the individual energy-loss timescales can be found in Appendix A.

Figure 1. Refer to the following caption and surrounding text.

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 microparsecs, whereas the starburst ring is located at about 1 kpc distance from the nucleus.

Standard image High-resolution image

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 vw; 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 avis and the Keplerian velocity . 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 by 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 SNR via 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 vsh equals the wind speed vw. Thus, we use the acceleration timescale (e.g., Romero et al. 2018; Murase et al. 2020)

and set the maximal CR energy according to that value where the acceleration timescale exceeds the competing total loss timescale, i.e., .

Note that due to the huge differences in the physical parameters, such as the gas density or the magnetic field strength, the different timescales differ significantly between the starburst and the corona zone, as shown in Figure 2. With respect to the maximal CR energy it is shown that in the 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 a few keV. 8

Figure 2. Refer to the following caption and surrounding text.

Figure 2. The different timescales of relativistic protons (upper panel) and electrons (lower panel) for the starburst zone (left) and the AGN corona zone (right). The assumed parameters are given in 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 processes are as follows: hadronic pion production, 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.

Standard image High-resolution image

2.1. 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 nonthermal electron and proton density n(T) in the starburst as well as the coronal zone, as shown in Figure 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 the case of the coronal region. This is due to significantly smaller energy-loss timescales of the relativistic electrons according to synchrotron and IC losses (see Figure 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 cutoff introduced by the acceleration. In Figure 3 we consider a parameter scenario that results from a fit to the data (as introduced in Section 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 will be introduced in the following.

Figure 3. Refer to the following caption and surrounding text.

Figure 3. The resulting CR spectra of the starburst zone (left) and the AGN corona zone (right) using the best-fit parameters as introduced in Section 3, hereby s = 2.2 for the starburst and s = 1.7 for the corona.

Standard image High-resolution image

2.2. Nonthermal Emission

Using the spectral energy distribution of CR protons and electrons we determine the nonthermal emission of those particles from radio to γ-ray energies. Hereby, both emission regions are considered to be spherically symmetric, with the inner coronal region constrained to a radius , where typically rcor ∈ [1, 100] and denotes the Schwarzschild radius, and the surrounding starburst ring extends between an inner radius Rico-Villas et al. (2021) and an outer radius . In the following we include the nonthermal photon emission by synchrotron radiation, inverse Compton scattering, nonthermal bremsstrahlung from the CR electrons according to a spectral emissivity epsilonsyn, epsilonic, and epsilonbrems, respectively, as well as hadronic and photopion production by the CR protons with a spectral emissivity epsilonpp and epsilonpγ , respectively. These hadronic processes also introduce high-energy neutrinos as well as secondary electrons and positrons with a source rate , where the last two terms introduce electrons/positrons by γ γ and Bethe–Heitler pair production.

In addition to these nonthermal emission processes, we also account for the free–free emission (epsilonff) 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 γ-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 Appendices B and C.

For a homogeneous source for which epsilonm 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 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 . 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 Multimessenger Data

3.1. Data and Fit Procedure

In the following we use the previously introduced model to explain the multimessenger data of NGC 1068. Hereby, we quantify the goodness of the fit by the chi-squared value , where Pj (Oj ) 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 9 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 vw 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 LX = 1043 L43 erg s−1 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 modeling 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 , 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 L43 = 0.9. In general, it is expected (e.g., Mayers et al. 2018) that a higher X-ray luminosity corresponds to a higher black hole mass MBH than a low X-ray luminosity. As shown by different observations (see][and references therein]GRAVITY2020 the expected range of the black hole mass yields MBH ∈ [0.8, 1.7] × 107 M and since a higher X-ray luminosity generally indicates a higher value of MBH (e.g., Mayers et al. 2018), we use MBH(L43 = 0.9) = 0.8 × 107 M and MBH(L43 = 7) = 1.7 × 107 M, respectively, yielding an Eddington luminosity of LEdd(L43 = 0.9) = 1.0 × 1045 erg s−1 and LEdd(L43 = 7) = 2.1 × 1045 erg s−1, respectively. And based on the relation by Hopkins et al. (2007) we obtain a bolometric luminosity of Lbol(L43 = 0.9) = 0.2 × 1045 erg s−1 and Lbol(L43 = 7) = 2.9 × 1045 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 as well as the mass accretion rate decrease for a decreasing LX value. The surrounding starburst region is dominant in the IR due to the 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 LIR 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).

Table 1. Summary of Fixed and Free Parameters of the Radiation Zones of NGC 1068.

Starburst
s ngas (cm−3) B (μG) η (kpc) (kpc)ϰ Tinj (MeV) vw (km s−1) LIR (1011 L)
1.5–2.50.01–0.2510–100010–5001–1001.3–1.80.75/30.011001.5
AGN corona
s finj ngas (109cm−3) B (kG) η Rcor ϰ Tinj (MeV) μvis LX (1043erg s−1)
1.5–2.50.00001–0.10.05–500.001–101–1005–5005/30.010.1(0.9, 7)

Note. 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 SN energy that gets converted to CRs, , the gas density, ngas, the magnetic field strength, B, the radius of the starburst ring, , 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 table.

Download table as:  ASCIITypeset image

But also for the other, nonfixed 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 (Merloni & Fabian 2001; Ricci et al. 2018) constraining the gas density according to ngas = ωT/(σT R) if the coronal plasma is not dominated by electron–positron pairs. And based on the virial gas temperature θgas = GMBH mp/(3Rcor kB) = mp c2/(6rcor kB) of the protons we are able to draw constraints on the coronal magnetic field strength according to , where a low plasma beta (β ∼ 0.1–3) is expected from numerical MHD simulations (Miller & Stone 2000; Jiang et al. 2014, 2019). 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:

  • (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) reobserved the large-scale emission of that source with VLA at 2, 6, and 20 cm using 6farcs5 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 —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 from 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 3 orders of magnitude larger than the inner corona, that we are modeling here. Thus, there is a chance that the so-called extended corona region at a distance of about a milliparsec to a parsec 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 ×106 K of this extended corona is actually able to explain the Very Long Baseline Array (VLBA) radio observation. Despite the recent objections (which we will discuss in more detail in Section 4) by Baskin & Laor (2021), we adopt this simple approach using a gas with an (electron/proton) density of 2.5 × 105 cm−3 and an electron temperature of 106 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 Figure 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 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 beam sizes of 20 and 60 mas, respectively (Garcia-Burillo et al. 2016; Impellizzeri et al. 2019). 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 size 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 parsecs. Hence, the flux prediction at 224 GHz is expected to become smaller than the one from Impellizzeri et al. (2019) at 256 GHz, if the mismatched coverage is taken into account properly. So, we account heuristically for this effect by introducing an additional, lower uncertainty of the observed flux Fobs dependent on the beam size δbs. 10
  • (iii)  
    The γ-ray data are taken up to 100 GeV from the fourth Fermi-LAT catalog of γ-ray sources (Fermi-LAT collaboration 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 modeling by Ricci et al. (2017) of the torus absorption 11 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 Figure 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, which 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.

3.2. 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 for an intrinsic coronal X-ray luminosity of LX = 7(0.9) × 1043 erg s−1. Even though the resulting best-fit spectra are almost equal for both of those two cases, except for the resulting neutrino flux at ≲1 TeV. The resulting best-fit parameter space of the corona shows some differences: For LX = 0.9 × 1043 erg s−1 the inner corona needs to extend about , 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 finj 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 LX = 7 × 1043 erg s−1, the Figure 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 evaluations if the algorithm converges toward 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 multidimensional 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 marginalize over the other dimensions, which is done by using its minimal chi-squared value . 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 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 . Note that the observed black hole mass of NGC 1068 is rather small compared to other AGN yielding a small Schwarzschild radius, so that . Still, in combination with a high gas density of ∼ 1010 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 the case of the starburst region a high gas density of ngas ∼ 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.

Figure 4. Refer to the following caption and surrounding text.

Figure 4. The chi-squared distribution of the starburst zone (left) and the AGN corona zone (right) dependent on the fit parameters, where we display for two different parameters and marginalize over the others. The red cross marks the best-fit parameter values that are also used in Figures 2, 3 and 5.

Standard image High-resolution image

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 LCR = 9.55(2.07) × 1043 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 PCR/Pgas = 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 CR-to-thermal gas density pressure ratios (of 12% and 6% in case of a high and low X-ray luminosity, respectively) for the corona. Figure 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 toward small frequencies in the far-infrared. 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 multiwavelength data. At about a 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 for 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 best-fit case of a low coronal X-ray luminosity (LX = 0.9 × 1043 erg s−1) the resulting neutrino flux at ≲1 TeV is about a factor of five smaller than what is shown in Figure 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.

Figure 5. Refer to the following caption and surrounding text.

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.01 arcsec, and black markers indicate a beam size of ≳1 arcsec. The dark gray 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 gray area indicates the thermal IR emission by dust grains of the starburst region. Note that the torus attenuation is not taken into account here.

Standard image High-resolution image

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 (finj ∼ 0.06) and a smaller radius () 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 (ngas ∼ 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 (PCR/Pgas ≃ 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.

4. Conclusions and Discussion

In this work, we introduced a spatially homogeneous, spherically symmetric, steady-state two-zone model for AGN-starburst composite galaxies. Using the multimessenger data of NGC 1068 from the radio up the γ-ray band as well as its recent indications of high-energy neutrino emission, we present the 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., Inoue et al. 2020b; Murase et al. 2020; 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., Miller & Stone 2000; Ricci et al. 2018). 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 the case of a strongly magnetized corona and a starburst ring with a high SN rate (∼0.5 yr−1). Such a full multimessenger 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 (Inoue et al. 2020b; Murase et al. 2020; 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 multimessenger data of NGC 1068 can only be explained if we account for the nonthermal 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 optical 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 106 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 × 105 K at 5 GHz. In addition, they showed that a hot free–free emitting gas (with electron temperatures of ≳107 K) also overproduces the observed X-ray luminosity of NGC 1068. Therefore, they exclude optically thin free–free emission in NGC 1068, on the subparsec scales. We noticed, that due to the given maximal extension of this region an electron temperature of at least 106 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 the 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 subtorus structures. As the coronal region is a perfect CR calorimeter, we do not expect any additional nonthermal 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 substructure 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. 2004, 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 Yoast-Hull et al. 2014; Eichmann & Becker Tjus 2016). 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 nonthermal 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 et al. 2017).

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:

where T stands for the relativistic kinetic energy of a CR particle with a rest energy E0. 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 nZ , 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 characterizing 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 the 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 ne 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 ultrarelativistic 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)

with and . The terms 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 . 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 particle energy-loss rate can be approximated in the range 1.2 GeV < E ≤ 108 GeV by Krakau & Schlickeiser (2015)

where ngas 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 toward 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 particle energy-loss rate of Zheng et al. (2016)

where 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 and correspond respectively to 1 MeV/mec 2 and 2γp E, while =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 is the fraction of energy lost by the ultrarelativistic proton (γp ≫ 1 and βp 1) in the interaction, therefore the inelasticity of the collision.

Appendix B: Emissivities

The emissivities (in units of cm−3 s−1eV−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 photohadronic pion production as well as Bethe–Heitler pair production of CR protons; (iii) free–free emission of the thermal gas; and (iv) γ γ pair production.

The synchrotron emission by CR protons is negligible (see Equations (A4) and (A7)) so that only synchrotron radiation by CR electrons is considered in the following. Its spectral power is given by (e.g., Blumenthal & Gould 1970)

with P0 = 2.65 × 10 1 0 eV s−1 Hz−1 and νs = 4.2 ×106 Hz. The isotropic spontaneous synchrotron emission coefficient of the relativistic electron distribution yields

where ne (γe) is the differential CR electron density.

The emissivity function for bremsstrahlung radiation produced by CR electrons is (Stecker 1971)

where σbrems = 3.38 × 10−26 cm2.

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 ultrahigh energetic cosmic ray proton spectra, is described by the following expression (Dermer & Menon 2009):

Hereby, for secondaries generated by single pion production we use μb, , , whereas for those produced in multi pion production we adopt μb, , . 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) and is the differential intensity of CR protons with a Lorentz factor . The dimensionless target photon energy Eγ is described in units of Ee,0.

Table 2. Multiplicities ζ and Mean Fractional Energies χ of Secondaries Formed in Photomeson Production

SpeciesSingle π Multi π
Neutrinos
Electrons
γ-rays

Download table as:  ASCIITypeset image

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 gff is the so-called Gaunt factor, which is given by Padmanabhan (2000):

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)

in the case of an homogeneous source, where ngas denotes the number density of target protons. The differential cross section d σa /dE is adopted from the Kamae et al. parameterization (Kamae et al. 2006) below a threshold energy of 4 GeV, and above that energy by the AAfrag parameterization (Koldobskiy et al. 2021).

For inverse Compton radiation the γ-ray emissivity is determined by Schlickeiser (2002)

where

and FKN(q, Γ) is previously defined in the inverse Compton energy-loss rate.

The Bethe–Heitler process of emissivity is given by Zheng et al. (2016)

where np is the differential CR proton number density and is the energy-loss rate as introduced in Equation (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 ) 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 Bottcher et al. (2013) the steady-state production rate of these e+/e pairs can be estimated by

where = γe/fγ and = γe/(1 − fγ ) denote the characteristic dimensionless γ-ray energies (in units of mec 2), and the absorption fraction is determined from the given optical depth due to this process (see Equation (C3)). 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.

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

Synchrotron self-absorption can be described by Schlickeiser (2002)

where ne(γe) denotes the differential number density of CR electrons, which are assumed to be isotropic, and is the total spontaneously emitted spectral synchrotron power of a single relativistic electron that has already been introduced in Equation (B1).

The thermal bremsstrahlung (free–free) radiation emitted by an electron moving in the field of an ion—such as given by Equation (B5)—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 ngas 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 re is the classical electron radius, denotes the dimensionless photon energy passing through a background of photons with energy , is the isotropic photon field, and is defined as

Here and the γ γ pair production cross section σγ γ dependent on the dimensionless interaction energy scm is given by

where and denotes the center-of-momentum frame Lorentz factor of the produced electron/positron.

Footnotes

  • 5  

    Supposing a SN rate (note that Condon1992 suggested a value of 0.04 instead of 0.02 for normal galaxies) where the star formation rate SFR ≃ 17[LIR/(1011 L)] Myr−1.

  • 6  

    Note that in these phenomena also stochastic diffuse acceleration may become relevant due to the presence of a turbulent plasma within the wind bubbles.

  • 7  

    A more accurate approach to the IR photon spectrum has been proposed by Casey (2012), where a coupled modified greybody plus a mid-infrared power law has been used, but these modifications have no impact on our results.

  • 8  

    Note that we do not account for the possible steepening of the turbulent power spectrum at energies below the thermal proton energy, which would lengthen the acceleration time considerably.

  • 9  
  • 10  

    Note that this additional uncertainty reduces the resulting χ2-value of the fit, but hardly affects the resulting best-fit parameters.

  • 11  

    Including the combined effect of photoelectric absorption and Compton scattering by neutral material as well as the absorption by ionized gas using the ZXIPCF model (Reeves et al. 2008).

10.3847/1538-4357/ac9588
undefined