Dynamics of Powerful Radio Galaxies

Abstract

Analytical models describing the dynamics of lobed radio sources are essential for interpretation of the tens of millions of radio sources that will be observed by the Square Kilometre Array and pathfinder instruments. We propose that historical models can be grouped into two classes in which the forward expansion of the radio source is driven by either the jet momentum flux or lobe internal pressure. The most recent generation of analytical models combines these limiting cases for a more comprehensive description. We extend the mathematical formalism of historical models to describe source expansion in non-uniform environments, and directly compare different model classes with each other and with hydrodynamic numerical simulations. We quantify differences in predicted observable characteristics for lobed radio sources due to the different model assumptions for their dynamics. We have made our code for the historical models analysed in this review openly available to the community.

1. Introduction

The first extragalactic radio sources were identified over seven decades ago by Bolton in the late-1940s [1]. Shortly after, in 1953, the first resolved image was captured of Cygnus A, now known as the archetypal “classical double” [2]. This breakthrough was followed by the pioneering efforts of radio survey groups in Australia [3] and the United Kingdom [4], which conducted the first large-scale radio surveys (for a comprehensive review, refer to [5]). The first quasar, 3C273, was discovered a few years later in 1963 (see [6,7,8] for a historical review). These pivotal developments laid the foundation for observational studies of radio galaxies. Subsequently, building upon Lynden-Bell’s ([9], 1969) proposal that black holes are responsible for the extreme luminosities observed in quasars, the 1970s saw the development of the first models describing the dynamical evolution of radio galaxies.

While varying in specific details, the majority of models in the literature share a similar overarching framework. These models in general consider two initially conical jets composed of particles that have been accelerated to relativistic velocities. The interaction between the jets and the intracluster medium (ICM) surrounding their host galaxy determines the subsequent evolution. Jets which retain sufficient forward ram pressure during their initial propagation phase (typically on galaxy scales) will be collimated by pressure from the ambient medium, or more likely, a build-up of plasma shed by the jet in the early stages of lobe formation [10,11]. Regardless, each collimated jet leads to the formation of a Mach disk (observable as a hotspot) and, as overpressured jet material flows back towards the equatorial plane, the inflation of a plasma lobe observable through synchrotron radiation; such objects are generally classified to have a Fanaroff and Riley [12] Type-II (FR-II) lobe morphology.

On the other hand, jets which suffer substantial entrainment (e.g., from stellar winds, [13,14] or the interstellar medium, [15]) will slow down to transsonic speeds (with respect to the internal lobe sound speed, of order 0.1c) and be disrupted. In this scenario, the jet momentum thrust is not important to the evolution of the lobe, and the role of the jet is simply to supply energy to the synchrotron-emitting lobe. The subsequent expansion of the lobes is determined by solving a set of fluid conservation equations; typically, the lobes undergo an initial momentum-dominated supersonic phase, followed by an adiabatic-expansion-driven coasting phase, and ultimately rising buoyantly in the later phases of evolution.

A large number of analytical and numerical models describing the evolution of active galactic nucleus (AGN) jets and lobes have been published since the first models were introduced over five decades ago. In this review, we summarise the different classes of lobed radio galaxy models, and provide a common framework to facilitate comparison both between the model classes and to more detailed hydrodynamic simulations. The dynamics of jetted Fanaroff and Riley [12] Type-I (FR-I) sources are not considered in this work; we refer the interested reader to the classical work of Bicknell [15].

The review is structured as follows. Section 2 presents early models by Rees [16] and Scheuer [17] (Model A), in which the forward thrust of uncollimated jets is balanced by the ram pressure from the ambient medium. We extend the formalism of Scheuer [17] to non-uniform environments, enabling direct comparison with more sophisticated modern analytical models for the first time. While these early models capture the fundamental aspects of jet termination and lobe formation, they neglect jet collimation by sideways ram pressure from the lobe (or ambient medium). In Section 3, we describe the analytical models proposed by Falle [18] and Kaiser and Alexander [19], which link jet collimation to subsequent lobe expansion. The past decade has seen the advent of environment-sensitive radio galaxy models beyond the self-similar solutions of Kaiser and Alexander [19] and related models (e.g., [20,21]). These models capture the evolution of lobe morphology in realistic environments [22,23] as well as the transition between jet- and lobe-driven expansion [23,24]. We describe these models in Section 4. In Section 5, we compare the consistency of predicted radio source dynamics between the main model classes, benchmark their evolutionary tracks against hydrodynamic simulations, and discuss their ability to generate synthetic AGN populations for parameter inversions. We conclude and suggest improvements to implement in the next generation of analytical models in Section 6.

2. Early Jet–Lobe Models

Rees [16] proposed a model in which conical jets emanating from the central engine of active galactic nuclei (AGNs) are pressure-balanced by the ram pressure of the ambient medium. In this model, the jets are a beam of low-frequency electromagnetic waves (LFEMW), the quantum field equivalent of a pair-plasma. The radiation pressure of this beam upon absorption by the ambient medium is $p_{rad}=Q/\left(\mathsf{\Omega }{R}^{2}c\right)$ for jet kinetic power Q and beam cross-sectional area $\mathsf{\Omega }{R}^2$ at radius R from the active nucleus. The pressure is increased if the interaction between the beam and ambient medium results in pair production, leading to a reaction pressure up to double that of the radiation pressure; the exact factor depends on the angle of reflected particles. The pressure contribution from the jet is, therefore, expressed as:

$$p_{jet}=\frac{{\kappa }_{1}Q}{\mathsf{\Omega }{R}^{2}c},$$

(1)

where $1⩽{\kappa }_1<2$ is a dimensionless constant describing both the fraction of the beam power that interacts with the ambient medium as radiation or particles, and the angle of reflection of those particles.

The forward jet thrust is balanced by ram pressure from the ambient medium $p_{ram}=\rho v^2$, where $\rho$ is the gas density of the assumed constant density ambient medium and $v=dR/dt$ is the advance speed of the jet head (see Figure 1). Scheuer [17] evaluated the resulting first-order differential equation in R assuming a constant jet half-opening angle ${\theta }_j$ (and thus, solid angle $\mathsf{\Omega }$) in their Model A. The jet length is related to source age t, and jet and environment parameters as:

$$R\left(t\right)={\left(\frac{{\kappa }_{1}Q}{\mathsf{\Omega }\rho c}\right)}^{1/4}{\left(2t\right)}^{1/2}.$$

(2)

The above approach assumes a constant density environment. However, the ambient medium on scales exceeding several kiloparsecs—typical of extended radio sources—is well represented by a symmetric power-law density profile of the form $\rho =k{r}^{-\beta }$, where the density parameter $k\equiv {\rho }_0{r}_0^{\beta }$ is a constant (e.g., [18,22]). The ram pressure applied by the ambient medium onto the expanding jet consequently weakens with distance from the central nucleus. In this review, we extend the Scheuer [17] Model A for the more general case of a power-law density profile, yielding:

$$R\left(t\right)={\left(\frac{{\kappa }_{1}Q}{\mathsf{\Omega }kc}\right)}^{1/(4-\beta )}{\left(\frac{(4-\beta )t}{2}\right)}^{2/(4-\beta )},$$

(3)

which converges to Scheuer’s [17] original constant density form when $\beta =0$, noting that $\rho =k$ in this limiting case. We use this more complete version of the model in the remainder of this work.

The synchrotron-emitting lobes inflated by the jets are typically assumed to have ellipsoidal morphology, with the ratio of major (aligned with the jet) to minor axes defined by the axis ratio $A=R/R_{\perp }$; we note that this differs (by a factor of 2) from the axial ratio $R_T=A/2$ of Kaiser and Alexander [19]. Scheuer [17] derive the volume of the ellipsoidal radio lobe associated with their LFEMW jets by considering the work done in inflating the cavity. The total energy in the cavity, U, increases over the time interval $\delta t$ due to the input kinetic power Q as:

$$\delta U=Q\delta t-p\delta V,$$

(4)

where $\delta V$ is the differential increase in volume, and the lobe pressure is given as follows (see, e.g., [19], specifically their Equation (15)):

$$p=\frac{U({\mathsf{\Gamma }}_c-1)(q+1)}{V},$$

(5)

where ${\mathsf{\Gamma }}_c$ is the polytropic index (or adiabatic index for an adiabatic equation of state; EoS) of the lobe plasma, and $q\ll 1$ is the ratio of energy in the magnetic field to that in the particles. This equation assumes the energy density (and thus, pressure) is approximately uniform throughout the lobe, a reasonable assumption given the high (∼$0.1c$) sound speeds in the lobes.

