Simulating the synchrotron emission of AGN with grid based relativistic hydrodynamics

Grid based numerical hydrodynamic codes have been used to simulate the evolution of a variety of astrophysical environments, ranging from stellar winds to accretion disks in binary systems. These simulations can be used as valuable tools to correlate the theoretical models and observational data. In this study a grid based relativistic hydrodynamic simulation is created, using the freely available numerical code PLUTO ver 4.2, to study possible causes of variability in AGN. The code uses shock-capturing Godunov-type methods to solve the fluid dynamic conservation equations on a structured mesh grid. To simulate the production of a jet similar to those observed in AGN a three dimensional static grid was created containing 256 × 256 × 512 computational cells. A uniform background medium was assigned to this grid, while a nozzle was defined on the z = 0 boundary to inject a pressure matched relativistic jet, with Γ = 10 (v = 0:995c), into the medium. The simulation showed the formation of a central collimated relativistic beam along with a surrounding cocoon region. In addition, a post-processing emission modelling code is being developed, in order to determine the emission produced by the relativistic jet. The emission model utilizes the data produced by the hydrodynamic code to calculate intensity maps due to synchrotron self absorption. The emission code takes into account the effects produced by the relativistic motion of the emitting regions (excluding light travel time) which can lead to effects such as the Doppler boosting observed in blazars. Synchrotron based emission models can be applied to other astrophysical objects, such as X-ray binary systems, which means that the application of the emission code can be adapted to model a variety of astrophysical sources.


