A new method for extending solutions to the self-similar relativistic magnetohydrodynamic equations for black hole outflows

Abstract

1 INTRODUCTION

Since their discovery, relativistic collimated outflows of matter have been observed in many astrophysical objects and they are known to be associated with accretion flows. Jets reveal themselves at different scales and redshifts, showing an extreme diversity in energetics, shapes and emission. Jets are found to be characteristic features of black holes (BH) systems, such as X-ray binaries (XRBs) and active galactic nuclei (AGN), as well as of young stellar objects (YSOs), explosive transients such as tidal disruption events (TDEs) and gamma-ray bursts (GRBs). Observations suggest that jets are an energetically important component of the system that hosts them, because the jet power appears to be comparable to the accretion power (see e.g. Rawlings & Saunders 1991; Nemmen & Tchekhovskoy 2015, for a more recent discussion). Significant evidence has been found of the effect of jets not only on the immediate proximity of the central object but also on their surrounding environment, where they deposit the energy extracted from the accretion flow (e.g. Gallo et al. 2005; Fabian 2012). To launch, accelerate and collimate a relativistic outflow over such large distances, magnetic fields need to be invoked. Observational evidences, such as polarization measurements both in the radio (Martí-Vidal et al. 2015) and in the hard X-rays (Laurent et al. 2011), support the idea that ordered magnetic fields are a key ingredient in the jet phenomena, and have a significant effect on their emission as well. Understanding what causes the jet to be launched from an accretion disc, which mechanisms determine its shape and extension and where the radiation is produced in these systems is one of the most fundamental questions in astrophysics and needs to be addressed promptly.

Multiwavelength continuum emission from jets and discs is observed for all accreting, jet-launching sources, XRBs in particular can provide essential pieces of information because extensive monitoring of these sources shows that they go through a duty cycle multiple times during their lives (Fender, Belloni & Gallo 2004). During an outburst, they spend most of their time in the low-luminosity hard state where they exhibit a mildly relativistic steady-state jet launched from a likely recessed disc. In this state, the jet dominates the total power and shows a characteristic power-law spectrum that may extend to high energies. Many outbursts show a rapid increase in luminosity, bringing the system close to its Eddington limit. The standard paradigm has the disc inner radius moving closer to the BH while the jet becomes ballistic, emitting superluminal knots. Eventually the jet switches off and the emission is dominated by the disc. The spectrum becomes softer and loses its non-thermal high-energy component (high-luminosity soft state). Finally, the system slowly decays into a hard state and the cycle restarts (see e.g. Belloni & Motta 2016, for a recent review). When the emission is dominated by the compact steady-state jet, its characteristic synchrotron emission can span several orders of magnitude in frequency. At the wavelength range where the synchrotron transitions from the optically thin (τ < 1) to the optically thick regime (τ > 1), the spectrum shows a break and it becomes a flat/inverted power law, characteristic of the self-absorbed synchrotron (see Romero et al. 2017, for a recent review). The region in the jet corresponding to the break frequency is believed to be the site where particles are first accelerated (Markoff 2010), potentially by internal shocks (Malzac 2014). The jet break can be inferred from observations to occur over a fairly large range of distances from the BH ∼ 10–10⁴ gravitational radii (hereafter r_g), and it has been seen to span 4 orders of magnitude in frequency during a state transition in a single object (Russell et al. 2014). Similar spectral features and a duty cycle 10⁷–10⁸ times longer are seen in AGN as well, and a power-unification scenario has been proposed independently by Merloni, Heinz & di Matteo (2003) and Falcke, Körding & Markoff (2004). They showed that AGN and XRBs of the same relative luminosity, rescaled by the mass of the BH, can be explained with the same physics framework. While XRBs provide unique constraints on the emission of jets and accretion discs thanks to the multiwavelength monitoring of the activity of the source during state-transition episodes, AGN are ideal for studying the structure of jets and the dynamical processes that shape them.

Indeed, in the case of nearby AGN, high spatial resolution observations are now possible with very long baseline interferometry (VLBI) in both cm and mm bands. These data provide unprecedented constraints on the geometry and the dynamics of jets, such as the jet opening angle, the height and the width of knots associated with standing shock features, such as Hubble Space Telescope (HST)-1 in M87's jet. Recently, high-resolution VLBI observations of M87 by Hada et al. (2016) resolved and imaged the inner core of the galaxy down to ∼10 r_g, revealing the innermost structure of the jet. VLBI/VLBA observations (Asada et al. 2014; Mertens et al. 2016) constrained the bulk Lorentz factor of the jet of M87 to be mildly relativistic, i.e. γ_j ∼ 1–3. In the near future, with the beginning of the Event Horizon Telescope (EHT) era, observations will resolve the nearest BH (Sagittarius A* and M87) down to the event horizon scale (see e.g. Doeleman et al. 2008). This unprecedented resolution will shed new light on the immediate proximity of BH, possibly unveiling the jet/disc connection and the mechanisms responsible for the acceleration and collimation in the first stages of jet formation. Finally, Meyer et al. (2013) measured the proper motions of the knots downstream of HST-1 with the HST, finding significant evidence of transverse and parallel motion with respect to the jet axis. This is evidence of a helical magnetic field beyond HST-1 and it brings important constraints on the modelling of jets at larger distances from the BH.

Using both XRBs and AGN to obtain insight on the apparently similar jet phenomena is extremely important as demonstrated by the activity in this field of research. However, an adequate modelling of jets, including a detailed treatment of both the radiative processes and the magnetohydrodynamics (MHD), is still far from being achieved.

Thanks to the dramatic improvement of computational power, accretion discs and jets can be modelled with general relativistic magnetohydrodynamic (GRMHD) simulations (e.g. Koide et al. 2002; McKinney 2006; Hawley & Krolik 2006; Tchekhovskoy, Narayan & McKinney 2011; Tchekhovskoy & Bromberg 2016). Full three-dimensional simulations allow detailed study of the stability of jets under different sets of initial conditions, and they provide a unique overview of how jets are launched, how the disc and jet interact during this phase and how they approach a stable configuration. However, the lack of crucial ingredients such as non-ideal processes or self-consistent radiative processes makes a direct test against observational data still a challenge. It is, however, worth noting that efforts are currently made to incorporate simplified treatments of electron microphysics in GRMHD simulations (e.g. Ressler et al. 2015; Mościbrodzka, Falcke & Shiokawa 2016).

A complementary method to simulations is given by semi-analytical models for multiwavelength emission from disc/jet systems, where the geometry and the dynamics are generally fixed and simplified. Many of such models have been proposed, for instance by Romero et al. (2003), Markoff, Nowak & Wilms (2005), Yuan, Cui & Narayan (2005), Potter & Cotter (2012), Zdziarski, Lubiński & Sikora (2012), Zdziarski et al. (2014), Pepe, Vila & Romero (2015) and they have been successful in reproducing the spectral energy distribution of accretion discs and jets. However, the treatment of magnetic fields is also greatly simplified and its orientation is usually not considered. Finally, they all present a certain degree of degeneracy between combinations of input parameters that give statistically equivalent fits to the same data set. Introducing a self-consistent treatment of MHD via semi-analytical models can help reduce the freedom in the parameter space, allowing for better constraining the fits to the observations.

