© The Authors 2026

1. Introduction

A significant fraction of active galactic nuclei (AGNs) launch collimated, relativistic jets that can extend from parsec to megaparsec scales, transporting enormous amounts of energy, up to ∼10⁴⁶ erg s⁻¹ in kinetic power (Blandford et al. 2019), from the vicinity of the supermassive black hole into the host galaxy and beyond. These jets play a central role in galaxy evolution and feedback, regulating star formation and heating cluster gas (McNamara & Nulsen 2007; Cavagnolo et al. 2010; Fabian 2012; Blandford et al. 2019). Among these, radio galaxies host jets misaligned with respect to the observer’s line of sight (i.e., seen at large angles to the jet axis), which are classified according to the Fanaroff-Riley (FR) scheme (Fanaroff & Riley 1974). This classification divides sources into FRI and FR II classes, based on radio morphology and luminosity. Fanaroff-Riley I jets are less luminous and show centrally peaked radio emission, whereas FR II jets are brighter and terminate in compact hot spots at the outer edges of their radio lobes (e.g., Fanaroff & Riley 1974; Bridle & Perley 1984). This was later interpreted in terms of FR I jet deceleration via mixing (Bicknell 1984, 2002). Unified schemes further link these morphological classes to orientation and accretion states, with FR I and FR II radio galaxies as the misaligned counterparts of BL Lacertae (BL-Lac) type objects and radio-loud quasars, respectively. Despite this phenomenological framework, the physical process responsible for the FR I-FR II dichotomy remains uncertain (e.g., Bicknell 2002; Wold et al. 2007; Tchekhovskoy & Bromberg 2016; Gourgouliatos & Komissarov 2018; Perucho 2019, 2020), with key questions about jet energetics and composition still unanswered.

One critical factor in distinguishing FR I and FR II sources is jet kinetic power. Fanaroff-Riley II jets typically carry higher power (≳10⁴⁵ erg s⁻¹) than FR I jets (e.g., Blandford & Konigl 1979; Godfrey & Shabala 2013). However, this is not the only intervening parameter; stellar composition and host galaxy properties also play a role, mainly at intermediate powers, as pointed out by, for example, Mingo et al. (2022). In practice, AGN jet power spans many orders of magnitude and is only indirectly inferred, with large uncertainties implied. Power and composition together govern how a jet interacts with its environment: high-power jets can remain collimated and relativistic to large distances, whereas lower-power jets are prone to entrainment and deceleration. However, a stability analysis of relativistic flows shows that thermodynamically cold (specific enthalpy dominated by rest-mass energy) and fast (kinetically relativistic) jets are intrinsically more stable (Perucho et al. 2005, 2010), and whereas powerful jets have more inertia (or mass flux), these properties do not necessarily need to be correlated, mainly at compact (sub-kiloparsec) scales. Still, all these quantities remain poorly constrained by observations.

As mentioned above, estimates of the AGN jet kinetic power remain challenging due to large systematic uncertainties (e.g., Ineson et al. 2017; Hardcastle & Croston 2020). Common methods either use scaling relations between radio synchrotron luminosity and kinetic power (Godfrey & Shabala 2015; Croston et al. 2018) or derive the work done in inflating X-ray cavities (McNamara & Nulsen 2012; Pasini et al. 2021). Radio-based estimates depend on assumptions about particle content, magnetic field strength, spectral shape, and source age and may vary by orders of magnitude depending on the presumed energy partition or proton content (Godfrey & Shabala 2015). X-ray cavity measurements offer a more direct probe of the jet energy budget but require well-resolved structures. They assume that the observed cavities fully trace jet energy with minimal leakage (McNamara & Nulsen 2012), in addition to an estimate of the jet age, which may itself be uncertain by more than one order of magnitude. Both methods provide time-averaged powers and may miss effects such as episodic activity, which is increasingly observed in remnant and restarted radio sources as array sensitivity increases (Brienza et al. 2017; Shabala et al. 2020).

The composition of jets is equally elusive but critically important for jet dynamics. Jets might be made of electron-positron pairs or pair-proton plasmas, and the presence of heavy particles greatly increases flow inertia (Scheck et al. 2002; Blandford et al. 2019) if they represent a significant fraction of the bulk flow mass flux. Observational probes of composition (such as Faraday rotation or lobe pressures) suggest a range of possibilities (e.g., Croston et al. 2018) at large scales but are often compatible with both leptonic and hadronic solutions. The question of jet composition is also intertwined with the potential emission of very energetic neutrinos (Mannheim 1993; Murase et al. 2014; Murase & Stecker 2023), as relativistic protons are necessary to produce them. Stars in the jet could be the missing link between mass-loading and PeV neutrino production (Fichet de Clairfontaine et al. 2025a).

Interestingly, FR I radio galaxies require a larger thermal population in their lobes than FR IIs, which is in agreement with the former being significantly mass-loaded (Croston et al. 2018). Theoretical studies show that light, pair-dominated jets decelerate strongly as soon as they entrain even a relatively small baryonic load. By contrast, jets that are already proton-rich carry far more kinetic power per unit radio luminosity and are much less affected by additional mass-loading (e.g., Komissarov 1994; Bicknell 1995). Importantly, jets interact with stars and interstellar medium (ISM) gas as they propagate: stellar winds, interstellar clouds, and small-scale instabilities at the jet surface can load jets with baryonic material (e.g., Komissarov 1994; Bowman et al. 1996; Perucho et al. 2014; Matsumoto et al. 2017; Gourgouliatos & Komissarov 2018; Perucho 2019, 2020; Anglés-Castillo et al. 2021). These processes were originally proposed to explain the FR dichotomy due to low-power jets being transonic and turbulent through mass-load, whereas high-power jets remain supersonic (Bicknell 1995). In this picture, mass-load and composition are intimately linked to the large-scale morphology of radio galaxies.

Furthermore, mass loading by stars has also been suggested to be connected to intriguing observational features recently identified. High-precision astrometry reveals small but systematic offsets between the apparent radio core (from VLBI) and the optical centroid (from Gaia measurements) in many AGNs (Plavin et al. 2019; Petrov et al. 2018). In particular, Plavin et al. (2019) found that in most cases, the optical position is displaced downstream along the jet direction by milliarcseconds (mas). This “radio-optical offset” implies that the optical emission centroid is often dominated by extended jet light rather than by the nucleus, in contrast to the expected transition controlled by opacity in the jet. Extended parsec-scale optical jets shifting the Gaia photocenter downstream, contamination by host galaxy optical emission, and accretion disk dominance producing upstream offsets have been proposed as possible explanations (Kovalev et al. 2017; Petrov et al. 2018; Kovalev et al. 2020b).

However, Fichet de Clairfontaine et al. (2025b) (hereafter FdC25) showed that stellar mass loading can naturally generate positive radio-optical offsets through the dissipation it induces in the jet if a fraction of the dissipated energy is invested in particle acceleration. In their numerical experiments, jets with FR I-like powers (L_j ∼ 10⁴³ erg s⁻¹) entrained mass from realistic stellar distributions, which forced a gradual conversion of kinetic energy into internal energy along the flow. Radiative-transfer calculations with the code RIPTIDE (Fichet de Clairfontaine et al. 2021) showed that dissipation can increase the minimum Lorentz factor of the nonthermal electron population downstream, causing the optical synchrotron emission to peak farther from the jet base than the radio core. As a result, offsets of a few parsecs naturally arise. Under this assumption, namely that the dissipation-driven increase in $\gamma {\prime }_{min}$ shifts the optical centroid, FdC25 found that even moderate stellar mass-loading rates (∼10⁻⁹ M_⊙ yr⁻¹ pc⁻³) reproduce the typical offset magnitudes and alignments reported observationally (see also Kovalev et al. 2017). In essence, the offset is a direct imprint of where dissipation becomes most efficient: the optical emission traces the downstream dissipation zone, while the radio peak remains anchored to the optically thick jet base.

The recent mass-load models connect two previously separate ideas. Historically, stellar entrainment was proposed as a possible explanation for the slow and diffuse nature of FR I jets (e.g., Komissarov 1994; Bowman et al. 1996; Perucho et al. 2014). We invoke the same physical mechanism to explain radio-optical position shifts. The region in which these interactions could be more intense would coincide with the strong dissipation region identified by Laing & Bridle (2014). FdC25 suggest that these positional offsets could act as a novel probe of jet physics: they report that jets with the most pronounced positive offsets correspond to moderately low-power FR I jets, whereas the most powerful jets would show no offset or even negative offset (optical upstream of radio) because they remain relativistic and do not show indications of strong, continuous dissipation at scales of tens to hundreds of parsecs. In other words, the sign and magnitude of the offset depend on the balance of jet power and entrainment. If this interpretation holds, then measuring radio-optical offsets could place independent constraints on jet energetics, complementing traditional methods.