Introduction
Relativistic astrophysical jet-like outflows have been associated with many different sources such as Active Galactic Nuclei (AGN), microquasars and gamma-ray bursts (GRB). In the case of AGN the jets can stretch over Mpc scales, carrying large amounts of energy away from the central engine and producing radiation over a wide range of the electro-magnetic spectrum. The spectral energy distribution of AGN are characterized by a double hump structure consisting of a low energy and a high energy component. The emission in the low energy component is dominated by synchrotron emission from non-thermal electrons, while the high energy component has been modelled by a combination of inverse-Compton and hadronic processes, see for e.g. [1].
Many radio surveys have been undertaken to study the synchrotron emission produced in the relativistic jets and how this relates to their dynamic morphology. For example, the MOJAVE program has studied the movement of superluminal components in AGN as well as the change in ejection direction and precession, over year-long time periods [ There are, however, still many unanswered questions about the processes that occur in jets and how they relate to the received emission. Numerical fluid dynamics is a powerful tool to study these dynamic structures. The macroscopic motion of relativistic jets can be described as a fluid and it must, therefore, adhere to the set of partial differential equations which describes the laws of conservation of mass, energy and momentum. By solving the conservation equations numerically, and evolving them with time, we can evolve the jet morphology to future states. This is a good method for investigating the variable nature of relativistic jets.
In order to compare the evolution of structures that are obtained using numerical simulations to the observed variability, we need to determine the synchrotron emission from the physical characteristics of the jet, such as internal energy, magnetic field strength, pressure and density.
In this proceedings we present a numerical simulation of an ideal relativistic jet. In this model the magnetic field was assumed to be dynamically negligible compared to the kinetic energy of the particles. A detailed discussion of the numerical simulation is given in section 2 while the emission modelling follows in section 3. Section 4 contains our results followed by a conclusion in section 5.

Numerical models
In hydrodynamical simulations a numerical jet environment is constructed and evolved with time according to the conservation equations of fluid dynamics which can be expressed as where U is a column vector consisting of conserved quantities, T ( U ) is a tensor containing the flux vectors as a function of U , and S( U ) is a tensor containing source terms that can be used to introduce effects such as viscosity and gravitational forces [3]. Since AGN jets are have a high bulk Lorentz factor it is necessary to perform the simulations in the relativistic hydrodynamic (RHD) regime. In this case, the constituents of equation (1) have the following form: where ρ is the proper density, P is the pressure, h is the specific enthalpy, Γ is the Lorentz factor, I is a 3x3 unit matrix and v is the three velocity [3]. All other fluid properties, such as the pressure and the internal energy, can be determined from the above mentioned variables [4]. For example, the pressure can be solved by combining the ideal caloric equation of state, where e is the internal energy density of the fluid, with the definition of the sound speed The complex nature of this set of equations (equation 1) makes it impossible to find analytical solutions for these types of environments. We, therefore, need to numerically solve these equations. For this study the numerical hydrodynamical code PLUTO was used to evolve the numerical simulation. PLUTO  uses high resolution shock capturing algorithms to solve the conservation equations and evolve them with time [5]. For our simulation we set up a three dimensional Cartesian mesh grid consisting of 256 × 256 × 512 computational cells. The grid was initially assigned a uniform background medium at rest. A nozzle with a radius of 4 computational cells, was created on the z = 0 axis at the bottom of the simulation and jet material was injected into the environment at a constant rate. For this simulation the density of the injected jet was chosen to be less than the ambient density (η = ρ b ρam = 10 −6 ). The simulation was performed using arbitrary units in order to avoid numerical errors. A pressure matched model was used for the simulation in order to help collimate the jet. The parameters used for the simulation were based on previous papers such as [4,6]. All parameters used in the simulation are listed in table 1.
The simulation was run using the HLLC Riemann solver, piecewise parabolic interpolation and characteristic tracing time stepping [7]. The environment was evolved until the bow shock created by the injection of jet material reached the edge of the computational domain.

Synchrotron emission
In order to compare the numerical simulations with the observed emission in radio frequencies we need to translate the three dimensional environment of physical parameters onto a two dimensional image plane. In order to do this we first calculate the synchrotron emission and absorption coefficients for each cell in the simulation, using a post-processing code. In this code we assumed that the electrons inside the jet consist of particles with a power-law energy distribution with the number density, n(γ), given by, where the normalization factor, n 0 , is determined by [8] n 0 = e(p − 2) Here m p is the mass of a proton, p is the power-law index and C E is the ratio of the maximum and minimum energies. For our initial estimates we chose p = 1.8 and C E = 10 3 . The magnetic field energy density is assumed to be an equipartition fraction of the internal energy density, where B = 10 −3 based on values fitted by [9]. An analytical delta-approximation model was used to estimate the synchrotron emission (j sy v ) and absorption (α sy v ) coefficients for each cell [10], given by α sy ν = 2 9 p + 2 mν 2 where, represents the electron gyro-frequency. Here q is the charge and m is the mass of the nonthermal electrons in the jet, c is the speed of light, u B is the magnetic field energy density given by equation 6 and ν is the frequency of emission in the co-moving frame [1].
Since the jet material is moving at relativistic velocities we need to apply Lorentz transformations to each cell. The synchrotron coefficients of each cell transform as, where the primed quantities are in the observer frame and the unprimed quantities are in the co-moving frame, β is the magnitude of the velocity in units of c and µ is the cosine of the angle between the observer and the velocity of the fluid in the observer frame. For our initial calculations additional effects such as the light travel time and the expansion of the universe were neglected. The post-processing code finally integrates the coefficients along a line of sight, s, defined by an azimuth angle φ and a polar angle θ , to create a 2D intensity map. This is illustrated in figure 1.

Results
The xz-slices of the density, pressure and velocity distributions for the simulation are shown in figure 1. From these figures we can note five distinct regions in the simulation. First we have the outer uniform ambient medium. A bow wave propagates through this medium due to the injection of the supersonic jet material through the nozzle on the initial z boundary. The bow wave increases the pressure of the medium as it propagates creating a high pressure region between the jet and the ambient medium. The centre of the jet contains collimated jet material moving at a high Lorentz factor through the environment. Surrounding this collimated beam we find a turbulent backflow of material which is generated as the jet material interacts with the shock front, called the working surface, which forms between the jet material and the ambient medium. Figure 2 shows 1 dimensional plots of the pressure and density in the z direction along the centre of the beam. We note initial oscillation in both pressure and density. The oscillations indicate the formation of stationary re-collimation shocks. The shocks are caused by the overpressured cocoon acting on the central beam. These shocks are damped as the jet material propagates along the z-direction. At larger distances from the injection site turbulence causes the beam of the jet to become unstable. Figure 2 shows this as a rapid change in the density and pressure of the beam. Finally we see a discontinuity in pressure and density at the working surface where the jet material collides with the ambient medium. at 15 GHz illustrating the time dependence of the main emission regions, or hot spots, in the emission. The hot spot at the head of the jet is shown to change its location, size and maximum intensity. As time progresses we note that the hotspot breaks apart to form several extended emission regions.

Conclusion
In this paper a numerical model of a relativistic jet was set up and evolved with time using the PLUTO hydrodynamical code. The results showed the formation of a relativistic beam surrounded by a cocoon of backflowing material. Plots of the pressure and density along the central beam revealed the presence of periodic re-collimation shocks within 40 jet radii. Beyond this region the beam becomes turbulent with rapid variations in pressure and density, ending in a shock front at the head of the jet. A delta-approximation synchrotron model was applied to the simulation to produce a 2D intensity map of an edge-on system. The resulting intensity maps showed time dependent emission regions varying in intensity, position and size. This result indicates that variations in the emission of AGN jets can occur even if the system is fed by a constant injection.
In order to make a precise comparison between the synthetic intensity maps and real observation, however, additional corrections to emission model must be applied for effects such as light travel time and expansion of the universe and radiative losses. The future aim for this project is to include an inverse Compton model into the emission code to investigate the correlation of variability at multi-wavelengths. Based on this we will attempt SED modelling of the simulations in order to produce a better qualitative comparison with observational data. We also plan to include a variable injection rate in the simulation to model how variable injection influences the morphology and stability of the resulting jet. Finally the project can be extended to a much larger scope, including many different astronomical sources such as X-ray binaries and gamma-ray bursts.