Hydrodynamical Simulations of Recollimation Shocks within Relativistic Astrophysical Jets

Astrophysical jets launched from active galactic nuclei can remain tightly collimated over large distances due, in part, to recollimation shocks. Formed within the jets due to their supersonic nature, recollimation shocks are predicted to leave signatures in the observed radio emission due to magnetic flux freezing and the geometric relationship between magnetic fields and the polarization of synchrotron radiation. In the course of this work, we will compare how predictions of emission from recollimation shocks change when the flow is modelled using a hydrodynamical code, as opposed to semi-dynamical and magnetohydrodynamical codes. Jets generally exhibit low levels of polarization, which implies a substantially disordered magnetic field. It is difficult to model such fields using magnetohydrodynamics, hence this work uses hydrodynamical code and a statistical treatment of the magnetic field (c.f. Scheuer and Matthews, 1990). It should then be possible to assess whether certain radio jet phenomena, such as knots and radio-cores, may be modelled as singular or multiple recollimation shocks. To date, the hydrodynamical code has been successfully built and executed on UCLan’s supercomputer cluster, and parallelepiped vector triads have been included to monitor the fluid deformation within the simulation, so that the emergent flux and polarization may be calculated. The parallelepiped advection is currently being verified and some results are discussed. Code for radiative transfer throughout the jet is also being implemented, in order to simulate images for comparison with previous works and observations.


Introduction
Astrophysical jets are produced by a variety of sources, but the biggest and most luminous examples originate from active galactic nuclei (AGNs). These jets have relativistic bulk motion (with Lorentz factors γ ∼ 10), and are likely composed of an electron-positron dominant plasma [1].
It has long been recognised that jets must interact with their surroundings in order to remain tightly collimated over distances that span up to the megaparsec scale [2]. As a jet propagates away from the AGN, it must displace the denser intergalactic medium (IGM) surrounding the galaxy. The initially over-pressured astrophysical jets will expand laterally as they propagate away from the AGN; because the flows are supersonic, the jet will over-expand and become under-pressured in comparison with the IGM, and thus be driven radially inwards. A series of recollimation shocks will be set up within the jet as it oscillates about an equilibrium pressure.
Radio emission from these jets results almost exclusively from synchrotron radiation, the polarization of which allows us to infer some characteristics of the magnetic field across the jet.
It is often assumed that magnetic flux-freezing tethers the jet plasma and magnetic field together; hence, a transverse shock would suppress the longitudinal field component and order the magnetic field perpendicularly to the flow velocity. Conversely, a shearing flow would enhance the longitudinal magnetic field component, giving rise to a field ordered parallel to the jet; this is generally correct [3], and Kahn [4] produced a fairly successful model with flux-freezing and shearing. It has been shown that moderate compression of a randomly tangled magnetic field can drastically increase the degree of order of the field, as seen from angles close to the plane of compression e.g., [5][6][7]. This compression also increases the plasma emissivity [8] (Section 3), meaning that the observed synchrotron radiation has a higher intensity and degree of polarization e.g., [9,10].
Although other phenomena may have greater influence on the overall jet structure e.g., [11], recollimation shocks should obviously leave polarization signatures in the emission of jets [12], so it may be possible to use recollimation shocks as a diagnostic in determining the some properties of their host AGN [13,14].

Fluid Dynamics and Simulations
Euler's equations ensure conservation of mass, momentum and energy for inviscid fluids [15], and are written as where the number density N, momentum density M, and energy density E are the conservation properties in the observer's frame of reference, pressure P is Lorentz invariant, and the velocity v is relative to the grid. The conserved variables relate to the rest frame values of number density n and energy density e by the following Lorentz transforms: for Lorentz factor γ = (1 − v · v/c 2 ) −1/2 . The Euler equations are presented here in conservation form: ∂ t U + ∇ · F = J, where the solution vector U contains the conservation quantities, F is the flux vector and J is the source term. This is vital for codes utilising shock-capturing methods (including the code used for this work); non-conservative equations may result in misplaced shocks and unstable solutions, depending on the simulation type [16] (Chapter 2.10).

Semi-Dynamical Simulations
Semi-dynamical models represent the action of a shock through jump conditions in the parameters, which gives a sharp edge to the shock. The shocked plasma is also assumed to cool rapidly and only emit near the shock, thus ignoring emission contributions from upstream and inter-shock regions.
Early examples of this approach modelled jet knots with planar propagating shock waves, to some degree of success [6,17]. Fixed, conical recollimation shocks with a randomly tangled magnetic field upstream were simulated [12], and Cawthorne [18] built upon by this model by adding a parallel upstream magnetic field component, and convolving the results with a circular beam, for comparison with observations. These models strongly suggest that knots in radio jets may result from conical shocks; later simulations of stationary shocks in jets 3C 120 [19] and 1803+784 [20] replicated characteristics observed in these sources.
Instead of treating the magnetic field disorder as unresolved on the scale of the simulation e.g., [6], Marscher [21] used a conical recollimation shock and a grid of cells, each with a randomly-oriented uniform magnetic field. The observations were well replicated, except that the model underestimated the γ-ray to X-ray luminosity ratio.