Equation (4) is a first-order differential equation describing the evolution of the total energy of the cavity. In Appendix A.1, we solve this differential equation assuming the cavity volume expands with increasing jet length as $V\left(R\right)={\kappa }_2{R}^{\alpha }$; here $\alpha ,{\kappa }_2>0$ are constants. This yields an expression for the lobe pressure in terms of the source age:

$$p\left(t\right)=\frac{Q({\mathsf{\Gamma }}_c-1)(q+1)}{{\kappa }_2[\alpha ({\mathsf{\Gamma }}_c-1)(q+1)+(4-\beta )/2]}{\left(\frac{\mathsf{\Omega }kc}{{\kappa }_{1}Q}\right)}^{\alpha /(4-\beta )}{\left(\frac{(4-\beta )t}{2}\right)}^{(4-\beta -2\alpha )/(4-\beta )},$$

(6)

where the constants $\alpha$ and ${\kappa }_2$ are evaluated below by considering the lobe volume evolution.

A major limitation of the Scheuer [17] model concerns the sideways expansion of the lobe, which is assumed to occur at the same velocity at every point on the lobe surface at any given time t. This expansion rate is derived by equating the lobe pressure to the ram pressure presented by the ambient medium as the lobe widens, i.e., $\rho v_{\perp }^2=p\left(t\right)$, where the ambient gas density is reasonably approximated as $\rho \sim k{R}^{-\beta }$. This sidewards expansion can only commence at locations already reached by the jet material. The half-width of the lobe at some location r along the jet axis is, thus, given by:

$$R_{\perp }\left(r\right)={\int }_{t\left(r\right)}^{t\left(R\right)}{v}_{\perp }\left(t^{*}\right)d{t}^{*},$$

(7)

where $t\left(r\right)$ is the time when the jet head reached the location r along the jet axis and $t\left(R\right)$ is the source age when the jet has its present length R. This integral is evaluated in Appendix A.2.

The lobe volume at the source age $t\equiv t\left(R\right)$ is found by integrating over all locations r along the jet axis:

$$\begin{array}{cc}\hfill V\left(R\right)\equiv {\kappa }_2{R}^{\alpha }& \hfill ={\int }_0^R\pi R_{\perp }{\left(r\right)}^{2}dr\hfill \\ \hfill & \hfill \propto R^{(14-5\beta -2\alpha )/2},\hfill \end{array}$$

(8)

Dimensional analysis shows that the only possible solution is for $\alpha =(14-5\beta )/4$. The lobe expansion is, therefore, not self-similar (which would require $V\propto R^3$, and hence, $\alpha =3$) unless $\beta =0.4$, i.e., in a gently declining density profile representative of cluster cores. For a uniform medium as considered by Scheuer [17], the exponent converges to $\alpha =\frac{7}{2}$, while in a steep environment ($\beta =2$) representative of the outer regions of groups or clusters, the lobe volume scales linearly with jet length leading to a rapid increase in the lobe axis ratio.

The constant of proportionality, ${\kappa }_2$, is similarly found by comparing terms not involving R, yielding:

$${\kappa }_2=\frac{16{\pi }^{1/2}{\left(\mathsf{\Omega }c\right)}^{3/4}{k}^{1/4}}{{\left[(14-5\beta )(18-5\beta )\right]}^{1/2}{\kappa }_1^{3/4}{Q}^{1/4}}{\left[\frac{({\mathsf{\Gamma }}_c-1)(q+1)}{(14-5\beta )({\mathsf{\Gamma }}_c-1)(q+1)+2(4-\beta )}\right]}^{1/2},$$

(9)

which converges to the expression found by (Scheuer [17], specifically their Equation (10)) in the limit of a uniform ambient medium and assuming ${\mathsf{\Gamma }}_c=\frac{4}{3}$.

This simple model neglects the sidewards ram pressure of the ambient medium (or lobe at later times) acting on the jet, which will lead to reconfinement shocks and ultimately the collimation of the jet. Scheuer [17] proposed a second model (their Model B) in which the jet is smoothly compressed into a collimated beam by the ambient medium; however, this assumption leads to unphysically narrow jets and consequently significantly faster jet-head advance speeds, which scale with jet cross-section y as $R\left(t\right)\propto y^{-4/9}{t}^{7/9}$. Scheuer [25] subsequently proposed that, in some sources, the jet may precess on a timescale which is short compared to the evolutionary timescale of the lobe. The time-averaged momentum flux of the jet is effectively spread over a larger cross-sectional area (equivalent to a larger jet opening angle), resulting in a slower growth rate along the jet axis (e.g., [26]). We return to this point in Section 5, when we compare the predictions of different models.

3. Lobe Expansion Models

The self-similar expansion model for the growth of quasar winds by Dyson et al. [27] spurred a new generation of analytical models based on the adiabatic expansion of the lobe bubble along a power-law ambient gas density profile. In particular, Falle [18] related the geometry and internal pressure of the expanding lobe to the dynamics of the jet (Section 3.1), enabling the Dyson et al. [27] model to be modified to consider the evolution of radio sources (Section 3.2). The Falle [18] model forms the basis for several models in the literature, including those by Kaiser and Alexander [19], Blundell and Rawlings [20], and Manolakou and Kirk [21].

3.1. Jet Collimation

Falle [18] revisited the dynamical modelling of jet collimation in their work in 1999 [18], considering an initially conical jet that reflects a strong shock of the surrounding medium upon reaching lateral pressure equilibrium. This reconfinement shock bounces between each side of the jet cavity preventing any further decay of the lateral thrust against the surrounding medium; this leads to a constant width, collimated jet with repeated cross-shaped structures of enhanced pressure and synchrotron emissivity. Falle [18] assumed that lobe formation occurs prior to jet collimation, and thus, that it is the lobe pressure which opposes the sidewards component of the jet thrust, not the ambient medium. Alexander [10] showed that jets may in fact be collimated by the ambient medium prior to the formation of lobes if $\sqrt{{\mathsf{\Gamma }}_x}sin{\theta }_j{M}_x<1$, where ${\mathsf{\Gamma }}_x$ is the polytropic index of the ambient medium, ${\theta }_j$ is the half-opening angle of the conical jet, and $M_x$ is the Mach number of the jet-head advance with respect to the sound speed of the ambient medium. Jet-head advance speeds are initially relativistic (e.g., VLBI observations of [28]), and thus, the external Mach number of the jet, is expected to be of an order of several hundred; only for very small opening angles of <1 degree would the jet be expected to be collimated by the ambient medium rather than the lobe.

The location of the initial reconfinement shock, z, is found by applying the Rankine–Hugoniot jump conditions for a plane-parallel shock to the lateral component of the flow that travels along the jet edge. The pressure of the lobe plasma is related by these conditions to the lateral component of the jet thrust as:

$$p\left(t\right)=\frac{2}{{\mathsf{\Gamma }}_j+1}{\rho }_j\left(z_1\right)v_j^2{sin}^2{\theta }_j,$$

(10)

where the jet plasma has a bulk velocity $v_j$, polytropic index ${\mathsf{\Gamma }}_j$, and density ${\rho }_j\left(z_1\right)$ at the critical radius $z_1$ at which the jet begins to collimate (i.e., location where the ram pressure first matches the lobe pressure). The density of the jet plasma will remain constant after this point until it reaches the jet head. We can, thus, derive an expression for the pressure acting on the contact discontinuity between the jet head and surrounding shocked gas using the shock jump conditions. That is:

$$\begin{array}{cc}\hfill p_h\left(t\right)& \hfill =\frac{2}{{\mathsf{\Gamma }}_j+1}{\rho }_j\left(z_1\right)v_j^2\hfill \\ \hfill & \hfill =\frac{p\left(t\right)}{{sin}^2{\theta }_j},\hfill \end{array}$$

(11)

where the second equality is obtained upon substitution of Equation (10). This expression has the same form as found by Kaiser and Alexander [19] (specifically their Equation (36)), and yields comparable pressure ratios to the numerically informed value obtained by Komissarov and Falle [29].

The lobe and jet-head region are surrounded by a shell of swept-up ambient medium that has been overrun by the bow shock generated by the expanding jet. This shocked gas is in approximate pressure equilibrium with the proximate lobe/jet-head plasma, but has significantly higher density, and thus, lower temperature (e.g., [30,31]), i.e., $p_s\left(t\right)\sim p_h\left(t\right)$, where $p_s\left(t\right)$ is the pressure just inside the bow shock along the jet axis. Together with conservation equations, these relationships between the conditions in the lobe cavity, the bow shock, and the ambient medium are sufficient for describing the evolution of the expanding radio source.

