© The Authors 2025

1. Introduction

Extragalactic jets emerging from actively accreting supermassive black holes (SMBHs) located in active galactic nuclei (AGNs) have been traced through their radio synchrotron emissions since the 1950s (Jennison & Das Gupta 1953; Hazard et al. 1963). Currently, such AGN outflows are being widely investigated through multi-wavelength and multi-messenger approaches, with the origins of such jets linked to the interplay between magnetic fields and rotation, either of the SMBH or of the accretion disk, through which matter is funneled into the SMBH (e.g., Blandford et al. 2019). The multi-wavelength properties of jets are now generally understood to be primarily governed by non-thermal processes such as synchrotron and inverse-Compton mechanisms (e.g., Hardcastle & Croston 2020). However, these observed emission processes are greatly influenced by a complex set of intrinsic microphysical processes (e.g., particle re-acceleration, matter entrainment, turbulence), many of which are still highly debated in astrophysical flows (e.g., Pe’er 2014).

With the recent advancements in radio observations using high-resolution and sensitive telescopes, there is now a promising expansion in exploring the lower-energy hump of the jets’ emission spectra dominated by the synchrotron mechanism (e.g., Saikia 2022; Mahatma 2023). The observations have not only detected new features associated with jet evolutions (Ramatsoku et al. 2020; Brienza et al. 2021; Brüggen et al. 2021; Knowles et al. 2022), but have also revealed detailed insights into the micro-scale processes governing such emission characteristics (Turner et al. 2018; Mahatma et al. 2020; Condon et al. 2021). In parallel with these investigations, ongoing advancements in numerical techniques (e.g., the Lagrangian-Eulerian particle-fluid technique) have paved the way for testing these intricate processes in high-resolution simulation grids, thereby enhancing our understanding of the impact of modeling such phenomena (Vaidya et al. 2018; Winner et al. 2019; Ogrodnik et al. 2021; Dubey et al. 2023, 2024).

With this progress, longstanding inquiries regarding radio galaxies have been revitalized, particularly the question of how far AGN jets can travel after their generation at the launching sites. Recent discoveries of distinct extragalactic jets extending nearly 5 and 7 Mpc in linear sizes necessitate a closer examination of their evolution mechanisms (Machalski et al. 2008; Oei et al. 2022, 2024a), especially when compared to the existence of AGN jets that appear to saturate in their evolution on galactic scales (O’Dea & Saikia 2021; Baldi 2023). The question is particularly relevant given the energy-dissipating effects of processes such as radiative and adiabatic cooling, along with the entrainment of ambient matter, on the radiating non-thermal particles of these sources across different redshifts (O’Dea & Baum 1997; Schoenmakers et al. 2000; Rossi et al. 2008). To further complicate the context, the detection of persistent emission in the radio lobes of distinct “giant” radio galaxies, which have ceased their active phase by halting jet ejection (Subrahmanyan et al. 2008; Bagchi et al. 2014; Cotton et al. 2020; Bagchi et al. 2024), calls for studies that explore the correlation between the dynamical evolution of these sources and their conversion into radiative fingerprints.

Giant radio galaxies (GRGs) are a minor subclass of the radio galaxy population characterized by their extensive linear sizes, reaching 700 kpc or more with the current Hubble constant value, H₀ = 70 km s⁻¹ Mpc⁻¹ (some studies suggest a revised threshold of 1 Mpc, e.g., Dabhade et al. 2023). Alongside efforts to link morphological traits with hypothesized formation mechanisms (Subrahmanyan et al. 1996; Malarecki et al. 2015; Cotton et al. 2020), current research on GRGs aims to reconstruct their intrinsic three-dimensional dynamical properties from their projected emission maps and assess the feasibility of the observed states within the framework of established jet evolution theories (Oei et al. 2023). For instance, properties such as the source’s spectral behavior and the magnetic-field estimates provide insights into the evolutionary timescales of jet-cocoon structures (Machalski et al. 2009; Machalski 2011). Meanwhile, the fractional polarization (FP) values and the magnetic-field distribution reveal zones of internal and external stress layers, highlighting regions of shear and compression within and in the vicinity of the cocoon (Laing 1981; Mack et al. 1997; Cotton et al. 2020). Studies on Faraday rotation measure emphasize the value of using GRGs as probes to estimate ambient-particle density and the correlation length of small-scale environmental turbulence (Cantwell et al. 2020; Stuardi et al. 2020). A key finding, as illustrated in Stuardi et al. (2020), suggests that larger GRGs are more likely to be detected in polarized emission, indicating that their extensive lobes extend beyond the dense central environment and evolve within a rarified medium in the peripheral regions. This observation is also consistent with several continuum studies of GRGs (Subrahmanyan et al. 2008; Oei et al. 2024b), highlighting their prevalence within the relatively low-density warm-hot intergalactic medium (WHIM; Davé et al. 2001), typically situated beyond the boundaries of galaxy groups or clusters.

Despite extensive numerical efforts to simulate AGN jets including the effects of microphysical processes such as particle re-acceleration and cooling (Borse et al. 2021; Mukherjee et al. 2021; Kundu et al. 2022; Dubey et al. 2023), no such simulations to date have specifically addressed GRGs. This gap is largely due to the substantial computational demands required to resolve small-scale physics (∼100 s of pc) while allowing jets to evolve over extended scales of a few megaparsec. In fact, numerical modeling to test the genesis hypotheses of GRGs (i.e., recognizing the dominant factors contributing to their growth: (a) longer temporal evolution, (b) efficient, collimated jets, (c) low-density environment to jet propagation, or (d) rapidly propagating restarting jets) has been notably scarce until recently. A 2D axisymmetric model has been introduced recently to explore the evolution of symmetric jets within a dense cluster environment, offering a preliminary framework for studying a subset of GRGs (Duan et al. 2025). However, it is now well established that the imposed symmetry and reduced dimensionality in these models often lead to discrepancies in results when compared to more accurate 3D simulations (Mignone et al. 2010; Bodo et al. 2016; Massaglia et al. 2016; Kundu et al. 2022).

Recognizing the necessity of investigating the influence of factors in producing extended GRGs (≳1 Mpc) within a full 3D framework, only recently has a set of five key models been developed (Giri et al. 2025), spanning a diverse range of jet-environment configurations in relativistic magnetohydrodynamic (RMHD) settings. These models employ symmetry-breaking mechanisms intrinsic to three-dimensional configurations, designed to more closely mimic realistic jet-medium interactions. The resulting dynamical features and their evolutionary implications are documented in Giri et al. (2025) (hereafter referred to as Paper I). While this study marks the first effort to develop 3D RMHD models for GRGs, the resolution achieved was necessarily limited. The entire jet injection region is defined over a cylindrical volume of ten cells, with the jet diameter corresponding to two compute cells – a compromise to balance physical accuracy and computational efficiency at this stage. Nonetheless, efforts were made to minimize the impact of this constraint, particularly by leveraging the natural expansion of the jet following its injection. A detailed discussion of these mitigation strategies and the effects of resolution enhancement can be found in Paper I.

Given the significance of transforming 3D dynamical data into insightful 2D projected emission maps, as discussed earlier, we intend to spearhead the radiative modeling of simulated GRGs. This study aims to comprehensively examine the correlation between the intrinsic dynamical behavior of (giant) jetted sources and their continuum and polarization properties. Additionally, it explores the underlying physical mechanisms that contribute to the suppression of the radiative signatures from such sources.

This paper is structured as follows: In Section 2, we provide a concise overview of the relevant details and conclusions from Paper I that form the foundation of this study. Section 3 presents the transformation of dynamic maps from Paper I into their corresponding radiative signatures at 1 GHz, which are examined in detail in this paper. This is followed by a discussion in Section 4 on the matter-transport processes responsible for generating the GRG cocoon. In Section 5, we explore the emission characteristics of the simulated GRGs, providing several quantitative insights into their properties. In Section 6, we examine the polarization characteristics of the simulated structures, linking them to their underlying dynamical properties. We present a summary of our key findings in Section 7.

2. Insights from previous work

This section provides a concise overview of the simulation framework and results from Paper I, revisiting its numerical formulation designed to test the GRG hypothesis. The section also summarizes key findings that underpin the present analysis.

2.1. Revisiting the simulation setup of Paper I

The five 3D simulations presented in Paper I feature an ambient medium characterized by a King’s β-profile, with the selected parameter values ensuring a stable galaxy group medium, thus maintaining a balance between gravity and pressure forces. In Cartesian coordinates, the density follows a profile

$$\begin{array}{cc}\hfill \rho (x^{\prime },y^{\prime },z)& ={\rho }_0{(1+\frac{{(x^{\prime }-x_0)}^2}{a^2}+\frac{{(y^{\prime }-y_0)}^2}{b^2}+\frac{{(z-z_0)}^2}{c^2})}^{-3\beta /2}\hfill \end{array}$$(1)

$$\begin{array}{cc}\hfill x^{\prime }& \equiv x\cos {10}^{°}-y\sin {10}^{°}\hfill \end{array}$$(2)

$$\begin{array}{cc}\hfill y^{\prime }& \equiv x\sin {10}^{°}+y\cos {10}^{°},\hfill \end{array}$$(3)

where, ρ₀ represents the core density (0.001 amu/cc, i.e., in atomic mass/cm³), while x′ and y′ denote axes rotated by 10° in order to break the symmetry of the jet flow within the ambient medium (insights on other parameters reveal as we discuss further). In all simulations, a cylindrical under-dense jet with a density of ρ_j = 10⁻⁵ × ρ₀, a radius of 1 kpc (corresponding to the resolution of the computational domain), and an initial length of 3 kpc is injected into the medium. The setup considers one-sided jet injection and propagation along the negative x-direction within a simulation domain extending ∼700 kpc in the jet-flow direction.

The choice of a 10 degree offset (in Eq. (1)) is introduced to mimic jet’s propagation in a more realistic environmental configuration rather than a perfect symmetry. Immediately after jet injection (after 3 kpc; length of the injection cylinder), the jet encounters a medium whose major/minor axis is misaligned with the jet-flow axis by approximately one cell at 5 kpc – consistent with the resolution of the computational grid. This asymmetry becomes increasingly significant with distance, reaching offsets of ∼9 cells at 50 kpc and ∼44 cells at 250 kpc.

The only varying conditions across the five runs are (a) the configuration of the ambient medium, which is tri-axial, with jets propagating either from its center or its edges, and (b) the classification of the jets into two types based on power, Fanaroff & Riley (1974) I and II types, injected into the computational domain. This formulation has been designed to test two hypothesized models for GRG formation (Dabhade et al. 2023): one suggesting that the dominant factors contributing to their growth is the more rarefied ambient medium, and the other proposing that the high jet power characteristics may have enabled the jet to overcome the frustration imposed by the surrounding environment.

A schematic diagram illustrating the five models is presented in Fig. 1. The simulations are labeled such that their nomenclature indicates the direction of the jet flow within the ambient environment and the associated jet power¹. We inject jets with low power, L_j = 2.3 × 10⁴⁴ erg/s, denoted by “lp” (Lorentz factor 3), and jets with high power, L_j = 7.2 × 10⁴⁴ erg/s, denoted by “hp” (Lorentz factor 5), corresponding to powerful FR-I jets and typical FR-II jets, respectively (Schoenmakers et al. 2000; Machalski et al. 2008). The acronyms “min”, “maj”, and “edge” refer to the jet’s injection position and flow direction within the tri-axial environment. Specifically, “min” indicates a jet originating at the center and propagating along the minor axis (encountering less obstacles in its path), “maj” represents a jet flowing towards the major axis (facing significant obstacles in its path), and “edge” denotes a jet injected from and propagating along the edges (x₀ = 0, z₀ = 0, y₀ = −6; referring jet’s injection at 600 kpc from the center) of the tri-axial medium (experiencing minimal obstacles in its path). The tri-axiality and axes configuration of the ambient medium are characterized by the parameters a, b, and c (the effective core radii), with values varying between 33 kpc and 66 kpc (“min”: a,  c ≡ 33 kpc and b = 66 kpc, “maj”: b,  c ≡ 33 kpc and a = 66 kpc, “edge”: a,  b,  c ≡ 66 kpc; refer to Paper I).

