MHD SIMULATIONS OF ACTIVE GALACTIC NUCLEUS JETS IN A DYNAMIC GALAXY CLUSTER MEDIUM

P. J. Mendygral; T. W. Jones; K. Dolag

doi:10.1088/0004-637X/750/2/166

1. INTRODUCTION

Observations of powerful outflows from active galactic nuclei (AGNs) show that they have a considerable effect on the energy budget and morphology of their host galaxy cluster. The jets from the supermassive black hole that powers an AGN inflate low density bubbles in the host galaxy cluster's intracluster medium (ICM) that are observed as X-ray cavities. The minimum energy required to inflate these cavities can be as high as 10⁶¹ erg (e.g., McNamara et al. 2005; Wise et al. 2007) and typically imply a power on the order of 10⁴¹–10⁴⁴ erg s⁻¹ averaged over the age of the cavity (e.g., Bîrzan et al. 2004, 2008; McNamara & Nulsen 2007; Diehl et al. 2008). Most X-ray cavities are paired as a result of the twin AGN jets that formed them, and most are filled with radio synchrotron emission from relativistic, cosmic-ray (CR) electrons; that is, "radio lobes." Absent thermal conduction models of the ICM without any heating source predict that galaxy clusters should establish "cooling-flows," where radiatively cooled ICM plasma falls toward the central galaxy driving up star formation to very high values (Fabian 1994). Observations do not show evidence of these classical cooling-flows, but in contrast, they show that mass accretion onto the central galaxy is at least an order of magnitude below cooling-flow predictions (Peterson & Fabian 2006). One popular explanation for regulating cooling in clusters is that energy output from the AGN in the central galaxy heats the ICM sufficiently to limit accretion. Indeed, the enthalpy content estimated in X-ray cavities appears to be sufficient to offset cooling of the ICM (e.g., Bîrzan et al. 2004, 2008). While the details of the physical processes involved in the transfer of energy from AGN outflows to the ICM are not fully understood, it is generally agreed that these outflows contribute to the dynamics and thermodynamics of the ICM.

At the same time there is also good evidence that the ICMs significantly affect AGN outflows. For example, several radio lobes that have been observed, such as Hydra A (Wise et al. 2007) and 3C 31 (Laing et al. 2008), show complex morphologies with various kinks and turns. The cited authors suggested that, for these objects, the jets and lobes may be reacting to bulk motions of the ICM.⁶ Further, there are groups of AGN-produced radio sources specifically classified by their bent morphologies; namely, narrow-angle tails (NATs) and wide-angle tails (WATs), depending on the apparent degree of bending. Owen & Rudnick (1976) proposed that these objects obtain their bent shape from the relative motion through the ICM of the galaxy hosting the AGN engine. In the case of WATs, the low peculiar velocity of the associated galaxies in combination with reasonable assumptions for the density and velocity of the plasma within the AGN outflows suggest that the motion of the ICM within the cluster is more responsible for their characteristic C-shape than the motion of the individual host galaxy within the cluster (e.g., Burns et al. 1994; Douglass et al. 2008). Burns et al. (1994) argued that cluster mergers can produce sufficient ram pressure in the ICM to cause this shape, and simulations have shown that motions produced by merger events with sufficient ram pressure to bend WATs are persistent for timescales exceeding 10⁸ yr (Loken et al. 1995; Roettiger et al. 1996; Ricker & Sarazin 2001). Similarly, recent simulations show large-scale ICM "sloshing" several gigayears after a gravitational encounter with a passing subcluster (e.g., Zuhone et al. 2010).

Direct detection of such bulk ICM flows became feasible with the Advanced Satellite for Cosmology and Astrophysics, and the Chandra X-ray Observatory has further improved such measurements. Dupke & Bregman (2006), for example, used Chandra observations to explore the velocity structures in the Centaurus Cluster (A3526), and they found that significant departures from hydrostatic equilibrium are possible even in a cool cluster that shows little evidence of a strong merger. ICM sloshing motions have also been observed in the central regions of several cool-core clusters (see Markevitch & Vikhlinin 2007 for a review).

There is limited simulation literature addressing interactions between AGN outflows and dynamical ICM flows, or "weather." For example, Balsara & Norman (1992) and Porter et al. (2009) carried out three-dimensional simulations of AGN outflows into supersonic cross winds that produced NAT-like structures. Loken et al. (1995) used three-dimensional hydrodynamic simulations of a low density jet obliquely crossing a plane shock to demonstrate that merger-produced shocks can decollimate a jet and cause it to bend. In a study of AGN energy deposition to the ICM, Heinz at al. (2006) and Morsony et al. (2010) used output from a cosmological simulation as initial conditions for hydrodynamical simulations that included a pair of AGN jets. To differentiate behaviors fully from hydrostatic clusters the latter study specifically selected a cluster that was far from being relaxed and that included a net rotation and a dynamically induced cold front. Not surprisingly, those authors found significant effects on jet propagation and lobe morphology due to the rather severe ICM weather.

Especially in the context of clusters exhibiting cavities produced by central dominant galaxies, which more typically associate with apparently more relaxed clusters, however, it is important to understand if such effects can take place in those types of clusters. Except for the Porter et al. (2009) idealized NAT study, there are no published simulations of AGN interactions with dynamical ICMs that included magnetic fields. However, since magnetic fields both introduce potential dynamical influences and are also essential to modeling the generation and transport of diagnostic nonthermal emissions, it is important to extend such simulations into the MHD regime. O'Neill & Jones (2010; hereafter OJ10) recently conducted a series of three-dimensional MHD simulations of steady and intermittent jets in an analytically defined, magnetized ICM close to hydrostatic equilibrium. The magnetic field was disordered on scales similar to the size of the ICM core and on average produced a pressure roughly 1% of the gas pressure, so probably characteristic of cluster strength fields. They found that for localized regions, those magnetic fields can have a significant effect on the morphology of the jet lobes. Whether this remains true for a dynamically self-consistent cluster evolved through a cosmological simulation has not been examined. Another important result from the OJ10 simulations was the lack of any significant deflections of the jets or lobes. The steady and intermittent jets in their study both maintained relatively straight-line paths, suggesting that a dynamical ICM is required to produce bent jet and/or bent lobe morphologies.

To address some of these open issues we present here a pair of three-dimensional simulations of MHD jets in a dynamic, magnetized ICM extracted from a high-resolution cosmological simulation. The selected cluster was chosen from among the more relaxed systems in the simulation. To expand the insights into these interactions and the influence of ICM weather on radio lobe morphology the two simulations directed the AGN jet pairs in two mutually orthogonal directions that probed distinct regions of the ICM. The paper is organized as follows. Section 2 describes the numerical setup for the calculations, Section 3 describes the initial conditions in the cluster, Section 4 qualitatively describes the results of the simulations, Section 5 discusses the effects of weather on the jets and lobes, and Section 6 summarizes the major conclusions. We assume below that H₀ = 70 km s⁻¹ Mpc⁻¹, Ω_M = 0.3, and Ω_Λ = 0.7.

2. NUMERICAL DETAILS