3.2. Lobe Adiabatic Expansion

Falle [18] presented their model in terms of the volume and pressure of the lobe; however, the surrounding shocked gas shell also receives a non-negligible fraction of the input energy from the central nucleus. Therefore, we express their equations in a more complete form considering both the lobe and shocked gas shell consistent with later work by Hardcastle [23], Turner and Shabala [32], and Turner et al. [24]. These models assume that it is primarily the thermal pressure of the lobe plasma that drives the source expansion; this pressure is uniform throughout the lobes due to high internal sound speed. The first law of thermodynamics relates the jet kinetic energy input to the thermal pressure p in the lobe and shell, and shocked gas volume $V_s$ (see Equation (6) of [10]):

$$V_s\frac{dp}{dt}+{\mathsf{\Gamma }}_{c}p\frac{d{V}_s}{dt}=({\mathsf{\Gamma }}_c-1)Q,$$

(12)

where ${\mathsf{\Gamma }}_c$ is the adiabatic index of the shocked gas and lobe plasma ($\frac{5}{3}$ for a non-relativistic fluid), and Q is the power injected into the shocked shell by the jet. This equation can be rewritten in terms of the pressure at the interface of the shocked shell and ambient medium (along the jet axis) using Equation (11). That is:

$$V_s\frac{d{p}_s}{dt}+{\mathsf{\Gamma }}_c{p}_s\frac{d{V}_s}{dt}=\frac{({\mathsf{\Gamma }}_c-1)Q}{{sin}^2{\theta }_j}.$$

(13)

Falle [18] showed that the shocked shell expands in a self-similar manner, leading to a constant scaling between the volume and cube of the jet length. That is:

$$V_s\left(R_s\right)={\kappa }_3{R}_s^3.$$

(14)

where $R_s$ is the radius of the shocked gas shell along the jet axis and ${\kappa }_3$ is a constant of proportionality. Kaiser and Alexander [19], and some subsequent authors (e.g., [33]), modelled the lobe as cylindrical (see Figure 2), with the major axis of the shocked gas shell being a factor of $A_s=1/sin{\theta }_j$ longer than the minor axis. The volume of the shocked shell is then:

$$V_s\left(R_s\right)=\pi {sin}^2{\theta }_j{R}_s^3\left(\mathrm{cylinder}\right).$$

(15)

In this approach, expansion in the jet direction is driven by the ram pressure acting on the jet-head region, $p_h\left(t\right)$, whilst sideways expansion is driven by shocked shell thermal pressure, $p\left(t\right)$. The more realistic assumption of an ellipsoid shocked shell adds a factor of $\frac{2}{3}$ to the above equation; however, the sidewards ram pressure then becomes a strong function of distance along the jet axis. The resulting average lobe and shocked shell pressures must be calculated numerically (see Section 4.1) rather than using the simple relation derived in Section 3.1.

The pressure of the shocked shell along the jet axis, $p_s\left(t\right)\sim p_h\left(t\right)$, is related to the density of the ambient medium, ${\rho }_x=k{R}_s^{-\beta }$, using the Rankine–Hugoniot shock jump conditions:

$$p_s\left(t\right)=\frac{2}{{\mathsf{\Gamma }}_x+1}k{R}_s^{-\beta }{\left(\frac{d{R}_s}{dt}\right)}^2,$$

(16)

where ${\mathsf{\Gamma }}_x$ is the adiabatic index of the ambient medium surrounding the shocked gas shell.

Substituting Equations (15) and (16) for the shocked shell volume and pressure along the jet axis into Equation (13) yields a second-order non-linear differential equation for the shell radius:

$$2R_s^{3-\beta }\frac{d{R}_s}{dt}\frac{d^2{R}_s}{d{t}^2}+(3{\mathsf{\Gamma }}_c-\beta )R_s^{2-\beta }{\left(\frac{d{R}_s}{dt}\right)}^3=\frac{({\mathsf{\Gamma }}_c-1)({\mathsf{\Gamma }}_x+1)Q}{2k\pi {sin}^2{\theta }_j}.$$

(17)

This equation can be solved by trialling another power-law solution, with the exponent again constrained by dimensional analysis. This yields:

$$R_s\left(t\right)={\left[\frac{({\mathsf{\Gamma }}_c-1)({\mathsf{\Gamma }}_x+1){(5-\beta )}^{3}Q}{18(9{\mathsf{\Gamma }}_c-4-\beta )k\pi {sin}^2{\theta }_j}\right]}^{1/(5-\beta )}{t}^{3/(5-\beta )}.$$

(18)

Kaiser and Alexander [19] included an additional correction in the denominator of their equivalent expression due to the energy associated with the higher pressure jet-head region; the $(9{\mathsf{\Gamma }}_c-4-\beta )$ term in Equation (18) becomes $(9[{\mathsf{\Gamma }}_c+({\mathsf{\Gamma }}_c-1)/{sin}^2{\theta }_j]-4-\beta )$. We note that this correction assumes that the differential increase in volume of the lobe and hotspot are equal as the source grows; this assumption is not particularly realistic for an ellipsoidal lobe geometry.

4. Semi-Analytic Models

Improved computation has recently enabled a new generation of analytical models with added complexity. These models typically solve systems of differential equations which lack an analytic solution to describe the evolutionary history of the radio source. Below, we summarise the main developments, including atmospheres beyond power-law density profiles (Section 4.1), considering both the ram and thermal pressure contributions to the expansion along the jet axis (Section 4.2), and modelling the relativistic jet in a distinct expansion phase prior to the onset of lobe formation (Section 4.3).

4.1. RAiSE (Version 2015)

Turner and Shabala [22] developed a semi-analytic model, Radio AGN in Semi-analytic Environments (RAiSE), based on the theory of the Falle [18] class of models (Section 3). These authors extended existing analytic approaches by considering piece-wise solutions to the governing differential equations in two dimensions. The RAiSE model included three key improvements over the earlier models: (1) ambient medium consistent with X-ray observations of clusters and semi-analytic galaxy formation models, (2) angular dependence of expansion velocity across the ellipsoidal contact surface, and (3) modelling of the morphological transition from supersonic to subsonic lobe expansion by using complete differential equations rather than limiting cases, which yield analytic expressions.

The lobe and shocked shell in their model are constructed from an ensemble of small angular volume elements in assumed pressure equilibrium. Each element of fixed angular width $d\theta$ is assumed to receive a constant fraction of the jet power as the cavity expands. This assumption yields self-similar expansion at early times when the shocked shell is expanding in the strong-shock supersonic limit, as in the earlier models of Kaiser and Alexander [19]. The volume of each small angular element of the shocked shell, $[\theta -\delta \theta /2,\theta +\delta \theta /2)$, is given by:

$$\delta V_s\left(\theta \right)=\frac{2\pi R_s^3\left(\theta \right)}{3}sin\theta \delta \theta ,$$

(19)

where $\theta$ is the angle between some location on the surface of the shocked shell and the jet axis and $R_s\left(\theta \right)$ is the radius of the initially ellipsoidal shell at that location (see Figure 3). Importantly, the shocked gas shell does not expand self-similarly as the steepness of the ambient gas density profile encountered by the non-spherical shell will in general differ across its surface, leading to different growth rates; this prediction is consistent with the higher axis ratios observed in the largest radio sources Mullin et al. [34].

The initial radius of each volume element is related to that along the jet axis by a geometric factor as $R_s(\theta ,t\to 0)={\eta }_s\left(\theta \right)R_s(\theta =0,t\to 0)$, where $\theta =0$ is aligned along the jet axis and ${\eta }_s\left(\theta \right)$ is defined as:

$${\eta }_s\left(\theta \right)=\frac{1}{\sqrt{({sin}^2\theta /{sin}^2{\theta }_j)+{cos}^2\theta }},$$

(20)

where we have assumed the same relationship between the jet half-opening angle, ${\theta }_j$, and axis ratio of the shocked shell, $A_s$, as discussed in Section 3.

