Dynamics of relativistic radio jets in asymmetric environments

A key feature of this modelling is the observationally motivated environments used as the initial conditions for the simulations. Our focus is on low-redshift clusters, which are reasonably well described by the isothermal beta profile (Cavaliere & Fusco-Femiano, 1978) for gas density as a function of radius from the core $\rho(r)$, where

This profile is parameterised by the core density $\rho_{0}$, core radius $r_{\textrm{c}}$, and power-law slope $\beta$. In this work, we use a density profile typical of a cluster environment, with $M_{\textrm{halo}}=3\times 10^{14}\,\textrm{M}_{\sun}$ and $\beta=0.38$. We set the core radius to be $0.1R_{\textrm{vir}}$, or $r_{c}=144\,\textrm{kpc}$, following Yates et al. (2018); observationally the core radii of clusters vary between tens and hundreds of kpc (Vikhlinin et al., 2006). These parameters give a central density of $\rho_{0}=1.4\times 10^{-27}\,\textrm{g cm}^{-3}$. In Figure 1 we show the environment density profile as a function of distance along the jet axis, for the different offset and inclination angles used in this work. Jets launched offset from the cluster centre experience significant asymmetry in the environment density profile; in Sect. 3.1 we argue that this asymmetry has a large impact on the resulting jet dynamics. We also note that the impact of inclination angle on density is less pronounced, however, an inclined jet does experience a gravitational acceleration asymmetry which affects remnant morphology.

The pressure profile is calculated from the density profile using $P(r)=\frac{k_{\textrm{B}}T\rho(r)}{\mu m_{\textrm{H}}}$. We take $\mu=0.60$ and $T=3.46\times 10^{7}\,\textrm{K}$, consistent with observed temperatures in clusters. Cooling is not included in our simulations. This is a valid assumption provided the environment cooling time is long; the central cooling time $t_{\textrm{cool}}$ for the environment presented here is on the order of $1.5\,\textrm{Gyr}$ for a metallicity of $Z=-1.0$ (Sutherland & Dopita, 1993), more than an order of magnitude longer than the jet lifetimes considered in this work.

The radial gravitational acceleration necessary to keep the environment in hydrostatic equilibrium is

where $c_{\textrm{s}}=\left(\frac{\Gamma P}{\rho}\right)^{1/2}$ is the environment sound speed. This is implemented as a source term in pluto. The stability of the environment is confirmed by evolving it for over $100\,\textrm{Myr}$, twice as long as the maximum jet time.

We offset the environment in both the Y and Z directions to study the parameter space detailed in Sect. 2.3. This is achieved through a transformation of coordinates with respect to the new environment centre, resulting in an offset radius given by

Version 4.3 of pluto is used to carry out the simulations presented here. The relativistic hydrodynamics physics module is used, along with the hllc Riemann solver, linear reconstruction, second-order Runge-Kutta time-stepping, the Taub-Mathews equation of state, and a Courant-Friedrichs-Lewy (CFL) number of $0.33$.

At each timestep pluto integrates the system of conservation laws described as

with conservative state variables $D,\boldsymbol{m},E$ (the laboratory density, momentum density, and total energy density respectively), fluid three-velocity $\boldsymbol{v}$, and pressure $p$.

The simulations were carried out on a static three-dimensional Cartesian grid centred at the origin with a side length of $l_{\textrm{x,y,z}}=400\,\textrm{kpc}$, cell count of $(n_{\textrm{x}},n_{\textrm{y}},n_{\textrm{z}})=(550,550,960)$, and reflective boundary conditions. A central uniform grid patch of $100$ cells is defined around the injection region in all three dimensions ($-2.5\rightarrow 2.5\,\textrm{kpc}$, a resolution of $0.05\,\textrm{kpc}$) to ensure that injection of the initially conical jet and its hydrodynamic collimation by ambient pressure are sufficiently resolved. Either side of the central grid patch is a geometrically stretched grid consisting of $430$ cells in Z, and $225$ cells in X and Y; giving a minimum resolution of $1.7\,\textrm{kpc}$ and $2.8\,\textrm{kpc}$ per cell respectively. Reflective boundary conditions were applied to both the lower and upper grid boundaries.

