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

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.

In this work, we explore a possible causal relation between the 2017 gamma-ray flare and the delayed radio flare of TXS 0506+056. We consider that the jet accelerates a population of electrons at the miliparsec scale, emitting a short gamma-ray flare, and subsequently meets a configuration of electron clouds at the parsec scale, leading to fresh particle acceleration that describes the delayed radio flare. We apply the model to the 2017 flare of blazar TXS 0506+056, fitting the parameters to the time-domain multi-wavelength emission, including radio data. We show that the data are well described by propagation of the initially loaded jet material and subsequent interaction with a configuration of three clouds of cold electrons at a distance of several parsecs. By fitting the gamma-ray flare and the multi-frequency time-domain radio data, we constrain the particle density evolution along the jet and the spatial configuration of the parsec-scale clouds.

Our time-domain model consists of a series of blobs launched from the inner jet, at a distance x0,BH from the SMBH. At its original location x0,BH, each blob has a radius r0′ and is permeated by a magnetic field of local strength B0′. We assume the blob to be initially located at a distance of a few times the BLR radius, x0,BH ≳ RBLR, consistent with surveys of AGN jet formation and with radiative source models. In order to estimate the continuously varying time-domain emission, we simulate a series of blobs launched sequentially from x0,BH in temporal increments of ∆tBH = (2r0′ /c)/Γ. In each blob, we assume that pre-accelerated electrons are injected with a power-law energy spectrum. We model the interactions of these electrons and the resulting radiative cooling using the publicly available time-dependent numerical code AM3. The code treats all particle species and magnetic fields as homogeneously and isotropically distributed in the plasma rest frame, and solves the time- and energy-dependent equations describing their cooling and respective electromagnetic emission, including non-linear cascades triggered by pair production. The blob propagates downstream with a constant Lorentz factor Γ, expanding at speed ηc. This expansion results in a reduction in the particle densities and the adiabatic cooling of the electrons and positrons. We assume the magnetic field strength to evolve as a simple power-law in xBH: B ′ (x) = B0′ × (xBH /x0,BH)−p. For a magnetic field structure lying between purely poloidal and purely toroidal, the index p lies in the range 1 ≤ p ≤ 2. As the results will show, our best-fit scenario favors a value closer to p ≳ 1. Since simulating the downstream propagation of all of these blobs throughout the entire duration of the flare would make the model fitting unnecessarily computation-intensive, we explicitly simulate the particle interactions in only a small number of representative blobs, and we assume that neighboring blobs emit identically to their nearest representative. For our ∼3.4 year injection period (as measured in the observer’s frame), we found that explicitly simulating seven representative blobs, equally spaced time over this time period, provides sufficient temporal resolution to adequately sample the flare profile, given the rate of variation of the parameters along the jet. To estimate the total time-dependent emission spectrum in the observer’s frame, we sum the contributions of each of the blobs (the representative ones as well as the chain of contiguous copies), shifted by a time interval in the observer’s frame of ∆t = (1 + z)(2r0′ /c)/δD compared to the previous blob in the chain, thus capturing their sequential launching. Finally, we account for the effect of photon interactions with the extragalactic background light (EBL), which leads to an additional attenuation of the observed gamma-ray flux. We adopt the frequency- and redshift-dependent attenuation factors as tabulated in the gammapy database assuming the EBL model by Franceschini et al. Following the above pipeline, we obtain the total multi-wavelength flux from the simulated inner jet flare as a function of photon frequency and time on the observer’s frame.

As the power in non-thermal electrons injected in the inner jet decreases exponentially, further particle injection in the parsec-scale jet is required to describe the radio flare. Reproducing the non-trivial behavior of the multi-frequency radio light curves requires a complex description of this particle injection profile. Our model explains the radio flare in the following way: when the first (leading) blob in the chain propagating from the inner jet reaches the parsec-scale jet, it runs into clouds of higher particle density, stationary in the radial direction in the SMBH frame. In the process, the jet accelerates a fraction of the electrons in these clouds to a non-thermal spectrum. We model this by injecting an additional electron power-law spectrum in the leading blob. Given that the acceleration mechanism and its efficiency in the parsec-scale jet is likely to differ from that in the inner jet, we allow for a different power-law index of the non-thermal electrons, γcloud, compared to the inner jet. We found empirically that a fine-grained description of the multi-frequency radio data requires a minimum of three clouds of different sizes, partially overlapping in the radial direction. We parameterize the size of the clouds in the direction transverse to the jet motion as a fraction of the jet’s local cross-sectional radius: rcloud,0 = ξ0 r, ξ0 < ξ1; rcloud,1 = ξ1 r, ξ1 < 1; rcloud,2 = r. Clouds 0 and 1 can be smaller than the jet’s cross-section, while cloud 2 is defined as the largest one, with a size of at least the jet’s cross-section. Modeling for simplicity the cross-section of each cloud as circular, while crossing cloud 0(1) the jet sees a cross-section of π[ξ0(1) r(t)]2, while for cloud 2 it is simply πr2 (t). Assuming the jet accelerates a constant fraction of the cloud particles, and given that the electron luminosity is dominated by the highest-energy electrons, the luminosity of non-thermal electrons at each point x (in the SMBH rest frame) follows approximately, Le,cloud i (x) ∝ Γβc (π rcloud squared ne (x)) Ee max, where ne is the cloud density at location x. As for simplicity we consider a constant value of Γ, the non-thermal electron luminosity injected in the leading blob at each point x is simply proportional to the local cloud particle density. We model this profile as an asymmetric Gaussian for each cloud i. The parameters L0,cloud,i, ρ, and ω(l)r,i correspond to the peak value, width, and asymmetry scale of the electron luminosity distribution in the leading blob caused by the crossing of each of the clouds. To simulate the multi-wavelength emission from the leading blob as it crosses our three-cloud configuration, we simulate independently the particle interactions in the overlapping volume between the blob and each of the clouds, evolving them independently in time and summing their contributions to the total time-dependent emission.