Following Turner and Shabala [22], and later work by Turner and Shabala [32] and Turner et al. [24], the adiabatic expansion of each angular volume element is related to the pressure imparted on that element at the surface, $p_s\left(\theta \right)$, its volume $\delta V_s\left(\theta \right)$, and the fraction of the input jet power associated with that element, $Q\delta \lambda \left(\theta \right)$. The function $\delta \lambda \left(\theta \right)$ is defined in Equation (20) of Turner et al. [24]. The first law of thermodynamics in Equation (13) gives:

$$\frac{d{p}_s\left(\theta \right)}{dt}\delta V_s\left(\theta \right)+{\mathsf{\Gamma }}_c{p}_s\left(\theta \right)\frac{d\left[\delta V_s\left(\theta \right)\right]}{dt}=({\mathsf{\Gamma }}_c-1)Q\delta \lambda \left(\theta \right).$$

(21)

Away from the contact surface, the pressure in the lobe is calculated as the spatial average of the surface pressures $p_s\left(\theta \right)$.

Turner and Shabala [22] used a similar expression for the pressure at the contact surface to Equation (16), but additionally considered: (1) the orientation of the expanding surface as it impacts the ambient medium and (2) terms describing evolution in the transonic and subsonic expansion regimes. That is, for expansion in the supersonic and transonic phases, we have:

$$p_s\left(\theta \right)=\frac{2}{{\mathsf{\Gamma }}_x+1}k{R}_s^{-\beta }\left(\theta \right){\left(\frac{{\zeta }_s\left(\theta \right)}{{\eta }_s\left(\theta \right)}\frac{d{R}_s\left(\theta \right)}{dt}\right)}^2-\frac{{\mathsf{\Gamma }}_x-1}{{\mathsf{\Gamma }}_x+1}\left(kl\right)R_s^{-\beta }\left(\theta \right),$$

(22)

where ${\zeta }_s\left(\theta \right)$ is a further geometric factor defined in Equation (13) of Turner et al. [24]. The radial temperature profile of the ambient medium is defined by Turner and Shabala [22] as $T=(\overline{m}/k_b)l{r}^{-\varsigma }$ for Boltzmann constant $k_B$, and average particle mass $\overline{m}\sim 0.6m_p$ (where $m_p$ is the proton mass). Below, we present their results assuming an isothermal medium for consistency with other authors, i.e., adopting $\varsigma =0$. The pressure in the subsonic regime is equal to the ambient pressure (see Equation (4) of [22]).

The second-order non-linear differential equation that results from substituting Equations (19) and (22) into Equation (21) cannot in general be solved to yield an analytic solution. Turner and Shabala [22], instead, rewrote the resulting equation as a system of two coupled first-order ordinary differential equations. These differential equations describe the velocity and acceleration at the contact surface of a given volume element $\delta V_s\left(\theta \right)$:

$$\begin{array}{cc}\hfill {\dot{R}}_s\left(\theta \right)& =v_s\hfill \\ \hfill {\dot{v}}_s\left(\theta \right)& =\frac{3({\mathsf{\Gamma }}_x+1)({\mathsf{\Gamma }}_c-1)Q{R}_s^{\beta -3}\delta \lambda }{8\pi v_s{({\zeta }_s/{\eta }_s)}^{2}ksin\theta \delta \theta }+\frac{(\beta -3{\mathsf{\Gamma }}_c)v_s^2}{2R_s}\hfill \\ & \hspace{1em}+\frac{({\mathsf{\Gamma }}_x-1)(3{\mathsf{\Gamma }}_c-\beta )l}{4R_s{({\zeta }_s/{\eta }_s)}^2},\hfill \end{array}$$

(23)

where $v_s$, $R_s$, $\delta \lambda$, ${\zeta }_s$, and ${\eta }_s$ are explicit functions of $\theta$, whilst the properties of the ambient medium, k, l and $\beta$, are implicit functions of $\theta$ as different sections of the contact surface reach a given distance into the spherically symmetric environment at different times. Turner and Shabala [22] used a standard fourth-order Runge–Kutta method to solve this system of equations, providing the analytic solution for the strong-shock limit as an initial condition.

4.2. Hardcastle Model

Hardcastle [23] improved on the lobe-dominated expansion model of Turner and Shabala [22] by explicitly considering the momentum flux of the jet plasma. This model considers expansion for two angles across the surface of the shocked shell, $\theta =0$ and $\frac{\pi }{2}$ (see Figure 4). The choice of these two angles is sufficient to model changes to the axis ratio of the lobe and shocked shell in lobed FR-IIs, with the larger number of angles considered by Turner and Shabala [22] only important in the transonic and subsonic expansion phases when the lobe deforms from its initial ellipsoidal shape.

The shocked shell pressures are derived from the ram pressure component along the jet axis (Equation (1)), and a component due to the internal energy of the relativistic lobe plasma acting along both axes (Equation (5), with $U=Qt$ and $q\ll 1$). The lobe and shocked gas are assumed to be in pressure equilibrium, as discussed in Section 3. The shocked shell pressures along the major and minor axes are then given by:

$$\begin{array}{cc}\hfill p_s(\theta =0)& \hfill =\frac{\epsilon Q{R}_s(\theta =0)}{2cV}+\frac{({\mathsf{\Gamma }}_c-1)\xi Qt}{V}\hfill \\ \hfill p_s(\theta =\frac{\pi }{2})& =\frac{({\mathsf{\Gamma }}_c-1)\xi Qt}{V}.\hfill \end{array}$$

(24)

where $V/R_s(\theta =0)$ is the cross-sectional area of the jet-head region for a cylindrical lobe of volume V, $\epsilon \sim 4$ here acts as a geometric correction factor reflecting more realistic lobe shapes, and $\xi \sim \frac{1}{2}$ is the fraction of the input jet kinetic power found in hydrodynamic simulations to be stored as internal energy of the relativistic lobe plasma; the remainder is stored as thermal and kinetic energy in the shocked gas shell.

Hardcastle [23] related the volume of the lobe, V, to that of the shocked shell, $V_s$ (including the interior lobe), by considering the ratio of total internal energies:

$$\frac{V}{V_s}=\frac{({\mathsf{\Gamma }}_c-1)\xi Qt}{[\xi {\mathsf{\Gamma }}_c+(1-\xi ){\mathsf{\Gamma }}_s-1]Qt+f(N,T,v_s)},$$

(25)

where $f(N,T,v_s)$ is a function describing the internal energy of the ambient medium swept up by the bow shock. This function depends on the total number of swept-up particles, N, their temperature, T, and their bulk velocity due to the expansion of the shocked shell, $v_s=d{R}_s/dt$ ([23], specifically their Equations (5) and (7)). For young sources, when the thermal energy of these particles is lower than the energy supplied by the jet to the forming shocked gas shell, the ratio of lobe to shocked shell volumes tends to $\xi$ for ${\mathsf{\Gamma }}_c={\mathsf{\Gamma }}_s=\frac{5}{3}$, or $\xi /(2-\xi )$ if the lobe is assumed to have a relativistic plasma (${\mathsf{\Gamma }}_c=\frac{4}{3}$). The volume of the shocked gas shell is of course also directly related to the lengths of the major and minor axes; for an ellipsoidal geometry this gives:

$$V_s=\frac{2\pi R_s(\theta =0)R_s^2(\theta =\frac{\pi }{2})}{3}.$$

(26)

We can, therefore, express the pressure along the major and minor axes of the shocked shell (Equation (24)) as a function of both axis lengths upon substitution of Equations (25) and (26).

Hardcastle [23] rewrote the Rankine–Hugoniot shock jump conditions (Equation (22)) in terms of the velocities along the major and minor axes of the lobe. However, unlike the previously discussed models, Hardcastle [23] expressed the jump conditions in terms of the sound speed, $c_s$, and adiabatic index, ${\mathsf{\Gamma }}_s$, of the shocked gas shell surrounding the lobe. The expansion rates along the major and minor axes are then given by:

$$\begin{array}{cc}\hfill \frac{d{R}_s(\theta =0)}{dt}& =c_s\sqrt{\frac{({\mathsf{\Gamma }}_s+1)[p_s(\theta =0)/p_x\left(R_s(\theta =0)\right)]-({\mathsf{\Gamma }}_s-1)}{2{\mathsf{\Gamma }}_s}}\hfill \\ \hfill \frac{d{R}_s(\theta =\frac{\pi }{2})}{dt}& =c_s\sqrt{\frac{({\mathsf{\Gamma }}_s+1)[p_s(\theta =\frac{\pi }{2})/p_x\left(R_s(\theta =\frac{\pi }{2})\right)]-({\mathsf{\Gamma }}_s-1)}{2{\mathsf{\Gamma }}_s}},\hfill \end{array}$$

