Bayesian parameter study of the Seyfert-starburst composite galaxies NGC 1068 and NGC 7469

Abstract: Multimessenger observation of the Seyfert-starburst composite galaxies NGC 1068 and NGC 7469 indicate a characteristic feature in the radio band (the so-called mm-bump) as well as indication of high-energy neutrinos by the AGN corona. Moreover, also the starburst ring of these sources is bright in the radio and hence, a potential source of $γ$-rays and neutrinos. We aim to explain the non-thermal features of these two sources with our homogeneous steady-state Seyfert-starburst composite model, which we refined in this work. Hereby, we account for stochastic diffuse acceleration and energy losses within the corona and $γγ$-pair attenuation of the escaping $γ$-rays. Since the non-thermal features of Seyfert sources contribute only marginally to the electromagnetic spectrum, only few data points can be assigned to the starburst ring or the AGN corona. Hence, prior information on the physical parameters is incorporated within a Markov Chain Monte Carlo approach to avoid overfitting. Based on this Bayesian parameter study we show, that the non-thermal features of NGC 1068 can be explained well. Still a more detailed treatment of the spatial inhomogeneities in the central region of the AGN could further improve the fit results. This manifests itself even more clearly in the case of NGC 7469, where the mm-bump needs to emerge from a coronal size $R_{\rm c}&gt;100\,\mathcal{R}<em>{\rm s}$, whereas (TeV-PeV)-neutrino emission requires $R</em>{\rm c}&lt; 10\,\mathcal{R}_{\rm s}$. Similar to what has previously been shown in other wavebands, our analysis highlights that the spatial extension of the so-called AGN corona depends the considered energy of the messenger. Hence, it seems that there is not a unique edge of the corona and a substantial progress in the understanding of these phenomena is expected if future analysis account for these spatial inhomogeneities.

• The paper develops a Bayesian MCMC model combining stochastic particle acceleration in AGN coronae with shock acceleration in starburst rings, including secondary-pair cascades and extended gamma-ray attenuation.
• For NGC 1068, the model reproduces the broadband spectrum and favors a coronal radius near 35 Schwarzschild radii, but typically underpredicts the observed TeV–PeV neutrino flux by about a factor of two when fitting the millimeter emission tightly.
• For NGC 7469, the results show that a homogeneous corona cannot simultaneously match the millimeter bump and candidate neutrinos above 100 TeV, indicating that energy-dependent or radially structured emission regions are likely required.

Motivation and scope

The two Seyfert-starburst composite galaxies NGC 1068 and NGC 7469 are among the most compelling candidates for high-energy neutrino production in AGN coronae. IceCube reported a $4.2\sigma$ neutrino signal at TeV energies from NGC 1068, while for NGC 7469, Sommani et al. found two IceCube real-time track alerts with a chance coincidence probability of only 0.08% (3.2$\sigma$), at energies significantly above those of the NGC 1068 events. This paper updates a homogeneous steady-state composite model — originally introduced to explain NGC 1068's full non-thermal spectrum — by refining the treatment of the AGN corona and applying a Bayesian Markov Chain Monte Carlo (MCMC) analysis to both sources simultaneously.

The central physical picture treats each source as two quasi-independent zones: a compact AGN corona on sub-milliparsec scales near the supermassive black hole, where stochastic diffusive acceleration (SDA) operates in the same environment as radiative losses; and a kiloparsec-scale circumnuclear starburst ring, fed primarily by supernova remnant shocks via diffusive shock acceleration (DSA). The model accounts for $\gamma\gamma$-pair attenuation of escaping gamma-rays against coronal X-ray targets.

Model refinements relative to prior work

Three key upgrades distinguish this treatment from the earlier E+22 formulation:

1. Self-consistent SDA in the corona: The momentum diffusion term is retained alongside continuous and catastrophic losses, yielding a substantially harder non-thermal proton spectrum than the simplified injection-only approach of E+22. This resolves a known tension: E+22 required an unusually soft proton spectrum inconsistent with expectations from stochastic acceleration.
2. Nonlinear secondary $e^\pm$ cascades: Secondary pairs from hadronic processes and from leptonic radiation (synchrotron, inverse Compton) are included iteratively until convergence. These become dynamically relevant below several hundred MeV.
3. Extended $\gamma\gamma$ attenuation: Beyond internal attenuation within the corona, the authors compute external attenuation using the radially decreasing X-ray photon field outside the corona, plus a "shielded" scenario where gamma-ray production is confined to $(1-r_{\rm sh})R_{\rm c}$ inside an X-ray-emitting shell.

A notable constraint on the turbulence strength parameter ($\eta^{-1} \lesssim 0.9$) arises from consistency with the quasilinear approximation used in the transport equation.

Bayesian methodology

Because the non-thermal signatures attributable specifically to either the starburst ring or the corona occupy only narrow spectral windows — radio/mm emission and GeV–TeV gamma rays (NGC 1068 only), plus TeV–PeV neutrinos — the likelihood is severely underconstrained without informative priors. Truncated normal priors in log-space encode physical or observational bounds on parameters such as the gas-to-CR pressure ratio $f_{\rm gas}$, plasma beta $\beta_{\rm B}$, Thomson optical depth $\omega_{\rm T}$, coronal radius $\sigma$0 (in units of $\sigma$1), starburst magnetic field $\sigma$2, gas density $\sigma$3, and CR injection spectral index $\sigma$4. A uniform prior applies to $\sigma$5 over $\sigma$6.