The simulations of AGN jets for this study were performed on a 1008×1008×1008 three-dimensional Cartesian grid with a uniform resolution of Δx = Δy = Δz = 1 kpc per zone using WOMBAT (P. J. Mendygral et al. 2012b, in preparation), which is a versatile, efficient, and highly scalable MHD code. Here, WOMBAT evolved the equations for non-relativistic ideal MHD with adiabatic index, γ = 5/3, using the second-order, total variation diminishing (MHDTVD) method of Ryu & Jones (1995). The solenoidal constraint of the magnetic field was maintained to machine accuracy with the "constrained transport" (CT) method of Ryu et al. (1998), which solves the induction equation with a directionally unsplit update using electromagnetic fields (EMFs) derived from the Riemann solver. The simulations included the total energy conserving method for gravity available in WOMBAT and described in detail in P. J. Mendygral et al. (2012a, in preparation). This method adds gravitational source terms to MHDTVD in a manner that maintains second-order code accuracy. The evolution of a passive CR electron momentum distribution, $f({\bm x}, p, t)$ , was also calculated in the simulations with a restricted form of the "Coarse Grained Momentum Volume" (CGMV; Jones & Kang 2005) implementation in WOMBAT. CGMV solves the diffusion–convection equation for a population of relativistic particles. For these simulations spatial CR diffusion outside shocks was neglected, while adiabatic changes in momentum, test-particle diffusive acceleration at shocks and radiative losses by synchrotron emission and inverse Compton scattering of cosmic microwave background (CMB) photons were included (see further details below and, e.g., Tregillis et al. 2001).

The boundary conditions for these simulations were based on a technique used by OJ10. To limit any influence from the boundaries, we used a very large domain (over a Mpc along a side) and modified continuous boundaries that limit the fluxes through those surfaces. In the boundary zones, a hydrostatic atmosphere with a constant sound speed was derived by extrapolating conditions from the last physical zone. Since the mass distribution was known on a much larger domain (see below), the gravitational potential was known out to radii far beyond the computational domain. Consequently, the exact value for the potential described below was used through the boundaries.

2.1. Setting Up the Cluster

The initial conditions for the simulations described here were extracted from a high-resolution simulation of a galaxy cluster performed with an MHD implementation of the smoothed particle hydrodynamic (SPH) code GADGET-3 (Dolag & Stasyszyn 2009). The selected cluster is g676 from Dolag et al. (2009), re-simulated at extreme high-resolution with MHD enabled (for details see F. Stasyszyn et al., in preparation). In short, the cluster used for this study has been extracted from a re-simulation of a Lagrangian region selected from a cosmological, lower resolution dark matter (DM) only simulation (Yoshida et al. 2001). This parent simulation has a box size of 684 Mpc and assumed a baryon fraction of f_bar = 0.13 and σ₈ = 0.9 for the normalization of the power spectrum. Using the "zoomed initial conditions" technique (Tormen et al. 1997), this region was re-simulated with higher mass and force resolution by populating the Lagrangian volume with a larger number of particles, while appropriately adding additional high-frequency modes drawn from the same power spectrum. The initial unperturbed particle distribution (before imprinting the Zeldovich displacements) was realized through a relaxed glass-like configuration (White 1996). Gas was then added to the high-resolution regions by splitting each parent particle into a gas and a DM particle. The gas and the DM particles were displaced by half the original mean inter-particle distance, such that the center of mass and the momentum of the original particle are conserved. The final mass resolution of the DM and gas particles in our simulations is m_dm = 12.1 × 10⁶ M_☉ and m_gas = 2.2 × 10⁶ M_☉, respectively. Thus, the clusters within its virial radius are resolved with 11 × 10⁶ DM and a similar number of gas particles. In the simulations, the gravitational softening length was chosen to be epsilon = 1.4 kpc which corresponds roughly to the mean particle separation in the center of the cluster at z = 0. For the initial magnetic field we choose a space filling, homogeneous, primordial magnetic field of B = 10⁻¹¹ G comoving as in Dolag et al. (2009).

For this work, we were interested in assessing the effects of cluster "weather" specifically for a cluster that observationally appears relaxed yet harbors significant bulk ICM flows. We selected the state of g676 at a redshift of z ∼ 0, which was the least dynamic period for this cluster. As seen in Figure 1, the last major merger (mass ratio of 0.27) occurred at z ≈ 0.8 (a little less than 7 Gyr before our jet simulations start), and the DM mass continued to grow slowly after that. The general morphology and environment of the cluster at z = 0 can be seen in Figure 1 of Dolag et al. (2009). Volume-averaged values for baryonic density, temperature, velocity, and the magnetic field were all computed from the SPH particles and mapped to a uniform grid matching the jet simulation grid. The simulation assumed a hydrogen fraction of X = 0.76 and a helium fraction of Y = 0.24, resulting in a mean molecular weight of μ = 0.59. An ideal equation of state was used to convert between pressure and temperature. Radiative cooling of the thermal plasma was not included in the GADGET calculations that formed cluster g676, so neither was it used in the jet simulations. A summary of relevant cluster ICM properties is given in Section 3.

**Figure 1.** Evolution of velocity dispersion (dark matter) and total mass as a function of redshift for the selected galaxy cluster *g676*. The vertical lines mark identified merger events with a mass ratio indicated by the labels.
Download figure:
Standard image High-resolution image

For the jet simulations, we used a static gravitational potential defined from the total gravitating mass (DM and baryonic matter) azimuthally averaged from the SPH values at the time used for the initial conditions. Figure 2 shows the gravitating mass as a function of radius. The total enclosed mass out to the virial radius of 1.4 Mpc is 1.53 × 10¹⁴ M_☉. The decision to use a static potential as opposed to evolving a population of DM particles and maintaining an updated full potential was based on simplicity, the relaxed nature of the cluster, and the relatively brief duration of the jet activity (200 Myr) compared to the dynamical timescale for the cluster mass distribution. Obviously over longer times or in a cluster involved in significant active restructuring it would be important to follow the evolution of the gravitational potential.

**Figure 2.** Gravitating mass profile used at all times to define the gravitational potential. The mass includes both dark matter and baryonic matter from the GADGET simulation and was derived from an azimuthal average at each radius.
Download figure:
Standard image High-resolution image

2.2. Setting Up the Jets

For each simulation, a pair of oppositely directed jets was launched from a cylinder centered in the grid (X = Y = Z = 0). The launching cylinder had a radius of r_j = 3 kpc and length on each side of the jet origin of l_j = 8 kpc. An additional 3 kpc thick collar and cap surrounded the cylinder to transition from jet conditions to ICM conditions and to contain the return electric current associated with the jet magnetic field (Section 2.2.2). Special consideration was given to how the cylinder interacted with the ICM. In these simulations, ambient flows could in principle enter the cylinder, even from the initial conditions. This can be particularly problematic if ICM magnetic fields enter the cylinder, where they can blend with the jet field and produce undesirable field geometries.

Two steps were taken to alleviate this. First, at the start of the simulation the ICM material in the vicinity of the jet cylinder location was gently pushed aside before the jet cylinder was initialized. In addition, reflecting boundaries were defined as defaults for the cylinder. The initial ICM evacuation was accomplished through a purely solenoidal velocity field imposed at the start of the simulation on a sphere with radius r_j centered on the grid with v_ϕ = v_e(r/r_j), and v_e set to the local sound speed. This sphere was overpressured by 90% as well. These conditions were allowed to evolve for 3.28 Myr, which was the sound crossing time over r_j. This removed ≈25% of the ICM plasma and magnetic energy from the jet cylinder region. The jet cylinder was then established at t = 3.28 Myr. To divorce the ICM and jet magnetic fields, any remaining ICM field penetrating the jet cylinder or collar was set to zero at this time. This last step effectively introduced some magnetic monopoles around the cylinder, but they had negligible effect on the jet flows and were maintained to low levels by the CT scheme in WOMBAT. The collar region around the active jet cylinder acted as an impenetrable electrically conducting surface because of the imposed reflecting boundaries for both the normal momentum and magnetic field. This prevented the loss of any conserved quantity into the launching region. The caps of the cylinder transitioned between reflecting boundaries and outflowing jet properties as the jet cylinder transitioned between active and inactive states.

