AGN variability from 3D relativistic jets

Radio-loud AGN variability from three-dimensional propagating relativistic jets

Yutong Li Shandong Provincial Key Laboratory of Optical Astronomy and Solar-Terrestrial Environment, Institute of Space Sciences,
Shandong University, Weihai, 264209, China
Department of Physics, The College of New Jersey, 2000 Pennington Rd., Ewing, NJ 08628-0718, USA
Paul J. Wiita Department of Physics, The College of New Jersey, 2000 Pennington Rd., Ewing, NJ 08628-0718, USA Terance Schuh Department of Physics, The College of New Jersey, 2000 Pennington Rd., Ewing, NJ 08628-0718, USA Geena Elghossain Department of Physics, The College of New Jersey, 2000 Pennington Rd., Ewing, NJ 08628-0718, USA Shaoming Hu Shandong Provincial Key Laboratory of Optical Astronomy and Solar-Terrestrial Environment, Institute of Space Sciences,
Shandong University, Weihai, 264209, China
10 July 201830 August 2018ZZZ
10 July 201830 August 2018ZZZ
10 July 201830 August 2018ZZZ

The enormous sizes and variability of emission of radio-loud AGNs arise from the relativistic flows of plasma along two oppositely directed jets. We use the Athena hydrodynamics code to simulate an extensive suite of 54 propagating three-dimensional relativistic jets with wide ranges of input jet velocities and jet-to-ambient matter density ratios. We determine which parameter sets yield unstable jets that produce jet dominated FR I type radio galaxy morphologies and which tend to produce stable jets with hot-spots and FR II morphologies. Nearly all our simulations involve jets with internal pressures matched to those of the ambient medium but we also consider over-pressured jets and discuss differences from the standard ones. We also show that the results are not strongly dependent on the adiabatic index of the fluid. We focus on simulations that remain stable for extended distances (60 – 240) times the initial jet radius. Scaled to the much smaller sizes probed by VLBI observations, the fluctuations in such simulated flows yield variability in observed emissivity on timescales from months. Adopting results for the densities, pressures and velocities from these simulations we estimate normalized rest frame synchrotron emissivities from individual cells in the jets. The observed emission from each cell is strongly dependent upon its variable Doppler boosting factor. We sum the fluxes from thousands of zones around the primary reconfinement shock. The light curves and power-spectra, with red-noise slopes between and , so produced are similar to those observed from blazars.

galaxies: active – BL Lacertae objects: general – galaxies: jets – quasars: general

Paul J. Wiita

thanks: E-mail:, liy@tcnj.eduthanks: E-mail:

1 Introduction

Active galactic nuclei (AGNs) are powered by accretion onto supermassive black holes in the centers of their host galaxies. Around 10 percent of AGNs are classified as radio-loud (e.g. Jiang et al., 2007), and are characterized by relativistic plasma jets, the extended linear structures that transport energy and particles from the compact central region out to kpc or even Mpc-scale extended regions (Begelman et al., 1984; Peterson, 1997). Variability in observed emission can be considered a defining characteristic of AGNs, and for these radio-loud AGN the majority of this variable emission is understood to arise from synchrotron emission from relativistic flows of plasma along two oppositely directed jets (Urry & Padovani, 1995). The combination of multi-band radio observations and theoretical modeling has provided extremely strong evidence for the presence of both moving and standing shocks in these jets (e.g. Aller et al., 1985; Marscher & Gear, 1985; Hughes et al., 1991; Lister et al., 2001, 2009; Marscher et al., 2008, 2010). It is now generally accepted that the largest flares arise from the production and relativistic propagation of new components seen as radio knots (e.g. Hughes et al., 1991) but smaller observed variations can arise from changes in the overall directions, and hence Doppler factors, of the jets (e.g. Camenzind & Krockenberger, 1992; Gopal-Krishna & Wiita, 1992), turbulence within the jets (e.g. Marscher et al., 2008; Marscher, 2014; Calafut & Wiita, 2015), or differential velocities between small regions within the jets (e.g. Giannios et al., 2009; Pollack et al., 2016). Changes in the overall direction of the inner portions of the jet have clearly been demonstrated using very long baseline interferometry (VLBI) (e.g. Biretta et al., 1986; Ros et al., 2000; Piner, 2003; Caproni & Abraham, 2004; Lister et al., 2013). Blazars, the set of AGNs that is the union of BL Lacertae objects and Flat Spectrum Radio Quasars (or at least the high-polarization subset of FSRQs) are highly variable throughout the electromagnetic spectrum; their enhanced brightness and strong variability are understood to arise through Doppler boosting when one of the jets is pointing close to our line-of-sight (e.g. Urry & Padovani, 1995).

Extended extragalactic radio sources have long been classified into two categories (Fanaroff & Riley, 1974) based upon their radio morphology. Fanaroff-Riley I (FR I) sources haves jet-dominated emission, with the majority of their fluxes arising from the inner halves of their total extent and are weaker. The FR II, or classical double sources, have emission dominated by lobes containing terminal hot-spots; frequently only one of the jets feeding those lobes are detected. They can be distinguished based on power: objects below at 178 MHz were typically found to have FR I morphologies (where is the Hubble constant in units of 70 km s Mpc). An examination of the radio luminosity at 1.4 GHz against the optical absolute magnitude of the host galaxy led to the finding (Owen & Ledlow, 1994; Ledlow & Owen, 1996) that the radio flux boundary line between FR I and FR II sources correlates as , which means that more radio power is required to form a FR II radio source in a more luminous galaxy. Furthermore, some hybrid-morphology radio sources (HYMORS) have been discovered that show FR I structure on one side of the radio source and FR II morphology on the other (Gopal-Krishna & Wiita, 2000). These sources are important in understanding the basic origin of the FR I and FR II dichotomy, where the different morphologies may be induced by intrinsically different jet properties, interactions with different environments on either side of the source (e.g. Bicknell, 1994; Gopal-Krishna & Wiita, 2001; Massaglia et al., 2016), or long-term temporal variations combined with the time-lag in the observer’s frame between evolving approaching and receding lobes (e.g. Gopal-Krishna & Wiita, 2000, 2004; Massaglia, 2003; Wold et al., 2007; Kawakatu et al., 2009; Kapińska et al., 2017).

Hydrodynamical simulations of propagating jets are of critical importance to the understanding of radio galaxies and quasars and now have a history spanning four decades (Rayburn, 1977; Wiita, 1978; Norman et al., 1982). These simulations give fundamental support to the idea of the twin-jet models for radio galaxies dating to Scheuer (1974) and Blandford & Rees (1974). Two basic parameters, the internal beam Mach number, , and the jet to ambient density contrast, , were found to govern the morphology and propagation properties of these jets (Norman et al., 1982). Good numerical and analytical understandings of 2D and 3D non-relativistic hydrodynamical flows were achieved by the mid-1990s (e.g. Hardee & Norman, 1988, 1990; Burns et al., 1991; Hooda et al., 1994; Bassett & Woodward, 1995; Bicknell, 1995; Bodo et al., 1995). Since then, the complexity of the simulations has increased in parallel with growing computational power and algorithm development, leading to better understanding of the jet phenomenon, focussing on the study of the morphology, dynamics and non-linear stability of jets at kpc and larger scales (e.g. Norman, 1996; Ferrari, 1998; Carvalho & O’Dea, 2002a, b; Krause & Alexander, 2007; Donohoe & Smith, 2016; Massaglia et al., 2016).

Early simulations of special relativistic jets were carried out by a few groups. Some focused on the propagation on nuclear regions corresponding to the parsec-scale structures revealed by VLBI (e.g Martí et al., 1995, 1997; van Putten, 1996; Mioduszewski et al., 1997), while others focused on propagation out to larger scales and production of the structures seen on the multi-kiloparsec scale of powerful radio galaxies (e.g. Rosen et al., 1999; Hardee et al., 2001; Martí & Müller, 2003; Leismann et al., 2005; Massaglia et al., 2016).

Of course, extragalactic jets emit by the synchrotron mechanism (at least from radio through ultraviolet bands) and equipartition arguments have long indicated that the energy density in magnetic fields should be substantial. Nonetheless, there have been a more modest number of simulations that include them dynamically in magnetohydrodynamical (MHD) codes, particularly in 3D (e.g. Clarke et al., 1986; Nishikawa et al., 1997; Leismann et al., 2005; O’Neill et al., 2005; Keppens et al., 2008; Gaibler et al., 2009; Mignone et al., 2010; Huarte-Espinosa et al., 2011; Hardcastle & Krause, 2014; Martí, 2015; Martí et al., 2016). The relative paucity of work in this area can be attributed to both the difficulty of maintaining conservation of magnetic flux and energy in MHD codes (especially special relativistic ones) and to the early indications that simulations of jets with dynamically important magnetic fields did not resemble observed radio source morphologies (Clarke et al., 1986). We shall not discuss MHD jet simulations further in this paper, but note that the Athena code we use is designed to handle them and that we are in the process of producing a suite of such simulations to be compared to the relativistic HD ones presented here.