To prevent ALMA mm data (with very small quoted uncertainties) from dominating the likelihood, a floor of 15% relative uncertainty is imposed on all data points unless otherwise stated. New observational inputs include:

• Revised NGC 1068 starburst flux: Using FIRST survey data at 1.43 GHz after removing jet contamination from an ellipse of $\sigma$7, the authors find that only about 7% of the total 1.43 GHz flux originates from the starburst ring — roughly a factor of three less than adopted in E+22.
• New VLA L-band measurement of NGC 7469 at 1.5 GHz ($\sigma$8 Jy).
• Fermi-LAT upper limits for NGC 7469: Analysis of ~15.5 years of data yields TS consistent with zero, confirming no significant GeV detection and providing 95% CL spectral limits.

Results: NGC 1068

The updated model reproduces the broadband non-thermal emission of NGC 1068 well. The GeV gamma-rays originate predominantly from the starburst ring, with increasing coronal contribution below $\sigma$9GeV energies. The posterior distribution for the starburst-ring CR injection index yields a strikingly hard value of $\gamma\gamma$0, far harder than the DSA expectation near $\gamma\gamma$1 and even harder than the prior mean of $\gamma\gamma$2. This suggests either energy-independent escape from accelerators or co-spatial acceleration and radiation within the ring itself. The inferred coronal radius $\gamma\gamma$3 (i.e., $\gamma\gamma$4) matches previous independent findings from Mutie et al., lending confidence to the spatial scale.

However, a systematic tension remains: for most best-fit scenarios the predicted TeV–PeV neutrino flux falls about a factor of two short of the IceCube measurements when the observed mm-bump is accurately reproduced. Relaxing the minimal uncertainty on the mm-data to $\gamma\gamma$5 allows higher values of $\gamma\gamma$6 and lower $\gamma\gamma$7, shifting the synchrotron self-absorption peak upward and enabling perfect agreement with IceCube — but this comes at the cost of loosening constraints on the coronal structure.

Results: NGC 7469

For NGC 7469, the situation is more constrained. Two independent MCMC fits were performed (coronal and starburst components decoupled, justified by negligible cross-talk given the assumed coronal origin of the neutrinos). Key outcomes:

• The starburst ring explains the large-scale radio data with $\gamma\gamma$8 — closer to canonical DSA expectations than in NGC 1068 — but requires a magnetic field $\gamma\gamma$9G, slightly above prior expectation, driven by the need to match primary synchrotron emission while respecting IC gamma-ray limits.
• The predicted GeV gamma-ray flux sits about a factor of two below current Fermi-LAT upper limits, leaving room for future detection or tighter constraint.
• Most critically, a homogeneous spherical corona cannot simultaneously explain the mm-bump and the candidate high-energy neutrino signal: matching the former requires $e^\pm$0, whereas producing neutrinos at $e^\pm$1TeV demands $e^\pm$2 absent shielding. Even relaxing uncertainties to 30% does not reconcile this discrepancy under current gamma-ray limits.

This size mismatch mirrors a known pattern across wavebands: X-ray reverberation studies favor $e^\pm$3 while SED-based mm modeling often yields $e^\pm$4. The paper argues the gap exceeds what time variability alone could account for, since even the most extreme reported mm variability implies only a factor-of-two size change.

Limitations and open questions

Several caveats qualify these conclusions:

• Spatial homogeneity assumption: The model assumes uniform density, magnetic field, and photon fields in each zone. The persistent tension in NGC 1068's neutrino normalization and NGC 7469's radius conflict likely stem from unmodeled spatial gradients in the central engine. The authors explicitly state that accounting for such inhomogeneities should improve fits, particularly for NGC 7469.
• Neglect of turbulent cascade damping: The analysis fixes the magnetic turbulence level and assumes non-thermal particles draw negligible energy from the flow. While verified self-consistently for the best-fit scenarios via the criterion $e^\pm$5, if damping were significant the resulting proton spectra would soften considerably.
• Simplified coronal geometry: The spherical, isotropic treatment ignores possible disk-like or stratified coronal structure, and the shielded scenario is treated as a limiting case rather than a physically motivated configuration.
• Neutrino spectral assumptions: The IceCube flux estimate adopts a power-law with index $e^\pm$6 for NGC 1068; alternative spectra peaking at different energies could shift the required normalization.
• Jet contribution excluded: For NGC 1068, jet-related extended emission is deliberately removed from the starburst flux extraction but not modeled explicitly.

An open question left by the work concerns whether a radially structured corona — with charged-particle confinement varying with Larmor radius so that the radio-emitting zone extends well beyond the $e^\pm$7-ray/neutrino zone — can quantitatively resolve the size discrepancy in both sources simultaneously.

Conclusion

This Bayesian refinement of a Seyfert-starburst composite framework demonstrates that NGC 1068's multimessenger phenomenology is broadly reproducible with moderate CR pressure ($e^\pm$8) and a coronal radius near $e^\pm$9, though a residual neutrino deficit persists under strict mm-bump fitting. For NGC 7469, the analysis establishes that no single homogeneous coronal scale satisfies both the mm-radio constraint and the candidate TeV-scale neutrino signal, pointing toward energy-dependent spatial extent of the emitting region as a necessary feature of any successful model.