An important feature of the jet cylinder for this investigation was that it could be oriented in any direction on the grid. This allowed jets to be launched at an arbitrary angle with respect to the large-scale flows and pressure structures in the ICM. To achieve this, we defined all vector quantities in the jet cylinder according to cylindrical coordinates. The jet cylinder is azimuthally symmetric, so only two coordinates are needed, z and r, where z is along the cylinder symmetry axis and r is radially outward. By default the jet z-axis was aligned with the grid Z-axis, and Euler angles were used to reorient the cylinder. After the rotation, all cylinder coordinates were converted into grid coordinates so that quantities could be mapped to the Cartesian grid. The specific orientations of the two jets in this study are discussed in section 4.

The active jet density and pressure were set to ρ_j = 4.165 × 10⁻²⁸ g cm⁻³ and P_j = 2.499 × 10⁻¹⁰ dyne cm⁻² (approximately 1/200 the initial local ICM density) and in approximate pressure equilibrium with the initial surrounding ICM. For the assumed adiabatic index of γ = 5/3, this gave a sound speed inside the jet of c_{s, j} = 10⁴ km s⁻¹ or 0.03c. A passively advected "color" (or mass fraction) tracer, C_j, was also introduced in the jet cylinder that was used to identify AGN and ICM plasma. The color was set to 1 in the jet cylinder, while all ICM material was initially assigned C_j = 0.

2.2.1. Jet Velocities

When it was active the jet velocity in the jet cylinder was given a spatial ramp from its origin in z to allow a smooth transition between the oppositely directed flows of the two jets. In particular, the velocity had the form, $v_{z,j}(z) =v_j \mathcal { F}(z)$ , where

$\begin{eqnarray} &&\mathcal {F}(z) = {\rm TANH}\left(1.8 \left[\frac{z}{l_{j}}\right]^{2}\right). \end{eqnarray} \tag{ 1 }$

The value 1.8 was determined empirically by minimizing the automatic application of diffusive, positive-definite Harten, Lax, and van Leer (HLL) fluxes by WOMBAT in the jet cylinder, which are used to eliminate unphysical pressures during peak flow times. The velocity in the jet cylinder was defined with this ramp to be v_{z, j}(z, t) = v_j(t)F(z), where v_j(t) = 1.2 × 10⁴ km s⁻¹, when the jet was fully on (see Section 2.2.3). This yields a pair of supersonic emergent jets with the flow at Mach number, M_j = 1.2, with respect to the jet sound speed. The peak speed of v_j corresponds to a Lorentz factor of 1.0008 and is sufficiently low that a non-relativistic MHD solver is applicable. The jet velocity was uniform in r across the jet core (r ⩽ 3 kpc), with a transition to the ambient medium velocity ∝r^−1/3.

The core luminosity of each of the two jets, ignoring small contributions from magnetic energy, is given by

$\begin{eqnarray} L_{j} &=& \pi r_{j}^{2} v_{j}\left(\frac{1}{2}\rho _{j}v_{j}^{2} + \frac{\gamma }{\gamma - 1}P_{j}\right)\nonumber\\ &=& \pi r_{j}^{2}\frac{1}{2}\rho v_{j}^3\left(1 + \frac{2}{M_j^2(\gamma - 1)}\right), \end{eqnarray} \tag{ 2 }$

where γ = 5/3. When the jets are fully active, the preceding parameters give a total luminosity (including both jets) of L_j = 6 × 10⁴⁴ erg s⁻¹. The transitional collar region adds roughly another 50% to the actual jet luminosity. The ratio of kinetic energy flux to enthalpy flux for these jets was ρ_jv²_j/5P_j = 0.48. The actual total energy injected onto the grid by the jet cylinder in these runs was 2.2 × 10⁶⁰ erg.

For later discussion it is also useful to know the thrust of the jets, F_j. The thrust is easily computed from the luminosity (with γ = 5/3) as

$\begin{equation} F_j =L_j \frac{2}{v_j \big(1 + 3/M_j^2\big)}, \end{equation} \tag{ 3 }$

so that the individual core jet thrusts at peak luminosity are F_j ≈ 1.6 × 10³⁵ dyne.

2.2.2. Jet Magnetic Field

The magnetic field in the jet launching cylinder was derived with a method that maintained $\nabla \cdot \mathbf {B}$ to machine accuracy. The jet field was purely toroidal inside the jet cylinder with B_ϕ = B_j(r/r_j) when r < r_j, where B_j = 7.92 μG (plasma β = 8πP_j/B²_j = 100). Outside of that radius, the field dropped off as r⁻³ and was zero outside of the jet cylinder collar. This kept the jet field separated from the ICM field and made the net electric current in each jet zero. The arbitrary orientation of the jet cylinder was an additional complication for defining the jet field. We implemented this by defining an electric field that determined the rate of change of the magnetic field in the cylinder according to Faraday's Law, $-\partial \mathbf {B}/\partial t = \nabla \times \mathbf {E}$ . Solving this equation allowed the magnetic field to be replaced as it was advected out of a zone at an arbitrary angle. Given the desired form of the magnetic field and the jet velocity as a function z, this electric field was

$\begin{eqnarray} &&E_{z}(r,z) = \left\lbrace \begin{array}{@{}l}B_{j}v_{j}(z)\frac{r_{j}}{2\Delta z}\left\lbrace \left(\frac{r}{r_{j}}\right)^{2} - \left[\frac{5}{4} - \frac{1}{4}\left(\frac{r}{r_{j}}\right)^{4}\right]\right\rbrace \\ \ms\ms \quad\quad {\rm for }\; r < r_{j}\; {\rm and }\; |z| \le l_{j} \\ \ms\ms B_{j}v_{j}(z)\frac{r_{j}}{2\Delta z}\left\lbrace \left[\frac{5}{4} - \frac{1}{4}\left(\frac{r_{j}}{r}\right)^{4}\right] - \left[\frac{5}{4} - \frac{1}{4}\left(\frac{r_{j}}{r_{b}}\right)^{4}\right]\right\rbrace \\ \ms\ms \quad\quad {\rm for }\; r_{j} < r < r_{b}\; {\rm and }\; |z| \le l_{j} \\ \ms\ms 0 \\ \ms \quad\quad {\rm otherwise}, \end{array}\right. \nonumber\\ \end{eqnarray} \tag{ 4 }$

where r_b = r_j + 3Δz. The magnetic field added to the zones in the jet cylinder and collar was simply $\delta \mathbf {B} = \Delta t(\nabla \times \mathbf {E})$ , where Δt is the current time step. The value of $\nabla \times \mathbf {E}$ only needed to be computed once at the beginning of the simulation.

2.2.3. Jet Source Time Evolution

AGN activity is likely to be unsteady on the ⩾10⁸ yr timescales required to form large-scale cavities, while the intermittency of jet activity has significant impact on the interactions with the ICM (see, e.g., OJ10). In these simulations, the activity of the jets was toggled on and off over approximately 26 Myr as a simple way to mimic such cyclic activity. In particular we applied to the characteristic jet velocity, v_j, and magnetic field, B_j, a ramp function, R(t), defined as

$\begin{eqnarray} R(t) &=& \eta (t) {\rm MAX}\left[0,{\rm TANH}\left(\frac{2(t - t_{{\rm on}})}{t_{{\rm blend}}}\right)\right]\nonumber\\ && + \left(1 - \eta (t)\right){\rm MAX}\left[0,-{\rm TANH}\left(\frac{2(t - t_{{\rm off}})}{t_{{\rm blend}}}\right)\right],\quad\quad \end{eqnarray} \tag{ 5 }$

