© The Authors 2025

1. Introduction

A limited fraction of active galactic nuclei (AGN) with accreting supermassive black holes (SMBHs; with mass ≥10^6.5 M_⊙; Baldi et al. 2021) at their centers exhibit the ejection of highly collimated plasma flows with bulk Lorentz factors reaching seven or higher (Kellermann et al. 1989; Merloni & Heinz 2007; Padovani et al. 2011; Saikia et al. 2016). Therefore, considerable focus has been placed on understanding the extent to which such high-speed jets can propagate given the significant jet hindrance caused by the denser ambient medium through which the jets travel. Studies have revealed a wide range of jet lengths (cf. Saikia 2022, for a recent review), with some ending in the sub-kiloparsec scales (compact radio cources, or CRSs; Baldi 2023), some extending to galactic scales (galaxy-scale jets, or GSJs; Webster et al. 2021; O’Dea & Saikia 2021), some propagating to hundreds of kiloparsecs (extended kiloparsec-scale jets; Hardcastle & Croston 2020), and a small fraction evolving on megaparsec scales (giant radio galaxies, or GRGs; Dabhade et al. 2023). When determining the availability of such sources, large surveys have reported that the longer the jet extends, the rarer these sources become (Hardcastle et al. 2016; Ishwara-Chandra et al. 2020; Gürkan et al. 2022). Being rare, these extended radio sources attract significant attention because of their extreme physical conditions, encouraging modeling efforts on various aspects such as the diverse ambient environments they propagate through and different jet propagation models that help them cover such distances. An important aspect to consider when studying extended radio jets is that they can also provide energetic feedback to the surrounding medium (Fabian 2012; McNamara & Nulsen 2012; Heckman & Best 2023). For an example, if kinetic energy of the outflow is efficiently deposited in the environment, it can significantly affect the star formation rate (Nesvadba et al. 2021; Drevet Mulard et al. 2023; Roy et al. 2024).

This paper focuses on understanding GRGs, which are recognized as sources with jetted structures extending to at least 700 kpc. Initially, the length cutoff for GRGs was fixed at ≥1 Mpc when the value of the Hubble constant was accepted to be H₀ = 50 km s⁻¹ Mpc⁻¹ (Willis et al. 1974). Many authors of recent works have tried to reset the length cutoff for GRGs again as ≥1 Mpc with the present-day estimation of the Hubble constant of H₀ = 70 km s⁻¹ Mpc⁻¹ (Andernach et al. 2021; Hardcastle et al. 2016; Brüggen et al. 2021). Regardless of the length threshold, the origin of such extended sources appears to be ambiguous, which is supported by the fact that they constitute a relatively small fraction of radio-loud AGNs (Dabhade et al. 2020a,b) – there are only a handful of sources that have linear extents exceeding 3 Mpc (Andernach et al. 2021). Giant radio sources reaching a total extent of 5 Mpc are extremely rare (Willis et al. 1974; Machalski et al. 2008; Oei et al. 2022), raising the question of whether there is a physical size limit to how far a jet can propagate. Notably, the recent discovery of an AGN jet spanning approximately 7 Mpc has introduced additional challenges to this topic (Oei et al. 2024a).

With advancements in the sensitivity and resolution of radio telescopes, the discovery of GRGs has increased, yet they still represent only a relatively small fraction of the extended-jetted population (Andernach et al. 2021; Bhukta et al. 2022; Oei et al. 2023a; Mostert et al. 2024). With the advancements in source statistics, the hypothesized formation mechanisms generating GRGs have been revised to incorporate new data and insights. Currently, two models are being considered to explain the origin of GRGs: one attributes their growth to propagation in a less dense ambient medium, while the other links their formation to complex nuclear activity (for a review, see, Dabhade et al. 2023).

The former model suggests that GRG growth occurs in the low-density warm-hot intergalactic medium typically found in the galaxy filaments of the Cosmic Web (Davé et al. 2001; Oei et al. 2024b) or that GRGs evolve near the edge of a galaxy group or cluster medium. Significant support for this model has been found in statistical studies, where available cluster and filament catalogs have been crossmatched or where polarization measures have been conducted (Stuardi et al. 2020; Oei et al. 2023a; Simonte et al. 2024), as well as in focused source studies that use the galaxy distribution around the hosts of such sources to infer their large-scale environments (Subrahmanyan et al. 2008; Safouris et al. 2009; Pirya et al. 2012; Malarecki et al. 2015; Oei et al. 2022). Additional evidence comes from the observation that GRGs also exhibit off-axis lobe distortions, which are linked to the presence of an underlying asymmetric ambient medium (Subrahmanyan et al. 1996; Malarecki et al. 2013). Furthermore, the lengths of the two lobes are typically dissimilar, and the ratio of bigger to smaller lobe lengths are marginally greater than those of their smaller radio galaxy (SRG) counterparts (< 0.7 Mpc), which cannot be fully explained by orientation effects alone (Saripalli et al. 1994; Ishwara-Chandra & Saikia 1999; Schoenmakers et al. 2000).

This model is, however, not free from caveats, as studies have found that ∼10% of GRGs evolve within a group or cluster medium, with many of them acting as the brightest cluster galaxy (Komberg & Pashchenko 2009; Dabhade et al. 2020a; Andernach et al. 2021; Sankhyayan & Dabhade 2024). The studies by Lan & Prochaska (2021) and Simonte et al. (2022) found no preference for ambient environment factors, such as the location of the host galaxy with respect to the surrounding large-scale structure or the properties of nearby galaxies that would favor the growth of GRGs over SRGs. Moreover, one needs to also confirm that the galaxy distribution used as a proxy for the gas distribution is in general true; an exception to this may be seen in Palma et al. (2000).

The latter model, which involves activity near the host AGN, includes possibilities such as a high-power jet ejection, a subsequently longer jet ejection time, or restarting activity that helps the lobes grow at a considerable rate. The possibility that recurrent activity of AGN contributes to the growth of GRGs seems unlikely to be a common occurrence given that the availability of such cases among GRGs is low (∼5%, Dabhade et al. 2020a,b). The fraction of extended jets exhibiting restarting activity is low in general (13–15%; Jurlin et al. 2020). A longer jet evolution time appears to be an important aspect and has been considered in several studies, which have indicated evidence of a greater jet evolution age for GRGs compared to their smaller counterparts (e.g., Jamrozy et al. 2008; Komberg & Pashchenko 2009; Machalski 2011; Dabhade et al. 2023). However, several studies have also estimated relatively younger (spectral) ages for GRGs (Mack et al. 1998; Konar et al. 2004) comparable to SRGs (Harwood et al. 2017; Turner et al. 2018).

The distribution of radio power (and thus jet power) of GRGs indicates that they predominantly lie above the classical Fanaroff & Riley (1974) I/II boundary (Dabhade et al. 2020a; Andernach et al. 2021; Simonte et al. 2024). This is often cited as the reason why the jets can overcome jet frustration (hindrance to jet motion) imposed by the ambient medium. However, if this is the case, it must also be explained why only ∼10% of FR II radio jets evolve into GRGs (Komberg & Pashchenko 2009). Despite the lower occurrence of FR I GRGs, they also appear to overcome the opposing force of the ambient environment to the jet flow to reach such distances. Notably, the majority of FR I sources reside in dynamically more evolved galaxy clusters or in rich galaxy groups (Zirbel 1997; Mingo et al. 2019). The discovery of GRGs exceeding 1 Mpc in length originating from spiral or disk-type host galaxies has further complicated efforts to disentangle the connection between GRG origin and central AGN activity (Hota et al. 2011; Bagchi et al. 2014; Clarke et al. 2017; Oei et al. 2023b). Clues from host AGN black hole masses have also been inconclusive, as an extended GRG of ∼5 Mpc size (Alcyoneus) hosts a black hole of ∼4 × 10⁸ M_⊙, while a much more massive black hole of mass ∼5 × 10¹⁰ M_⊙ (J103129.54+502959.1) produces a jet with a length of 0.81 Mpc (Oei et al. 2022, 2023a). As also reported by Kuźmicz & Jamrozy (2012), the median values in black hole masses in giant radio quasars (GRQs) are not considerably different from those of smaller radio quasars. On the other hand, when compared to GRGs, GRQs are reported to have a greater black hole mass and a higher Eddington ratio (Mahato et al. 2022).

A numerical approach to test these models and their caveats has thus become essential in order to enable examination of jet evolution phases across a broad spectrum of ambient configurations and jet powers. This is the primary focus of this work, as we aim to determine whether GRGs are similar to their smaller counterparts or if they represent a distinct class of sources. By formulating different jet ambient medium settings, we aim to explore the (thermo)dynamical states of the evolving intertwined structures. Given the length scales involved for GRGs, the simulations are computationally intensive, and to our knowledge, this work is the first of its kind. We intend to simulate sources with lengths extending to 1 Mpc and beyond, as it is expected that the larger the source, the more prominently the growth factors contributing to its expansion will appear (Oei et al. 2022).

This paper is structured as follows: In Sect. 2, we discuss the formulation of different simulation settings to test GRG formation models. Section 3 elaborates on the results obtained from the simulations, including the morphological distinctions of the grown structures (Sect. 3.1), their temporal evolution (Sect. 3.2), and dynamical characteristics (Sect. 3.3), and it includes a discussion on their relevance to observations. We then summarize our work in Sect. 4.

2. Numerical setup

Considering the crucial importance of modeling magnetic fields in the radio lobes of GRGs (Jamrozy et al. 2008; Machalski et al. 2009), as well as the presence of relativistic effects in the jets of GRGs (Saripalli et al. 1997; Giovannini et al. 1999; Hocuk & Barthel 2010), we formulated our setup using a relativistic magnetohydrodynamical framework with the PLUTO astrophysical code (Mignone et al. 2007). The relativistic magnetohydrodynamical equations to be solved are defined as follows:

$$\begin{array}{cc}\hfill {\partial }_{\kappa }(\rho u^{\kappa })& =0\hfill \\ \hfill {\partial }_{\kappa }(w{u}^{\kappa }{u}^{\delta }-b^{\kappa }{b}^{\delta }+p{g}^{\kappa \delta })& =0\hfill \\ \hfill {\partial }_{\kappa }(u^{\kappa }{b}^{\delta }-u^{\delta }{b}^{\kappa })& =0,\hfill \end{array}$$(1)

where κ, δ =(0, 1, 2, 3) (Mignone & Bodo 2006; Mignone et al. 2009). The symbols in Eq. (1) represent the following: ρ is the rest mass density; u^κ is the four velocity, Γ(1, v): v is the three velocity represented through Γ, the jet Lorentz factor; and b^κ is the covariant magnetic field prescribed through the lab field B as (b⁰, B/Γ + b⁰v). The total enthalpy and total pressure of the system are defined as w and p, respectively. The metric g^κδ takes the form diag(−1, 1, 1, 1).

2.1. Configurations of the ambient medium