A two-dimensional spherical axisymmetric grid is used in a resolution study to verify that jet recollimation is captured correctly. In these simulations the jet is aligned with the axis of symmetry, and the simulation domain is taken to be $(r,\theta)=(0.2\,\textrm{kpc},0\degr)\rightarrow(800\,\textrm{kpc},180\degr)$. A high-resolution uniform radial grid patch of $256$ cells covers the injection region ($r=0.2\rightarrow 2\,\textrm{kpc}$, a resolution of $7\,\textrm{pc}$ per cell), while the remainder of the domain is uniformly covered by $1650$ cells. Azimuthally the jet injection region is uniformly resolved for each side by $300$ cells over the first $5\degr$, and $250$ cells over the next $12\degr$. A lower resolution uniform grid patch of $400$ cells is placed from $\theta=17\degr\rightarrow 163\degr$, as for the adopted jet half-opening angle of $15\degr$ this region will mainly contain sideways expansion of the jet lobe and cocoon. At a distance of $200\,\textrm{kpc}$ this provides a resolution of $58\,\textrm{pc}$ across the jet head. The lower radial boundary has an inflow condition within the jet nozzle; reflective boundary conditions were applied outside the jet nozzle and at the upper radial boundary. Axisymmetric boundary conditions were applied to both the lower and upper azimuthal boundaries.

The jet injection region is defined in three dimensions as a sphere centred at the origin with radius $r_{0}=1\,\textrm{kpc}$ (Figure 2). Within this sphere the jet is injected as momentum and energy fluxes by overwriting the density, pressure, and velocity of these cells with the corresponding injection zone values, $(\rho,P,\boldsymbol{v})=(\rho_{\textrm{i}}(r),P_{\textrm{i}}(r),\boldsymbol{v}_{\textrm{i}}(r))$. The jet density and pressure are defined in the injection sphere as

$P_{\textrm{j}}$ is tuned to keep perturbations from the injection zone small. The velocity in the injection region is defined along the jet propagation axis as $v_{r}=v_{\textrm{j}}$ if $\theta\leq\theta_{\textrm{j}}$, and $0$ otherwise; here, $\theta_{\textrm{j}}$ is the half-opening angle of the jet. A passive tracer fluid is used to trace jet material; it is initialized to $0.0$ throughout the environment while given a value of $1.0$ inside the injection region. This fluid is then advected along with the jet flow and is used to quantify mixing between the jet and ambient material. At a radius of $r_{0}$, the edge of the jet injection zone is resolved by $6$ cells. In all simulations, the jet expands before collimation, so the number of cells across the collimated jet is larger than that of the injection region, providing sufficient resolution across the jet beam to capture recollimation dynamics.

For each one of the two relativistic jets of a bipolar radio source, the total power is given (Mukherjee et al., 2020) by

where $\gamma=1/\sqrt{1-v_{\textrm{j}}^{2}/c^{2}}$ is the bulk Lorentz factor of the flow, $c$ is the speed of light in a vacuum, $\Gamma$ is the adiabatic index, and $A_{\textrm{j}}$ is the cross-sectional area of the jet inlet. The jet density $\rho_{\textrm{j}}$ is calculated using Eq. 7 for a given combination of $Q,v_{\textrm{j}},P_{\textrm{j}},A_{\textrm{j}}$ as

When discussing relativistic fluid flows, it is useful to introduce a temperature parameter, $\Theta=P/(\rho c^{2})$ (Mignone & McKinney, 2007), and in the case of relativistic jets specifically, $\chi$, the ratio between rest-mass energy and the thermodynamic part of the enthalpy (Bicknell, 1994, 1995),

Cold jets with $\Theta\ll 1$ ($\chi\gg 1$) are dominated by kinetic energy, while hot jets with $\Theta\gg 1$ ($\chi\ll 1$) are dominated by thermal energy. When calculating the injection parameters of cold jets ($\chi\gg 1$), as is the case for the jets presented in this work, the ideal equation of state can be used (Mignone & McKinney, 2007). However, initially cold jet material will not remain kinetically dominated throughout the jet evolution due to shock-heating along the jet and at the hotspots. To account for this we use the Taub-Mathews (TM) (Taub, 1948; Mathews, 1971; Mignone & McKinney, 2007) equation of state, which is built into pluto and handles both non-relativistic and relativistic gases. In the TM equation of state, the relativistic specific enthalpy is given as

