# Simulations of Relativistic Jets with GENESIS

## 1 Introduction

Astrophysical jets are continuous channels of plasma produced by some active galactic nuclei that are currently observed in radio frequences. The relativistic nature of the plasma has been inferred from (esentially) two observational evidences: (i) the existence of superluminal motions of some radio components and, (ii) the high flux variability (even smaller than one day for some sources). Since several years the dynamical and morphological properties of axisymmetric relativistic jets are investigated by means of relativistic hydrodynamic simulations (e.g.,). In addition, relativistic MHD simulations have been performed in 2D () and 3D (). In their 3D simulations,  and , have studied mildly relativistic jets (Lorentz factor, $`W=4.56`$) propagating both along and obliquely to an ambient magnetic field.

In this work we report on high-resolution 3D simulations of relativistic jets with the largest beam flow Lorentz factor performed up to now (7.09), the largest resolution (8 cells per beam radius), and covering the longest time evolution (75 normalized time units; a normalized time unit is defined as the time needed for the jet to cross a unit length). These facts together with the high performance of our hydrodynamic code allowed us to study the morphology and dynamics of 3D relativistic jets for the first time.

The calculations have been performed with the high–resolution 3D relativistic hydrodynamics code GENESIS , which is an upgraded version of the code developed by . GENESIS integrates the 3D relativistic hydrodynamic equations in conservation form in Cartesian coordinates including an additional conservation equation for the beam-to-external density fraction to distinguish between beam and external medium fluids. The code is based on a method of lines which first discretizes the spatial part of the relativistic Euler equations and solves the fluxes using the Marquina’s flux formula . Then the semidiscrete system of ordinary differential equations is solved using a third order Runge-Kutta algorithm . High spatial accuracy is achieved by means of a PPM third order interpolation . The computations were performed on a Cartesian domain (X,Y,Z) of size $`15R_b\times 15R_b\times 75R_b`$ ($`120\times 120\times 600`$ computational cells), where $`R_b`$ is the beam radius. The jet is injected at $`z=0`$ in the direction of the positive $`z`$-axis through a circular nozzle defined by $`x^2+y^2R_b^2`$. Beam material is injected with a beam mass fraction $`f=1`$, and the computational domain is initially filled with an external medium ($`f=0`$).

We have considered a 3D model corresponding to model C2 of , which is characterized by a beam-to-external proper rest-mass density ratio $`\eta =0.01`$, a beam Mach number $`M_b=6.0`$, and a beam flow speed $`v_b=0.99c`$ ($`c`$ is the speed of light) or a beam Lorentz factor $`W_b7.09`$. An ideal gas equation of state with an adiabatic exponent $`\gamma =5/3`$ describes both the jet matter and the ambient gas. The beam is assumed to be in pressure equilibrium with the ambient medium.

The evolution of the jet was simulated up to $`T150R_b/c`$, when the head of the jet is about to leave the grid. The scaled final time $`T\mathrm{4.6\hspace{0.17em}10}^4\left(R_b/100\mathrm{pc}\right)\mathrm{yr}`$ is about two orders of magnitude smaller than the estimated ages of powerful jets. Hence, our simulations cannot describe the long term evolution of these sources.

Non–axisymmetry was imposed by means of a helical velocity perturbation at the nozzle given by

$$v_b^x=\zeta v_b\mathrm{cos}\left(\frac{2\pi t}{\tau }\right),v_b^y=\zeta v_b\mathrm{sin}\left(\frac{2\pi t}{\tau }\right),v_b^z=v_b\sqrt{1\zeta ^2},$$
(1)   
where $`\zeta `$ is the ratio of the toroidal to total velocity and $`\tau `$ the perturbation period (i.e.,$`\tau =T/n`$, $`n`$ being the number of cycles completed during the whole simulation). This velocity field causes a differential rotation of the beam. The perturbation is chosen such that it does not change the velocity modulus, (i.e., mass, momentum and energy fluxes of the beam are preserved).

## 2 Morphodynamics of 3D relativistic jets

Here we consider two models: A, which has a $`1\%`$ perturbation in helical velocity ($`\zeta =0.01`$) and $`n=50`$ and B, with $`\zeta =0.05`$ and $`n=15`$. Figure 1 shows various quantities of the model A in the plane $`y=0`$ at the end of the simulation. Two values of the beam mass fraction are marked by white contour levels. The beam structure is dominated by the imposed helical pattern with amplitudes of $`0.2R_b`$ and $`1.2R_b`$ for A and B, respectively.

The overall jet’s morphology is characterized by the presence of a highly turbulent, subsonic cocoon. The pressure distribution outside the beam is nearly homogeneous giving rise to a symmetric bow shock (Fig. 1b) in model A. Model B shows a very inhomogeneous pressure distribution in the cocoon. As in the classical case , the relativistic 3D simulation shows less ordered structures in the cocoon than the axisymmetric models. As seen from the beam mass fraction levels, the cocoon remains quite thin ($`2R_b`$) in A and widens ($`4R_b`$) in B.