Motivated by these considerations, we undertook a systematic study of jet dynamics, composition, and resulting emission with a focus on radio-optical offsets. In particular, we extended the work of FdC25, limited to a single jet power, by exploring a range of jet kinetic powers and host stellar distributions. We parameterized the stellar density profile (using a core radius r_c, s) to represent different galactic environments. By varying the jet power and mass-loading profile, we aimed to quantify the interplay between these factors and their impact on jet deceleration and the relative positioning of radio and optical emission. Our goal was to assess whether observed radio-optical offsets can help constrain the kinetic power and composition of AGN jets, thereby shedding light on the long-standing questions surrounding the FR I and II classification and jet feedback.

The paper is organized as follows. Section 2 describes our numerical tools and setup. Then, in Sects. 3 and 4, we present our results related to jet dynamics and radiative output, respectively. In Sect. 5, we discuss these results within the framework of radio-optical offset observations and other observational evidence before presenting our conclusions in Sect. 6. Throughout this paper, quantities given in the jet co-moving frame are primed, and we assume a flat Λ Cold Dark Matter (ΛCDM) cosmology with H₀ = 69.6 km s⁻¹ Mpc⁻¹, Ω₀ = 0.29, and Ω_Λ = 0.71.

2. Numerical setup

2.1. Jet simulation method

The dynamical simulations were carried out with the quasi-1D relativistic magnetohydrodynamic (RMHD) code described Anglés-Castillo et al. (2021) and further used in FdC25. Under the approximations of a narrow jet (jet radius much smaller than its length) and a relativistic flow, these models can be built by solving the time-dependent (magneto-)hydrodynamical equations for the transversal flow with the time coordinate acting as the axial coordinate in the steady flow (quasi-one-dimensional approximation; Komissarov et al. 2015), with the appropriate boundary conditions at the jet–ambient medium interface. The code uses a Synge equation of state to describe a mixture of ideal, relativistic gases of protons, electrons, and positrons. Jets are injected with leptonic composition (leptonic fraction η ≡ ρ_e/ρ = 1), in pressure equilibrium with the ambient medium at the grid boundary (see Martí 2015; Martí et al. 2016), with constant initial flow density ρ_j, Lorentz factor γ_j, and axial ($B_{\mathrm{j}}^z$). The toroidal magnetic field is also fixed, and its profile is shown in Anglés-Castillo et al. (2021). The total jet power is defined as

$$\begin{array}{c}\hfill L_{\mathrm{j}}=\pi R_{\mathrm{j}}^2({\rho }_{\mathrm{j}}{h}_{\mathrm{j}}{\gamma }_{\mathrm{j}}^2+{\left(B_{\mathrm{j}}^{\varphi }\right)}^2)v_{\mathrm{j}},\end{array}$$(1)

where R_j is the jet radius, h_j is the specific enthalpy, and v_j is the flow velocity. For reference, we adopted the injection parameters of model J4_C from Anglés-Castillo et al. (2021): a cold leptonic jet with rest-mass density ρ_j, 0, Lorentz factor γ_j, 0 = 6, and injection radius R_j = 1 pc. In contrast to FdC25, where jet power was fixed at L_j = 10⁴³ erg s⁻¹ for an FR I-type value, we explored a wide range of jet powers from 10⁴¹ to 10⁴⁶ erg s⁻¹ (including both FR I and FR II regimes).

Throughout this work, we varied the jet kinetic power solely by modifying the rest-mass density at injection. All other jet parameters (particularly the Lorentz factor, specific enthalpy, and magnetic field strength) were kept fixed. This choice was motivated by the expectation that, at the scales studied, these parameters would remain close to the values adopted here, and that reaching jet powers of 10⁴⁶ erg s⁻¹ by increasing them would result in unrealistic values. Although mixed solutions (involving, e.g., simultaneous increases in rest-mass density and Lorentz factor) are also possible, we adopted the simplest approach to isolate the effect of the rest-mass density on the jet power, and therefore on the radio–optical offset. This choice was further supported by the results of Anglés-Castillo et al. (2021), who show that varying the injection Lorentz factor within the range γ_j = 5 − 10 (consistent with the typical values observed in parsec-scale jets; e.g., Lister et al. 2021) does not significantly alter either the deceleration distance or the global impact on the dynamics (see their Fig. 5). Regarding the magnetic field, Zamaninasab et al. (2014) showed that the jet magnetic flux correlates with accretion power over seven orders of magnitude across the AGN population, implying $B(1\mathrm{pc})\propto L_{\mathrm{j}}^{1/2}$ for a fixed jet cross-section; field strengths inferred from parsec-scale VLBI core-shift measurements (Lobanov 1997) span roughly one to two orders of magnitude across the population. However, by the parsec scales at which our simulations begin, jets have already converted the bulk of their Poynting flux into kinetic energy through magnetic acceleration (Vlahakis & Konigl 2004; Komissarov et al. 2007; Chatterjee et al. 2019), transitioning from a magnetically dominated flow near the black hole to a kinetically dominated one. At the tens-of-parsec to kiloparsec scales where stellar mass loading and dissipation operate, the magnetic contribution to the energy flux is most likely negligible compared to the bulk kinetic flux. Hence, fixing B_j at injection is thus a simplification, but it is not expected to alter our qualitative conclusions.

The computational grid covers 20 pc × 2000 pc in cylindrical coordinates (r, z) with 800 × 5000 cells, spanning along several stellar core radii (see Anglés-Castillo et al. 2021). This resolution is adequate to reach convergence, as shown in Appendix A. In total, we made 400 simulations (100 simulations with different jet powers L_j for two stellar core radii, r_c, s = 0.5, 1 kpc, and 100 simulations with different stellar core radii r_c, s for two jet powers L_j = 10⁴³, 10⁴⁴ erg s⁻¹).

We emphasize that the quasi-1D RMHD framework adopted here was designed to accurately capture the dynamics of relativistic, supersonic flows. In practice, the numerical description remains reliable as long as the jet bulk axial velocity satisfies v ≳ 0.9c (corresponding to Lorentz factors γ_j ≳ 2.3). Below this threshold, the assumptions underlying the steady, axisymmetric treatment become progressively less accurate, as the flow is expected to enter a strongly decelerated and potentially turbulent regime that cannot be faithfully captured by our present approach.

2.2. Ambient medium and mass load

The galactic ISM is modeled by a King-like profile, as in Anglés-Castillo et al. (2021),

$$\begin{array}{c}\hfill p_{\mathrm{a}}(z)=p_{\mathrm{a},0}{[1+{\left(\frac{z}{r_{\mathrm{c}}}\right)}^2]}^{\alpha },\end{array}$$(2)

with α = −1.095, the core radius r_c, base pressure p_a, 0, and density ρ_a, 0 chosen to ensure pressure balance at the jet base, following Anglés-Castillo et al. (2021). In practice, we set r_c = 1.2 kpc, as in previous studies (Perucho et al. 2014). Mass loading from stellar winds was introduced via a source-term profile that traces the stellar density. We took

$$\begin{array}{c}\hfill Q(z)=Q_0{[1+{\left(\frac{z}{r_{\mathrm{c},\mathrm{s}}}\right)}^2]}^{{\alpha }_{\mathrm{s}}},\end{array}$$(3)

with α_s = −1.095 and where r_c, s is the core radius of the host’s stellar distribution. Here, Q(z) (in units of g yr⁻¹ pc⁻³ Bowman et al. 1996) gives the mass injection rate per unit volume at height z. In our simulations, we varied r_c, s over the range 0.1 − 1.5 kpc to assess its effect on jet dynamics and emission. This represents another extension of the study of FdC25, where it was defined in the range r_c, s ∈ [1.0, 1.5] kpc. The base loading rate Q₀ was chosen to reflect a typical old stellar population in an elliptical galaxy, Q₀ ∼ 5 × 10²⁴ g yr⁻¹ pc⁻³ (which corresponds to an average mass-loss rate of 10⁻⁹ − 10⁻¹⁰ M_⊙ yr⁻¹ for a stellar density n_★ ∼ 1 pc⁻³), within the range studied in FdC25. The stellar winds were assumed to be initially cold (proton-electron plasma) and to be entrained into the jet flow (see, e.g., Bosch-Ramon et al. 2012; Perucho et al. 2017). Table 1 summarizes the key physical parameters for the jet, ambient medium, and mass loading in our simulation suite.