with an equivalent adiabatic index

The relativistic generalisation of density contrast (Martí et al., 1997; Krause, 2005; Bromberg et al., 2011) is defined as

with relativistic specific enthalpy $h=1+\Gamma e/\rho c^{2}$ for a cold jet given the equation of state for an ideal gas, $p=(\Gamma-1)e$. $\eta_{\textrm{r}}$ can be compared with the non-relativistic jet density contrast $\eta$ such that for a light one-dimensional jet at velocity $v_{\textrm{j}}$, the hotspot would propagate at a velocity of $v_{\textrm{HS}}=\sqrt{\eta_{\textrm{r}}}v_{\textrm{j}}$.

We simulate a radio jet pair, each with a power of $Q=3\times 10^{38}\,\textrm{W}$, typical of moderate-power FR IIs (Turner et al., 2018b; Shabala et al., 2020). The main science simulations consist of a pair of relativistic jets on a three-dimensional grid. We vary the offset from the cluster centre ($0$ or $1$ core radii, corresponding to $0$ or $144\,\textrm{kpc}$), as well as the inclination angle with respect to the centre ($0\degr$, $15\degr$, or $45\degr$). The jets are relativistic, $\gamma=5$ ($v_{j}=0.98c$), have a half-opening angle $\theta_{\textrm{j}}=10\degr$, and are initially overpressured and underdense compared to the ambient medium. Finally, we conduct a resolution study with high-resolution two-dimensional relativistic simulations to verify that we are capturing the jet recollimation dynamics correctly in the three-dimensional simulations.

The simulations and their parameters are listed in Table 1. To aid with the discussion, in all offset cases we define the jet directed towards the cluster centre as the primary jet, while the jet directed away from the cluster centre as the secondary jet.

The simulations presented here were carried out on the kunanyi (Tasmanian Partnership of Advanced Computing) and Raijin (National Computational Infrastructure) high-performance computing facilities. Between $1400$ and $4800$ cores were used, and the three-dimensional relativistic simulations each took approximately $300,000$ CPU hours.

Emission due to synchrotron radiation is the main method by which the large-scale structure of radio sources can be observed. Synchrotron radiating particles are likely accelerated both along the jet and at the jet hotspots, continuing to radiate as they flow back into the radio lobes (Hardcastle et al., 2016; Matthews et al., 2020; Rieger & Duffy, 2021). Thus to compare observed radio source structure with simulations, the synchrotron emissivity of a simulated radio source needs to be calculated. While the simulations presented here do not contain the magnetic fields necessary for a complete treatment of synchrotron radiation, the emissivity can be calculated by assuming a constant departure from equipartition for the magnetic field and electron energy densities (Kaiser et al., 1997; Turner & Shabala, 2015). Magnetohydrodynamic simulations show that the ratio of particle to magnetic pressure (plasma $\beta$) remains on average constant over time (Hardcastle & Krause, 2014), but has large spread spatially due to lobe turbulence (Gaibler et al., 2009). Thus our assumption has the effect of smoothing out any local magnetic field variations across the lobes.

The synchrotron emissivity of a packet of electrons with a power-law distribution of electron energies $N(E)=\kappa E^{-q}$ emitting at an angular frequency $\omega$ is given by Longair (2011) as

The local magnetic field strength $B$ is related to the local electron density $u_{\textrm{e}}$ through a constant departure from equipartition, $f_{\textrm{B}}=u_{\textrm{B}}/u_{\textrm{e}}$. Kaiser et al. (1997) present a relationship between cocoon pressure and energy densities, $P=(\Gamma_{c}-1)(u_{\textrm{e}}+u_{\textrm{B}}+u_{\textrm{T}})$, which is used to relate the local pressure to the local electron density, assuming a negligible thermal energy component. The electron distribution normalisation $\kappa$ depends on observationally-informed distribution properties (electron Lorentz factor limits $\gamma_{\textrm{min}}$ and $\gamma_{\textrm{max}}$, energy power-law slope $q$), and local pressure $P$. The synchrotron emissivity for a given frequency $\nu$ and local pressure $P$ can then be written as