$\begin{equation} \eta (t) =\left\lbrace \begin{array}{@{}ll}0 & {\rm for }\; t \ge \frac{t_{{\rm on}} + t_{{\rm off}}}{2} \\ 1 & {\rm otherwise}, \end{array}\right. \end{equation} \tag{ 6 }$

where t_on and t_off are the start and end time for the current jet cycle and t_blend = 3.28 Myr. A duty cycle of ∼50% was used for both simulations with a period of 26.2 Myr. Thus, the jet velocity and added magnetic field were given a time dependence as v_j(z, t) = v_j(z)R(t) and $\delta \mathbf {B}(r,z,t) = \Delta t[\nabla \times \mathbf {E}(r,z)]R(t)$ .

2.3. Cosmic-ray Electrons

In order to enable realistic synthetic observations of nonthermal emissions we incorporated in the jet outflows a population of test-particle CR electrons with Lorentz factors from 10 to 1.6 × 10⁵ and evolved their distribution, $f({\bm r}, p, t)$ , using the CGMV routine in WOMBAT. This range in Lorentz factors is sufficiently large for the typical magnetic field strengths in these simulations to model synchrotron emission from a few MHz to tens of GHz. No CR electrons were included in the initial ICM, although any electrons introduced by the jets and then mixed into the ICM were followed in the CGMV scheme. For these simulations, CRs were subject to adiabatic changes in momentum as well as energy losses due to synchrotron emission and inverse Compton scattering of CMB photons. CRs were also subject to first-order Fermi acceleration at shocks (diffusive shock acceleration (DSA)).

As described by Tregillis et al. (2001), the characteristic scattering length for electrons in this energy range and in field strengths of a few μG, typical for these simulations, is on the scale of the solar system, and the acceleration timescale is on the order of years. These scales are both unresolved in these simulations, which resolves lengths of a kiloparsec and time in kiloyears. This mismatch in scales can be exploited by assuming that the acceleration process for CRs at a shock is completed in a time that is far less then a typical dynamical time step in the simulation. For this application, the CGMV implementation does not include diffusion terms and simply redistributes existing post-shock CRs into a power law according to test-particle DSA theory. The slope q_s of the resulting power law is related to the measured compression ratio of the shock, r, by q_s = 3r/(r − 1). To implement DSA WOMBAT includes a directionally unsplit shock detector that locates and characterizes shock strengths. Eight boundary zones along domain boundaries allowed accurate shock detection and characterization up to the edge of the local domain boundary and the physical grid boundary.

The CR population was introduced in the flow of AGN plasma out the jet cylinder. Since the CR electrons were passive, their number was somewhat arbitrary, except in the calculation of nonthermal emissions. For the synthetic observations reported here, we assume the emergent jets carried a population of CRs with a density of 5 × 10⁻⁷ cm⁻³ (about 0.2% of the jet bulk proton number density) and had a single power-law slope of q_j = 4.5. For consistency this population should not be expected to contribute significantly to the subsequent flow dynamics. In fact, with this density and slope, the nominal CR electron pressure at the jet source was 4 × 10⁻¹² dyne cm⁻², or approximately 2% of the "thermal" gas pressure in the AGN plasma. Especially, given the softer equation of state for these CRs, the downstream CR pressure influences would, indeed, have been minor. Since the nonthermal emissions of interest to us here are associated with the jets, we disabled injection of fresh CR electrons from the thermal pool for these simulations in order to avoid contamination by CRs that would have been injected in ICM shocks.

2.4. Synthetic Observations

Synthetic observations of the simulations were calculated following techniques similar to those described in Mendygral et al. (2011). We focus on two models for optically thin emission; thermal bremsstrahlung from the hot plasma and synchrotron radiation from the CR electrons. X-ray emission in an energy band typified by Chandra and low frequency radio emission similar to those obtainable by LOFAR are of particular interest. We note that, although emission due to inverse Compton scattering of CMB photons is simple to include, it is negligible for the energy/frequency range of interest with the CR electron numbers we assume here.

The thermal bremsstrahlung emissivity was calculated in every computational zone as

$\begin{eqnarray} j_{\nu _{{\rm local}}} &=& 5.4\times 10^{-39}\,g_{ff}(\nu _{{\rm local}},T_{e})\nonumber\\ && \times Z_{i}^{2}\frac{n_{e}n_{i}}{T_{e}^{1/2}}e^{-h\nu _{{\rm local}}/T_{e}}\, {\rm erg\; cm}^{-3}{\rm s}^{-1}\,{\rm sr}^{-1}, \end{eqnarray} \tag{ 7 }$

where ν_local = ν_obs(1 + z) for a redshift of z, T_e is in keV, and all other quantities are in cgs units. The free–free Gaunt factor, g_ff, was computed by interpolation from the values calculated for plasma with typical ICM properties in Table 1 of Nozawa et al. (1998). For simplicity, we assumed a fully ionized, Z_i = 1 pure hydrogen gas with an ideal equation of state. The pure hydrogen assumption, however, does not match the abundances that were used to define the composition of the baryons described in Section 2.1. Although formally the Gaunt factor does depend on the ion charge Z_i, the modifications to g_ff for typical temperatures and observations energies of these simulations for a Z_i = 2 gas were only a few percent from Table 1 of Nozawa et al. (1998). The remaining dependence is the Z²_i term in Equation (7). Since the abundance is uniform on the grid, this only sets a normalization for bremsstrahlung emission. Finally, only for the thermal emission computation, the mean molecular weight μ = 0.5 was used to define the temperature in a zone as T_e = T_i = T(keV) = μPm_H/(1.602 × 10⁻⁹ρ), where P, m_H, and ρ are all in cgs units. To simulate the bandwidth of real instruments, Equation (7) was numerically integrated over a range of frequencies. Unlike the technique used in Mendygral et al. (2011), where emission from zones with AGN plasma was clipped for numerical reasons, we do not eliminate emission from those zones in this study. Significant and complex mixing between ICM and AGN plasma in these simulations prevented a clear cutoff based on the passive color tracer, C_j.

Synchrotron emissivity was calculated using the local CR spectrum and was given by Jones et al. (1974) as

$\begin{eqnarray} &&\epsilon _{\nu _{{\rm local}}} = j_{\alpha }\frac{4\pi e^{2}}{c}f_{p_{c}}p_{c}^{q_{c}}\left(\frac{\nu _{B}^{\alpha + 1}}{\nu _{{\rm local}}^{\alpha }}\right), \end{eqnarray} \tag{ 8 }$

where $f_{p_{c}}$ is the CR phase space density at momentum p_c, q_c is the CR spectral slope at p_c, and j_α is a constant of order unity. The spectral index is related to the CR spectral slope as α = (q_c − 3)/2. The critical momentum is calculated in units of m_ec as $p_{c} = \sqrt{2 \nu _{{\rm local}} / (3 \nu _{B})}$ , where the gyrofrequency is given by ν_B = eB/(2πm_ec) and B is the magnitude of the magnetic field in the plane of sky.

A ray casting engine was used to convert volume emissivities into images with the source set at a user-defined distance. In the engine, rays were cast normal to an image plane that was rotated to a user-defined orientation with respect to the grid. Trilinear interpolation was used to determine the emissivity at regular intervals along the ray. The sum of emissivities along a line of sight determined the intensity, and the intensity was then converted into a flux by applying the appropriate redshift correction and multiplying it by the solid angle of a pixel determined by the luminosity distance D_L of the source. For both X-ray and radio observations presented here, the computational domain was set at a redshift of z = 0.0594 (D_L = 240 Mpc), which is approximately the same as the Hydra Cluster (Wise et al. 2007).