Fig. 1. Schematic representation of the jet-environment interaction framework developed to investigate the influence of factors on the development of GRGs. The setup simulates a one-sided relativistic jet propagating up to ∼700 kpc through a stratified medium representing a galaxy-group environment, capturing the varying influence of the environment on the jet’s motion. We also used two different jet powers to mimic Fanaroff-Riley type I and II jets, allowing us to examine the role of the jet’s primary thrust in the development of such radio galaxies. This parameter space yields a total of five simulation scenarios, discussed in detail in Paper I and also concisely in Section 2.1.

Figure 1 summarizes our simulation runs by a schematic representation, highlighting the environmental structure, the direction of the jet flow, the simulation domain, and the configuration of the coordinate axes. This figure includes all five simulation setups considered.

The simulations were performed in dimensionless units. A conversion to physical units can be conducted based on three fundamental scales: the length, L₀ = 100 kpc, the velocity, V₀ = c, i.e. the speed of light, and the density, ρ₀ = 0.001 amu/cc.

Other units can be derived using these fundamental scales, including the time, T₀ = L₀/V₀ = 0.326 Myr, the gas pressure, P₀ = ρ₀V₀² = 1.5 × 10⁻⁶ dyn/cm², and the magnetic field strength, $B_0=\sqrt{4\pi {\rho }_0{V}_0^2}=4.3\times {10}^{-3}$ Gauss.

On all sides of the simulation domain an outflow boundary condition is applied, except for the boundary layer hosting the jet injection nozzle (right boundary, along x-axis; Fig. 1), where we assign a reflective boundary condition. This setup is meant to mimic one half of a bipolar jet simulation, in particular for cases for which a substantial amount of backflow material reaches the jet injection region. The outflow boundaries apply a zero-gradient extrapolation which allows matter to freely leave the domain but does not artificially suppress inflow velocities. At earlier times, when the jet is still propagating well within the simulation volume, the surrounding environment near the boundaries of the simulation box remains effectively static, as governed by the combined balance of pressure and gravity set in the initial conditions for the ambient medium.

It is important to note that in all cases, the jet is injected (along negative x-axis) with a magnetic field in a purely toroidal configuration, which is defined as located in the plane perpendicular to the jet flow direction (i.e., in the y − z plane),

$$\begin{array}{c}\hfill B_y=B_{j}r\sin (\theta ),B_z=-B_{j}r\cos (\theta ),\end{array}$$(4)

where, (r,  θ) is the polar coordinate in the (y,  z) plane.

The strength of the injected magnetic field (B_j) is calibrated by assuming that the magnetic energy constitutes ∼1% of the enthalpy of the jet beam while injection (Mukherjee et al. 2020; Rossi et al. 2017). This ratio is referred to as σ, usually denoted as magnetization, and is defined as

$$\begin{array}{c}\hfill \sigma =\frac{B_j^2}{{\gamma }^2{\rho }_j{h}_j};{\rho }_j{h}_j=\frac{5}{2}{P}_j+\sqrt{\frac{9}{4}{P}_j^2+{({\rho }_j)}^2},\end{array}$$(5)

where ρ_jh_j is the matter-energy flux (enthalpy), P_j is the pressure, and γ the bulk Lorentz factor of the jet (subscript j for jet). Notably, the formulation of σ is based on the Taub-Matthews equation of state, commonly used to describe a relativistic gas (Taub 1948; Mignone et al. 2005). We choose, σ = 0.01, implying ⟨B_j⟩≃9 μG at the injection distance from the galactic center of 3 kpc.

We apply two tracer variables in order to be able to disentangle between material of the ambient medium (Q₁) and the jet (Q₂). Tracers are passive scalar variables that evolve according to the advection equation (Matthews et al. 2019)

$$\begin{array}{c}\hfill \frac{\partial (\rho Q_i)}{\partial t}+\mathrm{\nabla }·(\rho Q_{i}v)=0,\end{array}$$(6)

When tied to a fluid as an initial condition, their subsequent evolution is governed by the density and velocity, ρ,  v, of the associated fluid. This approach allows us not only to track the individual evolution of each fluid in a multi-fluid system, but also to determine the relative fraction of each fluid within a computational cell (value varying between 0 to 1). We emphasize, that this capability is essential for analyzing the matter entrainment processes that crucially influence the radiative appearance of GRGs.

2.2. Revisiting relevant conclusions of Paper I

For the reader’s convenience we repeat the following key conclusions from Paper I, as they are exceedingly relevant in the context of our current study.

1. We observed the emergence of distinct cocoon morphologies under different scenarios, evolving over varying dynamical timescales with jets of different power. This led to the conclusion that a combined information of morphology, jet power, and structural age of GRGs holds significant potential to unravel the underlying formation mechanisms.

2. The simulated lobes of giant radio galaxies are over-pressured across all models, indicating that lobe expansion remains active as long as the jet activity persists. The lateral (to jet-flow) variation in pressure showed a sharp jump in values at the cocoon-ambient medium interface irrespective of the models, following the overpressure nature of the lobes.

3. The axial ratio, defined as the ratio of cocoon length to width, reveals that powerful jets exhibit rapid, elongated expansion with arrowhead-like structures. In contrast, low-powered jets or jets encountering stronger environmental interaction produce wider lobes.

We keep these key insights in mind while further investigating the emission behavior of such radio galaxies and will return to these points when searching for correlations between the 3D dynamical evolution and the projected 2D emission maps.

3. Formulating radiative maps for the present investigation

Here, we outline the techniques adopted to convert the dynamical state of the simulated GRGs, as presented in Paper I, into emission and polarization maps. These maps are then used for further examination.

3.1. Post-processing I: Generating intensity maps

We converted the three-dimensional data of the simulated structures corresponding to the five simulation cases, as obtained at the evolved epochs detailed in Paper I, (see also Fig. 2, for a reference), into informative projected emission maps.

Fig. 2. Sliced velocity maps (x − y plane, z = 0) for the resultant GRG morphologies from five simulation cases (labels indicated in each plot) illustrating the underlying matter-transport mechanisms responsible for these structures. These maps are displayed at dynamical ages of ∼167 Myr (GRG_lp_min), 196 Myr (GRG_lp_maj), 49 Myr (GRG_hp_edge), 137 Myr (GRG_hp_maj), and 69 Myr (GRG_hp_min), as also highlighted in Paper I. The colormap represents the x-component of velocity, capturing the bulk-flow (in blue) and back-flow (in red) dynamics, while the contours depict the y-component of velocity, with solid lines indicating 0.005c and dashed lines indicating −0.005c, to highlight the behavior of secondary structures forming within the domain. To analyze the impact of jet power on morphology and the influence of the tri-axial environment on the resulting structure, comparisons are made column-wise and row-wise, respectively.

To evaluate the synchrotron emissivity of the simulated GRGs, we build upon the framework outlined in Vaidya et al. (2018), while adopting a more simplified dynamic-to-emission conversion methodology as described in Meenakshi et al. (2023). The foundational approach employed here excludes the effects of non-thermal particle re-acceleration and radiative cooling but accounts for naturally evolving dynamical processes such as matter entrainment and adiabatic cooling. The choice of this methodology, i.e. adopting a simplified form of Bessel function in emissivity evaluation, stems from the resolution requirements necessary to accurately incorporate all micro-scale processes (Dubey et al. 2023), which are nearly unfeasible with most contemporary computational resources, particularly for giant radio galaxy simulations (as also discussed in detail in Appendix A.1. of Paper I). This limitation has led to the recommendation of utilizing 2D models (Kundu et al. 2022). However, as noted earlier, reduced dimensionality and imposed symmetry often result in deviations from the outcomes of fully 3D models (e.g., Mignone et al. 2010), which justifies the adoption of this approach for the present study.

The emissivity equation for an ensemble of electrons (mass: m_e, charge: e) distributed in a power-law fashion over an energy spectrum (E_i to E_f), defined through Lorentz factors of γ_i (10²) to γ_f (10⁶) with a fixed power-law index p ( = 2.2; Machalski 2011), can be expressed in a co-moving frame as

$$\begin{array}{c}\hfill {\mathcal{J}}_{\mathrm{syn}}^{\prime }({\nu }^{\prime },\widehat{{\mathit{n}}^{\mathbf{\prime }}},{\mathit{B}}^{\mathbf{\prime }})=\frac{\sqrt{3}{e}^3}{4\pi m_e{c}^2}|{\mathit{B}}^{\mathbf{\prime }}\times \widehat{{\mathit{n}}^{\mathbf{\prime }}}|{\displaystyle {\int }_{E_i}^{E_f}}{\mathcal{N}}^{\prime }(E^{\prime })F(x)d{E}^{\prime },\end{array}$$(7)

where ν′ is the emitting frequency, and vectors B′ and $\widehat{\mathit{n}}\prime$ are the system’s magnetic field and the line of sight direction, respectively. Assuming isotropic conditions for the particle distribution 𝒩′(E′) and incorporating the Bessel integrals for the Bessel function, F(x) with x ≡ ν′/ν′_cr, the ratio of the emitted to the critical synchrotron frequencies, the expression for the emissivity simplifies to

$$\begin{array}{c}\hfill \begin{array}{c}\hfill {\mathcal{J}}_{\mathrm{syn}}^{\prime }=N_0\frac{3^{p/2}{e}^2{\nu }^{\prime -(p-1)/2}{|{\mathit{B}}^{\mathbf{\prime }}\times \widehat{{\mathit{n}}^{\mathbf{\prime }}}|}^{(p+1)/2}}{2c(p+1)}\\ \hfill \mathrm{\Gamma }(\frac{p}{4}+\frac{19}{12})\mathrm{\Gamma }(\frac{p}{4}-\frac{1}{12}){\left(\frac{e}{2\pi m_{e}c}\right)}^{(p+1)/2}.\end{array}\end{array}$$(8)

Here, N₀ (the constant in particle spectral distribution) can be derived by assuming that the energy stored in the non-thermal particles is equal to a fraction of energy carried by either magnetic field, E_B (Giri et al. 2022a,b; Upreti et al. 2024), or by the internal energy, E_th (Mukherjee et al. 2021; Meenakshi et al. 2023), of the simulated system. Adopting the latter, N₀ can be expressed as follows,

$$\begin{array}{c}\hfill N_0=\frac{ϵE_{\mathrm{th}}}{m_e{c}^2{\displaystyle {\int }_{{\gamma }_i}^{{\gamma }_f}}{\gamma }^{(-p+1)}d\gamma },\end{array}$$(9)

where ϵ indicates the fraction of internal energy (E_th = gas pressure/(Adiabatic gas constant − 1)) carried by the non-thermal particles, adopted to be 0.1, following Meenakshi et al. (2023).

We note that all primed variables are defined in the co-moving frame and require conversion to the observer’s frame as follows,