Following the above jet modeling procedure, we obtain the total time-domain photon flux, dNtotal (t)/dEγ dt = dNjet (t)/dEγ dt + dNcloud (t)/dEγ dt + dNsteady /dEγ dt, where we adopt the constant steady-state component, dNsteady /dEγ dt, from a previous work describing the quiescent-state spectrum from the extended jet. We then fit this model by comparing this flux with the RATAN-600 and Fermi-LAT light curves for TXS 0506+056. We start with an initial parameter guess based on a by-eye fit of the broad data features, such as the approximate two-and-a-half year delay between the gamma and radio flares and the order-of-magnitude luminosity in both wavelengths. We then use this initial guess as a seed for MCMC sampling. We use 500 independent sampler chains, each performing 80 steps, evolved with the Metropolis-Hastings algorithm. We evaluate approximately ∼40,000 model realizations. We treat the first 30% of steps as a burn-in and parameter exploration. We explore the parameter ranges listed in Tab. I, informed by general principles: for example, we assume a bulk Lorentz factor Γ ≲ 10, as determined by previous studies of radio morphology, a spectral index p ≳ 1, as expected generally from Fermi processes, and a magnetic field strength in the inner jet of order of magnitude B ∼G, as determined in previous source models. Because we restrict the sampling to the neighborhood of a single high-likelihood solution rather than the entire parameter space allowed by the priors, we can interpret the reported spread as a local likelihood-based uncertainty rather than a formal Bayesian credible interval derived from a global posterior sampling. We assume measurement uncertainties are Gaussian and independent.

We show the best-fit multi-wavelength light curves in the observer’s frame in Fig. 3: the upper five panels show the multi-frequency radio light curves, while the two lower panels show the gamma-ray light curve in the two analyzed frequency bands. As we can see in Fig. 3, the model captures the non-trivial frequency dependence of the radio flare, owing to the complex spatial profile of the clouds swept by the parsec-scale jet. At 2 GHz, shown in the upper panel, the model correctly describes a late peak, about 3.5 years after the peak of the gamma-ray flare. In contrast, higher radio frequencies display an earlier rise phase with a more temporally extended flare containing two sub-peaks; as we discuss below, this is due to the jet crossing the earlier clouds. The decay phase is shorter than the rise phase, with a duration of only about a year, and is roughly simultaneous across all frequencies, in agreement with the data.

To physically interpret the predicted radio light curves, we show in Fig. 4 the best-fit spatial distribution of the electron injection luminosity in the parsec-scale jet resulting from the crossing of each of the clouds. This quantity is proportional to each cloud’s spatial density profile, as the jet is assumed to have constant speed and accelerate a constant fraction of the particles that it picks up. We can see that the clouds lie at a radial distance of between 3 and 7 pc. All clouds impinge on the jet a peak non-thermal electron power of ∼ 2 × 1044 erg/s. The most upstream cloud is the most extended one in the radial direction, starting at about 2 pc and peaking just below 6 pc. At the same time, it has a relatively small transversal radius (cf. ξ0 Tab. II), making the emission region relatively compact. As our leading blob crosses this first cloud, the synchrotron emission from the freshly accelerated electrons describes the slow early rise observed between 5 and 22 GHz. This is controlled by the parameter ωl,0. Owing to the compactness of the emitting region, the escaping photons are partially self-absorbed; this explains the fact that the 2 GHz flux does not initially increase as fast as at higher frequencies. The jet blob then crosses the second and third clouds, both peaking just beyond 6 pc. While the second cloud is also compact (cf. ξ1 Tab. II), the third cloud is larger, as described in Sec. III B, and the corresponding emitting region spans the entire jet cross section. This makes the region more optically thin, leading to the late 2 GHz flare. The simultaneous emission at higher frequencies contains a contribution from both cloud 2 and cloud 1. Finally, as the blob moves away from the cloud complex, the radio flux decreases simultaneously at all frequencies. This decay time is determined by the parameters ωr,1 and ωr,2, which describe the decrease of the spatial density profile of the last two clouds. Altogether, the model predicts that the bulk of the radio emission during the flare is entirely confined to the innermost ∼7 pc, roughly compatible with measurements of the VLBI core.