(27)

where $p_x\left(r\right)$ is the spherically symmetric ambient gas pressure profile. This coupled system of non-linear ordinary differential equations can be solved using a standard fourth-order Runge–Kutta method, with initial conditions of $R_s(\theta =0)=c{t}_0$ and $R_s(\theta =\frac{\pi }{2})=c{t}_0$ for some small time $t_0$.

4.3. RAiSE (Version 2023)

The RAiSE [22] model discussed in Section 4.1 was first extended to make predictions for the spatial distribution of emission at radio [35] and X-ray wavelengths [32], and was most recently used to incorporate important changes to jet and lobe dynamics. The Turner et al. [24] model includes: (1) a relativistic jet expansion phase modelled prior to the formation of a lobe, (2) formation of lobes within a surrounding bow shock, and (3) a separation of the ram and thermal components of the jet and lobe pressure.

The Turner et al. [24] model is implemented using the same computational framework as the original RAiSE model discussed in Section 4.1, specifically, using coupled differential equations which are solved for small angular volume elements of the lobe and shocked shell. Below, we present a concise derivation of their solution for the expansion of the relativistic jet (Section 4.3.1), summarise their methodology to model the subsequent lobe formation/inflation (Section 4.3.2), and finally present their method to separate the ram and thermal components of the jet and lobe pressure (Section 4.3.3).

4.3.1. Relativistic Jet Expansion

The relativistic hydrodynamic conservation equations relate the properties of fluids upstream and downstream of a shock discontinuity via the stress-energy tensor. The conservation equations for a relativistic fluid are expressed in terms of comoving quantities including gas density $\rho$, gas pressure p, dimensionless specific enthalpy h (i.e., enthalpy divided by $c^2$), and the non-zero spatial component of the four-velocity $u=\gamma v/c$ (hereafter shortened to four-velocity) relative to the shock front. The conservation equations for a relativistic fluid are (e.g., [24], and references therein):

$$\begin{array}{c}\rho \gamma vs.={\rho }_1{\gamma }_1{v}_1\left(\mathrm{continuity}\right)\end{array}$$

(28a)

$$\begin{array}{c}\rho h{\gamma }^2{v}^2+p={\rho }_1{h}_1{\gamma }_1^2{v}_1^2+p_1\left(\mathrm{momentum}\right)\end{array}$$

(28b)

$$\begin{array}{c}\rho (h\gamma -1)\gamma vs.={\rho }_1(h_1{\gamma }_1-1){\gamma }_1{v}_1\left(\mathrm{energy}\right).\end{array}$$

(28c)

where the fluid downstream of the shock is represented by the subscript ‘1’; no subscript refers to the upstream fluid.

The conservation of energy expression in Equation (28c) is related to the rate of energy input by the jet, Q, by multiplying by the cross-sectional area of the jet (cf. [36], specifically their Equation (26)). That is:

$$Q={\rho }_j(h_j{\gamma }_j-1){\gamma }_j{v}_j{c}^2\mathsf{\Omega }{r}^2,$$

(29)

where the factor of $c^2$ is added to convert the dimensionless enthalpy, $h_j$, to the specific enthalpy, $v_j$ is the bulk velocity of the jet plasma, and ${\gamma }_j$ is the corresponding Lorentz factor. We can, therefore, obtain an expression for the density of the jet plasma some distance r along the jet in terms of the dimensions and energetics of the jet as follows:

$${\rho }_j\left(r\right)=\frac{Q}{{\gamma }_j{v}_j{c}^2(h_j{\gamma }_j-1)\mathsf{\Omega }{r}^2},$$

(30)

where the density at the jet head, ${\rho }_j\equiv {\rho }_j\left(R_s\right)$, is of particular interest for the radio source dynamics.

We now derive the Rankine–Hugoniot jump conditions relating the density and velocity of both the jet plasma and the ambient medium. The bulk velocity of the ambient medium in the observer frame is zero at all times for random particle motions. As a result, the bulk velocity of these particles in the frame of the shock front, $v_1$, is exactly equal to the expansion rate of the shock in the observer frame, $v_s$, i.e., $v_1\equiv -v_s$. By contrast, the bulk velocity of the upstream fluid particles in the jet is non-zero, defined as $v_j$ in the observer frame. Following [24], the conservation of momentum flux equation can, therefore, be rewritten as:

$${\rho }_j{h}_j{\gamma }_j^2{\gamma }_s^2{(v_j-v_s)}^2={\rho }_x{h}_x{\gamma }_s^2{v}_s^2,$$

(31)

where $h_j$ is the dimensionless specific enthalpy of the jet, and ${\rho }_x$ and $h_x$ are the density and dimensionless specific enthalpy of the (external) ambient medium, respectively. Rearranging yields a relationship between the jet-head advance speed and the bulk velocity of the jet (cf. [37,38,39], and subsequent authors):

$$v_s\equiv \frac{d{R}_s}{dt}=\frac{v_j}{1+{[{\rho }_j{h}_j{\gamma }_j^2/\left({\rho }_x{h}_x\right)]}^{-1/2}},$$

(32)

where the dimensionless quantity ${\eta }_R={\rho }_j{h}_j{\gamma }_j^2/\left({\rho }_x{h}_x\right)$ is a function of properties of the jet and ambient medium. That is:

$${\eta }_R\left(r\right)=\frac{Q{h}_j{\gamma }_j}{k{h}_x{v}_j{c}^2(h_j{\gamma }_j-1)\mathsf{\Omega }{r}^{2-\beta }}.$$

(33)

where we have made use of the power-law approximation for the local density of the ambient medium, $\rho =k{r}^{-\beta }$, and Equation (30) for the jet plasma density.

The jet length is found by integrating Equation (32) with respect to time; however, an analytical solution is only possible in the limits ${\eta }_R\to 0$ and ${\eta }_R\to \infty$ (e.g., [40,41]). Turner et al. [24] solved this integral numerically using a fourth-order Runge–Kutta method on the following system of three ordinary differential equations:

$$\begin{array}{cc}\hfill {\dot{R}}_s& =v_s\hfill \\ \hfill {\dot{v}}_s& =\frac{(\beta -2)v_j{v}_s}{2R_s{\eta }_R^{1/2}{[1+{\eta }_R^{-1/2}]}^2}\hfill \\ \hfill {\dot{\gamma }}_s& =\frac{{\gamma }_s^3{v}_s{\dot{v}}_s}{c^2}.\hfill \end{array}$$

(34)

We note that in the interests of clarity, we have omitted the transverse density and velocity structures of the flow along the jet from the above derivation; we refer the reader to Section 2.2.2 of Turner et al. [24] for a complete description.

4.3.2. Lobe Formation

The energy supplied by the central nucleus is initially focussed over a small range of angles within the half-opening angle of the jet. Beyond some lobe formation length scale, the energy must be distributed across the $2\pi$ steradians of the shocked shell. The source expansion in these two phases is described by the differential equations for the relativistic jet (Section 4.3.1) and lobe and shocked shell (Section 4.1). Turner et al. [24] combined these frameworks by modelling the expansion of the radio source as a two-phase fluid, where each angular volume element is assumed to comprise a fraction $\mathsf{\Lambda }\left(t\right)$ of lobe plasma at any given time t. Turner et al. [24] related the acceleration of the ellipsoidal bow shock surrounding the lobe to the acceleration in the jet- and lobe-dominated expansion phases, ${\dot{v}}_{s,jet}$ (Equation (34)) and ${\dot{v}}_{s,lobe}$ (Equation (23)), respectively, as follows:

$${\dot{v}}_s\left(\theta \right)=[1-\mathsf{\Lambda }]{\dot{v}}_{s,jet}\eta \left(\theta \right)+\mathsf{\Lambda }{\dot{v}}_{s,lobe}\left(\theta \right),$$

(35)

where $\mathsf{\Lambda }$ is the fractional contribution of the lobe plasma to the acceleration of the bow shock at a given time. The other two coupled ordinary differential equations (for the velocity and derivative of the Lorentz factor) are identical for both fluids, and thus, do not require any modification.

Turner et al. [24] defined the transition from a jet-dominated to a lobe-dominated flow based on the length scale at which lobe formation commences. This length scale is calculated by equating the densities of the jet plasma and ambient medium (e.g., [10,11]). Turner et al. [24] parametrised the transition from jet- to lobe-dominated expansion by using the ratio of these densities:

$$\mathcal{L}\left(t\right)=\frac{{\rho }_j}{{\rho }_x}=\frac{{\eta }_R\left(R_s(\theta =0,t)\right)}{{\gamma }_j^2}$$

(36)

where ${\eta }_R\left(r\right)$ is defined in Equation (33) and is evaluated for the length of the jet at time t. Turner et al. [24] used this ratio to calculate the fractional contribution of the lobe to source expansion:

$$\mathsf{\Lambda }\left(t\right)=e^{-{\mathcal{L}}^2\left(t\right)/(2log2)},$$

(37)

where $\mathsf{\Lambda }\left(t\right)\to 0$ in the jet-dominated expansion phase and $\mathsf{\Lambda }\left(t\right)\to 1$ in the lobe-dominated phase.

This two-phase fluid model describes the evolution of the bow shock across the transition from a jet-dominated to lobe-dominated flow. A more complete description requires consideration of lobe formation inside the shock front. We refer the interested reader to Section 2.4.2 of Turner et al. [24] for these details.

4.3.3. Thermal Pressure

The Turner et al. [24] relativistic jet model (Section 4.3.1) and their earlier lobe-dominated expansion model (Section 4.1) derive the jet and lobe length evolution by considering conservation of momentum flux (Equation (31)); however, the relative magnitudes of the ram and thermal pressure components after the interaction are not explicitly calculated. These pressure components are difficult to separate directly using the conservation equations; however, we know the lobe evolution is driven entirely by the thermal component in the limit $t\to \infty$. Turner et al. [24], therefore, found the thermal pressure at earlier times by iteratively solving (backwards in time) the relevant differential equations with the initial condition at $t\to \infty$. We refer the interested reader to Section 2.3.4 of Turner et al. [24] for a complete description of the separation of the ram and thermal components of the lobe internal pressure.

5. Discussion

In preceding sections, we have presented the theory underpinning the key classes of analytical models describing the dynamics of kiloparsec-scale radio AGN jets and lobes. The same physical principles are considered in each of these models, notably ram pressure against the ambient medium and an adiabatic equation of state; however, their implementation between model classes differs greatly—as we discuss in Section 5.1. We compare the accuracy of predictions for each model type relative to the outputs of a three-dimensional relativistic hydrodynamic simulation in Section 5.2. We then assess for which regions of parameter space the different model classes yield comparable results, and conversely, those where large differences are expected, in Section 5.3.

5.1. Similarity of Key Model Classes

The four key classes of analytical models examined in this review share common physical principles to explain the dynamics of kiloparsec-scale radio sources (see

Table 1. Summary of the model assumptions for each of the four key classes of analytical models considered in this review. First column: analytical model. Second column: type of solution (analytical or semi-analytic), together with a note on the assumed ambient gas density profile. The third and fourth column: physical process driving the jet and lobe length expansion at early and late times, respectively. The fundamental equations used to derive the source evolution, and the time dependence of the length expansion, are also listed.

Model	Type	Early-Time Evolution	Late-Time Evolution
Scheuer (1974; Model A)	analytical (constant density)	momentum flux (Bernoulli equation) $R\left(t\right)\propto t^{2/(4-\beta )}$	–
Falle (1991)	analytical (power-law density profile)	–	internal pressure (first law of thermodynamics) $R\left(t\right)\propto t^{3/(5-\beta )}$
Hardcastle (2018)	semi-analytic (spherically symmetric density profile)	momentum flux (non-relativistic shock-jump conditions) $R\left(t\right)\propto t^{2/(4-\beta )}$	internal pressure (non-relativistic shock-jump conditions) $R\left(t\right)\propto t^{3/(5-\beta )}$
Turner et al. (2023; RAiSE)	semi-analytic (spherically symmetric density profile)	momentum flux (relativistic hydrodynamics) $R\left(t\right)\propto t$ ${}^{†}$	internal pressure (first law of thermodynamics) $R\left(t\right)\propto t^{3/(5-\beta )}$

† Solution for jet-head length expansion prior to lobe formation length scale.

). Scheuer [17] modelled the forward expansion of the source based on the momentum flux of the jet and invokes internal energy to calculate the sidewards expansion of the lobe. Falle [18], instead, modelled the forwards expansion by considering the adiabatic expansion of the lobe due to an increase in internal energy while using the jet momentum flux to relate the shape of the lobe to the opening angle of the jet. Meanwhile, the Hardcastle [23] and Turner et al. [24] models smoothly transition their dynamics between the jet- and lobe-dominated expansion phases predicted by these earlier models.

5.1.1. Early-Time Evolution

The source length evolution in the Scheuer [17] model, and the (jet-dominated) early-time expansion phases of the Hardcastle [23] and Turner et al. [24] models, are derived considering the relative amplitudes of the momentum flux of the jet and the thermal pressure of the ambient medium (Equation (28b)). The relativistic hydrodynamic equations used in the theory of Turner et al. [24] can be simplified to obtain the expressions found by the Scheuer [17] class of models. In particular, the jet-head advance speed derived in Equation (32) is integrated to yield the source length as follows:

$$\begin{array}{cc}\hfill R\left(t\right)& ={\int }_0^t\frac{d{R}_s}{dt}dt\hfill \\ \hfill & \approx {\left(\frac{Q{h}_j{\gamma }_j{v}_j}{\mathsf{\Omega }k{h}_x{c}^2(h_j{\gamma }_j-1)}\right)}^{1/(4-\beta )}{\left(\frac{(4-\beta )t}{2}\right)}^{2/(4-\beta )},\hfill \end{array}$$

(38)

where the second equality is valid for ${\eta }_R\left(r\right)\ll 1$, which coincides with the formation of lobes on a length scale of order 1 kpc [24]. This equation converges to that proposed by Rees [16] and Scheuer [17] (their Model A) by taking their limit of massless particles with velocity c (i.e., $h_{j,x}\to 1$ and ${\gamma }_j\to \infty$):

$$R_s\left(t\right)={\left(\frac{Q}{\mathsf{\Omega }kc}\right)}^{1/(4-\beta )}{\left(\frac{(4-\beta )t}{2}\right)}^{2/(4-\beta )},$$

(39)

where for their assumption of a uniform ambient medium, we set $\beta =0$ and $k=\rho$.

The similarity of the early-time evolution predicted by the Hardcastle [23] model and Scheuer [17] class of models is immediately apparent by comparing their expression for the jet-head pressure (Equation (24)). In the limit $t\to 0$, their expression converges to the ram pressure component as follows:

$$\begin{array}{cc}\hfill p_s(\theta =0)& =\frac{\epsilon Q{R}_s(\theta =0)}{2cV}\left(\mathrm{Hardcastle}\right)\hfill \\ & =\frac{{\kappa }_{1}Q}{\mathsf{\Omega }{R}_s^2(\theta =0)c}\left(\mathrm{Scheuer}\right),\hfill \end{array}$$

(40)

where the geometric factor assumed by Hardcastle [23] is identified as $\epsilon =2{\kappa }_1/\mathsf{\Omega }$, using the terminology of Scheuer [17]; in other words, this geometric factor largely corresponds to the solid angle of the jet and the jet-head region. Hardcastle [23] derived the jet-head advance speed using the non-relativistic Rankine–Hugoniot shock jump conditions (Equation (27)), largely equivalent to the ram pressure argument employed by Scheuer [17]. Because of this, the source length evolution predicted by the Hardcastle [23] model early in the source lifetime matches that of the Scheuer [17] class of models for the same input parameters, i.e., $R\propto t^{1/2}$ for a constant density ambient medium. By contrast, the jet-dominated expansion phase of the Turner et al. [24] model yields a very different evolutionary history to these two model types, with $R\propto t$ prior to lobe formation, as those authors do not take the non-relativistic limit of the hydrodynamic equations.

5.1.2. Late-Time Evolution

The lobe length expansion in the Falle [18] class of models and the (lobe-dominated) late-time expansion phases of the Hardcastle [23] and Turner et al. [24] models are calculated by considering the internal energy of the relativistic lobe plasma evolving under an adiabatic equation of state. The first law of thermodynamics (Equation (12)) applied to a self-similar ellipsoidal shell yields lobe length growth of the form $R\propto t^{3/5}$ (Equation (18), constant gas density form is stated here for simplicity), consistent with findings for supernova remnants. The Falle [18] class of models use this dependence directly, while the lobe-dominated expansion phase of the Turner et al. [24] model also permits an evolving lobe axis ratio due to a non-power law ambient medium. The late-time evolution of the Turner and Shabala [22] model (and subsequent versions of RAiSE) ultimately transitions from the supersonic to subsonic regime as the jet-head advance speed slows ($R\propto t^{1/3}$). This is not critical when modelling powerful lobed radio sources which drive strong shocks—the focus of this review—but is essential when considering the coasting (inactive) phase of remnant radio sources.