$$\begin{array}{cc}\hfill {\mathcal{J}}_{\mathrm{syn}}(\nu ,\widehat{\mathit{n}},\mathit{B})& ={\mathcal{D}}^2{\mathcal{J}}_{\mathrm{syn}}^{\prime }({\nu }^{\prime },\widehat{{\mathit{n}}^{\mathbf{\prime }}},{\mathit{B}}^{\mathbf{\prime }})\hfill \end{array}$$(10)

$$\begin{array}{cc}\hfill {\nu }^{\prime }& =\frac{1}{\mathcal{D}}\nu \hfill \end{array}$$(11)

$$\begin{array}{cc}\hfill \widehat{{\mathit{n}}^{\mathbf{\prime }}}& =\mathcal{D}[\widehat{\mathit{n}}+(\frac{{\gamma }^2}{\gamma +1}\mathit{\beta }·\widehat{\mathit{n}}-\gamma )\mathit{\beta }]\hfill \end{array}$$(12)

$$\begin{array}{cc}\hfill {\mathit{B}}^{\mathbf{\prime }}& =\frac{1}{\gamma }[\mathit{B}+\frac{{\gamma }^2}{\gamma +1}(\mathit{\beta }·\mathit{B})\mathit{\beta }],\hfill \end{array}$$(13)

where 𝒟 is the Doppler factor and β is the velocity vector in units of speed of light. We thereby generated intensity maps (I_ν) at a frequency of 1 GHz, which is relevant for several contemporary radio telescopes, by integrating the emissivity-data for the line of sight along the z-axis, in order to capture the giant radio galaxy’s behavior in the plane of the sky, as has been the case for a majority of GRGs (Oei et al. 2023).

3.2. Post-processing II: Generating polarization maps

Building on the definitions defined in Section 3.1, we utilized the expression for polarized emissivity as reported in Vaidya et al. (2018),

$$\begin{array}{c}\hfill {\mathcal{J}}_{\mathrm{pol}}^{\prime }({\nu }^{\prime },\widehat{{\mathit{n}}^{\mathbf{\prime }}},{\mathit{B}}^{\mathbf{\prime }})=\frac{\sqrt{3}{e}^3}{4\pi m_e{c}^2}|{\mathit{B}}^{\mathbf{\prime }}\times \widehat{{\mathit{n}}^{\mathbf{\prime }}}|{\displaystyle {\int }_{E_i}^{E_f}}{\mathcal{N}}^{\prime }(E^{\prime })G(x)d{E}^{\prime },\end{array}$$(14)

where the Bessel function is G(x) = xK_2/3. Utilizing the simplified form of the Bessel integral (Rybicki & Lightman 1979), the polarized emissivity in the observer’s frame simplifies to (Meenakshi et al. 2023),

$$\begin{array}{c}\hfill {\mathcal{J}}_{\mathrm{pol}}={\mathcal{D}}^2{N}_0\frac{3^{p/2}{e}^2{\nu }^{\prime -(p-1)/2}{|{\mathit{B}}^{\mathbf{\prime }}\times \widehat{{\mathit{n}}^{\mathbf{\prime }}}|}^{(p+1)/2}}{8c}\\ \hfill \mathrm{\Gamma }(\frac{p}{4}+\frac{7}{12})\mathrm{\Gamma }(\frac{p}{4}-\frac{1}{12}){\left(\frac{e}{2\pi m_{e}c}\right)}^{(p+1)/2},\end{array}$$(15)

where the conversion of ν′, B′, $\widehat{\mathit{n}}\prime$ follows the equations highlighted in Section 3.1. The Stokes parameters for linear polarization (U_ν and Q_ν) can subsequently be computed along the line of sight (z-axis) as

$$\begin{array}{cc}\hfill U_{\nu }(\nu ,X,Y)& ={\displaystyle {\int }_{-\infty }^{\infty }}{\mathcal{J}}_{\mathrm{pol}}\text{sin}2\mathcal{X}dZ\hfill \\ \hfill Q_{\nu }(\nu ,X,Y)& ={\displaystyle {\int }_{-\infty }^{\infty }}{\mathcal{J}}_{\mathrm{pol}}\text{cos}2\mathcal{X}dZ,\hfill \end{array}$$(16)

where the local polarization angle 𝒳 is defined as (following, Del Zanna et al. 2006),

$$\begin{array}{c}\hfill \begin{array}{c}\hfill \text{sin}2\mathcal{X}=-\frac{2q_X{q}_Y}{q_X^2+q_Y^2};\text{cos}2\mathcal{X}=\frac{q_X^2-q_Y^2}{q_X^2+q_Y^2}\end{array}\end{array}$$(17)

and

$$\begin{array}{c}\hfill \begin{array}{c}\hfill q_X=(1-{\beta }_Z)B_X-{\beta }_X{B}_Z;q_Y=(1-{\beta }_Z)B_Y-{\beta }_Y{B}_Z.\end{array}\end{array}$$(18)

The X and Y axes formulate the sky-plane in the observer’s frame, where the projected maps are generated. The angle 𝒳 is measured clockwise from the Y-axis, given the axis is aligned towards the north in the sky-plane. Magnetic field lines are then overlaid on the Stokes-I image (intensity), with the polarization angle 𝒳 being rotated by 90°. The length of the field lines is made proportional to the fractional polarization values, denoted as ∏, at each point. The fractional polarization, ∏, at the locations of synchrotron emissions is defined in terms of the Stokes parameters as,

$$\begin{array}{c}\hfill \prod =\frac{\sqrt{U_{\nu }^2+Q_{\nu }^2}}{I_{\nu }}.\end{array}$$(19)

4. Exploring velocity patterns in GRGs

In this section, we extend the analysis from Paper I by examining the velocity patterns of the generated giant radio galaxy structures to understand how different model assumptions result in different cocoon structures.

The velocity plots shown in Fig. 2 are comprehensive yet straightforward to interpret. The x-component of velocity is represented in a colormap, where the v_−x component, indicating the active jet’s bulk flow, is shown in blue shades, while the backflowing matter is displayed in red shades via v_+x, both on a logarithmic scale for enhanced visualization.

In addition, in order to capture the significance of matter flowing along the y-direction in certain cases – critical for revealing the full structural details – we have overlaid contours of the y-component of velocity at values of 0.005c (v_+y; solid lines) and −0.005c (v_−y; dashed lines), where c being the light speed.

We may follow two approaches in order to comparing the distinct evolutionary features of these runs, as we discuss below.

4.1. Impact of jet power

A comparison between the low jet power models in the left column and the high jet power cases in the right column highlights how the differences in jet power may influence the general dynamical behavior. For example, we observe that lower-powered jets (FR-I-like) tend to form wider lobe structures. This can be attributed to the jets’ inability to maintain collimation over longer distances. This effect is further aggravated by slight bending (along y-axis) caused by the external ambient medium, as the medium is rotated by 10° to break symmetry. An illustration of the initial configuration of the environment (for the ‘maj’ cases) is presented in Fig. 3, highlighting the asymmetry introduced between the major axis and the jet propagation axis.

Fig. 3. Initial density configuration of the ambient environment, illustrating a 10° rotation of the ambient medium’s principal axis, in order to introduce an asymmetry to the jet-flow, while the jet injection axis remains fixed along the negative x-direction. The density distribution is shown using colormap, with the corresponding log ρ[amu/cc] values annotated along the contour lines. This represents one of the simulation scenarios employed, as shown in Fig. 1.

Together, these effects lead to a de-collimated, broad flow that distributes material over a wider opening angle (similar to e.g., Monceau-Baroux et al. 2015; Nawaz et al. 2016). As a result of this divergence, the bulk flow speed of the fluid significantly decreases along the −x axis. This reduction in speed also limits the formation of backflow material, causing the material to propagate more slowly, which thereby struggles to reach back the location where the jet originated. The combination of a wider bulk flow and a slowed backflow ultimately contributes to the formation of significantly wider lobes (refer to the v_y contours). In this context, we note that the presence of mild back-flowing plasma in FR-I jets has been documented by Laing & Bridle (2012), exhibiting characteristics consistent with our findings.

In contrast, the cases in the right column, characterized by higher flow speeds (FR-II-like), consistently exhibit well-collimated jet beams up to the jet head. This indicates rapid propagation, as the jet thrust is concentrated over a smaller area (Bromberg et al. 2011). Consequently, this results in a significant pressure imbalance at the jet head–ambient medium interface, generating a higher volume of backflowing material, even at greater speeds (Cielo et al. 2017).

As a result, the backflowing material stretches out predominantly along the jet flow direction. To note, the rapid expansion we observe for the high-powered jets is simply evident from the fact that they are able to cover a linear extent of nearly 700 kpc within 50–140 Myr, whereas the low-powered jets reach a smaller linear extent of only up to 600 kpc, despite evolving over a longer duration of 170–200 Myr (Paper I).

4.2. Impact of the environment

We now investigate the influence of the environment. In Fig. 2, this means comparing row-wise. We observe distinctions arising from the influence of the ambient medium’s configuration. In the minor-axis cases shown in the top row, we observe that for a low-powered jet, decollimation leads to the lobe expanding laterally, aided by the pressure gradient force that is strongest along the minor axis (similar to, e.g., Giri et al. 2023; Bruno et al. 2024). As a result, the lobe’s lateral size is largest in this case (GRG_lp_min) among all the models.

Conversely, for the high-powered case (GRG_hp_min), the jet remains collimated and maintains high thrust, enabling rapid propagation (again, aided by the pressure gradient). Consequently, the backflow material struggles to reach all the way back to the launching region, ultimately producing an arrowhead-shaped structure.

In the middle row (Fig. 2), the jet propagates along the major axis of the environment. With greater obstruction to the jet flow due to the lower rate of change in pressure values (low pressure gradient) along this direction, the low jet power run exhibits lobes that are not as wide as those in the ‘GRG_lp_min’ case. However, as the jet enters into the computing domain with the same power as ‘GRG_lp_min’ and evolves for a longer period (196 Myr versus 167 Myr), subsequently failing to spread laterally, generates substantial backflow.

This backflow returns to the zone of jet-injection and then follows the path of maximum pressure gradient, i.e., along the minor axis of the environment. This expansion of backflowing material in the shape of a wing is further captured in the y-component of velocity (Fig. 2) and results in a structure known as X-shaped radio galaxies (XRGs; Giri et al. 2024). This mechanism of XRG formation is referred to as the classic backflow scenario (Capetti et al. 2002; Rossi et al. 2017) and is believed to occur in many (giant) radio sources with X-shape (Saripalli & Subrahmanyan 2009; Cotton et al. 2020; Bruni et al. 2021).

For high-powered jets flowing along the major axis, a similar X-shaped structure is formed. However, the prominence of the secondary lobes, or wings, is reduced since the jet propagates rapidly (less decollimation), leading to a decrease in the wing-to-lobe length ratio (ratio < 0.8; Cheung 2007). Notably, due to the jet’s fast propagation, even though a significant amount of backflow is generated – as expected for higher-powered jets – the bulk fraction of matter struggles to return quickly to the jet injection zone. Consequently, a substantial amount of material accumulates within the active lobe cocoon (along −x axis), creating a morphology known as a nose-cone shape (Stone & Hardee 2000). The overlaid contours illustrate the y-velocity of gas-flow, showing explicitly the channeling of matter along the minor axis of the environment to form the wings.

The asymmetry in the v_y contours with respect to the jet-flow direction for ‘maj’ cases can be attributed to the effect of tilted ambient medium, which is inclined by 10° to break symmetry to the jet flow (Fig. 3). Furthermore, for the high-power cases, in general, the contour distribution reveals the formation of strong bow shocks in the environment due to the rapid expansion of the jet. In contrast, for the low-power cases, the bow shock is less prominent and remains almost attached to the entire frontal part of the lobe, reflecting slower progression of lobe expansion along the x-axis.