With a best-fit Doppler factor of δD = 7.8 for a viewing angle of 5◦ (both assumed constant), the cloud distribution responsible for the radio emission is constrained to lie between 2 and 7 pc. The first cloud, lying between 2 and 5 pc, has the smallest cross-sectional overlap with the jet, yielding a compact emission zone. This explains the early rise in the radio flux at and above 5 GHz, while at 2 GHz the emission is absorbed via SSA. As the jet crosses the second and third clouds, the emission at and above 5 GHz undergoes a double peak in the time domain, about 2.6 and 3.2 years after the gamma-ray flare. As the third cloud intersects with the entire jet cross-section, the emission zone is the most optically thin. This describes the long delay observed in the 2 GHz band, where the peak occurs only about 3.4 years after the gamma-ray flare. While the model does not directly capture neutrino emission, it constrains the conditions necessary to produce a non-trivial multi-frequency radio flare. We thus probe source properties that may be key to unlocking a broader and more consistent picture of multi-messenger emission from blazars and jetted AGN at large. More broadly, our results show that fitting time-domain multi-wavelength data allows us to increase the complexity of source models, pushing them beyond the state-of-the-art single-zone framework. Including more detailed observational data, such as time-domain VLBI imaging, may eventually allow the inclusion of even more complex geometries, enabling new probes of the physics of AGN jets.

Improvements to AI Systems:

1. Spatio-Temporal Multi-Zone Radiative Modeling

• Extend current single-zone radiative codes (e.g., AM3) to handle spatially resolved, time-evolving jet geometries with multiple interacting particle populations (inner jet blobs + parsec-scale clouds).

• Implement dynamic coupling between particle injection, adiabatic expansion, magnetic field decay, and external cloud interactions in a moving-frame simulation.

• Enable automatic generation of multi-frequency light curves (radio to gamma) with self-consistent synchrotron self-absorption (SSA) and EBL attenuation.

2. Bayesian Inference with Physics-Informed Priors for Jet Parameters

• Build an MCMC-based fitting engine that jointly constrains jet kinematics (Lorentz factor, Doppler factor, viewing angle), magnetic field profile (power-law index p), particle injection spectra, and cloud spatial density profiles (asymmetric Gaussians with position, width, asymmetry).

• Incorporate priors from VLBI morphology, Fermi-LAT spectra, and radio core size to reduce degeneracies.

• Provide local likelihood uncertainties and posterior predictive checks for time-domain data.

3. Automated Cloud Configuration Discovery

• Develop an AI module that iteratively proposes and tests multi-cloud configurations (number, size, radial position, density profile) to match observed delayed radio flares with sub-peaks.

• Use gradient-free optimization (e.g., differential evolution) to explore cloud parameters, then refine with MCMC.

• Output physical constraints (e.g., cloud distances 2–7 pc, density contrasts) that can be compared with VLBI core measurements.

4. Time-Delay Causality Inference

• Train a model to automatically identify causal links between high-energy flares (gamma/neutrino) and delayed lower-frequency emissions (radio) by fitting the full time-domain multi-wavelength dataset.

• Quantify the lag distribution across frequencies and infer the acceleration site (e.g., parsec-scale cloud crossing) from the light-curve morphology.

5. Multi-Messenger Event Association

• Integrate neutrino event data (e.g., IceCube) as additional constraints: given a gamma-ray flare, predict the probability of a delayed radio flare and the required hadronic content, using the fitted electron/proton ratios and radiation densities.

• Enable real-time classification of neutrino-blazar associations based on predicted radio light-curve signatures.

---

What the Improved AI System Can Do:

• Simulate the full electromagnetic evolution of a flaring blazar from sub-parsec to parsec scales, reproducing complex multi-frequency light curves (e.g., TXS 0506+056) with delayed radio peaks and sub-peaks.

• Fit observational data (radio, optical, X-ray, gamma) simultaneously, yielding posterior distributions for jet and cloud parameters (e.g., magnetic field index p≈1, cloud distances 2–7 pc, Doppler factor ≈7.8).

• Predict the timing and spectral shape of future radio flares following a gamma-ray flare, enabling targeted multi-wavelength follow-up campaigns.

• Constrain the physical conditions (particle densities, magnetic fields, acceleration efficiencies) in parsec-scale jet regions, directly testable with VLBI and future CTA observations.

• Assess the likelihood of neutrino production in specific blazar flares by linking the fitted electron distributions to hadronic models, improving multi-messenger source identification.

• Automatically discover new jet-cloud interaction geometries from time-domain data, potentially revealing hidden structures in AGN jets beyond current single-zone assumptions.