A number of such models have been developed since the 80s-90s, with pioneering works by Blandford & Payne (1982), Lovelace & Contopoulos (1990), Li, Chiueh & Begelman (1992), Contopoulos (1994), Sauty & Tsinganos (1994), Bogovalov & Tsinganos (1999). The MHD system of equations describing an accelerating flow is a highly non-linear system that changes nature from elliptical to hyperbolic several times across the interval of integration. Independently of the geometry of the system, it has been shown that three critical surfaces exist in correspondence to such transitions and they determine as well the onset of magnetosonic waves of different type. These are called fast, slow magnetosonic and Alfvén waves. The slow magnetosonic singular surface (SMSS) appears close to the central object and upstream of the Alfvén surface. The fast magnetosonic singular surface (FMSS) is located downstream of the Alfvén point (AP) and it is suspected to be linked with a recollimation of the streamlines describing the flow (see e.g. Meier 2012, for a complete discussion and derivation). The FMSS, therefore, could tentatively be identified with the jet break seen in observations of AGN and XRBs, while the properties of the flow at SMSS could instead provide important constraints on the input parameters of radiative transfer models, such as e.g. the magnetic-to-thermal energy ratio at the base and the initial bulk velocity. However, the MHD equations that exhibit all the three singular surfaces must include gravity and the effect of thermal pressure in the total internal energy of the flow to properly describe the region close to the BH. Moreover, to describe typical astrophysical jets, the equations need to allow the flow to become relativistic. Vlahakis et al. (2000, herefter VTST00) and Vlahakis & Königl (2003, hereafter VK03) derived the MHD system of equations under the assumption of radial self-similarity, first including a simplified gravity term (kinetic term) and, later removing it, to include relativistic effects. Only recently, Polko, Meier & Markoff (2010, 2013, 2014, hereafter, respectively, PMM10, PMM13 and PMM14) derived the equations for a self-similar relativistic MHD flow including enthalpy and gravity, and found solutions crossing smoothly all three singular surfaces. Under the assumption of self-similarity, the singular surfaces are cones and intersect a streamline only in one point. They are usually referred to as ‘modified’ slow point (MSP) and ‘modified’ fast point (MFP) and AP. We will adopt this terminology from here on. For a more detailed history of the derivation of the model, we address the interested reader to the papers cited above and references therein.

Since only at the AP analytical formulae can be written explicitly, locating the MFP and MSP and performing the integration of the equations through them is not an easy numerical problem to solve. These unknowns add a large degree of complexity to the problem and common integration techniques fail to retrieve solutions in most cases. Although successful in solving such highly non-linear set of MHD equations, PMM14 were limited in the range of jet solutions that they could retrieve due to their numerical approach. Following VK03, PMM14 developed an algorithm where the integration starts from the AP and ‘shoots’ towards the other two singular points (MFP and MSP). The AP can be regularized analytically by using the De L'Hôpital rule on the terms that are in indefinite form 0/0. The other two singular points are later found by extrapolation (for more details about the shooting method, see e.g. Press et al. 1993 and Stoer & Bulirsch 2013). When the equations are integrated towards a singularity, however, the error with respect to the exact solution can be large. As a consequence of the limitations intrinsic to the adopted numerical scheme, PMM14 were limited to a small portion of the parameter space.

In this paper, we present a new method for exploring the parameter space and finding solutions to the MHD equations building on the work of the PMM papers. Using this method, we are able to explore more efficiently a much larger fraction of the parameter space. This will allow the model to be applicable to many other difficult flow solution problems in astrophysics, as well as potentially other fields. Moreover, this class of models can be coupled with fairly accurate radiative transfer models (see Markoff et al. 2005; Maitra et al. 2009; Connors et al. 2017; Crumley et al. 2017), which can handle spectral fitting in a reasonable computational time.

This paper follows the following structure. In Section 2, we describe the system of equations that we solve and discuss modifications in the equations compared to PMM14 with respect to their dependence on the gravitational potential. In Section 3, we present our numerical scheme, and in Section 4, we compare our results with PMM14. We show that the solutions are extremely sensitive to the gravity terms and when the corrected functions of the gravitational potential are used, the self-similarity assumption is more easily broken. In Section 5, we present a partial study of the parameter space, unaccessible to previous studies. We included the details of our method in a series of Appendixes. In Appendix A, we define the equations that we solve in explicit form with the corrected gravity terms and the new Alfvén regularity condition (ARC). In Appendix B, we describe the derivation of the functions of the pseudo-potential in the gravity terms. In Appendix C, we discuss in detail our approach in finding the locations of the unknown singular points, MFP and MSP, and how we perform the integration. Finally, in Appendix D, we give the conversion of the most relevant quantities into physical units.

2 SET-UP OF THE EQUATIONS

Our goal is to describe an outflow launched from an accretion disc in the presence of a magnetic field by the Blandford–Payne mechanism (Blandford & Payne 1982). We use cylindrical coordinates (ϖ, ϕ, z) and spherical coordinates (r, ϕ, θ) and impose axial symmetry, i.e. ∂/∂ϕ = 0 and assume the flow to be stationary, ∂/∂t = 0. For the electric field E, we impose the freeze-in condition, |$\boldsymbol {E} = -\boldsymbol {v} \times \boldsymbol {B}$|⁠, for the Ohm's law in ideal MHD and due to the axial symmetry E_ϕ = 0. In what follows, we adopt the notation of VK03 and PMM14.

As in VK03, we assume that the dependence of each variable on r and θ is separable and that the radial dependence can be expressed as a power law of r (see VK03). This assumption leads to self-similar solutions that obey a system of ordinary differential equations plus algebraic constraints. The four unknowns are the specific relativistic enthalpy ξc², the poloidal Alfvén Mach number, M = (γV_p/B_p)(4πρ₀ξ)^1/2 (V_p and B_p are the poloidal components of the velocity and the magnetic field, ρ₀ is the baryon rest mass density, γ is the bulk Lorentz factor), the cylindrical radius x in units of the light cylinder, x = ϖΩ/c (c is the speed of light, Ω is the angular velocity of the flow) and the angle ψ between the streamline and the horizontal axis. These are functions of θ only. We further substitute x = x_AG, such that for each streamline, G is the cylindrical radius scaled by its value at the AP.

Equation (1) is the so-called wind equation for the poloidal Alfvén Mach number, M, and can be derived from the Euler equation. Equation (2) describes the evolution of the dimensionless cylindrical radius and can be derived from the Euler equation. Equation (3) is the Bernoulli equation describing the energy conservation along the poloidal component of the magnetic field line. From this equation, we derive ψ, once M² and G² are known at each integration step. Finally, equation (4) is the relation between the poloidal M, the enthalpy ξ and the dimensionless adiabatic parameter q (Table 1 ), which we use to derive ξ at each step.

However, we noticed that the rate of change of ψ with respect to θ as derived from the system of equations (1)–(4) was not consistent with the prediction of equation (5), with a substantial discrepancy upstream of the AP (see the right-hand panel of Fig. 1).

Figure 1.

Open in new tabDownload slide

Left-hand panel: System of coordinates we adopt to describe a solution of equations (1)–(4), which is typically the ‘reference’ streamline identified with the label |$\alpha \equiv \varpi _A^2/\varpi ^2 = 1$|⁠. The four unknowns in the system (1)–(4), together with all the other quantities, are functions of θ, which is the angle between a point in the streamline and the z-axis. The angle that the tangent to the streamline makes with the horizontal axis is ψ, while the distance from a point of the streamline to the z-axis is defined by its cylindrical radius in units of the light cylinder x. Right-hand panel: Discrepancy between the derivative of ψ, calculated directly from equation (5), dψ/dθ, and inferred numerically at each step from equation (3), Δψ/Δθ with the corrected functions of the gravitational potential (top panel) and with the one used by PMM14 (bottom panel). The MSP and the BH are on the right. [A colour version of this figure is available in the online version.]