Numerical Simulations
Shocks may be more accurately represented by numerical simulations because their shapes, positions and strengths are calculated, instead of assumed. However, the use of ordered magnetic fields in many numerical simulations means that plasma compression due to a relatively weak shock, such as a recollimation shock, will have little effect on the field polarization.
An early, non-relativistic numerical simulation of the jet in M87 predicted the spacing of the knots in the jet fairly accurately. Later, a more sophisticated numerical simulation was able to passively advect a randomly tangled magnetic field with the flow of the jet, and determine the spectral ageing of the plasma's relativistic electron population [22]. In the model, the jet propagates into an ambient gas of equal pressure, with marker particles advected from the injection region. Each marker particle has an associated parallelepiped vector triad that is initially an orthogonal set of unit vectors. The distortion of the magnetic field at the particle's position is monitored through parallelepiped deformation. The dynamics of the flow were simulated first, and then the synchrotron emission was determined using the parallelepipeds. This approach gave results that quantitatively agreed with the observed polarization and brightness distributions, but strongly over-predicted fractional polarization.
Many more numerical recollimation shock and/or jet simulations have been produced, with notable contributions from the group including J. L. Gómez and J. M. Martí e.g., [23][24][25]. Magnetohydrodynamical (MHD) simulations have also been used to produce jet simulations; these treat the magnetic field as dynamically important, e.g., [26][27][28]. However, the MHD approach requires extremely high resolutions to represent the disordered field that observations lead us to expect [29], resulting in a prohibitively large simulational run-time if MHD code were to be used in this work.

Methodology
A three-dimensional, numerical, parallelized Godunov code for relativistic hydrodynamics [30,31] is being used to produce simulations of astrophysical jets. The code uses RHLLE (an approximate Riemann solver [32,33]) to implement the Godunov approach in the simulations, by evolving the initial discontinuities at cell-boundaries with time.
It was determined that the code could be modified in such a way that the parallelepiped vector triads advect through the simulation in a manner to the conserved variables. By doing this, it can be seen whether randomly tangled magnetic fields advected through recollimation shocks can produce the features observed in knots and radio-cores. The parallelepipeds are initially an orthogonal set of unit vectors; as the fluid experiences shear, compression or expansion, the parallelepiped representing the enclosed volume should be appropriately deformed in shape, orientation and/or size. Parallelepiped deformation may then be used to determine the magnetic field configuration throughout the flow.
This adaptation of the method used by Matthews and Scheuer [22] (discussed in Section 2.2) has the advantage of using a fixed grid of markers that remain evenly distributed throughout the simulation, regardless of fluid flow. The other major endeavour is to evaluate the radiative transfer in a separate code. This will involve implementing the findings of Cawthorne and Hughes [7] in order to calculate the observed emission.
Due to complications detailed in Section 4, the method of determining parallelepiped deformation is currently under scrutiny. When the necessary modifications have been implemented, the simulations will be compared to observations, which are expected to show corresponding features in polarization and emission characteristics.

Results
In order to verify that the parallelepipeds accurately monitor their volume elements, a simulation was produced in which a relativistic, planar shock impinges on ambient medium of equal density. As the density of the jet and ambient medium are initially equal, it may be assumed that the mass density at any point in the simulation is the inverse of the parallelepiped volume V = (c · (a × b)), where a, b and c are the vectors that monitor the parallelepiped's shape. It was found that the parallelepipeds do monitor the compression of the plasma quite well, as can be seen in Figure 1. Clearly, as the inverse of both the parallelepiped volume and the shock-normal vector magnitude give reasonable estimates of the density, we are dealing with a special case. When dealing with recollimation shocks, we must expect shocks that are not normal to the local velocity field. Currently, work is being done to confirm that the parallelepipeds are correctly advecting in three-dimensional simulations, and radiative transfer code is in the process of being written.
For a more complex, three-dimensional case, a jet was initialized as a cylinder of uniform radius and velocity v 0 = 0.9798c, embedded within the stationary ambient medium, which is overdense by a factor of 10, and underpressured by a factor of 0.33, with respect to the initialized jet. These values are fixed for the inlet plane (z ≤ 2 cells), while the rest of the simulation (3 ≤ z ≤ 300) is able to evolve with time. It was found that if the parallelepipeds were prevented from evolving for an extended region of z ≤ 10 cells (effectively having separate boundary conditions for the hydrodynamic variables and the parallelepipeds), the amount of shear they demonstrated was greatly reduced (see Figure 2). This stopped the large inlet plane discontinuity between the jet and ambient medium affecting the parallelepipeds, and the effect was exacerbated by the shear layer. However, it is also clear from Figure 2 that a steady growth of the parallelepiped component parallel to the jet is still present, to a degree greater than expected from the modest variations in the component of velocity parallel to the jet axis.
As an alternative approach, code is currently being written which determines the streamline between each grid-point and the jet inlet; advecting slightly displaced particles from each streamline origin should enable the calculation of parallelepiped deformation for each grid-point. Author Contributions: This manuscript was primarily written by C.D.K., and revised by both T.V.C. and P.A.H.; likewise, the analysis was mostly performed by C.D.K., with suggestions for methodology, discussion and interpretation provided by the other two authors.