The Athena code, developed by J. Stone and colleagues (Gardiner & Stone, 2005; Stone et al., 2008; Beckwith & Stone, 2011) is an excellent choice for a wide range of astrophysical computational fluid dynamics problems. It is a highly efficient, grid-based, code for astrophysical MHD that was developed primarily for studies of the interstellar medium, star formation, and accretion flows. Athena has the capability to include special relativistic hydrodynamics (RHD) and MHD, static (fixed) mesh refinement and parallelization using domain decomposition and MPI. The discretization is based on cell-centered volume averages for mass, momentum, and energy, and face-centered area averages for the magnetic field. In order to solve a series of discretized partial differential equations expressing conservation laws, the rest density, pressure, velocity, internal energy, and magnetic field are calculated though in these RHD simulations the magnetic field was set to zero.

Here we simulate propagating jets with a very extensive suite of both medium-power and high-power jets that are run for substantial problem times using larger computational resources. These allow us to better distinguish between parameters likely to produce jets that remain stable for extended distances (FR II type) from those that do not (FR I type). The other key aspect of this work is that we examine the light curves and power spectra from approximated jet emission taking into account the different velocities of the cells within the jet. Although the literature on simulations of jets for extragalactic radio sources is now extensive, such computations of light curves (e.g. Marscher, 2014) and the resulting power-spectra produced by the jet motions (Calafut & Wiita, 2015; Pollack et al., 2016) are exceedingly rare. Here we consider emission variations that are produced by changing Doppler factors, densities and pressures from a fixed region in a propagating jet.

However, it has been noted (Pollack et al., 2016), that faster sub-grid variability is likely to involve relativistic turbulence (e.g. Zrake & MacFadyen, 2013) which is almost certainly present with these jets, at least in some portions (Marscher et al., 2008, 2010). Marscher (2014) has created a sophisticated “turbulent extreme multi-zone” model that posits the presence of turbulence and computes beautiful multi-band light curves for flows through a standing-shock within a jet of fixed bulk velocity. However, neither this TEMZ model, nor the simpler turbulent model of Calafut & Wiita (2015), take into account the variability arising from the bulk jet motions. Our group has previously performed relativistic jet simulations, albeit in only two dimensions, using the Athena code in Pollack et al. (2016), where wiggles in slab jets were computed in conjunction with a simplified turbulence model. The jet velocity and density changes produced slower variations while the turbulence produced faster ones that connected with them; these models had the advantage of producing variations over five orders of magnitude in time, and with power spectral slopes similar to those usually observed in AGN. There the changes in observed synchrotron emission were estimated by summing the fluxes from a vertical strip of zones behind the reconfinement shock and power spectral densities (PSDs) were calculated from the light curves for both turbulent and bulk velocity origins for variability; these yielded PSDs with slopes in the ranges to and to , respectively (Pollack et al., 2016). Here we extend the bulk calculations to 3D and include orders of magnitude more zones in the light curve estimations.

In Section 2 we provide a description of our RHD simulations. Our basic results are given in Section 3 for an exceptionally large number of jet simulations with a substantial range of jet velocities and densities; there we consider typical jet evolutions and compare higher and lower power inputs that respectively give rise to stable and unstable jets. In Section 4 we describe how we compute reasonable estimates for the light curves, focusing on the portion of the jet near the first reconfinement shock, where both intrinsic emissivity and Doppler boosting are high. In Section 5 we consider differences between pressure-matched and over-pressured jets, the effects of different adiabatic indices and give a brief comparison of 2D and 3D jets computed with Athena. In Section 6 we present and discuss light curves and power spectral densities (PSDs) for representative simulations. Our conclusions are given in Section 7.

2 Simulations of 3D RHD jets

We use the Athena code (Gardiner & Stone, 2005; Stone et al., 2008; Beckwith & Stone, 2011) for special relativistic hydrodynamics to produce 3D simulations of jets propagating through initially uniform external, or ambient, media, with wide ranges of powers. A relativistic hydrodynamic jet problem in three-dimensions is available in the problem files bundled with Athena111 Our jets are computed as launched along the -axis in cartesian coordinates, and the circular jet inlet is taken to have a radius unit (so that all lengths can be scaled to that unit) in the plane. To model a hydrodynamic jet, the initial physical parameters of inlet jet velocity, (assumed constant across the inlet cross-section for our initially cylindrical jets), proper (rest) ambient and jet densities ( and , respectively), ambient and jet pressures ( and ) and adiabatic index, , must be specified. In our production runs we normally use but we have performed a few runs with and will discuss the differences so produced. However, as has been shown in earlier work, the dominant variables are and . The great majority of our simulations are performed for pressure-matched jets () but we also consider a few over-pressured jets, with . Boundary conditions are set to outflow everywhere except at the inflowing jet inlet. We used a Courant number of 0.4 in our runs except when smaller values were needed for the code to function successfully. The constant density ambient media we model here are good approximations in both the nuclear region of a galaxy (scales pc) and on extended motions in intracluster media (ICM, on scales kpc). However, studies of such jets in media with decreasing densities, as would be appropriate for propagation through galactic halos (on 1 to 10 kpc scales) are also important (e.g. Gopal-Krishna & Wiita, 1991; Hooda & Wiita, 1998; Carvalho & O’Dea, 2002b; Jeyakumar et al., 2005; Massaglia et al., 2016) and we also neglect the likely important impact of large-scale motions in the ICM on propagating jets (e.g. Loken et al., 1995; Bliton et al., 1998). The Athena code permits the user to optimize the simulation by allowing a choice of the order of reconstruction (always taken as second-order piecewise parabolic in our production runs), as well as providing options such as static mesh refinement which can change the resolution of a grid by a factor of two per refinement level. This can allow for more focus on parts of the grid that contain interesting phenomena, such as the jet or shocks, although we did not employ mesh refinement.

With substantial experimentation involving different code parameters, our best overall results for faster jets came from simulations of higher resolution 3D RHD jets with zones with the length along the axis extending to 120 , and on for grids extending only out to 60 (or 10 zones in each length unit, hence 20 zones across the jet diameter). However, many preliminary runs and some production runs at those lengths were performed at a medium resolution (5 zones in each unit, or 10 across the jet diameter) including one with the jet propagating an extremely long distance (for 3D simulations, anyway) of 240 . For pc, as is reasonable for jets propagating on extragalactic scales, one time unit equals 1630 yr. However, if we consider these simulated jets to be propagating on VLBI scales, a reasonable value for might be 1 pc and a time unit in the source rest frame then would be 3.26 yr. These domain extensions, , and numbers of cells in each dimension, , for our models are given in Table 1 along with the corresponding values of jet velocities , Lorentz factor , and jet-to-ambient density ratios . Also given in Table 1 are the equivalent Mach number , jet luminosity , and the imposed limit time of each run . Results given in this table are the times that the jets start to become unstable , the distances that jets retain stability , whether or not the jets reach the ends of the grids (“End”), and the expected FR type. We also computed three overpressured simulations with ambient (dimensionless) pressure while the jet pressure (); the default pressure values are ). The three simulations that were also performed with overpressured jets are marked with asterisks in Table 1 and those also performed with a different adiabatic index are noted with daggers.

A total of 54 RHD simulations have been performed with the Athena code with different combinations of jet velocities () and jet-to-ambient matter density ratios (). The simulations contain a range of from 0.0005 to 0.0316 and a range of initial from 0.7c () to 0.995c (), where, as usual, the Lorentz factor is with . A summary of the results of these simulations are shown in Figure 1. The circles in Figure 1 represent runs with jets that eventually go unstable before the end of the grid at 60 or 120 is reached, and can be considered to correspond to FR I radio sources if scaled to extragalactic dimensions. Triangles show parameters of runs with jets that remain stable enough to retain terminal shocks throughout the simulation (to even 240 ) and thus are plausibly representative of FR II sources. The boundary between stable and unstable jets corresponds nicely to their nominal powers, as will be discussed further in Section 5.1.

Figure 1: Summary of the stability (and likely morphology) of our 54 simulations. Circles represent FR I runs which have unstable jets; triangles are FR II runs with stable jets.

One can compare Newtonian jets and special relativistic jets through defining and as the equivalent classical density ratios and Mach numbers, respectively (Rosen et al., 1999). In the formulation of Carvalho & O’Dea (2002a),


Here is a factor related to the energy density of relativistic electrons, and is defined as


where for a pair plasma jet (Case 1, Carvalho & O’Dea (2002a)), while , neglecting the thermal component, is typically taken for a proton-electron jet (Case 2); here is the average Lorentz factor of the relativistic electrons, and , the ratio of the number of electrons to total particles within the jet, is taken to be 0.5. For a limiting fully relativistic adiabatic index , is 0.106 (Case 1 in Carvalho & O’Dea (2002a)). For our calculations, we have nearly always employed a value of appropriate for an ionized monatomic gas; then , in approximate accordance with . We include in Table 1.

The following expression gives the energy flux of a jet where relativistic particles are important (Bicknell, 1994, 1995)


where is the jet pressure and is the ratio of cold matter energy density to enthalpy and is the internal energy density; when the internal energy and pressure are dominated by relativistic particles then (Bicknell, 1995). We use this expression, taking a fully-ionized gas mixture of hydrogen and helium so the mean molecular weight, and then , where is the number density of the ambient gas. In the limit where the cold material dominates over the relativistic particles, Equation 4 reduces to