The reason can be understood as follows: Meier (2012) obtained the Bernoulli and the transfield equations neglecting terms |$\propto P_{\rm g}^2$| (see equations F.16 and F.18 in Meier 2012). Our equation (3) corresponds to Meier's equation F.16. In order to obtain equation (1), we take derivatives of the Bernoulli equation (equation 3) and therefore of P_g (see Appendix B). However, we also keep the Bernoulli equation in its original form. If we were to approximate the term 1/(1 − P_g) in the derivatives, the resulting equation would not be consistent with equation (3) anymore, because we would obtain a different version of the Bernoulli equation when integrating them back. Of course, this difference would be more pronounced upstream of the AP where P_g is greater, explaining the discrepancy we measured.

We use equation (3) in order to evaluate quickly ψ and therefore we need to retain in C₁ the full term 1/(1 − P_g) and calculate full derivatives of P_g with respect to θ when needed.¹ This approach brings back self-consistency between equations (1)–(4) and equation (5).

Analogously, the gravity correction term to C₂ comes from derivatives of P_g in the transfield equation written to the same order of approximation as the Bernoulli equation (equation F.18 of Meier 2012). For keeping the same consistency as mentioned before, we keep the full terms of the derivatives of P_g in C₂ as well (Appendix B).

3 NUMERICAL METHOD

The system of equations (1)–(4) allows for different families of solutions that may or may not cross the singular point(s), similar to the case for hydrodynamic (Parker 1958) or magnetohydrodynamic (Weber & Davis 1967) solar winds. In the case of the system of equations (1)–(4) and under the assumption of self-similarity, there can be up to three singular points: the AP, the modified fast magnetosonic point MFP and the modified slow magnetosonic point MSP (see Fig. 2). There is only one family of solutions that crosses all the points while it accelerates away from the BH (see figs 1 and 2 of Weber & Davis 1967, solution u_α1). Therefore, we shaped our approach in such a way that we automatically select solutions having this topology.

Figure 2.

Open in new tabDownload slide

Three-dimensional plot of the streamlines corresponding to model (6) in Table 2 [see fig. 1a in Contopoulos (1995) for a comparison]. Left: Zoom at the base (MSP and AP black circles, starting from the bottom, respectively). Right: MFPs. [A colour version of this figure is available in the online version.]

Each solution is characterized by the 10 parameters described in Table 1. In particular, F and Γ describe the magnetic field configuration and the kind of plasma in the system; therefore, they can be considered as defined for a given class of outflows. The magnetization parameter, σ_M, gives the efficiency with which matter is pulled out of the rotating plasma at the base of the jet. As shown in Table 1, it is a function of the angular velocity of the streamlines, the magnetic field distribution and mass load. Hence, for a jet to be efficiently launched, it has to be strongly magnetized, rapidly rotating and/or have little mass load, as pointed out by Fendt & Ouyed (2004). The σ_M parameter could be used as a fitted parameter, like we did to reproduce PMM14's results, but in general we keep this fixed, assuming again that this parameter will be defined by specific applications. Following a similar logic, we usually fix the position and the inclination of the field line, and the radial distance from the axis of the AP with θ_A, ψ_A and ϖ_A. This last parameter is related to the strength of the gravitational potential as described in Appendix A. We will discuss in detail the role of the fixed parameters in Section 5. In principle, the rest of the parameters should follow from the integration of the system of equations (1)–(4) once a suitable set of initial conditions is given for M², G² and θ (see VK03).

The integration is not simple, however, because the initial location of the MSP and MFP is not known. Therefore, some of the unknown parameters must be fixed, while the remaining ones are being retrieved as part of the solution process. The constraint we have on the solution, i.e. that all three singular points must be crossed, determines how many parameters will be found and how many need to be fixed a priori. The procedure that we follow is then an iterative one: (1) we fix a set of parameters (F, Γ, σ_M, ϖ_A, θ_A ψ_A), (2) we make an educated guess for the remaining ones (θ_MSP, θ_MFP, q, |$x_{\rm A}^2$|⁠), (3) we derive initial conditions for the integration and integrate the equations, and finally (4) we evaluate the ‘goodness’ χ of the solution, improve the guesses and integrate again until a high enough χ is achieved. We now justify and describe in more detail our approach for each of these steps.

3.1 Fixed and guessed parameters

As a basic requirement we impose that all our solutions cross smoothly all the three critical points: MSP, AP, MFP (see Fig. 2). At the AP analytical conditions are known (see VTST00). Given all the parameters,³ these conditions allow the determination of |$G^2_{\rm A}$|⁠, |$M^2_{\rm A}$|⁠, dG²/dθ|_A, and dM²/dθ|_A through the ARC (see Appendix A), thus allowing integration away from the AP.

The regularity conditions at the modified magnetosonic points MFP and MSP are not known analytically. They are both of the same kind: the denominator |$\mathcal {D}$| and the the numerators |$\mathcal {N}_1$| and |$\mathcal {N}_2$| of equations (1) and (5) must vanish, while dM²/dθ and dψ/dθ remain finite. Other authors, like VTST00, VK03, PMM13, PMM14, found their solutions by using a shooting method to integrate from the AP upstream towards the MSP and downstream towards the MFP. This method has major caveats, however, since it is very hard to numerically integrate towards singular points, in which the numerators and the denominator approach zero simultaneously, while keeping the accuracy of the solution. It is however numerically stable to integrate away from a singularity. This inspired us to explore a different approach.

We first guess the locations of the critical points, θ_MFP and θ_MSP and derive values for M², G² and their derivatives with respect to θ based on the regularity conditions |$\mathcal {N}_1= \mathcal {N}_2= \mathcal {D}= 0$| evaluated at the θ_MFP/MSP of choice (see Appendix C for our numerical technique). We then are able to integrate away from each initial guess for the modified magnetosonic points and avoid numerical inaccuracies. At the same time, we integrate away from the AP towards both magnetosonic points and consider how good the match is for the values of M² and G² of the various solutions at the mid-points θ_mid, MFP = (θ_A + θ_MFP)/2 and θ_mid, MSP = (θ_A + θ_MSP)/2 (see Fig. 3). These are in total four conditions, which imply four free parameters. θ_MFP and θ_MSP are two necessary ones and we are left with the freedom to choose two other parameters. We chose to leave |$x_{\rm A}^2$| and q free and fix the others. This choice is the most natural and convenient one: |$x_{\rm A}^2$| immediately determines |$M^2_{\rm A}$| as per equation A15, while knowing q and the position of the AP allows our algorithm to derive dM²/dθ|_A very quickly.

Figure 3.

Open in new tabDownload slide

Schematic of the method to derive the form of the functions M²(θ) or G²(θ). The solid black lines are the branches of integration with the corresponding direction shown with arrows. A typical situation where the given set of input parameters plus the guesses for the θ_MFP, θ_MSP, |$x_{\rm A}^2$| and q at these points do not provide a smooth solution through all the singular points and we see large offsets at the mid-points (red solid arrowed lines), while the closest good solution looks like the red dotted line. Note that the last integration branches towards the disc and downstream of MFP are not used in the evaluation of the fitness, but calculated at later times. [A colour version of this figure is available in the online version.]

3.2 Initial conditions at the singular points

As mentioned above, the AP is completely determined. As for the two modified magnetosonic points, there is no analytical condition that can be used to regularize the equations. We use a combination of root-finding techniques that allow us to find the values of M² and G² that give |$\mathcal {N}_1=\mathcal {N}_2=\mathcal {D}=0$|⁠. A similar procedure can be applied to both the MFP and the MSP. Once the values of M² and G² at the singular points are found, the last step before starting the integration is finding the derivative of M². By making use of equation (2) and of M² and G² at the points of interest, we determine dM²/dθ by finding the root of the function C5. The integration from all the three singular points can now start. This step is very important, serving as a numerical regularity condition, but it is somewhat laborious, so that for a more detailed discussion, we refer the interested reader to the Appendix C.