The bottom row highlights an intriguing case of GRG formation, where the jet flows along the edges of the environment with minimal obstacles, producing an arrow-head structure. The well-collimated jet propagates rapidly (covering ∼650 kpc in 49 Myr), generating substantial and fast backflows that return toward the jet origin. These backflows, influenced by buoyancy forces, create asymmetry in the cocoon’s lateral extent as plasma follows the pressure gradient outward (from the center of the ambient medium; Fig. 1). This supplies a significant clue for guessing underlying GRG formation, as many are suspected to evolve outside galaxy groups or cluster media. Observationally, GRGs mapped near voids or along large-scale density gradients often show bending in their lobes (e.g., Safouris et al. 2009; Malarecki et al. 2015), consistent with this result.

From these discussions, it is evident that the combined effects of the ambient medium and jet power drive the appearance of GRGs. Also noteworthy is the observation that in cases of high-powered jets, the colormap shows the presence of v_+x components in the immediate vicinity of the jet beam, suggesting that material may have leaked from the jet-column due to instabilities, and subsequently flow alongside the jet. This is likely caused by the formation of a shearing layer between the sheath material and the jet-spine, which facilitates the leakage through Kelvin-Helmholtz (KH) instabilities (Borse et al. 2021; Wang et al. 2023). This phenomenon could be an intriguing factor for further investigation, as it may provide additional insights into the survival and stability of jet flows over even greater distances.

5. Radio continuum maps at 1 GHz

This section provides a detailed analysis of the emission maps derived from the techniques outlined in Section 3. It is divided into three subsections: the first (Sect. 5.1) discusses the morphological features observed in the maps for each GRG case; the second (Sect. 5.2) presents a quantitative analysis of these maps, leading to inferences relevant to observational studies; and the final subsection (Sect. 5.3) explores the correlation between the obtained radiative results and the underlying three-dimensional dynamical state.

5.1. Understanding GRG dynamics through emission maps

The generated emission maps effectively encapsulate the internal dynamical behavior of the GRG systems, presenting a more concise representation by reducing dimensionality. These maps integrate the combined influence of multiple three-dimensional properties into a single two-dimensional framework. Since the maps are integrated along the line-of-sight vector, any features distinctly observable in the emission maps indicate the prominence of those layers within the underlying 3D structure. Below, we discuss the projected sky maps classified according to their formation processes.

5.1.1. Jet propagation through regions of moderate obstructions

The 1 GHz emission maps shown in Fig. 4(A) reveal that extrinsic properties – such as the cocoon length, its geometric shape, and the overall projected axial ratio (length-to-width ratio) – exhibit consistent similarities with the dynamical maps (Paper I) for the chosen line-of-sight. Fig. 4(A) contains the evolved structures at ages of nearly 167 Myr and 69 Myr for a low-powered jet (“GRG_lp_min”, top panel) and a high-powered jet (“GRG_hp_min”, bottom panel), respectively, propagating along the environment’s minor axis direction¹.

Fig. 4. A. Intensity maps (log Iν) at 1 GHz radio frequency for two simulation cases: a low-power jet (“GRG_lp_min”; top) and a high-power jet (“GRG_hp_min”; bottom). Both jets are propagating along the minor axis of the ambient medium. While the overall structures resemble the dynamical maps from Paper I, the emission maps offer deeper insights, highlighting crucial aspects of the observability of these simulated giant radio structures at their evolved stages (∼167 Myr: top; ∼69 Myr: bottom).

The decollimation of the jet flow, combined with the counterclockwise bending of the bulk flow, induces lobe formation in the “GRG_lp_min” case. This suggests that GRG lobes may not be passive regions of the jet-cocoon system but rather zones of intensified turbulence (e.g., Kundu et al. 2022), fostering the development of decollimated filamentary structures and cascading multi-scale vortices. The decollimated flow extends to the cocoon head along the linear growth direction, indicating ongoing lobe expansion. This is corroborated by the presence of a mild bow shock adjoining the frontal lobe, as evidenced in the dynamical evolution outlined in Paper I. The well-defined boundary of the lobe in this case is particularly notable, suggesting that, in addition to the structure being 1.9 times over-pressured relative to the surrounding medium, Kelvin-Helmholtz (KH) instabilities have a lesser impact on disrupting the cocoon in such low jet power scenarios (low-magnetization cases; Mukherjee et al. 2020).

The differences in cocoon disruption become particularly evident when comparing with the high-power jet case (“GRG_hp_min”). Despite having an over-pressure of four times that of the environment, the cocoon edges in this case are more susceptible to matter entrainment by KH instabilities. This phenomenon arises from the combination of high jet power and the substantial density contrast between the jet and the surrounding medium, which together ensure the survival of the jet spine but compromise the stability of the cocoon (Rossi et al. 2008). These findings underscore the critical role of instabilities and entrainment processes in shaping the jet-cocoon dynamics, particularly in the long-term evolution and survival of giant radio galaxies.

In the “GRG_hp_min” case, the jet spine exhibits a bifurcated structure, attributed to the concentration of magnetic field lines along the jet beam’s edges, followed by the formation of recollimation zones. This behavior aligns with high-resolution numerical studies (Dubey et al. 2023) and is consistent with recent high-resolution, sensitive observations of extragalactic jetted radio sources (Condon et al. 2021; Velović et al. 2023; Weżgowiec et al. 2024).

The “GRG_hp_min” case is particularly notable for two key reasons. (a) It exhibits the transition of the jet spine into a decollimated, turbulent state with tendril-like extensions, providing insight into the jet-lobe transition process in giant radio galaxies (with relevance to smaller RGs, e.g., Condon et al. 2021; Velović et al. 2023; Weżgowiec et al. 2024). (b) The decollimation of the jet spine at ∼450 kpc, despite the cocoon extending to 700 kpc, suggests that hotspots (or enhanced emission zones) located well within the cocoon of observed GRGs may not always indicate a restarting jet phase (e.g., Sebastian et al. 2018; Andernach et al. 2021; Oei et al. 2022). This observation may challenge the conventional interpretation of hotspot positions as definitive indicators of restarting activity in GRGs (see the numerical models, e.g.: Horton et al. 2023; Young et al. 2024). Decollimation in this high-power jet case yields warmspot (wide emission zone) instead of a distinct hotspot, along with tendril-like filamentary structures, further analyzed in Sections 5.1.4 and 5.3.

5.1.2. Jet propagation through regions of maximal obstructions

Due to the jets propagating along a path with a lower pressure gradient (i.e., the major-axis direction) and encountering greater resistance to free flow in both the low jet power (top panel, Fig. 4(B)) and high jet power cases (bottom panel, Fig. 4(B)), the formation of lobe structures follows. This obstruction causes the jet to decelerate faster, resulting in earlier decollimation into tendril-like extensions compared to the jet flow along the minor axis situations, where collimation is preserved over a greater distance (∼300 kpc vs. 400 kpc). Both scenarios exhibit noticeable wing emission, highlighting the influence of lateral deflection of backflow material from the active lobe. Additionally, the interaction of matter, likely originating from the opposite jet arm, contributes to this effect (modeled using a “reflective” boundary condition near the jet injection boundary to replicate this natural phenomenon; Capetti et al. 2002). These effects cause an increase in internal energy within the wings highlighting the fact that wings are not passively evolving structures in XRGs. Instead, they function as regions of enhanced turbulence and random shock sites. This finding has important implications for understanding the persistent emission and anomalous spectral behavior observed in the wings of many XRGs, which is believed to be caused by particle re-energization, as indicated by both numerical studies and observations (Giri et al. 2022b; Gopal Krishna & Dabhade 2022; Patra et al. 2023).

Focusing on the evolved structure (at ∼196 Myr) of the low jet power case (“GRG_lp_maj”), we observe notable behavioral features, such as the jet’s bending in response to the ambient medium’s 10° rotational setup to break symmetry (Fig. 3). This interplay further induces the formation of a large vortex-like structure in the lobe, with jet material rotating in a clockwise direction (as also inferred from velocity streamlines; not shown here). Additionally, this motion leads to the creation of smaller-scale caldera-like structures, characterized by emission valleys surrounded by enhanced emission ridges. The compression of jet-matter (consequently the magnetic field; Section 6) in the vicinity of lobe edges is particularly prominent (e.g., Condon et al. 2021). These dynamic processes ultimately contribute to the development of the lobe structure observed along the active jet axis.

In the high-powered jet case (“GRG_hp_maj”), the jet resists bending due to its higher thrust, and upon decollimation, the jet beam disperses the jet material over a wider angle forming spread-out emission ridges (e.g., Cotton et al. 2020). A slight asymmetry in the channeling of back-flow material has been observed in the lobe, with stronger channeling occurring towards the southern side with respect to the jet axis. Despite the jet’s decollimation, the thrust appears to persist up to the frontal region, forming a top-hat-like structure (see Section 4). The sustained thrust of the jet results in an intensified interaction with the ambient medium in the top-hat region of the jet, which is further discussed in Section 5.1.4.

As the jet travels larger distance compared to “GRG_lp_maj”, it loses internal energy and magnetic field strength more due to expansion, leading to a reduced contribution to the emission map from the frontal lobe part (compared to “GRG_lp_maj”; Fig. 4(B)). This emission map also reiterates that enhanced emission regions, deeply embedded within the structure, do not necessarily imply a jet restart (Young et al. 2024; Fanaroff et al. 2021). This underscores the need for additional diagnostics, such as pressure jumps at lobe boundaries (see Section 5.2.1), to reliably differentiate between jet phases. The cocoon is susceptible to greater disruption in this case due to increased shear with the surrounding environment.

5.1.3. Jet propagation through regions of minimal obstructions

This scenario involves a high-powered jet propagating through the edges of the ambient medium (“GRG_hp_edge”, Fig. 4(C)), where minimal flow resistance allows for rapid propagation, as evidenced by the jet covering ∼670 kpc (one-sided) in ∼49 Myr. At the same time, this also results in a higher shear between the cocoon and the surrounding environment, significantly affecting the visibility of the cocoon’s extent. Due to the high thrust, the jet decollimates gradually (compared to “GRG_hp_maj”), generating finger-like filamentary extensions. Additionally, the increased jet speed leads to the leakage of matter from the jet spine, a feature observed as a result of enhanced shear between the jet sheath and spine, which fosters instability growth (see, Section 4). The combination of these effects results in the mixing of matter from different regions, such as thermal and non-thermal particles, as well as freshly injected and aged particles (Turner et al. 2018; Rossi et al. 2024). Such mixed flow, along with backflow and buoyancy effects, creates greater turbulence within the cocoon, which explains the enhancement of emission closer to the jet-launching zone rather than near the jet head (e.g., PKS 2356-61; Ursini et al. 2018). The prominent emission zone originating from the jet also diminishes midway along the total linear extent of the cocoon in this case. Despite the jet’s bulk flow experiencing minimal decollimation among the high jet power runs, its propagation through a rarefied environment results in a mini warm spot formation at the jet head position (see Section 5.1.4).

5.1.4. Hotspots and jet dynamics

It is worth noting that our emission maps lack prominent hotspots at the jet head in high-power cases. Instead, peak emission appears midway along the jet path, despite sustained thrust to the edge. Hotspots, typically marking the termination shock where powerful jets interact strongly with the ambient medium (Bromberg et al. 2011), form as jets decelerate from near-light-speed at launching sites to ≲0.1c at lobe fronts (for GRGs; Jamrozy et al. 2008). However, their radiative prominence also critically depends on the distribution and strength of the magnetic field (Black et al. 1992; Pyrzas et al. 2015; Baghel et al. 2023).