The simulations are carried out using dimensionless variables. However, it is sensible to present the results with scalings to physical units in order to facilitate comparison with astrophysical radio sources. Thus, if we consider large-scale jets on scales kpc propagating through a uniform intracluster medium (ICM) we can take kpc and can take ICM number densities, temperatures and pressures to have typical values around , K and thus dyne (Sarazin, 1988). The values of for these parameters are given in Table 1, where a range from erg s to erg s is spanned. In Section 4 we consider variability for which jets propagating through the nuclear core of the galaxy (which is also nearly constant in density inside the core radius) are relevant instead; then, a sensible scale is set by taking  pc and for the ambient gas. The equivalent jet powers for the nuclear region are essentially reduced by a factor of when compared to the large-scale extragalactic ones for the same parameters.

0.7 1.40 0.01 1000 0.4 1.98 8.05(44) 700 24.0 yes FR I
0.8 1.67 0.00316 1100 0.4 2.70 5.82(44) 550 18.5 no FR I
0.8 1.67 0.01 1300 0.4 2.70 1.82(45) 325 25.8 yes FR I
0.8 1.67 0.0316 800 0.4 2.70 5.71(45) 275 36.8 yes FR I
400 0.4 3.27 5.90(45) 240 32.7 yes FR I
0.9 2.29 0.005 1000 0.4 4.17 2.75(45) 350 58.7 no FR I
0.9 2.29 0.01 1000 0.4 4.17 5.45(45) 425 64.2 yes FR I
0.9 2.29 0.02 600 0.4 4.17 1.09(46) 400 111.5 yes FR I
0.9 2.29 0.0316 300 0.4 4.17 1.72(46) yes FR II
0.93 2.72 0.00316 500 0.3 5.12 2.83(45) 100 16.6 no FR I
0.93 2.72 0.005 500 0.4 5.12 4.45(45) 225 31.6 no FR I
0.93 2.72 0.01 500 0.4 5.12 8.88(45) 240 51.2 no FR I
0.95 3.20 0.001 1100 0.3 6.15 1.35(45) 300 31.1 no FR I
0.95 3.20 0.00316 600 0.4 6.15 4.35(45) 170 31.5 no FR I
0.95 3.20 0.005 400 0.4 6.15 6.85(45) 350 57.7 no FR I
0.95 3.20 0.01 500 0.4 6.15 1.37(46) yes FR II
0.95 3.20 0.02 500 0.4 6.15 2.73(46) yes FR II
0.95 3.20 0.0316 500 0.4 6.15 4.30(46) yes FR II
0.96 3.57 0.0075 400 0.4 6.93 1.35(46) yes FR II
0.96 3.57 0.01 400 0.4 6.93 1.80(46) yes FR II
200 0.4 6.93 3.58(46) yes FR II
0.96 3.57 0.0316 300 0.4 6.93 6.15(46) yes FR II
0.97 4.11 0.00316 800 0.2 8.06 7.90(45) 350 61.5 yes FR I
0.97 4.11 0.005 550 0.4 8.06 1.27(46) 380 80.8 yes FR I
0.97 4.11 0.0075 500 0.4 8.06 1.90(46) yes FR II
0.97 4.11 0.01 500 0.4 8.06 2.53(46) yes FR II
0.97 4.11 0.02 1000 0.4 8.06 5.05(46) yes FR II
0.97 4.11 0.0316 500 0.4 8.06 7.98(46) yes FR II
0.98 5.03 0.00316 600 0.3 9.97 1.27(46) 400 80.7 yes FR I
0.98 5.03 0.005 600 0.2 9.97 2.01(46) yes FR II
0.98 5.03 0.0075 800 0.4 9.97 3.03(46) yes FR II
400 0.4 9.97 4.03(46) yes FR II
200 0.4 9.97 8.05(46) yes FR II
0.98 5.03 0.0316 400 0.4 9.97 1.27(47) yes FR II
0.985 5.80 0.0075 300 0.4 11.56 4.18(46) no FR II
0.985 5.80 0.02 500 0.4 11.56 1.11(47) yes FR II
0.985 5.80 0.0316 500 0.3 11.56 1.75(47) yes FR II
0.99 7.09 0.0005 600 0.1 14.20 4.32(45) 270 52.5 no FR I
0.99 7.09 0.001 600 0.1 14.20 8.63(45) 280 57.9 no FR I
0.99 7.09 0.00316 500 0.3 14.20 2.73(46) yes FR II
200 0.4 14.20 6.53(46) yes FR II
0.99 7.09 0.01 150 0.4 14.20 8.70(46) yes FR II
0.99 7.09 0.02 400 0.3 14.20 1.73(47) yes FR II
0.99 7.09 0.0316 400 0.3 14.20 2.73(47) yes FR II
0.992 7.92 0.0075 300 0.4 15.89 8.30(46) yes FR II
0.992 7.92 0.01 300 0.4 15.89 1.11(47) yes FR II
0.993 8.47 0.001 500 0.2 17.01 1.27(46) yes FR II
0.993 8.47 0.00316 400 0.3 17.01 4.23(46) yes FR II
0.993 8.47 0.0075 300 0.4 17.01 9.58(46) yes FR II
0.993 8.47 0.02 300 0.3 17.01 2.55(47) yes FR II
0.994 9.14 0.0075 300 0.3 18.38 1.12(47) yes FR II
0.994 9.14 0.0316 300 0.2 18.38 4.73(47) yes FR II
0.995 10.01 0.01 300 0.3 20.15 1.81(47) yes FR II
0.995 10.01 0.0316 300 0.1 20.15 5.73(47) yes FR II
Table 1: Simulation parameters
  • 1. is the prescribed time limit of the simulation.

  • 2. gives the extent of the domain in units of .

  • 3. means numbers of grid points in the three directions.

  • 4. The Courant, Friedrichs, & Lewy (CFL) Number

  • 5. in erg s

  • 6. is the time that jet becomes noticeably unstable, if it does.

  • 7. is the distance out to which the jet remains stable, if it goes unstable.

  • 8. indicates whether or not the bow shock reaches (or goes beyond) the end of the grid before .

  • The run was also computed for an overpressure: .

  • The run was also computed with .

Figure 2: Snapshot slices of fluid densities (on a logarithmic scale) for , = 0.0075, at different stages: (left) ; (center) ; (right) . There are 10 zones per jet radius and zones to 60 along the -axis and to along the - and -axes.

3 Results

Snapshots of fluid density slices (on a logarithmic scale) in the plane for a stable jet with and at three different time steps are shown in Figure 2. All of the ambient densities are taken as fixed and the same (a non-dimensional value of 10.0, denoted by a red color), in our simulations. The length of this simulation, and roughly half of our others, is 60 , and the time steps shown are 60, 120, and 190. For the three images the scale is such that the maximum density (purple) is in the bow shock and is times that of the ambient density (red) and the minimum density (white) is of the ambient density, or a little less than half that of the inflowing jet matter.

The highest density region is the bow shock advancing into the ambient medium, and behind it is the shocked ambient medium which appears in shades of red, orange and yellow in these log plots. A contact discontinuity separates this material from the cocoon of shocked jet material which is left behind as the jet advances and appears in shades of green and blue. Substantial vortices form within the cocoon and these, along with Kelvin-Helmholtz instabilities, facilitate the mixing of shocked ambient and shocked jet fluids. A series of strong internal shocks and rarefactions form within the jet, as can be seen by the changes from the blue color corresponding to the injected jet density to the turquoise and white colors, respectively. The first such shock is called the reconfinement shock (even though subsequent ones also help keep the jet confined) and it maintains a relatively constant location after the jet achieves a length some ten times its radius. In this case, the jet remains well collimated with a terminal shock (corresponding to an observed radio hot spot) until the end of the problem grid is reached. It is thus classified as a stable source, with a cocoon similar in shape to many FR II RGs. For these simulations we examine 1D plots of the density or pressure along the jet axis and for the ones we call stable the strongest shock in the jet itself occurs at this termination point, which remains very close behind the bow shock. At early times the jet is essentially axially symmetrical, but as it propagates outward, small numerical fluctuations grow and modest asymmetries have clearly arisen before the jet leaves the grid. When we compare simulations with the same input parameters made at moderate resolutions (10 zones across the jet diameter) with our higher resolution runs (20 zones per diameter) we see few significant differences in the distances traveled at given times, the stability, or the gross structure of the shocks; however, the higher resolution runs certainly resolve more features of the internal jet shocks and better display the eddies within the cocoon.

