An Evolving Leptonic Jet Model for Delayed Radio Flares in Neutrino Blazars

Abstract: The jets of blazar active galactic nuclei (AGN) are promising sites of hadron acceleration and subsequent neutrino production, owing to their extreme intrinsic power and high radiation density. Potential associations of IceCube neutrinos with blazars displaying delayed radio flares, such as TXS~0506+056 and PKS~1424+240, may support this scenario. However, the mechanisms and location of particle acceleration in the jet remain unclear. The 2017 IceCube event associated with TXS 0506+056 was concurrent with the initial peak of a flare in gamma rays, followed by a longer flare at radio frequencies peaking three years later. State-of-the-art single-zone radiative source models focus on the compact high-energy-emitting region, and thus fail to describe this subsequent radio enhancement. In this work, we model the time-domain multi-frequency evolution of TXS 0506+056 by considering electron emission along the physically extended jet. We couple a numerical particle interaction framework with a dynamic description of the temporal and spatial jet evolution down to the parsec scale. We fit the model to multi-wavelength data, including multi-frequency radio light curves. The results suggest that the parsec-scale jet meets a population of older, cooled electrons near the edge of the observable radio core, establishing a causal relation between the initial high-energy emission and trailing radio activity. This spatially and temporally resolved modeling approach can offer new insights into the physical conditions along the jet of flaring neutrino candidate blazars.

• The paper develops a spatially and temporally resolved leptonic jet model that jointly fits the 2017 Fermi-LAT flare and multi-frequency RATAN-600 radio light curves of TXS 0506+056.
• The model requires at least three partially overlapping electron clouds between 2 and 7 parsecs, reproducing frequency-dependent radio peaks from roughly 2.6 to 3.4 years after the gamma-ray flare and implying a Doppler factor of 7.8.
• The fit favors a magnetic-field decline close to B′ ∝ x⁻¹, while its under-predicted very-high-energy gamma rays and leptonic-only treatment leave proton acceleration, neutrino production, and the absolute radio-emission location unresolved.

The 2017 coincidence of IceCube-170922A with a gamma-ray flare from TXS 0506+056 established that blazar jets are plausible neutrino sources, but it also exposed a modeling gap: the radio flare that followed peaked roughly three years after the neutrino and gamma-ray event, at frequencies down to 2 GHz, while standard single-zone radiative models—tuned to compact, high-energy-emitting regions—are structurally incapable of reproducing such delayed, frequency-structured radio emission. In "An Evolving Leptonic Jet Model for Delayed Radio Flares in Neutrino Blazars" (2608.12696), Kochocki, Rodrigues, and Whitehorn address this gap with a spatially and temporally resolved leptonic jet model applied to TXS 0506+056, fitting both the 2017 Fermi-LAT gamma-ray light curve and multi-frequency RATAN-600 radio light curves within a single causal framework.

Observational motivation

The paper assembles two time-domain datasets. The first consists of RATAN-600 radio light curves of TXS 0506+056 at 2.3, 4.7, 8.2, 11.2, and 22.3 GHz [Sotnikova_2022], binned coarsely (167-day bins) for visualization, together with derived spectral indices between the lowest pair, highest pair, and all bands. Two features drive the modeling: lower frequencies peak later than higher ones, and the spectral index hardens during the rise and steepens during decay while never reaching the fully self-absorbed value $\alpha_{\mathrm{SSA}} = 5/2$—behavior the authors attribute to a superposition of regions with different optical depths rather than any single optically thick zone.

The second dataset is a reanalysis of public Fermi-LAT data in two energy bands (100 MeV–5 GeV and 5–500 GeV), using 121-day bins over December 2009 to November 2024, with a standard Pass 8 P8R3 binned likelihood analysis and test-statistic-based upper limits. The authors note that the nearby sub-GeV flaring source PKS 0502+049 could cause confusion before 2016 but not during the period of interest.

