Kink Instability of Force-Free Jets: a Parameter Space Study
E. Sobacchi, Y. E. Lyubarsky, M. C. Sormani

TL;DR
This study investigates the kink instability in force-free relativistic jets, identifying key parameters influencing growth rates and suggesting the instability's role in energy conversion in jets like M87.
Contribution
The paper provides a comprehensive parameter space analysis of kink instability in force-free jets, deriving a simple formula for growth rates and applying it to real astrophysical jets.
Findings
Growth rate depends on magnetic field gradient and Lorentz factor.
Kink instability becomes nonlinear near the jet's acceleration cessation.
Supports kink instability as a mechanism for energy conversion in jets.
Abstract
In the paradigm of magnetic acceleration of relativistic jets, one of the key points is identifying a viable mechanism to convert the Poynting flux into the kinetic energy of the plasma beyond equipartition. A promising candidate is the kink instability, which deforms the body of the jet through helical perturbations. Since the detailed structure of real jets is unknown, we explore a large family of cylindrical, force-free equilibria to get robust conclusions. We find that the growth rate of the instability depends primarily on two parameters: (i) the gradient of the poloidal magnetic field; (ii) the Lorentz factor of the perturbation, which is closely related to the velocity of the plasma. We provide a simple fitting formula for the growth rate of the instability. As a tentative application, we use our results to interpret the dynamics of the jet in the nearby active galaxy M87. We…
Peer Reviews
No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.
Videos
No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.
Kink Instability of Force-Free Jets: a Parameter Space Study
E. Sobacchi1,2, Y. E. Lyubarsky1 & M. C. Sormani3
1 Physics Department, Ben-Gurion University, P.O.B. 653, Beer-Sheva 84105, Israel
2 Department of Natural Sciences, The Open University of Israel, 1 University Road, P.O.B. 808, Raanana 4353701, Israel
3 Institute for Theoretical Astrophysics, Zentrum für Astronomie der Universität Heidelberg, Albert-Überle-Str. 2, 69120 Heidelberg, Germany E-mail: [email protected]
Abstract
In the paradigm of magnetic acceleration of relativistic jets, one of the key points is identifying a viable mechanism to convert the Poynting flux into the kinetic energy of the plasma beyond equipartition. A promising candidate is the kink instability, which deforms the body of the jet through helical perturbations. Since the detailed structure of real jets is unknown, we explore a large family of cylindrical, force-free equilibria to get robust conclusions. We find that the growth rate of the instability depends primarily on two parameters: (i) the gradient of the poloidal magnetic field; (ii) the Lorentz factor of the perturbation, which is closely related to the velocity of the plasma. We provide a simple fitting formula for the growth rate of the instability. As a tentative application, we use our results to interpret the dynamics of the jet in the nearby active galaxy M87. We show that the kink instability becomes non-linear at a distance from the central black hole comparable to where the jet stops accelerating. Hence (at least for this object), the kink instability of the jet is a good candidate to drive the transition from a Poynting-dominated to a kinetic-energy-dominated flow.
keywords:
Magnetohydrodynamics (MHD) – Instabilities – Galaxies: jets – Galaxies: individual: M87
1 Introduction
Astrophysical jets are ubiquitous in a wide variety of events, ranging from small-scale protostellar objects to large-scale extragalactic jets. The jets from microquasars (e.g. Mirabel & Rodriguez 1999), Active Galactic Nuclei (AGN; e.g. Urry & Padovani 1995) and Gamma Ray Bursts (GRBs; e.g. Piran 2004) are accelerated to relativistic speeds.
One of the most promising explanations for jet launching is energy extraction from a rotating, magnetised source (e.g. Blandford 1976; Lovelace 1976; Blandford & Znajek 1977). In the simplest scenario, magnetic field lines anchored to the central object act as sliding wires for the plasma that is accelerated by the magnetic tension. In this scenario, one of the key points is the fate of the magnetic fields at large distances from the source. In the context of a steady, axisymmetric, ideal MHD flow, both analytical and numerical works have shown that the magnetic energy could be converted into the plasma kinetic energy up to equipartition (i.e. corresponding to a magnetisation ), but achieving further acceleration is generally difficult (e.g. Komissarov et al. 2007, 2009; Lyubarsky 2009, 2010, 2011; Tchekhovskoy et al. 2008, 2009, 2010).
Moreover, even for a relatively low magnetisation (), only weak shocks are possible, which makes the jet too radiatively inefficient to be consistent with observations of GRBs (Zhang & Kobayashi, 2005; Mimica et al., 2009a, b; Mimica & Aloy, 2010; Narayan et al., 2011). Spectral fitting of AGNs also require the plasma to be matter-dominated in the emission region, typically located at hundreds/thousands of gravitational radii from the black hole (Ghisellini et al., 2010; Tavecchio et al., 2011). Hence, it is crucial to identify some mechanism for efficient conversion of the magnetic energy into kinetic energy well beyond equipartition.
A possibility is that the instabilities in the MHD flow eventually destroy its regular structure and cause the release into the plasma of the energy stored in the magnetic fields (e.g. Lyubarsky 1992; Eichler 1993; Spruit et al. 1997; Begelman 1998; Giannios & Spruit 2006). In force-free jets, the most dangerous helical modes are indeed unstable if the poloidal magnetic field has a non-vanishing gradient (Istomin & Pariev, 1996; Lyubarsky, 1999).
In the context of relativistic jets, most of the analytical works focused on non-rotating flows (e.g. Appl et al. 2000; Bodo et al. 2013), or considered the limit of long/short wavelengths (e.g. Lyubarsky 1999; Tomimatsu et al. 2001; Nalewajko & Begelman 2012). Narayan et al. (2009) studied a two-parameter family of cylindrical, force-free equilibria, assuming a rigid impenetrable wall at the outer cylindrical radius; here we do not make this assumption, and we get complementary results. For the more relevant case of rotating jets, in this paper we address (i) the dependence of the maximum growth rate of the instability on the gradient of the poloidal field (only the limit of long wavelengths has been studied so far); (ii) the relation between the group velocity of the perturbation and the velocity of the plasma. As a first step, we consider force-free jets.
In general, analytical results have proven extremely useful to interpret numerical simulations that explored the evolution of the kink instability in the non-linear regime and revealed its fundamental effect on the jet’s structure after tens/hundreds of light crossing times (e.g. Nakamura et al. 2007; Mizuno et al. 2009, 2012, 2014; Mignone et al. 2010; O’Neill et al. 2012; Singh et al. 2016). Some consensus has eventually been reached on the fact that the kink instability can play an important role in the transition to a matter-dominated flow, at least for causally connected jets. Recently, different authors (Porth & Komissarov, 2015; Tchekhovskoy & Bromberg, 2016) even proposed that whether the jet becomes kink-unstable may explain the FRI/FRII dichotomy of radio galaxies (Fanaroff & Riley, 1974).
Since the detailed internal structure of real jets is unknown, it is important to identify the most fundamental physical parameters controlling how the instability develops. We attack this problem by selecting a large class of background solutions (Mizuno et al., 2012), and we explore a wide (three-dimensional) parameter space to get robust conclusions. We provide a simple analytic formula that approximately reproduces the growth rate of the kink instability and depends on (i) the gradient of the poloidal magnetic field; (ii) the Lorentz factor of the perturbation, which is closely related to the drift velocity of the plasma. As a tentative application, we use our result to interpret the dynamics of the jet in the active galaxy M87, recently resolved down to hundreds of gravitational radii from the central black hole (Mertens et al., 2016).
The paper is organised as follows. In Section 2 we present the background solution and the linearised equations for helical perturbations. In Section 3 we find the dispersion relation for the kink modes and we study in detail the group velocity and the growth rate of the instability. In Section 4 we discuss how our results can be applied to the jet of the active galaxy M87. Finally, in Section 5 we summarise our conclusions.
2 Fundamental equations
The fundamental equations governing relativistic jets in the ideal MHD approximation are the Maxwell’s equations
[TABLE]
where and j are the charge and current densities. These are coupled with the condition of infinite conductivity
[TABLE]
where v is the velocity of the flow. In the case of force-free flows, Euler’s fluid equation reduces to
[TABLE]
2.1 Unperturbed solution
The equation for the steady-state equilibrium configuration of a cylindrical, force-free jet can be derived from Eq. (1) - (4). It is
[TABLE]
where
[TABLE]
is the angular velocity of the field lines (with this definition ). Note that, since the toroidal magnetic field is generated by the rotation of the base of the jet, the sign of is opposite to the sign of (for jets propagating in the positive direction).
We consider a poloidal magnetic field of the form
[TABLE]
where is the typical scale of the jet core and the parameter defines the poloidal field profile (e.g. Mizuno et al. 2012). Note that is also equal to minus the logarithmic derivative of calculated at the core scale, .
With this poloidal field, Eq. (5) has an analytical solution for the toroidal field
[TABLE]
where
[TABLE]
is the ratio between the toroidal and poloidal components of the magnetic field in a non rotating jet. In the following we parametrize the angular velocity as
[TABLE]
where defines its slope at large radii.
Both the current density and the Poynting flux corresponding to Eq. (7) - (10) are peaked at . Nevertheless, this solution may formally give a finite/diverging total current, depending on the asymptotic (i.e. ) profile of the fields. The current flowing through a circular section of the jet with radius is . Since from Eq. (8) one can show that is never decaying faster than , the total current cannot formally vanish at infinity. However, if the fields sharply decline outside of some radius (which we may associate with the typical size of the accretion disk, whose magnetic field is confining the jet), the total current vanishes. Such a scenario does not significantly change the fields/currents in the most relevant region within a few , but produces small currents at large radii (i.e. ) which balance the total flux. Since the perturbations are well localised at the core, this slight modification does not affect their time evolution.
2.2 Linearised equations for the perturbed flow
Stability of cylindrical equilibria can be investigated considering perturbations on the magnetic surfaces of the form
[TABLE]
In this paper we use the method developed by Solov’ev (1967) for non-relativistic MHD flows. Lyubarsky (1999) has adapted this method to force-free jets, showing that the time evolution of the perturbation is described by a second order linear differential equation, namely
[TABLE]
The functions and can be expressed in terms of the unperturbed solution as
[TABLE]
where
[TABLE]
Here we have defined and . In general, the solution of Eq. (12) needs to satisfy two homogeneous boundary conditions, namely (i) an arbitrary normalisation; (ii) vanishing at infinity. Hence, the value of is automatically determined by these requirements.
In the following we focus on the mode, corresponding to the fastest growing instability (e.g. Bateman 1978). For completeness, we include a brief discussion of modes with higher in Appendix A. When is complex (the relevant case since we aim to study jet instabilities) the function does not vanish at any radius . For , at we have ; therefore, expansion at the first order close to the origin implies . With an arbitrary normalisation (e.g. ) we are left with the initial conditions of a standard Cauchy problem, and Eq. (12) can then be integrated numerically to find the solution for any .111Numerical integration is started from a positive (where is the typical scale of the jet core). Taylor expansion of the solution is used to find the proper initial conditions.
However, the asymptotic behaviour of the solution poses an additional constraint. For a given mode, when we have . Hence, in this limit . Since we are looking at small perturbations (i.e. ), the physical solution is the decaying exponential and we need the outer boundary condition .222We use a finite for the outer boundary condition, namely . We have checked that the choice of and does not affect the results. The value of matching this further condition can be found with the standard shooting method for eigenvalue problems (e.g. Press et al. 2002).
3 Results
This section is dedicated to study the solution of Eq. (12) for different configurations of the background fields. In Figure 1 we show the solution of Eq. (12) when , , (solid/dashed lines correspond to the real/imaginary parts of respectively). The wavelength is comparable to the size of the core (we use ). Note that the perturbation is well concentrated at the core of the jet, and practically vanishes at .
3.1 Dispersion relation
In Figure 2 we show the dispersion relation in a force-free jet for the imaginary and real parts of (/ in left/right panels respectively). As a fiducial model, we take , and . Then, we keep two of these parameters fixed and we study the effect of changing the other one (, and from top to bottom).
First of all note that the most unstable wavelengths are comparable but longer than the scale of the jet core (), and the typical growth rate of the instability is some fraction of the light crossing frequency ( for the parameter range considered here).
For a given triplet the right panels in Figure 2 show that the group velocity of the perturbation is nearly constant over almost the entire unstable region. However, it shows significant variations at the shortest unstable wavelengths (largest ), also exceeding by in few cases. This is due to the fact that, when , the solution of Eq. (12) has a resonance which is difficult to treat numerically. Since these complications arise only in a limited region (where the instability is weak), they do not affect the general behaviour of the perturbation.
The maximum growth rate of the instability, , decreases with , while the group velocity of the perturbation, (calculated for the most unstable ), increases. For example, and when . This result is consistent with the analytical study of Lyubarsky (1999), who found a relativistic suppression of the instability in the limit of long wavelengths. Note that the most unstable increases with in the same range.
When the poloidal field is nearly flat (i.e. small ), the instability is growing slowly, and completely disappears when (Istomin & Pariev, 1996; Lyubarsky, 1999). The poloidal field profile parameter has a strong impact on the dispersion relation: there is one order of magnitude difference in the growth rate () when . In the same range, the group velocity is decreasing, though less significantly ().
Both the maximum growth rate and the group velocity of the instability are almost independent on the asymptotic slope of the angular velocity, . In the range , the variation is for and for (note that the difference in is larger, but the slope for the most unstable is almost constant).
3.2 Group velocity of the perturbation
The Lorentz factor of the perturbation (calculated for the most unstable wavelength, ) is fundamental to understand its time evolution. Therefore, one would like to connect it to the Lorentz factor of the plasma, which is easier to visualise physically. The problem is that in force-free jets the Lorentz factor is generally undetermined.
However, sufficiently far from the source the velocity of the plasma approaches a pure drift motion in the electromagnetic fields, i.e. (Tchekhovskoy et al., 2009). In this case, for our setup one can easily calculate the Lorentz factor of the plasma as
[TABLE]
where and are given by Eq. (9) and (10) respectively. Since depends on the radius, we take the Lorentz factor at the core scale, , as a proxy for the typical Lorentz factor of the plasma. It is possible to show that peaks around this scale, and that the drift velocity eventually vanishes at infinity (i.e. when ).
To study the connection between and , we select random points in our three-dimensional parameter space (in the range ; ; ) and we calculate the corresponding dispersion relations. The upper limit on and the lower limit on are selected in order to avoid the case , which is difficult to treat numerically.
In Figure 3 we plot the Lorentz factor of the perturbation () versus the typical Lorentz factor of a pure drift motion () for the 15 different choices of the parameters. As we can see, and are clearly correlated; this is a natural result since the instability develops in the plasma comoving frame. In particular, for all the combinations of parameters we used, we found (dashed lines in the figure). Hence, one can identify the two Lorentz factors within reasonable accuracy.
3.3 Growth rate of the instability
In general, it would be useful to have a quick way to estimate the growth rate of the instability, without the need to recalculate the dispersion relation for each different jet structure we are interested to study. This is particularly true while considering realistic jets, where the detailed structure of the fields is unknown and we necessarily rely on order-of-magnitude estimates of the relevant physical parameters.
In Figure 4, the maximum growth rates (in the rest frame of the perturbation, ) for the 15 background solutions considered above are shown as dots as a function of . Since there is a relatively small scatter at a given , it is possible to estimate the growth rate of the instability as
[TABLE]
where the solid line in Figure 4 corresponds to
[TABLE]
For comparison, the dashed lines show the same Eq. (20) with different coefficients ( and for the lower/upper curves). Note that the growth rate vanishes when (i.e. for a flat poloidal field). This is particularly important while interpreting the results of numerical simulations, since the stability of the jet crucially depends on the gradient of the poloidal field at the core.
Using Eq. (19) we can calculate the characteristic growth time of the instability as
[TABLE]
where we have substituted from Eq. (20). In general, one would expect the kink instability to become non-linear after few . This result is in general agreement with numerical simulations with a similar setup, typically finding that the kink instability significantly develops for a time , before saturating in the fully non-linear regime (e.g. Mizuno et al. 2009, 2012, 2014; O’Neill et al. 2012; Singh et al. 2016). Moreover, Mizuno et al. (2012) explicitly considered the effect of the poloidal field gradient on the perturbation, showing that for small the instability is severely suppressed also in the non-linear regime.
4 Implications for the jet of M87
In the paradigm of magnetic launching, the kink instability is often invoked to explain how a jet can transfer its energy from the Poynting flux to the kinetic energy of the plasma. For this mechanism to be efficient, the jet needs to be strongly causally connected, i.e. , where and are the opening angle and the Lorentz factor of the jet (e.g. Komissarov et al. 2009; Lyubarsky 2009; Tchekhovskoy et al. 2009; Granot et al. 2011; Porth & Komissarov 2015). This condition can be verified in AGN (though with a large scatter; e.g. Pushkarev et al. 2009; Clausen-Brown et al. 2013), while it is violated by GRBs (Kumar & Zhang 2015 and references therein).333McKinney & Blandford (2009) performed a 3D simulation of a rapidly rotating black hole producing an approximately conical jet, with and . They found this jet to be stable, retaining a high magnetisation out to gravitational radii despite small wiggles interpreted as a signature of the kink modes. Interestingly, this jet is at the boundary for strong causality (they have ). Moreover, the kink instability typically becomes non-linear around the largest scale they simulated (see below).
A fundamental prototype of collimated jet is that in the active galaxy M87. Mertens et al. (2016) recently resolved the dynamics of this jet down to hundreds of gravitational radii from the central black hole. This galaxy is therefore an ideal case to study the mechanisms acting during the launching of the jet. The jet has an approximately parabolic shape, and the plasma accelerates linearly out to (where is the gravitational radius) from the source. For the central black hole we assume a mass of (e.g. Walsh et al. 2013), corresponding to a gravitational radius . During this phase, the typical transverse scale of the jet spans a range (see their Figure 6). After the linear acceleration phase, the Lorentz factor slowly increases for several orders of magnitude in distance.
Even in perfectly collimated jets, the kink instability becomes non-linear only at a distance from the source. Using , our Eq. (21) gives
[TABLE]
where, according to the discussion above, we have identified the Lorentz factors of the perturbation and of the plasma.444Formally, our Eq. (21) was derived considering force-free jets. Of course, the effect of a finite magnetisation (even if ) and the presence of a confining external medium can affect the result. However, we believe that at least the order of magnitude is preserved. Interestingly, this length scale is comparable with the end of the linear acceleration regime, i.e.
[TABLE]
This suggests the following scenario: (i) close to the central engine, the flow accelerates while dominated by the Poynting flux, or at most in a state of equipartition; (ii) at a typical distance the kink instability enters its non-linear regime, eventually transferring the energy of the jet from the Poynting flux to the plasma (note that also blazar observations require the energy conversion to be completed around this scale; e.g. Ghisellini et al. 2010; Tavecchio et al. 2011); (iii) farther away, the jet is dominated by the kinetic energy of the plasma and the acceleration almost stops.
Of course, additional (potentially large) uncertainties are due to the unknown value of . However, at least when the jet is accelerated from down to , the poloidal flux is concentrated in the vicinity of the axis (Beskin & Nokhrina, 2009; Lyubarsky, 2009), and a typical seems a reasonable description for the core of the jet (Tchekhovskoy et al., 2009).
In general, one would expect a strong dissipation in the region where the global structure of the jet is destroyed by the kink modes. Unfortunately, resolving the jet down to these scales in wavebands different than the radio is not feasible with current facilities. However, the time variability of the light curves provides interesting constraints on the size of the emitting region. For M87, variabilities on extremely short time scales () have been detected in the TeV region, a factor faster than in other bands (Aharonian et al., 2006). This puts an upper limit on the size of the emitting region, (where is the Doppler factor of the jet and ). Using and a viewing angle (Mertens et al., 2016),555Note that for M87 we have . Hence, it is possible to observe the emitted radiation despite beaming. one eventually finds . Since this upper limit is comparable with the typical transverse scale of the jet, , the bulk of the emission may come from the same region where the kink instability becomes non-linear.
In the first regime of the outlined scenario (i.e. Poynting-dominated flow), theoretical models predict oscillations of the jet cross section (Lyubarsky, 2009); alternatively, some oscillations may be due to the kink modes while still in the linear regime. Also in the third regime, the Kelvin-Helmholtz instability of the kinetic-energy-dominated plasma can result in similar patterns (Hardee 2000; Lobanov et al. 2003; see also Lobanov & Zensus 2001). These regimes may correspond to the oscillations of the instantaneous opening angle detected in the M87 jet, at and respectively (Mertens et al. 2016, their Figure 6). In this interpretation, the intermediate case (i.e. , where oscillations are not clear), would correspond to the transition from a Poynting to a kinetic-energy dominated jet, driven by the kink instability in its non-linear regime.
5 Conclusions
We have explored the effect of the kink instability on cylindrical, force-free jets in the linear regime. To get robust conclusions, we have considered a large class of background solutions (Mizuno et al., 2012). In principle, the growth rate can depend on all the parameters describing the background fields, namely: (i) the angular velocity, ; (ii) its slope at large radii, ; (iii) the logarithmic derivative (calculated at the core scale) of the poloidal magnetic field, .
Calculating the dispersion relation for different combinations of the parameters, we have found that the group velocity of the perturbation, , is closely related to the velocity of the plasma in the background fields. Hence, one can estimate
[TABLE]
where is the typical Lorentz factor of the plasma and . This is a natural result since the instability develops in the plasma comoving frame.
The growth rate of the instability (corresponding to the most unstable wavelength) can be expressed in terms of and , while it is insensitive to . We have provided a simple equation reproducing our results:
[TABLE]
where is the scale of the core (note that narrow jets are more unstable). In particular, the growth rate is suppressed due to time dilation from the rest frame of the perturbation (see the factor in the equation above). We also confirm previous results (Lyubarsky, 1999; Mizuno et al., 2012), finding that the growth rate of the perturbation is severely suppressed for a nearly flat poloidal field (i.e. ).
Applying these results to the well resolved jet of the active galaxy M87 (Mertens et al., 2016), we have shown that the kink instability becomes non-linear at a distance from the central black hole comparable to where the jet stops accelerating. This scenario is broadly consistent with both (i) the oscillations of the instantaneous opening angle of the jet, which may be due to the kink/Kelvin-Helmholtz modes, in the Poynting/kinetic-energy dominated regimes respectively; (ii) the variability of the light curve (at all wavelengths shorter than the radio, and in particular at TeV energies; Aharonian et al. 2006), suggesting that the size of the region where bright emission is expected due to dissipation is comparable to the transverse scale of the jet at the relevant distance from the source (i.e. where the kink modes become non-linear). Hence (at least for this object), we have suggested that the kink instability of the jet may be the mechanism driving the transition from a Poynting-dominated to a kinetic-energy-dominated flow.
Acknowledgements
ES and YEL acknowledge support from the Israeli Science Foundation under Grant No. 719/14. MCS acknowledges support from the Deutsche Forschungsgemeinschaft in the Collaborative Research Center (SFB 881) “The Milky Way System” (subprojects B1, B2, and B8) and in the Priority Program SPP 1573 “Physics of the Interstellar Medium” (grant numbers KL 1358/18.1, KL 1358/19.2). MCS furthermore thanks the European Research Council for funding in the ERC Advanced Grant STARLIGHT (project number 339177).
Appendix A Instability of the high-m modes
When , solving Eq. (12) one needs to pay attention to the fact that , but when . Hence, we have to take the proper initial conditions, i.e. and (here the value of the derivative corresponds to an arbitrary normalisation).
In Figure 5 we show the dispersion relation for , where different peaks correspond to modes with different . At high , the peak shifts to large and the instability is suppressed as expected. Interestingly, this suppression is stronger for the dashed () than for the solid curves (); hence, at least in the linear regime, the high- modes may be even less important in the relativistic case. In particular, we were unable to find the dispersion relation when and , probably because of numerical issues when the jet becomes almost stable (i.e. ).
The reference list from the paper itself. Each links out to its DOI / PubMed record.
- 1Aharonian et al . (2006) Aharonian F. et al., 2006, Science, 314, 1424
- 2Appl et al . (2000) Appl S., Lery T., Baty H., 2000, A&A, 355, 818
- 3Bateman (1978) Bateman G., 1978, MHD Instabilities. MIT Press, Cambridge, MA
- 4Begelman (1998) Begelman M. C., 1998, Ap J, 493, 291
- 5Beskin & Nokhrina (2009) Beskin V. S., Nokhrina E. E., 2009, MNRAS, 397, 1486
- 6Blandford (1976) Blandford R. D., 1976, MNRAS, 176, 465
- 7Blandford & Znajek (1977) Blandford R. D., Znajek R. L., 1977, MNRAS, 179, 433
- 8Bodo et al . (2013) Bodo G., Mamatsashvili G., Rossi P., Mignone A., 2013, MNRAS, 434, 3030