3. CLUSTER PROPERTIES: INITIAL CONDITIONS

To set the stage for analysis, we outline now and illustrate in Figures 3–5 the initial properties of the ICM for our cluster before the jets were launched. Figure 3 shows from the simulation initial conditions azimuthally averaged radial profiles for gas density and pressure (along with the derived temperature and sound speed) as well as the ICM flow speed with respect to the DM core and the magnetic field strength. The variance and extremes of each quantity are also shown. The cluster has a dense core with a central ICM density of ρ_{ICM, 0} = 1.05 × 10⁻²⁵ g cm⁻³ and pressure P_{ICM, 0} = 2.69 × 10⁻¹⁰ dyne cm⁻². This gives a core temperature of T_{ICM, 0} = 1.6 keV, for the assumed mean molecular weight of μ = 0.6, and a sound speed of c_{ICM, 0} = 653 km s⁻¹. The azimuthally averaged temperature falls off slowly with radius, which is typical for non-cool-core clusters (Hudson et al. 2010). This is not surprising, as radiative cooling was not included in the GADGET simulation that evolved the cluster to this time. The average magnetic field magnitude in the core is B_{ICM, 0} = 4.2 μG, which corresponds to a relatively high plasma β_{ICM, 0} = 8πP_{ICM, 0}/B²_{ICM, 0} = 382. Figure 4 shows a radial profile of the ICM entropy for the initial conditions. This cluster did have a low entropy core, which is a standard property of SPH simulations like the one that produced the cluster (e.g., Mitchell et al. 2009).

**Figure 3.** Radial profile of ICM properties from the initial conditions. The mean with standard deviations are the solid lines surrounded by the colored band, and the minimum and maximum values at each radial sample are the dashed lines. Note that the lower variation values for the magnetic field magnitude were consistent with zero for some radii and were then truncated.
Download figure:
Standard image High-resolution image

**Figure 4.** Radial profile of ICM entropy from the initial conditions. The mean with standard deviations are the solid line and color band, and the minimum and maximum values at each radial sample are the dashed lines.
Download figure:
Standard image High-resolution image

Although this cluster was selected especially because its last major merger was several gigayears in the past and at first appearance it seems relatively relaxed, closer inspection reveals significant ICM "weather" characteristic of sloshing behaviors, thus, reminding us that g676 has continued to interact with smaller DM halos and to accumulate mass slowly. Note first from Figure 3 that the mean magnitude of the ICM velocity in the inner core is v_{ICM, 0} = 142 km s⁻¹ (Mach 0.2) and that it increases to ≈400 km s⁻¹ by 100 kpc from the cluster center. For a typical sound speed of 450 km s⁻¹ outside 100 kpc, these motions approach Mach 1. They are primarily relatively large-scale flows across the cluster, since by contrast Zhuravleva et al. (2011) report isotropic, random velocities with respect to the mean of only ∼50 km s⁻¹ in the core of this cluster, increasing to around ∼100 km s⁻¹ outside a 100 kpc radius. The turbulent pressure in the core region (P_turb ∼ 3 × 10⁻¹² dyne cm⁻²), therefore, is less than 1% of the core thermal pressure (Figure 3 shows profiles of density and thermal pressure).

Further indications of the sloshing behavior are found in Figure 5, which shows orthogonal slices of gas entropy, S = P/ρ^γ. The initial cluster conditions were rotated such that the jet axis from R1 (see Section 4) is in two of the planes on the left panel, and the panel on the right was created in the same fashion for R2 (see Section 4). An inverse color map was used, where lighter colors correspond to lower entropy. The jet axis is shown as a dark line from the cluster center. Overlaid are ICM velocity vectors whose length corresponds to the gas speed. Both views show complex flow patterns and highly asymmetric entropy distributions, which are classic signatures of sloshing (Zuhone et al. 2010). Despite these flow patterns this cluster does still qualify by common observational measures as an example of a relaxed system. A synthetic X-ray observation of the initial cluster in the 0.5–8 keV band is shown in Figure 6 along a line of sight roughly orthogonal to the mean ICM velocity in the cluster core. This line of sight would also be along the jet axis from the R1 simulation (see Section 4). Although there is some asymmetry of the cluster core for this particular line of sight, other observational orientations lead to more symmetric images. As a quantitative measure we construct from this image the concentration parameter, c, sometimes used to separate relaxed clusters from those that have recently undergone mergers (e.g., Santos et al. 2008):

$\begin{equation} c = \frac{S(R < 100 {\rm \,kpc})}{S(R < 500 {\rm \,kpc})}, \end{equation} \tag{ 9 }$

where S is the X-ray flux in the aperture. A relaxed cluster would have a more concentrated core, so a higher value of c, than a cluster that was recently disturbed in a merger. Cassano et al. (2010) find that a value c > 0.2 separates clusters without radio halos (not recently disturbed) from those with radio halos and a recent major dynamical disturbance. From the X-ray observation in Figure 6, the initial cluster conditions used for this study have a concentration parameter of c = 0.66, so the X-ray profile would not be considered dynamically disturbed by this observational method. This measurement was repeated for two other orthogonal lines of sight through the cluster with variation in c of less than 1%. For reference we note that it has been pointed out (Cassano et al. 2010) that the concentration parameter is less sensitive to projection effects than the so-called power ratio based on multipole expansions of the X-ray brightness distribution (e.g., Buote & Tsai 1995; Jeltema et al. 2005; Ventimiglia et al. 2008; Böhringer et al. 2010; Cassano et al. 2010).

**Figure 6.** Synthetic X-ray observation from 0.5 to 8 keV of the initial cluster conditions. The image is scaled logarithmically over ≈3 orders of magnitude. The line of sight for this observation was along the jet axis from the R1 simulation described in Section 4.1.
Download figure:
Standard image High-resolution image

4. JET SIMULATION RESULTS

Two AGN outflow simulations for the parameters defined above in Section 2.2 were run with the jet cylinder at different orientations on the grid. Both AGNs were placed at rest in the cluster center in order to distinguish the influences of the dynamic ICM from host galaxy motions within the cluster. The first simulation, R1, had the jet cylinder orientation chosen arbitrarily, but so that the jet axis was not along a primary grid axis. As it turned out the jets were directed roughly orthogonal to the mean ICM velocity inside a 200–300 kpc cluster radius, which represents the final scale of the cavities formed in the ICM by the AGN outflows. The second simulation, R2, had the jet cylinder in the orthogonal plane to R1 and in a direction more nearly aligned with the mean interior ICM flows. The total elapsed time for each simulation was 200 Myr. The jets cycled on and off with a 50% duty cycle for 6 periods of 26.2 Myr per cycle and then left off after 160.5 Myr of this activity, allowing the AGN inflated bubbles to evolve without more energy injection for about an additional 40 Myr. We note in passing that a steady jet of the same properties over this same time interval would propagate off the grid. A steady jet would create a structure roughly twice as long as those produced in the two simulations we present. Additionally, the X-ray cavity produced by a steady jet would be significantly less round in shape as compared to those produced by the intermittent jets in the simulations shown here.

4.1. Simulation 1 (R1) Evolutionary Summary