where the full derivation for $j_{0}$ in units of W Hz^-1, the volume emissivity coefficient, is given in Yates et al. (2018).

Spectral ageing is not included in this work, so the synthetic radio spectrum has an identical shape for all frequencies; frequency serves only to change the normalisation $j_{0}$. The emissivity is weighted by the jet tracer fluid, to account for mixing of jet and ambient material. We note that this method for calculating the emissivity does not include any losses, which alter the observed morphology of the radio source (Turner et al., 2018a) and the evolution of luminosity with time (Kaiser et al., 1997; Shabala & Godfrey, 2013; Turner & Shabala, 2015; Hardcastle, 2018). In this work, we focus on the dynamics of the radio source, and so neglect both radiative and adiabatic losses, in addition to Doppler boosting due to relativistic fluid motions. English et al. (2016) showed that the effect of Doppler boosting on total luminosity is small due to both the small amount of emission associated with the jets and the relatively low bulk Lorentz factors found on the grid, consistent with our simulations.

Using Eq. 14, the emissivity per unit volume for each simulation grid cell can be calculated. We take the electron energy power-law slope to be $q=2.2$ as in (Hardcastle & Krause, 2013), corresponding to a spectral index $\alpha=-0.6$ typical of radio lobes; the ratio of magnetic to particle energy densities as $f_{\textrm{B}}=0.1$; and the minimum and maximum electron Lorentz factors as $\gamma_{\textrm{min}}=10$ and $\gamma_{\textrm{max}}=10^{5}$. The total lobe luminosity $L_{\nu}$ is an integral of emissivity over the lobe volume,

The two-dimensional surface brightness can be obtained by integrating along a chosen line of sight $r$,

Following Yates et al. (2018), we place our simulated radio sources at $z=0.05$, interpolate the surface brightness maps onto a pixel size of $1.8\,\textrm{arcsec}^{2}$, and approximate the effects of an observing beam by convolving the surface brightness maps with a two-dimensional circular Gaussian beam with a full width at half maximum of $5\,\textrm{arcsec}$. The observing frequency is chosen to be $\nu=1.4\,\textrm{GHz}$. These parameters are characteristic of the VLA FIRST (Becker et al., 1995) survey, and similar to ongoing and future mid-frequency surveys by SKA pathfinders. The radio sources are placed in the plane of the sky so that the line-of-sight integral is performed along the x-axis. Our synthetic emissivity model captures the hotspot and recently shocked lobe emission with reasonable accuracy, but also produces equatorial emission which is not observed in real radio sources due to electron aging; a detailed discussion of this point is deferred to Yates-Jones et al. (in preparation).

We show logarithmic density maps of our main simulations at the end of the simulation time in Figure 3. The jets re-collimate successfully before entraining gas later via three-dimensional instabilities and transitioning to turbulence. The jet disruption is very similar to the mechanism described in Massaglia et al. (2016), and is essentially due to the low Mach number in the collimated jet. When the jet first collimates, the internal relativistic Mach number is above $10$. Along the jet, it then oscillates between approximately $6$ and $10$, with a declining trend of the minima. Where it drops below $5$ (shown as the white contour in Figure 3), the jet disrupts. We find that the jet disruption radius is roughly equal for a primary/secondary jet pair in a given environment over the simulated time. Additionally, no significant differences are present between environments.

In Figure 4 we show midplane slices of velocity along the jet axis. The white contours (which denote a bulk velocity of $0.5c$) highlight the deceleration of jet material after the disruption point. Clear backflow is present in all simulations, and slight asymmetries are introduced in the inclined jets due to the environment (the lower two panels of Figure 4). The overall lobe shape is pinched in the case of the primary jet (see, e.g., top-right panel of the same figure).