In Fig. 5, we present dynamical maps from our high-powered GRG jet simulations. The top panels show the 3D tracer distribution, revealing the evolution of the jet spine and the shape of the head region that primarily drives the backflowing plasma. The bottom panels show 2D slices focusing on the jet head, displaying the B-field magnitude and pressure distribution. Notably, around the midpoint of the jet’s extent, the jet spine begins to lose its collimation and develops mild undulations or wiggling. This induces matter leakage into the surrounding cocoon, generating a spine-sheath morphology beyond mid-distance. Concurrently, the jet head in these high-power cases produces substantial backflowing material.

Fig. 5. B. Continued. Intensity maps (log Iν) at 1 GHz radio frequency depicting the X-shaped radio galaxy in their giant stages, resulting from jet propagation along the direction of highest jet frustration (i.e., along the major axis direction). The low jet power case (“GRG_lp_maj”; top) and the high jet power case (“GRG_hp_maj”; bottom) exhibit distinct decollimation patterns, with both cases generating internal regions of emission valleys and enhanced emission ridges. The wings are luminous, indicating their active evolutionary stages (∼196 Myr: top; ∼137 Myr: bottom).

This continuous leakage and backflow transport B-field lines into the cocoon, which leads to a dilution of the magnetic field strength at the jet head by an order of magnitude or more in comparison with the jet injection zone (weakening the emission) – even while the thrust remains largely intact, as seen in the pressure maps (Fig. 5). For instance, in the GRG_hp_maj case, where the jet experiences greater obstruction by the environment and forms a broader cocoon, the magnetic field becomes more evenly distributed across the cocoon. This explains why its emission map shows lobes illuminated over a wider region near the jet head, accompanied by the formation of a warm spot (Leahy et al. 1997) in the top-hat region. By contrast, in the other two high-power runs, the B-field remains concentrated around the jet sheath, leading to emission that is largely confined along the jet beam and gradually fades as it approaches the jet head. These results point to the possibility that some fraction of GRGs are likely larger than current observations indicate. As radio surveys advance in sensitivity and resolution, they could expose the full extents of these sources and identify additional GRGs now misclassified as smaller radio galaxies (Hardcastle et al. 2019; Koribalski 2025).

To further assess the weakening of magnetic field strength along the jet, particularly at the terminal regions, we computed the plasma beta (β = 8πP_jet/B_jet²) in the vicinity of the jet head. Our analysis reveals β ≳ 20 in this region, indicating that the magnetic pressure plays a subdominant role compared to the thermal pressure at the jet terminus. Notably, in the cocoon enveloping the head, β rises even further–reaching values ∼500 or higher–suggesting significant magnetic field dilution. This trend, along with the increasingly tangled nature of the magnetic field topology, is further illustrated in Fig. 13 in Section 6.

Such features recovered in our GRGs – including the emergence of top-hat–shaped warm spots, and peak emission zones forming well within the lobe boundaries – have also been documented in several other observed GRG cases, such as: Clarke et al. (2017), Chen et al. (2018), Sebastian et al. (2018), Oei et al. (2024a). However, to better understand the typical occurrence of hotspots at the jet termini in powerful GRGs, further numerical investigations are needed (e.g., through high-resolution simulations of jet flows capable of modeling flow bifurcation, branching, and complex spine–sheath dynamics; Laing et al. 2015; Horton et al. 2023). As suggested in Paper I, a comprehensive understanding of these complex emission features requires systematic exploration across a broader parameter space – such as higher jet-to-ambient density contrasts (ρ_j/ρ₀) and stronger injected magnetic fields – awaiting future study. In this context, we note that we have initiated a related investigation focused on the early evolutionary stages of GRGs (evolving on the galactic scales), examining the role of magnetic field strength in jet collimation and consequent hotspot formation (Giri et al., in prep.). Additionally, the adoption of advanced emission mapping techniques, such as a Lagrangian-Eulerian particle-fluid framework on high-resolution grids, is planned in future work, facilitated by recent GPU-accelerated code developments.

5.2. Quantitative exploration of the emission maps

5.2.1. Insights from the lateral intensity variation

We begin with the intensity profiles, as shown in Fig. 6, illustrating the lateral emission variation of the simulated structures at their evolved stages. The variation is assessed at the maximum cocoon width, representing a general demonstration of the underlying radio galaxy (Subrahmanyan et al. 1996; Malarecki et al. 2013; Hardcastle 2018). To ensure a fair comparison within the cocoon and exclude the intensity variations around its immediate edges, we apply a threshold of 10²³ erg s⁻¹ Hz⁻¹ kpc⁻² sr⁻¹ to the intensity values. In all cases, the intensity values are notably higher around the jet-spine structure, reflecting its active flow and thrust. A shift in the peak emission position from the central axis or mid-point is observed in low jet power scenarios, attributed to the bending of the jets, as can be seen in the dynamical maps shown in Fig. 2.

Fig. 6. C. Continued. Intensity map (log Iν) at 1 GHz radio frequency revealing an intriguing case where the jet flows with minimal hindrance due to its propagation along the edges of the surrounding environment (“GRG_hp_edge”). The faster propagation causes increased shear between the cocoon and the environment as well as between the jet spine and sheath layers and results in the generation of instabilities, leading to the formation of pronounced turbulent structures at 49 Myr.

The most notable feature in this lateral intensity analysis (Fig. 6) is the intensity variation observed in regions extending outward from the jet spine. For both high and low jet power cases, the emission decreases sharply from a relatively high value to the threshold intensity as one moves laterally from the jet spine (Condon et al. 2021). This suggests that as long as the jet remains active, the cocoon structure retains its dominance across a substantial lateral extent of the radio galaxy. Conversely, as the radio galaxy slows in its evolution (e.g., through adiabatic expansion or cessation of jet activity) or undergoes significant entrainment, the emission gradually diminishes outward from the primary flow, as suggested by the observational findings of Subrahmanyan et al. (2008). This gradual reduction in lateral emission may therefore serve as an indicator of intermittent activity or lower active stages of radio galaxies, complementing other diagnostic parameters such as core dominance or hotspot position estimation (e.g., Subrahmanyan et al. 2008; Cotton et al. 2020).

5.2.2. Insights into the radio power–jet power correlation

Motivated by the growing interest in the radio power versus linear size (P-D) distribution for GRGs (Andernach et al. 2021; Simonte et al. 2022), we analyzed this trend (see Fig. 7) using our five simulated outputs at their evolved stages (one-sided length exceeding 350 kpc). Considering the evolution track of each simulations, for the lower jet power scenarios (depicted in blue), we observe an increase in radio power with increasing source age. The high jet power situations (depicted in red), however, exhibit a variable trend, with total radio power sometimes increasing or decreasing with age. A closer examination suggests that this behavior correlates with the extent and shape of the radio galaxy, i.e., sources with wider and more extended structures, also indicative of prominent lobe formation, tend to exhibit higher radio power for a given jet power. This trend may not be universally applicable to all GRGs (Ishwara-Chandra & Saikia 1999), as the influence of radiative cooling becomes increasingly significant with the source’s growing size and redshift evolution – an effect not explicitly accounted for in our calculations. In the following sections, we briefly examine the relative dominance of various cooling mechanisms in shaping the radiative properties of GRGs.

Fig. 7. Dynamic maps illustrating the jet spine and head structures, elucidating the absence of prominent hotspots in our synthetic emission maps for high-power jet simulations. The top panel presents a 3D volume-rendered visualization of the jet tracer, highlighting the jet beam – jet head morphology at their corresponding evolutionary ages. The bottom panel shows 2D slices focused on the jet head: the top row maps the magnetic field magnitude, while the bottom row maps the pressure distribution. These maps reveal that although the jet sustains its thrust up to the terminal shock, the magnetic field strength drops sharply at this region – largely advected away by the backflowing plasma – thereby suppressing the formation of bright hotspots.

Our observation of Fig. 7 also suggests that sources with greater linear extents do not exhibit any trend in radio power values, with some showing comparable or even lower radio power estimates (e.g., Dabhade et al. 2020; Simonte et al. 2024). This observation hints that the combination of radio power and linear extent alone is insufficient to fully capture the complexities of processes within a GRG system. It underscores the importance of considering additional metrics, such as the relationship between radio power and the projected emission area of the cocoon (P-A), to gain deeper insights (also reiterating from the above paragraph).

The P-D analysis (Fig. 7) also emphasizes a key conclusion: using a simplistic analytical formulation to relate radio power and jet power (e.g., Hardcastle 2018) for statistical GRG samples (Dabhade et al. 2020) may lead to inaccurate estimates. This is due to the formulation’s dependence on various factors, including the complex geometry and associated parameters of the ambient environment, which significantly influence the jet-induced morphology and radio power (Mingo et al. 2019). These factors, which vary depending on the parametric space where the GRG formed, make it unsuitable for straightforward application to individual GRGs.

The radio power values obtained at 1 GHz for the simulated GRGs are slightly on the lower side. Adopting a redshift of 0.05, the estimated flux values for these simulations are ∼135 mJy (“lp_min”), 67 mJy (“hp_min”), 73 mJy (“lp_maj”), 86 mJy (“hp_maj”) and 21 mJy (“hp_edge”), respectively. It is important to note that these values correspond to the one-sided cocoon only, and thus, the total values are likely to be higher. Additionally, we note that the emission maps were generated under the assumption that 10% of the internal energy is converted into emission (Section 3.1). The mentioned ranges of the estimated radio power and flux values, however, are relevant with those reported for several GRGs in similar frequency bands, as referenced in, Machalski et al. (2009), Machalski (2011), Andernach et al. (2021), Delhaize et al. (2021), Charlton et al. (2024). For reference, in Fig. 7, a representative one-sided radio power and length estimation of three GRGs from Charlton et al. (2024) is shown.

5.2.3. Magnetic field estimates and their implications

To infer magnetic field strengths critical for estimating GRG parameters such as age, we calculated and compared different B-field estimates, as shown in Fig. 8. The first estimate, shown as purple points, represents the median values of dynamically simulated magnetic fields (B_dyn) at the evolved times of the five simulations, observed to vary around 0.15 μG (see also, Paper I).

The second estimate represents the median magnetic field values calculated using the equipartition assumption (B_eq), which assumes equal energy distribution between radiating non-thermal particles and the magnetic field. Following Hardcastle et al. (2002), we deduce the B_eq values as follows,

$$\begin{array}{c}\hfill m_e{c}^2{\displaystyle {\int }_{{\gamma }_{\min }}^{{\gamma }_{\max }}}\gamma N(\gamma )d\gamma =\frac{B_{\mathrm{eq}}^2}{8\pi },\end{array}$$(20)

where N(γ) is the non-thermal electron distribution, with energy represented by the Lorentz factor γ, ranging from 10² to 10⁶, as discussed in Section 3.1. The equipartition field strengths range from 0.65 μG to 1.94 μG (Schoenmakers et al. 2000; Andernach et al. 2021; Machalski et al. 2008), with the ratio of dynamical to equipartition values as 0.15 (“lp_min”), 0.17 (“hp_min”), 0.18 (“lp_maj”), 0.08 (“hp_maj”), and 0.2 (“hp_edge”) (e.g., see the conclusion on powerful radio galaxies; Mahatma et al. 2020). For reference, we highlighted the B_eq values for a sample of GRGs from Schoenmakers et al. (2000) in a gray zone in Fig. 8 (see also, Tamhane et al. 2015; Charlton et al. 2024), illustrating their similarity to our estimated values, while the true intrinsic values are nearly an order of magnitude lower. The review by Dabhade et al. (2023) summarizes a range of B_eq values from 1 to 16 μG, with a median of 5 μG. However, when compared to X-ray measurements (independent of the equipartition assumption), this appears to be nearly an order of magnitude lower, further emphasizing our reported situation.