2.3. Radiative transfer with the RIPTIDE code

We computed synthetic synchrotron emission maps by post-processing our relativistic jet simulations with the RIPTIDE code (Fichet de Clairfontaine et al. 2021, 2022). In this model, the nonthermal electron population in the jet’s co-moving frame is prescribed as a power law,

$$\begin{array}{c}\hfill \frac{\mathrm{d}{n}_{\mathrm{e}}^{\prime }}{\mathrm{d}{\gamma }_{\mathrm{e}}^{\prime }}=K{\gamma }_{\mathrm{e}}^{\prime -p}\mathrm{for}{\gamma }_{\mathrm{e},\min }^{\prime }<{\gamma }_{\mathrm{e}}^{\prime }<{\gamma }_{\mathrm{e},\max }^{\prime },\end{array}$$(4)

with a spectral index of p ≃ 2.2 (a typical value for mildly relativistic shock acceleration; see e.g. Ostrowski & Bednarz 2002; Lemoine & Pelletier 2003). The normalization constant K and the low-energy cutoff $\gamma {\prime }_{\mathrm{e},min}$ were determined from the local nonthermal electron number and energy densities ($n{\prime }_{\mathrm{e}},e{\prime }_{\mathrm{e}}$) via the relations of Gómez et al. (1995):

$$\begin{array}{cc}& K={\left(\frac{e_{\mathrm{e}}^{\prime }(p-2)}{m_{\mathrm{e}}{c}^2(1-C_{\mathrm{E}}^{2-p})}\right)}^{p-1}{\left(\frac{1-C_{\mathrm{E}}^{1-p}}{n_{\mathrm{e}}^{\prime }(p-1)}\right)}^{p-2},\hfill \end{array}$$(5)

$$\begin{array}{cc}& {\gamma }_{\mathrm{e},\min }^{\prime }=\frac{1}{m_{\mathrm{e}}{c}^2}\frac{e_{\mathrm{e}}^{\prime }}{n_{\mathrm{e}}^{\prime }}\frac{p-2}{p-1}\frac{1-C_{\mathrm{E}}^{1-p}}{1-C_{\mathrm{E}}^{2-p}},\hfill \end{array}$$(6)

where $n_{\mathrm{e}}\prime$ is the number density of the nonthermal electrons and $e_{\mathrm{e}}\prime$ its corresponding energy density, and C_E is the fixed ratio $C_{\mathrm{E}}=\gamma {\prime }_{e,max}/\gamma {\prime }_{\mathrm{e},min}={10}^3$. We considered radiative cooling to be dynamically negligible at the simulated scales of hundreds of parsecs to kiloparsecs (see FdC25).

We assumed that a fixed fraction of the local thermal energy and leptonic number is transferred to the nonthermal electrons: $e{\prime }_{\mathrm{e}}={ϵ}_{\mathrm{e}}e{\prime }_{\mathrm{th}}$ and $n{\prime }_{\mathrm{e}}={\xi }_{\mathrm{e}}n{\prime }_{\mathrm{e},\mathrm{th}}$, with ϵ_e = ξ_e = 0.1 (e.g., Gómez et al. 1995; Mimica & Aloy 2012; Fromm et al. 2016). Here $e{\prime }_{\mathrm{th}}$ and $n{\prime }_{\mathrm{e},\mathrm{th}}$ are the co-moving internal energy and lepton density from the simulation. Hence, $\gamma {\prime }_{\mathrm{e},min}\propto e{\prime }_{\mathrm{th}}/n{\prime }_{\mathrm{e},\mathrm{th}}=\overline{ϵ}{\prime }_{\mathrm{e}}$, i.e., the mean internal energy per particle. As a consequence, $\gamma {\prime }_{\mathrm{e},min}$ increases in regions of strong dissipation (e.g., shocks or turbulent mixing regions) in which the internal energy is enhanced and efficient particle acceleration can occur.

This dependence is central to reproducing the observed radio-optical positional offsets (FdC25). In our model, higher-energy electrons are preferentially produced in regions where jet-star interactions enhance local dissipation, shifting the optical emission centroid downstream relative to the opacity-defined radio core. In contrast, adopting a fixed $\gamma {\prime }_{\mathrm{e},min}$ prescription (e.g., $\gamma {\prime }_{\mathrm{e},min}=1$) suppresses this spatial differentiation: the optical emissivity either remains co-spatial with the radio emission or falls below the observational detection threshold, depending on the local physical conditions. We note that even in the case of intrinsically co-spatial emission, systematic offsets may arise observationally because optical positions correspond to flux-weighted centroids, whereas VLBI radio positions trace the brightest compact core component. In our analysis, centroids are extracted from the point source function (PSF) convolved synthetic maps without full instrumental forward modeling and therefore represent idealized observational proxies.

The synchrotron emissivity j′_ν′ and absorption coefficients α′_ν′ were computed in the co-moving frame at each cell using standard formulae (e.g., Rybicki & Lightman 1979; Katarzyński et al. 2001). These coefficients were then transformed to the observer frame via the Doppler factor δ = [γ_j(1 − β_j cos θ_obs)]⁻¹ (with the bulk Lorentz factor γ_j at the cell and line-of-sight angle θ_obs), according to

$$\begin{array}{c}\hfill j_{\nu }={\delta }^2{j}_{{\nu }^{\prime }}^{\prime }\mathrm{and}{\alpha }_{\nu }={\delta }^{-1}{\alpha }_{{\nu }^{\prime }}^{\prime },\mathrm{with}\nu =\frac{\delta }{1+z}{\nu }^{\prime },\end{array}$$(7)

so that the emission is boosted in the direction of motion. We solved the radiative transfer equation dI_ν/dτ_ν = S_ν − I_ν (assuming a constant synchrotron source function, S_ν = j_ν/α_ν, within the cell) along each line of sight to accumulate the emission and absorption through all the cells (Rybicki & Lightman 1979). Relativistic beaming makes the approaching jet appear significantly brighter and more compact to the observer. We accounted for the counter-jet by mirroring the emission and deboosting it. However, for the range of observational angle study, the effect of these counter-jets on the radio and optical centroids is expected to be limited (the intensity ratio between the jet and counter-jet is expected to be ∼10³ and ∼10⁵ for θ_obs = 30° and 5°, respectively; see Ghisellini et al. 1993).

Finally, we generated synthetic radio (43 GHz, chosen as a representative VLBI frequency for parsec-scale AGN imaging) and optical (5 × 10¹⁴ Hz) synchrotron emission maps by projecting the emergent intensity onto the sky plane. In post-processing, we convolved the maps with representative Gaussian point-spread functions and computed image-plane flux-weighted centroids, from which we defined the apparent angular separation d_app between radio and optical positions (Plavin et al. 2019). This procedure provides an idealized, centroid-based estimate of the positional offset; we neither performed visibility-based VLBI forward modeling nor simulated group-delay astrometry. To mimic observational conditions, we applied an optical flux threshold consistent with Gaia DR3 sensitivity (F_AB ≳ 10⁻⁴ Jy at 3σ; Brown et al. 2021) and adopted an astrometric precision of ∼0.1 mas. A millijansky-level radio flux cut has a negligible impact. The resulting d_app values were compared with the observed radio-optical offsets following FdC25. The jet geometry and notation used throughout this work are illustrated in Fig. 1; an example of the resulting radio and optical surface-brightness maps, with the extracted centroid positions and the corresponding d_app, is shown in Fig. 2.

Fig. 1. Schematic of the jet geometry. The approaching jet propagates from the nucleus at angle θobs to the line of sight. The radio and optical centroids (zradio, zopt) are projected onto the sky plane, defining the apparent offset dapp into the sky plane from the observer’s point of view. The solid line in the jet marks z(v = 0.9c); the stellar distribution (rc, s) is shown in orange (not drawn to scale).

