{'fraction_non_alphanumeric': 0.03277469041263729, 'fraction_numerical': 0.012383494508379644, 'mean_word_length': 4.629064039408867, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Numerical simulations of dynamics and emission from relativistic astrophysical jets

Broadband emission from relativistic outflows (jets) of active galactic nuclei (AGN) and gamma-ray bursts (GRBs) contains valuable information about the nature of the jet itself, and about the central engine which launches it. Using special relativistic hydrodynamics and magnetohydronamics simulations we study the dynamics of the jet and its interaction with the surrounding medium. The observational signature of the simulated jets is computed using a radiative transfer code developed specifically for the purpose of computing multi-wavelength, time-dependent, non-thermal emission from astrophysical plasmas. We present results of a series of long-term projects devoted to understanding the dynamics and emission of jets in parsec-scale AGN jets, blazars and the afterglow phase of the GRBs.


Introduction
Relativistic jets are found in a variety of astrophysical scenarios such as active galactic nuclei (AGN), microquasars, gamma-ray bursts (GRBs) and those resulting from tidal-disruption events (TDEs). While jet temporal and spatial timescales can vary from relatively small and short-lived (GRB jets) to very large and long-lived (AGN jets), what they have in common is that they are relativistic, seem to be very collimated and appear to be launched from a region very close to the accreting black hole (see e.g., [1] for a recent review and introduction to relativistic jets).
In this paper we review a series of projects whose aim is to understand the dynamics and emission from relativistic jets by focusing on specific scenarios appearing in relativistic jets of AGNs, blazars and GRBs. We begin by describing the numerical method in section 2, where we describe the equations of special relativistic magnetohydrodynamics and the code MRGENESIS (section 2.1), and then give some details about the treatment of non-thermal particles (section 2.2) and the radiative transfer (section 2.3), both of which are handled by the SPEV. In section 3 we present three of the applications of MRGENESIS +SPEV method, and we conclude in the section 4.

Numerical method
The numerical method we use consists of three parts, each dealing with a different aspect of the jet dynamics and emission. In the section 2.1 we discuss the equations of the special relativistic magnetohydrodynamics (RMHD). In sections 2.2 and 2.3 the treatment of non-thermal particles and emission, as well as that of the radiative transport is explained. Finally, the structure of the SPEV code is explained in 2.3.1.

Special relativistic magnetohydrodynamics
We solve the system of equations of ideal special RMHD [2], which consists of the mass and the total-energy momentum conservation 1 where the covariant derivative is denoted by ∇ α . ρ is the fluid rest-mass, u α is the fluid 4-velocity, and the energy-momentum tensor is defined as The Minkowski metric is η αβ = diag(−1, +1, +1, +1), and the total pressure is p * := p + |b| 2 /2, where p is the fluid pressure. b α is the magnetic field four-vector, Here B is the magnetic field 3-vector and v the fluid velocity 3-vector, whose Lorentz factor is defined as Γ := 1 − v 2 −1/2 . The total specific enthalpy is defined as where ǫ is the specific internal energy, which depends on the equation of state p = p(ǫ, ρ). We use the TM approximation to the Synge equation of state, where the specific enthalpy h is defined as [3] h = (5/2)(p/ρ)+ (9/4)(p/ρ) 2 + 1. Finally, the magnetic field has to satisfy the induction equation and the divergence constraint ∇ · B = 0. We use the code MRGENESIS [4,5,6,7] to numerically solve the system of equations of RMHD. MRGENESIS uses operator-splitting, finite-volumes, explicit method and employs approximate Riemann solvers for the computation of intercell fluxes, as well as total variation diminishing second and third-order Runge-Kutta methods for the time integration (for a recent overview of jet simulation codes see e.g. [8]). MRGENESIS has been massively parallelized using OpenMP and MPI libraries and has achieved reasonable scaling up to 10 4 cores on the MareNostrum machine at the Barcelona Supercomputing Centre.

Non-thermal particles and emission
As discussed in the section 1, a small fraction of very energetic particles ("non-thermal electrons" in the following) is responsible for the observed jet emission. In [9] we discuss in detail the method which deals with the spatial and temporal evolution of the non-thermal electrons in a magnetized relativistic fluid. Here we summarize the most important properties of the method.
2.2.1. Non-thermal electron evolution We assume that electrons are injected at relativistic shocks and use the phenomenological injection term which denotes the number of particles injected per unit volume and unit time, where γ is the electron Lorentz factor γ = E/(m e c 2 ), E being its energy and m e being its mass. Q 0 is the normalization of the injection energy spectrum, q is the power-law index and the function S(x; a, b) has the value 1 if a ≤ x ≤ b and 0 otherwise 2 . The temporal evolution of the electron energy distribution in the fluid comoving frame 3 is governed by the kinetic equation [10]: where n(t; γ) is the number density of electrons in an interval dγ around the Lorentz factor γ at the time t. Energy gains and losses are described by the termγ := dγ/dt. For sufficiently small intervals of time we can write the energy losses as where the adiabatic gains or losses are described by the quantity k a := d ln ρ dt , which denotes the compression or expansion of the fluid element, while the synchrotron losses are denoted by 3m 2 e c 2 where B ′ is the magnetic field strength in the fluid comoving frame and σ T is the Thomson cross section.
As described in section 3.2 of [9], we discretize the electron energy space in a number of logarithmically spaced bins and use a semi-analytic solver to accurately solve the equation (7) with a minimum computational effort.

Non-thermal radiation processes
We consider the following radiation processes in our simulations: synchrotron emission, synchrotron-self Compton (SSC) scattering, and the external inverse-Compton (EIC) scattering.
Synchrotron emission is produced when high-energy electrons gyrate around the magnetic field lines, producing a non-thermal, broadband emission spectrum. It is computed for each energy bin of the electron distribution using the interpolation method described in the sections 2.1.3 and 4.3.1 of [11].
Inverse-Compton scattering is a scattering of a low-frequency photon off a high-energy electron, whereby the frequency of the outgoing photon can be many orders of magnitude larger than that of the incoming one. In our numerical method it is computed analytically (see section 2.2.2 of [11]) assuming that the incoming emission spectrum is a power-law, as is the energy spectrum of the electrons off which the photons scatter. This method has been used in [12] to compute the EIC scattering of the ultraviolet stellar, photons, as well as in [13] to compute both the SSC and EIC emission in the blazar jets.

Radiative transfer
To compute the time-and frequency-dependent synthetic image we use the algorithm described in the Appendix A of [9]. Both this algorithm and the one described in the previous section are building blocks of the code SPEV. SPEV solves the radiative transfer equation where I ν , j ν and α ν are the specific intensity, emissivity and absorption coefficient at frequency ν. s is the distance along the line of sight which is seen by the observer at the time T . The relation between the time of observation and the time t in the laboratory frame (the frame in which the jet evolution is simulated) is given by where D is the distance of the jet injection point from the Earth and we choose that s = 0 at that point. As can be seen, in order to solve the equation (9), it is not enough to process each saved simulation snapshot at time t snapshot . Instead, it is necessary to process all snapshots simultaneously because of the fact that distance from the observer and the time at which the snapshot has been taken are not independent, but must satisfy the equation (10). In practice, this means that the radiative transfer equation (9) cannot be solved before all the information along the line of sight has been gathered. This makes the problem tightly coupled and highly non-local, complicating the code paralellization and increasing the CPU and memory costs.
2.3.1. Structure and paralellization of the SPEV code The evolution and emission from the nonthermal particles (section 2.2) and the computation of the resulting emission (section 2.3) are done by the SPEV code whose structure can be sketched as follows: after creating the initial Lagrangian particles (representing the pre-existing non-thermal electrons), SPEV iteratively runs through the three phases: evolution of existing particles, injection of new particles, and calculation of the emission (i.e., contribution to the pixels of a virtual detector on which the emission is computed). SPEV has two important characteristics very useful for its paralellization: the evolution of each Lagrangian particle during the simulation does not affect the evolution of the rest of the particles, and the contribution to each pixel of the detector is independent from the contribution to any other pixel. Therefore, SPEV is parallelized over particles and detector pixels. The data is partitioned into different sets and distributed among different computing units (e.g., multiple multicore-nodes connected by means of the network or by means of buses). Each set of particles can evolve separately without the need of any kind of communications during the iterative process. When the iterative process ends, data is transferred to the master node and then the radiative transfer equation (9) is solved for each pixel. More details about the parallel implementation are provided in [14,15]. This parallel approach has shown very good scalability on different HPC-systems, especially on medium scale multi-sockets multicores of up to 50 cores.

Applications
In this section we outline a number of application of MRGENESIS + SPEV approach to relativistic jets. In section 3.1 we discuss the radio emission from parsec-scale jets. The section 3.2 is devoted to address the highly variable and highly energetic emission from the blazars, while section 3.3 presents results of a long-term study of gamma-ray burst afterglows.

Parsec-scale jets
Parsec-scale jets are radio features on a scale of several light years. They are seen emerging from AGNs over periods of months to years. It is widely accepted that the radio maps are not direct observations of the physical quantities (density, pressure, velocity, composition, etc.) in the jet, but are influenced by a number of relativistic effects (Doppler shift, beaming) as well as by the degradation of the image due to a finite resolution of the radio telescopes. By using numerical simulations of the jet dynamics and, subsequently, by computing the synthetic radio maps we can test theoretical jet models and directly compare them to the observations, as well as study some of the events which are expected to occur in them 4 . In a previous work [9] we studied one such event, which is an injection of a perturbation into a steady jet shown on Figure 1. This jet is over-pressured and under-dense with respect to the surrounding medium and has a number of recollimation shocks, which can be seen as knots in the density profile shown on Figure 1.
The perturbation has a form of a temporary increase in the fluid velocity at the jet injection point located at the left grid boundary. This results in the development of shock and rarefaction waves which propagate down the jet and interact with the recollimation shocks, resulting in a complex, non-linear evolution of the perturbation and several trailing components (see section 7 of [9] or more details). Figure 2 shows the radio maps taken at six different epochs of the jet evolution. In the first three panels we can see the main component (perturbation) as a dark region propagating through the jet. In the last three panels we see a gradual reestablishment of the steady jet.
The contours show the image degraded to the resolution available with today's radio telescopes. As can be appreciated, many of the small-scale details of the jet evolution and emission are lost in the relatively strong blurring of the image due to the low resolution of the actually observable radio maps. Therefore it is important to use the hydrodynamic simulations and possess the ability to compute time-dependent radio maps from those simulations because this enables us to compute both the full-resolution and the degraded images and thus evaluate which of the observed features of the jet dynamics and emission are intrinsic to the jet and which are the artefact of the finite observational resolution.

Blazar jets
Blazars are jets pointing almost directly towards us. They are characterised by high variability of their emission, often in the form of flares observed by X-ray satellites 5 . It is thought that the cause of the flares are the collisions of parts of the jet with different velocities ("shells" in the following), whereby a fraction of the shell kinetic energy is dissipated by the shocks which propagate through the shells as a result of their collision, and a fraction of the dissipated energy is radiated and observed as a flare. This is called the internal shock scenario [18,19,20,21,22,23,6,24,3,13]. Figures 3 and 4 illustrate the internal shocks model and the typical distance scales at which shells collide in blazar jets.
We study the dynamics and emission from single shell collisions, using one-dimensional [23] and two-dimensional [21] hydrodynamic simulations, as well as one-dimensional    magnetohydrodynamic simulations [6]. We found that the lateral shell spreading is negligible on the length scales considered and that the one-dimensional approximation is accurate [21]. Furthermore, we find that the magnetization of the flow plays a crucial role in the efficiency of the conversion of kinetic into thermal energy (and subsequently into radiation). From the computational point of view we found the hydrodynamic simulations to be rather expensive and unsuitable for parameter space studies.
In the follow-up works [3,13] we simplify the hydrodynamic approach in order to be able to devote more resources to the computation of the emission and to be able to cover larger parameter space. We study a large number of shell collisions where both shells have the same energy and size, but their velocity and the degree of magnetization can vary. We denote the shell Lorentz factor by Γ := (1 − v 2 /c 2 ) −1/2 , where v is the shell velocity. The magnetization for cold and relativistic shells (e.g. p/ρc 2 ≪ 1 and Γ ≫ 1, where ρ and p are the shell rest-mass density and the thermal pressure) is defined as σ := B 2 /(4πΓ 2 ρc 2 ) . We compute the exact solution of the Riemann problem 6 for each shell collision and determine the strength and velocity of the shocks which are formed by the shell interaction. By covering a large parameter space of possible shell magnetization (we studied 10 6 shell collisions with σ ranging from 10 −6 to 10 3 for each shell independently) we find that the "sweet spot" for the efficiency of the kinetic energy conversion is not for σ ≪ 1, as might be expected (the smaller the σ, the easier it is to shock the fluid), but rather for (σ L = 1, σ R = 0.1), where subscripts L and R denote the faster (left) and the slower (right) shell, respectively [3]. This result makes the magnetized internal shock model viable. In the follow-up work [13] we compute the multi-wavelength, time-dependent emission from ≈ 10 3 shell collisions. 7   Figure 6. Average spectrum of the flare shown in figure 5. The total spectrum is shown in black, while the synchrotron, SSC and EIC contributions are shown in blue, red and green, respectively. The seed photons for the EIC are assumed to be monochromatic (frequency 2 × 10 14 Hz). The viewing angle (angle between the jet axis and the line of sight) is 5 • , and the jet is assumed to be located at redshift z = 0.5. Figures 5 and 6 show an example of the light curve and spectra produced by a collision of two moderately magnetized shells in a typical blazar jet (Rueda, Mimica & Aloy, 2012, in preparation). In figure 5 can see that the optical and γ-ray light curves peak somewhat earlier than the X-ray light curve. The reason is that the SSC emission (section 2.2.2) is produced by the scattering of the synchrotron emission which needs a finite time to travel from the place where it is produced to the site of its scattering. The sharp drop in the emission is related to the moment when shocks cross the shells and cease dissipating energy. Average spectra are shown in figure 6. The low-frequency emission is dominated by the synchrotron emission, while the high-energy spectrum is dominated by the SSC component. As can be seen, the correct computation of the SSC component is essential, and since this component is computationally most expensive, it puts an upper limit on a number of models which can be computed on today's machines (see footnote 7).

Gamma-ray burst afterglows
Gamma-ray bursts (GRBs) are produced by an internal dissipation of energy in the relativistic outflow. The GRB afterglow emission is produced by the interaction of the relativistic jet with the external (circumburst) medium. It is not known to what extent the jet is magnetized. If the outflow is initially dominated by the thermal energy (fireball model [26,27]) then it is expected that its σ ≪ 1. If, on the other hand, the outflow is initially dominated by the Poynting-flux [28,29], then it is possible to find a significantly magnetized jet at the onset of the afterglow phase.
We study the interaction of an arbitrarily magnetized GRB jet with the external medium. We model the interaction by considering the deceleration of a dense, magnetized shell as it moves through and compresses a constant density medium in front of it. One of the crucial differences between the evolution of a non-magnetized and a magnetized shell is in the existence of a reverse shock (a shock propagating through the shell): if σ is high enough, then the interaction between the shell and the external medium is not strong enough to shock the shell material itself. By analyzing timescales of the shell deceleration and the propagation of magnetohydrodynamic waves through the shell we find the conditions for the existence of a reverse shock depending on the GRB magnetization, energy, Lorentz factor, duration and the density of the external medium [30]. Figure 7 shows the parameter space spanned by the ejecta magnetization and a dimensionless parameter ξ (see figure caption for definition). As we can see, there is a substantial portion of the parameter space (shaded region), where the reverse shock does not form. Since it is absent, the particles are not accelerated and we expect to fail to observe an emission from this shock during the early afterglow, which would be consistent with the absence of such observations in most GRB optical afterglows [31].
In the follow-up works [7,32] we confirmed these results by means of numerical simulations. Using MRGENESIS we performed high-resolution, long-term one-dimensional relativistic magnetohydrodynamic simulation of a shell-medium interaction, and then use SPEV to compute the resulting light curves. Figure 8 shows an example of the optical light curve produced by nonmagnetized, weakly-magnetized and strongly-magnetizeed shells with the same initial energy, Lorentz factor and width. As can be seen, the strongly magnetized shell light curve does not possess the optical flash feature that the other two possess. Therefore, strong magnetization might be required for the paucity of the observed optical flashes to be explained. Since most of the observed GRBs have ξ in the range [0. 3,3], this would require σ > 0.1 for almost all GRBs [30,32]. This would imply that GRB jets could be strongly magnetized even at such late stages of their evolution as the afterglow phase, which puts constraints on the amount of magnetic dissipation that can happen at previous stage (prompt emission).
The simulations needed to correctly compute afterglow light curves are very computationally demanding due to the high resolution required in the direction of the jet motion (e.g. more , where E, ∆ 0 and Γ 0 are the ejecta initial energy, width and the Lorentz factor, and ρ ext is the ext. medium density. The shaded region shows the ejecta which do not form a reverse shock, while the vertical boundary at σ 0 = 0.3 indicates those ejecta whose radial spreading becomes significant and makes the formation of a reverse shock inevitable [30].  Figure 8. Optical (5 × 10 14 Hz) light curve produced by the non-magnetized (σ 0 = 0), weakly magnetized (σ 0 = 0.01 and stronglymagnetized (σ 0 = 1) shells (black, blue and green lines, respectively) whose initial conditions are E = 3 × 10 47 J, Γ 0 = 300 and ∆ 0 = 3 × 10 9 m. The non/weakly magnetized ejecta have a strong reverse shock which produces an optical flash, while the strongly magnetized ejecta light curve does not have this feature [32]. than 10 4 zones in the initial shell) and the long timescales needed to compute the light curve of sufficient duration (approximately 10 7 numerical iterations). Therefore, performing this type of calculations to cover a full parameter space of feasible GRB afterglows is not possible at the moment. However, the existence of re-scaling laws (see section 4.4 of [7]) can help us to perform a small number of simulations and then extend their results to cover a somewhat larger parameter space.

Conclusions
We have presented a framework consisting of a numerical relativistic magnetohydrodynamic code MRGENESIS and the non-thermal radiative transfer code SPEV. We show how its modular nature has enabled it to be successfully been applied to a number of current problems and issues in the relativistic jets physics, and how it enables us, via the computation of synthetic observations, to directly compare the results of our simulations with the observations from ground-and space-based telescopes. In the future we plan to improve the scaling of MRGENESIS to more than 10 4 cores, and also to add new radiative processes to SPEV (e.g., polarization, improved inverse-Compton, radiative transfer in curved space-times, etc.). The code scaling improvements will allow us to perform full 3D magnetohydrodynamic simulations and radiative transfer calculations, relevant for astrophysical scenarios where axial symmetry cannot be assumed.