Figure 7 shows color, C_j (= jet mass fraction), volume renderings from the R1 simulation at five representative times (t = 34.4, 49.2, 98.4, 147.6, and 196.8 Myr). The views are orthogonal to the jet axis. All the panels project the same view except the lower right, as indicated by the stick diagrams in the two bottom panels. A line segment in the lower right panel indicates a projected 100 kpc length. To aid structural registration the location of the jet cylinder is shown independent of jet activity.

At 34.4 Myr (upper left), the AGN is active for the second time in the simulation, so the jets are collimated, and there are small lobes formed from the first active AGN phase. As mentioned, these jets drove through a strong cross flow in the ICM associated with sloshing motions. Already at this time the lobes are obviously asymmetric and no longer aligned with the jet axis. At 49.2 Myr the AGN is inactive, and it is apparent that the AGN plasma is concentrated at the ends of the lobes. The panels at 98.4 Myr and 147.6 Myr both show times when the jets are on. In both cases, the visible jets clearly are not as well collimated as they were at 34.4 Myr. The ICM pressure near the launching region is approximately a factor of two lower at these times than the value at 34.4 Myr. Recall from Section 2.2 that the jets were launched with a fixed pressure related to the initial ICM, so the jets are now overpressured as they emerge, causing them to expand. Since the internal jet Mach number is low (M_j ≈ 1.2), re-collimation is slow. In addition, bow shocks have formed ahead of the new outflows, as well as jet termination shocks, further enhancing the dispersion of jet plasma into the cavity associated with the previous activity. The asymmetry between the lobes and their misalignment with the jet outflows are quite obvious, especially in the final views at 196.8 My, sr. The lobes extend roughly the same distance from the cluster center (300 kpc), but there are different offsets of the lobes from the jet axis. The most significant differences between the two lobes are the deflection in roughly the $-\hat{z}^{\prime \prime }$ -direction, so that viewed roughly along the $\hat{y}^{\prime \prime }$ -axis they present a "C-shaped" morphology, and from another, orthogonal perspective, a flattened shape of the left lobe compared to the right lobe, as seen in the bottom right panel of Figure 7. Finally, we note that the AGN plasma is very nonuniformally distributed inside the lobes, with some significant mixing of ICM and AGN plasma in several regions, made evident by the gray (turquoise for color figures) regions in all the figure panels. Such mixing should have various observational consequences, including the Faraday rotation behaviors along sight lines through the lobes, as we will discuss in a subsequent publication.

To that point, however, we mention here that magnetic fields, while not dynamically dominant in these simulations, are essential in the production and transfer of radio emissions. In those roles they can help to illuminate dynamical properties of the AGN/ICM interactions. Figure 8 shows volume renderings from the R1 simulation of the color, C_j, along with selected magnetic field lines colored dark to light by magnitude during an active AGN phase at 137.8 Myr. In the left panel, the jet to the right is coming toward the observer and the left jet is directed away. From this perspective, the AGN outflows are deflected away from the observer. Also in the left panel, ICM field lines can be seen stretched and dragged outward with AGN plasma into the leftward lobe formed from previous activity. Those ICM field are actually draped around the AGN plasma. The field lines colored white show where the stretching from the jet outflow has most intensified the field. OJ10 reported a similar result for their intermittent jet simulations. There is a stray magnetic field line seen in the left panel of Figure 8. Despite attempting to isolate the field within the jet, this line is rendered due to the blending of the jet and ICM magnetic fields in the collar region around the jet cylinder. In the right panel, the magnetic field inside this same lobe is shown from a slightly different perspective wrapping around regions consisting primarily of AGN plasma. Field lines are also evident tracing the outline of the lobe.

**Figure 8.** Volume rendering of passive tracer *C_j* and selected magnetic field lines from R1 during the active AGN phase at 137.8 Myr. Field lines are colored by magnitude with brighter field line colors representing stronger fields. Left: a view showing the tightly coiled jet field near the launching region. ICM field lines are seen stretched over the AGN plasma from a previous AGN cycle. Right: field lines in this view within the leftward jet lobe in the other panel are seen following the outline of the lobe.
Download figure:
Standard image High-resolution image

4.2. Simulation 2 (R2) Evolutionary Summary

The evolution of the R2 simulation is shown in Figure 9 again through volume renderings of the color C_j. Many of the same morphological characteristics of R1 visible in Figure 7 can be seen in this rendering of R2 as well. Just as in the R1 case the R2 jets appear well collimated early in the simulation, but again they overexpand as the propagate through the lower pressure cavities. The AGN plasma in the lobes has a qualitatively similar distribution to R1. Again the AGN plasma is primarily deposited near the outer edges of the lobes. The most striking difference between the R1 and R2 morphologies is that while the lobes in the R1 simulation extended about the same distance from the AGN, the lengths of the two lobes in R2 are of very different length. The lobe on the right side of Figure 9 does not extend nearly as far into the cluster as the lobe on the left side in these images. This is evident even in the early 34.4 Myr image. The asymmetry continues to increase throughout the simulation. At 196.8 Myr the difference in extent is approximately 100 kpc or roughly 50%. Recall that these jets are directed in an orthogonal direction to those in the R1 simulation. In fact, the jet axis is roughly aligned parallel and anti-parallel to the prevailing ICM velocity in the inner cluster mentioned previously. The relationship between cavity morphology and cluster weather will be discussed in more detail in Section 5.

The structure of several magnetic field lines during jet activity at 137.8 Myr from R2 are seen in Figure 10. Similar to the structures in Figure 8, the toroidal jet field is shown in the left panel, along with ICM field lines that have been stretched and wrapped over AGN plasma from the last period of activity. Magnetic field lines inside the leftward lobe are visible in the right panel. Those lines also follow the outline of lobe structures, as they did in the R1 simulation. The fields seen here appear to be less tangled than those seen in Figure 8. This same lobe is the one on the left-hand side of the images in Figure 9. Its distortions are not as extreme as the deflected jet and lobe shown from R1 in Figure 8.

4.3. Synthetic X-Ray Observations

Figure 11 shows 0.5–8 keV X-ray observations of R1 at various epochs taken along a line of sight that was perpendicular to the jet axis. The resolution is 1 arcsec pixel⁻¹ at the assumed 240 Mpc distance, and each observation was divided by the best-fit double β-profile (see Mendygral et al. 2011) to emphasize the X-ray cavities. In these images, the lower X-ray cavity is associated with the left-hand lobe in Figure 7. The observation of the initial conditions shows significant departures in surface brightness from the double β-profile, particularity at large distances from the cluster center (see also Figures 5 and 6) reminding us, once again, that the ICM is not truly relaxed and the cluster is not isolated. At 65.6 Myr, after two full AGN periods, the jets are active and the inner pair and outer pair of cavities are seen separated by a thin bright rim. These features are similar to the multi-cavity system in A2052 (Blanton et al. 2009). The observation at 131.2 Myr shows the same morphology with the pairs of cavities on larger scales. Five complete AGN cycles have been executed by this time. Ripples are seen emanating from the outer pair of cavities, which are outlined by bright rims that are ∼30% hotter than the surrounding gas. The ripples were seen for the intermittent simulation presented in Mendygral et al. (2011), and they are similar to the ripples observed in the Perseus Cluster (Fabian et al. 2003). Additionally, bright rims hotter than their surroundings have been observed in Centaurus A (Kraft et al. 2003, 2007) and NGC 3801 (Croston et al. 2007), but for other objects, the bright rims are, in fact, cooler than their surroundings (Diehl et al. 2008).