We measure the lobe morphology using a tracer threshold of $10^{-6}$. Gaibler et al. (2009) have shown that lobe morphologies do not strongly depend on the exact value of this parameter. Figure 5 shows the lobe length evolution as a function of time. We also plot the lobe length evolution as modelled by RAiSE (Turner & Shabala, 2015), using the same jet parameters and environment. The jets undergo a fast expansion phase for the first few Myrs, before beginning to slow; this initial ‘breakout’ phase is not captured by the RAiSE model. However, at later times, the rate of jet length evolution is consistent with the analytical model.

In all offset cases, the primary jet is shorter—as predicted by analytic theory (e.g. Kaiser & Alexander, 1997) for a rising density profile—and the lobe is pinched compared to the secondary jet. This causes the observed lobe length asymmetry; initially, the primary and secondary lobes are similar, but they diverge as the source ages. Inclining the jet with respect to the cluster centre has no significant effect on lobe length.

In Figure 6 we plot the evolution of the length and width ratios for all three-dimensional simulations. All offset simulations show significant departures from the baseline, demonstrating that the secondary lobe is becoming both longer and wider than the primary with time.

Following Hardcastle & Krause (2013) we define the lobe axial ratio as the ratio between the lobe length and the lobe width, as measured at the midpoint of the lobe; the axial ratio for all simulations is plotted in Figure 7. The lobes produced by relativistic jets are initially elongated with high axial ratios, during the breakout phase of the jet. As the jet head propagation slows, the lobes begin to inflate because they are overpressured compared to the environment, causing the axial ratio to approach the range $[1.5,2]$. This is significantly lower than in the non-relativistic simulations of Hardcastle & Krause (2013). Since a large fraction of observed axial ratios is between $1.5$ and 2 (Mullin et al., 2008, Figs 6-9), our relativistic jets with self-consistent hydrodynamic collimation may explain many radio sources better than non-relativistic ones, although projection effects non-trivially affect this.

In contrast to the length and width ratio evolution, the difference in axial ratio between the primary and secondary lobes for the $1r_{\textrm{c}}$-offset relativistic simulations is small. At later times this difference becomes more pronounced, with the primary lobes having a consistently higher axial ratio compared to the secondary lobes. The width ratio evolution is a function of both the length ratio (as the lobe width is measured at the midpoint between the core and the hotspot), as well as the pinching occurring in the primary lobes.

Relativistic jets produce initially narrow lobes. The length evolution of the lobe is driven by balancing both the jet head pressure and momentum flux against the ambient density and pressure, while the transverse expansion is determined solely by the lobe pressure (Hardcastle & Krause, 2013). All simulations have axial ratios that evolve with time, indicating that they are not expanding self-similarly, consistent with analytical models of Turner & Shabala (2015); Hardcastle (2018).

In addition to the purely morphological differences, the primary and secondary jets have different lossless lobe luminosities and size-luminosity tracks, as shown in Figure 8 and Figure 9. The secondary jet has consistently higher lobe luminosities for a given time, due to the increased volume. Lobe pressure can be considered homogeneous across the entire radio source due to the high sound speed, except for the overpressured hotspots. This results in a strong correlation between luminosity and lobe volume for the lossless calculations used in this work, and so the secondary jets produce lobes with higher luminosities due to their larger volume. The luminosity evolution of all our primary lobes is very similar; the same is true for the luminosity evolution of all our secondaries. Hence, we conclude that our results are robust for variations of the angle the jet makes with the environment symmetry axis.

Figure 9 shows the impact of environment profiles on size-luminosity tracks. The primary PD tracks of the offset simulations begin to flatten out after the initial $\sim 80\,\textrm{kpc}$, as the environment itself starts to flatten out. The secondary PD tracks on the other hand continue rising for the $1r_{\textrm{c}}$-offset simulations inclined at $0$ and $15\degr$.

We show synthetic radio images in Figure 10. We note once more that radiative losses would make the central parts of the lobes less prominent and increase the surface brightness contrast between the hotspot and lobe. Despite the jets becoming unstable and undergoing entrainment, the general source structure is FR II-like and clear hotspots are produced. Subtle differences are consistently visible between the primary and secondary lobes. The primary jets exhibit enhanced surface brightness at the hotspots as expected from theory (Kaiser et al., 1997).