Fig. 8. Lateral intensity variation measured at the maximum cocoon width for five simulation scenarios (Fig. 4(A−C)), providing a quantitative assessment of intensity distribution and identifying active zones within the cocoon. Subplots correspond to different simulation runs categorized by ‘min’, ‘maj’, and ‘edge’ cases, enabling a comparison of jet power on cocoon activity and evolution.

With the assumption that the spectral age follows the relation, t_spec ∝ B^−3/2 (van der Laan & Perola 1969), the differences in equipartition values, which vary by nearly an order of magnitude from the dynamical values, can lead to significant discrepancies in the spectral age compared to the true dynamical age (with the spectral age being smaller). Thus, the inherent (dynamical) age of many GRGs may need to be re-evaluated, which might highlight discrepancies among various age estimation methods (see, e.g., Machalski et al. 2009; Machalski 2011). This re-assessment is crucial for testing the proposed formation mechanisms of GRGs under the standard jet evolution hypothesis and clarifying the inconsistencies observed across studies concerning the GRG ages, such as Jamrozy et al. (2008), Komberg & Pashchenko (2009), Sebastian et al. (2018), Marecki et al. (2021).

The third estimate is the CMB-equivalent magnetic field strength (B_CMB) at a redshift (z) of 0.05, calculated using the formula,

$$\begin{array}{c}\hfill B_{\mathrm{CMB}}=3.18\times {(1+z)}^2\mathrm{\mu }\mathrm{G}.\end{array}$$(21)

The equivalent level is shown with a dashed line in Fig. 8, indicating higher values than both B_dyn and B_eq, suggesting the dominant role of the inverse-Compton cooling process in GRGs, which remains significant even at this low-redshift regime (aligning with, for example, Ishwara-Chandra & Saikia 1999; Schoenmakers et al. 2000; Subrahmanyan et al. 2008).

5.2.4. Influence of matter entrainment on morphological appearance

The entrainment of ambient matter into the cocoon primarily arises from the significant velocity shear and density contrast between the relatively static ambient medium and the dynamic jet-cocoon system. This process is facilitated by the growth of instabilities, including Kelvin-Helmholtz (Borse et al. 2021) and Rayleigh-Taylor instabilities (Abolmasov & Bromberg 2023), combined with the effects of large-scale cocoon rotation (Gourgouliatos & Komissarov 2018). These mechanisms drive turbulent mixing, enabling the transport (beside momentum exchange) of ambient material deep into the cocoon and into the vicinity of the jet-beam structure, significantly impacting the radiative properties of the system.

In our simulated systems, we estimated the amount of total entrained ambient matter using the tracer Q₁ (utilizing the location and corresponding Q₁ fraction) within the boundary of the jet-induced cocoon, defined by Q₂ ≥ 10⁻⁷ (Mukherjee et al. 2020), as shown in Fig. 9. We reiterate that Q₁ and Q₂ are passive scalars used to trace different media: Q₁ for ambient gas and Q₂ for jet gas (see Section 2.1). The results reveal an initial rapid phase of entrainment, attributed to the jet’s propagation through the denser core region, followed by a steady increase in the entrained mass as the cocoon length increases. The encapsulated mass represents a considerable fraction of the underlying ambient medium (baryonic mass of the ambient medium ∼3.3 × 10¹² M_⊙) at the evolved times.

Fig. 9. Radio power at 1 GHz of the (one-sided) cocoon versus (one-sided) cocoon length, presented at the evolved stages (dynamical ages attached) of the simulated GRGs. Two jet power values used in the simulations are also highlighted, in blue and red. The evolutionary track of the simulated GRGs in the Power-Distance (P-D) diagram is also shown from the point they exceed 350 kpc, qualifying as GRGs in their total extent. For reference, we include measurements of three observed GRGs from Charlton et al. (2024) in similar observing band, representing a typical one-half estimation. We note that the above calculation does not explicitly account for the cooling effects.

In the low-powered scenarios (Fig. 9), the jet decollimation and wider lobe formation results in higher entrainment. This behavior is expected and is suspected to be a key factor in the jet’s rapid deceleration, leading to its transition from a collimated jet flow to an FR-I-like lobe morphology (Massaglia et al. 2016; Abolmasov & Bromberg 2023; Rossi et al. 2024). Similarly, in the high jet power case with increased jet opposition (“GRG_hp_maj”), where the formation of a mini-lobe structure is observed in the frontal part of the flow, the accumulation of ambient matter remains comparable. The other high-powered situations show lower levels of matter accumulation, as the high-speed flow helps maintain the collimation of the jet spine for a longer duration. However, due to the increased velocity shear between the ambient medium and the cocoon, and also being significantly less dense (with the injected jet density being 10⁻⁵ times that of the ambient density), the cocoon begins to experience disruption, although the central spine remains relatively protected (Rossi et al. 2008). A direct consequence of this phenomenon is illustrated in Fig. 4 (A: lower-panel), which shows the described jet-spine and cocoon behavior. The lowest mass entrainment occurs when the jet flows along the edges of the ambient medium (Fig. 9), resulting from its propagation through a lower-density region and the collimation of the flow due to rapid propagation. However, following the processes noted earlier, cocoon disruption is evident in this case also.

Interestingly, for the “GRG_hp_edge” case (Fig. 4(C)), the jet-spine is observed to leak matter into the cocoon-sheath structure (around ∼200 kpc), offering insight into the turbulent mixing of non-thermal, older back-flowing particles with freshly injected ones. The enhanced shear between the jet’s spine-sheath structure, along with instabilities from recollimation shocks, likely contributes to these features. While understanding these phenomena is crucial for spectral studies of radio galaxies (Turner et al. 2018), we have not focused on them in detail due to the need for significantly higher resolution for convergence, warranting a separate study. Previous works have highlighted the importance of such investigations in greater detail (Borse et al. 2021; Abolmasov & Bromberg 2023; Wang et al. 2023), which should be extended to GRGs, as jet stability is vital for their survival over large scales (since greater matter entrainment leads to enhanced jet deceleration; Rossi et al. 2024).

5.3. Correlating emission features with dynamical properties: Filamentary extensions

Observations of extragalactic radio jets with contemporary telescopes have unveiled novel features in a growing number of sources, broadly categorized as ‘filamentary threads’, whose origins remain a topic of ongoing debate (Rudnick et al. 2022). Prominent and extended filaments are often observed emanating from the cocoon structures of these jets, extending outward into the surrounding intracluster medium where the radio source is embedded (Ramatsoku et al. 2020; Rudnick et al. 2022; Koribalski et al. 2024). In addition, application of advanced edge-detection techniques have revealed a network of intricate filamentary structures within the cocoon as well (Condon et al. 2021; Velović et al. 2023). These internal filaments exhibit complex geometries, often intertwining or displaying irregular patterns, suggesting the interplay of various dynamical processes, such as turbulence, shear flows, and magnetic field compressions (Rudnick et al. 2021; Upreti et al. 2024; Weżgowiec et al. 2024).

Recently, observations of two GRGs with the MeerKAT radio telescope² have unveiled intricate internal (Cotton et al. 2020) and extended external filamentary structures (Cotton et al. 2025), emphasizing the need for a deeper investigation into their properties and origins in such large-scale jetted sources. These structures may offer critical insights into questions such as their potential role in non-thermal particle transport and re-acceleration. In this context, we have sought to identify correlations between the filamentary structures observed in our formulated emission maps (Figs. 4(A−C)) and the underlying three-dimensional dynamics responsible for their formation.

The emission maps for the five simulated cases¹ were analyzed for this purpose. A detailed examination of these maps (Figs. 4(A−C)) summarized that the bulk flow of the jet-beam undergoes substantial widening, particularly with respect to its injection radius (i.e., 1 kpc). This behavior is especially prominent in the low jet-power cases and occurs after the jet propagates a certain distance within the cocoon, well before reaching its full extent. This widening gives rise to filamentary extensions with different geometries, including tendril-shaped or finger-like structures, caldera-like formations characterized by emission valleys surrounded by enhanced emission ridges, and large vortex ring-like features spanning the lobe.

Despite such complex geometric appearances observed in the projected 2D maps, we found striking similarities to these emission layers in their three-dimensional counterparts data, when analyzing the variation of the Alfvénic speed within the jet-cocoon geometry. The Alfvén speed (v_Alfvén) is a key indicator of magnetized flow, defined as follows (given the sub-relativistic evolution of the cocoon material; Paper I),

$$\begin{array}{c}\hfill v_{\mathrm{Alfvén}}=\frac{B}{\sqrt{4\pi \rho }},\end{array}$$(22)

where B and ρ represent the magnetic field strength and the plasma density at each location within the cocoon, respectively. To ensure a fair comparison with the line-of-sight integrated emission map (where the line of sight is along the z-axis), we project the 3D distribution of v_Alfvén onto the x − y plane (by averaging and collapsing all values along the z-direction onto the z = 0 plane), as the dominant 3D layers of the jet-cocoon system will also reflect themselves in the projected map. This is shown in Fig. 10, where the projected v_Alfvén maps were subsequently overlaid with emission contours, with the contour levels carefully selected for each case to prominently highlight the filamentary structures noticed in the intensity maps (Figs. 4(A−C)). For instance, in the “GRG_lp_min” case, the contour levels in log(I_ν) (in units of erg s⁻¹ Hz⁻¹ kpc⁻² sr⁻¹) were set to [23.8, 24.45, 25.0, 25.25], while for the “GRG_lp_maj” case, the levels were chosen as [24.08, 24.45, 24.6, 24.9].

Fig. 10. Estimation of the strengths of various magnetic field values (in μG) for our five simulated scenarios, as shown on the horizontal axis. The magnetic fields included are the dynamical or intrinsic simulated B-field (Bdyn), the equipartition magnetic field (Beq), and the estimated CMB-equivalent magnetic field strength (BCMB) at a redshift of 0.05. The gray zone highlights the range of Beq values adopted for a sample of GRGs; details in Section 5.2.3.

The striking resemblance between the projected v_Alfvén maps and the emission contours, which trace the filamentary geometries, reveals two critical insights. First, the widening of the jet into the broader cocoon flow appears to be the primary driver of filament formation, with jet instabilities playing a comparatively minor role. Second, these filaments correspond to magnetized flows, effectively acting as magnetic threads within the cocoon. The first finding underscores the importance of turbulence-driven particle transport and acceleration within the cocoon, suggesting that this process could be a significant contributor to the observed radio emission in giant radio galaxies (e.g., see, Mukherjee et al. 2021; Kundu et al. 2022, for relevance in smaller radio galaxies). The second insight is particularly compelling, as it provides, to our knowledge, one of the most recent numerical indications of the magnetized nature of filaments (consistent with the observed B-field pattern of external filaments; Rudnick et al. 2022; Velović et al. 2023). The magnetization of these filaments, which is anticipated to manifest in polarization maps, has been observed and is discussed in detail in Section 6.