Hardcastle [23] did not make any explicit assumptions about the lobe axis ratio; however, in the limit $t\to \infty$, their expressions for the expansion rate along both the major and minor axes converge, leading to a spherical lobe (at least in the special case of a constant density ambient medium). The expression for their growth rate at late times (Equation (27)) becomes:

$$\frac{d{R}_s}{dt}=\sqrt{\frac{({\mathsf{\Gamma }}_s+1)([\xi {\mathsf{\Gamma }}_c+(1-\xi ){\mathsf{\Gamma }}_s-1]Qt+f(N,T,v_s))}{2{\rho }_x{V}_s}},$$

(41)

where we have assumed that the pressure along both axes is dominated by the internal energy component (Equation (24)) and that the sound speed in the shocked gas is comparable to that of the ambient medium (i.e., $c_s\sim \sqrt{{\mathsf{\Gamma }}_s{p}_x/{\rho }_x}$), and applied the relationship for the ratio of the lobe and shocked shell volumes (Equation (25)). This differential equation again yields a solution of the form $R\propto t^{3/5}$ for a constant density ambient medium, if the internal energy contribution from the shocked gas, $f(N,T,v_s)$, is ignored. Such an assumption is reasonable for high jet powers or moderate ambient densities; however, the late-time evolution diverges significantly from this prediction for weak jets or high density environments—we examine this point further in Section 5.2.2.

5.1.3. Source Morphology

The lobe volume and, hence, axis ratio are calculated using the remaining equations not previously invoked in the calculation of the source length expansion. For example, in the Scheuer [17] model, source expansion along the jet axis is modelled based on ram pressure arguments, while the sidewards expansion is derived considering the internal energy of the injected lobe plasma; Hardcastle [23] made similar arguments, while Turner et al. [24] invoked internal energy to describe the formation and inflation of their lobe within the confines of a surrounding shocked gas shell. By contrast, the self-similar expansion assumed in the Falle [18] class of models, and the lobe-dominated expansion phase of the Turner et al. [24] model, implicitly sets the lobe volume based on its length, which is calculated from internal energy. These models use ram pressure arguments to relate the internal conditions of the lobe to those of the ambient medium (Equation (11)), and to relate the jet half-opening angle to the source axis ratio.

5.2. Comparison to Hydrodynamic Simulations

To test the analytical models, we compare the dynamics of the four key model classes introduced above to the results of hydrodynamic simulations run using the PLUTO code [42]. Below, we describe the existing simulations run by Yates-Jones et al. [43] for powerful FR-II radio galaxies and test the analytical models by comparing the predicted evolution of the source length, lobe axis ratio, and jet-head pressure throughout the source lifetime to the simulation results.

5.2.1. Hydrodynamic Simulation Dynamics

The hydrodynamic simulations of Yates-Jones et al. [43] consider an initially conical jet of half-opening angle ${\theta }_j=10$ degrees and bulk flow with Lorentz factor ${\gamma }_j=5$. Their high-powered jet ($Q=3\times {10}^{38}$ W) expands into a spherically symmetric King profile with a core density of ${\rho }_c=2.41\times {10}^{-24}$ kg m${}^{-3}$, core radius of $r_c=144$ kpc, and slope described by the coefficient ${\beta }^{\prime }=0.38$ (for details, see [31]). Those authors considered both jets located at the centre of the gas distribution, as well as offset jets. In this review, we will only consider their cluster-centred jet simulation, as the theory underpinning the analytical models assumes a spherically symmetric environment. Their simulations resulted in lobed Fanaroff and Riley [12] Type-II sources forming on a timescale of ∼1 Myr, and considered the late-time evolution up to an age of 35.1 Myr.

We extract time series for the evolution of the source length, the lobe axis ratio, and the jet-head pressure from the hydrodynamic simulation outputs (for details, see [24]). These are critical dynamical quantities in calculation of both the source evolution and the radio-frequency synchrotron emission, and thus, should be considered to assess the accuracy of the analytical models. The potentially more informative lobe volume and volume-weighted pressure are poorly constrained prior to lobe formation as large regions near the core remain partially occupied by ambient gas; the calculation of these quantities in the hydrodynamical simulation is highly dependent on the threshold used to separate ambient gas from the jet plasma.

5.2.2. Accuracy of Analytical Models

The critical intrinsic parameters characterising radio source evolution in analytical models are the jet kinetic power, source age, and properties of the ambient medium. The spherically symmetric King profile used by Yates-Jones et al. [31] for the ambient medium in their simulations is readily modelled by both the Hardcastle [23] and Turner et al. [24] models, but not the older models. The original form of the Scheuer [17] model assumes a constant density ambient medium, while the Falle [18] model employs a slightly more general power-law of the form $\rho =k{r}^{-\beta }$; in Section 2, we similarly derived the Scheuer [17] model for a power-law gas density profile. To facilitate a meaningful comparison between all four model classes, we modify the Falle [18] and Scheuer [17] model following the approach of Turner and Shabala [22], by approximating the gas density profile as a series of contiguous power-laws with piecewise solutions. The implementation of these models for a general ambient medium is available on our GitHub online repository1.

Source Length

The source length evolution for the four analytical models is shown in the top panel of Figure 5. The jet power, source age, and ambient gas density profile are in all cases identical to the inputs to the hydrodynamic simulation of Yates-Jones et al. [31]; however, some of the more minor model parameters are varied to obtain the best representation of each model class. Specifically, the Scheuer [17] and Falle [18] model evolutionary tracks are shown for three plausible values of the jet half-opening angle ${\theta }_j$, while the Hardcastle [23] model is shown for three values of their equivalent geometric factor $\epsilon$. The free parameters in the Turner et al. [24] model have previously been calibrated based on this hydrodynamic simulation, and thus, results for only a single set of parameters for this model are shown. The resulting evolutionary tracks are consistent with the discussion in Section 5.1: the Scheuer [17] and (jet-dominated) early-time Hardcastle [23] models follow an approximately $R\propto t^{1/2}$ growth rate (the dependence expected for a flat atmosphere), while the Falle [18] and (lobe-dominated) late-time Hardcastle [23] models follow an $R\propto t^{3/5}$ expansion rate (also see Figure 6). By contrast, the relativistic hydrodynamic equations used by Turner et al. [24] yield quite different evolutionary tracks in the jet-dominated expansion phase, and are more consistent with the hydrodynamic simulation. At late times (>10 Myr), the Turner et al. [24] model predicts slower expansion close to an $R\propto t^{1/2}$ relationship, converging towards the same limit as the other models at later times. We explore this in more detail in Section 5.3.

Lobe Axis Ratio

The evolution of the lobe axis ratio is shown in the middle panel of Figure 5. The Turner et al. [24] model captures axis ratio evolution prior to lobe formation, as it considers a two-phase fluid with an initially low, but non-zero, fraction of jet plasma in the region between the jet and bow shock. By contrast, the other analytical models initially disagree with the hydrodynamic simulation, as they assume a plasma-filled lobe structure throughout the source lifetime. Upon lobe formation, the Hardcastle [23] model agrees well with the hydrodynamic simulation assuming a compact jet-head region with a geometric factor $\epsilon \gtrsim 40$ (corresponding to a jet head of radius 12.6 kpc for a typical 100 kpc jet). The Falle [18] class of models, of course, yields a constant lobe axis ratio (which is why these models are often referred to as “self-similar”) throughout the source evolutionary history, while the Scheuer [17] model predicts a rapidly increasing axis ratio as the lobe encounters the steeper sections of the ambient gas density profile ($A\propto t^{(5\beta -2)/(16-4\beta )}$; see Section 2). These last two models do not explicitly separate the lobe and shell material, so for a fairer comparison, we are guided by the numerical results of Turner and Shabala [32] in assuming that the lobe axis ratio scales to that of the shell as $A=A_s^{1.7}$.

Jet-Head Pressure