At this stage, we have both initial values and derivatives for M² and G², evaluated at θ_MFP/MSP, and we can finally start the integration inwards towards the AP. Since we integrate away from singular points, a standard adaptive step Runge–Kutta scheme is sufficient to integrate equations (1) and (2) (giving M² and G²), while we retrieve ψ from equation (3) and ξ from equation (4). Thanks to this procedure it is also possible to integrate downstream from the MFP further away from the BH and upstream from the MSP towards the equatorial plane of the disc: we simply repeat the same procedure on the other side of the critical points.

3.3 Goodness of solution and solution finding

In order to find the best parameters (θ_MSP, θ_MFP, q, |$x_{\rm A}^2$|⁠) for a given problem set (F, Γ, σ_M, ϖ_A, θ_A ψ_A) that minimize the offsets at each mid-point simultaneously (see Fig. 3), we make use of the open-source Bayesian inference algorithm multinest (Feroz & Hobson 2008; Feroz, Hobson & Bridges 2009; Feroz et al. 2013). multinest is a very robust software package which is also fast due to MPI parallelization. It also has the advantage of not requiring derivatives of the fitted function, like the Newton-Raphson method for example. This makes it even better suited for our case, because the derivatives can only be calculated numerically and such an approach would increase the numerical error and would prevent most attempts at finding a good fit. We only use the fitting algorithm of multinest. In this case, the only information that needs to be given to the algorithm is the goodness of fit of a given calculation. It then proceeds to maximize such a function by exploring the given domain of the free parameters.⁴

We accept a solution and stop the iterations when we reach a fitness value⁵ of χ ≥ 10⁹. During the process of finding the solutions, we noted that there is a highly non-linear relation between the amount of change of the different parameters q, |$x_{\rm A}^2$|⁠, θ_MFP and θ_MSP and the resulting change in the fitness function equation (12), such that each parameter should be known with at least six significant digits to make a good fit.

4 COMPARISON WITH PREVIOUS SOLUTIONS

We start our parameter study by recovering the solutions presented in table 1 of PMM14. We list our parameters corresponding to PMM14's reference and first solutions in Table 2, as models I and II. We identify two factors that can explain the discrepancies in the parameters: (1) differences in the numerical scheme and (2) definition of the functions of the gravitational potential in the gravity terms (equation A13 and Appendix B). Comparing then the parameters published in PMM14 (first solution) with model IIa, both listed in Table 2, we noticed that whilst there is not an appreciable difference in q and σ_M, which are the fitted parameters for PMM14, the locations of both MFP and MSP change by a few  per cent.

As discussed in Section 2, if we remove the approximation of a small potential in the gravity terms, we find substantial differences in our solution parameters (see models Ia–Ib, Id–Ie and IIa–IIb), which are clearly seen when drawing the corresponding streamlines: in Fig. 4, we plot the projected streamlines in the zϖ-plane for the models Ia–If. Model Ia is our reference solution, which differs from PMM14 only for the numerical method used to solve the equations (1)–(4) (black solid line). Model Ib is the reference solution with the corrected potential functions (purple dashed line) and model Ic is the reference solution with the corrected potential functions, but with Newton potential (green dotted line). It is evident that introducing the approximation of small P_g in the derivatives of P_g itself (see discussion in Section 2 and Appendix B) reduces the effect of gravity close to the BH (the right-hand panel of Fig. 4). Larger discrepancies at MSP appear when we consider a relativistic gas with adiabatic index 4/3 [see models Id (blue dot–dashed line), Ie (orange dot–dot–dashed line) and If (yellow long-dashed line)]. Generally, the effect of a small Paczyńsky–Wiita potential is very similar to a Newtonian potential. We define the last recollimation point (LRP) as the last point of the streamline, where the integration downstream of MFP stops. In the left-hand panel of Fig. 4, we note that the positions of the MFP and the LRP lie far apart from each other when comparing solutions with corrected or approximated functions of the gravitational potential.

Figure 4.

Open in new tabDownload slide

Comparison of streamlines of models Ia–If from Table 2. Colours correspond to the same model in both panels: black solid line is model Ia, purple short-dashed line is model Ib, dotted-green line is model Ic, dot–dashed blue line is model Id, orange dot–dot–dashed line is model Ie, and yellow long-dashed line is model If. The models Ia–Ic have adiabatic index of 5/3 and differ for the potential, P_g, used (N or PW) and whether the functions f_i(P_g) are approximated (A) or corrected (C). The models Id–If are the same but with Γ = 4/3. Left: MFP region, where the black dots mark the height of the MFP for the corresponding streamline. Right: Zoom on the AP and MSP. The AP (red dot) is fixed in all the solutions; therefore, the streamlines converge to the same point with the same derivative. The black dots define the location of the MSP of each streamlines. [A colour version of this figure is available in the online version.]

4.1 The self-similarity assumption

As pointed out by PMM14, the self-similarity assumption is a serious limitation intrinsic in the derivation of the equations. The inclusion of gravity further complicates the matter. Radial self-similarity, even without the inclusion of gravity, introduces a few geometrical constraints. Ultimately, to properly quantify these issues, our solutions will be benchmarked against GRMHD simulations in future works. In this section, we focus on determining the effect of gravity on the self-similarity assumption. In Fig. 5, we show two families of solutions obtained varying ϖ_A from the reference models Ia using equation (A14) and Ib (Table 2) and equation (A13). All the solutions found with the corrected functions of the gravitational potential (equation A13) can be integrated down to the disc mid-plane (z = 0). Each field line will have a characteristic angular velocity, Ω(ϖ, z), which is a constant of motion along the field line, but it will vary from one field line to another (see Appendix D). As expected from relativistic self-similar models, e.g. Li et al. (1992) and VK03, Ω(ϖ, z = 0) follows a profile which is ∝ R⁻¹. However, the gravity terms add some perturbation in the proximity of the BH. The maximum deviation between the angular velocity of a single field line at z = 0 compared to Ω_K = (R/R_g)^−3/2 (G = c = M = 1) is a factor of 2, with our field lines rotating at lower speed than the Keplerian one. All streamlines should look as scaled copies of the same shape; however, where gravity is strong enough, for example upstream of the MSP, the field lines bend and cross each others until ϖ_A is larger than ∼19r_g [panel (d) of Fig. 5]. The MFP is heavily affected as well for smaller values of ϖ_A (≲ 19). It is worth noting that for the family of solutions with model Ib as reference, we could not obtain solutions for ϖ_A < 9 due to the stronger gravitational potential, while for solutions found around model Ia, we were able to go down to ϖ_A = 5.

Figure 5.

Open in new tabDownload slide

Family of solutions obtained by varying the cylindrical radius of the AP, ϖ_A. The reference solution (ϖ_A = 15) for the upper panels is model Ia and for the bottom panels is Ib in Table 2. Upper panels: We vary ϖ_A in the interval 5–25, as the lines go from blue to red. On the left, the black dots on the curves mark the position of the MFP, while on the right they are showing the MSP. The bigger black dot at the origin of the axis is the BH. Lower panels: Same but the interval is restricted to 9–25. [A colour version of this figure is available in the online version.]