Substantial work aimed at understanding the ambient medium in which GRGs thrive indicates that these giant radio jets prefer to reside in cosmic filaments, smaller groups, and low-mass clusters, the density of which can be modeled using the general King’s beta-profile (Machalski et al. 2009; Oei et al. 2023b). This tri-axial density profile (Cavaliere & Fusco-Femiano 1976) can be expressed (in the Cartesian system) as

$$\begin{array}{c}\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}.\end{array}$$(2)

Eq. (2) is characterized by ρ₀, which indicates the core density of the medium, and β, which dictates the density gradient. By setting ρ₀ to 0.001 amu/cc (∼200 times the critical cosmic density of the local universe) and β to 0.55, the particle density ranges from 10⁻³ cm⁻³ at the core to 10⁻⁶ cm⁻³ at the outskirts (∼700 kpc), consistent with typical environmental density estimated around GRG lobes (Mack et al. 1998; Stuardi et al. 2020). Given such a profile, the total baryonic mass evaluated within one megaparsec is ∼3.3  ×  10¹² M_⊙, and including the dark-matter the total mass approaches 1.9  ×  10¹³ M_⊙ with a virial radius (Gaspari et al. 2020) of 747 kpc, representing a poor galaxy group or warm-hot intergalactic medium (Sun et al. 2009; Oei et al. 2023b).

To break any symmetry that may arise while the jet propagates along the symmetry axes of the 3D configuration, we rotated the ambient medium by 10° in the x − y plane, keeping the z-axis as the rotation axis. Therefore, x′ and y′ in Eq. (1) are defined as follows:

$$\begin{array}{c}\hfill x^{\prime }=x\cos {10}^{°}-y\sin {10}^{°},y^{\prime }=x\sin {10}^{°}+y\cos {10}^{°}.\end{array}$$(3)

The values of a, b, and c control the tri-axial spheroid shape of the medium and represent the effective core radii, whereas (x₀,  y₀,  z₀) represents the center of the ambient medium. To examine jet growth across a range of ambient conditions, we selected three different configurations (given that the jet always flows along the −x direction) by varying these values as follows:

1. Jet flow along the minor axis (“min”): Values of b = 66 kpc, while a and c = 33 kpc to reduce jet frustration.

2. Jet flow along the major axis (“maj”): Values of a = 66 kpc, while b and c = 33 kpc to increase jet frustration.

3. Jet flow along the edges (“edge”): Values of a,  b, and c = 66 kpc to create a spherical environment and then initiate jet propagation from the edge (600 kpc above host’s center).

The aforementioned configurations are illustrated in Fig. 1 through a schematic diagram and 2D slices of the initial ambient medium setup. Such choices of core radii, β, and ρ₀ for a galaxy group or cluster medium are typically adopted in analytical formulations of these environments (Hardcastle 2018; Musoke et al. 2020; Giri et al. 2023). The tri-axial nature of the ambient medium is also not uncommon and has been observed in statistical X-ray studies of (unrelaxed) cluster sources (e.g., Hodges-Kluck et al. 2010; Brüggen et al. 2021), as well as in many individual X-ray studies of groups and cluster radio sources (as inferred through preliminary assessment) (Hodges-Kluck & Reynolds 2011; Bogdán et al. 2014; Randall et al. 2015; Ramatsoku et al. 2020; Pandge et al. 2022). Additionally, the 3-dimensional shape and dynamics of dark matter halos around galaxies and galaxy groups or clusters is a frontier topic of ongoing research, but numerical simulations (Chua et al. 2019; Emami et al. 2021; Baptista et al. 2023) and limited observations (e.g., Buote et al. 2002) suggest that they are not perfectly spherical but may assume a triaxial spheroidal geometry, with long axis preferably oriented parallel to longer axis of the central galaxy or group. This is consequence of the fact that the gravitational force that shapes dark matter halos are not perfectly symmetrical, which governs their 3D equilibrium figures in a rotating frame, much like a Jacobi ellipsoid (Chandrasekhar 1969). Significantly, for the extended warm-hot (temperature ≲1 keV) circum-galactic halo surrounding the galaxy 2MASX J23453268-0449256 (a spiral), which hosts a (∼1.6) Mpc scale GRG (Bagchi et al. 2014), the shape of the X-ray halo was found to be ellipsoidal rather than spherical (Mirakhor et al. 2021; Bagchi et al. 2024) with jet axis at a steep angle to the major axis of the halo.

Fig. 1. Concise overview of the simulations conducted in this work. Left column: Brief description of the jet ambient medium configurations used to explore a broad spectrum of jet environment parameters. Middle column: Ambient configuration in relation to the Cartesian axes, with a 3D box representing the schematic simulation domain. Right column:x − y slice of density (at z = 0) illustrating the initial variation of (ρ/ρ0) as a color map, with overplotted pressure (P/P0) contours at values of [4.5,  2.2,  1.2,  0.67]×10−7, outward. Jet injection zones are marked with black and white dots in the density slices. The values of ρ0 and P0 are 0.001 amu/cc and 1.5 × 10−6 dyn/cm2 respectively, and lengths are represented in units of 100 kpc. To break the symmetry with respect to the jet flow, the ambient medium has been slightly rotated.

The presence of GRGs within rich groups, clusters, and even super clusters appears to be more prevalent than previously assumed, with many serving as the brightest galaxies within their environment (Sankhyayan & Dabhade 2024), which influences the “min” and “maj” cases. In the case of “edge”, the jet is propagating at the periphery of the ambient environment, considering the virial radius of less evolved groups and poor clusters to be varying around ∼800 kpc (Hota et al. 2011; Stuardi et al. 2020; Oei et al. 2023b).

We assigned the pressure distribution P(x′,y′,z) of the ambient medium following the density distribution (correlated with the adiabatic gas constant γ), which then generates an isothermal atmosphere of temperature 1.6 keV (Oei et al. 2023b). A representation of pressure values through contours has also been highlighted in Fig. 1, indicating the drastic, slow, and relatively constant changes in pressure gradient along the minor, major, and edge configurations, respectively. To maintain the static equilibrium of such a medium, a gravitational profile has been assigned, derived from the hydrostatic equilibrium equation as

$$\begin{array}{c}\hfill \mathrm{\nabla }P=\rho \mathbf{g},\end{array}$$(4)

where g is the acceleration due to gravity. The three components of g are represented as follows:

$$\begin{array}{cc}\hfill g_{x^{\prime }}& =-\frac{3\beta }{\gamma }\left(\frac{x^{\prime }-x_0}{a^2}\right){(1+\frac{{(x^{\prime }-x_0)}^2}{a^2}+\frac{{(y^{\prime }-y_0)}^2}{b^2}+\frac{{(z-z_0)}^2}{c^2})}^{-1}\hfill \\ \hfill g_{y^{\prime }}& =-\frac{3\beta }{\gamma }\left(\frac{y^{\prime }-y_0}{b^2}\right){(1+\frac{{(x^{\prime }-x_0)}^2}{a^2}+\frac{{(y^{\prime }-y_0)}^2}{b^2}+\frac{{(z-z_0)}^2}{c^2})}^{-1}\hfill \\ \hfill g_z& =-\frac{3\beta }{\gamma }\left(\frac{z-z_0}{c^2}\right){(1+\frac{{(x^{\prime }-x_0)}^2}{a^2}+\frac{{(y^{\prime }-y_0)}^2}{b^2}+\frac{{(z-z_0)}^2}{c^2})}^{-1}.\hfill \end{array}$$(5)

2.2. Configurations of the propagating jets

We have introduced the jet into the computational domain by constructing a cylindrical injection zone with a radius r_j of 1 kpc and a length of approximately 3 kpc. The jets are set to propagate along the negative x-axis in all simulations. To model the growth factors responsible for generating GRGs that extend beyond 1 Mpc in total length, and considering the extensive computational demands of each run, our simulations are designed to extend from +10 kpc to −710 kpc along the x-direction. Therefore, we modeled only the one-sided propagation of the jet. We inject it into the domain from (x,  y,  z)≡(0,  0,  0) in all cases (see, Fig. 1). This setup also required us to maintain the jet radius consistent with the resolution of the domain, resulting in 720 grid points in the x-direction (a resolution enhancement effect is discussed in Appendix A). Along the other two axes (y,  z), the domain was varied around the extent of [ − 180,  180] kpc to optimally capture the jet morphology (see, Table 1).

To discuss the thermodynamical properties of the jet, we incorporated an underdense, approximately pressure-matched, magnetic, relativistic jet into the domain. The jet is underdense with respect to the ambient medium, with ρ_j = 10⁻⁵ρ₀, and is slightly overpressured with the ambient environment (see Appendix B) at the launching zone (P_j ∼ 5.36 × 10⁻¹² dyn/cm²), following: Rossi et al. (2017), Musoke et al. (2020), Mukherjee et al. (2020), Giri et al. (2023). A toroidal magnetic field has been introduced along with the jet, configured as follows:

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

where (r,  θ) is the polar coordinate in the (y,  z) plane, perpendicular to the jet propagation axis, and so (B_y,  B_z) represent the components of the toroidal magnetic field in the (y,  z) plane. The term B_j is a constant and assigns the strength of the magnetic field, which is determined through the jet magnetization parameter σ estimated as the ratio of Poynting flux to the matter energy flux as