We investigated the typical lifetimes of these filaments within the cocoon and found them to be highly dynamic, with their geometry and appearance evolving on timescales significantly shorter than the lifetime of the GRG age. In Fig. 11, we present intensity maps from three simulations at an earlier epoch, ∼10 Myr before (top-row) the snapshots shown in Figs. 4(A-C) (bottom-row), corresponding to the data-saving interval of our simulations. The comparison between the two snapshots demonstrates that the most prominent filaments within the cocoon, formed directly by the channeling of the jet’s bulk flow, undergo significant topological changes on timescales shorter than 10 Myr. To constrain this variability further, we conducted a higher temporal-resolution simulation with data saved every 0.3 Myr intervals for “GRG_hp_min”¹ (a computationally intensive task, hence performed for one case only). Analysis of these better time-resolved emission maps reveals that prominent large-scale filamentary structures (distinctly identifiable through visual inspection) by their formation via jet flow branching and direct matter channeling, have estimated lifetimes of approximately 3–10 Myr. Intermediate-scale filaments, arising from jet’s bulk flow widening and decollimation, are comparatively shorter-lived, with lifetimes of ∼0.9–3 Myr. In contrast, smaller-scale filaments, identified qualitatively and less prominent in nature, are predominantly located farther from the bulk jet beam, contributing to the dynamic cocoon formation, which are highly transient, with lifetimes shorter than 0.3 Myr. To better constrain the survival age of these transient structures, higher-resolution data and significantly faster data-saving intervals are required.

Fig. 11. Amount of ambient medium matter entrained into the cocoon for the five simulated scenarios and with the evolution plotted against the linear growth of the cocoon (also mimicking the temporal growth).

6. Polarization measures at 1 GHz

In Fig. 12, we show the projected magnetic-field (B-field) lines (as inferred from the Stokes parameters; Section 3.2) overlaid on the emission maps at their evolved times, with fractional polarization values proportional to the line lengths. An important observation from the maps is that, despite of our jet injection model (i.e., a purely toroidal B-field injected with the jet flow, initially perpendicular to the x-axis; Section 2.1), the B-field lines eventually align with the flow direction (broadly along the x-axis), with added intricacies as the flow shows additional variations (e.g., bending, branching). This observed behavior aligns with the model proposed by Laing (1981), emphasizing the critical role of shock compression at the cocoon-ambient medium interface. Supplemented by the bulk flow, this process of compression increases shear and stress on the field lines, causing them to stretch and align tangentially to the lobe boundaries (see, e.g., Roberts et al. 2008).

Fig. 12. Projected color maps depicting the distribution of Alfvén velocity (vAlfvén; normalized by light-speed) on a logarithmic scale, overlaid with emission contours derived from the maps presented in Figs. 4(A−C). These combined maps enable the identification of a strong correlation between dominant emission features, such as the emission filaments, and their probable dynamical origins (e.g., the magnetized flow) across the five simulation cases, as indicated in the individual sub-figures. The 3D maps of Alfvén speed is projected onto the x − y plane to allow direct comparison with the emission maps with their line-of-sight observation aligned along the z-axis.

In the low jet power cases in Fig. 12, a notable feature of the frontal lobes is the presence of lobe-sized vortices, where the bulk flow bends (counterclockwise in “GRG_lp_min” and clockwise in “GRG_lp_maj”)¹. The projected B-field morphology closely follows the flow topology in these regions, reinforcing the conclusion that the jet material drags the magnetic field lines along with it as the flow progresses (e.g., Baidoo et al. 2023). The FP values remain high along the flows until vortex motion becomes dominant, inducing turbulence in both the flow and the magnetic field. This turbulence reduces the FP values in the projection, as reflected in the shorter line lengths (Fig. 12).

To correlate this observation further with the dynamical B-field distribution, we present 2D slices of magnetic fields from the 3D simulated data, in Fig. 13. In the figure, the contours represent the poloidal components of the B-field (B_p), while the colormap illustrates the toroidal component (B_ϕ). The slices are chosen from different planes along the z-axis to showcase individual layer’s contribution to the sky projected polarization map (note: the polarization maps in Fig. 12 considers contribution of all cocoon layers, and so, the sliced dynamical maps can be treated only as a simplified interpretation), as indicated in the figure. In the low jet power cases, we observe that the dynamical B-field lines predominantly align along the flow direction. However, as smaller-scale vortices and turbulence layers dominate, both the strength and alignment of the field lines are significantly impacted (transforming into more chaotic arrangements).

Fig. 13. Emission maps for three simulations (labeled in the sub-figures) at times ∼10 Myr earlier than those presented in Figs. 4(A-C), reiterated in the bottom row. A comparison between these temporal snapshots underscores the dynamic and transient nature of filamentary structures within the cocoon, while also aiding in constraining the typical lifetimes of these multi-scaled internal features. For further details, refer to Section 5.3.

A closer examination of the polarization maps also reveals the intricate filamentary structures within the lobes, as discussed in detail in Section 5.3. The projected B-field lines are observed to closely follow these filamentary trails, reaffirming that such internal filaments represent magnetic flows, resembling intertwined strands (e.g., Sebokolodi et al. 2020; Andati et al. 2024).

In the high jet power cases, a consistent pattern is observed where the B-field lines align strongly with the jet spine, extending nearly up to the jet head. However, as the distance increases, the stretching of the field lines gradually weakens the magnetic field. By comparing Figs. 12 and 13, it becomes clear that the B-field strength remains significant as long as the jet maintains its collimation. As the distance increases from jet injection, the backflow material becomes more dominant, carrying away part of the magnetic field. The generated backflow material also shows alignment of the field lines along the flow, but their turbulent nature becomes increasingly significant (similar to, Mukherjee et al. 2020), particularly beyond the jet sheath structure, as prominently seen in the “GRG_hp_edge” case in both Figs. 12 and 13.

For the winged sources (the “maj” simulations), the magnetic field lines in the wings tend to align with the flow (fluid motion around the y-axis). The flow within these secondary lobes originate from the channeled backflowing plasma from the active lobe into the wing cavity, further reinforced by the backflowing matter from the opposing jet arm (discussed also in Section 5.1.2). The B-field lines within the wings, observed to approximately aligning towards the y-axis in the projected polarization map, indicate the deflection and channelling of backflowing plasma into the wing. This configuration is particularly evident in the sliced dynamic maps in Fig. 13. At the head location of the wings, the compression of the B-field lines creates a bow-like arrangements.

Building on the above discussion, which highlights the development of turbulence in the cocoon and its significant impact (diminishing) on FP estimates, we further analyzed the variation of FP values with the increasing length of the GRGs. This is presented in Fig. 14, which shows a general decrease in the median FP values with growing length, suggesting that GRGs become increasingly turbulent in their cocoons as they evolve. Initially, the FP values are relatively high, although still below the theoretical upper limit of ∼71%, indicating that the jets carry an ordered magnetic field due to their high thrust (see also Fig. 13). However, as the GRG grows, the FP values gradually decline, showing a potential saturation at larger lengths, stabilizing around ∼40% (see Willis & Strom 1978; Cotton et al. 2020, for relevance of these FP values in the context of GRGs).

Fig. 14. Intensity variation maps (as background) for five simulated scenarios at their evolved stages presented as color maps, overlaid with projected magnetic field distributions. The line lengths indicate proportionality to the corresponding FP values at each point. Median FP values are noted for ‘GRG_lp_min’ (41.9%), ‘GRG_hp_min’ (45.8%), ‘GRG_lp_maj’ (38.4%), ‘GRG_hp_maj’ (39.2%), and ‘GRG_hp_edge’ (35.1%), and are annotated in each subfigure alongside a scale bar for reference.

Fig. 15. The x − y slices of the magnetic field configuration, taken at different z-plane positions. In each slice, the contours represent the poloidal magnetic field component (Bp), while the toroidal component (Bϕ) is shown using a colormap (in μG). Together, these components provide insight into (a) the dynamic characteristic of magnetic field structures within the jet–cocoon system and (b) their individual contribution to the features seen in polarization maps (Fig. 12). Additionally, the variation across z demonstrates the non-axisymmetric bulk motion of the jet beam in three dimensions, emphasizing the limitations of employing a 2D model. To aid comparison of the field configurations within each group, red dashed lines have been added as visual guides to the reader.

Fig. 16. Evolution of the median FP values in the simulated emission maps for the five situations, plotted against their growing length (mimicking also their aging), illustrating a gradual decrease in FP values followed by a near-settling phase. The FP values are calculated at 1 GHz at a spatial resolution of 1 kpc, corresponding to an angular resolution of 0.9 arcsecond for a redshift of 0.05.

Whether this saturation indicates a settling phase in the evolution of GRGs or represents a brief transitional phase requires further investigation. Future studies, incorporating higher-resolution simulations and exploring the correlation length of the turbulent magnetic field distribution, will be essential to address this question. In that case, the effect of resolution on FP estimations can also be explored by convolving high-resolution maps to the beam sizes of different telescopes.

7. Summary

The numerical modeling of giant radio galaxies and their long-term evolution has recently advanced, with the first 3D relativistic magnetohydrodynamic simulations presented in Paper I (Giri et al. 2025). Paper I focused on the intrinsic dynamical properties of GRGs, examining how their thermodynamic and magnetic characteristics evolve over age through a series of five simulations with varying jet powers and differing levels of resistance from the ambient medium to jet propagation. However, it is crucial to validate these results against observable quantities, particularly those measurable by contemporary radio telescopes. Such observations are key to addressing enduring questions in GRG physics, including the persistence of their emission despite the diminishing influence of several physical processes over time. Using a standard emission treatment to transform the 3D data of dynamical entities into sky-projected continuum and polarization maps at 1 GHz, we uncover the role of a few of these mechanisms intricately linked to the evolution of GRGs, as detailed below:

1. In studying the matter transport processes within the GRG cocoon, we found that decollimation and bending are characteristic features of low-power jets, which lead to lobe formation. In contrast, high-power jets maintain their collimation, resulting in arrowhead- or nose-cone-shaped structures. This morphological distinction is additionally influenced by tri-axiality of the surrounding environment, giving rise to features such as winged sources and asymmetric cocoon extents. These findings thus suggest that GRG morphology can serve as a valuable proxy for disentangling the contributions of jet power and environmental factors.

2. A prominent feature observable from the emission maps is the formation of large-scale vortex motion in low-powered GRGs and the branching or bifurcation of jet flows in high-powered GRGs. Both phenomena lead to decollimation and the generation of filamentary structures with varying geometries, such as tendril- or finger-like extensions and caldera-like formations, characterized by emission valleys surrounded by enhanced emission ridges, beside the formation of other smaller-scaled non-linear topologies.

3. In the high-powered jet cases, prominent enhanced emission region often appears where the jet begins to lose its collimation, located approximately midway along the total cocoon extent. Despite identifying (weak) warm spots at the ends of two out of three high-powered simulations, the presence of such mid-cocoon enhanced emission zones could lead to misclassifications, suggesting that the GRG is restarting in nature (similar to concerns highlighted in recent numerical studies of smaller radio galaxies). In such scenarios, we propose studying the lateral intensity variations in GRGs (perpendicular to the bulk flow). An inactive or low-activity cocoon would exhibit a gradual decline in lateral emission, whereas active, high-powered GRGs would display abrupt variations in lateral intensity measures.

4. A quantitative analysis of the five synthetic continuum maps reveals a set of key findings, as follows. First, applying a simple jet power–radio power conversion relationship to a statistical sample of GRGs can be misleading, as it varies with the parameters shaping individual GRG formation. Second, different estimates of magnetic field values show significant deviations from equipartition assumptions. Inverse Compton cooling by CMB photons appears to be a dominant cooling mechanism for GRGs, even at low redshifts. Lastly, matter entrainment has a diminishing role in sustaining GRG emission; low-powered jets entrain more environmental matter due to lobe formation, while high-powered jets retain collimation and thrust, limiting entrainment in the jet spine but heavily impacting the cocoon’s persistence.