The degree in which self-similarity is affected by gravity depends on the other parameters, such as σ_M and ψ_A: we show this in Fig. 6. The four panels in the figure show the values of the fitted parameters of various families of solutions, as a function of ϖ_A. In the absence of gravity, self-similarity would ensure that |$x_{\rm A}^2$|⁠, q, θ_MFP and θ_MSP be constant within the same family of solutions, independent of any variation in the value of the cylindrical radius ϖ_A. In Fig. 6, we show families that correspond to four models of Table 2: Ia (purple), Ib (green), 1 (light blue), 2 (orange).

Figure 6.

Open in new tabDownload slide

Fitted parameters (⁠|$x_{\rm A}^2,q,\theta _{\rm MFP},\theta _{\rm MSP}$|⁠) of the same families of solutions as function of the Alfvén cylindrical radius, ϖ_A. The purple line with crosses and the green line with stars represents the families of solutions corresponding to the reference parameters of the models Ia and Ib in Table 2, respectively. The blue line with hollow squares and the orange line with filled squares are derived from the reference solutions given in Table 2 as models (1) and (2), respectively. The vertical dashed lines mark the reference ϖ_A. [A colour version of this figure is available in the online version.]

Solutions corresponding to models Ia and Ib have lower ψ_A (=40°) and lower σ_M (=0.02) with respect to model (1) and (2) that have ψ_A = 45°, σ_M = 0.5 and ψ_A = 46°, σ_M = 0.1, respectively. When we introduce gravity, we are introducing a disturbance in the radial scaling of the streamlines that becomes more pronounced as we move closer to the BH, i.e. ϖ_A ≲ 15. Gravity breaks self-similarity for all the parameters at play, but its effect is more dramatic on θ_MSP. Nonetheless, the median of the absolute deviations from the results’ median (mad = median|X_i − median(X)|) for the parameters |$x_{\rm A}^2,\ q,\ \theta _{\rm MFP}, \ \theta _{\rm MSP}$| over the full interval of ϖ_A is at most ∼ 10 per cent. This gives us the confidence in the use of this method.

5 PARAMETER SPACE STUDY

We present here a more detailed analysis of an ensemble of solutions found in a (ψ_A, σ_M) slice of the parameter space, for a specific choice of a subset of fixed parameters, namely F = 0.75, θ_A = 60°, Γ = 4/3 and ϖ_A = 15. Although we fixed it for all runs, ϖ_A has been changed when testing self-similarity for models Ia, Ib, (1) and (2) of Table 2. Changing ϖ_A changes the radius of the streamline, allowing one to check self-similarity of a given outflow with the parameters F, θ_A, σ_M and Γ. We populate a two-dimensional grid of solutions by varying ψ_A and σ_M, while fitting for |$x_{\rm A}^2,\ q,\ \theta _{\rm MFP}$| and θ_MSP and keeping all the other parameters fixed. In Table 2, models (3)–(6) belong to this grid. We use them to present features also found in other solutions in the grid. Models (7) and (8) are solutions found in other areas of the parameter space and we include them into the discussion to show their peculiar characteristics. Our primary goal in such exploration is to expand the pool of solutions found by PMM14 and search for a wider range of relativistically boosted jet solutions, i.e. with bulk Lorentz factors in the range of 2–10, in order to find suitable solutions for future applications to astrophysical objects. Considering the difficulties of starting from a random initial position, we started by changing the polytropic index to 4/3 and keeping the other fixed parameters as in the reference solution in PMM14 (the first row in Table 2), then slowly increasing σ_M and ψ_A until it was possible to find solutions. As we can see from Fig. 7, the parameter space is not a continuous volume and has patches where no solution exists. Due to the high dimensionality of this space, it is often difficult to know a priori where solutions can be found. We will focus on the nature of these boundaries later on in this section.

Figure 7.

Open in new tabDownload slide

Two-dimensional plots (⁠|$x_{\rm A}^2$| versus q in the left-hand panel and ψ_A versus |$x_{\rm A}^2$|⁠) of the parameter space we explored by making discrete steps in the parameters ψ_A and σ_M, while fitting for |$x_{\rm A}^2,q,\theta _{\rm MFP}$| and θ_MSP. The third dimension (σ_M in the left-hand panel and S/ξγ with S = −ϖΩB_ϕ/Ψc² at the MSP in the right-hand panel) is expressed with the colour scheme on the right of each panel. The two crosses in the left-hand panel are the first and reference solutions of PMM14 with corrected gravity terms. In the right-hand panel, we show which solutions of the grid are magnetically or thermally dominated at the base. The black dots are the solution that have roughly equipartition between magnetic and thermal energy at the base. On the left of these all solutions are thermally powered jets at the base, while on the right there are solutions that are Poynting dominated at the base. [A colour version of this figure is available in the online version.]

The solutions found in the (ψ_A, σ_M) slice presented here, lie on a three-dimensional curve that extends on a limited range of the interested parameters. Although solutions with σ_M < 0.05 could be recovered, we dedicated little time to the exploration of this class of solutions since they appear to have very slow bulk velocities, γ_MFP ≳ 1.

Just from the solutions in this sparse two-dimensional grid, we see already a variety of flow shapes and dynamics. For instance, we show in Fig. 8 six examples of the dynamical evolution of the energy terms (equation 6) along the jet streamlines corresponding to the solutions (3)–(8) in Table 2. The panels (a)–(c) differ in the value of σ_M, and indeed the Poynting energy at the base increases going from (a) to (c), while the enthalpy and kinetic energy remain unchanged. At the launching site, the jet can be powered by a different source of energy depending on the parameters, being first thermally dominated [0.05 ≥ σ_M < 0.2, panel (a)], crossing equipartition between magnetic and thermal energy at σ_M = 0.2 (panel (b)) and later becoming magnetically dominated [0.2 < σ_M ≤ 0.75, panel (c)]. However, we note that all these solutions have very little thermal energy overall, so they are still relatively ‘cold’ jets. When σ_M increases even further [panel (d)], the Poynting energy completely overtakes all the other energy contributions until it converts entirely into kinetic energy at ∼1000 r_g. A common feature of all our jet solutions is the conversion of the primary source of energy (enthalpy or Poynting energy) at the base into kinetic energy at some distance from the BH between the AP and the MFP, as it is to be expected for relativistic flows (Komissarov, Vlahakis & Königl 2010). The jet is always kinetically dominated by the time it approaches the MFP. There is another channel for the energy exchange, however, which we discuss below.

Figure 8.

Open in new tabDownload slide

Evolution of the energy components along the jet streamlines. The black dotted lines indicate the position of the three singular points (from left to right, MSP, AP and MFP). The black solid line is the total energy μ^΄, the green line is the specific enthalpy ξγ, the light blue line is the energy carried by the magnetic field −ϖΩB_ϕ/Ψ, the purple line is the kinetic energy γ and the dashed red line is the function (1 − P_g) of the gravitational potential. As described in Section 2 (equation 6), each component is multiplied by this function to account for gravity in the energy balance. Panel (a) shows the energy balance for model (3) with σ_M = 0.05, F = 0.75, panel (b) is model (4) with σ_M = 0.20, F = 0.75, panel (c) is model (5) with σ_M = 0.75, F = 0.75, panel (d) is model (6) with σ_M = 1.45, F = 0.75, panel (e) is model (7) with σ_M = 1.45, F = 0.85 and panel ( f) is model (8) with σ_M = 0.05, F = 0.75. [A colour version of this figure is available in the online version.]

5.1 Counter-rotation in hot jets