The flow field outside the beam shows that the high velocity backflow is restricted to a small region in the vicinity of the hot spot (Fig. 1e), the largest backflow velocities ($`0.5c`$) being significantly smaller than in 2D models. The flow annulus with high Lorentz factor found in axisymmetric simulations is also present, but it is reduced to a thin layer around the beam and possesses sub-relativistic speeds ($`0.25c`$) in model A and mildly relativistic ($`0.7`$) in B. The size of the backflow velocities in the cocoon do not support relativistic beaming in case of small perturbations but such possibility is open in larger ones.

Within the beam the perturbation pattern is superimposed to the conical shocks at about 26 and 50 $`R_b`$. The beam of A does not exhibit the strong perturbations (deflection, twisting, flattening or even filamentation) found by other authors  for 3D classical hydrodynamic jets;  for 3D classical MHD jets). This can be taken as a sign of stability, although it can be argued that our simulation is not evolved far enough. For n15p01, the beam is about to be disrupted at the end of our simulation. Obviously, the beam cross section and the internal conical shock structure are correlated (Figure 1).

The helical pattern propagates along the jet at nearly the beam speed which could yield to superluminal components when viewed at appropriate angles. Besides this superluminal pattern, the presence of emitting fluid elements moving at different velocities and orientations could lead to local variations of the apparent superluminal motion within the jet. This is shown in Fig. 1f, where we have computed the mean (along each line of sight, and for a viewing angle of 40 degrees) local apparent speed. The distribution of apparent motions is inhomogeneous and resembles that of the observed individual features within knots in M87 .

The jet can be traced continuously up to the hot spot which propagates as a strong shock through the ambient medium. Beam material impinges on the hot spot at high Lorentz factors in A case, but the beam Lorentz factor strongly decreases for B. We could not identify a terminal Mach disk in the flow. We find flow speeds near (and in) the hot spot much larger than those inferred from the one dimensional estimate. This fact was already noticed for 2D models by  and suggested by them as a plausible explanation for an excess in hot spot beaming.

We find a layer of high specific internal energy (typically more than a tenfold larger than that of the gas in the beam core, see Fig. 1d) surrounding the beam like in previous axisymmetric models . The region filled by the shear layer is defined by $`0.2<f<0.95`$. It is mainly composed of forward moving beam material at a speed smaller than the beam speed (Fig. 1e). The intermediate speed of the layer material is due to shear in the beam/cocoon interface, which is also responsible for its high specific internal energy. The shear layer broadens with distance from 0.2$`R_b`$ near the nozzle to 1.1$`R_b`$ (in A) or 2.0$`R_b`$ (in B) near the head of the jet. The diffusion of vorticity caused by numerical viscosity is responsible for the formation of the boundary layer. Although being caused by numerical effects (not by the physical mechanism of turbulent shear) the properties of PPM–based difference schemes are such that they can mimic turbulent flow to a certain degree . The existence of such a boundary layer has been invoked by several authors  to interpret a number of observational trends in FRI radio sources. Such a layer will produce a top-bottom asymmetry due to the light aberration , and additionally, it can be used to explain the rails in polarization found by . Other authors  have found evidence for these boundary layers in FRIIs (3C353) radio sources.

The jet’s propagation proceeds in two distinct phases. First it propagates according to a linear of 1D phase, and then the behavior depends on the strength of the perturbation: it accelerates to a propagation speed which is $`20`$% larger than the corresponding 1D estimate in model A or it deccelerates up to $`0.37c`$. The second result partially agrees with the one obtained by .

The axial component of the momentum of the beam particles (integrated across the beam) along the axis decreases by more than a 30% within the first 60 $`R_b`$. Neglecting pressure and viscous effects, and assuming stationarity the axial momentum should be conserved, and hence the beam flow is decelerating. The momentum loss goes along with the growth of the boundary layer.

In model A, although the beam material decelerates, its terminal Lorentz factor is still large enough to produce a fast jet propagation. On the other hand, in 3D, the beam is prone to strong perturbations which can affect the jet’s head structure. In particular, a simple structure like a terminal Mach shock will probably not survive when significant 3D effects develop. It will be substituted by more complex structures in that case, e.g., by a Mach shock which is no longer normal to the beam flow and which wobbles around the instantaneous flow direction. Another possibility is the generation of oblique shocks near the jet head due to off–axis oscillations of the beam. Both possibilities will cause a less efficient deceleration of the beam flow at least during some epochs. At longer time scales the growth of 3D perturbations will cause the beam to spread its momentum over a much larger area than that it had initially, which will efficiently reduce the jet advance speed.