Figure 3: 3D jet simulation with parameters and =0.0316 (high resolution, at ): (top) longitudinal slice at of the logarithms of the fluid densities; (bottom) pressure (also on a logarithmic scale) shown between , and magnitude of the velocity and density (also on logarithmic scales), shown for , with velocity shown for (including a quarter-cross section at and density for ; in addition, most of a cross-section of for is shown.

Another high resolution simulation of a somewhat slower, but denser, jet (with ) is shown in Figure 3; in this case the simulation ran to a distance of 120, and the plots are made shortly before it reaches that distance, at . The upper panel shows a log slice similar to that shown earlier for the other jet in Figure 2. The lower-panel of Figure 3 presents different cross-sections of the pressure, density and velocity distributions. The top plot shows that the jet in the center remains clearly stable for over 70% of its length, but it starts to wiggle a bit near the end of this 120 grid. This means it may eventually transition into an FR I morphology at an even later time, though it might also appear as an FR II with hot-spot offset from the current jet direction; such morphologies are common, including in the canonical FR II, Cygnus A (e.g. Carilli & Barthel, 1996). The left portion of the bottom combined plot illustrates the typical pressure structures for these strong, pressure-matched, jets, where both the ambient and initial jet pressures show as a light blue shade while a much higher pressure (in red) is seen in the bow shock and much of the shocked ambient gas (red to yellow). The shocked material that passed through the jet and fills the radio emitting cocoon contains only modest pressure variations (shades of blue) and does not differ greatly from the average internal jet pressure. Within the jet itself, pressure variations corresponding to shocks and rarefactions are visible. The logarithm of the magnitude of the velocity is given in the upper-right portion of the lower panel of Figure 3. The unshocked ambient medium, which doesn’t move, is shown as white and the shocked ambient medium, which is moving very slowly () appears as yellow and green. The cocoon moves at modest speeds (up to ), but the only strongly relativistic motions are seen within the jet itself, where the average velocity stays very close to the injection speed of (purple). The cross section of density shows a nearly axisymmetric distribution, particularly in the shocked ambient material, but modest perturbations have appeared within the cocoon material; the apparent dominance of an mode at this stage is almost certainly attributable to the use of cartesian coordinates to model a cylindrical jet.

Figure 4: High resolution simulations shown at or very close to 60 : (top left) for at ; (top right) , at ; (middle left) for at ; (middle right) for at ; (bottom left) for at ; (bottom right) for , also at .

Four different high resolution simulations extending out to 60 are shown in Figure 4. The top four plots are the snapshots of slices of the fluid densities (on a logarithmic scale) with the input parameters: , (top left); , (top right); , (middle left); and , (middle right). The bottom panels give the distribution of pressure (bottom left) and velocity (bottom right) for that last simulation. There are clearly different morphology evolutions displayed. For , , the jet starts to become unstable at the distance of 24.0 which has been listed in Table 1 and as it continues to develop it would certainly produce an FR I morphology on extragalactic scales, particularly when ambient material is entrained in the jet and renders it transonic (e.g. Bicknell, 1994, 1995). We define the position at which the jet terminates within the lobe as the strongest density enhancement measured along the jet axis. In these unstable simulations, any additional density enhancements further out toward the hotspots are weaker and visual inspection of 2D density plots shows they involve swinging or effectively fragmented flows. In such relatively weak jets the bow shock actually can fade into a subsonic pressure wave as is seen in the top left panel of Fig. 4. The second case, with and , is essentially stable out to this distance, with the strongest jet shock still close behind the bow shock, but some wiggles are developing in the jet and it is possible that it will become unstable shortly. The other two simulations remain fully stable while propagating past the end of the 60 grid and would exhibit FR II morphology; these differences will be discussed further in Section 5.1.

The pressure distribution at the bottom left of Figure 4 shows that the pressure at the head of bow shock is much higher than that in the jet material but a high pressure terminal shock is also present at the end of the jet just behind the bow shock, along with three internal shocks, which would likely be viewed as knots within the jet. The velocity structure shown in the lower right panel makes it clear that the jet velocity remains highly relativistic until the terminal shock and that the flow in the cocoon is much slower. Hence, even though the cocoon volume is much higher than that of the jet, if the source is viewed close to the line of sight the much greater Doppler boosting of the jet would mean that emission from the jet would dominate the observed brightness.

An even more powerful jet, with , is shown to remain beautifully stable all the way out to 240 in Figure 5. At this high power the cocoon is very narrow and it would not resemble the typical FR II extragalactic radio source morphology. However, if scaled down to galactic nuclear scales it is reminiscent of the morphology of many VLBI jets (e.g. Martí et al., 1995, 1997; Mioduszewski et al., 1997).

Figure 5: Density slice for a simulation extending to 240 performed at moderate resolution for input parameters of at .

4 Estimation of Emission Brightness

Any radiating material traveling at relativistic speed will be affected by Doppler boosting, which strongly beams the emission in the direction of its velocity, providing an amplification in brightness and a shortening of the observed time over which a variation occurs (e.g. Urry & Padovani, 1995; Gopal-Krishna et al., 2003). The flow in the entire jet wiggles back and forth by a small amount and there are certainly differential motions within the jet. Hence, the angle of the velocity of each cell with respect to the observer is always changing, resulting in the apparent rise and fall of the received flux from individual cells within the jet region. In addition, the density and pressure of each cell varies with time. The Doppler factor for a cell characterized by coordinates and numerated by is given by


where is the Lorentz factor of the flow in that voxel, is its speed in units of , and is the angle of the velocity of that cell with respect to the observer. Within our simulations we have the vector velocity and values of density and pressure in every cell and for an assumed value of for the line-of-sight angle to the jet axis we can find .

In principle it would be optimal to include boosted emissivity from all zones in the jet and cocoon, but this is computationally overwhelming as it must be done for thousands of time intervals in post-processing to produce a light curve. The regions we used for computing the source brightness were determined visually after examining most of the the entire suite of simulations, since only zones moving at high velocities will obtain strong boosting and only zones with higher densities and pressures are expected to have high rest-frame emissivities. A location close to the (first) reconfinement shock is optimal (Pollack et al., 2016) and is illustrated in Figure 6. The explicit range in which we chose to compute emission is 12 in length, and 2 in width and height (). Therefore, the total number of zones used for computing the emission are 6000 in medium resolution runs (5 zones in each length unit) and 48000 for high resolution runs (10 zones in each unit). The input parameters of the example shown in Figure 6 are , =0.0075, at (the density for moderate resolution simulation is illustrated). We remove much, but not all, of the initial transient burst in emission produced by the formation of the jet and the passage of the bow shock by only considering light-curves to begin once the recollimation shock had formed and had a nearly stable location (taken to be from in our simulations) and we end the computations once the simulation was recomputed with the dimensionless time extended to 1000 or more units. Doing so meant that the jets usually propagated beyond the right-hand boundary of the grid, but any loss of material from the grid would have essentially no impact on variations occurring in the inner quarter or eighth of the entire simulation. We exported data at intervals of 0.25 time-units, thus typically producing 3400 output sets of velocity, density and pressure for the relevant zones between , which was sufficient for use in the calculation of power spectra spanning about three decades in frequency.

The basic emitted flux per unit volume can be considered arbitrary for our purposes as we are only interested in the relative variations, so it was set to unity for an individual zone. The physical mechanism producing the jet emission from the radio through (at least) soft X-ray bands is synchrotron emission, the intensity of which scales with the number of relativistic particles and the square of the strength of the magnetic field. As we have not incorporated the microphysics necessary to compute the fractions of ultrarelativistic electrons in each zone and their distributions in energy, we assume that the zone rest mass density can be taken to be proportional to the number of relativistic electrons that produce synchrotron emission in that zone. We also scaled the emission by the pressure in that zone, as the pressure, , can be reasonably taken to be proportional to the square of the magnetic field strength, though of course we have no information on the orientation of that “field”, so we can say nothing about the polarization of the radiation. For each zone at each output temporal step, we use the velocity data to calculate a Doppler boosting factor based upon an assumed viewing angle to the -direction along the initial jet direction. The resulting observed flux was then summed along with the observed fluxes of the other zones within the range from that time-step to yield a total observed estimated flux for that time. We do not take into account the temporal lags for receipt of emission from different zones in the jet as has been done in models or simulations that involve many fewer zones (e.g. Marscher, 2014; Calafut & Wiita, 2015; Pollack et al., 2016), as this would demand much greater computational resources and is not expected to produce significant differences to the magnitude of amplitude fluctuations nor to the shape of the power-spectra. With these approximations, the observed flux from this key portion of the jet is estimated as


where is the Doppler boosting factor, is the pressure in each zone, and is the density. The observed flux is proportional to the () power of the Doppler boosting factor, where is the slope of the synchrotron spectrum (with ) and is 2 for continuous flow and 3 for a shock in a flow (Begelman et al., 1984). We took to be equal to a typical value of 0.5 for a relativistic jet, and was set to 2 in this work.

Figure 6: The selected emission region, chosen to be close to the reconfinement shock. =0.99c, =0.0075, t=130.

5 Discussion

We have carried out a series of simulations to explore the effects of the different parameters on the jet propagation, evolution, and resulting morphologies. As shown in Section 3, snapshots of density, pressure and velocity retrain almost axisymmetric distributions for some time but eventually numerical perturbations seed a range of instabilities even without any prescribed variations in jet power or direction. For weaker sources these instabilities basically disrupt the jet and while the bow shock continues to advance without a coherent jet making it to the outer part of the lobe and this can produce a FR I type morphology. Most of our simulations were sufficiently powerful so that they indicated essentially stable jets which retain a strong terminal jet shock for extended times; these look like FR II radio sources, particularly for the lower values of . We also computed a few simulations that were over-pressured, with and in this section we also consider the differences in morphology between 3D and 2D (slab-like) simulations.

5.1 Distinction between FR I and FR II types

Although our focus has been on the propagation of powerful, stable jets, we did consider a range of jet powers that spanned a factor of 985 in order to study the distinction between FR I and FR II morphologies. To illustrate some of the key differences, plots of the distances against time of both the bow shock position along the x-axis (or peak of the wave front separating uncompressed ambient medium from that affected by the jet) and jet terminus (the strongest shock within the jet) for nine different runs have been shown in Figure 7. These are plotted with different colors for different input parameters, marked as c1 – c9 on the figure, with properties listed in Table 2. Solid lines are the spread of bow shocks and triangles are jet termini, (as defined in Sect. 3) which always start out close to the bow shock. The first two belong to FR I morphologies and their slower growth puts them on the right side of Figure 7, and the others are FR II, as shown on the left 7 curves.

Figure 7: Solid curves are the leading edges of the bow shocks, and the triangles give the jet termini behind the bow shocks for 9 simulations, labeled c1 to c9; parameters for these runs are listed in Table 2.
Run c1 c2 c3 c4 c5 c6 c7 c8 c9
0.7 0.8 0.96 0.97 0.98 0.985 0.985 0.99 0.993
1.40 1.67 3.57 4.11 5.03 5.79 5.79 7.09 8.47
0.01 0.01 0.02 0.0075 0.0316 0.02 0.0075 0.0075 0.0075
8.05(44) 1.82(45) 3.58(46) 1.90(46) 1.27(47) 1.11(47) 4.18(46) 8.70(46) 9.58(46)
Table 2: Parameters for bow-shock and jet terminus distance plots
Figure 8: Comparison between standard (left, ) and overpressured (right, ) simulations close to 60 . Top row: , =0.02 (); middle: =0.02 (); bottom: =0.0075 ().

It is quite clear that the jets in cases c1 and c2 terminate further and further behind the bow shock once they go fully unstable; at this time ( for c1 and for c2) the terminal location of the jet, which is the strongest shock within the flow, essentially stalls. Nonetheless additional energy, momentum and mass continue to flow through the jet into the cocoon, which still has sufficient pressure for some time to inflate a bow shock. However, because those inputs are spread over a widening cocoon the rate of shock advance slows substantially and may eventually decay into a weaker bow wave. For runs c3–c9 the terminal jet shocks which continue to be the strongest within the jets (and would be observed as hot spots) remain very close to the bow shocks throughout the simulations. The rate of advance of the bow shock remains remarkably constant for nearly all of these FR II like runs, although the least powerful of these (c4) actually accelerates at , where the bow shock narrows and allows a more focused thrust.

In summary, unstable FR I type morphologies are produced by slower and/or lower density jets while faster/higher density ones yield FR IIs. The critical variable is the jet power and from Table 1 we see that (at least out to the distances we simulated) and for the particular assumptions about constant jet radii and ICM properties we considered, the boundary between FR I and FR II sources is around erg s.

Unfortunately, there are no iron-clad ways to use observations to constrain values from the measured radio fluxes (for FR IIs) or X-ray cavity measurements (particularly helpful for FR Is) but attempts to make such connections yield lower values than those given in our Table 1. For instance, Godfrey & Shabala (2013) estimate jet power values for a substantial sample of FR II sources to range between and erg s while Cavagnolo et al. (2010) estimate those for FR I sources to range between and erg s. A reasonable boundary between these types appears to lie around erg s (Godfrey & Shabala, 2013). Our simulations have focused on the higher power sources that might yield FR IIs, so even lower and values than those we considered would certainly yield FR I types of much lower power, even if we kept the physical quantities to which we were scaling our simulations the same. In addition, we could have certainly chosen lower, but still reasonable, values for or with which to scale our simulations and doing so would have improved the match with these observational estimates. We note that there is a large spread, up to a factor of , in the jet powers that are estimated to produce the same observed radio flux for FR I sources, while the relation for FR II sources is tighter, but still spans nearly a decade (Cavagnolo et al., 2010; Godfrey & Shabala, 2013).

Massaglia et al. (2016) recently performed several 3D HD simulations of supersonic jets propagating through media with declining densities and pressures to investigate the morphologies of low power jets. They showed that the jet energy in the lower power sources, instead of being deposited at the terminal shock, was gradually dissipated through turbulence. The jets spread out while propagating, and they smoothly decelerate while mixing with the ambient medium and produce the plumes characteristic of FR I objects. These simulations confirm the early analytical work that showed that propagation on modestly relativistic jets through a declining atmosphere (usually on scales of kpc) should induce a jet slow-down to transsonic speeds and a substantial increase of the jet opening angle (Gopal-Krishna & Wiita, 1991; Bicknell, 1994, 1995). While these effects are almost certainly important and might help explain the discrepancy between the power boundary for our FR I and FR II simulations and the lower one estimated from observations, we have here shown that FR I type sources nonetheless can form for relatively weak relativistic jets, even in ambient media taken to have constant density and pressure.

Figure 9: Density slices for of simulations when the bow-shocks reach for in the upper row and in the lower row: (left) , =0.02 (); (right) , =0.02 ().

5.2 Overpressured simulations of powerful jets

In all of our simulations mentioned so far (as well as those in the vast majority of the jet simulation literature cited in Section 1) it has been assumed that the jet is initially in pressure equilibrium with the ambient medium, with the ratio of jet pressure and ambient pressure . This assumption is motivated by the fact that an overpressured jet develops a strong recollimation shock that leads to the development of equal pressures between the internal and external gas (e.g. Belan et al., 2010). This equilibration probably occurs at radii within the initial 100 pc of the jet length and so assuming that the jet is already in pressure equilibrium at injection is quite sensible for sources propagating through the ICM or an even lower density IGM. However, on smaller scales, the jets may well be overpressured and so we have considered three such runs which were simulated with , and have been marked with asterisks in Table 1.

Density snapshots at late times for these overpressured simulations are shown in the right column of Figure 8, where the equal pressure runs with the otherwise identical input parameters are displayed in the left column of Figure 8 at the same times. In all cases the morphologies of the simulation are similar, particularly for the bow shocks and cocoons. There are, however, a few modest differences that can be observed: the overpressured simulations propagate a slightly longer distance than do the equal pressure ones at the same times, their jets are wider close to the inlet and their first recollimation shock occurs further down the jet. All of these differences would be expected because of the higher jet pressures.

5.3 Effects of varying adiabatic index

An optimal treatment of the equation of state (EoS) in RHD simulations would allow for smooth transitions between fully relativistic jet fluids, where the adiabatic EoS, , is a good approximation, and the hot, but still essentially classical, fluid in the unshocked ambient medium where is highly accurate. Various computationally efficient approximations to the complex analytical (Synge, 1957) fully relativistic EoS for non-magnetized fluids have been proposed (e.g. Mignone et al., 2005; Ryu et al., 2006) and this transitional EoS approach was extended to RMHD flows by Mignone & McKinney (2007). Unfortunately, the Athena code is not designed to accommodate different EoSs for the different fluids, nor can it currently handle a transitional EoS. Hence we have computed three simulations where is assumed throughout to see if consistent differences arise in the results for a different adiabatic EoS.

The full Synge EoS was used in Scheck et al. (2002) where they analyzed three 2D RHD jet simulations with different compositions (leptonic or baryonic) and found very little difference between the morphologies and dynamic behaviors of those light ( or ) jets. Mignone & McKinney (2007) presented a different EoS comparison, also for 2D RHD simulations, for one very fast () and low density () set of simulations where they found the jet to propagate substantially faster than the one. They also computed a transitional EoS simulation which, unsurprisingly, led to results in between those of the two fixed values.

In Figure 9 we display late-time density slices of two of our 3D RHD high-resolution runs where we used in the upper row and in the lower row. The first run () has become unstable by the time the jet goes somewhat beyond , while the second run () remains stable out to at least . There are no huge differences between the and simulations in terms of basic shape and stability, as was also seen for a simulation with which we do not display. However, there are interesting distinctions. These 3D simulations with move faster across the grid by modest factors of than those with . This greater advance speed implies stronger and narrower bow-shocks and slightly narrower cocoons when compared at equal distances (as shown in Figure 9), but when compared at equal times, while the bow-shocks remain narrower, the cocoons, at least in the outer-halves of the sources, are of similar width. It is worth noting that the differences between jets with alternative adiabatic indices in our 3D simulations are in the same senses as those seen earlier for 2D jets (Scheck et al., 2002; Mignone & McKinney, 2007), but are not as pronounced in 3D. It also appears that the jets in the weaker simulations may remain stable for modestly longer distances and times.

5.4 Comparison of RHD jets in 3D and 2D

Our work here on 3D RHD simulations follows our previous 2D RHD relativistic jet simulations (Pollack et al., 2016) which were focused on the variability induced by jet motions. Even though the 3D runs are clearly superior and were enabled by The College of New Jersey’s (TCNJ’s) recent acquisition of a large computer cluster, it should be worthwhile to compare the morphology of jets produced in 2D and 3D RHD simulations using the Athena code at the same high resolutions. As shown in Figure 10, there are major differences between 3D (left) and 2D-slab (right) simulations regardless of input parameters. The top row shows a weaker jet that goes unstable rather quickly and the lower one a more powerful jet that remains stable at the times when their bow shocks reach the end of the grid at . This critical distinction is seen in both 3D and 2D simulations, with the jet going unstable at roughly the same distance in both of the lower-power cases. The morphology of 2D simulations is more symmetric than 3D ones as fewer instabilities can be excited in the former. The big difference is that 2D simulations inflate much wider bow shocks and cocoons and thence lose material off the grid along the upper and lower boundaries. The jet widths are similar, though slightly wider for the 2D runs at large distances. Concomitantly, 2D simulations take longer times to cross the entire grid and their jet termini are much further behind the bow shock than they are in the 3D simulations.

Figure 10: High resolution RHD runs for 3D and 2D slabs at the same distances: (top-left) 3D simulation with , = 0.01 at ; (top-right) 2D with the same parameters for ; (bottom-left) 3D with , = 0.0075, ; (bottom-right) 2D with the same parameters, .

6 Light curves and power density spectra

A key new aspect of this work is our consideration of the variations of emission produced by 3D propagating relativistic jets. To encompass variability time-scales of practical interest, i.e., those that are shorter than decades, we must scale the jets downward so that they correspond to jets emerging from the vicinity of the central engine where the length scales are parsec-like and not kpc-like and where Doppler boosting shortens the observed times of variations when compared to those occurring in the rest frame of the jets. Our viewing angles to the relativistic jets in AGN we call blazars are quite small and so we have usually taken in computing crude light curves as described in Section 4. Observed light curves for an extended run (to in the source rest frame) with are shown in Figure 11 for different viewing angles: and . Smaller viewing angles produce higher observed fluxes over the shorter periods in the observer’s frame. The emission for this case when viewed at has an amplitude about 25 times higher than that seen at thanks to the higher boosting factor.

Using the nuclear scaling ( pc) the rest frame time unit is 3.26 yr and the interval at which the light curves were sampled is yr. For the range of bulk Doppler factors considered in Figures 11 and 12, of 2.3 to 9.2, this means that time intervals down to 0.089 yr (essentially one month) are being simulated. The nominal intrinsic powers for such jets moving through the galactic nuclear region lower by a factor of than those quoted in Table 1, where they were scaled for propagation through the ICM, yielding a rescaled range from erg s to erg s. However, here we are fundamentally concerned with seeing if the fractional levels of variability produced by changes in the Doppler boosting of the different cells of the jet in the region close to the first reconfinement shock, and the resulting power-spectra, are in reasonable agreement with those seen in blazars. We note that the values of and , which set the nominal through Eqn. (5), can all have significant ranges about our nominal values.

Figure 11: Light curves for , = 0.02 at different viewing angles: (left) ; (middle) ; (right) = .

For several emission light curves, all explicitly evaluated at a viewing angle of , the power spectral density (PSD) were computed using a Lomb-Scargle periodogram, which is a method to estimate the PSD generated by a physical process that can sensibly deal with unevenly sampled time series data (Horne & Baliunas, 1986), although they are computed at equal intervals here. The periodogram is evaluated using the algorithm presented in Press & Rybicki (1989) in order to achieve fast computational speeds. The observed PSDs of AGNs have either been measured in X-ray bands for Seyfert galaxies or in optical bands for quasars and blazars. In either case, they are usually found to be approximated by red noise at lower frequencies, with the power depending on frequency roughly as with , while at high frequencies, measurement white noise () dominates; in a small minority of cases, a flatter slope, , may be present below a relatively low break frequency (Uttley & McHardy, 2005; González-Martín & Vaughan, 2012; Edelson et al., 2013; Wehrle et al., 2013; Revalski et al., 2014; Smith et al., 2018).

Four examples of our simulated emission light curves (left) and their corresponding PSDs (right) are shown in Figure 12. For the first two cases, both with the same value of jet velocity () but different jet densities, these light curves show overall substantial variations (55–75%) with some rapid fluctuations. During the first 150 (proper) time units of the 1000 for which these were run, the bow-shock and terminal shock move past the portion of the grid that we analyze (which is shown in Fig. 6), so we do not consider that portion of the light curve. The decreasing fluxes seen in most of the early portions of the simulations as plotted arise from the transient effect of the establishment of the the reconfinement shock, whose position and strength often takes a while to stabilize. After that phase, the variations are more modest (on the several percent level) and reflect the changes in the velocities, densities and pressures of the relevant jet zones once a quasi-steady state is established in the inner portion of the jet. But there are stronger longer-term variability trends for the faster jets shown in the lower two panels (, and ), both for the same . In both of these cases the baselines can be considered to be composed of two parts: a declining trend for the first 20–30%, and a constant one for the remaining portions of the light curves. Before proceeding to compute the PSDs of the emissions, it probably makes sense to subtract those baseline trends, but we have computed the PSDs with and without subtracting the trends and give the results both ways in Table 3. The top panels of the and simulated light curves show the full light curves and the baselines while the bottom panels show the light curves after subtraction of the best-fit line trends for each simulation.

Figure 12: Observed light curves (left) and PSDs (right) at viewing angles of , where parameters and slopes are marked on the panels. The light curves for the bottom two simulations ( and ) are shown in two ways: upper sub-panels with full observed emission variability in black and the long-term trend baseline components underlying the variability in red; the bottom sub-panels show the light curves after subtraction of these best-fit trend lines. In the right column images the (red) circles are binned data points and (green) lines are the best fits yielding the quoted low-frequency slopes; PSDs computed from the post-subtracted trends are shown in the bottom two panels.
0.9 2.29 0.005 2.1 1.8
0.9 2.29 0.01 2.2 1.5
0.9 2.29 0.02 2.2 1.8
0.9* 2.29 0.02 2.4 1.5
0.95 3.20 0.01 2.2 1.4
0.95* 3.20 0.01 2.4 3.3
0.95 3.20 0.02 2.2 1.7
0.97 4.11 0.01 2.2 1.4
0.97* 4.11 0.01 2.5 2.2
0.97 4.11 0.02 2.2 1.9
0.97* 4.11 0.02 2.4 2.8
0.985 5.80 0.02 2.2 1.6
0.985* 5.80 0.02 2.3 1.5
0.99 7.09 0.01 2.1 1.9
0.99* 7.09 0.01 2.3 2.4
  • A linear trend was subtracted from the first portion of the emission light curve before the PSD was computed.

Table 3: PSD results

The right column of Figure 12 shows the corresponding PSDs of these light curves with those for which the trends were subtracted being displayed for the bottom two panels. To compute the PSDs, we adopted Hanning window functions which reduce red-noise leakage in Fourier transform filtering (Harris, 1978; Martinis & Geller, 2014). There are many fewer data points on the low frequency end of the PSDs than on the high end. To mitigate this inequality, when computing the slopes of the PSDs we binned the data points at 0.1 dex frequency intervals and computed an average for each bin, which have been marked with red circles on the images. These averaged points for each bin were used to fit a power-law line to the data from the log–log plot of power against frequency to calculate a final slope up until the white noise contribution becomes important. Although each light curve is unique and has significant differences, we found a rather narrow range in the best fitting slopes of the PSDs. The slopes and reduced values for the straight-line fits to the red-noise portions of the PSDs that we have computed are listed in Table 3 for 9 distinct runs, of which 5 were also calculated after trend subtractions were performed. The best fitting slopes of the resulting PSDs had the range . We note that the subtraction of a linear baseline from parts of the light curve produces a PSD slope that is steeper by 0.1–0.3 compared to the unsubtracted PSDs. Regardless of approach taken, our PSDs are in accord with observations of optical and X-ray variations from AGN, which have been seen to range between and but do seem to cluster between and (MacLeod et al., 2010; González-Martín & Vaughan, 2012; Edelson et al., 2013; Revalski et al., 2014; Kasliwal et al., 2015; Mohan et al., 2016; Smith et al., 2018; Wehrle et al., 2013, 2018). Similar results () were obtained from our earlier and cruder 2D RHD simulations that did include the modest effects of taking the time-travel lags into account (Pollack et al., 2016).

It is certainly gratifying that these simulations produce light curves that are reminiscent of those of many blazars, in terms of having both active and quiescent periods and that the PSD slopes are also in accord with observations; however, we make no claim that this approach provides a unique explanation for those phenomena. In particular, we note that very fast variations cannot plausibly be obtained by the changes in the properties of jet zones as we have modeled them. Activity on scales much smaller than the grid cells in our jet simulations appear to be necessary as well. These would involve turbulence (Marscher, 2014; Calafut & Wiita, 2015; Pollack et al., 2016) or ultrarelativistic mini-jets within the jets, most likely produced by magnetic reconnection (e.g. Giannios et al., 2009; Biteau & Giebels, 2012; Kagan et al., 2016; O’Riordan et al., 2017).

7 Conclusions

We have simulated an exceptionally large suite of over 50 3D RHD propagating jets using the Athena code, with an emphasis on investigating their suitability for modeling observed blazar variability. Our simulations of propagating jets have spanned a significant range of velocities for the initial bulk flows () that cover the great majority of the velocities deduced for radio galaxies and blazars (e.g. Lister et al., 2009). These flows are light, as is appropriate to radio jets, with jet density to ambient medium density ratios between and . Both medium resolution (5 zones per jet radius) and high resolution (10 zones per radius) simulations have been completed extending out to 60, 120, or even 240 jet radii along the direction of motion; in all cases our simulations had widths of 50 jet radii in the two perpendicular directions so there was no loss of matter out of the grid along those transverse boundaries or the need to worry about waves reflecting off those boundaries and unphysically distorting the jet flow.

These simulations span a sufficient range in power so that the weaker ones go unstable before they pass through our simulation volumes. When scaled to extragalactic dimensions and parameters expected for intracluster (or intergalactic) ambient media these cases yield FR I type morphologies; if scaled to dimensions and properties probed by VLBI, such dying jets would represent failed radio sources. The majority of our simulations took more advantage of the relativistic velocities computable with Athena and correspond to powerful sources that remain stable for very extended distances and times. On large scales these would be FR II radio galaxies and on the small scales they would be young radio galaxies, probably seen as compact steep spectrum sources if not aimed close to our line of sight. Although bow shocks continue to be driven by the weaker sources for some time after the jets have clearly gone strongly unstable, the rate of propagation slows down, while our powerful relativistic jets remain extremely well collimated for very extended times and the rate of advance of both the bow shocks and the jet tips are close to constant for them.

The main suite of simulations was performed for jets injected in pressure balance with the ambient medium and where an adiabatic index of 5/3 was taken for both jet and ambient fluids. However, a few overpressured simulations () have been performed to compare with equal pressure simulations. Their morphologies are quite similar, though unsurprisingly we saw that otherwise identical overpressured simulations propagate a little faster than do the standard ones and they have wider jets, at least at the early stages of the simulations. We also performed a few simulations where instead of 5/3 was assumed: the former propagated somewhat faster than the latter and may stay stable for slightly longer times but do not display major differences. Comparison between our standard 3D RHD and 2D (slab-like) RHD simulations were also made and they show that 2D simulations have wider bow shocks and cocoons and take longer to propagate across the entire grid than do 3D ones.

We approximate the synchrotron light curves of these propagating jets by summing up the fluxes emitted from zones with different densities and pressures, where the intrinsic emissivity is taken to be proportional to the product of and . As we deal with propagating jets, these physical quantities vary with time but are higher for more powerful jets. Each of those rest frame “emissivities” is then differentially Doppler boosted by the changing velocities each zone evinces at each time-step downloaded from the simulations. We then sum up these contributions over a substantial portion of the jet close to the first reconfinement shock and then compute observed total brightnesses and their power spectra. These do seem to be able to reproduce key aspects of blazar light curves and yield PSD slopes a little steeper than , which also are found to characterize a significant majority of observations of the optical and X-ray variations from AGN.

The biggest shortcoming of these simulations is that they have been restricted to flows where the magnetic fields have not been incorporated. We are in the process of computing MHD simulations which we will then analyze in a similar fashion to see if those extra degrees of freedom (and complication) yield interestingly different variability characteristics.


We thank the anonymous referee for useful criticisms and suggestions that have led to a substantially improved manuscript. The authors acknowledge use of the ELSA high performance computing cluster at The College of New Jersey for conducting the research reported in this paper and are most grateful to Shawn Sivy for assistance in porting the Athena code to the cluster and for his efforts in constructing and maintaining this facility. This cluster is funded by the National Science Foundation under grant number OAC-1828163. ParaView software was used for the production of most of the figures. YL was supported in part by the China Scholarship Council (file No. 201706220272). This work is partly supported by Natural Science Foundation of China under grant No. 11873035, the Natural Science Foundation of Shandong province (No. JQ201702) and the Young Scholars Program of Shandong University (No. 20820162003).


  • Aller et al. (1985) Aller, H. D., Aller, M. F., & Hughes, P. A. 1985, ApJ, 298, 296
  • Bassett & Woodward (1995) Bassett, G. M., & Woodward, P. R. 1995, ApJ, 441, 582
  • Beckwith & Stone (2011) Beckwith, K., & Stone, J. M. 2011, ApJS, 193, 6
  • Begelman et al. (1984) Begelman, M. C., Blandford, R. D., & Rees, M. J. 1984, RvMP, 56, 255
  • Belan et al. (2010) Belan, M., de Ponte, S., & Tordella, D. 2010, Phys. Rev. E, 82, 026303
  • Bicknell (1994) Bicknell, G. V. 1994, ApJ, 422, 542
  • Bicknell (1995) Bicknell, G. V. 1995, ApJS, 101, 29
  • Biretta et al. (1986) Biretta, J. A., Moore, R. L., & Cohen, M. H. 1986, ApJ, 308, 93
  • Biteau & Giebels (2012) Biteau, J. & Giebels, B. 2012, A&A, 548, A123
  • Blandford & Rees (1974) Blandford, R. D., & Rees, M. J. 1974, MNRAS, 169, 395
  • Bliton et al. (1998) Bliton, M., Rizza, E., Burns, J. O., Owen, F. N., & Ledlow, M. J. 1998, MNRAS, 301, 609
  • Bodo et al. (1995) Bodo, G., Massaglia, S., Rossi, P., et al. 1995, A&A, 303, 281
  • Burns et al. (1991) Burns, J. O., Norman, M. L., & Clarke, D. A. 1991, Sci, 253, 522
  • Calafut & Wiita (2015) Calafut, V., & Wiita, P. J. 2015, JApJ, 36, 255
  • Camenzind & Krockenberger (1992) Camenzind, M., & Krockenberger, M. 1992, A&A, 255, 59
  • Caproni & Abraham (2004) Caproni, A., & Abraham, Z. 2004, ApJ, 602, 625
  • Carilli & Barthel (1996) Carilli, C. L., & Barthel, P. D. 1996, A&ARv, 7, 1
  • Carvalho & O’Dea (2002a) Carvalho, J. C., & O’Dea, C. P. 2002a, ApJS, 141, 337
  • Carvalho & O’Dea (2002b) Carvalho, J. C., & O’Dea, C. P. 2002b, ApJS, 141, 371
  • Cavagnolo et al. (2010) Cavagnolo, K. W., McNamara, B.R., Nulsen, P. E. J., et al. 2010, ApJ, 720, 1066
  • Clarke et al. (1986) Clarke, D. A., Norman, M. L., & Burns, J. O. 1986, ApJ, 311, L63
  • Donohoe & Smith (2016) Donohoe, J., & Smith, M. D. 2016, MNRAS, 458, 558
  • Edelson et al. (2013) Edelson, R., Mushotzky, R., Vaughan, S., et al. 2013, ApJ, 766, 16
  • Fanaroff & Riley (1974) Fanaroff, B. L., & Riley, J. M. 1974, MNRAS, 167, 31P
  • Ferrari (1998) Ferrari, A. 1998, ARA&A, 36, 539
  • Gaibler et al. (2009) Gaibler V., Krause, M., & Camenzind, M. 2009, MNRAS, 400, 1785
  • Gardiner & Stone (2005) Gardiner, T. A., & Stone, J. M. 2005, JCoPh, 205, 509
  • Giannios et al. (2009) Giannios, D., Uzdensky, D. A., & Begelman, M. C. 2009, MNRAS, 395, L29
  • Godfrey & Shabala (2013) Godfrey, L. E. H., & Shabala, S. S. 2013, ApJ, 767, 12
  • González-Martín & Vaughan (2012) González-Martín, O., & Vaughan, S. 2012, A&A, 544, A80
  • Gopal-Krishna & Wiita (1991) Gopal-Krishna, & Wiita, P. J. 1991, ApJ, 373, 325
  • Gopal-Krishna & Wiita (1992) Gopal-Krishna, & Wiita, P. J. 1992, A&A, 259, 109
  • Gopal-Krishna & Wiita (2000) Gopal-Krishna, & Wiita, P. J. 2000, A&A, 363, 507
  • Gopal-Krishna & Wiita (2001) Gopal-Krishna, & Wiita, P. J. 2001, A&A, 373, 100
  • Gopal-Krishna & Wiita (2004) Gopal-Krishna, & Wiita, P. J. 2004, arXiv: astro-ph.0409761
  • Gopal-Krishna et al. (2003) Gopal-Krishna, Stalin, C. S., Sagar, R., & Wiita, P. J. 2003, ApJ, 586, L25
  • Hardcastle & Krause (2014) Hardcastle, M. J., & Krause, M. G. H. 2014, MNRAS, 443, 1482
  • Hardee & Norman (1988) Hardee, P. E., & Norman, M. L. 1988, ApJ, 334, 70
  • Hardee & Norman (1990) Hardee, P. E., & Norman, M. L. 1990, ApJ, 365, 134
  • Hardee et al. (2001) Hardee, P. E., Hughes, P. A., Rosen, A., & Gomez, E. A. 2001, ApJ, 555, 744
  • Harris (1978) Harris, F. J. 1978, Proceedings of the IEEE, 66, 51
  • Hooda et al. (1994) Hooda, J. S., Mangalam, A. V., & Wiita, P. J. 1994, ApJ, 423, 116
  • Hooda & Wiita (1998) Hooda, J. S., & Wiita, P. J. 1998, ApJ, 493, 81
  • Horne & Baliunas (1986) Horne, J. H., & Baliunas, S. L., 1986, ApJ, 302, 757
  • Huarte-Espinosa et al. (2011) Huarte-Espinosa, M., Krause, M., & Alexander, P. 2011, MNRAS, 417, 382
  • Hughes et al. (1991) Hughes, P. A., Aller, H. D., & Aller, M. F. 1991, ApJ, 374, 57
  • Jeyakumar et al. (2005) Jeyakumar, S., Wiita, P. J., Saikia, D. J., & Hooda, J. S. 2005, A&A, 432, 823
  • Jiang et al. (2007) Jiang, L., Fan, X., Ivezić, Z., et al. 2007, ApJ, 656, 680
  • Kagan et al. (2016) Kagan, D., Nakar, E., & Piran, T. 2016, ApJ, 826, 221
  • Kapińska et al. (2017) Kapińska, A. D., Terentev, I., Wong, O. I. et al. 2017, AJ, 154, 253
  • Kasliwal et al. (2015) Kasliwal, V. P., Vogeley, M. S., & Richards, G. T. 2015, MNRAS, 451, 4328
  • Kawakatu et al. (2009) Kawakatu, N., Kino, M., & Nagai, H. 2009, ApJ, 697, L173
  • Keppens et al. (2008) Keppens, R., Meliani, Z., van der Holst, B., & Casse, F. 2008, A&A, 486, 663
  • Krause & Alexander (2007) Krause, M., & Alexander, P. 2007, MNRAS, 376, 465
  • Ledlow & Owen (1996) Ledlow, M. J., & Owen, F. N. 1996, AJ, 112, 9
  • Leismann et al. (2005) Leismann, T., Antón, L., Aloy, M. A., et al. 2005, A&A, 436, 503
  • Lister et al. (2001) Lister, M. L., Tingay, S. J., & Preston, R. A. 2001, ApJ, 554, 964
  • Lister et al. (2009) Lister, M. L., Cohen, M. H., Homan, D. C., et al. 2009, AJ, 138, 1874
  • Lister et al. (2013) Lister, M. L., Aller, M. F., Aller, H. D., et al. 2013, AJ, 146, 120
  • Loken et al. (1995) Loken, C., Roettiger, K., Burns, J. O., & Norman, M. 1995, ApJ, 445, 80
  • MacLeod et al. (2010) MacLeod, C. L., Ivecić, Z., Kochaneck, C. S., et al. 2010, ApJ, 721, 1014
  • Marscher & Gear (1985) Marscher, A. P., & Gear, W. K. 1985, ApJ, 298, 114
  • Marscher et al. (2008) Marscher, A. P., Jorstad, S. G., D’Arcangelo, F. D., et al. 2008, Natur, 452, 966
  • Marscher et al. (2010) Marscher, A. P., Jorstad, S. G., Larionov, V. M., et al. 2010, ApJL, 710, L126
  • Marscher (2014) Marscher, A. P. 2014, ApJ, 780, 87
  • Martí (2015) Martí, J. M. 2015, MNRAS, 452, 3106
  • Martí et al. (1995) Martí, J. M., Müller, E., Font, J. A., et al. 1995, ApJL, 448, L105
  • Martí et al. (1997) Martí, J. M., Müller, E., Font, J. A., et al. 1997, ApJ, 479, 151
  • Martí & Müller (2003) Martí, J. M., & Müller, E. 2003, LRR, 6, 7
  • Martí et al. (2016) Martí, J. M., Perucho, M., & Gómez, J. L., 2016, ApJ, 831, 163
  • Martinis & Geller (2014) Martinis, J. M. & Geller, M. R., 2014, Phys. Rev. A90, 022307
  • Massaglia (2003) Massaglia, S. 2003, Ap&SS, 287, 223
  • Massaglia et al. (2016) Massaglia, S., Bodo, G., Rossi, P., Capetti, S., & Mignone, A. 2016, A&A, 596, 12
  • Mignone et al. (2005) Mignone, A., Massaglia, S., & Bodo, G. 2005, ApJS, 160, 199
  • Mignone & McKinney (2007) Mignone, A., & McKinney, J. C. 2007, MNRAS, 378, 1118
  • Mignone et al. (2010) Mignone, A., Rossi, P., Bodo, G., Ferrari, A., & Massaglia, S. 2010, MNRAS, 402, 7
  • Mioduszewski et al. (1997) Mioduszewski, A. J., Hughes, P. A., & Duncan, C. G. 1997, ApJ, 476, 649
  • Mohan et al. (2016) Mohan, P., Gupta, A. C., Bachev, R., & Strigachev, A. 2016, MNRAS, 456, 654
  • Nishikawa et al. (1997) Nishikawa, K.-I., Koide, S., Sakai, J.-I, Christodoulou, D. M., Sol, H., & Mutel, R. L. 1997, ApJ, 484, 794
  • Norman et al. (1982) Norman, M. L., Winkler, K.-H., Smarr, L., & Smith, M. D. 1982, A&A, 113, 285
  • Norman (1996) Norman, M. L. 1996, in ASP Conf. Ser. Vol. 100, Energy Transport in Radio Galaxies and Quasars, ed. P. E. Hardee, A. H. Bridle & J. A. Zensus (San Francisco: ASP), p. 405
  • Owen & Ledlow (1994) Owen, F. N., & Ledlow, M. J. 1994 in ASP Conf. Proc. 54, The First Stromlo Symposium: The Physics of Active Galaxies, ed. G. V. Bicknell, P. J. Quinn, & M. A. Dopita (San Francisco: ASP), p. 319
  • O’Neill et al. (2005) O’Neill, S. M., Tregillis, I. L., Jones, T. W., & Ryu, D. 2005, ApJ, 633, 717
  • O’Riordan et al. (2017) O’Riordan, M., Pe’er, A., McKinney, J. C. 2017, ApJ, 843, 81
  • Peterson (1997) Peterson, B. M. 1997, An Introduction to Active Galactic Nuclei, (Cambridge: Cambridge U. Press)
  • Piner (2003) Piner, B. G., Unwin, S. C., Wehrle, A. C., et al. 2003, ApJ, 588, 716
  • Press & Rybicki (1989) Press, W. H., & Rybicki, G. B. 1989, ApJ, 338, 277
  • Pollack et al. (2016) Pollack, M., Pauls, D., & Wiita, P. J. 2016, ApJ, 820, 12
  • Rayburn (1977) Rayburn, D. R. 1977, MNRAS, 179, 603
  • Revalski et al. (2014) Revalski, M., Nowak, D., Wiita, P. J., Wehrle, A. E., & Unwin, S. C. 2014, ApJ, 785, 60
  • Ros et al. (2000) Ros, E., Zensus, J. A., & Lobanov, A. P. 2000, A&A, 354, 55
  • Rosen et al. (1999) Rosen, A., Hughes, P. A., Duncan, G. C., & Hardee, P. E. 1999, ApJ, 516, 729
  • Ryu et al. (2006) Ryu, D., Chattopadhyay I., & Choi, E. 2006, ApJS, 166, 410
  • Sarazin (1988) Sarazin, C. L. 1988, X-Ray Emission from Clusters of Galaxies (Cambridge: Cambridge University Press)
  • Scheck et al. (2002) Scheck, L., Aloy, A. A., Martí, J. M., Gómez, J. L., & Müller, E. 2002, MNRAS, 331, 615
  • Scheuer (1974) Scheuer, P. A. G. 1974, MNRAS, 166, 513
  • Smith et al. (2018) Smith, K. L., Mushotzky, R. F., Boyd, P. T., Malkan, M., Howell, S. B., & Gelino, D. M. 2018, ApJ, 857, 141
  • Stone et al. (2008) Stone, J. M., Gardiner, T. A., Teuben, P., Hawley, J. F., & Simon, J. B. 2008, ApJS, 178, 137
  • Synge (1957) Synge, J. L. The Relativistic Gas (Amsterdam: North Holland Publ.)
  • Urry & Padovani (1995) Urry, C. M., & Padovani, P. 1995, PASP, 107, 803
  • Uttley & McHardy (2005) Uttley, P., & McHardy, I. M. 2005, MNRAS, 363, 586
  • van Putten (1996) van Putten, M. H. P. M. 1996, ApJL, 467, L57
  • Wehrle et al. (2013) Wehrle, A. E., Wiita, P. J., Unwin, S. C., et al. 2013, ApJ, 773, 89
  • Wehrle et al. (2018) Wehrle, A. E., Carini, M. C., & Wiita, P. J. 2018, submitted to ApJ
  • Wiita (1978) Wiita, P. J. 1978, ApJ, 221, 436
  • Wold et al. (2007) Wold, M., Lacy, M., & Armus, L. 2007, A&A, 470, 531
  • Zrake & MacFadyen (2013) Zrake, J., & MacFadyen, A. I. 2013, ApJ, 763, L12
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description