Panel (e) of Fig. 8 shows the energy components of model (7). This jet is roughly at equipartition at its foot point, with both enthalpy and Poynting energy being large. This is a hot, magnetized jet. We see that, starting from the MSP, the Poynting energy increases as a consequence of the transfer of a fraction of the thermal energy into the magnetic field, while the kinetic energy increases at lower pace. The fraction of energy transferred between the different components is regulated by the conservation of the total energy μ^΄. This third channel of energy transfer from the kinetic and thermal components to the magnetic component has never been seen before in a semi-analytical model, although present in simulations (see model B2H and section 5.5 in Komissarov et al. 2009) and discussed analytically by Sauty et al. (2012) and Cayatte et al. (2014). This result is thus important because it demonstrates that our semi-analytical framework produces the full range of flows seen also in MHD simulations, and is not limited to the simplest scenarios. The increase of the Poynting energy also results in an increase in the magnetic component of the angular momentum |$L^{\prime }_{\rm M}$| in equation (7) that corresponds to |$L^{\prime }_{\rm HD}$| becoming negative (see the right-hand panel of Fig. 9 and fig. 2 of Cayatte et al. 2014). The change in sign of |$L^{\prime }_{\rm HD}$| is due to V_ϕ becoming negative (see the left-hand panel of Fig. 9). In the cold regime, we normally see both components of the velocity, V_p and V_ϕ, always above zero. The jet starts off with a larger toroidal velocity that then decreases in correspondence to the poloidal component taking the lead. In the case of hot magnetized solutions, the toroidal velocity can become negative before the canonical behaviour of a cold jet is restored at larger distances from the BH. This inversion of sign of the toroidal component of the velocity is interpreted as a counter-rotation of the jet with respect to the disc. Sauty et al. (2012) and Cayatte et al. (2014) claim that the counter-rotation in jets is the signature of the magnetization of the jet and it can be due to several effects: deceleration of the flow, steep gradients of the magnetic field and/or energy transfer from enthalpy to the magnetic field. Model (7) belongs to the latter case. In Fig. 9, we start from model (7) and progressively decrease F. In the left-hand panel, we see how for higher F the solutions become hotter and their toroidal velocity minima decrease and they exhibit counter-rotation.

Figure 9.

Open in new tabDownload slide

Left-hand panel: Toroidal and poloidal velocities in units of c for a series of solutions obtained for increasing F. The solution of this series with higher F is model (7) in Table 2. Right-hand panel: Angular momentum and energy components defined in equations 8 and 9 for model (7). The vertical solid lines are the MSP, AP and MFP of this solution, while the dashed line marks the height of the jet where |$L^{\prime }_{\rm HD}$| becomes negative. The shaded grey areas are shown as a comparison with the range considered in fig. 2 of Cayatte et al. (2014). [A colour version of this figure is available in the online version.]

Similar interplay between energy terms can also result in almost oscillatory behaviours in the evolution of the energy components, seen for instance when the initial amount of enthalpy in the system is close to its minimum (=1, i.e. q ∼ 10⁻⁵–10⁻⁹) and σ_M < 1 [panel (f) of Fig. 8, depicting solution (8) of Table 2]. Here, we see two episodes where the energy is transferred to the magnetic field from the thermal energy and the flow is decelerating. However, the toroidal component of the velocity does not reverse sign, therefore no counter-rotation is established. This behaviour is usually associated with a highly wound-up streamline below MSP. The radial distance of the streamline during this phase is not constant, inducing a small change in the gravitational potential. Usually, after an initial acceleration powered by the leading force, either thermal or magnetic, the flow starts decelerating while approaching the jet axis, then the gravitational potential becomes relevant again, leading to an acceleration of the flow. Such interplay of forces can happen a few times and it results in complete winding up of the field line before the jet gains enough poloidal velocity to be launched outwards after crossing MSP. If both q and σ_M are small, upstream of the MSP the jet has little poloidal component in its velocity; therefore, it keeps slowly spiralling upwards, until it gains enough poloidal velocity to be slung out of the disc. Solutions with lower ψ_A could not be found because the jet will not have enough initial energy to acquire sufficient poloidal velocity to then be launched out.

5.2 On the effect of σ_M and ψ_A

We show two examples of the evolution of the toroidal and poloidal components of the velocity for two sets of solutions in Fig. 10 obtained by varying σ_M and ψ_A. We note that at small z (upstream of the MSP) as σ_M increases (the left-hand panel of Fig. 10) the toroidal component of the velocity increases, while the poloidal component first becomes larger and then it starts decreasing. For increasing z the toroidal component decreases towards the MFP, but, as σ_M increases, a maximum appears around the MSP. The poloidal velocity increases with distance from the BH and accelerates, rapidly approaching the speed of light close to the MFP/LRP. Indeed, we cannot find solutions with a larger magnetization parameter for the given set of input parameters (ψ_A = 46°, θ_A = 60°, F = 0.75, Γ = 4/3 and ϖ_A = 15). In the right-hand panel of Fig. 10, we show a similar plot of the velocity components for a series of solutions obtained by increasing ψ_A and fixing σ_M = 0.55, θ_A = 60°, F = 0.75, Γ = 4/3 and ϖ_A = 15. We see that varying either ψ_A or θ_A individually restricts the search to a small range in such parameters. From the ARC (see equation A20 of Appendix A), we can put a constraint on the sum of these two angles in order to select solutions that have a negative derivative of the poloidal Mach number, dM²/dθ (equation 1), at the AP (90° < θ_A + ψ_A < 180°). This interval for the sum of θ_A and ψ_A ensures that the flow is accelerating while moving away from the BH. From the evolution of the velocity components with respect to this sum, we see that the range where we can find solutions, for a given set of fixed parameters, is much smaller than the one inferred from the sign of dM²/dθ at Alfvén and that the solution can be radically different between the lower end and the upper end of the range. We also note that the larger ψ_A becomes, the closer the MFP moves towards the BH. We cannot find solutions for ψ_A > 47° probably because the LRP happens before the MFP, while our method is focused on finding solutions that pass through all the three singular points. When ψ_A decreases, we see that suddenly the jet has a much larger toroidal velocity, while the poloidal component is almost zero. The last solution is therefore the one for which the jet can still be launched from beyond MSP, while it keeps circulating upstream of this point.

Figure 10.

Open in new tabDownload slide

Left-hand panel: Toroidal and poloidal velocities in units of c for a series of solutions obtained increasing σ_M at fixed ψ_A = 46°. Right-hand panel: Same figure as the left-hand panel, but for increasing ψ_A at fixed σ_M = 0.55. The other fixed parameters for both plots are θ_A = 60°, F = 0.75, Γ = 4/3 and ϖ_A = 15. [A colour version of this figure is available in the online version.]

5.3 Exploring other regions of the parameter space

We already presented here a few solutions that were found outside the initial two-dimensional grid found by varying ψ_A and σ_M. Driven by the need to better understand which parameters lead to a shift in the observable values, such as the Lorentz factor of the jet at the MFP/collimation region, the height of the MFP, the energy balance at the launching site, we chose to try to explore different regions of the parameter space, i.e. varying previously fixed parameters. In the left-hand panel of Fig. 11, we show some of the alternative directions that we pursued. We note that each of the parameters used has a substantial effect in determining the height of the MFP compared to the relatively small steps we adopted in this search. Also, the variation of some of them (particularly σ_M) result in large changes in the bulk Lorentz factor of the jet at MFP, e.g. a Δ σ_M ∼ 1.2 around σ_M = 0.70 corresponds to Δ γ_MFP ∼ 4.

Figure 11.

Open in new tabDownload slide