Fig. 2. Synthetic synchrotron surface-brightness maps at 43 GHz (left) and optical 5 × 1014 Hz (right) for a representative run with Lj ≃ 1043 erg s−1, rc, s = 500 pc, θobs = 30°, at redshift z = 1. The dashed lines mark the intensity-squared centroid positions zradio (blue) and zopt (red); the bracket to the right of the optical panel indicates the resulting dapp = 6.17 mas.

3. Jet dynamics

In this section, we investigate how the jet kinetic power, together with the spatial distribution of stellar mass loading, affects jet dynamics. For all simulated jets, we fixed the average mass-load rate Q₀ = 5 × 10²⁴ g pc⁻³ yr⁻¹, which is a typical value producing observable radio-optical shifts, as shown in FdC25.

To ensure the physical consistency of our analysis, all dynamical quantities were evaluated only up to the axial position z(v = 0.9c), defined as the distance from the jet base at which the bulk flow velocity first drops below 0.9c. This position sets the upper boundary of the region where the relativistic quasi-1D approximation adopted in our simulations remains valid. The position z(v = 0.9c) depends strongly on the jet power and on the stellar core radius: for jet powers above L_j ≃ 3 × 10⁴³ erg s⁻¹ the jet velocity remains above 0.9c throughout the simulation box for any value of r_c, s. Below L_j ≃ 3 × 10⁴³ erg s⁻¹, z(v = 0.9c) moves upstream rapidly for decreasing values of jet power and increasing stellar core radius. For example, at L_j ≃ 10⁴³ erg s⁻¹, we have z(v = 0.9c)∼500 pc, for any value of r_c, s. Below L_j ≃ 10⁴² erg s⁻¹, z(v = 0.9c)< 100 pc, and the jet opens abruptly near the base. That is why in the rest of the paper, we distinguish between simulations in which z(v = 0.9c)≥2 kpc, corresponding to the full extent of the computational domain, and those in which z(v = 0.9c)< 2 kpc. In the latter case, physical quantities were evaluated within the relativistic region z ≤ z(v = 0.9c) and, when required, extrapolated up to z = 2 kpc. These are indicated with triangles in the figures. This approach allows us to compare average jet properties over a common axial distance while preserving the internal consistency of the relativistic flow approximation. When extrapolation was applied, each quantity was treated separately, as detailed below. Specifically, over the fitting interval [0.2 z(v = 0.9c), z(v = 0.9c)], the jet radius R_j was extrapolated by fitting a conical profile R(z) = s z forced through the origin by least squares, and evaluating R(2 kpc) = 2000 s pc. A conical morphology is physically motivated for mildly relativistic, decelerating jets (Laing & Bridle 2014; Porth & Komissarov 2015). The mean electron energy $\overline{ϵ}{\prime }_{\mathrm{e}}$ and proton fraction $1-\overline{\eta }$ were each fitted with a power law A z^α in log-log space over the same interval, with the slope α determined by the fit, then averaged over the full 0 − 2 kpc range. In a conical jet without mass loading, adiabatic losses drive power-law scalings of macroscopic quantities with distance (Kaiser 2006); in the presence of stellar mass entrainment, the exact index was modified, but we adopted this functional form as a first-order approximation.

To facilitate the interpretation of our results, we studied the role of the jet power and the stellar spatial distribution separately. We analyzed the impact of these two parameters on three representative quantities: the jet radius at a height of z = 2 kpc (end of simulation box, top panel in Fig. 3), the average energy per lepton within the jet (central panel in Fig. 3), and the average jet baryonic (proton) mass fraction (third panel in Fig. 3), both averaged within the jet region where z(v > 0.9c), which provide information about the macroscopic structure and microphysical state of the flow. This analysis complements and extends the study run by FdC25 for jets with L_j = 10⁴³ erg s⁻¹.

Fig. 3. Left: Jet radius Rj at z = 2 kpc (top), average energy per lepton $\overline{ϵ}{\prime }_{\mathrm{e}}$ (second), average proton fraction within the jet (third), and axial deceleration distance z(v = 0.9c) (bottom) as a function of the jet power Lj, for two stellar core radii rc, s (see legend). Each symbol corresponds to one simulation. The circles indicate simulations where the jet remains relativistic throughout the full 2 kpc box; the triangles indicate simulations where the jet decelerates before 2 kpc, for which the jet radius is obtained by extrapolation (see Sect. 3). The dashed line in the bottom panel marks the 2 kpc box limit. Right: Same as left but shown as a function of the stellar distribution core radius rc, s for fixed jet powers (see legend).

The top panel of Fig. 3 (left) shows the jet radius at z = 2 kpc as a function of jet power for two values of r_c, s. For L_j ≤ 10⁴³ erg s⁻¹, jets undergo strong lateral expansion, reaching transverse radii in the range ∼10 − 300 pc, with only a weak dependence on r_c, s. In this regime, the flow is strongly affected by pressure imbalance and stellar mass loading, leading to efficient deceleration and pronounced decollimation (Perucho et al. 2014; Anglés-Castillo et al. 2021). Around L_j ∼ 10⁴³ erg s⁻¹, the jet radius at 2 kpc decreases rapidly with increasing jet power, marking a transition from a mass-loaded, decollimated regime to a more inertia-dominated flow. For higher jet powers, jets remain significantly collimated and exhibit characteristic radii of ∼2 − 3 pc at kiloparsec scales.

The second panel of Fig. 3 (left) shows the average energy per lepton within the jet as obtained from our post-processing modeling. For the lowest-power jets, it remains around ∼1 GeV, consistent with dissipation resulting from heavy mass entrainment. The average lepton energy decreases steadily when the baryonic loading becomes progressively less efficient as a dissipation channel for jet powers larger than 10⁴³ erg s⁻¹. For instance, the value is ∼10⁻² GeV at L_j = 10⁴⁵ erg s⁻¹ and approaches the electron rest-mass energy ($\overline{ϵ}{\prime }_{\mathrm{e}}\sim m_{\mathrm{e}}{c}^2$) around L_j = 10⁴⁶ erg s⁻¹. This trend reflects the decreasing ability of stellar mass loading to convert a significant fraction of the bulk jet kinetic energy into internal energy of the leptonic population as jet inertia increases.

The third panel of Fig. 3 (left) shows the corresponding entrained average proton fraction, $1-\overline{\eta }$, within the jet (η being the electron fraction, set to unity at the jet base). At low jet powers, the jet becomes heavily mass-loaded, with baryonic fractions approaching unity. As the jet power increases, the proton content drops sharply, reaching values ≲10% for L_j ≥ 10⁴⁴ erg s⁻¹. This reduction reflects both the reduced relative importance of entrainment and the enhanced inertia provided by the original leptonic jet content. The bottom panel shows the position z(v = 0.9c) where the jet begins to decelerate significantly inside the simulation box. An unperturbed jet shows z(v = 0.9c) = 2 kpc, which is the axial size of our simulation box.

The right panels of Fig. 3 show the same quantities as those on the left as a function of the stellar distribution core radius, for two intermediate jet powers. As in the left panels, the circles indicate whether the jet radius or the other quantities are directly measured at 2 kpc or averaged over the full jet, while the triangles indicate extrapolation of the radius, and averaged quantities within the jet region where z(v = 0.9c). The abrupt change in slope seen in the top panel around r_c, s ≃ 400 pc marks the transition from simulations where the jet remains relativistic throughout the full 2 kpc box (circles, direct measurement) to those that are significantly decelerated, and reach v = 0.9c before z = 2 kpc (triangles, extrapolated radius). This transition is primarily physical: a more extended stellar distribution, i.e., larger r_c, s, sustains mass loading over a larger axial range, increasing the cumulative momentum transfer and triggering earlier deceleration. This steep drop is visible in the bottom panel of Fig. 3, which shows z(v = 0.9c) dropping sharply from 2 kpc to ∼450 pc within a narrow interval of ∼80 pc in r_c, s. The jet radius data points shown as triangles should therefore be interpreted with care, as they represent an extrapolation from the relativistic region assuming a conical geometry.

The top right panel of Fig. 3 shows that the influence of r_c, s on the jet radius is not as strong as that of jet power. We see that jets with powers above 10⁴⁴ erg s⁻¹ are only weakly affected by entrainment from stellar winds, independently of the stellar distribution. At L_j = 10⁴³ erg s⁻¹, an increase in r_c, s from 100 to 1500 pc results in an increase in the jet radius from ∼2.5 to ∼8 pc. In contrast, significant changes can be observed in both average lepton energy (second panel) and proton fraction (third panel). The former rises from ∼10⁻¹ to ∼1 GeV, whereas the latter increases from ∼25% to nearly unity for a L_j = 10⁴³ erg s⁻¹ jet. Both magnitudes reflect sustained dissipation (with the consequent efficient particle acceleration) over larger interaction scales. The plots also indicate that the effects are much weaker for a jet power of 10⁴⁴ erg s⁻¹.

