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.
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.,  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.
2 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.
2.1 Special relativistic magnetohydrodynamics
We solve the system of equations of ideal special RMHD, which consists of the mass and the total-energy momentum conservation111In this section we assume a system of units where the speed of light and the factor is included in the definition of the magnetic field.
where the covariant derivative is denoted by . is the fluid rest-mass, is the fluid 4-velocity, and the energy-momentum tensor is defined as
The Minkowski metric is , and the total pressure is , where is the fluid pressure. is the magnetic field four-vector,
Here is the magnetic field 3-vector and the fluid velocity 3-vector, whose Lorentz factor is defined as . The total specific enthalpy is defined as
where is the specific internal energy, which depends on the equation of state . We use the TM approximation to the Synge equation of state, where the specific enthalpy is defined as  . Finally, the magnetic field has to satisfy the induction equation
and the divergence constraint .
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. ). MRGENESIS has been massively parallelized using OpenMP and MPI libraries and has achieved reasonable scaling up to cores on the MareNostrum machine at the Barcelona Supercomputing Centre.
2.2 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  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 , being its energy and being its mass. is the normalization of the injection energy spectrum, is the power-law index and the function has the value if and otherwise222We fix and determine the rest of the parameters by assuming a proportionality between the number and energy density of the non-thermal electrons and those of the thermal fluid (see equations 36 and 37 of ).. The temporal evolution of the electron energy distribution in the fluid comoving frame333We assume that non-thermal electrons are advected with the thermal fluid. The spatial distribution of non-thermal electrons is represented by Lagrangian tracer particles. is governed by the kinetic equation :
where is the number density of electrons in an interval around the Lorentz factor at the time . Energy gains and losses are described by the term . For sufficiently small intervals of time we can write the energy losses as
where the adiabatic gains or losses are described by the quantity
which denotes the compression or expansion of the fluid element, while the synchrotron losses are denoted by
where is the magnetic field strength in the fluid comoving frame and is the Thomson cross section.
2.2.2 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 .
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 ) 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  to compute the EIC scattering of the ultraviolet stellar, photons, as well as in  to compute both the SSC and EIC emission in the blazar jets.
2.3 Radiative transfer
To compute the time- and frequency-dependent synthetic image we use the algorithm described in the Appendix A of . 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 , and are the specific intensity, emissivity and absorption coefficient at frequency . is the distance along the line of sight which is seen by the observer at the time . The relation between the time of observation and the time in the laboratory frame (the frame in which the jet evolution is simulated) is given by
where is the distance of the jet injection point from the Earth and we choose that 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 . 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 non-thermal 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.
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.
3.1 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 them444Our work on blazar jets (section 3.2) as well as studies of the stability of initially very magnetized () jets (e.g. ) indicate that at parsec-scales distances AGN jets are expected to be at most very weakly magnetized (e.g. ). In a previous work  we studied one such event, which is an injection of a perturbation into a steady jet shown on Figure 2. 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 2.
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  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.
3.2 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 satellites555For an example of X-ray flares from a well known blazar Mrk 421 see e.g., . 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 4 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  and two-dimensional  hydrodynamic simulations, as well as one-dimensional magnetohydrodynamic simulations . We found that the lateral shell spreading is negligible on the length scales considered and that the one-dimensional approximation is accurate . 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 , where is the shell velocity. The magnetization for cold and relativistic shells (e.g. and , where and are the shell rest-mass density and the thermal pressure) is defined as . We compute the exact solution of the Riemann problem666We use the exact Riemann solver  for the one-dimensional ideal relativistic magnetohydrodynamic Riemann problem with magnetic fields perpendicular to the flow velocity. 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 shell collisions with ranging from to for each shell independently) we find that the “sweet spot” for the efficiency of the kinetic energy conversion is not for , as might be expected (the smaller the , the easier it is to shock the fluid), but rather for , where subscripts and denote the faster (left) and the slower (right) shell, respectively . This result makes the magnetized internal shock model viable. In the follow-up work  we compute the multi-wavelength, time-dependent emission from shell collisions.777 We could not compute the emission from all of the shell collisions studied in  because it was not feasible with present supercomputers (the computation time for the emission from shell collisions exceeds thousand hours).
Figures 6 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 6 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).
3.3 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 . 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 . Figure 8 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 .
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 non-magnetized, 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 , this would require 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 than zones in the initial shell) and the long timescales needed to compute the light curve of sufficient duration (approximately 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 ) can help us to perform a small number of simulations and then extend their results to cover a somewhat larger parameter space.
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 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.
PM, MAA and CA acknowledge the support from the ERC grant CAMAP-259276 and the grants AYA2010-21097-C03-01 and PROMETEO-2009-103. JMRB acknowledges the support from the Grisolia fellowship GRISOLIA/2011/041. PM, MAA, ST and CA acknowledge the support from the Consolider grant CSD2007-00050. The authors thankfully acknowledge the computer resources, technical expertise and assistance provided by the Barcelona Supercomputing Center - Centro Nacional de Supercomputación.
-  Böttcher M, Harris D E and Krawczynski H 2012 Introduction and Historical Perspective; in Relativistic Jets from Active Galactic Nuclei (eds M. BÃ¶ttcher, D. E. Harris and H. Krawczynski) (Wiley-VCH Verlag GmbH & Co. KGaA) chap 1, pp 1–16 ISBN 9783527641741
-  Anile A M 1989 Relativistic fluids and magneto-fluids: With applications in astrophysics and plasma physics (Cambridge University Press) ISBN 9780521304061
-  Mimica P and Aloy M A 2010 Mon. Not. R. Astron. Soc. 401 525–532 (Preprint 0909.1328)
-  Aloy M A, Ibáñez J M, Martí J M and Müller E 1999 Astrophysical J. Suppl. Ser. 122 151–166 (Preprint arXiv:astro-ph/9903352)
-  Leismann T, Antón L, Aloy M A, Müller E, Martí J M, Miralles J A and Ibáñez J M 2005 Astron. Astrophys. 436 503–526
-  Mimica P, Aloy M A and Müller E 2007 Astron. Astrophys. 466 93–106 (Preprint arXiv:astro-ph/0611765)
-  Mimica P, Giannios D and Aloy M A 2009 Astron. Astrophys. 494 879–890 (Preprint 0810.2961)
-  Aloy M A and Mimica P 2012 Simulations of Jets from Active Galactic Nuclei and Gamma-Ray Bursts; in Relativistic Jets from Active Galactic Nuclei (eds M. BÃ¶ttcher, D. E. Harris and H. Krawczynski) (Wiley-VCH Verlag GmbH & Co. KGaA) chap 10, pp 297–339 ISBN 9783527641741
-  Mimica P, Aloy M A, Agudo I, Martí J M, Gómez J L and Miralles J A 2009 Astrophysical J. 696 1142–1163 (Preprint 0811.1143)
-  Kardashev N S 1962 Soviet Astron. 6 317
-  Mimica P 2004 Numerical Simulations of Blazar Jets and their Non-thermal Radiation Ph.D. thesis Max-Planck-Institut für Astrophysik
-  Mimica P and Giannios D 2011 Mon. Not. R. Astron. Soc. 418 583–590 (Preprint 1106.1903)
-  Mimica P and Aloy M A 2012 Mon. Not. R. Astron. Soc. 421 2635–2647 (Preprint 1111.6108)
-  Tabik S, Mimica P, Plata O, Zapata E L and Romero L F 2011 HiPC 2011 1 1–10
-  Tabik S, Romero L F, Mimica P, Plata O and Zapata E L 2012 Computer Physics Communications 183 1937–1946
-  Giannios D and Spruit H C 2006 Astron. Astrophys. 450 887–898 (Preprint arXiv:astro-ph/0601172)
-  Brinkmann W, Papadakis I E, Raeth C, Mimica P and Haberl F 2005 Astron. Astrophys. 443 397–411 (Preprint arXiv:astro-ph/0508433)
-  Rees M J and Meszaros P 1994 Astrophysical J. Lett. 430 L93–L96 (Preprint arXiv:astro-ph/9404038)
-  Kobayashi S, Piran T and Sari R 1997 Astrophysical J. 490 92 (Preprint arXiv:astro-ph/9705013)
-  Daigne F and Mochkovitch R 1998 Mon. Not. R. Astron. Soc. 296 275–286 (Preprint arXiv:astro-ph/9801245)
-  Mimica P, Aloy M A, Müller E and Brinkmann W 2004 Astron. Astrophys. 418 947–958 (Preprint arXiv:astro-ph/0401266)
-  Kino M, Mizuta A and Yamada S 2004 Astrophysical J. 611 1021–1032 (Preprint arXiv:astro-ph/0404555)
-  Mimica P, Aloy M A, Müller E and Brinkmann W 2005 Astron. Astrophys. 441 103–115 (Preprint arXiv:astro-ph/0506636)
-  Bošnjak Ž, Daigne F and Dubus G 2009 Astron. Astrophys. 498 677–703 (Preprint 0811.2956)
-  Romero R, Martí J M, Pons J A, Ibáñez J M and Miralles J A 2005 Journal of Fluid Mechanics 544 323–338 (Preprint arXiv:astro-ph/0506527)
-  Goodman J 1986 Astrophysical J. Lett. 308 L47–L50
-  Paczynski B 1986 Astrophysical J. Lett. 308 L43–L46
-  Usov V V 1992 Nature 357 472–474
-  Meszaros P and Rees M J 1997 Astrophysical J. Lett. 482 L29 (Preprint arXiv:astro-ph/9609065)
-  Giannios D, Mimica P and Aloy M A 2008 Astron. Astrophys. 478 747–753 (Preprint 0711.1980)
-  Gomboc A, Kobayashi S, Mundell C G, Guidorzi C, Melandri A, Steele I A, Smith R J, Bersier D, Carter D and Bode M F 2009 Optical flashes, reverse shocks and magnetization American Institute of Physics Conference Series (American Institute of Physics Conference Series vol 1133) ed Meegan C, Kouveliotou C and Gehrels N pp 145–150 (Preprint 0902.1830)
-  Mimica P, Giannios D and Aloy M A 2010 Mon. Not. R. Astron. Soc. 407 2501–2510 (Preprint 1004.2720)