Left-hand panel: Distribution of the Lorentz factors at MFP of all solutions, as a function of the height of the MFP. The blue crosses are solutions belonging to the two-dimensional search varying ψ_A for given values of σ_M. The other symbols are for solutions found by varying θ_A (orange hollow circles), σ_M (black hollow downward triangles), F (red hollow triangles) and PMM14 with correct gravity terms (yellow hollow diamonds). The arrows mark the direction of increasing values of the given parameter. Right-hand panel: Distribution of the same collection of solutions shown in the left-hand panel with respect to cylindrical radius versus height of the MFP (black crosses) and of the LRP (red crosses). The shaded light grey areas correspond to the radial coordinates of the jet at the jet break with error bars reported by Russell et al. (2014) for the XRB MAXI J1836-194. Solutions are present within each region where the jet break has been seen (grey areas). The areas with the oblique-line pattern denotes the predicted heights for the jet break. The dark grey rectangle is the area where the knot HST-1 in the jet of M87 is estimated to be by Asada & Nakamura (2012). The authors pointed out that the jet cross-section at HST-1 is smaller than the one predicted by conical and parabolic jet models. We note that HST-1 could be tracing an intermediate location between the MFP and LRP, where the jet is rapidly collimating towards the axis. [A colour version of this figure is available in the online version.]

We find that increasing F, thereby changing the scaling of the magnetic field (B ∝ r^F − 2), is the most efficient way of moving into a different area of the parameter space, letting us touch terminal Lorentz factors of about 11. However, as discussed e.g. in Blandford & Payne (1982), Contopoulos (1995) and VK03, the higher F the more the MFP is moving to larger distances from the BH, eventually to infinity. Although in principle there is no constraint in our algorithm, with the exception of loss of numerical accuracy, on how high the MFP can lie with respect to the BH, we indeed have not yet found solutions with F ≳ 0.9 for the parameters explored so far due to the aforementioned numerical accuracy issues; therefore, all our solutions are in the so-called return-current regime, i.e. the current decreases with radius.

As a further general feature shared by all solutions, we point out that they are suddenly terminated soon after the MFP, while rapidly recollimating towards the jet axis (see also Fig. 2). From a numerical point of view that happens because the denominator of the wind equation (equation 1) and the equation for ψ (equation 5) goes to zero, D → 0. We see that the streamlines are rapidly becoming vertical (ψ → π/2). At the moment, it is difficult to say whether this is due to physical effects such as compression of the gas or just due to the ‘polar axis singularity’ imposed by the self-similarity assumption, which makes the equations degenerate for θ → 0.

6 COMPARISON WITH OBSERVATIONS

The goal of the parameter search presented here has been primarily to populate an initial portion of the parameter space where solutions that are good candidates for the application to real sources reside. As shown in Fig. 11, we retrieve solutions that have bulk Lorentz factors at the MFP ranging from 1to11 and MFP's height that spans 4 orders of magnitude (∼10³–10⁷r_g). Moreover, we find jets with different initial conditions exhibiting a large range of magnetic-to-thermal energy ratios and different degrees of winding at the base. The set of initial conditions, the positions of the singular points, velocity, density and magnetic field profiles can be either directly compared to observations or can be used to provide constraints for the calculation of the emission and the polarization degree expected from a particular jet configuration.

For instance, the MFP is the location where the jet streamlines start to recollimate and the flow downstream of this region loses causal contact with the upstream flow. These are favourable conditions for the onset of a shock, while maintaining the structure upstream intact (PMM10). If, furthermore, we assume that there is a correspondence between the MFP and the frequency of the jet break, the self-absorption turnover in the synchrotron spectrum from the region where particle acceleration initiates in the jet (Markoff, Falcke & Fender 2001; Markoff et al. 2005; Markoff 2010), we can use observational constraints on the position of the jet break to determine the best jet solution within the set that we have.

In AGN, the equivalent jet break can be at large offsets from the BH (>10⁴r_g) and there could be substantial contamination in the determination of the jet break by the interaction of the jet with the environment (Russell et al. 2015) or the scenario could be further complicated by a multiple-flow dynamics (see Meier 2003 for a theoretical motivation of structured jets and e.g. Giroletti et al. 2004; Harris & Krawczynski 2006; Meyer et al. 2013 for observational evidence on the radial structure of the jets of M87). High-resolution observations can now resolve the details of the jet structure in AGN down to a few gravitational radii. Radio interferometry and high-resolution optical and X-ray observations unveiled bright subfeatures, usually referred to as knots, which are identified with shocks and particle acceleration. Knots can occur everywhere along the jet. Some of them are attributed to internal shocks caused, for example, by the interactions of the jet with a dense external cloud (Mendoza & Longair 2001), variations in the injected mass and velocity (Rees & Meszaros 1994; Malzac 2014) or magnetic recollimation (Markoff 2010). However, knots that occur at larger distances from the BH can also be produced by self-collimation of the streamlines (PMM14) or a change in the ISM pressure profile (Nakamura & Asada 2013).

The jet of M87, for example, is one of the main targets for this type of study. It exhibits a complex pattern of knots with different values and orientations of the proper motion. In particular, the HST-1 knot has received significant interest over time. HST-1 is located approximately at 2 − 5 × 10⁵r_g, which also corresponds to the region where the Bondi radius resides (r_B ∼ 3.5 × 10⁵r_g for M_BH = 6 × 10⁹ M_⊙). The Bondi radius marks the volume where the gravitational potential of the BH dominates over the thermal energy of the gas contained within this volume. The density profile within and outside the Bondi radius can be substantially different, hence the jet could experience a steep gradient in the external pressure which could induce a recollimation shock (for observational evidence, see Russell et al. 2015, and for recent simulations, see Barniol Duran, Tchekhovskoy & Giannios 2017).

Asada & Nakamura (2012) and Asada et al. (2014) show that the M87 jet maintains a parabolic profile up to HST-1. The authors show that this is evidence of the extension of the acceleration and collimation region of the jet of M87 up to the end of the sphere of influence of the BH.

If we assume then that the HST-1 knot in the jet of M87 coincides with a recollimation shock and acceleration region in the jet produced by self-collimation, we can identify it with the MFP or LRP. In all the solutions presented in this paper, we see that the jet continues recollimating downstream of the MFP and we possibly see conditions that could lead to a shock due to the compression of the gas in correspondence of the LRP. The flow could become causally disconnected but still remain smooth through the MFP, then undergoing a shock only when it is overcompressed. The distance between the MFP and LRP is relatively large in AGN so that it can be resolved by current high-resolution observations.

Using the interval given by Asada & Nakamura (2012) for the width of the jet at HST-1 (dark grey area in Fig. 11, right-hand panel), we see that the estimated position of HST-1 falls in between the MFP and the LRP of the same jet solutions. HST-1 has a smaller cross-section with respect to a canonical paraboloidal jet, probably indicating that the streamlines are converging towards the jet axis and a shock could take place soon after the MFP.

Although the identification of HST-1 with a self-collimation shock is tempting, also the change in the external pressure could influence the collimation of the jet and the shock formation due to the proximity of the Bondi radius. There is also increasing evidence that a first recollimation occurs within the first few tens of r_g (Prieto et al. 2016; Hada et al. 2016), suggesting a more complex scenario.

In XRBs, the position of the jet break has been identified in a few cases in the range of of 10–1000 r_g. The turnover has been directly observed only in three cases (GX339-4 by Corbel & Fender 2002, MAXI J1836-194 by Russell et al. 2014 and 4U 0614+091 by Migliari et al. 2006). Russell et al. (2014) show that the jet break shifts by ∼3 orders of magnitude in frequency during state transition in the source MAXI J1836-194. The estimated width of the acceleration region also presents large variations during the transition of the source, decreasing as the jet break moves closer to the BH. With the solutions that we have collected so far, we cover the full range of the jet widths seen in MAXI J1836-194 (see the right-hand panel of Fig. 11). Solutions found within the grey areas in Fig. 11 not only can reproduce the inferred widths of the jet at the location of the shock, but they give a full set of predicted physical properties of the corresponding jet, such as the height of the shock (identified with the MFP), the bulk velocity of the flow and the magnetic field morphology at any point of the jet.