We finally compare the predicted evolution of the jet-head pressure to that measured from the hydrodynamic simulation (Figure 5, bottom panel). The majority of the evolutionary history probed by the hydrodynamic simulation (≲10 Myr) is associated with a significant ram pressure component at the jet head, in addition to a thermal component that scales approximately in proportion to the ram pressure component. The Scheuer [17] and Hardcastle [23] models directly consider the ram pressure component, and hence, predict a power-law evolution of jet head pressure ($p\propto t^{-1}$ for a constant gas density ambient medium; see Section 2 and Section 4.2), which is broadly similar to the thermal jet head pressure component measured by the hydrodynamic simulation. The self-similar expansion model of Falle [18] yields a flatter relationship of the form $p\propto t^{-4/5}$. Turner et al. [24], instead, explicitly modeled both the thermal and ram pressure components throughout the source evolutionary history; this model accurately captures the steepening in the rate of change of jet-head pressure during the transition between these two limiting cases.

5.3. Parameter Space Exploration

In this section, we compare the consistency of observable predictions between the four model classes. We select three parameters for our comparison. Source length and axis ratio are directly measurable model predictions, while synchrotron luminosity integrated over the entire lobe is a quantity which can be approximated from source dynamics—noting that detailed consideration of particle acceleration and loss processes is required for a full calculation (e.g., [44]). These radio source attributes are critical for estimating the jet energy budget through a parameter inversion of observables [22,45].

We investigate the behaviour of the different model classes for a range of input parameters, specifically the single-jet kinetic power Q, core density ${\rho }_0$, and scale radius $r_c$ of the ambient gas density profile. The base-case set of parameters is chosen to match the hydrodynamic simulation of Yates-Jones et al. [31], but with ${\rho }_0={10}^{-23}$ kg m${}^{-3}$ and $r_c=100$ kpc. The jet half-opening angles for the Scheuer [17] and Falle [18] models are set to ${\theta }_j={10}^{\circ }$ and ${\theta }_j={20}^{\circ }$, respectively, as these closely match the hydrodynamic simulation evolution for the three key dynamical parameters of source length, lobe axis ratio, and jet-head pressure (cf. Figure 5). For the same reasons, we set the geometric factor to $\epsilon =40$ in the Hardcastle [23] model. The evolutionary tracks for each model class are considered both for the base set of parameters and for a factor of ten variation to one of either the jet power ($Q=3\times {10}^{37}$, $3\times {10}^{38}$ or $3\times {10}^{39}$ W), core density (${\rho }_0={10}^{-24}$, ${10}^{-23}$ or ${10}^{-22}$ kg m${}^{-3}$), or scale radius ($r_c=10$, 100 or 1000 kpc).

Source Length

The source length evolution for the four model classes are shown in Figure 6. Changes in jet power and core density largely result in a constant offset (in logarithmic space) to the source length, consistent with the expected scalings in flat atmospheres of $R\propto Q^{1/4}$ and $R\propto {\rho }_0^{-1/4}$ for the jet-dominated [17] model, and $R\propto Q^{1/5}$ and $R\propto {\rho }_0^{-1/5}$ for the lobe-dominated [18] model class. The more sophisticated models of Hardcastle [23] and Turner et al. [24] capture the transition between these limiting cases. For the same input parameters, the Hardcastle [23] and Turner et al. [24] models predict similar dynamical evolution for the majority of simulated sources. However, the jet-dominant phase lasts longer for high-powered jets and in denser environments in the Turner et al. [24] model, making this model more sensitive to changes in these parameters than the other models. Variations to the scale radius of the ambient gas density profile yield qualitatively similar behaviour between the model classes, with faster expansion occurring for steeper atmospheres.

Lobe Axis Ratio

The evolution in the lobe axis ratio is shown in Figure 7. The Hardcastle [23] and Turner et al. [24] models both predict the lobe axis ratio evolutionary tracks shift horizontally along the source age axis in response to variations in the jet power and core density, with lobe formation occurring earlier for less powerful jets and/or denser ambient media (e.g., see [24]). The Scheuer [17] model predicts a similar response to changes in the jet power and core density for the approximately constant ambient gas density section of the lobe axis ratio evolutionary tracks (i.e., ≲1 Myr), but is inconsistent at later times once the ambient density profile begins to steepen (see Section 5.2). By contrast, variations to the scale radius produce lobe axis ratio evolution that is largely inconsistent between the different model classes. The self-similar model of Falle [18] yields a constant lobe axis ratio for all input parameters by construction.

Synchrotron Luminosity

Radio galaxies are detectable through their synchrotron emission. In analytical radio source models, this emission is typically calculated by assuming a scaling between the lobe pressure and magnetic field, acceleration of emitting particles at the jet termination shock, and subsequent losses due to adiabatic, synchrotron, and inverse-Compton losses due to upscattering of cosmic microwave background photons.

A full calculation of the radio-frequency luminosity of each lobe is beyond the scope of this review (see, e.g., [22,23,44,45]). However, a useful estimate can be obtained by considering only adiabatic losses, which are directly related to the evolution of lobe pressure. The “lossless” luminosity calculated in this approach is related to the lobe volume and pressure as $L_{\nu }\propto p^{(\alpha +3)/2}V$, where $\alpha \sim 0.7$ is the spectral index of the non-aged radio spectrum. The total internal energy of the lobe, $U\sim pV$, represents the fraction of the jet energy that is transferred to the lobe, and hence, is very similar in all models for a fixed time and jet power. In this case, the synchrotron luminosity scales with pressure as $L_{\nu }\propto p^{(\alpha +1)/2}U$.

The evolution in the lossless synchrotron luminosity for each of the model classes is shown in Figure 8 at a rest-frame frequency of 151 MHz. The predictions of the Scheuer [17] model are once again inconsistent with the later models, with radio luminosities up to a factor of 100 greater than the other models at late times due to their higher lobe axis ratios, and subsequently, smaller volumes and higher pressures. The remaining three models are consistent at late times—when their lobe pressures are largely derived based on changes in lobe internal energy—for all sets of input jet powers, core densities, and scale radii. The Falle [18] and Hardcastle [23] model classes also agree at early times, in contrast to the Turner et al. [24] model, which predicts higher luminosities as a result of significantly higher pressures prior to lobe formation when the system is dominated by the momentum flux of the relativistic jet. At these early times, the luminosities are a factor of 100–1000 higher than in the lobe-only models of Falle [18] and Hardcastle [23].

6. Concluding Remarks

We have summarised and compared the main classes of analytical models describing the dynamics of kiloparsec-scale lobed radio galaxies. These models can be separated into two main classes, depending on whether the expansion of the radio source is driven by the momentum flux from the jet or by the internal lobe pressure. We presented the Scheuer [17] and Falle [18] models, respectively, to describe the general characteristics of other literature models in either of these two classes. We also examined separately the more recent models proposed by Hardcastle [23] and Turner et al. [24], which combine aspects of both jet momentum flux and lobe pressure.

We compared the different model classes against each other, and with high-resolution hydrodynamic simulations, for a range of realistic input parameters. Our key findings are as follows:

• Jet momentum flux and lobe internal pressure dominate the early- and late-time radio source evolution, respectively. Both must be considered for a complete radio source model describing source dynamics after the lobe formation phase (∼1 Myr; Section 5.2.2).
• Realistic ambient gas density profiles (i.e., not constant or power-law) produce radio sources which are inconsistent with the self-similar lobe evolution predicted by the Falle [18] class of models (Section 4.1 and Section 4.2). This naturally explains the large axis ratios seen in giant radio galaxies.
• Relativistic jet dynamics is important for an accurate description of early source evolution, before the lobe formation phase (Section 5.2.2).

We make three of the four models considered in this review openly available. The code extending the Scheuer [17] and Falle [18] models to general atmospheres,2 and RAiSE3 (version 2023), are available in our GitHub online repository. Hardcastle [23] has also made their code available; we refer the interested reader to their paper.

We conclude with a brief reflection on the next generation of analytical models. Existing analytical models neglect the interaction between the jet and the multi-phase interstellar medium of the host galaxy. Hydrodynamic simulations that model this interaction predict that the jet can spend $\gtrsim 1$ Myr in the galaxy (e.g., [46]). The majority of observed radio sources are compact and short-lived [47,48], and hence, these processes are likely to be relevant to the bulk of the radio source populations. At the other end of the radio galaxy evolution scale, Hardcastle [23] pointed out that for extreme losses, such as expected in large sources at high redshift, it is possible for the majority or even all of the jet energy to be radiated away. This mechanism can potentially limit the maximum size to which a radio source can grow. Existing analytical models decouple source dynamics from the synchrotron and inverse Compton radiative loss mechanisms, and hence, are not currently capable of tackling this issue.