$$\begin{array}{c}\hfill \sigma =\frac{B_j^2}{{\mathrm{\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}$$(7)

where ρ_jh_j is the relativistic enthalpy, estimated for a Taub–Matthews equation of state (Taub 1948; Mignone et al. 2005), and Γ is the bulk Lorentz factor of the jet. We set the value of σ to 0.01, a value typically assumed at the early stages of jet propagation (Rossi et al. 2017; Mukherjee et al. 2020). To ensure that ∇ ⋅ B = 0 is always maintained in the computational domain, we employed the divergence cleaning method in our computations, as designed by Dedner et al. (2002). We conducted these simulations with second-order spatial accuracy, utilizing linear reconstruction and the HLLC Riemann solver (Mignone & Bodo 2006).

The jet power at the time of injection is controlled (in our runs) by the parameter Γ, the bulk Lorentz factor of the jet. For our simulations, Γ is set to values of 3 and 5 to represent two different scenarios. The jet power can then be calculated as follows (Matthews et al. 2019; Giri et al. 2023):

$$\begin{array}{c}\hfill Q_j=\pi r_j^2{\mathrm{v}}_j[\mathrm{\Gamma }(\mathrm{\Gamma }-1){\rho }_j{c}^2+\frac{\gamma }{\gamma -1}{\mathrm{\Gamma }}^2{P}_j+\frac{B_j^2}{4\pi }].\end{array}$$(8)

Here, v_j is the jet propagation speed during its injection into the domain (estimated from Γ values). The adiabatic index γ varies between 5/3 and 4/3 depending on the temperature of the medium, as prescribed by the Taub-Matthews equation of state.

The evaluated jet power is approximately 2.3 × 10⁴⁴ erg/s for Γ = 3 and 7.2 × 10⁴⁴ erg/s for Γ = 5. It is important to note that the radio power (at 1.4 GHz) of an evolved jet system is typically lower than the mechanical power of the jet by a factor of 3–300, as highlighted in studies by Bîrzan et al. (2004) and Rossi et al. (2017). Based on the models of Hardcastle (2018), the radio power at 150 MHz can also be estimated from the jet’s kinetic power. Relying on these estimates, the two different jet scenarios can then be enumerated as follows (see also, Mingo et al. 2019; Dabhade et al. 2020b):

1. Low jet power (“lp”):

   • Q_j erg/s; Expected radio power resembling powerful FR I sources

2. High jet power (“hp”):

   • Q_j erg/s; Expected radio power resembling standard FR II sources.

A combination of jet power and ambient medium configuration results in five distinct simulation settings, as highlighted in Table 1 and schematized in Fig. 1. Although GRGs are predominantly of FR II nature, the existence of ∼5 − 8% GRGs with FR I characteristics (Dabhade et al. 2020b; Simonte et al. 2024), necessitates verification of the applicability of hypothesized models in this regime of low-powered jet propagation, which is why the choice of “lp” above. These simulations therefore possess the potential to evaluate the effectiveness of growth factors, such as jet power or ambient medium hindrance, in forming the GRGs through comparative analysis. We note that due to the limited availability of tailed GRGs (Andernach et al. 2021), and the tendency of FR I sources to thrive in rich cluster or dynamically evolved group mediums (Zirbel 1997; Mingo et al. 2019), we have not explicitly considered simulating an FR I jet propagating through the outskirts of the environment. Additionally, the interstellar medium (ISM) of the host galaxy has not been specifically modeled here, as our study focuses on sources with jetted structures extending beyond one megaparsec, representing a late-stage evolution. The impact of the ISM on the early phases of such jetted structures will be addressed in a separate study.

All of our simulations have been performed in dimensionless units, with conversions to physical scales defined by the following units: length (L₀ = 100 kpc), velocity (c = 2.99 × 10¹⁰ cm/s), and density (ρ₀ = 0.001 amu/cc). Units for other variables can be derived from these three fundamental units. For instance, the unit of time can be calculated as the ratio of the unit length to the unit velocity, resulting in T₀ = 0.326 Myr. Finally, we note that two standard passive scalar tracers were injected into the domain: tr1 to identify jet material and tr2 to identify ambient material, with values ranging between 0 and 1.

3. Giant radio galaxy properties and observational relevance

This section provides a detailed overview of the findings of this work, addressing several key aspects: the morphological characteristics of the evolved structures and their relevance to observed sources, the implications for the ages of such giant jetted structures, and an examination of the (magnetohydro)dynamical properties to determine whether GRGs are similar to or distinct from their smaller counterparts.

3.1. Morphological distinction

3.1.1. Low-powered jet along the minor axis of the environment

The cocoon topology resulting from the propagation of a low-powered jet along the minor axis of the ambient medium at the late stage of its evolution is presented in Fig. 2 (“GRG_lp_min”). After an extended evolution period of approximately 166 Myr, a fat lobe has formed, measuring nearly 600 kpc in length and 360 kpc in width. The broad lobe is best recognized in the tracer plot (tracer values, tr1 ≥10⁻⁷; Mukherjee et al. 2020), shown in the top-right panel of Fig. 2. This panel also clearly illustrates the velocity distribution of cocoon material within the lobe via contours, which vary between 0.01c and 0.1c, indicating that the lobe constituents are evolving sub-relativistically. The jet ridge-line appears to be mildly relativistic even after propagating over such large distances. Since we have simulated only one-sided jet propagation, the total extent of the structure is expected to reach 1.2 Mpc, thus easily qualifying as a giant radio source.

Fig. 2. Simulation “GRG_lp_min” (low-powered jet propagating along the minor axis of the environment) showcasing the structure of the evolved GRG at the highlighted age. Top-left panel: x − y variation of (ρ/ρ0), accompanied by contours of velocity (0.1c for white, 0.01c for black). Top-right panel: Tracer distribution showing regions with jet material and their associated fraction. The highest color bar value indicates a grid cell completely filled with jet material (flipped for better visualization), and the contours represent the same velocities as earlier. Bottom-left panel: Three-dimensional representation of (ρ/ρ0) of the lobe, providing finer details with an extra dimension. Bottom-right panel: Contours of pressure values (Plobe/P0) indicating a minimal variation inside the cocoon (flipped for better visualization). The cocoon also appears to be overpressured compared to the ambient medium, as indicated by the ratio above the subplot. Here, ρ0 = 0.001 amu/cc, P0 = 1.5 × 10−6 dyn/cm2, and length is represented compared to L0 = 100 kpc. Details are given in Sect. 3.1.1.

The morphology of the cocoon is well captured in the 3D volume-rendered plot (bottom-left of Fig. 2), with density thresholds and color transparency applied as represented by the color bar. This diagram illustrating how the propagating jet transitions into a decollimated, broad flow (around 250 kpc), resulting in the formation of such an extended lobe. Similar transition zones in low-powered jets have also been recognized in: Monceau-Baroux et al. (2015), Nawaz et al. (2016), Giri et al. (2022a), showcasing numerically the process of generating plumes by a decollimated jet.

The pressure plot in the bottom-right of Fig. 2 shows contours of pressure values (P/P₀), indicating minimal variation of thermal pressure inside the cocoon. This is due to the lobes being overpressured relative to the ambient medium by a factor of 1.9, which effectively confines the lobe boundary against the ambient material (similar to, e.g., Matthews et al. 2019; Giri et al. 2023). This overpressured phase results in the generation of a bow shock, for example, at the forward-end of the lobe, which is clearly visible in the density slice plot, where the ambient material appears compressed (relevant to a low jet power case of Simionescu et al. 2009).

From Fig. 2, it is evident that a significant portion of the jet plasma accumulates in the frontal region of the lobe, likely observable in emission maps, as depicted by the 3D density map. The relevance of such structures in giant radio sources can be found in statistical studies such as: Subrahmanyan et al. (1996), Mack et al. (1997), Malarecki et al. (2013), Andernach et al. (2021), Oei et al. (2023a). Individual studies on giant lobes have also been conducted, with Oei et al. (2022) suggesting a lower resistance of the ambient medium to the jet, Cantwell et al. (2020) and Dabhade et al. (2022) indicating morphological transitions between jet and the fat-lobe, beside reporting a significantly higher time of evolution for such structures, Saripalli & Subrahmanyan (2009) and Hota et al. (2011) suggesting a preference for the GRGs to propagate along the minor axis of the host medium, and Oei et al. (2023b) pointing to the lowest lobe-to-ambient medium pressure ratio for such extended configurations, which are also the properties relevant to the present simulation (“GRG_lp_min”).

3.1.2. High-powered jet along the minor axis of the environment

In this scenario (“GRG_hp_min”), a high-powered jet propagating along the minor axis of the ambient medium forms a structure shown in Fig. 3 at nearly 68 Myr. The jet cocoon structure appears thin, with an axial ratio (length to width) of 5.5, similar to the GRG cases highlighted in Subrahmanyan et al. (1996). We emphasize the significantly shorter evolution time for the jet to reach 700 kpc compared to a low-powered jet along the same axis. With a jet power increased by a factor of 3, compared to the earlier case, the time required for the lobe to extend to 700 kpc is reduced by ∼2.5 times, indicating a substantial impact on the jet evolution time. This suggests that lobe age, in conjunction with radio power and morphology, can serve as a proxy for understanding GRG formation models.

Fig. 3. Morphological collage similar to Fig. 2 but for simulation “GRG_hp_min” (high-powered jet propagating along the minor axis of the environment) showcasing the structure of the evolved GRG at the highlighted age. We emphasize the much shorter evolution time for the jet to reach ∼700 kpc compared to the low-powered jet along the minor axis. Here, ρ0 = 0.001 amu/cc, P0 = 1.5 × 10−6 dyn/cm2, and length is represented compared to L0 = 100 kpc. Details are given in Sect. 3.1.2.

The frontal part of the lobe still moves rapidly, with speeds ≥0.1c, as seen from the contour plots in the top-rows of Fig. 3. This high speed causes the bow shock to remain attached to the jet head. Despite the high power of the jet, which would typically result in more back-flowing material from the jet head (Cielo et al. 2017), the rapid expansion of the jet leads to a thinner cocoon morphology. However, the combination of substantial back-flow and the lower axial ratio of the cocoon geometry results in higher cocoon pressure, increasing the lobe-to-ambient pressure difference by a factor of 4.0. This is illustrated in the pressure plot in the bottom-right of Fig. 3, which also shows the formation of back-flowing plasma from the jet ridgelines, likely around the recollimation shocks.

The bottom-left subplot of Fig. 3 showcases a 3D volume-rendered configuration, highlighting the formation of prominent recollimation shocks (see also the discussion in Appendix B). The cocoon develops Kelvin-Helmholtz instabilities at the jet ambient medium interface, a typical characteristic of kinetically dominated jets (Acharya et al. 2021; Giri et al. 2022a). However, these instabilities only marginally hinder the growth of the structure, and the jet head yet remains well-collimated, enabling rapid propagation through the denser environment. It is noteworthy that this jet flows from the center of the ambient medium, suggesting that GRG formation from Brightest Cluster Galaxies is theoretically plausible. This corroborates the findings of a recent study by Oei et al. (2024b), which suggests that luminous GRGs are more likely to be found in denser environments. However, the lower occurrence of such cases (BCG-GRGs) may be due to the general tendency of FR II jets to avoid residing in evolved clusters or rich groups (Gendre et al. 2013; Joshi et al. 2019). With the one-sided lobe length of the simulated structure exceeding 700 kpc, the total extent is expected to reach 1.4 Mpc, thereby comfortably meeting the definition of a GRG in the literature.

The availability of GRG morphologies such as this is commonly observed in statistical studies. A powerful jet with thin cocoon geometry and hotspots at the lobe ends can be observed in studies such as: Stuardi et al. (2020), Dabhade et al. (2020b), Simonte et al. (2022), Oei et al. (2023a). There have been individual studies on similar morphological sources, such as Subrahmanyan et al. (2008), indicating FR II GRGs expecting to have well-bounded lobes if not a relic; Jamrozy et al. (2008), highlighting a lower spectral age for such sources; Delhaize et al. (2021), showing the recollimation shock knots, and Machalski et al. (2011), suggesting ballistic propagation of the jet with a younger age of evolution.

3.1.3. Low-powered jet along the major axis of the environment

When a low-power jet propagates along the major axis of the ambient environment (“GRG_lp_maj”), it encounters higher jet disruption compared to other scenarios in our simulations, resulting in a shorter jet length over a longer evolution time. The one-sided lobe length reaches approximately 500 kpc (see, Fig. 4), indicating a total extent of 1 Mpc, qualifying it again as a GRG. Propagating along a path of lowest pressure gradient (Fig. 1), the higher opposing force to the jet propagation leads to a complex morphology, known as X-shaped radio galaxies (e.g., Mahatma et al. 2023). The formation of such sources is primarily believed to result from the back-flow of plasma from the jet head, which gets diverted along the minor axis of the ambient medium, a path with the highest pressure gradient that aids the backflowing plasma to follow that path, forming the structure known as a wing (Capetti et al. 2002; Hodges-Kluck & Reynolds 2011; Rossi et al. 2017; Giri et al. 2022b). Therefore, the wing typically consists of the older jet plasma and lies at a large angle to the active lobe (for a review, see Giri et al. 2024). Giant X-shaped radio galaxies are recognized as GRG-XRGs, examples of which can be found in the recent past as well (Saripalli & Subrahmanyan 2009; Cotton et al. 2020), which have been found to support the classic back-flow scenario.

Fig. 4. Collage similar to Fig. 2 but for simulation “GRG_lp_maj” (low-powered jet propagating along the major axis of the environment) showcasing the structure of the evolved GRG at the highlighted age. The structure has evolved to only 500 kpc on one side (∼1 Mpc in total extent) over a longer period of evolution, indicating higher jet frustration imposed by the environment while traveling along the major axis. The emerged structure is known as an X-shaped radio galaxy, and given its giant size, it is typically recognized as a GRG-XRG. Here, ρ0 = 0.001 amu/cc, P0 = 1.5 × 10−6 dyn/cm2, and length is represented compared to L0 = 100 kpc. Deatils are given in Sect. 3.1.3.

The active lobe evolves mostly with sub-relativistic speed, as expected (Fig. 4). As the jet grows in size, it appears to become disconnected from the back-flowing material supply to the wing, which is evident in the 3D volume rendered plot. This phenomenon may indicate the observation highlighted by Saripalli & Subrahmanyan (2009) that the larger the extent of the sources along the active lobe, the smaller the wing to lobe length ratio becomes. The lobe edges are susceptible to Kelvin-Helmholtz instabilities, leading to the mixing of ambient material with cocoon material. The jet flows along the major axis of the ambient medium with similar internal energy input as the case “GRG_lp_min”, which is now stored within a smaller active lobe volume, causing higher internal pressure. Consequently, the lobe becomes overpressured relative to the ambient environment by a factor of 3.3. It is interesting to mention that the wing structure has marginally lesser internal pressure than the active lobe, but still remains overpressured compared to the ambient medium, promising its slower but further growth (see the back-flowing case of Giri et al. 2023). Given the substantial backflow reaching the center, the hot plasma is anticipated to thermally influence both the interstellar medium and the nuclear region (e.g., Landt et al. 2010).

It is worth noting that the time required for the jet to form a lobe extending 500 kpc on one side is 196 Myr for the highlighted simulation setting. This represents a considerably longer dynamical evolution time, several times higher than the spectral ages typically associated with GRGs (see, for a review, Dabhade et al. 2023). While the dynamical age can vary by a factor of 1 to 5 compared to the spectral age for GRGs (Machalski et al. 2009), the obtained dynamical age of the formed structure here is substantially long. Therefore, as previously mentioned, combining the age of evolution of the GRG with its radio power and morphology may allow for the disentanglement of formation models.

The occurrence of X-shaped radio galaxies themselves comprises only 5 − 10% of known radio galaxies (Leahy & Parma 1992; Yang et al. 2019), suggesting that GRG-XRGs may be even less common. However, this fraction may be influenced by the sensitivity limits of radio telescopes, as demonstrated by Yang et al. (2019). With significant improvements in current and upcoming radio telescopes, the discovery rate of such systems is expected to increase dramatically (see, e.g., Mahatma 2023). At present, examples of GRG-XRGs can be found in statistical studies conducted by: Saripalli & Subrahmanyan (2009), Dabhade et al. (2020a), Bruni et al. (2021). While the later studies focuses on modeling the morphological appearances of these systems, the former showcases how the propagation axis of the jet in GRGs with wings tends to align along the major axis of the host galaxy. Given that the large-scale environment of the majority of XRGs follows the geometry and configuration of the host galaxy (Hodges-Kluck et al. 2010), it is expected that GRGs with wings would also propagate along the major axis of their large-scale environment, similar to what our simulation result indicates. In general, off-axis lobe distortions forming wings in GRGs are not rare, as observed in studies by Subrahmanyan et al. (1996) and Malarecki et al. (2013).

It is important to note that the formation of X-shaped radio galaxies may not solely result from classical back-flow scenarios. Other mechanisms, such as binary SMBH mergers, unstable mass accretion by the SMBH, or galaxy mergers causing rotation in the ambient environment that shifts the jet propagation direction, have also been reported to produce winged galaxies (see, review by: Gopal-Krishna et al. 2012, Giri et al. 2024). This suggests that, while not all (GRG-)XRGs may be attributed to back-flow processes, this mechanism could nontheless explain the majority of cases due to its naturally occurring nature (similar to the conclusions of, e.g., Saripalli & Subrahmanyan 2009, Joshi et al. 2019, Bruni et al. 2021, Dabhade & Krishna 2022).

3.1.4. High-powered jet along the major axis of the environment

A high-powered jet can overcome the increased jet frustration imposed by the environment when propagating along its major axis, as seen in the “GRG_hp_maj” case, resulting in the evolved structure shown in Fig. 5 at 137 Myr. This structure resembles the earlier scenario of “GRG_lp_maj” (Fig. 4), but with a shorter secondary lobe perpendicular to the active lobe, identified as mini-wings (primary to secondary lobe length ratio < 0.8; Cheung 2007). The reason for the shorter wing is that, with the increase in jet speed, the growth rate of the cocoon length also accelerates. Consequently, a substantial amount of the back-flow material struggles to reach the center and the wing, despite the faster production of back-flowing plasma for a high-speed jet. This phenomenon is numerically demonstrated by Rossi et al. (2017) and supported by observations (e.g., Gillone et al. 2016), indicating that only low-powered jets (∼ FR I/II boundary) propagating along the lowest pressure gradient path of the ambient medium (around the major axis) can generate wings from the overpressured back-flow model. The observational existence of such mini-winged GRGs has also been reported by Saripalli & Subrahmanyan (2009), who inferred that the more extended the source, the shorter the wings it possesses. A handful of similar sources have been identified in statistical studies by Malarecki et al. (2013) and Stuardi et al. (2020). Considering that the one-sided active lobe length reaches 700 kpc, the total extent of this source is expected to be 1.4 Mpc, thereby qualifying as a prominent GRG with lobe-to-ambient medium thermal pressure ratio of 2.9 (bottom-right of Fig. 5).

Fig. 5. Collage similar to Fig. 2 but for simulation “GRG_hp_maj” (high-powered jet propagating along the major axis of the environment) showcasing the structure of the evolved GRG at the highlighted age. This morphology exhibits the formation of a mini-wing perpendicular to the active lobe. Additionally, it demonstrates the phenomenon of rapid jet propagation (top-hat jet front), leaving the cocoon behind at the given time. Such topology are not uncommon for GRGs and have been observed in many instances. Here, ρ0 = 0.001 amu/cc, P0 = 1.5 × 10−6 dyn/cm2, and length is represented compared to L0 = 100 kpc. Details are given in Sect. 3.1.4.

The morphology also reveals a top-hat-like feature at the lobe front, best captured in the volume rendered image, indicating that the jet has broken out of the cocoon due to its propagation speed (Fig. 5). While this prominent feature is not consistently present throughout the course of its evolution, the existence of an arrowhead-like structure in the lobe front is evident, representing the rapid propagation of this high-powered jet, despite it being a major-axis source.

Such extended jet head features have been observed in GRGs as well, exemplified by: Saripalli et al. (1994), Safouris et al. (2009), Bruni et al. (2021), who also have documented the presence of multiple recollimation shocks, similar to ours (see, volume-rendered plot; bottom-left of Fig. 5; and also Appendix B). Chen et al. (2018) have also shown the evidence of such jet head extension along with the formation of fat lobe in the frontal part of the cocoon, due to the encounter of higher jet obstruction along its propagation direction as inferred from the galaxy distribution (similar to what we report). Another example could be found in the work of Oei et al. (2022), which illustrates a fat frontal cocoon attached to the top-hat jet structure. However, the estimated jet power derived from radio observations appears to be on the lower side, possibly due to non-negligible radiative cooling given the source’s size and redshift. These subtle features, combined with the lobe configuration, are challenging to identify in statistical studies unless the GRG morphology is highly resolved; however such structures are not rare (Simonte et al. 2024; Dabhade et al. 2020b).

3.1.5. High-powered jet along the edges of the environment

In the “GRG_hp_edge” scenario (Fig. 6), the jet propagates in a low-density environment along the edges of the tri-axial ambient medium. This results in the jet encountering less resistance compared to the other models considered here. Combined with the high power of the jet, the encountered drag against jet travel is negligible, allowing for rapid jet growth. Consequently, the jet reaches a distance of approximately 650 kpc (with a total expected extent of 1.3 Mpc) in just 49 Myr, a dynamical age comparable to that of SRGs (Harwood et al. 2017).

Fig. 6. Collage similar to Fig. 2 but for simulation “GRG_hp_edge” (high-powered jet propagating along the edges of the environment) showcasing the structure of the evolved GRG at the highlighted age. The dynamical age of 49 Myr is considerably low for a GRG and comparable to SRGs, hinting at the potential formation of even larger-scale GRGs within a satisfactory time frame compared to radiative ages. Here, ρ0 = 0.001 amu/cc, P0 = 1.5 × 10−6 dyn/cm2, and length is represented compared to L0 = 100 kpc. Details are given in Sect. 3.1.5.

Considering the arrowhead-like structure at the jet head, it can be inferred that the jet is propagating at a faster rate even after covering such distances. If the jet propagation continues unhindered for a considerably longer period, it could cover an enormous length-scale, explaining the rarer giant sources with lengths extending beyond 3 Mpc (Andernach et al. 2021). As mentioned earlier for the case “GRG_hp_min”, the combination of a high-power jet producing substantial back-flowing material, along with a thinner cocoon due to faster jet propagation, generates significantly higher lobe pressure. Compared to the pressure of the ambient medium at the edges, the lobe-to-ambient environment pressure difference becomes a factor of 7.9. Being high-powered, the recollimation shocks are quite evident in the pressure plot (bottom-right of Fig. 6; see also the discussion in Appendix B).

A detailed examination of the lobe morphology reveals interesting features such as the formation of a broad cocoon near the jet ejection zone, with the cocoon extending more toward the positive y-axis, that is, toward the outskirts of the ambient medium. This asymmetry is caused by two factors. First, the high-speed back-flowing plasma (as shown in the velocity contours in the top panel of Fig. 6) allows the jet material to reach the center rapidly. Then, as the material starts accumulating, the buoyancy force (acting outward from the ambient medium, in the direction of the pressure gradient) generates this asymmetric extension. Such buoyancy effects have been observed in smaller radio sources (e.g., Leahy & Williams 1984). However, for GRGs, we suspect observing such cocoon structures may be difficult as the extended zone consist of older, cooled jet plasma (see, PKS 2356-61 for a contrary; Ursini et al. 2018). Nonetheless, morphological structures similar to the frontal part of the lobe have been recognized in many statistical studies (Brüggen et al. 2021; Dabhade et al. 2020b; Stuardi et al. 2020; Simonte et al. 2022; Oei et al. 2023a). Individual studies, such as Machalski et al. (2008) referred the existence of an exceptionally low ambient medium density, Delhaize et al. (2021) identifying a similarly shaped but bent source likely evolving in the outskirts of a poor group, and Konar et al. (2004), Sebastian et al. (2018) suggesting a considerably low age and similar jetted morphology, implying our simulated scenario at play.

3.2. Temporal signatures

3.2.1. Evolution of axial ratios (length to width)

The ratio of the length to the width of the cocoon has been used for GRGs to investigate whether this parameter (axial ratio) can provide insights into distinguishing GRGs from their smaller counterparts. For example, studies by Subrahmanyan et al. (1996), Saripalli et al. (1996) have shown that there is no significant difference in axial ratios between GRGs and SRGs. Furthermore, a detailed understanding of the axial ratio geometry of the radio lobes aids in accurately estimating dynamical quantities (Hardcastle 2018), particularly concerning GRGs (Subrahmanyan et al. 2008; Machalski et al. 2008, 2009). Therefore, we have plotted the axial ratio distribution along the length of the cocoon for evolved systems as tracked by the tracer plots in Fig. 2–6. Additionally, we have plotted the temporal evolution of this parameter, providing better insights into the correlation with their earlier phases.

This is depicted in Fig. 7, where the lobe of the GRG is delineated by tracer values, tr1 ≥10⁻⁷ (Mukherjee et al. 2020). For cases where the jet propagates rapidly (“GRG_hp_min” and “GRG_hp_edge”), the cocoon exhibits similar evolutionary patterns in both their giant and smaller phases, indicating a self-similar expansion (bottom two panels of Fig. 7). As mentioned earlier, morphological similarities with GRG structures akin to these two cases have been observed, and a few studies have also reported axial ratio values for such cases that resemble our findings (e.g., Machalski et al. 2008). In this context, the study by Andreasyan & Sol (1999) emphasized that elongated radio sources with FR II-type jet power tend to preferentially orient themselves along the minor axis of the host galaxy, which aligns with the jet ambient medium configuration applied for these two (low-jet frustration) cases.

Fig. 7. Demonstration of temporal evolution of the structural topology originated from five simulations conducted in this work. Each row presents a simulation case, with labels attached to the rightmost part of that row. The evolution of the sources starts from the right column and progresses leftward, with the corresponding ages mentioned in the top corner of each box. The red-lines represent lateral extent of the radio lobe, that is, the width of the lobe, and the numbers attached to the lines denote the axial ratios (ratio of total length to width at that time) spread across the cocoon.

Considering the cases “GRG_lp_maj” and “GRG_hp_maj”, both produce similar morphology with wings. The geometry of the active lobe, inferred from the axial ratios, resembles a typical cylinder with minimal variation along its length. Due to higher jet frustration in the low-powered jet, the axial ratio values are slightly smaller, indicating a fatter lobe (with respect to length covered) compared to the high-powered jet. The jet with low power also produces relatively longer wings than the high-powered case, as reflected by the axial ratio. Compared to their earlier phases with smaller lengths, a mild deviation from self-similar evolution is observed, where increase of width near the frontal part of the lobe has been observed over time. Such axial ratio values, however, are not uncommon in GRGs, as highlighted by Subrahmanyan et al. (1996).

The run “GRG_lp_min” represents the case with the lowest jet power along the minor axis of the ambient medium, generating a fat lobe with corresponding axial ratio values as highlighted in Fig. 7. The morphology also deviating considerably from a cylinder of radius similar to the maximum width of the lobe. Furthermore, compared to its earlier phases, a deviation from self-similar evolution has been observed. Considering their axial ratio values, Machalski et al. (2009) noted resembling numbers for several grown GRGs.

Thus, the parameter “axial ratio” serves as a supplementary tool in deciphering the formation processes of GRGs, complementing the conclusions acquired from age, morphology, and radio power of GRGs. A concise overview of this is well represented in Fig. 7.

Finally, to assess the deviation in cocoon volume compared to geometric shapes like a cylinder or sphere, we analyzed the cocoon volume fraction. Initially, we determined the volume occupied by the cocoon material for all simulations at the specified final times outlined for each case in Fig. 7. Then, we calculated the length of the cocoon and the maximum width of the (primary) lobe to compare the actual cocoon volume (V_c) relative to a cylinder (V_cyl) with the length as cocoon length and maximum width as its diameter, and to a sphere (V_sphr) with a diameter equivalent to the length of the cocoon. The estimations are reported in Table 2. We note that these values are estimated for the one-sided cocoon-structure simulated in our cases.

3.2.2. Evolution of lobe lengths (comparison to theory)

The continuous variation of the length growth of cocoon over the course of its evolution time has been tracked and presented in Fig. 8. As, we are tracking the cocoon geometry, we have used tracer values as, tr1 ≥10⁻⁷, similar to Mukherjee et al. (2020), and then find out the largest cocoon length along the −x direction. Figure 8 showcases a diverse set of behaviors consistent with the earlier discussion, with a comparison to theoretical estimates.

Fig. 8. Development of cocoon length over the course of evolution for the five simulation cases considered in this work. The shaded region represents the theoretically allowed zones in the length-age diagram for the extended parameter space considered across the simulations. A diverse set of behavior compared to theoretical predictions has been observed here, as discussed in Sect. 3.2.2. Notably, the blue dashed line represents the case where the jet grows at the outskirts of the ambient medium, thereby, propagating faster than the theoretical predictions.

To deduce the theoretical framework of cocoon length versus age, we have adopted the formulation of Kaiser & Alexander (1997), which has been utilized in deducing the ages of GRGs corresponding to their length by, for example, Machalski et al. (2009) and Machalski (2011). They used the geometry of the bow-shock formed from mild to high-powered jets while propagating in a β-profiled environment, similar to our adopted models. Their analytical calculations highlighted the relationship between the length of the jet, consequently cocoon (L_j), and the time of evolution (t) as

$$\begin{array}{c}\hfill L_j=c_1{a}_0(\frac{t}{\tau }{)}^{3/(5-\beta )};\tau =(\frac{a_0^5{\rho }_0}{Q_0}{)}^{1/3}.\end{array}$$(9)

Here, c₁ is a dimensionless constant, a₀ represents the core radius of the ambient profile that varies with the values of (a,  b,  c) in our cases, ρ₀ is the density of the ambient medium at the core, and Q₀ primarily indicates the kinetic power of the jet. The formulation of c₁ is complicated and is dependant on the jet opening angle ϑ (assumed to be 5°; Horton et al. 2020), adiabatic indexes for ambient medium, sub-relativistic cocoon and relativistic jet material as γ_x = 5/3, γ_c = 5/3 and γ_j = 4/3 respectively (following Taub–Matthews equation of state; Mignone et al. 2005), and two other constants c₂,  c₃ as

$$\begin{array}{cc}\hfill c_1& ={\left[\frac{c_2}{c_3{\vartheta }^2}\frac{({\gamma }_x+1)({\gamma }_c-1){(5-\beta )}^3}{18\{9[{\gamma }_c+({\gamma }_c-1)\frac{c_2}{4{\vartheta }^2}]-4-\beta \}}\right]}^{1/(5-\beta )}\hfill \end{array}$$(10)

$$\begin{array}{cc}\hfill c_2& ={[\frac{({\gamma }_c-1)({\gamma }_j-1)}{4{\gamma }_c}+1]}^{{\gamma }_c/({\gamma }_c-1)}\frac{{\gamma }_j+1}{{\gamma }_j-1}\hfill \end{array}$$(11)

$$\begin{array}{cc}\hfill c_3& \approx \frac{\pi {\vartheta }^2}{c_2}.\hfill \end{array}$$(12)

Given the diverse set of parameter choices in our simulation cases, we evaluated the theoretically allowed range of values and highlighted them with a shaded region in Fig. 8.

Figure 8 shows that jets with lower kinetic power lag significantly behind theoretical estimates throughout their evolution, taking substantially longer to reach comparably shorter distances. This discrepancy arises because such jets tend to produce fat lobes by loosing the jet’s collimation or by bending, whose lateral extent grows faster over time (consequently growth in jet hindrance), deviating from self-similar evolution (see Sect. 3.2.1). In contrast, the theoretical evaluation assumes a self-similar expansion of the cocoon structure, where the jet is assumed to remain collimated. The evolution of such cases has been previously recognized by Rossi et al. (2017) and Giri et al. (2022b).

The cases with high-powered jets are found to reach the theoretically allowed rate of expansion. As indicated earlier, the faster evolving jet of “GRG_hp_min”, which develops self-similar growth, mostly follows the theoretical gray zone. However, in the beginning phases of it, the cocoon’s geometric evolution in the dense core of the ambient medium is almost spherical rather than elongated causing the lag. While travelling in the outskirts, the jet’s ram pressure and the environment’s pressure gradient force accelerate the jet propagation. For “GRG_hp_maj”, being a major-axis source, the core’s influence remains effective for a longer time, leading to the extended lag in the early phase of its propagation. Once out of the core region, the continuous ram pressure of the jet increases its speed. However, at later stages, the lobe expansion is starting to decreases again, as the cocoon is found to deviate from a self-similar expansion. Such evolutionary behavior of high-powered jets has been recognized previously by Hardcastle (2018).

For the case of “GRG_hp_edge”, the evolving length of the jet or cocoon exceeds the theoretical prediction, as expected. This occurs as the jet propagates at the outskirts of the ambient medium with minimal jet obstruction compared to the theoretical prediction, which assumes the jet always propagates from the center of the ambient medium. We again note the significantly shorter evolution time required for this jet to cover giant scales, indicating why GRGs prefer to avoid dense environments and why large-scaled GRGs (≥3 Mpc) evolve in underdense environments (Sankhyayan & Dabhade 2024).

3.3. (Magnetohydro)dynamical characteristics

3.3.1. Lobe pressure and expansion speed

Since GRGs also experience smaller phases during their early evolutionary stages, it is worthwhile to explore the changes in thermodynamical quantities throughout their evolution. In our simulations, the different GRG structures evolve over varying timescales. To ensure a consistent comparison across all cases, we plotted the lobe expansion speed and the fractional change in lobe pressure against the length of the cocoon, which serves also as a proxy for temporal evolution (see, Fig. 9).

Fig. 9. Diagrams showing variations of lobe expansion speed (left) and the fractional change in lobe pressure (right) versus the one-sided lobe lengths as they evolve over time. Since the change in lobe pressure and lobe expansion speed are intricately connected variables, we plotted them alongside each other to investigate whether the giant and smaller phases of the jet exhibit any distinct or dissimilar characteristics. The shaded region extends from 350 kpc onward. For more details, refer to Sect. 3.3.1.

The lobe speed was plotted from 10 Myr onward, and since the jets cover different distances by then, the initial starting points for simulation cases changes. As the jet cocoon structures evolve further, in the case of lobe expansion speed (left panel of Fig. 9), it is found to vary marginally. However, after the one-sided jet covers a distance of about 350 kpc, the expansion speed of the lobe decreases rapidly (for four cases) or remains nearly constant (for one case). The smaller-scale fluctuations, primarily appearing at the late stages, result from a combination of the jet’s ram pressure fluctuations due to its decollimation or bending and propagation along the higher pressure gradient in the outskirts of the environment.

The fractional change in pressure has been plotted starting from 20 Myr, since it represents a differential quantity (right panel of Fig. 9), defined as follows:

$$\begin{array}{c}\hfill \frac{\mathrm{\nabla }P}{P}=\frac{P_i-P_{i-1}}{P_i},\end{array}$$(13)

where i is the data file number as the simulation output, saved every 9.8 Myr. As with the previous analysis, the initial length values differ across the simulation cases. The lobe pressure of all simulated structures decreases with increasing length or time, as also indicated by the fractional pressure change plot. A detailed examination of which reveals that the values of ∇P/P exhibit either marginal changes (in one case) or more rapid decreases (in four cases) after the jet reaches a length of approximately 350 kpc, compared to the earlier phases of their evolution.

These insights suggest a possible phase transition for GRGs around a particular length scale, where they transition from their smaller evolutionary stages into their giant class. This could help determine whether GRGs differ in their lobe-evolution properties compared to SRGs. A similar phase transition behavior was observed in jet simulations by Rossi et al. (2024) while analyzing the evolution of maximum pressure in the jet environment system. These scale-free simulations further support our findings. However, to draw definitive conclusions for GRGs, the result requires more statistical input. Future studies should explore various other parameter spaces, including the evolution of GRGs in voids (Oei et al. 2022), in super clusters (Sankhyayan & Dabhade 2024), and the formation of GRGs in high-powered FR IIs and low-powered FR Is to investigate GRGs in the non-preferable zones of the 𝒫 − 𝒟 diagram (Delhaize et al. 2021; Simonte et al. 2024). These derived findings, however, could guide observational studies in statistically exploring the differences between SRGs and GRGs.

3.3.2. Comparison of lobe pressure to environmental pressure

To depict the pressure variation across the cocoon-ambient medium configuration more accurately, we have plotted a sliced pressure variation along the y-axis, perpendicular to the jet propagation, at the position of maximum cocoon width. This choice is made because the pressure variation across this direction reasonably represents the formed lobe, similar to findings by Hardcastle & Krause (2013) and Malarecki et al. (2013).

The result is showcased in Fig. 10, where the subplots present the same situation. However, the right panel represents cases where the maximum width of the cocoon fully occupies the y-axis, and tentative relaxed ambient pressure levels are drawn with disconnected lines in the bottom as obtained from the initial conditions when the jet has not entered the domain.

Fig. 10. Distribution of normalized pressure across the y-axis (lateral to jet propagation), indicating pressure variation within the lobe, immediate surroundings, and unperturbed ambient medium. The variation is shown at the position of maximum lobe width along the y-axis. The left and right panels represent the same situation for different runs. However, in the right panel, the computational domain is almost entirely occupied by the lobe material in the lateral direction. So, we use the ambient medium pressure values from the unperturbed initial condition for a tentative representation.

Both panels of Fig. 10 illustrate a general pattern of pressure variation. The pressure within the lobe remains fairly constant, with minor fluctuations indicating the turbulent nature of the cocoon. Large spikes in pressure within the lobes are typically associated with the jet and are centered around y = 0, except for the green dashed line representing “GRG_lp_maj”, where the active jet is significantly bent (see Fig. 4). Immediately beyond the lobe, significant increase in pressure are observed, corresponding to bow shocks that effectively confine the lobe boundaries. Several features in this zone are particularly noteworthy, such as jet bending, ambient matter entrainment at the lobe interface, and buoyancy effects for sources evolving at the medium’s edge, which make the bow shock strength asymmetric in both sides of the lobe. Further out, the unperturbed ambient medium exhibits lower pressure values. From the pressure jump values in Fig. 10 around the bow shock and the undisturbed ambient medium, one can estimate the Mach number (right panel is not recommended for this estimation, as the ambient values are tentative).

It is noteworthy that the pressure difference between the bow shock and the lobe is higher for high-powered jet cases, whereas for low-powered jets, it is not as prominent. This further suggests that as the power of the jet decreases or during phases of jet activity cessation, the lobe pressure is anticipated to dissipate and gradually align with the ambient pressure (see, e.g., Subrahmanyan et al. 2008; Cotton et al. 2020). Thus, the pressure disparity between the lobe and its immediate surroundings could serve as an indicator to discern a less active or relict phase of GRGs.

3.3.3. Evolution of magnetic field and total energy

Other parameters that require special attention include the evolution of the (dynamical) magnetic field and the total energy (comprising kinetic, thermal, and magnetic components).

The evolution of the magnetic field (left panel of Fig. 11) does not exhibit a clear pattern throughout its evolutionary track, showing nearly random phases of growth and decay. Despite these fluctuations, the values are not distinctly different across the simulations. Regardless of the varying jet ambient medium settings and evolutionary paths, all runs seem to converge toward a final value of around 0.15 μG. This suggests that, regardless of the different evolution history GRGs take, if the initial conditions for the magnetic field of the jet are similar, the cocoon is expected to have a comparable final average magnetic field strength. Determining whether this pattern is common to all GRGs requires further testing with a broader parameter space and longer simulations reaching greater extents than those we have currently simulated.

Fig. 11. The left panel illustrates the length (and consequently time) evolution of the dynamical magnetic field, while the right panel shows the evolution of total energy. Despite the seemingly random fluctuations in the magnetic field, the values tend to converge around a specific range, approximately 0.15 μG. In contrast, the total energy exhibits a relatively smooth evolution. The lines represent the five simulation cases considered in this work, all of which have a one-sided lobe length exceeding 500 kpc. Detail in Sect. 3.3.3.

We note that the estimated magnetic field here represents the dynamical magnetic field, which is expected to be lower than the equipartition field strength by a factor of ∼2 − 10, as suggested by observational and numerical studies with varying jet power, such as those by: Croston et al. (2005), Mahatma et al. (2020), Giri et al. (2022a,b). Such discrepancies in magnetic field values further contribute to the inconsistency in age estimation of extended radio sources between dynamical and spectral ages. The observation that these discrepancies are found in GRGs further supports the deviation in dynamical and spectral B-field strengths (Machalski et al. 2009; Machalski 2011). Based on these factors, the estimated dynamical B-field therefore appears relevant to GRGs (Schoenmakers et al. 2000; Machalski et al. 2008; Sebastian et al. 2018; Andernach et al. 2021).

Concerning the evolution of total energy (right panel of Fig. 11), the observed pattern is an expected trend relevant to extended-jetted sources (Giri et al. 2023). The total energy (E_T) has been evaluated as

$$\begin{array}{c}\hfill E_T=V_c[\frac{1}{2}{\rho }_c{v}_c^2+\frac{1}{{\gamma }_c-1}{P}_c+\frac{B_c^2}{8\pi }],\end{array}$$(14)

where the subscript c is attached to parameters to represent their values associated with the cocoon determined through the tracer (tr1) values. The 1st, 2nd, and 3rd parts in Eq. (14) refer to the kinetic, thermal, and magnetic components, respectively.

The total energy continuously increases as the jet injects energy into the lobe unhindered. The lobe is primarily dominated by kinetic energy, with the fraction of thermal energy being ≲0.001 (Oei et al. 2022), and the magnetic energy being almost negligible in comparison to the kinetic energy values (fraction ≲10⁻⁶). As a result, simulations with a higher jet power exhibit a higher total energy than those with a lower jet power. Interestingly, the low jet power cases show a slight decrease in the rate of energy increase, supporting the idea that the formation of the lobe, which deviates from self-similar expansion, allows the kinetic energy of the lobe material to relax slightly. This effect is particularly prominent in the case of “GRG_lp_min”. The estimated values align with what has typically been reported for GRGs, such as by Machalski et al. (2008) and Andernach et al. (2021).

4. Summary

Understanding the extent to which relativistic jets can propagate has been a topic of significant interest. This paper focuses on GRGs with megaparsec scale extents, which form a relatively small subset of radio-loud AGNs. Due to their rarity and the diverse range of jet environment configurations they are associated with, the origin of such extended structures has been a topic of considerable debate. This paper investigates the jet evolution phases of GRGs across a broad spectrum of ambient medium configurations and jet powers with the aim of numerically testing the plausible parameter space of the models and identifying their limitations. By adjusting the ambient environment to increase (major-axis sources) and decrease (minor-axis sources) jet frustration and varying the jet power between high and low values, we developed five distinct setups. These include a scenario where the jet propagates along the edges of the environment with minimal resistance to the jet flow. The findings of this study can be summarized as follows:

1. We observed the emergence of distinct cocoon morphologies for GRGs in varying jet environment settings, indicating the critical role played by these two physical aspects. We demonstrated that knowledge of morphology and jet power, when supplemented by the estimated age of the cocoon, seems to aid in the understanding of their likely formation history. Although precedents for such characteristics have been found in observational studies, they are still subject to verification with current and upcoming sensitive radio observations.

2. This study highlights the notably accelerated growth of the jet cocoon structure when FR II jets advance along the peripheries of an ambient medium and extend to megaparsec scales within a time frame comparable to SRGs on a scale of hundreds of kiloparsecs. This elucidates the prevalence of GRGs linked with FR II jets and their expansion in a low-density setting. This phenomenon likewise explains the development of giant radio structures on even larger scales within the standard jet propagation models.

3. Regardless of the models employed, the (active) GRG lobes consistently appear overpressured relative to the surrounding medium and are well confined by a shocked shell of this (ambient) medium. However, no consistent pattern of overpressure values has been observed; for example, there is no definitive correlation indicating that a high-power jet or a jet propagating along the path of maximum obstruction yields higher lobe pressures. Observing the lateral pressure variation in the lobe and its immediate surroundings, depending on whether the pressure decrease is sharp or gradual, may help distinguish whether a GRG is in a less active or relic phase.

4. Measuring the axial ratio (ratio of length to width) seems promising for understanding whether GRGs adhere to self-similar expansion patterns. While high-powered jets generally exhibit such expansion, low-powered jets diverge from this behavior primarily because of lobe formation, which significantly decelerates lobe expansion. However, this assessment may be influenced by the projection effect, a factor we aim to validate in a separate study.

5. While comparing GRGs with their smaller phases, a shift in the variation of lobe expansion speed and fractional change in lobe pressure has been observed around a length scale of ∼350 kpc (one-sided) for our formulated simulation settings, indicating the likely presence of a transition between SRGs and GRGs. However, this claim needs to be rigorously tested with even more extended parameter sets, supplemented by observational input (discussed below), in order to confirm its general appearance.

6. The (dynamical) magnetic field within the cocoon of GRGs exhibits no discernible pattern in its temporal variation. However, regardless of the models employed, consistent convergence to similar magnetic field strengths occurs when the initial magnetic field strength of the launched jet remains the same for all cases. The lobes of GRGs typically remain kinetically dominated.

7. The formation of off-axis lobes in GRGs resembling wings is a distinct feature that arises due to the influence of the tri-axial dynamics of the large-scale medium, and the lobes grow alongside the jet’s giant phases. Detection of these passively evolving lobes by present-day and forthcoming sensitive radio telescopes holds the potential to shed light on whether winged sources constitute a minor subset of radio-loud AGNs or if the prevalence of winged sources among radio galaxies is limited by the sensitivity thresholds of current radio telescopes.

This study focuses on the dynamical behavior of GRGs, emphasizing the morphology formed and variations in magnetohydrodynamical properties based on a set of jet environment configurations. While shedding light on GRG growth, an extension of this work could involve investigating the micro-physical processes occurring inside the jet cocoon, such as particle reacceleration or cooling effects (Vaidya et al. 2018; Mukherjee et al. 2021; Kundu et al. 2022), which are believed to have significant impacts on GRG appearances, particularly in high-redshifts (Hintzen et al. 1983; Schoenmakers et al. 2000; Konar et al. 2004). It is also crucial to extend numerical exploration of GRG properties to include examples such as GRGs found in voids (Oei et al. 2022) or in super clusters (Sankhyayan & Dabhade 2024) in order to assess the general applicability of the highlighted models. Moreover, examining the formation of GRGs in high-powered FR IIs and low-powered FR Is is essential to understanding their scarcity in less favored regions of the 𝒫 − 𝒟 diagram (Delhaize et al. 2021; Simonte et al. 2024). We also underscore the necessity of statistically modeling these giant jetted systems within expansive cosmological contexts (e.g., Davé et al. 2019; Nelson et al. 2019) to understand their feedback effects on the environment and their potential utility as probes for characterizing the warm-hot intergalactic medium.

Expanding upon this discussion, it is worth noting that a group of parameters in our simulated models, although carefully chosen to initially mimic standard jet environment systems, are held constant at single values. These include, for example, the density contrast between the jet and the ambient medium, the configuration and strength of the jet magnetic field, the eccentricity of the environment and parameters related to the King’s β-profile. Previous studies on various astrophysical jet systems have underscored their critical role in stabilizing and sustaining a jet’s evolutionary path. These studies have demonstrated how initial magnetic field strength and configuration can influence the development of Kelvin-Helmholtz or kink instabilities, ultimately affecting entrainment, jet expansion speed, and long-term evolution (Mukherjee et al. 2020; Rossi et al. 2024). Work by Rossi et al. (2017) has shown that different jet-to-ambient density contrast values impact jet flow speed, which in turn affects the linear growth of a jet system over time (see also its effect on FR-I type jets; Rossi et al. 2008). Focusing on the role of environmental eccentricity, Hodges-Kluck & Reynolds (2011) emphasized its impact on the development of prominent wings in X-shaped systems, where it facilitates the formation of notable XRGs (with a high wing-to-lobe length ratio) as eccentricity increases. A discussion on the influence of environmental parameters on the evolution of powerful jets can also be found in, for example, Hardcastle (2018).

It can thus be speculated that understanding the development of giant phases in straight jets is complex, especially from a theoretical perspective, given the combined influence of multiple underlying mechanisms. Consequently, more focused studies on GRGs are essential, and some aspects of these parameter studies will be addressed in our forthcoming works.

Acknowledgments

GG is a postdoctoral fellow under the sponsorship of the South African Radio Astronomy Observatory (SARAO). The financial assistance of the SARAO toward this research is hereby acknowledged (https://www.sarao.ac.za/). GG also acknowledges the pivotal discussions with William D. Cotton at the MeerKAT@5 conference, which sparked this study. JB acknowledges the support from Department of Physics and Electronics, Christ University, Bangalore. RPD acknowledges funding by the South African Research Chairs Initiative of the Department of Science and Innovation and National Research Foundation (Grant ID: 77948). 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/). 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 Inter-University 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). We extend our gratitude to Jeremy Smith (Inter-University Institute for Data Intensive Astronomy) and Anna Bosman (University of Pretoria) for their assistance in this regard. GG acknowledges financial support from the IDIA for a research visit to the University of Cape Town.

References

Appendix A: Effects of resolution enhancement

A.1. Setting up the resolution runs

In light of the lack of numerical studies available to validate models explaining the growth of GRGs, this work has been developed within a fully 3D relativistic magnetohydrodynamic framework in realistic jet environment settings, representative of typical GRG formation conditions. Furthermore, symmetry-breaking operations inherent to 3D dynamics have been incorporated, lending this study a distinctive methodological approach.

However, as with all simulation models, it is essential to assess the effect of resolution on the results, as achieving convergence in terms of resolution in 3D is often computationally challenging. For instance, there is general consensus that simulating jetted sources at high resolution requires approximately 8 grid zones within the jet radius to capture relevant details accurately (Stone & Norman 1994; Aloy et al. 1999). However, studies such as Mignone et al. (2010), Hodges-Kluck & Reynolds (2011), and Giri et al. (2022a) have demonstrated the use of even higher grid resolution within the jet radius (∼15 grid cells) to better capture transverse structures, such as shock formation along the jet beam and the influence of pressure gradients on evolving structures. A recent study by Dubey et al. (2023) employed an even greater number of grid cells within the jet radius (∼25), noting that while the overall external morphology may appear similar with lower resolution, key micro-physical processes within the jet cocoon, may vary significantly.

Achieving the grid resolution necessary to capture the majority of underlying physics is ideal but challenging, particularly for simulating jets at giant scales with the computing resources currently available. For instance, using 20 grid cells per jet radius would require a 3D grid of [14400 × 8000 × 8000] zones for a single run, which poses significant computational demands. Thereby, the highlighted earlier studies have either omitted several physical aspects to focus on a broader overview or have limited their simulations to the jet’s evolution at smaller scales. Driven by the goal of resolving both the jet and the large-scale lobe structures of GRGs with equal importance, and with our computational resources distributed across five models, we achieved only 2 grid cells per jet diameter. Consequently, a comparison with a slightly higher-resolution run has become essential.

We conducted our experiment on the case "GRG_hp_min," considering it a powerful source that is more susceptible to instabilities and capable of generating more shocks within the cocoon. The model presented in Table 1 for this case is designated as "Reference_Res." We then conducted two additional sets of simulations: one with 4 grid zones per jet diameter, resulting in a grid of [1440 × 600 × 600] to cover a domain of [720 × 300 × 300] kpc³, and the other performed in a zoom-in manner focusing on the early phases of jet evolution, with a grid of [1050 × 700 × 700] to cover a domain of [210 × 140 × 140] kpc³, raising ten grid cells per jet diameter. These simulations were designated as "High_Res (extended)’ and "High_Res (zoom-in)," respectively. A brief overview of the simulations are presented in Table A.1.

A.2. Insights from the resolution runs

From Fig. A.1, the immediate picture that arises is the lesser expansion rate of the jet cocoon structure in the higher-resolution case compared to the reference one. This result is not unexpected, given that higher-resolution runs can better resolve the effect of pressure gradient forces in the ambient medium (responsible for lateral expansion of the cocoon); additionally, the jet beam is now susceptible to instabilities, which affect the jet expansion speed. Regardless, it is clear that a high-powered jet propagating along the minor axis of the environment (i.e., in a direction of lower jet frustration) undergoes rapid lengthwise growth (∼300 kpc of one-sided jet length covered in 39 Myr of dynamical age). Consequently, as noted earlier, such a jet would not require much time to evolve into a standard GRG.

Fig. A.1. Comparison of x − y (z = 0) jet structures (based on tracer values) obtained from the "Reference_Res" and "High_Res (extended)" runs at ∼39 Myr (dynamical age). These two cases follow the longer-term evolution of the jet.

To illustrate the emergence of instabilities in the jet beam, which likely impact the jet expansion speed, we plotted a projected map of the highest 2000 Lorentz factor values for the "Reference_Res" case and 4000 values for the "High_Res (extended)" case (as it has twice the resolution). These points, expected to align with the primary jet beam, are projected onto a plane (mimics a typical 3D configuration) where Lorentz factor values are represented through a color bar, formulating a 2D histogram. This visualization is shown in Fig. A.2, where the two jet beams are slightly offset to fit within a single plot. The bending of the jet is significantly more evident in the higher-resolution case than in the reference case, where it is only mildly observed. This highlights that jet instabilities linked to powerful jets in GRGs are a natural occurrence that can impact the linear motion of the jets. Given the importance of these instabilities, a dedicated study is necessary to understand the growth history of these micro-features, as their evolution is also influenced by multiple factors, including the configuration and strength of the magnetic field, the density contrast between the jet and its surrounding medium. Recent high-resolution observational studies have also documented instabilities of this nature, reinforcing the relevance of these findings in understanding jet behavior in GRGs (Dabhade et al. 2022).

Fig. A.2. Projected visualization of the 3D jet beam structures for the "Reference_Res" and "High_Res (extended)" cases, represented by Lorentz factor values (highest few thousand values) shown in the color bar. The jet positions are shifted vertically to accommodate both cases within a single plot. Notably, the bending of the jet beam in the higher-resolution case is more pronounced than in the reference case, suggesting the development of instabilities with the jet’s motion.

Moving forward, we aimed to investigate the impact of resolution on external properties (jet length and cocoon volume) and internal property (internal energy) across the three runs, including the "High_Res (zoom-in)" case. First, we examined the evolution of jet length with dynamical age for the three cases, as illustrated in Fig. A.3 (left). Notably, the initial development of the jet structure around the denser inner core is consistent across all three cases. However, as time progresses, instabilities in the jet beam emerge, resulting in differences in jets’ linear evolution with age. Despite these variations, the evolution ultimately aligns with the theoretically permitted zones, suggesting that the results obtained at the reference resolution do not introduce any significant deviations from the established conclusions. This is further supported by the observation that the evolution of cocoon volume for all three cases closely follows one another over time. The deviations are not markedly different; in fact, the values eventually align with the primary evolution pattern (Fig. A.3 (middle)). This indicates that the 3D morphological characteristics, in terms of topology, remain consistent among the cases. In the case of internal energy variation, differences in values are observed (the patterns of evolution remain consistent though; Fig. A.3 (right)). As anticipated, the highest resolution run, which focuses on small-scale jet evolution, yields elevated energy values due to its ability to resolve more shock structures generated by the jet flow, resulting in increased heating of the cocoon material. In contrast, the reference resolution case exhibits lower heating, as expected from this perspective. However, if we consider the "High_Res (extended)" case, it ultimately converges with the "Reference_Res" case over a longer evolution time frame. Thus, it may be premature to assert that increasing resolution will consistently lead to higher internal energy detection within the cocoon of simulated GRGs. As the system ages, shock strength likely diminishes, referring that the observed internal energy increase could be an initial feature that subsides over time. Additionally, the cocoon volume evolution shows similar converging behavior at later stages, indicating that any early energy enhancement may smooth out as the system progresses.

Fig. A.3. Evolution patterns of the internal and external behavior of the simulated jet cocoon system over time for the three resolution cases discussed in Table A.1. The parameters plotted here are the lengthwise cocoon growth (left), cocoon volume growth (middle), and internal energy growth (right) across the three systems simulated at varying numbers of grid cells inside the jet diameter.

Appendix B: Observability of recollimation shocks

Our main simulations (Table 1) began with only 2 grid cells across the jet diameter. However, due to the initial jet pressure being ∼1.4 to 3.1 times higher than that of the surrounding environment near the injection zone, the jet rapidly expands to occupy multiple grid cells as propagation initiates. This expansion is evident in Fig. A.2, where the Lorentz factor values at the launching site quickly spread across several grid cells. Additional hints of this expansion appears in the 3D volume-rendered images in Fig. 2 - 6. However, to investigate this matter in detail, we have plotted the pressure values within the jet structures for our main runs at specific stages of their development (ages as mentioned in Fig. 2 - 6). Each plot shows a 2D zoom-in slice of the developed jet structure in the x − y (z = 0) plane, as illustrated in the Fig. B.1. The figure clearly illustrates the expansion across all cases, likely diminishing the impact of resolution limitations as the expanded jet occupies multiple grid cells. A comprehensive analysis of these patterns is still pending, but the influence of magnetohydrodynamic parameters on this behavior is evident. Our planned parameter studies will explore these effects in greater depth in future work.

Fig. B.1. Two-dimensional slices of pressure values (indicated by the color bar) are plotted in the x − y (z = 0) plane, with a zoom-in of the jet beam region, showing the developed structures of the five primary simulations from Table 1 at their respective evolutionary stages, as seen in Fig. 2 - 6 (with each case labeled in the subplots). Pressure contours in black solid-lines have been overlaid to illustrate the complex pressure distribution within the jet beam, highlighting zones of compression and rarefaction. This approach helps visualize the formation of recollimation shocks, particularly prominent in high-powered jet cases. In lower-powered jets, these shocks are noticeable in early evolution length-scales, though compression zones continue to appear as the jets undergo bending. This analysis also reveals that the jet beam expands across multiple grid zones shortly after injection, effectively reducing the impact of resolution constraints.

We further examined small-scale pressure behavior within the jet beam by overlaying pressure contours on the pressure map. This approach highlights zones of compression and rarefaction within the jet beam structure for each case. As shown in Fig. B.1, high-powered jets display distinct zones of pressure compression, with oscillatory variations in the Lorentz factor (not shown here) suggesting the formation of recollimation shocks. Hints of this oscillatory behavior can also be observed in the 3D volume-rendered plots of density presented in Figs. 3, 5, and 6. In contrast, low-powered jets exhibit pressure compression close to the launching sites, which gradually diminishes over distance. However, as low-powered jets are more susceptible to instabilities under the same magnetohydrodynamic parameters compared to high-powered jets, bending occurs further along their length, leading to renewed compression zones at these bend locations (e.g., Condon et al. 2021; Weżgowiec et al. 2024). Testing the visibility of these compressed zones in sky-projected emission maps is an area we plan to explore in future studies.

All Tables

All Figures

Fig. 1. Concise overview of the simulations conducted in this work. Left column: Brief description of the jet ambient medium configurations used to explore a broad spectrum of jet environment parameters. Middle column: Ambient configuration in relation to the Cartesian axes, with a 3D box representing the schematic simulation domain. Right column:x − y slice of density (at z = 0) illustrating the initial variation of (ρ/ρ0) as a color map, with overplotted pressure (P/P0) contours at values of [4.5,  2.2,  1.2,  0.67]×10−7, outward. Jet injection zones are marked with black and white dots in the density slices. The values of ρ0 and P0 are 0.001 amu/cc and 1.5 × 10−6 dyn/cm2 respectively, and lengths are represented in units of 100 kpc. To break the symmetry with respect to the jet flow, the ambient medium has been slightly rotated.

In the text

Fig. 2. Simulation “GRG_lp_min” (low-powered jet propagating along the minor axis of the environment) showcasing the structure of the evolved GRG at the highlighted age. Top-left panel: x − y variation of (ρ/ρ0), accompanied by contours of velocity (0.1c for white, 0.01c for black). Top-right panel: Tracer distribution showing regions with jet material and their associated fraction. The highest color bar value indicates a grid cell completely filled with jet material (flipped for better visualization), and the contours represent the same velocities as earlier. Bottom-left panel: Three-dimensional representation of (ρ/ρ0) of the lobe, providing finer details with an extra dimension. Bottom-right panel: Contours of pressure values (Plobe/P0) indicating a minimal variation inside the cocoon (flipped for better visualization). The cocoon also appears to be overpressured compared to the ambient medium, as indicated by the ratio above the subplot. Here, ρ0 = 0.001 amu/cc, P0 = 1.5 × 10−6 dyn/cm2, and length is represented compared to L0 = 100 kpc. Details are given in Sect. 3.1.1.

In the text

Fig. 3. Morphological collage similar to Fig. 2 but for simulation “GRG_hp_min” (high-powered jet propagating along the minor axis of the environment) showcasing the structure of the evolved GRG at the highlighted age. We emphasize the much shorter evolution time for the jet to reach ∼700 kpc compared to the low-powered jet along the minor axis. Here, ρ0 = 0.001 amu/cc, P0 = 1.5 × 10−6 dyn/cm2, and length is represented compared to L0 = 100 kpc. Details are given in Sect. 3.1.2.

In the text

Fig. 4. Collage similar to Fig. 2 but for simulation “GRG_lp_maj” (low-powered jet propagating along the major axis of the environment) showcasing the structure of the evolved GRG at the highlighted age. The structure has evolved to only 500 kpc on one side (∼1 Mpc in total extent) over a longer period of evolution, indicating higher jet frustration imposed by the environment while traveling along the major axis. The emerged structure is known as an X-shaped radio galaxy, and given its giant size, it is typically recognized as a GRG-XRG. Here, ρ0 = 0.001 amu/cc, P0 = 1.5 × 10−6 dyn/cm2, and length is represented compared to L0 = 100 kpc. Deatils are given in Sect. 3.1.3.

In the text

Fig. 5. Collage similar to Fig. 2 but for simulation “GRG_hp_maj” (high-powered jet propagating along the major axis of the environment) showcasing the structure of the evolved GRG at the highlighted age. This morphology exhibits the formation of a mini-wing perpendicular to the active lobe. Additionally, it demonstrates the phenomenon of rapid jet propagation (top-hat jet front), leaving the cocoon behind at the given time. Such topology are not uncommon for GRGs and have been observed in many instances. Here, ρ0 = 0.001 amu/cc, P0 = 1.5 × 10−6 dyn/cm2, and length is represented compared to L0 = 100 kpc. Details are given in Sect. 3.1.4.

In the text

Fig. 6. Collage similar to Fig. 2 but for simulation “GRG_hp_edge” (high-powered jet propagating along the edges of the environment) showcasing the structure of the evolved GRG at the highlighted age. The dynamical age of 49 Myr is considerably low for a GRG and comparable to SRGs, hinting at the potential formation of even larger-scale GRGs within a satisfactory time frame compared to radiative ages. Here, ρ0 = 0.001 amu/cc, P0 = 1.5 × 10−6 dyn/cm2, and length is represented compared to L0 = 100 kpc. Details are given in Sect. 3.1.5.

In the text

Fig. 7. Demonstration of temporal evolution of the structural topology originated from five simulations conducted in this work. Each row presents a simulation case, with labels attached to the rightmost part of that row. The evolution of the sources starts from the right column and progresses leftward, with the corresponding ages mentioned in the top corner of each box. The red-lines represent lateral extent of the radio lobe, that is, the width of the lobe, and the numbers attached to the lines denote the axial ratios (ratio of total length to width at that time) spread across the cocoon.

In the text

Fig. 8. Development of cocoon length over the course of evolution for the five simulation cases considered in this work. The shaded region represents the theoretically allowed zones in the length-age diagram for the extended parameter space considered across the simulations. A diverse set of behavior compared to theoretical predictions has been observed here, as discussed in Sect. 3.2.2. Notably, the blue dashed line represents the case where the jet grows at the outskirts of the ambient medium, thereby, propagating faster than the theoretical predictions.

In the text

Fig. 9. Diagrams showing variations of lobe expansion speed (left) and the fractional change in lobe pressure (right) versus the one-sided lobe lengths as they evolve over time. Since the change in lobe pressure and lobe expansion speed are intricately connected variables, we plotted them alongside each other to investigate whether the giant and smaller phases of the jet exhibit any distinct or dissimilar characteristics. The shaded region extends from 350 kpc onward. For more details, refer to Sect. 3.3.1.

In the text

Fig. 10. Distribution of normalized pressure across the y-axis (lateral to jet propagation), indicating pressure variation within the lobe, immediate surroundings, and unperturbed ambient medium. The variation is shown at the position of maximum lobe width along the y-axis. The left and right panels represent the same situation for different runs. However, in the right panel, the computational domain is almost entirely occupied by the lobe material in the lateral direction. So, we use the ambient medium pressure values from the unperturbed initial condition for a tentative representation.

In the text

Fig. 11. The left panel illustrates the length (and consequently time) evolution of the dynamical magnetic field, while the right panel shows the evolution of total energy. Despite the seemingly random fluctuations in the magnetic field, the values tend to converge around a specific range, approximately 0.15 μG. In contrast, the total energy exhibits a relatively smooth evolution. The lines represent the five simulation cases considered in this work, all of which have a one-sided lobe length exceeding 500 kpc. Detail in Sect. 3.3.3.

In the text

	Fig. A.1. Comparison of x − y (z = 0) jet structures (based on tracer values) obtained from the "Reference_Res" and "High_Res (extended)" runs at ∼39 Myr (dynamical age). These two cases follow the longer-term evolution of the jet.
In the text

Fig. A.2. Projected visualization of the 3D jet beam structures for the "Reference_Res" and "High_Res (extended)" cases, represented by Lorentz factor values (highest few thousand values) shown in the color bar. The jet positions are shifted vertically to accommodate both cases within a single plot. Notably, the bending of the jet beam in the higher-resolution case is more pronounced than in the reference case, suggesting the development of instabilities with the jet’s motion.

In the text

	Fig. A.3. Evolution patterns of the internal and external behavior of the simulated jet cocoon system over time for the three resolution cases discussed in Table A.1. The parameters plotted here are the lengthwise cocoon growth (left), cocoon volume growth (middle), and internal energy growth (right) across the three systems simulated at varying numbers of grid cells inside the jet diameter.
In the text

Fig. B.1. Two-dimensional slices of pressure values (indicated by the color bar) are plotted in the x − y (z = 0) plane, with a zoom-in of the jet beam region, showing the developed structures of the five primary simulations from Table 1 at their respective evolutionary stages, as seen in Fig. 2 - 6 (with each case labeled in the subplots). Pressure contours in black solid-lines have been overlaid to illustrate the complex pressure distribution within the jet beam, highlighting zones of compression and rarefaction. This approach helps visualize the formation of recollimation shocks, particularly prominent in high-powered jet cases. In lower-powered jets, these shocks are noticeable in early evolution length-scales, though compression zones continue to appear as the jets undergo bending. This analysis also reveals that the jet beam expands across multiple grid zones shortly after injection, effectively reducing the impact of resolution constraints.

In the text