5. Striking similarities between the filamentary structures within the cocoon (as detected in the emission maps) and the Alfvénic velocity map, supported by aligned projected B-field lines and higher fractional polarization values along such threads in the polarization maps, confirm the magnetized nature of these internal filaments. These prominent filamentary structures within the cocoon, often identified in observed maps through simple edge-detection techniques, persisting for several Myr (dynamical age).

6. Our findings suggest that the reduction in fractional polarization values is linked to the development of turbulence within the cocoon of GRGs. As the jetted sources age and grow in length, the fractional polarization exhibits a diminishing trend, indicating a transition of the magnetic fields into more chaotic configurations. Such a phenomenon can further be modeled with varying resolutions (or beam sizes) and frequencies to determine the correlation length of the underlying turbulence. Regardless of our jet injection models, the projected B-field lines align nearly parallel to the bulk flow of the jet-cocoon system.

While this study serves as a foundational step to correlate the simulated dynamic and synthetic emission characteristics of GRGs, it provides a basis upon which future advancements and refinements can be achieved. First, a broader parameter set needs to be explored numerically to model the observed characteristics of GRGs in a statistical sample, including their prevalence as FR II GRGs with prominent hotspots at the jet head. This would involve exploring scenarios with a lower jet-to-environment density contrast, such as incorporating heavier jets, simulating jet propagation in a denser ambient medium, or injecting jets with higher power. Furthermore, examining diminishing effects, including matter entrainment in high-resolution simulation grids, its dependence on magnetic field strength, and the incorporation of radiative cooling into simulations, is crucial for understanding the persistence of emission during the longer-term evolution of GRGs. Several of these aspects were highlighted in Paper I and are planned for further investigation in our future work. Additionally, emphasis can be placed on translating these dynamical simulations into pragmatic radio maps, using more accurate emission conversion methods, such as the Lagrangian-Eulerian approach (e.g., Vaidya et al. 2018), which accounts for particle re-acceleration processes and the impact of turbulence on polarization behavior.

Acknowledgments

GG and KT acknowledge support from the South African Department of Science and Innovation’s National Research Foundation under the ISARP RADIOMAP Joint Research Scheme (DSI-NRF Grant Number 150551). CF acknowledges support by the German Research Foundation DFG via the research unit FORS 5195. JB acknowledges the support from the Department of Physics and Electronics, Christ University, Bangalore. JD and GG acknowledges partial research-travel support by the National Research Foundation of South Africa (Ref Number CSUR240426216203). GG acknowledges the travel support from Gianluigi Bodo and Paola Rossi (INAF Torino) through Ob. Fu. 1.05.01.89.01 and Ricerca di base gruppo extragalattico, where discussions played a pivotal role in shaping this project. RPD acknowledges funding by the South African Research Chairs Initiative of the Department of Science, Technology, and Innovation and National Research Foundation (Grant ID: 77948). We acknowledge the use of the ilifu cloud computing facility – https://www.ilifu.ac.za/, a partnership between the University of Cape Town, the University of the Western Cape, Stellenbosch University, Sol Plaatje University, the Cape Peninsula University of Technology and the South African Radio Astronomy Observatory. The ilifu facility is supported by contributions from the InterUniversity Institute for Data Intensive Astronomy (IDIA – a partnership between the University of Cape Town, the University of Pretoria and the University of the Western Cape), the Computational Biology division at UCT and the Data Intensive Research Initiative of South Africa (DIRISA). The authors acknowledge the Centre for High Performance Computing (CHPC), South Africa, for providing computational resources to this research project (https://www.chpc.ac.za/).

References

All Figures

Fig. 1. Schematic representation of the jet-environment interaction framework developed to investigate the influence of factors on the development of GRGs. The setup simulates a one-sided relativistic jet propagating up to ∼700 kpc through a stratified medium representing a galaxy-group environment, capturing the varying influence of the environment on the jet’s motion. We also used two different jet powers to mimic Fanaroff-Riley type I and II jets, allowing us to examine the role of the jet’s primary thrust in the development of such radio galaxies. This parameter space yields a total of five simulation scenarios, discussed in detail in Paper I and also concisely in Section 2.1.

In the text

Fig. 2. Sliced velocity maps (x − y plane, z = 0) for the resultant GRG morphologies from five simulation cases (labels indicated in each plot) illustrating the underlying matter-transport mechanisms responsible for these structures. These maps are displayed at dynamical ages of ∼167 Myr (GRG_lp_min), 196 Myr (GRG_lp_maj), 49 Myr (GRG_hp_edge), 137 Myr (GRG_hp_maj), and 69 Myr (GRG_hp_min), as also highlighted in Paper I. The colormap represents the x-component of velocity, capturing the bulk-flow (in blue) and back-flow (in red) dynamics, while the contours depict the y-component of velocity, with solid lines indicating 0.005c and dashed lines indicating −0.005c, to highlight the behavior of secondary structures forming within the domain. To analyze the impact of jet power on morphology and the influence of the tri-axial environment on the resulting structure, comparisons are made column-wise and row-wise, respectively.

In the text

Fig. 3. Initial density configuration of the ambient environment, illustrating a 10° rotation of the ambient medium’s principal axis, in order to introduce an asymmetry to the jet-flow, while the jet injection axis remains fixed along the negative x-direction. The density distribution is shown using colormap, with the corresponding log ρ[amu/cc] values annotated along the contour lines. This represents one of the simulation scenarios employed, as shown in Fig. 1.

In the text

Fig. 4. A. Intensity maps (log Iν) at 1 GHz radio frequency for two simulation cases: a low-power jet (“GRG_lp_min”; top) and a high-power jet (“GRG_hp_min”; bottom). Both jets are propagating along the minor axis of the ambient medium. While the overall structures resemble the dynamical maps from Paper I, the emission maps offer deeper insights, highlighting crucial aspects of the observability of these simulated giant radio structures at their evolved stages (∼167 Myr: top; ∼69 Myr: bottom).

In the text

Fig. 5. B. Continued. Intensity maps (log Iν) at 1 GHz radio frequency depicting the X-shaped radio galaxy in their giant stages, resulting from jet propagation along the direction of highest jet frustration (i.e., along the major axis direction). The low jet power case (“GRG_lp_maj”; top) and the high jet power case (“GRG_hp_maj”; bottom) exhibit distinct decollimation patterns, with both cases generating internal regions of emission valleys and enhanced emission ridges. The wings are luminous, indicating their active evolutionary stages (∼196 Myr: top; ∼137 Myr: bottom).

In the text

Fig. 6. C. Continued. Intensity map (log Iν) at 1 GHz radio frequency revealing an intriguing case where the jet flows with minimal hindrance due to its propagation along the edges of the surrounding environment (“GRG_hp_edge”). The faster propagation causes increased shear between the cocoon and the environment as well as between the jet spine and sheath layers and results in the generation of instabilities, leading to the formation of pronounced turbulent structures at 49 Myr.

In the text

Fig. 7. Dynamic maps illustrating the jet spine and head structures, elucidating the absence of prominent hotspots in our synthetic emission maps for high-power jet simulations. The top panel presents a 3D volume-rendered visualization of the jet tracer, highlighting the jet beam – jet head morphology at their corresponding evolutionary ages. The bottom panel shows 2D slices focused on the jet head: the top row maps the magnetic field magnitude, while the bottom row maps the pressure distribution. These maps reveal that although the jet sustains its thrust up to the terminal shock, the magnetic field strength drops sharply at this region – largely advected away by the backflowing plasma – thereby suppressing the formation of bright hotspots.

In the text

	Fig. 8. Lateral intensity variation measured at the maximum cocoon width for five simulation scenarios (Fig. 4(A−C)), providing a quantitative assessment of intensity distribution and identifying active zones within the cocoon. Subplots correspond to different simulation runs categorized by ‘min’, ‘maj’, and ‘edge’ cases, enabling a comparison of jet power on cocoon activity and evolution.
In the text

Fig. 9. Radio power at 1 GHz of the (one-sided) cocoon versus (one-sided) cocoon length, presented at the evolved stages (dynamical ages attached) of the simulated GRGs. Two jet power values used in the simulations are also highlighted, in blue and red. The evolutionary track of the simulated GRGs in the Power-Distance (P-D) diagram is also shown from the point they exceed 350 kpc, qualifying as GRGs in their total extent. For reference, we include measurements of three observed GRGs from Charlton et al. (2024) in similar observing band, representing a typical one-half estimation. We note that the above calculation does not explicitly account for the cooling effects.

In the text

Fig. 10. Estimation of the strengths of various magnetic field values (in μG) for our five simulated scenarios, as shown on the horizontal axis. The magnetic fields included are the dynamical or intrinsic simulated B-field (Bdyn), the equipartition magnetic field (Beq), and the estimated CMB-equivalent magnetic field strength (BCMB) at a redshift of 0.05. The gray zone highlights the range of Beq values adopted for a sample of GRGs; details in Section 5.2.3.

In the text

	Fig. 11. Amount of ambient medium matter entrained into the cocoon for the five simulated scenarios and with the evolution plotted against the linear growth of the cocoon (also mimicking the temporal growth).
In the text

Fig. 12. Projected color maps depicting the distribution of Alfvén velocity (vAlfvén; normalized by light-speed) on a logarithmic scale, overlaid with emission contours derived from the maps presented in Figs. 4(A−C). These combined maps enable the identification of a strong correlation between dominant emission features, such as the emission filaments, and their probable dynamical origins (e.g., the magnetized flow) across the five simulation cases, as indicated in the individual sub-figures. The 3D maps of Alfvén speed is projected onto the x − y plane to allow direct comparison with the emission maps with their line-of-sight observation aligned along the z-axis.

In the text

Fig. 13. Emission maps for three simulations (labeled in the sub-figures) at times ∼10 Myr earlier than those presented in Figs. 4(A-C), reiterated in the bottom row. A comparison between these temporal snapshots underscores the dynamic and transient nature of filamentary structures within the cocoon, while also aiding in constraining the typical lifetimes of these multi-scaled internal features. For further details, refer to Section 5.3.

In the text

Fig. 14. Intensity variation maps (as background) for five simulated scenarios at their evolved stages presented as color maps, overlaid with projected magnetic field distributions. The line lengths indicate proportionality to the corresponding FP values at each point. Median FP values are noted for ‘GRG_lp_min’ (41.9%), ‘GRG_hp_min’ (45.8%), ‘GRG_lp_maj’ (38.4%), ‘GRG_hp_maj’ (39.2%), and ‘GRG_hp_edge’ (35.1%), and are annotated in each subfigure alongside a scale bar for reference.

In the text

Fig. 15. The x − y slices of the magnetic field configuration, taken at different z-plane positions. In each slice, the contours represent the poloidal magnetic field component (Bp), while the toroidal component (Bϕ) is shown using a colormap (in μG). Together, these components provide insight into (a) the dynamic characteristic of magnetic field structures within the jet–cocoon system and (b) their individual contribution to the features seen in polarization maps (Fig. 12). Additionally, the variation across z demonstrates the non-axisymmetric bulk motion of the jet beam in three dimensions, emphasizing the limitations of employing a 2D model. To aid comparison of the field configurations within each group, red dashed lines have been added as visual guides to the reader.

In the text

	Fig. 16. Evolution of the median FP values in the simulated emission maps for the five situations, plotted against their growing length (mimicking also their aging), illustrating a gradual decrease in FP values followed by a near-settling phase. The FP values are calculated at 1 GHz at a spatial resolution of 1 kpc, corresponding to an angular resolution of 0.9 arcsecond for a redshift of 0.05.
In the text