The outer R1 cavities at 131.2 Myr are distorted and no longer aligned with the axis of the inner pair of cavities. These distortions are largely due to cluster "weather," as mentioned above and discussed more fully in Section 5. The approximately Mach 2 bow shock is well separated from the cavities at this time. By the time of the final observation, at 196.8 Myr, the AGN has been inactive for 40 Myr, and there is no longer an apparent pair of inner cavities. The cavity system appears to include only one pair that extend ∼300 kpc from the cluster center and are ∼150 kpc wide. The cavities are still outlined by bright rims, and the ripples can still be seen. The lower cavity shows the most distinct distortions. The bow shock, now Mach ∼1.3, is still visible.

Figure 12 shows 0.5–8 keV observations of the R2 simulation, also divided by best-fit double β-profiles. In these images the line of sight is again perpendicular to the jet axis. Many of the same feature seen in the R1 observations are present, such as the ripples from the intermittent AGN activity, bright rims around the large cavities, and a bow shock with a similar Mach number to the R1 case. A unique and very interesting feature of the R2 observations, however, is that the 65.6 Myr and 131.2 Myr observations show an incomplete pair of inner cavities, whereas R1 showed a complete pair. The rim around the upper-inner cavity is visible in these two images, but the lower-inner cavity is minimally outlined with a faint rim. It is unlikely that this cavity would be seen in a real observation of such an object with finite sensitivity, and only three of the four cavities might be identified. In the last image at 196.8 Myr the inner cavities are absent, and only a single pair of cavities with bright rims, that are less noticeable than those from R1, are visible. The lower cavity is significantly shorter and more circular in appearance than the upper cavity. As mentioned previously and discussed in more detail in Section 5, this asymmetry also results from large-scale ICM flows (weather) present in the cluster initial conditions. As mentioned previously, the R2 jets are more or less aligned with this flow, whereas the R1 jets were roughly orthogonal to the mean interior ICM flow.

4.4. Synthetic Radio Observations

Figures 13 and 14 show at 131.3 Myr and 198.2 Myr the 178 MHz synthetic synchrotron intensity observations of the R1 and R2 simulations, respectively. These observations are of the same line of sight as the respective X-ray observations discussed above, but they are from a smaller volume to show more detail. We note to avoid confusion that in all the images the jet launch cylinder remains artificially illuminated. The left-hand panels in both figures were configured with 1.7 arcsec beams (pixels) at the assumed 240 Mpc distance, and the right-hand panels are the same observations but Gaussian smoothed over 6 pixels to simulate a 10 arcsec beam. The line of sight through the grid in Figure 13 (Figure 14) is the same as the X-ray image in Figure 11 (Figure 12). Such low frequency radio observations are of particular interest because they highlight low energy CRs that provide a reasonable tracer for AGN plasma, since their losses are not large. For a typical field value in the lobes of ∼1.5 μG much of the emission at 178 MHz is from ∼3 GeV electrons, whose radiative lifetimes against inverse Compton and synchrotron losses in this context are ∼ 200 Myr, assuming a small cluster redshift.

**Figure 13.** Synthetic observations of synchrotron emission from R1 at 178 MHz. The units on the color bar are Jy beam⁻¹ with a 1.7 arcsec beam (pixel). Left: full resolution observation. Right: same observation smoothed with a Gaussian kernel over 6 beams (pixels).
Download figure:
Standard image High-resolution image

**Figure 14.** Same as Figure 13 for R2.
Download figure:
Standard image High-resolution image

For both simulations at 131.3 Myr, the bright active jets are seen in the figures penetrating their associated lobes. Emission from previous AGN activity extends further out, but is considerably dimmed because of radiative aging of the CR population. These images closely resemble a "double–double" radio galaxy (DDRG; Schoenmakers et al. 2000) such as 3C 293, which also shows a misalignment between the inner and outer doubles (Joshi et al. 2011). At the end of the simulations, most of the DDRG resemblance is gone because roughly 50 Myr has passed since the AGN was last active. R1 shows a distinctly bent morphology, similar to a WAT radio source. The R2 observation in Figure 14 has a second pocket of emission in the lower lobe, which has not extended as far from the cluster center as the upper lobe. The radio emission shown in Figures 13 and 14 was contained within the boundaries of the X-ray cavities seen in Figures 11 and 12, respectively.

In addition to qualitative morphological evaluations that reveal simple dynamical relationships, the above synthetic observations can be examined quantitatively to good effect, as well. For example, Tregillis et al. (2002) used such observations at multiple frequencies to test observational methods for measuring magnetic field strengths and properties of the CR population. In addition to total intensity maps, our synthetic observation tool can also take advantage of the vector magnetic field and CR electron plasma information in the simulations to construct polarization maps that include Faraday rotation. These can be used, for example, to test techniques for measuring the structure of magnetic fields in the AGN and in the ICM surrounding it. We defer those studies to subsequent publications.

5. THE ROLE OF CLUSTER "WEATHER"

The asymmetries between the jet lobe and X-ray cavity pairs in these simulations clearly demonstrate that asymmetries in the properties of the ICM, that is, cluster "weather patterns," can strongly influence the morphology of radio lobes and ICM cavities. Since most radio galaxies in cluster environments show some degree of distortion, it is useful to understand how to establish the dynamical links to observed AGN outflow distortions as metrics of the ICM structure. In our jet simulations we can exclude direct gravitational forces as contributors to nonradial distortions of the AGN plasma and entrained ICM, since the gravitational potential in action during the simulations was spherically symmetric. Nonspherical gravity was essential, of course, in the development of the ICM conditions present at the start of our simulations, especially those we describe as sloshing.

A simple model can illuminate the main character of the AGN/ICM–weather interactions at work during our simulations. For this purpose we can picture the AGN outflows as being ejected into a wind at some angle to the jet axis. It is sufficient for our purposes to look at the two limiting cases having the winds aligned with the jet axis (i.e., "head" and "tail" winds) and orthogonal to the jets ("cross" winds). For simplicity, we assume a uniform ICM and steady jets. In doing this, we should remember that our simulated environment is nonuniform and the simulated jets are intermittent. Still, the model can provide some insights and consistency checks.

In the head and tail (aligned) wind case that applies approximately to our R2 simulation we expect the length of the jet structure, ℓ, on the tail wind side (v_w > 0 assuming v_j > 0) to be extended in comparison to a stationary ICM, while the head wind side (v_w < 0) will be shortened. A modified version of the classic calculation to balance the jet thrust against ram pressure from the ICM as the jet advances (Blandford & Rees 1974) provides a simple estimate of the wind's influence. The model estimates, if the local ICM density is ρ_ICM and the jet thrust is F_j, that the rate at which the jet "head" extends, v_h > 0, is

$\begin{eqnarray} v_h &\approx & v_{h0} \left[1 + \frac{v_w}{|v_w|}\left(\frac{\rho _{{\rm ICM}} v_w^2 \pi r_h^2}{F_j}\right)^{1/2}\right]\nonumber\\ &=& v_{h0}\left(1 + \frac{v_w}{v_j}\sqrt{\chi }\right), \end{eqnarray} \tag{ 10 }$

where r_h is the cross-sectional radius of the head, and

$\begin{eqnarray} &&v_{h0} = \frac{v_j}{1 + \sqrt{\chi }}, \end{eqnarray} \tag{ 11 }$

with χ = (r_h/r_j)²(ρ_ICM/ρ_j). In the applicable,"light jet" limit, $\sqrt{\chi }\gg 1$ this gives the intuitively obvious result,

$\begin{eqnarray} &&v_h \approx v_{h_0} + v_w. \end{eqnarray} \tag{ 12 }$