In summary, the obtained jet dynamics highlight the dominant role of jet kinetic power in setting the jet collimation and composition at kiloparsec scales, with a transition occurring over a narrow jet power range centered around L_j ∼ 10^43.5 erg s⁻¹ for the average mass-loading rate adopted here (Q₀ = 5 × 10²⁴ g yr⁻¹ pc⁻³). This transition shows only a weak dependence on the stellar spatial distribution. Nevertheless, different stellar populations and mass-loss histories are expected to shift the transition power, as discussed in FdC25.

4. Synthetic radiation and radio-optical offsets

From the synthetic radio (43 GHz) and optical (5 × 10¹⁴ Hz) surface-brightness maps obtained as described in Sect. 2.3, we extracted flux-weighted brightness positions, “centroid” positions, and computed their projected angular separation, hereafter the radio-optical offset d_app, following FdC25. Throughout this work, we adopt z = 1, corresponding to the median redshift of the Plavin et al. (2019) sample and retain only cases where the optical emission exceeds the Gaia detection threshold.

These maps also allowed us to quantify the apparent jet opening angle ${\theta }_{\mathrm{j}}^{\mathrm{app}}$ and to investigate how the spatial morphology of the emission varies with jet power, viewing angle, and stellar distribution. Consistent with the dynamical limitations discussed in the previous section, synthetic emission maps were computed over the full simulation domain, but radio and optical centroid positions (and derived quantities) were only considered physically meaningful when the corresponding centroid lies upstream of z(v = 0.9c). The axial position z(v = 0.9c) was converted into an apparent projected distance from the jet base using the viewing angle θ and the angular-diameter distance D_A = D_L/(1 + z)², which relates physical sizes to observed angles at cosmological distances. This conversion is purely geometrical: relativistic effects such as Doppler boosting modify the observed emissivity along the jet but do not shift the apparent position of the kinematic transition in the steady-state approximation adopted here. Results related to radio and optical centroids situated downstream z(v = 0.9c) are reported anyway and are only used for qualitative interpretations.

4.1. Viewing angle

The top and central panels of Fig. 4 show the positions of the radio and optical centroids (top) and the resulting offset d_app (central) as a function of jet power, for a viewing angle of 5° (left), and two stellar core radii, r_c, s = 0.5 and 1 kpc. At low powers, the jets undergo strong deceleration, and the radio and optical centroids are nearly coincident (small offsets; central panel). With increasing jet power, the separation between the two peaks grows, with the optical peak appearing ahead of the radio one (positive offsets). This behavior is common to all stellar distributions. At larger jet powers, the radio-optical offset reaches a similar maximum value for both stellar distributions (d_app ∼ 1 mas) but is shifted in jet power (∼7 × 10⁴² and ∼2 × 10⁴³ erg s⁻¹, for r_c, s = 0.5 and 1 kpc, respectively). This maximum value remains above the Gaia sensitivity, showing a window where the radio-optical core shift can be more easily detectable for jet powers around 10⁴³ erg s⁻¹, at redshift z = 1.

Fig. 4. Left: Radio and optical centroid positions (zradio and zopt; top), resulting radio-optical offset dapp (center, where the error bars correspond to the 1σ uncertainty on the centroid positions, estimated from the effective PSF widths), and apparent jet opening angle ${\theta }_{\mathrm{j}}^{\mathrm{app}}$ derived from the 43 GHz emission (bottom) shown as a function of jet power Lj for two stellar core radii rc, s (see legend). The dashed line in the central panel represents the Gaia astrometrical precision. These results are obtained for a viewing angle θobs = 5°. The gray symbols indicate cases for which the corresponding centroid lies downstream of z(v = 0.9c) and are therefore excluded from the analysis. Right: Same as left but for a viewing angle θobs = 30°.

Beyond the maxima, the radio–optical offset first decreases and then remains close to the minimum value (∼0.5 mas) for r_c, s = 0.5 kpc, whereas for r_c, s = 1 kpc it increases for L_j ≳ 6 × 10⁴³ erg s⁻¹. This growth in d_app reflects the increasing extent of the dissipation region in powerful jets, which shifts the optical peak downstream, while the radio peak remains unchanged because it is set by opacity. In principle, this behavior would make such shifts detectable even in powerful sources. However, the energy dissipation in these jets is too low, causing the optical emission to fall below the sensitivity of Gaia (L_j ≳ 4 × 10⁴⁴ erg s⁻¹, for r_c, s = 0.5 kpc; L_j ≳ 8 × 10⁴⁴ erg s⁻¹, for r_c, s = 1 kpc).

Figure 4 (right) shows the same quantities for a larger viewing angle of 30°. The top panel shows the location of the radio and optical centroids. Doppler (de-)boosting sets the maximum jet power at which the optical centroid remains visible. For instance, it is lost around 3 × 10⁴³ erg s⁻¹ for r_c, s = 0.5 kpc and at ∼10⁴⁴ erg s⁻¹ for r_c, s = 1 kpc. In this figure, the gray symbols indicate the cases in which the radio and/or optical centroid positions are downstream of the position where the velocity becomes smaller than 0.9c, which occurs for low-power jets. Despite this limitation, we expect similar trends as the ones observed for a smaller viewing angle (see Fig. 4, left): a plateau of the radio-optical offset due to strong dissipation closer to the jet base, before it falls to zero for the weakest jets.

The radio emission maps were used to estimate the apparent jet opening angle ${\theta }_{\mathrm{j}}^{\mathrm{app}}$ directly from the projected surface-brightness distribution. See Pushkarev et al. (2017) and Kovalev et al. (2020a) for details on opening angle measurements from VLBI images. Assuming a conical morphology, we measured the transverse jet width at multiple axial positions within the relativistic region (where v > 0.9c), using a fixed fractional intensity threshold of the peak brightness and derived ${\theta }_{\mathrm{j}}^{\mathrm{app}}=2arctan[(w/2)/z]$ (w being the local jet width, and z the associated axial position) from the maximum projected half-opening angle. This quantity represents a purely image-based, projected opening angle and therefore depends both on intrinsic jet geometry and on the viewing angle. The bottom panels of Fig. 4 show the resulting ${\theta }_{\mathrm{j}}^{\mathrm{app}}$ values derived from 43 GHz images at viewing angles of 5° (left) and 30° (right), respectively, as a function of jet power.

In the first case, we see how mass-load by stellar winds can contribute to decollimation of jets with powers below L_j ∼ 2 × 10⁴³ erg s⁻¹, while not significantly affecting more powerful ones. The resulting apparent opening angles reach relatively high values at a viewing angle of 5° (up to 50°) for low-power jets before a sharp decrease around 10⁴³ erg s⁻¹, while jets remain quasi-cylindrical for larger powers. At viewing angles of 30° a similar behavior is observed, although the values given for jets with larger opening angles are based on estimates at z(v = 0.9c).

4.2. Stellar distribution

To further examine the effect of the stellar core radius on the radio-optical offset and jet collimation, we show the centroid positions (top panel) and corresponding offset d_app (middle panel) as a function of r_c, s in Fig. 5, for a viewing angle of 5° (left), and for a viewing angle of 30° (right). The plots display the results for two jet powers: L_j = 10⁴³ and 10⁴⁴ erg s⁻¹.

Fig. 5. Left: Radio/optical centroid position (zradio/zopt, top), the resulting radio-optical offset dapp (center, where the error bars correspond to the 1σ uncertainty on the centroid positions, estimated from the effective PSF widths) for two jet powers Lj (see legend), and the apparent opening angle from the radio emission (bottom) as a function of the stellar radius rc, s and obtained with θobs = 5°. The dashed line in the central panel represents the Gaia astrometrical precision. If any, the gray points represent excluded radio/optical centroids, as their positions are downstream of the region where v > 0.9c. Right: Same as left but with θobs = 30°.