Model construction

The model couples the time-dependent lepto-hadronic code AM$^3$ (Klinger et al., 2023) to an explicitly evolving jet geometry. A chain of blobs is launched sequentially from $x_{0,\mathrm{BH}} = 3.5\times10^{17}$ cm—a few BLR radii outside the broad-line region—with each blob expanding at $\eta c$, propagating at constant Lorentz factor $\Gamma$, and embedded in a magnetic field declining as $B' \propto x_{\mathrm{BH}}^{-p}$. Electrons are injected with a power-law spectrum whose normalization follows an asymmetric double-exponential function of observer-frame time, peaking at $T_{\mathrm{peak}} = 1.35$ yr. Radiative processes include synchrotron emission, SSA, inverse Compton scattering on local electrons plus external BLR and dusty-torus photons, $\gamma\gamma$ pair production with cascades, adiabatic cooling, and EBL attenuation of the escaping gamma rays.

The central physical hypothesis concerns the parsec-scale jet: when the leading blob reaches several parsecs, it encounters three partially overlapping, radially stationary clouds of cold electrons. The jet accelerates a fraction of these electrons into a non-thermal power law (with index $\gamma_{\mathrm{cloud}}$ allowed to differ from the inner-jet index), producing fresh synchrotron emission that constitutes the delayed radio flare. Each cloud's density profile is parameterized as an asymmetric Gaussian in radius, with transverse sizes expressed as fractions $\xi_0 < \xi_1 < 1$ of the local jet cross-section; cloud 2 spans the full cross-section by construction. The authors found empirically that no fewer than three clouds are required to reproduce the fine structure of the multi-frequency light curves.

Nineteen parameters are sampled via Metropolis-Hastings MCMC (500 chains, ~40,000 model evaluations) around a by-eye seed solution, with uniform priors over physically motivated ranges. The steady-state baseline flux is taken from a prior extended-jet leptohadronic model of this source [Rodrigues:2025cpm], supplemented by an additional one-zone electron population fitted by eye to archival radio data.

Results

The best-fit solution achieves good agreement across all seven light curves. The key inferred quantities are:

Parameter	Best-fit value
Bulk Lorentz factor $^3$0	4.54
Initial blob radius $^3$1	$^3$2 cm
Expansion rate $^3$3	0.14
Inner magnetic field $^3$4	1.15 G
Magnetic field index $^3$5	1.01
Injected electron index $^3$6	1.94
Cloud electron index $^3$7	1.11
Cloud radial extent	2–7 pc

The corresponding Doppler factor is $^3$8 for the assumed 5° viewing angle. The fit implies a specific geometric interpretation of the flare's frequency structure. The most upstream cloud (peaking just below 6 pc) has a small transverse overlap ($^3$9), creating a compact, partially self-absorbed region whose synchrotron emission drives the early rise at ≥5 GHz while suppressing the initial 2 GHz response. The second and third clouds, peaking beyond 6 pc, produce a double-peaked structure at higher frequencies about 2.6 and 3.2 years after the gamma-ray peak. Because the third cloud spans the entire jet cross-section, its emission zone is the most optically thin, accounting for the late 2 GHz peak roughly 3.4 years after the gamma-ray maximum. The bulk of the radio emission is confined to the innermost ~7 pc, broadly consistent with VLBI core measurements [Kun_2018, Ros:2019bgo].

A notable structural result is the best-fit magnetic field index $x_{0,\mathrm{BH}} = 3.5\times10^{17}$0, near the poloidal limit of the allowed range—an outcome also obtained independently in the proton-synchrotron scenario of Rodrigues et al. [Rodrigues:2025cpm]. The implication is that the fit favors a jet in which the toroidal field component grows only weakly out to parsec scales, a constraint on jet magnetization structure that emerges purely from time-domain photometry.

Relation to prior frameworks