Empirically constraining some of the jet properties, such as the Lorentz factor of the flow, indirectly from observations is very difficult. For instance, XRBs are thought to have mildly relativistic jets with 2 < γ_j < 5, while the Lorentz factor of AGN could span a larger range. However, recent works have shown that the range of indirectly inferred Lorentz factors for XRBs could be biased towards the lower end (Fender 2003; Miller-Jones, Fender & Nakar 2006) and in fact be the same of AGN (γ_j ∼ 1–50). With the solutions that we present here, the velocity profile is determined as part of a jet solution, and can be constrained together with other observables.

With the increasing attention on polarization studies and the growing number of current and upcoming devoted facilities, particularly the first X-ray polarimeter IXPE (Imaging X-ray Polarimetry Explorer; Weisskopf et al. 2016), it is crucial to properly model the jet magnetic field and the emission produced by such configuration of fields and particles. Indeed, as shown by various works, e.g. Perlman et al. (1999), Russell & Shahbaz (2014) and Avachat et al. (2016), multiwavelength polarization studies can help constrain the magnetic field morphology over large scales. The parameter degeneracy in the evaluation of the spectral energy distribution of the given source is also greatly reduced by polarization measurements. The need of a more detailed MHD treatment of the modelling of jets, which could be used to infer their emission is therefore compelling.

7 SUMMARY AND CONCLUSION

We presented here a new numerical scheme to solve the equations for stationary, axisymmetric, radially self-similar, relativistic MHD jets. Such class of jet models have been studied extensively over the last three decades and adopted in numerical studies and simulations of the acceleration/collimation region of jets with various applications, from GRBs to YSOs. In the case of semi-analytical models of PMM10, PMM13, PMM14, the numerical approach adopted for the integration of the equations describing the outflow was heavily affecting the efficiency of the parameter space search. Moreover, in the most recent development, the gravity terms in PMM14 were calculated with excessive approximations, which led to extra inaccuracies in the results. We avoid inconsistent approximations and also use the correct derivatives of the gravitational potential with respect to θ: in this way, we use fully self-consistent gravity terms. With this new set-up, we are able to recover and compare the solutions found in PMM14 as a test for our algorithm. The consistency between the equations affects the parameters of the solutions and, consequently, the observables that we derive from them. Once the extra approximations are removed, the contribution of gravity is stronger and this leads to a more restricted range of radii within which self-similarity holds for a given set of outflow parameters.

We were able to find solutions corresponding to single streamlines in large parts of the parameter space. The multiple jet configurations retrieved from the initial search presented in this paper are shown to be diverse in geometrical and dynamical properties, such as magnetic-to-thermal energy ratio at the jet foot point, bulk velocity, morphology of the streamlines and the position of the MFP. However, the parameter space is not continuous and the exploration along a specific direction, i.e. varying the value of a given parameter, can be interrupted due to the approach of physical limits or the break-down of our assumptions. Solutions can cease to exist because of the jet velocity approaching the speed of light around the MFP or due to the impossibility of launching the jet because of a too small energy reservoir at the base. In the latter case, the flow keeps circulating around the jet axis and never gains enough poloidal velocity to be slung out. Alternatively, solutions cannot be found when the LRP takes place before the MFP, while our method is constructed so to cross all three singular points.

By studying the evolution of the energy and the velocity components for all solutions, we encounter both cold and hot relativistic jets, the latter configurations exhibiting counter-rotation at distances ∼0.1–10⁵ r_g. With the exploration of the parameter space that we have conducted, we are able to cover a large area of the parameter space and to move towards a desired direction, e.g. higher bulk Lorentz factors, by making use of the trends that we observed so far.

Taking into account the large variety of jet configurations and dynamics that we can obtain with our numerical scheme, the extension to other accreting objects, e.g. YSOs, is within the capabilities of this model. Future steps will be taken to increase the density of the solution grid and cover larger volumes of the parameter space. This mapping of the parameter space is ultimately aimed at providing a dynamic and flexible base for future coupling with a radiative code. That, in turn, will eventually allow data-fitting of both XRBs and AGN, which will be the subject of forthcoming papers. We are also aware of the need to benchmark the results of our method to fully numerical simulations to properly quantify the limits of our approach.

Acknowledgements

CC is supported by the Netherlands Organisation for Scientific Research (NWO) grant no. 614.001.209. SM acknowledges support from NWO VICI grant no. 639.043.513. YC is supported by the European Union's Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie Global Fellowship grant agreement no. 703916.

REFERENCES

APPENDIX A: DEFINITIONS OF THE EQUATIONS

We present here the system of equations (1)–(4) in its full form as given in PMM14 and we report the changes in the function of the gravitational pseudo-potential.

APPENDIX B: DERIVATION OF THE PSEUDO-POTENTIAL FUNCTIONS IN THE GRAVITY TERMS

APPENDIX C: INITIAL CONDITIONS

Once we find the initial conditions for M² and G² and their derivatives at the three critical points, we perform the integration from each point with an adaptive stepsize Runge–Kutta scheme as described in Press et al. (1993, hereafter NR93), until we reach the mid-points θ_mid, MFP and θ_mid, MSP. However, setting up the initial conditions is both the real challenge and the base for the robustness of our scheme.

As we mentioned in the Section 3, given the full set of parameters, the initial conditions at the AP are readily calculated. The only numerical step is finding dM²/dθ. The derivative of M² at Alfvén is commonly referred to as p_A and can be found from the ARC as described in Polko et al. (2010, 2014) and in equation (A20).

Now the integral of motion μ^΄ can also be calculated (see equation 6) and all the other quantities can be analytically calculated at the AP.

The situation is more complicated at the modified magnetosonic points, since no analytical condition is known. Our solution is a sequence of root finding routines to find the zeroes of specific functions. The method is almost the same at both MSP and MFP; therefore, we will describe only the case of the MFP, highlighting the few differences when they exist.

In order to find the minimum of the critical function, we create a grid in the M² and G² space with boundaries dictated by the location of the AP. We then use the three points that have the lowest values of the critical function |$\mathcal {C}$| as the starting nodes of a global simplectic minimizing algorithm such as amebsa (see NR93, for a description). When the minimum has been found, we try to polish the solution with a further local powell minimizing routine (NR93). At the end, we further verify that the absolute values of all the three numbers |$\mathcal {N}_1$|⁠, |$\mathcal {N}_2$| and |$\mathcal {D}$| are close enough to zero (<10⁻¹⁰). While this procedure of polishing the first solution and then cross-checking again seems redundant and CPU consuming, it is our experience that this process avoids having spurious starting points that would contaminate the results of the next steps in our solution finding scheme.

Although we now know the initial values of M² and G², we cannot yet start the integration, since the derivative dM²/dθ is still unknown: i.e. equation (1) would still be zero over zero. In order to recover the value of the derivative, we search the roots of another function.

In most cases, we find two roots for M²(θ): one smaller and one larger than M²(θ_MFP). Based on the physical picture of the accelerating jet and that dM²/dθ|_A < 0, we assume that the derivative is negative also at the critical points. Therefore at MFP, we choose M²(θ) < M²(θ_MFP) and at MSP we take M²(θ) > M²(θ_MSP). We cannot justify this assumption based on general principles, and cannot provide a definitive proof, but we performed various tests allowing for the alternative choices, and never managed to find a solution. Now that we have the initial conditions we can integrate.

APPENDIX D: CONVERSION TO PHYSICAL QUANTITIES