At 5° (Fig. 5, left column), only radio-optical offsets for r_c, s ≳ 0.2 − 0.4 kpc (depending on the jet power) are resolvable by Gaia – i.e., they appear over the Gaia angular resolution indicated by a dashed line in the plot. In contrast, at 30° (Fig. 5, right), all offsets are resolvable, but not all are detectable or lie within the region where v > 0.9c. For instance, for a 10⁴⁴ erg s⁻¹ jet, offsets can only be detected when r_c, s > 1 kpc. At a jet power of 10⁴³ erg s⁻¹, and r_c, s > 0.4 kpc, radio and/or optical centroids are situated downstream of our validity region and are therefore not considered. Nonetheless, we expect a similar trend to the one derived at θ_obs = 5°, i.e., an increase in the offset before reaching a local maximum followed by a plateau or a smooth decline.

The dependence of the apparent jet opening angles with r_c, s is shown in the bottom panels of Fig. 5. At low jet powers, the opening angle increases with r_c, s, reflecting stronger dissipation, in contrast to the case of more powerful jets, which remain basically cylindrical.

4.3. Summary of the results

In summary, there is a clear entanglement between the roles of the jet power, stellar distribution, and jet orientation on the values obtained for the radio-optical offsets. For small observation angles, in general, smaller jet powers and enhanced mass-load (via larger stellar distribution cores) are required for the offset to be observable. For larger viewing angles, offsets suffer from detectability rather than resolvability (see also FdC25). We observe that an extended stellar distribution core could compensate for an, in principle, limited degree of dissipation in the case of powerful jets. The radio-morphological transition from large to low opening angle jets occurs at comparable powers of about L_j ∼ 10⁴³ erg s⁻¹ and is reflected in the jet radii obtained at 2 kpc. Although we do not included all the possible dissipative mechanisms in our limited quasi-1D simulations, this result gives us the correct order of magnitude of jet powers at which large-scale morphologies transit from FR II to FR I. The inclusion of extra dissipative processes would probably enhance the differences between the two regimes found here, but we do not expect that this would change the jet power interval that we find for the transition.

Our results indicate that, at cosmological distances (e.g., z ∼ 1), the optical jet emission is unlikely to remain sufficiently bright relative to the nucleus to produce measurable radio-optical centroid offsets once the jet reaches the strongly dissipative regime characteristic of FR I jets for θ_obs > 5° (Laing & Bridle 2014). This does not contradict the presence of prominent optical jets in nearby sources such as M 87, where spatially resolved observations reveal extended emission powered by additional large-scale dissipation mechanisms (e.g., standing shocks such as HST-1) that are not included in our present modeling.

5. Discussion

Under the assumption that a fraction of the dissipated energy is invested into particle acceleration, our jet simulations show that detectable milliarcsecond radio-optical separations arise within a restricted range of jet kinetic powers, namely 10⁴² ≲ L_j ≲ 10⁴⁴ erg s⁻¹. At low jet powers or large viewing angles, a significant fraction of the simulated jets lie beyond the axial position z(v = 0.9c), where our quasi-1D treatment ceases to be reliable (see Fig. 3), thereby limiting the range of offsets that can be robustly interpreted. In such cases, the trends identified in this work should be regarded as indicative of the underlying physical behavior rather than as precise quantitative predictions.

Combined with the offset distributions modeled in FdC25 and considering only the impact of stellar mass loading within the jet, we can draw two main conclusions: (i) dissipation induced by stellar mass loading constitutes a viable physical mechanism for particle reacceleration and for producing positive radio-optical offsets; and (ii) variations in jet power, stellar distribution properties, and viewing angle systematically shift the location of the centroids along the jet and modulate both its amplitude and detectability. The evolution of these parameters could imprint redshift-dependent trends on the observed radio-optical offset population, making this phenomenon a sensitive probe of jet-galaxy interactions across cosmic time.

We stress that the specific jet-power ranges discussed in terms of core offset detectability are not universal but depend on the adopted mass-loading normalization Q₀ and on redshift. Varying either parameter primarily results in a systematic shift of the same qualitative behaviors toward higher or lower jet powers, rather than in a change of the underlying trends. In this sense, these results should be interpreted as a representative slice of a broader parameter space. In FdC25, which explicitly explored the dependence of radio-optical offsets on Q₀, redshift, and viewing angle, we provided a framework to extrapolate our findings beyond the specific choices adopted here and to connect them across cosmic time.

Finally, although we have established a common initial jet Lorentz factor for all our models, we do not expect significant differences in our results. Anglés-Castillo et al. (2021) showed that the deceleration (and dissipation) distance was very similar, regardless of the injection Lorentz factor, within the expected values in this region, i.e., 5 − 10 (see, e.g., Fig. 5 in that paper).

5.1. Implications for the FRI/FR II dichotomy

In Sect. 4 we saw a clear transition from basically cylindrical to large opening angle jets around L_j ∼ 10⁴³ erg s⁻¹. This probably reveals the jet powers at which direct mass-load by stellar winds may play a significant role in jet deceleration (see also Perucho et al. 2014; Anglés-Castillo et al. 2021).

If we assume that canonical FR I jets studied in, for instance, Laing & Bridle (2014) have jet powers around 10⁴³ erg s⁻¹, we see that the apparent opening angles inferred at a viewing angle of 30° (see, e.g., Fig. 4, right) are systematically smaller than those measured. However, we derive similar opening angles to those reported in the aforementioned paper for smaller jet powers.

This difference is expected for several reasons. First, our estimates are restricted to milliarcsecond scales and, by construction, only include regions of the flow that remain dynamically valid within our modeling framework, i.e., upstream of the axial position z(v = 0.9c). Furthermore, we see that for a viewing angle of 30° (more typical of radio-galaxies), a nonnegligible fraction of the optical peaks are located beyond z(v = 0.9c). As a result, the opening angles reported here characterize the still-collimated, mildly decelerated jet rather than the fully developed flaring region where velocities are reported to be significantly smaller (≃0.3 − 0.4 c Laing & Bridle 2014).

Second, our synthetic maps probe the jet emission at 43 GHz, which fades rapidly with distance due to adiabatic expansion losses, in contrast with the 5 − 8 GHz VLA observations used by Laing & Bridle (2014). At arc-second scales, the combined effects of an ambient pressure profile, jet deceleration, and enhanced dissipation lead to nonadiabatic evolution and substantial lateral expansion, producing the large opening angles characteristic of FR I flaring regions. These processes are largely associated with the entrainment of external material at the jet boundaries, driven by mechanisms such as small-scale instabilities or the direct interaction of stars with the jet surface (e.g., Matsumoto et al. 2017; Gourgouliatos & Komissarov 2018; Perucho 2020). Such boundary-driven entrainment and large-scale geometric flaring cannot be captured within the quasi-1D, steady-state framework adopted here, which is instead designed to follow the internal evolution of the jet spine under distributed mass loading. Accurately modeling the transition to the flaring region would therefore require fully 2- or 3D simulations, which we defer to future work.

5.2. Implications on the jet power and stellar properties

Our simulations show that each parameter of the model has a significant impact on the radio–optical centroid positions, as well as on the amplitude (d_app) and detectability of the resulting offsets. In particular, increasing the jet kinetic power L_j shifts both radio and optical centroids downstream, increases d_app, but reduces detectability due to less efficient dissipation. It should be noted that these trends correspond to a controlled exploration of the parameter space in which individual quantities are varied independently. In reality, the model parameters are intrinsically coupled, and the resulting observables reflect the combined effects of jet power, mass loading, and viewing geometry. As a consequence, these trends may develop more complex, non-monotonic behaviors when multiple parameters vary simultaneously, as shown in Figs. 4 and 5. In particular, the complex evolution of the radio–optical offsets with the jet power can be understood as a direct manifestation of these coupled dependencies, arising from transitions between regimes dominated by dissipation, jet inertia, and radiative or centroid weighting effects.

At the redshifts considered, the main stellar contribution to mass-load could be the thermally pulsing asymptotic giant branch (TP-AGB) stars, which are characterized by average mass-loss rates between Ṁ ∼ 10⁻⁹−10⁻⁵ M_⊙ yr⁻¹ and low wind velocities, ∼10 km s⁻¹ (Leitner & Kravtsov 2011; Höfner & Olofsson 2018; O’Shea et al. 2025). Their average mass is expected to be, on average, on the order of several solar masses (Herwig 2005). If we consider a conservative AGB stellar density of n_AGB ∼ 10⁻³ pc⁻³ (O’Shea et al. 2025), we derive Q₀ values between 2 × 10²¹ to 2 × 10²⁵ g yr⁻¹ pc⁻³, which includes the value used in this work and in FdC25.