The paper positions itself against three classes of models. Single-zone models fail categorically because the compactness required for high-energy emission renders the zone opaque to low-frequency radio. Continuous extended-jet treatments [Lucchini:2021scp, Zacharias:2022pea, Rodrigues:2025cpm] reach the radio band but under-predict flux below a few tens of GHz, which the authors interpret as evidence that a homogeneous flow with continuous acceleration misses either a discrete particle-loading event, a shock-driven acceleration mechanism more efficient at low magnetization, or an optically thin escape channel such as a slower sheath. A contemporaneous work [stathopoulos2026delayedradioflaresneutrinoassociated] also attributes the delayed flare to parsec-scale particle acceleration but invokes jet deceleration for the timing and predicts a single exponential profile peaking simultaneously across frequencies; the present model instead reproduces the observed frequency-dependent peak times through spatial superposition of distinct electron populations, which simpler frameworks cannot easily replicate.

For neutrino astrophysics, the result establishes a concrete causal chain: the inner-jet perturbation responsible for the 2017 gamma-ray flare propagates outward and triggers fresh acceleration at parsec scales. Since delayed radio flares have now been characterized in multiple neutrino-candidate blazars [kochocki2026characterizinggammaradiodelayedflaring], this framework offers a way to constrain the environments in which any associated hadronic acceleration would occur—even though the present model deliberately excludes protons.

Limitations and open questions

The authors are explicit about several caveats. The computational cost of evolving particle populations over years restricts the MCMC to a local exploration around a single seed solution; the reported spreads are local likelihood uncertainties, not globally representative credible intervals, and alternative solutions cannot be excluded. Standard degeneracies persist among $x_{0,\mathrm{BH}} = 3.5\times10^{17}$1, viewing angle, and the absolute distance between gamma-ray and radio emission sites, compounded by intrinsic leptohadronic degeneracies between $x_{0,\mathrm{BH}} = 3.5\times10^{17}$2, $x_{0,\mathrm{BH}} = 3.5\times10^{17}$3, and $x_{0,\mathrm{BH}} = 3.5\times10^{17}$4.

Two tensions deserve emphasis. First, the model under-predicts the observed 5–500 GeV gamma-ray flux during the flare, which the authors connect to hybrid scenarios in which protons dominate above ~100 GeV; verifying the findings under proton-inclusive assumptions remains open. Second, recent VLBI evidence for an ultra-fast spine surrounded by a slow sheath in TXS 0506+056 [Kovalev:2026fba] would place the radio-emitting plasma at distances comparable to the gamma-ray region but off-axis, whereas this one-dimensional coaxial model places it at 2–7 pc along the axis. The authors argue that the relative configuration and sizes of the emitting zones are robust because they are fixed by the multi-frequency timing and opacity evolution, while absolute positions are not. Additional unmodeled effects include possible jet curvature or precession, orientation-dependent Doppler boosting, and gravitational lensing amplification suggested for this and other neutrino-associated blazars [Britzen:2025bww].

Conclusion

This work demonstrates that time-domain, multi-frequency fitting can push blazar modeling beyond the single-zone paradigm in a controlled way: by propagating an electron-loaded jet from the BLR boundary to parsec scales and requiring interaction with a three-cloud density profile, the model reproduces both the 2017 gamma-ray flare and the full frequency-dependent structure of the delayed radio flare of TXS 0506+056, constraining the cloud complex to lie between 2 and 7 pc with individual non-thermal electron injection powers of order $x_{0,\mathrm{BH}} = 3.5\times10^{17}$5 erg/s. The cost is a thirty-parameter descriptive framework whose global uniqueness is untested, and whose leptonic restriction leaves the hadronic—and hence neutrino—content of the flare unconstrained. The specific open questions are whether proton-inclusive versions of this geometry remain viable given the GeV-band deficit, and whether time-domain VLBI imaging can break the distance degeneracies that currently limit absolute localization of the radio-emitting zones.