Setting the (observed) ratio of lengths, $\delta = \ell _+/\ell _- \sim (v_{h_0} + |v_w|)/(v_{h_0} - |v_w|)$ , one can find, noting that v_h0 is the average extension rate on the two sides, an estimate of the wind speed as, |v_w| ∼ v_h0(δ + 1)/(δ − 1).

At the end of our R2 simulation (t ≈ 200 Myr) one jet cavity extends roughly 300 kpc from the AGN, while the other reaches a distance only about 200 kpc from the AGN (see Figure 9). Thus, we can estimate δ ∼ 1.5. The average length is then about 250 kpc, which over the approximately 200 Myr life of the AGN activity corresponds to an average extension rate (∼v_h0) around 2500 km s⁻¹ (approximately Mach 5). Using the relation derived in the previous paragraph this would produce an estimate for the ICM wind speed, v_w ∼ 500 km s⁻¹. On these scales the actual ICM flow speed is around 400 km s⁻¹ (see Figure 3), so given the simplicity of the estimate, consistent with this result.

The orthogonal, cross wind case, related to our R1 simulation, is familiar in NAT and WAT contexts. In particular a transverse ram pressure from the wind will deflect the jet and its cocoon (e.g., Begelman et al. 1979; Jones & Owen 1979). In the absence of a wind the jets will deposit momentum, energy (including pressure), and AGN plasma in the head region. The energy and plasma will flow backward and expand in response to the added pressure to form an AGN plasma cocoon, which is represented by X-ray cavities in clusters. The cavities are reasonably close to pressure balance with the surrounding ICM (e.g., McNamara & Nulsen 2007, OJ10), so the rate of lateral expansion will be of the order of the ICM sound speed, c_ICM. The R1 cocoons at the end of the simulation have, on average, transverse radii, R_c ∼ 100 kpc, corresponding to mean expansion rates around 500 km s⁻¹. That is close to c_ICM and consistent with this expectation.

In a cross wind the cocoon plasma will experience a ram-pressure-based lateral force, F_w ∼ 2ρ_ICMv²_wℓR_c, where ℓ is the length of the cocoon and R_c is its transverse radius. The mass in the cocoon is roughly ρ_cℓπR²_c, with ρ_c the density in the cocoon. So, the time for F_w to accelerate the cocoon to the wind speed is roughly, t_wc ∼ (π/2)(ρ_c/ρ_ICM)(M_c/M_w)t_c, where M_w = v_w/c_ICM and t_c ≈ R_c/(M_cc_ICM) is the time for the cocoon to expand. With both M_c ∼ 1 and M_w ∼ 1 the cocoon will then reach terminal speed at rest in the ICM on a timescale less than it takes the coccon to expand. On the other hand, if the wind speed exceeds the expansion speed; namely, if M_w > M_c, the cocoon will be stripped from the jet, leaving a "naked" jet. Our R1 jets do possess cocoons, since M_w < M_c. Assuming that the cocoons move across the jet path at the wind speed, we can estimate the down wind displacement angle of the cocoon compared to the jet direction at its source, θ_d, as tan θ_d ∼ v_w/v_h. For the R1 structures at the end of the simulation, this gives tan θ_d ∼ 0.15–0.2, so θ_d ∼ 10° − −15°. That is consistent with the deformations visible on the left side of Figure 7.

For completeness we recall that a naked jet is deflected by a cross wind over a characteristic length, ℓ_d (e.g., Begelman et al. 1979; Porter et al. 2009):

$\begin{equation} \ell _d({\rm naked}) \sim \frac{\rho _j}{\rho _{{\rm ICM}}}\frac{v_j^2}{v_w^2} r_j \sim \frac{P_j}{P_{{\rm ICM}}}\frac{M_j^2}{M_w^2} r_j. \end{equation} \tag{ 13 }$

The expression on the far right emphasizes that the deflection length, ℓ_d of a naked jet in pressure equilibrium with its environment depends only on the jet radius and the relative Mach number in the jet flow to the Mach number of the wind. If the jet is cloaked it will be deflected by a transverse pressure difference across the cocoon resulting from F_w. The deflection length in this case can be estimated to be (Jones & Owen 1979),

$\begin{equation} \ell _d({\rm cloaked}) \sim \frac{\rho _j}{\rho _{{\rm ICM}}}\frac{v_j^2}{v_w^2} R_c = \frac{R_c}{r_j}\ell _d({\rm naked}), \end{equation} \tag{ 14 }$

so the bending will be more gradual for a given wind. In our case the expected ratio is roughly an order of magnitude consistent with the simulation results.

6. CONCLUSIONS

We have presented two MHD simulations of AGN jets inside of a galaxy cluster that was adapted from an SPH cosmological simulation to explore the effects of cluster "weather" on the jets and lobes. The important results from this work are as follows.

1.
Evidence that ICM flows in these simulations were sufficient to produce asymmetries between the opposing AGN jets and lobes was presented. Although pressure variations in the ICM (unrelated to the presence of a gravitational potential) and magnetic stresses were present, they were either not as large as the variations in the ICM motions or were very localized. The existence of significant "weather" in the cluster used for these simulations is interesting because it was specifically selected as a relatively relaxed system. Section 3 showed that typical flow speeds in the cluster were one the order of 400 km s⁻¹, which appear to be sufficient to modify the jets and lobes. This result is in agreement with the predictions of Hardcastle et al. (2005), who argued that flows as low as 100–300 km s⁻¹ could create bent WAT-like morphologies for jets similar to those in our model (high velocity and low density). These results also show that the AGN outflows act as probes of bulk motions in the ICM, even in a seemingly relaxed cluster. Simple dynamical models utilizing the influence of the ICM ram pressure on the AGN cocoons seem to give useful estimates of the ICM properties. This makes it possible to infer velocity structure information in environments where direct measurement is not possible.
2.
Magnetic field lines anchored in the ICM are stretched outward by the intermittent jets. This causes field lines to wrap around the jets and AGN plasma in the lobes, and in localized regions, the stresses associated with the magnetic field were on the order or even greater than the inertial stresses. Such regions serve as a reminder that even for the "high-β" plasmas modeled here, Maxwell stresses can significantly effect the evolution of a plasma and emphasize the importance of modeling these systems with MHD.
3.
Synthetic X-ray observations of these simulations show complex cavity structures associated with the intermittent jets and ICM motions. Multiple cavity pairs were "observed" as restarted jets pushed into cavities formed by previous outbursts. In one example given, a new cavity forming inside a larger cavity did not have an obvious companion cavity from the opposing jet. Bright rims were seen outlining the cavities, and they were associated with the injection of energy by restarting jets. Additionally, cavities were seen misaligned with the jet launching region due to the deflection of the jets by ICM flows.
4.
Synthetic radio observations made during periods of jet activity resembled DDRG sources with a bright pair of jets interior to fainter sources further out in the cluster. This state was temporary, however, as the CR electrons from the new outburst would eventually reach the end of the outer lobe. The ICM motions created a WAT morphology of the lobes in one simulation.

This work was supported at the University of Minnesota by the NSF grant AST0908668. P.J.M. was supported in part by the Graduate Dissertation Fellowship at the University of Minnesota. K.D. acknowledges support by the DFG Priority Programme 1177, the DFG Cluster of Excellence "Origin and Structure of the Universe," and the DFG Research Unit 1254. We thank the anonymous referee for help with improving the original manuscript.

MHD SIMULATIONS OF ACTIVE GALACTIC NUCLEUS JETS IN A DYNAMIC GALAXY CLUSTER MEDIUM

Article metrics

Permissions

Author affiliations

Dates

ABSTRACT

1. INTRODUCTION