Concerning the radial extent of the core of the gas and stellar distributions, r_c, s, it gives an idea of the characteristic interaction/deceleration scale. This value can be estimated using fits to Nuker profiles (see, e.g., Lauer et al. 2007), and it should be correlated with the characteristic deceleration distance in low-power jets. In fact, at these redshifts, early-type galaxies have a typical stellar mass ∼10¹¹ M_⊙ and a median effective radius (the radius enclosing most of the stellar mass) between 2 and 3 kpc (van der Wel et al. 2014), which is consistent with our assumptions and the comparison made between our results and Plavin et al. (2019).

In summary, our proposed model to explain the radio-optical offsets detected in the quasar population relies on the interaction between jets and the bulge stellar component of the host. Although quasars might be fed by powerful jets, if they propagate through a relatively extended and dense stellar distribution, i.e., maximizing both r_c, s and Q₀, an intense process of dissipation within that region, together with a small viewing angle, can bring the optical centroids close to the base of the radio jet. In this case, the effect of cosmological distance on the detectability is counterbalanced by the Doppler boosting, which matches with the fact that the sample in Plavin et al. (2019) naturally selects Doppler-boosted objects.

5.3. Local Universe (z < 1)

The predominance of quasars in radio-optical offset samples at z ∼ 1 could simply reflect the cosmic evolution of black-hole accretion, with luminous, high-Eddington AGN peaks populating intermediate redshifts (Shankar 2009), in contrast to low-accretion systems such as BL Lac objects and Seyfert galaxies dominating at late cosmic times (Heckman & Best 2014). Neither the angular resolution nor the received fluxes are, in principle, limiting factors for offset detection at low redshifts.

That being said, massive ellipticals at z < 1 are typically dominated by old stellar populations (Thomas et al. 2005). In such systems, stellar mass return is primarily driven by evolved low- and intermediate-mass stars in the red giant and TP-AGB phases. Although TP-AGB stars dominate the mass-load budget, their mass-loss rates decline significantly with age (O’Shea et al. 2025). Consequently, old elliptical galaxies would provide lower cumulative mass-loading into their jets. Regarding jets, Seyferts are also known to host low-power jets (typically < 10⁴⁴ erg s⁻¹) with relatively large viewing angles (θ_obs ≤ 30° for Seyfert II type), whereas BL Lacs show higher intrinsic kinetic power jets (typically between 10⁴³ to 10⁴⁵ erg s⁻¹) with smaller viewing angles (θ_obs ≤ 5°, see, e.g., Urry & Padovani 1995; Heckman & Best 2014; Ghisellini et al. 2014).

We can discuss the significant radio-optical offsets found in these two types of galaxies in the frame of our results. Seyfert galaxies frequently contain circumnuclear molecular gas reservoirs, as revealed, for example, by ALMA observations in NGC 1052 (Medeiros et al. 2020). This gas can either fuel new intermediate-age star formation or be directly entrained by the jet through cloud shredding, effectively increasing the mass-loading rate Q₀. In the context of our simulations, even a modest increase in Q₀ can move a source into the detectable offset regime. In Seyfert I galaxies, where the accretion disk is directly visible, cases in which mass loading is insufficient to trigger efficient energy dissipation will naturally result in disk-dominated optical emission, producing negative radio-optical separations, consistent with the upstream offsets reported by Plavin et al. (2019). This interpretation is further supported by Plavin (2026), who showed that optical flares in blazars systematically induce Gaia centroid motions: during flux rises, the optical centroid shifts upstream, while during flux decay, it drifts downstream. They localize the flaring region to within < 1 mas of the VLBI core, corresponding to a few parsecs in projection. In the context of our framework, this behavior naturally complements the negative radio-optical separations expected in Seyfert I systems when the accretion disk dominates, while also suggesting that compact jet flares near the radio core may temporarily reduce or even reverse the downstream offsets produced by distributed dissipation. In contrast, in Seyfert II galaxies, the accretion disk is obscured, and any detectable radio-optical offset must arise from downstream jet-related emission, leading exclusively to positive separations.

In the case of BL Lacs, Doppler boosting substantially enhances the optical synchrotron emission produced along the jet, making offsets detectable even when stellar winds are weaker. In contrast, our model predicts that higher-power objects such as FR II galaxies would rarely produce large downstream offsets, as expected from the irrelevance of dissipation in their case. For instance, Cygnus A (with L_j > 10⁴⁵ erg s⁻¹ Snios et al. 2018) displays no detectable VLBI-Gaia displacement. Such jets remain relativistic and emit only weak extended optical radiation outside the core.

5.4. Early Universe (z > 2)

Above redshifts z > 2, the incidence of radio-optical offsets rapidly decreases (Plavin et al. 2019). Galaxies may not have yet developed significant TP-AGB populations at these early epochs (Bruzual & Charlot 2003). Instead, mass-load could be dominated by winds of Wolf-Rayet (WR) stars, with terminal velocities of v ∼ 10³ km s⁻¹ (Mengel & Tacconi-Garman 2007) and mass-loss rates as high as Ṁ ∼ 10⁻⁴ M_⊙ yr⁻¹ (Barlow et al. 1981). However, WR stars are intrinsically rare (only a few thousand are expected in the whole Milky Way) and are preferentially found in galaxies with intense, localized star formation (Liang et al. 2020). Recent James Webb Telescope observations have confirmed their presence in galaxies at z ∼ 2.2 (Curti et al. 2026) and even up to z ∼ 6.1 (Berg et al. 2026), but these detections remain associated with compact star-forming regions rather than widespread stellar populations. Consequently, their volumetric number density is expected to be low, implying that WR-jet interactions could only sporadically occur and over limited spatial scales. Such interactions are governed by large-scale, strong shocks rather than gradual mixing from the distributed stellar population (e.g., Hubbard & Blackman 2006). The resulting coupling therefore leads primarily to dramatic, localized perturbations, while the efficiency of sustained volumetric mass loading into the bulk relativistic jet is reduced compared to red giants and/or TP-AGB winds (e.g., Wykes et al. 2015). Supernova explosions within the jet represent an even more extreme case: the interaction triggers reverse shocks and hydrodynamical instabilities that drive turbulent mixing of the ejecta with the jet material, temporarily decelerating the flow (Longo et al. 2025). At the same time, radio AGNs could be intrinsically more powerful at these epochs, with typical jet kinetic powers exceeding 10⁴⁵ erg s⁻¹ (Clark et al. 2014), which would limit the dissipation induced by mass loading.

Combined with a strong observational bias toward bright, disk-dominated quasars in Gaia samples, these physical conditions suppress significant offsets or favor negative values at high redshifts. Specifically, the optical jet emission typically falls below the Gaia sensitivity, causing the measured optical centroid to be dominated either by the inner jet (suppressing offsets) or by the accretion disk (resulting in a negative offset), rather than by synchrotron emission from the extended jet.

6. Conclusions

We investigated the physical origin and detectability of parsec-scale radio-optical offsets in AGN jets using axisymmetric RMHD simulations with continuous stellar mass loading, post-processed with the RIPTIDE radiative transfer code. We assumed that a fixed fraction of the dissipated energy, ϵ_e ∼ 0.1, is transferred to the acceleration of nonthermal particles. Our results indicate that detectable VLBI-Gaia separations of typically ∼0.1 − 4 mas (corresponding to projected distances of a few to a few tens of parsecs at z = 1) arise only within a restricted region of jet-environment parameter space, governed by the interplay between jet kinetic power, stellar mass-loading efficiency, and viewing geometry. Importantly, radio-optical offsets are flux-weighted astrometric signatures rather than purely geometric tracers of dissipation and therefore depend on the relative optical contributions of the jet, accretion disk, and host galaxy. The offsets reported in this work were derived by restricting the analysis to centroid positions located upstream of the axial position z(v = 0.9c), beyond which our quasi-1D treatment ceases to be reliable. As a consequence, the trends identified here characterize the inner, mildly decelerated jets, and the estimates of centroid positions should be regarded as lower limits when dissipation extends farther downstream, particularly at large viewing angles or for extended stellar distributions. Nevertheless, we find these trends to be physically meaningful and useful to derive clear conclusions from our work.

For the stellar environments explored here, intermediate-power jets (L_j ∼ 10^42.5 − 10⁴⁴ erg s⁻¹) embedded in hosts with substantial TP-AGB wind densities (Q₀ ∼ 5 × 10²⁴ g yr⁻¹ pc⁻³) produce the largest positive (downstream) radio-optical offsets, with apparent separations of ∼0.1 − 4 mas. Below this power range, dissipation occurs too close to the jet base for the radio and optical centroids to separate significantly, while above it the flow remains largely ballistic and the optical synchrotron emissivity is insufficient to displace the centroid unless aided by unusually dense stellar environments or opacity effects.

The jet-power window identified here as the plausible regime for offset detection is not universal but scales with the adopted stellar mass-loading normalization. As demonstrated in FdC25, variations in Q₀ or redshift shift this window without altering the qualitative trends identified in this work, and the semi-analytical framework of that study provides a natural basis for extrapolating our results across a broader parameter space and cosmic time. Viewing geometry further modulates offset detectability through Doppler boosting and projection effects. This naturally explains why flat-spectrum quasars, BL Lac objects, and partially aligned FR I or Seyfert jets dominate VLBI-Gaia offset samples (Plavin et al. 2019).

Finally, while stellar mass loading is unlikely to be the sole dissipation mechanism in relativistic jets, it constitutes an unavoidable and persistent source of jet perturbation in galactic nuclei. Adding to the results of previous studies (e.g. Perucho et al. 2014; Anglés-Castillo et al. 2021), we also conclude here that significant jet deceleration produced solely by stellar winds could only be possible for jets with powers L_j ≤ 10⁴³ erg s⁻¹. However, even when subdominant, the effect of stellar winds in jets sets a baseline level of dissipation that can be amplified by additional processes such as the development of instabilities or recollimation shocks.

Future improvements in both optical and radio astrometry will directly test the parameter regime identified in this work. In particular, astrometric precision at the tens-of-microarcsecond level will be crucial for distinguishing small offsets from noise. Moreover, improved multifrequency VLBI imaging will help disentangle compact core positions from extended jet structure by reducing structure-dependent systematic errors. The planned next-generation Very Large Array (ngVLA, Selina et al. 2018), with an order-of-magnitude increase in sensitivity and resolution relative to current centimeter-millimeter facilities (1.2 − 116 GHz) and long baselines for high-fidelity imaging, will enable a more accurate localization of compact jet cores as well as detailed mapping of jet morphology, thereby refining offset measurements and tests of the predicted dependence on jet power and viewing geometry.

Acknowledgments

We thank the anonymous reviewer for the relevant comments and suggestions that helped us improve the manuscript. We also thank Sasha Plavin and Elena Shablovinskaya for valuable comments on the manuscript. This work has been supported by Banco Santander and Universitat de València (International Mobility Grant 2025). This work has been supported by the Spanish Ministry of Science through Grant PID2022-136828NB-C43, and from the Generalitat Valenciana through grant CIPROM/2022/49. YYK was supported by the MuSES project, which has received funding from the European Union (ERC grant agreement No 101142396). Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or ERCEA. Neither the European Union nor the granting authority can be held responsible for them.

References

Appendix A: Numerical resolution convergence

To verify that our results are not significantly affected by the spatial resolution of the quasi-one-dimensional RMHD grid, we performed a resolution convergence study on a representative simulation in which a clear radio-optical offset is present. We selected one representative run (L_j ≃ 10⁴³ erg s⁻¹, r_c, s = 500 pc), which lies well within the power range where d_app is detectable and maximal (see Fig. 4). The fiducial grid for this run has 800 × 5000 cells over 20 × 2000 pc (r × z). We reran the simulation at resolution multipliers ranging from 0.1× to 2× the fiducial, corresponding to grids from 80 × 500 to 1600 × 10 000 cells, and recomputed the synthetic emission maps and centroid positions with RIPTIDE at each resolution.

The results are shown in Fig. A.1. At very low resolutions (≲0.25×), the grid is too coarse to resolve the jet structure and produces no detectable synchrotron emission (gray shaded region). At intermediate resolutions (0.3 × −0.4×), the optical flux falls below the Gaia threshold (open symbols), yielding unreliable centroid estimates. From 0.5× onward, both the radio and optical centroids stabilize rapidly (panel a), and d_app converges to the fiducial value of ≃1.1 mas (panel b, dashed line). The 2× run reproduces this value to within ≲2%, confirming that the fiducial resolution is adequate. We therefore conclude that our results are numerically converged and that the radio-optical offsets reported in this work are not resolution artifacts.

Fig. A.1. Numerical resolution convergence from one representative run (Lj ≃ 1043 erg s−1, rc, s = 500 pc, θobs = 5°). Top: Radio (zradio, blue) and optical (zopt, red) centroid positions as a function of the resolution multiplier relative to the fiducial 800 × 5000 grid. The filled circles pass the Gaia flux threshold; the open circles do not. The light-gray band marks the resolutions at which no synchrotron emission is produced. Bottom: Corresponding dapp values (dark-green circles). The dashed line indicates the fiducial (1×) reference.

All Tables

All Figures

Fig. 1. Schematic of the jet geometry. The approaching jet propagates from the nucleus at angle θobs to the line of sight. The radio and optical centroids (zradio, zopt) are projected onto the sky plane, defining the apparent offset dapp into the sky plane from the observer’s point of view. The solid line in the jet marks z(v = 0.9c); the stellar distribution (rc, s) is shown in orange (not drawn to scale).

In the text

	Fig. 2. Synthetic synchrotron surface-brightness maps at 43 GHz (left) and optical 5 × 1014 Hz (right) for a representative run with Lj ≃ 1043 erg s−1, rc, s = 500 pc, θobs = 30°, at redshift z = 1. The dashed lines mark the intensity-squared centroid positions zradio (blue) and zopt (red); the bracket to the right of the optical panel indicates the resulting dapp = 6.17 mas.
In the text

Fig. 3. Left: Jet radius Rj at z = 2 kpc (top), average energy per lepton $\overline{ϵ}{\prime }_{\mathrm{e}}$ (second), average proton fraction within the jet (third), and axial deceleration distance z(v = 0.9c) (bottom) as a function of the jet power Lj, for two stellar core radii rc, s (see legend). Each symbol corresponds to one simulation. The circles indicate simulations where the jet remains relativistic throughout the full 2 kpc box; the triangles indicate simulations where the jet decelerates before 2 kpc, for which the jet radius is obtained by extrapolation (see Sect. 3). The dashed line in the bottom panel marks the 2 kpc box limit. Right: Same as left but shown as a function of the stellar distribution core radius rc, s for fixed jet powers (see legend).

In the text

Fig. 4. Left: Radio and optical centroid positions (zradio and zopt; top), resulting radio-optical offset dapp (center, where the error bars correspond to the 1σ uncertainty on the centroid positions, estimated from the effective PSF widths), and apparent jet opening angle ${\theta }_{\mathrm{j}}^{\mathrm{app}}$ derived from the 43 GHz emission (bottom) shown as a function of jet power Lj for two stellar core radii rc, s (see legend). The dashed line in the central panel represents the Gaia astrometrical precision. These results are obtained for a viewing angle θobs = 5°. The gray symbols indicate cases for which the corresponding centroid lies downstream of z(v = 0.9c) and are therefore excluded from the analysis. Right: Same as left but for a viewing angle θobs = 30°.

In the text

Fig. 5. Left: Radio/optical centroid position (zradio/zopt, top), the resulting radio-optical offset dapp (center, where the error bars correspond to the 1σ uncertainty on the centroid positions, estimated from the effective PSF widths) for two jet powers Lj (see legend), and the apparent opening angle from the radio emission (bottom) as a function of the stellar radius rc, s and obtained with θobs = 5°. The dashed line in the central panel represents the Gaia astrometrical precision. If any, the gray points represent excluded radio/optical centroids, as their positions are downstream of the region where v > 0.9c. Right: Same as left but with θobs = 30°.

In the text

Fig. A.1. Numerical resolution convergence from one representative run (Lj ≃ 1043 erg s−1, rc, s = 500 pc, θobs = 5°). Top: Radio (zradio, blue) and optical (zopt, red) centroid positions as a function of the resolution multiplier relative to the fiducial 800 × 5000 grid. The filled circles pass the Gaia flux threshold; the open circles do not. The light-gray band marks the resolutions at which no synchrotron emission is produced. Bottom: Corresponding dapp values (dark-green circles). The dashed line indicates the fiducial (1×) reference.

In the text