© ESO 2018

1 Introduction

Active galactic nuclei (AGN) jets are now recognized as being multi-component outflows related to accretion onto a supermassive black hole (SMBH). The outer kiloparsec/megaparsec-scale lobes are fed by a powerful hadronic plasma, which most likely originates from the Keplerian accretion disk via the magnetocentrifugal launching mechanism, Blandford & Payne (1982). This hadronic population could be responsible for the second component peaking in the γ-ray band of the spectral energy distribution (SED) of some blazars. As explained in Böttcher et al. (2013), this is possible if a relativistic jet of protons contributes significantly to the radiative output through proton synchrotron emission or photo-pion production. However the purely hadronic synchrotron models for blazars present a major problem due to the requirement for having very high powers in the jets. Indeed, hadronic processes are very inefficient (Sikora 2011). They are only energetically favorable for high-frequency synchrotron peak (HSP) blazars, where the jet is likely to be highly magnetized (Petropoulou & Dermer 2016).

The alternative approach to explaining the origin of the high-energy emission is to consider leptonic models. There, the radiative output throughout the electromagnetic spectrum is assumed to be dominated by leptons both at low and high frequencies. Inverse Compton scattering of soft photons (IC) by relativistic nonthermal electrons in the jet is supposed to be the most probable mechanism for γ-ray production, at least for strong-line and low-frequency synchrotron peak – LSP – blazars. Until recently, leptonic models in quasi steady-state were very successful in modeling the SED for almost all classes of blazars (e.g., Celotti & Ghisellini 2008). However leptonic models with a one-zone component where only one portion of the jet dominates the emission are now questioned. Mostly, they fail to produce the extremely high bulk Lorentz factors in a very compact emission region required for ultra-fast variability of some TeV blazars (Begelman et al. 2008). Despite several proposed mechanisms mostly based on small-scale inhomogeneities in the jet (see references in Vovk & Babić 2015), it is not yet possible to know the location of the variable γ-ray emission; it may be at the base of the jet or up to parsec-scale distances from the central black hole. Moreover most of the sources in which ultra-rapid variability events have been observed are BL Lac objects emitting very-high-energy radiation.

For BL Lacs the application of the one-zone synchrotron-self Compton model to the SED implies that the jet is weakly magnetized and the emitting region is far from equipartition. Both conditions are required if magnetic-reconnection-powered emission is at the origin of the ultra-rapid variability (Tavecchio 2016). A most important constraint comes from the launching and the acceleration processes, which should lead to almost equipartition between the magnetic and the kinetic energy fluxes as can be seen in other blazars. In order to unify the BL Lac and radiogalaxy populations, Ghisellini et al. (2005) already proposed a structured jet model with two components; a faster core (the spine) surrounded by a slower sheath or layer. This two component model is able to explain peculiar features of TeV-emitting BL Lacs, such as the absence of fast superluminal components and the presence of a limb-brightened radio structure (Giroletti et al. 2004; Nagai et al. 2014). Tavecchio & Ghisellini (2016) showed that this structured jet model gives an extra source of soft photons intervening in the IC emission for BL Lacs. It is due to the radiative interplay between the two components and allows to reproduce the emission in equipartition conditions. This spine-sheathjet structure has also been explored in Sikora et al. (2016) for strong-line blazars in order to explain why γ-ray variations are often observed to have much larger amplitudes than the corresponding optical variations, as well as other observed features of γ-ray flares. For the nearby radiogalaxy M 87, a direct detection of different motions was possible and an extended mildly relativistic flow is seen surrounding a relativistic central jet (Mertens et al. 2016).

The possible seed photon population at the origin of IC scattering in AGNs is questionable. The dominant seed source depends critically on the location of the primary emitting region in the jet (Finke 2016). Radio array observations of individual blazars are used to localize the γ-ray emission site among the superluminal components seen along the jet. However, they cannot put strong constraints on localizing the emissionregions of the high-energy flares. More theoretical and numerical modelings and carefully analyzed polarization data are needed (see Kravchenko et al. 2016 and references therein). We note also the treatment used in Hada et al. (2011) for locating the central black hole relatively close to the jet base-radio core in M 87 using the dependence of the core position on frequency. There are cases where the spine-sheath model should work. The leptons in superluminal components after leaving the black hole corona interact with the sheath where IC scattering of photons occurs and produce the observed γ-ray activity (Marscher et al. 2010; Casadio et al. 2015).

For AGN jets, two-component models have been theoretically studied in Sol et al. (1989) and further developed by several authors (e.g., Fabian & Rees 1995). In those studies the spine jet is the inner leptonic component (electron-positron pairs), which is self collimated inside a hadronic disk wind, and should be separated from the wind by a force-free, empty region. On the other hand, it has been shown numerically that the high collimation efficiency of the outer-disk wind is sufficient to prevent decollimation from the inertia of the inner relativistic plasma, in the case where the initial magnetic field lines going out from the launching region, that is, the black hole corona, are radial (e.g., Gracia et al. 2009). Moreover synthetic synchrotron maps built from such magnetohydrodynamic (MHD) models can reproduce the opening angle of M 87 up to a large distance from the black hole, and the change between center-brightening near the core to limb-brightening further away.

The standard accretion-jet model includes an advection-dominated accretion flow (ADAF) which takes place at least in the inner part of the accretion disk in AGNs, except for very high accretion efficiency (Narayan & Yi 1994). The presence of this ADAF structure in the disk center does not prevent the launching of the internal jet as can be seen from observations near the black hole for radiogalaxies (see the Faraday rotation measures from observations of M 87 in Feng et al. 2016). In the high-accretion regime for AGNs, where a standard thin disk is found and powerful outflows are observed, both an internal jet and an outer disk wind coexist.

Numerical simulations in the framework of general relativity have been performed to model the inner jet formation assuming a large-scale magnetic field right from the beginning, mostly using monopolar configuration anchored in the ergosphere. Komissarov (2007) was able to ensure the viability of the Blandford–Znajek mechanism (Blandford & Znajek 1977) and the possibility of building a relativistic particle outflow. Similar relativistic simulations of jets were obtained by McKinney & Blandford (2009) where the central Poynting-flux-dominated force-free jet is self-consistently feeding the funnel of the accretion-disk wind created along the axis by the very high toroidal field in a magnetically arrested disk (MAD). McKinney et al. (2012) and Tchekhovskoy et al. (2011) have also derived a scaling law between the magnetic flux of the field threading the black hole and the mass-accretion rate on the horizon. This has been confirmed observationally by Zamaninasab et al. (2014). In addition, the jet power, directly related to the magnetic flux on the horizon, was calculated in their simulations, allowing for the derivation of a net flow efficiency between the accretion and the jet.

Beyond the launching region, the internal structure of magnetized and relativistic jets can be probed by relativistic MHD (RMHD) simulations, as has been done by Martí et al. (2016). They characterized the internal jet structure of overpressured, steady jets in connection with the dominant energy (internal, rest-mass, and magnetic). They showed that transverse equilibrium with a significant toroidal magnetic component implies a structure with a central spine and a surrounding layer with lower thermal and total pressures. Thus, we have strong clues, both from observations and simulations, that relativistic jets have a transverse structure. On the other hand, the jet launching mechanism associated with magnetic flux threading the black hole seems to produce only very light spine jets with rarified gas such that the inner outflow is force-free. The inner beam would be so light that it would be invisible at large distances. An alternative issue, in order to model the spine jet, is to extend the self-similar solutions for jet outflows in young stellar objects to general relativity. This model can be derived analytically and used as initial conditions in numerical simulations.

A self-similar model of the spine jet has been proposed by Meliani et al. (2006) for nonrotating black holes in Schwarzschild metric, extending the previous Newtonian self-similar model of Sauty & Tsinganos (1994). The Meliani et al. (2010) model aimed to explain the dichotomy between FRI and FRII sources in terms of magnetic collimation efficiency. Then, Globus et al. (2014), partially generalized this model to a Kerr metric. However they fail to give an exact generalization of the magnetic collimation criterion linked to the efficiency of the magnetic rotator. In both cases, the authors assume that the light cylinder is sufficiently far from the spine jet to have its effects neglected. Physically, this means that both models cannot produce solutions consistent across the light cylinder.

In the present paper, we present a new extension of the nonrelativistic meridional self-similar solutions of Sauty & Tsinganos (1994), both in Schwarzschild and Kerr metrics. This model can produce solutions for relativistic jets emerging from a spherical corona surrounding the central part of a Kerr black hole and its inner accretion disk.

Unlike previous models in a Schwarzschild or a Kerr metric, this new model for a nonrotating, as well as for a rotating black hole includes a self-consistent crossing of the light cylinder in the solution. The Alfvén Mach surfaces are not spherically symmetric. However, the Alfvén transition surface is spherical and includes the light cylinder.

In Sect. 2 we summarize the ideal MHD equations in the 3+1 formalism and in Sect. 3 construct the model and derive its equations. In Sect. 4 we describe how to solve the equations and illustrate this procedure with various new solutions in a Schwarzschild metric. Then, four different solutions in a Kerr metric are presented in Sect. 5 with high Lorentz factor and different geometries potentially applicable to AGNs and γ-ray bursts (GRB). Some characteristics of those solutions of the model are discussed as well as how the magnetic collimation evolves with the rotation of the black hole.

2 Steady axisymmetric relativistic MHD outflows

2.1 Kerr metric

The first step in building self-similar solutions in relativistic flows is to define the metric. In fact the central massive black hole dominates the gravitational field in the near regions and determines completely the metric field. Thus, in Kerr metric, the geodesics are defined as follows. $$\begin{array}{lll}\hfill & \text{d}{s}^2=\hfill & -\left(1-\frac{r_{s}r}{{\rho }^2}\right)c^2\text{d}{t}^2-\frac{2r_{s}rca}{{\rho }^2}{\sin }^2\theta \text{\hspace{0.17em}d}t\text{\hspace{0.17em}d}\varphi \hfill \\ \hfill & \hfill & +\frac{{\rho }^2}{\Delta }\text{d}{r}^2+{\rho }^2\text{d}{\theta }^2+\frac{{\Sigma }^2}{{\rho }^2}{\sin }^2\theta \text{\hspace{0.17em}d}{\varphi }^2.\hfill \end{array}$$(1)

We have the usual notations of the elements, $$\begin{array}{l}\Delta =r^2+a^2-r_{s}r,\\ {\rho }^2=r^2+a^2{\cos }^2\theta ,\\ {\Sigma }^2={(r^2+a^2)}^2-a^2\Delta {\sin }^2\theta ,\end{array}$$

where $a=\frac{\mathcal{J}}{\mathcal{M}c}$ and $r_s=\frac{2\mathcal{G}\mathcal{M}}{c^2}$.

We note that $\mathcal{J}$ is the angular momentum of the massive central object, $\mathcal{M}$ is its mass, h is the lapse function, ω is the angular velocity of zero angular momentum observers (ZAMO) and we use a for the length-scale related to the angular momentum of the black hole (Kerr scale). We can define the dimensionless spin of the black hole a_H in units of the gravitational radius r_s∕2 such that, a_H = 2a∕r_s. Furthermore, β is the shift vector.

The lapse function h, the angular velocity ω of zero angular momentum observers (ZAMO) and the shift vector coordinates can be written as: $$\begin{array}{l}h={\left(1-\frac{r_{s}r}{{\rho }^2}+{\beta }^{\varphi }{\beta }_{\varphi }\right)}^{1/2}=\frac{\rho }{\Sigma }\sqrt{\Delta },\\ \omega =\frac{ac{r}_{s}r}{{\Sigma }^2},{\beta }_{\varphi }=-\frac{\omega }{c}{\varpi }^2,{\beta }^{\varphi }=-\frac{\omega }{c},\end{array}$$

with $\varpi =\frac{\Sigma }{\rho }\sin \theta$. The corresponding line elements for the Kerr metric are given in Appendix A.

2.2 Maxwell’s equations

The next step is to define the electromagnetic field in this metric. Using covariant derivatives, we can write Maxwell’s equations in Kerr space, assuming stationarity and axisymmetry, $$\begin{array}{l}\nabla \cdot E=4\pi {\rho }_e,\\ \nabla \cdot B=0,\\ \nabla \times (hE)=\left(B\cdot \nabla \frac{\omega }{c}\right)\varpi {ϵ}_{\varphi },\\ \nabla \times (hB)=\frac{4\pi h}{c}J-\left(E\cdot \nabla \frac{\omega }{c}\right)\varpi {ϵ}_{\varphi },\end{array}$$

where ${({ϵ}_i)}_{i=\mathrm{1...3}}$ is the space orthonormal basis. We note that all quantities in the above equations are given in the ZAMO frame. We can split all vector fields in a poloidal component in the meridional plane and a toroidal one along the azimuthal direction. The poloidal magnetic field B_p can be expressed in terms of the magnetic flux function A, $$B_{\text{p}}=\nabla \times \left(\frac{A}{\varpi }{ϵ}_{\varphi }\right),$$(11)

and using Faraday’s law, we get the electric field, $$\nabla \times \left(hE-\frac{\omega }{c}\nabla A\right)=0.$$(12)

We can introduce the electric potential, Φ, but the electric field E is not directly proportional to the gradient of the electric potential, $$hE=\frac{\omega }{c}\nabla A-\nabla \Phi .$$(13)

The condition of ideal MHD for infinite electrical conductivity leads to $$E+\frac{V\times B}{c}=0.$$(14)

We note that ∇ is the covariant derivative on a space hyper-surface; see Appendix B for the expression of its coordinates.

2.3 Equations of motion

In the Kerr metric, the 3+1 formalism gives the equation for mass conservation, the Euler equation, and the energy conservation, respectively, $$\begin{array}{l}\nabla \cdot ({\rho }_0\gamma hV)=0,\\ {\rho }_0\gamma \left(V\cdot \nabla \right)\left(\gamma \xi V\right)+{\rho }_0\xi {\gamma }^2\left[c^2\nabla \ln h+\frac{\varpi \omega V^{\widehat{\varphi }}}{h}\nabla \ln \omega \right]\\ +\nabla P={\rho }_{e}E+\frac{J\times B}{c},\\ {\gamma }^2{\rho }_0\xi c\left[V\cdot \nabla \ln (\gamma \xi h)+\frac{\omega \varpi V^{\widehat{\varphi }}}{h{c}^2}V\cdot \nabla \ln \omega \right]=\frac{J\cdot E}{c}.\end{array}$$

Here $V^{\widehat{\varphi }}$ is the toroidal component of the bulk flow speed as seen by the ZAMO. The factor γ is the bulk Lorentz factor, ρ₀ is the mass density, and ξc² the specific enthalpy measured in the comoving frame of the outflow, that contains kinetic enthalpy of perfect relativistic gas ξ_K and some heating term Q∕c². $$\xi ={\xi }_K+\frac{Q}{c^2}.$$(18)

For the kinetic enthalpy we use the Taub–Matthews approximation of ideal fluid equation of state; for more details see Taub (1948), Meliani et al. (2004) and Mignone et al. (2005). $${\xi }_K=\frac{5}{2}\left(\frac{P}{{\rho }_0{c}^2}\right)+\sqrt{1+{\left(\frac{3P}{2{\rho }_0{c}^2}\right)}^2}.$$(19)

The energy conservation has been derived in the frame of the ZAMO and equivalently the first law of thermodynamics can be obtained by projecting the conservation of the energy-momentum tensor along the fluid 4-velocity but in the comoving frame. Assuming infinite conductivity, the contribution of the electromagnetic field is null and only the thermal energy affects the variation of the enthalpy of the fluid, giving, $${\rho }_0\left(V_{\text{p}}\cdot \nabla \right)(\xi c^2)=\left(V_{\text{p}}\cdot \nabla \right)P.$$(20)

2.4 Constants of motion

Under the assumptions of steadiness and axisymmetry, the MHD equations in general relativity can be partially integrated to yield severalfield/streamline constants (Beskin 2010). We already deduced those constants including the magnetic field with the same formalism in Cayatte et al. (2014). Here we present the derivation of the equations, following the notations of Tsinganos (1982), in order to compare with previous self-similar models, for example, Meliani et al. (2006) and give the choice of the first integrals.

Steady and axisymmetric flows are characterized by a function A that defines the geometry of the magnetic flux surfaces. In the poloidal plane, field lines are lines of constant magnetic flux A and first integrals will be functions of A, among which the mass flux Ψ. The poloidal velocity can be expressedin terms of Ψ, $$4\pi {\rho }_0\gamma h{V}_{\text{p}}=\nabla \times \left(\frac{\Psi }{\varpi }{ϵ}_{\varphi }\right).$$(21)

The frozen-in condition for ideal MHD flows, gives, in the toroidal direction, combined with Eq. (11), $$4\pi {\rho }_0\gamma h{V}_{\text{p}}={\Psi }_A{B}_{\text{p}},$$

where Ψ_A ≡dΨ∕dA is the magnetic to mass flux ratio.

The poloidal components of the law of flux freezing (Eq. (14)) give in turn the iso-rotation law, $$\Omega -\omega =\frac{h{V}^{\widehat{\varphi }}}{\varpi }-\frac{{\Psi }_A{B}^{\widehat{\varphi }}}{4\pi {\rho }_0\gamma \varpi },$$(22)

where Ω(A) ≡ cdΦ∕dA, which is the isorotation frequency, constant along each magnetic flux tube.

By integrating the Euler equation in the toroidal direction, we get the conservation of the angular momentum flux L(A), $$L=\varpi \left(\gamma \xi V^{\widehat{\varphi }}-\frac{h{B}^{\widehat{\varphi }}}{{\Psi }_A}\right).$$(23)

The last equation to integrate is the energy conservation. In other words, we may take the Euler equation projected along the time axis of the 3+1 decomposition, and integrate it under the hypothesis of steadiness, $$\mathcal{E}-L\omega =\gamma \xi h{c}^2-\frac{h\varpi (\Omega -\omega )}{{\Psi }_A}{B}^{\widehat{\varphi }}.$$(24)

2.5 Toroidal fields

Using the three last integrals of motion, we may express the toroidal components of the velocity and the magnetic fields and the enthalpy density as functions of these first integrals and the poloidal components. Using the standard procedure of inversion we get $$\begin{array}{l}\varpi \frac{h{B}^{\widehat{\varphi }}}{{\Psi }_A}=\frac{L\left[h^2{c}^2+{\varpi }^2\omega (\Omega -\omega )\right]-\mathcal{E}{\varpi }^2(\Omega -\omega )}{\left(M_{\text{Alf}}^2-h^2\right)c^2+{\varpi }^2{(\Omega -\omega )}^2},\\ \varpi \gamma \xi V^{\widehat{\varphi }}=\frac{M_{\text{Alf}}^{2}L{c}^2-(\mathcal{E}-L\Omega ){\varpi }^2(\Omega -\omega )}{\left(M_{\text{Alf}}^2-h^2\right)c^2+{\varpi }^2{(\Omega -\omega )}^2},\\ \gamma h\xi =\frac{M_{\text{Alf}}^2(\mathcal{E}-L\omega )-h^2(\mathcal{E}-L\Omega )}{\left(M_{\text{Alf}}^2-h^2\right)c^2+{\varpi }^2{(\Omega -\omega )}^2},\end{array}$$

where we have defined the poloidal Alfvén Mach number, $$M_{\text{Alf}}^2=h^2\frac{V_{\text{p}}{}^2}{V_{\text{Alf}}^2}=\frac{4\pi h^2{\rho }_0\xi {\gamma }^2{V}_{\text{p}}{}^2}{B_{\text{p}}{}^2}=\frac{\xi {\Psi }_A{}^2}{4\pi {\rho }_0}.$$(28)

This definition of the poloidal Alfvén Mach number is consistent with the definition used by Meliani et al. (2006) and includes the lapse function. This is also the definition taken by Breitmoser & Camenzind (2000) because the velocity, hV_p, calculated with the universal time is continuous across the event horizon.

The numerator and denominator of Eq. (25) are zero at the Alfvén transition surface, if the following two equations are satisfied, $${M_{\text{Alf}}^2|}_{\text{a}}=h_{\text{a}}^2\left[1-\frac{{\varpi }_{\text{a}}^2{(\Omega -{\omega }_{\text{a}})}^2}{h_{\text{a}}^2{c}^2}\right],$$(29) $$\frac{{\varpi }_{\text{a}}^2{(\Omega -{\omega }_{\text{a}})}^2}{h_{\text{a}}^2{c}^2}=\frac{L(\Omega -{\omega }_{\text{a}})}{(\mathcal{E}-L{\omega }_{\text{a}})}.$$(30)

The denominators of Eqs. (26) and (27) are identical to the one of Eq. (25) and there numerators are a linear combination of Eqs. (29) and (30). Hence, the numerators and the denominators of Eqs. (26) and (27) are also zero at the Alfvén transition surface.

We can reformulate the above equations by changing variables. We rescale the cylindrical radius with ch∕(Ω − ω) leading to the dimensionless cylindrical radius x and introduce the parameter x_MR, $$\begin{array}{ll}x=\frac{\varpi (\Omega -\omega )}{hc},\hfill & x_{\text{MR}}^2=\frac{L(\Omega -\omega )}{(\mathcal{E}-L\omega )}.\hfill \end{array}$$(31)

Hence we can write, $$\begin{array}{l}{B}^{\widehat{\varphi }}=\frac{-(\mathcal{E}-L\omega ){\Psi }_A}{cx}\frac{x^2-x_{\text{MR}}^2}{M_{\text{Alf}}^2-h^2\left(1-x^2\right)},\\ h\gamma \xi \frac{V^{\widehat{\varphi }}}{c}=\frac{(\mathcal{E}-L\omega )}{c^{2}x}\frac{M_{\text{Alf}}^2{x}_{\text{MR}}^2-(1-x_{\text{MR}}^2)h^2{x}^2}{M_{\text{Alf}}^2-h^2\left(1-x^2\right)},\\ \gamma h\xi =\frac{(\mathcal{E}-L\omega )}{c^2}\frac{M_{\text{Alf}}^2-h^2(1-x_{\text{MR}}^2)}{M_{\text{Alf}}^2-h^2\left(1-x^2\right)}.\end{array}$$

The second condition, Eq. (30), at the Alfvén transition surface becomes, keeping the first one unchanged, $${x^2|}_{\text{a}}={x_{\text{MR}}^2|}_{\text{a}}.$$(35)

In Kerr metric, the parameter $x_{\text{MR}}^2$ is an extension of $x_{\text{A}}^2$ defined by Meliani et al. (2006). It measures the amount of energy carried by the electromagnetic field. This is the energy flux of the magnetic rotator (MR) divided by the total energy flux of the outflow in the co-rotating frame, $\mathcal{E}-L\omega$. This new parameter x_MR, conversely to the previous x_A, is not any more constant along a field line, since ω is not an integral of the motion.

3 Model equations

3.1 Angular expansion

The MHD equations and the metric constitute a coupled set of highly nonlinear equations that cannot be solved analytically. The approach followed so far for Newtonian flows has been to look for solutions with separable variables in the frame of self-similarity. However, this technique cannot be applied in the frame of general relativity due to the complexity of the metric even for the simpler cases of a Schwarzschild or a Kerr metric. Instead, we may model the jet close to its symmetry axis, that is, to describe the spine jet, by expanding all variables with sin θ to second order.

Along the polar axis where ϖ and θ go to zero, we may define the spherical Alfvén radius to be the distance r_⋆ from the center where the Alfvén transition surface condition, $M_{\text{Alf},\theta =0}^2=h_{\star }^2$, applies. The subscript ⋆ denotes the value of a physical quantity at the Alfvén transition surface, along the polar axis. We shall use this location to write all our quantities in dimensionless form. Thus, the dimensionless spherical radius is $$R=\frac{r}{r_{\star }}.$$(36)

At R = 1, the velocity is V_⋆, the magnetic field B_⋆, the density ρ_⋆, the enthalpy ξ_⋆ and the lapse function h_⋆. Because of the Alfvén transition along the polar axis, we have, $$B_{\star }{}^2=4\pi {\gamma }_{\star }{}^2{\rho }_{\star }{\xi }_{\star }{V}_{\star }{}^2.$$(37)

Thus the dimensionless magnetic flux function α is defined as, $$\alpha =\frac{2}{r_{\star }{}^2{B}_{\star }}A.$$(38)

Moreover we can expand to the second order the metric of the system in dimensionless form using the characteristic dimensions of the system defined in the previous Section. This introduces the two following new parameters, $$\mu =\frac{r_s}{r_{\star }},l=\frac{a}{r_{\star }}=\frac{\mathcal{J}}{Mc{r}_{\star }}\Rightarrow \frac{2a}{r_{\star }}=\frac{2l}{\mu },$$(39)

which are respectively the Schwarzschild radius in units of the Alfvén radius and the dimensionless black hole spin.

Another dimensionless parameter is needed to describe the gravitational potential, as in the classical model. This parameter ν represents the escape speed at the Alfvén point along the polar axis in units of V_⋆. Then, the value of V_⋆ is fixed by the following condition, $$\nu =\frac{V_{\text{esc},\star }}{V_{\star }}=\sqrt{\frac{2\mathcal{G}\mathcal{M}}{r_{\star }{V}_{\star }^2}}\Rightarrow V_{\star }^2=\frac{\mu }{{\nu }^2}{c}^2.$$(40)

Thus, tosecond order in sin θ the ZAMO angular velocity and the lapse function are written as, $$\begin{array}{lll}\omega \hfill & =\hfill & \frac{lc\mu R}{r_{\star }{(R^2+l^2)}^2}\left(1+\frac{l^2{h}_z{}^2}{R^2+l^2}{\sin }^2\theta \right),\hfill \\ h_{}\hfill & =\hfill & \sqrt{1-\frac{\mu R}{R^2+l^2}}\left(1-\frac{\mu l^{2}R}{2{(R^2+l^2)}^2}{\sin }^2\theta \right).\hfill \end{array}$$

In order to simplify our notation, we define the lapse function along the polar axis, $$h_z(R)=h(R,\theta =0)=\sqrt{1-\frac{\mu R}{R^2+l^2}},$$(43)

and the polar shift of the metric, $${\omega }_z(R)=\omega (R,\theta =0)=\frac{lc\mu R}{r_{\star }{(R^2+l^2)}^2}.$$(44)

See Appendix A for details.

It will be useful to introduce the dimensionless polar shift function (see also Eq. (59)), $${\overline{\omega }}_z(R)=\frac{{\omega }_z{r}_{\star }}{V_{\star }{h}_{\star }}=\frac{l\sqrt{\mu }\nu R}{h_{\star }{(R^2+l^2)}^2}.$$(45)

We also expand the magnetic flux function to second order in sin θ. The magnetic flux is an even function which is zero along the polar axis due to axisymmetry and because of the symmetry around the equatorial plane. Thus all odd orders are zero and the first nonvanishing even order is the second order in colatitude. If we keep the lowest order in the expansion we get $$\alpha (R,\theta )=f(R){\sin }^2\theta ,$$(46)

where f is the inverse of the classical expansion factor for solar coronal holes (see Tsinganos & Sauty 1992). This expansion, similarly to the classical self similar model of Sauty & Tsinganos (1994), is equivalent to a hypothesis of separation of the variables in the magnetic flux function.

Thus from Eq. (A.17), the cylindrical radius can also be seen as an expansion in the magnetic flux. This is physically more meaningful as the magnetic flux is constant on a given mass flux tube. Moreover, several free integrals depend solely on this magnetic flux. We define the dimensionless cylindrical radius G in units of the polar Alfvén radius as, $$G(R)=\sqrt{\frac{R^2+l^2}{f(R)}}.$$(47)

The cylindrical radius can be written in the various following forms, $${\varpi }^2=r_{\star }^2(R^2+l^2){\sin }^2\theta =r_{\star }{}^2{G}^2\alpha =G^2{\varpi }_a^2.$$(48)

We can also write the metric as an expansion in α (see also Appendix A), $$\begin{array}{lll}\omega \hfill & =\hfill & \frac{lc\mu R}{r_{\star }{(R^2+l^2)}^2}\left(1+\frac{l^2{h}_z{}^2{G}^2}{{(R^2+l^2)}^2}\alpha \right),\hfill \\ h_{}\hfill & =\hfill & \sqrt{1-\frac{\mu R}{R^2+l^2}}\left(1-\frac{\mu l^{2}R{G}^2}{2{(R^2+l^2)}^3}\alpha \right).\hfill \end{array}$$

Of course we can always reverse our point of view and go back to the expansion in θ. This would be the case if we wanted to use the steady analytical solution as initial conditions for numerical simulations.

We can parametrize the geometry of the flux tubes with the logarithm derivative of f denoted F, $$F=\frac{\text{d\hspace{0.17em}ln\hspace{0.17em}}f}{\text{d\hspace{0.17em}ln\hspace{0.17em}}R}=2\left(\frac{R^2}{R^2+l^2}-\frac{\text{d\hspace{0.17em}ln\hspace{0.17em}}G}{\text{d\hspace{0.17em}ln\hspace{0.17em}}R}\right).$$(51)

The angle χ of the magnetic poloidal field line with the radial direction (see Sauty et al. 1999) is given in our metric by, $$\tan \chi =\frac{\sqrt{R^2+l^2-\mu R}}{2R}F\tan \theta .$$(52)

3.2 Choice of the Alfvén surface and pressure

We can expand all physical quantities to the first order in α. Thus the Alfvén number is given by, $$M_{\text{Alf}}=M(R)\left(1+M_1(R)\alpha \right).$$(53)

Contrary to previous self-similar models, the Alfvén number cannot be spherically symmetric because of the presence of the cylindrical radius in units of the “light cylinder” x, in the numerator and the denominator of Eqs. (32)–(34). This is induced by the regularity conditions, Eqs. (29) and (30), and the sphericity of the Alfvén surface. The surface x = 1 is the so-called outer “light cylinder”. Of course this surface may not be exactly cylindrical if x depends also on α, which may be the case for instance close to the black hole where ω has a strong dependence on α or if Ω is not constant with α. Therefore, this is rather a light surface, but for the sake of simplicity we refer to it as a “light cylinder” in the rest of the text.

Similarly, the pressure can be expanded to first order, $$P(R,\alpha )=P_0+\frac{{\gamma }_{\star }{}^2{\rho }_0{}_{\star }{\xi }_{\star }{V}_{\star }{}^2}{2}\Pi (R)\left(1+K(R)\alpha \right),$$(54)

where P₀ is a constant.

In order to simplify and as a first step, we assume for both equations that the radial dependence of the nonpolar component of the Alfvén number and the pressure are simply constant, M₁ (R) = m₁ = cst, K(R) = κ = cst. Thus, $$M_{\text{Alf}}=M(R)\left(1+m_1\alpha \right).$$(55)

We note that we have m₁ = 0 in previous models; see Meliani et al. (2006) and Globus et al. (2014).

3.3 Choice of the free integrals

Free integrals are also expanded to the first order in the magnetic flux. The mass to magnetic flux ratio is similar to the one in the classical case, expanded as, $${\Psi }_A{}^2(\alpha )=\frac{4\pi {\rho }_0{}_{\star }{h}_{\star }^2}{{\xi }_{\star }}(1+\delta \alpha ).$$(56)

where δ is a free parameter describing the deviations from spherical symmetry of the ratio number density/enthalpy as in Meliani et al. (2006) and not of the density itself, conversely to Sauty & Tsinganos (1994).

The total angular momentum loss flux density is given by $$J=\gamma {\rho }_{0}Lh{V}_{\text{p}}=\frac{L{\Psi }_A}{4\pi }{B}_{\text{p}}.$$(57)

Thus it is natural to expand the quantity LΨ_A rather than L itself. LΨ_A is also the poloidal current density along the polar axis and writes as $$L{\Psi }_A=\lambda h_{\star }{B}_{\star }{r}_{\star }\alpha .$$(58)

The isorotation law can be expanded to first order as well as the total energy, $$\Omega ={\Omega }_{\star }(1+w_1\alpha ),$$(59)

and $$E=E_{\star }(1+e_1\alpha ),$$(60)

where we see from Eq. (24) that $E_{\star }=h_{\star }{\gamma }_{\star }{\xi }_{\star }{c}^2$.

Although we have some freedom with the choice of w₁ and e₁, we could choose e₁ = 0 and w₁ = −δ∕2 to restrict ourselves to the values of the previous models, in particular in Schwarzschild metric; see Meliani et al. (2006) and Meliani et al. (2010). In fact, the isorotation function Ω does not need to be expanded beyond the zeroth-order term because Ω always appears multiplied by another quantity as in (Ω − ω)ϖ or LΩ. Thus, the value of w₁ is free and does not affect the solution. Conversely, the value of e₁ affects the whole dynamics and we shall study the effects of its variation in a future publication. We already discussed the fact that taking a weak dependence of Ω on α has the advantage of minimizing the variation of the “light cylinder” near the base of the jet. Thus for the sake of simplicity, we shall study here the case where e₁ = 0 and w₁ = 0.

3.4 Constraints on the Alfvén Mach number, the isorotation law, and the angular momentum flux

The value of m₁ is, in fact, determined by the prescription to cross the Alfvén transition surface. In order for the denominator in Eqs. (25)–(27) to vanish at the Alfvénic transition, the two following relations given in Eqs. (29) and (30) must be fulfilled. They can be expanded to first order. For the first regularity condition we get, $${M_{\text{Alf}}|}_{\text{a}}=h_{\star }(1+m_1\alpha )\text{\hspace{0.22em}\hspace{0.22em}\hspace{0.22em}with\hspace{0.22em}\hspace{0.22em}\hspace{0.22em}}{m}_1=-\frac{\mu }{2}\left(\frac{{\lambda }^2}{{\nu }^2}+\frac{l^2}{{(1+l^2)}^3}\right).$$(61)

The first term in the right part of Eq. (61) is due to the “light cylinder”, and the second one to the nonsphericity of the gravitational field in Kerr metric. m₁ is negligible whenever the rotational speed λV_⋆ is sub-relativistic and either the μ parameter or the angular momentum of the black hole are negligible too. Since m₁ < 0, there is a limiting field line where we have M_Alf = 0 since the magnetic flux increases going out from the polar axis.

To apply the second regularity condition we use the numerator of Eq. (25) and we get to the first order in α: $${\Omega }_{\star }-{\omega }_{\star }=\frac{\lambda V_{\star }{h}_{\star }}{r_{\star }}.$$(62)

Thus we can write $$\begin{array}{lll}\varpi (\Omega -\omega )\hfill & =\hfill & G(R)\sqrt{\alpha }\lambda V_{\star }{h}_{\star }\Lambda (R),\hfill \end{array}$$(63)

where, $$\begin{array}{lll}\Lambda (R)\hfill & =\hfill & \left[1+\frac{\sqrt{\mu }\nu l}{\lambda h_{*}}\left(\frac{1}{{(1+l^2)}^2}-\frac{R}{{(R^2+l^2)}^2}\right)\right].\hfill \end{array}$$(64)

The regularity conditions on the Alfvén surface fixes the value of m₁. Thus the critical Alfvén surface is a sphere, as in previous meridional self-similar models. We warn, however, that the Alfvén transition surface is a generalized or modified Alfvén surface as it takes into account the modification by the “light cylinder”.

Simultaneously, surfaces of constant Poloidal Alfvén Mach Number, M_Alf = const. (see Eq. (55)) are not spherical surfaces, conversely to the one defined by Meliani et al. (2006). Two effects modify it; first the “light cylinder” effect, which was neglected in Meliani et al. (2006) and Globus et al. (2014), and second the frame-dragging effect (Lense-Thirring).

3.5 Expansion of the velocity and magnetic fields

The model is obtained using an expansion to the second order for sin θ in the Euler equation. Due to axisymmetry, first-order terms are zero along r and ϕ while the antisymmetry along θ gives the zeroth and second orders as null along the colatitude.

Then, for the poloidal velocity field, this gives $$\begin{array}{lll}{V}^{\widehat{r}}\hfill & =\hfill & \frac{V_{\star }{M}^2}{h_{\star }^2{G}^2}\{1+{\sin }^2\theta [\frac{1}{2}\left(\frac{l^2{h}_z^2}{R^2+l^2}-1\right)\hfill \\ \hfill & \hfill & +\frac{R^2+l^2}{G^2}\left(\frac{{\lambda }^2\mu }{{\nu }^2}\left(\frac{{\Lambda }^2{N}_B}{D}+\frac{{\overline{\omega }}_z}{\lambda }\right)-e_1-\frac{\delta }{2}+2m_1\right)]\}\hfill \\ V^{\widehat{\theta }}\hfill & =\hfill & -\frac{V_{\star }{h}_z{M}^2\sqrt{R^2+l^2}F}{2h_{\star }^{2}R{G}^2}\sin \theta .\hfill \end{array}$$(65)

And for the poloidal magnetic field, we get $$\begin{array}{lll}{B}^{\widehat{r}}\hfill & =\hfill & \frac{B_{\star }}{G^2}\left[1+\frac{1}{2}(\frac{l^2{h}_z^2}{R^2+l^2}-1){\sin }^2\theta \right],\hfill \\ B^{\widehat{\theta }}\hfill & =\hfill & -\frac{B_{\star }{h}_{z}F\sqrt{R^2+l^2}}{2G^{2}R}\sin \theta .\hfill \end{array}$$

Now from Eqs. (25) and (26), we can calculate to the first order in sin θ the toroidal components of fields $$\begin{array}{l}{V}^{\widehat{\varphi }}=-\frac{\lambda V_{\star }{h}_z\Lambda N_V}{h_{\star }{G}^{2}D}\sqrt{R^2+l^2}\sin \theta ,\\ B^{\widehat{\varphi }}=-\frac{\lambda B_{\star }{h}_{\star }\Lambda N_B\sqrt{R^2+l^2}}{h_{z}D{G}^2}\sin \theta ,\end{array}$$

where the functions N_V, N_B, and D have been generalized, $$\begin{array}{lll}\hfill & N_V\hfill & =\frac{M^2}{h_{\star }^2\Lambda }-G^2,\hfill \\ \hfill & N_B\hfill & =\frac{h_z^2}{h_{\star }^2\Lambda }-G^2,\hfill \\ \hfill & D\hfill & =\frac{h_z^2-M^2}{h_{\star }^2}.\hfill \end{array}$$

3.6 Expansion of the enthalpy, densities, and electric field

We used Eq. (27) to deduce the enthalpy, $$\gamma h\xi c^2={\gamma }_{\star }{h}_{\star }{\xi }_{\star }{c}^2\left[1+\alpha \left(e_1-\frac{{\lambda }^2\mu }{{\nu }^2}\left(\frac{{\Lambda }^2{N}_B}{D}+\frac{{\overline{\omega }}_z}{\lambda }\right)\right)\right],$$(73)

and the mass density is given by $$\begin{array}{ll}{\gamma }^2{\rho }_0\xi \hfill & ={\gamma }_{\star }^2{\rho }_0{}_{\star }{\xi }_{\star }\frac{h_{\star }^4}{h_z^2{M}^2}[1+\sin {\theta }^2\{\frac{\mu l^{2}R}{{(R^2+l^2)}^2}\hfill \\ \hfill & +\frac{R^2+l^2}{G^2}\left(2e_1-2m_1+\delta -\frac{2{\lambda }^2\mu }{{\nu }^2}\left(\frac{{\Lambda }^2{N}_B}{D}+\frac{{\overline{\omega }}_z}{\lambda }\right)\right)\}].\hfill \end{array}$$(74)

In GRMHD, we also need the expressions of the electric field and the charge density. The electric field is a second-order term for the radial component and a first-order term for the θ-component, $$\begin{array}{lll}\hfill & E^{\widehat{r}}\hfill & =-\frac{\lambda V_{\star }{h}_{\star }{B}_{\star }}{2c}\frac{(R^2+l^2)F\Lambda }{R{G}^2}{\sin }^2\theta ,\hfill \\ \hfill & E^{\widehat{\theta }}\hfill & =-\frac{\lambda V_{\star }{h}_{\star }{B}_{\star }}{c}\frac{\Lambda \sqrt{R^2+l^2}}{h_z{G}^2}\sin \theta .\hfill \end{array}$$

Using Maxwell–Gauss Eq. (7), we calculate the charge density from the divergence of the above electric field, to zeroth orderonly, $${\rho }_{\text{e}}=-\frac{\lambda V_{\star }{B}_{\star }{h}_{\star }}{2\pi r_{\star }c}\frac{\Lambda }{h_z{G}^2}.$$(77)

With all these quantities we are able to expand the Euler equation. The radial component is expanded to the second order and the colatitude component to the first order. From the expansion of poloidal components in the Euler equation (Eq. (16)) and using Eq. (51), we can reverse the system to get the equations of the model (see Appendix C for details).

3.7 “Light cylinder”

The rescaling value ch∕(Ω − ω) for the cylindrical radius used in Eqs. (32)–(34) has been defined by Meliani et al. (2006) as the “light cylinder”. It is a surface of revolution Σ_LC where, $$x^2=\frac{{\varpi }^2{(\Omega -\omega )}^2}{h^2{c}^2}=1.$$(78)

On the “lightcylinder”, the electric field |E| is equal to the poloidal component of the magnetic field |B_p|. In the present publication, Σ_LC designed the external “light cylinder”, that is, x = +1, though, strictly speaking, this is not a cylinder as explained earlier, but a surface of revolution. This external “light cylinder” is outside the Alfvén surface since the denominator of Eqs. (32)–(34) is equal to $M_{Alf}^2$ on the “light cylinder”, is negative before crossing the Alfvén transition surface, and is positive after crossing it. At large distance in the jet, the lapse function h goes to unity and Ω − ω tends to Ω, which is assumed constant in our model. Thus, Σ_LC is located on a constant cylindrical radius along the z axis, becoming a real cylinder.

From the iso-rotation law, we get, $$\frac{V^{\widehat{\varphi }}}{c}=x+\frac{{\Psi }_A{B}^{\widehat{\varphi }}}{4\pi {\rho }_0\gamma hc}=x+\frac{V_{\text{p}}}{c}\frac{B^{\widehat{\varphi }}}{B_{\text{p}}}.$$(79)

As in special relativity, the second term of Eq. (79) cannot be neglected in the vicinity of the “light cylinder”. The sign of $B^{\widehat{\varphi }}$ is such that $V^{\widehat{\varphi }}$ always remains less than the speed of light (Vlahakis 2015). Moreover after crossing the “light cylinder” one of the two following conditions must be fulfilled. Either, we have $\mid B^{\widehat{\varphi }}\mid$ » B_p or V_p » $V^{\widehat{\varphi }}$, or both.

The term x was neglected in the equation of the previous relativistic meridional-self-similar models, Meliani et al. (2006); Globus et al. (2014). Hence, these models could not produce jets crossing the “light cylinder”. Conversely, in this model this quantity is taken into account. We assume an expansion in sin(θ) of this quantity.

Contrary to the two previous models we can choose the dependence of the isorotation frequency with the magnetic flux (see discussion on Eq. (59)) and this choice will not affect the solution. Thus, if Ω does not depend strongly on the magnetic flux A, even at the base of the jet, the ratio h∕(Ω − ω) will be nearly constant. The reason for this is that ϖ is larger than the Alfvén radius, which is at least a few times the Schwarzschild radius. As a consequence, the departure of Σ_LC from a real cylinder is unnoticeable.

3.8 Domain of validity

The equations of the model are the result of an inversion of the expanded conservation equations. Therefore, it will be useful to quantify therelative error of the expansion we made in order to properly analyze our results and to obtain the domain of validity of these results. To have an idea of the domain of validity, we quantify the rest in the expansion of the momentum equation, for each force Fⁱ (R, sin θ), $$\begin{array}{lll}{F}^i(R,\sin \theta )\hfill & =\hfill & F_0^i(R)+F_1^i(R)\sin \theta +F_2^i(R){\sin }^2\theta \hfill \\ \hfill & \hfill & +R^i(R,\sin \theta ){\sin }^3\theta ,\hfill \end{array}$$(80)

where F is one of the following forces: gravitational, centrifugal, inertial, electric or magnetic pressure, and so on. We define a new function in order to map the relative error: $$R^i(R,\sin \theta )\underset{\theta \to 0}{~}{g}^i(R).$$(81)

For example, in the case of the electric force, in Schwarzschild metric, we get, assuming solid rotation (w₁ = 0), $$\begin{array}{lll}\mid R_{\text{El}}(R,\sin \theta )\mid \hfill & =\hfill & \frac{B_{\star }^2}{4\pi r_{\star }}\frac{{\lambda }^2{h}_{\star }^2\mu }{{\nu }^2}\frac{R}{h_z^2{G}^4}\left(\frac{F^2{h}_z^2}{2}+\frac{F{h}_z^2}{2}-3+\frac{\text{d}F}{\text{d}R}\right)\hfill \\ \hfill & \hfill & \times \sqrt{1+{\sin }^2\theta \left(\frac{F^2{h}_z^2}{4}-1\right)}.\hfill \end{array}$$(82)

The relative error on the electric force which tends to zero in the asymptotic regime of cylindrical jets is defined as $$\text{err}=\frac{\mid R^{\text{El}}(R,\sin \theta ){\sin }^3\theta \mid }{\mid F^{\text{El}}(R,\sin \theta )\mid }.$$(83)

Even at the base of the jet this error can be reduced, as can be seen in Fig. 1 for the solution in Kerr metric presented in Sect. 5.1 when the co-latitude is less than 30^°. To get an estimate of the error in the expanded forces, we should add all relative error terms or take the largest one. This gives an estimate of the domain of validity of the solutions for a given set of parameters. We postpone the full error analysis for a future paper.

Fig. 1 Relative error on the electric force for a recollimating oscillating solution in Kerr metric (K1, see Sect. 5.1). Color isocontours correspond to the relative error in the electric force. Field lines anchored into the black hole magnetosphere and in the accretion disk are plotted in black solid lines. The limiting field line between the inner jet coming from black hole corona and the outflow outgoing from the accretion disk is plotted in red. The “light cylinder” is indicated by a green solid line. The cylindrical radius and the distance above the equatorial plane are in units of Schwarzschild radius.

3.9 The magnetic collimation efficiency, ϵ

By writing the first law of thermodynamics in the frame of the fluid along streamlines of an axisymmetric flow, we can construct a constant of the motion, as in the classical case. The first law of thermodynamics reduces to the adiabatic law if the heating is included in some effective enthalpy (see Eq. (20)). Thus (ξc²) is an effective specific enthalpy, like for polytropic flows where the enthalpy also hides the heating (cf. Sauty & Tsinganos 1994) but generalized for relativistic outflows (see Eq. (18)). Using Eq. (28), we can rewrite the first law of thermodynamics in the following form, $$\xi {\Psi }_A^2{c}^2{\frac{\text{d}\xi }{\text{d}R}|}_{\alpha =\text{cst}}=4\pi M_{Alf}^2{\frac{\text{d}P}{\text{d}R}|}_{\alpha =\text{cst}}.$$(84)

As the magnetic-to-mass-flux ratio and the total energy flux are constant along each streamline, this is equivalent to $${\frac{\text{d}({\Psi }_A^2{\xi }^2{c}^2)}{\text{d}R}|}_{\alpha =\text{cst}}=8\pi M_{Alf}^2{\frac{\text{d}P}{\text{d}R}|}_{\alpha =\text{cst}}.$$(85)

We note that Ψ_Aξc² is proportional to the thermal energy. If we write $${\Psi }_A^2{\xi }^2{c}^2={[{\Psi }_A^2{\xi }^2{c}^2]}_0(R)+\alpha {[{\Psi }_A^2{\xi }^2{c}^2]}_1(R),$$(86)

and using the expressions of the pressure and the Mach number, we get an equation of the form $$\frac{\text{d}{[{\Psi }_A^2{\xi }^2{c}^2]}_0}{\text{d}R}+\alpha \frac{\text{d}{[{\Psi }_A^2{\xi }^2{c}^2]}_1}{\text{d}R}=B_{\star }^2{M}^2\frac{\text{d}\Pi }{\text{d}R}\left[1+(\kappa +2m_1)\alpha \right].$$(87)

We see, as in the classical case, that the second term of the pressure is proportional to the first one such that, $$\frac{\text{d}{[{\Psi }_A^2{\xi }^2{c}^2]}_1}{\text{d}R}-(\kappa +2m_1)\frac{\text{d}{[{\Psi }_A^2{\xi }^2{c}^2]}_0}{\text{d}R}=0.$$(88)

We deduce from the previous equation that the quantity ϵ, which is defined by $$ϵ(R)B_{\star }^2={[{\Psi }_A^2{\xi }^2{c}^2]}_1-(\kappa +2m_1){[{\Psi }_A^2{\xi }^2{c}^2]}_0=\text{cst}.,$$(89)

is a dimensionless constant for all the field lines. To give explicitly ${[{\Psi }_A^2{\xi }^2{c}^2]}_0$ and ${[{\Psi }_A^2{\xi }^2{c}^2]}_1$, it may be useful to write $${\Psi }_A^2{\xi }^2{c}^2={\Psi }_A^2\frac{{(h\gamma \xi c)}^2}{h^2}\left(1-\frac{{(V^{\widehat{\varphi }})}^2}{c^2}\right)-\frac{M_{Alf}^4{B}_{\text{p}}^2}{h^2}.$$

Finally the calculation leads to $$\begin{array}{l}ϵ=\frac{M^4}{h_z^2{h}_{\star }^4{G}^2(R^2+l^2)}[\frac{h_z^2{F}^2(R^2+l^2)}{4R^2}-\frac{R^2}{(R^2+l^2)}\\ -(\kappa -2m_1)\frac{(R^2+l^2)}{G^2}]-\frac{{\nu }^2(2e_1-2m_1+\delta -\kappa )R}{h_z^2(R^2+l^2)}\\ -\frac{{\nu }^2{l}^{2}R{G}^2}{h_z^2{(R^2+l^2)}^3}+\frac{2{\lambda }^2}{h_z^2}\left(\frac{{\Lambda }^2{N}_B}{D}+\frac{{\overline{\omega }}_z}{\lambda }\right)+{\lambda }^2{\left(\frac{\Lambda N_V}{h_{\star }GD}\right)}^2.\end{array}$$(90)

This equation is similar to Eq. (71) in Meliani et al. (2006) and can be interpreted the same way. The parameter ϵ measures the efficiency of the magnetic rotator to collimate the flow. At the outflow base, ϵ is the relative difference of the transverse variation of internal energy that is simply the exchange of work done by the macroscopic forces. As this is perpendicular to the flow axis, this means that ϵ really measures the transverse force which collimates the flow and mainly its magnetic component.

We note that the quantity − 2m₁ appears twice in Eq. (90). First, it is associated with κ, having a similar effect to the nonspherically symmetric pressure in the term which is given as a factor of M⁴. Second, it is associated with 2e₁ in the term corresponding to the excess or the deficit of the gravitational energy not compensated by the thermal driving at the base of the jet.

To conclude we can also derive the magnetic collimation efficiency in a different form; after some calculations, we can write $$\begin{array}{lll}ϵ\hfill & =\hfill & -\frac{{\nu }^2{h}_{\star }^4}{\mu {\gamma }_z^2{h}_z^2}\frac{\partial }{\partial \alpha }\ln {\left(\frac{P-P_0}{{\rho }_0\xi }\right)|}_{\alpha =0}\hfill \\ \hfill & =\hfill & h_{\star }^2\frac{{\nu }^2}{\mu }\left(1-\frac{\mu }{{\nu }^2}\right)\frac{{\xi }_z^2}{{\xi }_{\star }^2}\left({\frac{\partial }{\partial \alpha }\ln ({\rho }_0\xi )|}_{\alpha =0}-\kappa \right).\hfill \end{array}$$(91)

This new relation brings a link between the total enthalpy on the axis and its logarithmic variation with α. In particular, the sign of ϵ seems to connect the balance between logarithmic variation of total enthalpy per unit of volume and the meridional increase of the pressure. The factor indicates that |ϵ| probably tends to decrease for solutions which reach ultra-relativistic speed.

4 Methodology for obtaining solutions

4.1 Numeral integration

Appendix C gives the coupled ordinary differential Eqs. (C.1), (C.5), (C.6) for the four quantities of the Alfvén number, dimensionless radius, expansion factor, and pressure, (M², G², F,  Π). Details on the method and numerical techniques for the integration of these differential equations of the model system can be found in Sauty & Tsinganos (1994) and Meliani et al. (2006). Briefly, by using a Runge–Kutta scheme we start integrating from the Alfvén surface using the continuity relations. Integrating upwind by adjusting the value of the slope of the derivative of the Alfvén number at the Alfvén transition, the unique value of F_⋆ can be found that allows the crossing of the modified slow magnetosonic surface, for a given value of the pressure at the Alfvén transition Π_⋆. Then, integrating downwind, the value of Π_⋆ can be further adjusted. After several iterations of upstream and downstream integrations, the program automatically finds the solution that crosses all critical points, a proxy that the asymptotic pressure converges towards the demanded value. In particular, in the following Sections we have selected only solutions with the minimum possible value of Π_⋆, the so-called limiting solutions (see Sauty et al. 2004). We either have a collimated jet where, at infinity, Π_∞ is minimum, or a conical wind where, at infinity, Π_∞ is zero. If necessary, we can always add a constant value to the pressure P₀ to ensure that the pressure is positive everywhere in the flow.

In the following, we only outline briefly the points which differ from previous related studies for building our model in the framework of general relativistic MHD in a Kerr metric. Some interesting spine jet solutions, both for their properties close to the black hole and at large distance, are presented in Sect. 5. The presented solutions depend on a number of parameters and a systematic parametric study of the model is postponed for a following paper.

4.2 Alfvén regularity conditions

The regularity conditions at the Alfvén transition (R = 1) for the azimuthal components and enthalpy have already been discussed in Sect. 2.5. In order for the field lines to not have a kink at the Alfvén transition (i.e., F is a continuous function across R = 1), we impose an extra regularity condition on the transfield equation giving the four physical quantities (M², G², F, Π) at R = 1 (Tsinganos & Trussoni 1991, Sauty & Tsinganos 1994). In other words, similarly to the classical model Sauty & Tsinganos (1994), we should take appropriately into account this regularity condition at R = 1, which gives the ratio τ which is involved in the magnetic toroidal component, $$\begin{array}{lll}\tau \hfill & =\hfill & \frac{N_B}{D}{|}_{\star }=\frac{\frac{\text{d}{N}_B}{\text{d}R}{|}_{\star }}{\frac{\text{d}D}{\text{d}R}{|}_{\star }}\hfill \\ \hfill & =\hfill & \frac{h_{\star }^2\left(\frac{2}{1+l^2}-F_{\star }\right)-\frac{\mu (1-l^2)}{{(1+l^2)}^2}-\frac{h_{\star }l\sqrt{\mu }\nu }{\lambda }\frac{l^2-3}{{(1+l^2)}^3}}{p-\frac{\mu (1-l^2)}{{(1+l^2)}^2}},\hfill \end{array}$$(92)

where p is the slope of the square of the Alfvén number at the Alfvén transition, and $$p\equiv \frac{\text{d}{M}^2}{\text{d}R}{|}_{\star }.$$(93)

More specifically, the expansion factor F is determined by Eqs. (C.1) and (C.3). Also in Eq. (C.3) are the ratios of N_V and N_B with D. Furthermore, D appears in the denominator with a power higher by one than the powers of N_V and N_B in the numerator. Hence, in order that dF∕dR does not diverge and the slope of F(R) is continuous across R = 1 we require that ${\mathcal{N}}_F\cdot D{\mid }_{\star }=0$ such that dF∕dR is finite at R = 1 wherein we have 0/0. Thus, near the Alfvén transition, we may define the function ${\mathcal{N}}_F\cdot D~\mathcal{P}(M^2,G^2,F,\Pi )$. Then, in order to avoid a singularity at the Alfvén transition, a necessary condition is to choose p, F_⋆, Π_⋆ such that this function is zero, that is, $\mathcal{P}(M^2=h_{\star }^2,G^2=1,F_{\star },{\Pi }_{\star })=0$, and thus ${\mathcal{N}}_F\cdot D{\mid }_{\star }=0$. After some algebra we finally get a second degree polynomial for F_⋆, namely, $$\mathcal{A}(p)F_{\star }^2+\mathcal{B}(p)F_{\star }+\mathcal{C}(p,{\Pi }_{\star })=0,$$(94)

with $$\begin{array}{l}\mathcal{A}(p)={\lambda }^2{h}_{\star }^4+\frac{h_{\star }^2}{4}{\left(p-\frac{\mu (1-l^2)}{{(1+l^2)}^2}\right)}^2\\ \mathcal{B}(p)=[\frac{1}{2}{\left(\frac{\mu (1-l^2)}{{(1+l^2)}^2}-p\right)}^3\\ -2{\lambda }^2{h}_{\star }^2(p+\frac{2}{1+l^2}-\frac{2\mu (1-l^2)}{{(1+l^2)}^2}\\ +\frac{l\sqrt{\mu }\nu h_{\star }}{\lambda }\frac{3-l^2}{{(1+l^2)}^3})]\\ \mathcal{C}(p,{\Pi }_{\star })={\lambda }^2(p+\frac{2h_{\star }^2}{1+l^2}-\frac{2\mu (1-l^2)}{{(1+l^2)}^2}\\ +{\frac{l\sqrt{\mu }\nu h_{\star }}{\lambda }\frac{3-l^2}{{(1+l^2)}^3})}^2\\ +(\kappa {\Pi }_{\star }-\frac{1}{{(1+l^2)}^2}-\frac{2{\lambda }^2\mu }{{\nu }^2}\\ -2{\lambda }^2-\frac{l^2(2\mu +{\nu }^2)}{{(1+l^2)}^3}){\left(p-\frac{\mu (1-l^2)}{{(1+l^2)}^2}\right)}^2.\end{array}$$

In this way, the regularity condition at the Alfvén transition is automatically satisfied and no more constraints are needed at R = 1.

4.3 Effect of a nonspherical Alfvén number in a Schwarzschild metric

To illustrate our model, we present the following two solutions, which are built in the framework of the Schwarzschild metric. The first one corresponds to a solution of a model presented in Meliani et al. (2006), in which m₁ = 0.

The chosen values of the other parameters are λ = 1.0, κ = 0.2, δ = 1.2, ν = 0.8, ℓ = 0, μ = 0.1, e₁ = 0. We compare this solution to a solution with the same parameters but by keeping the value of m₁ given by Eq. (61), m₁ = −0.078. In both solutions, the value of Π_⋆ is the minimum value of the limiting solution. Such solutions have the minimum amplitude of oscillations in the jet.

In Fig. 2 the radial velocity on the polar axis is compared for the two solutions, while field lines in the poloidal plane are plotted in Fig. 3. We note that in Eqs. (C.1) and (C.4) and Eqs. (C.1) and (C.3), giving the plasma acceleration and the variation of the expansion factor with the radius R, respectively, there are several terms proportionalto the factor (κ − 2m₁). As m₁ is always negative, it is evident that − m₁ effectively increases the transverse pressure gradient, which is proportional to κ. In other words, the effect of − m₁ is similar to the effect of κ which enforces collimation for κ > 0. Hence, taking into account a nonsphericalAlfvén number (m₁≠0) introduces an extra collimation force which explains why the second solution with a nonzero m₁ is more collimated. Indeed, the width of the jet at infinity over its value at the base, G_∞ ∕G₀, decreases from 13.74 to 9.42, which is similar to an increase of κ and this fact can be checked directly by looking at the poloidal field line’s shape shown in Fig. 3.

Since more tightly collimated solutions have a smaller super-Alfvénic acceleration, the second solution reaches a lower Lorentz factor asymptotically. Thus, similarly to κ, the introduction of a negative m₁ ≠ 0 leads to a decrease of the velocity because of the tighter collimation. The higher collimation reduces the pressure gradient along the axis, which in turn decreases the acceleration due to pressure driving on large distances (see Sauty et al. 2004).

Additionally, m₁ appears within the term (κ + 2m₁ − δ − 2e₁) which appears in the plasma acceleration function ${\mathcal{N}}_{M^2}$ in Eqs. (C.1) and (C.4) and the function ${\mathcal{N}}_F$ determining the expansion factor F in Eqs. (C.1) and (C.3). The first three terms (κ + 2m₁ − δ) arise from the variation across the field lines of the heat content P∕ρ with α, that is, $\partial /\partial \alpha \{[\Pi (R)M^2(R)](1+\kappa \alpha ){(1+m_1\alpha )}^2/(1+\delta \alpha ]\}\simeq [\Pi (R)M^2(R)](\kappa +2m_1-\delta )$, while the fourth term e₁ is proportional to the variation of the total energy $\mathcal{E}$ with α. The bigger this term is, the larger the initial acceleration (see Sauty & Tsinganos 1994), because it is linked with the distribution of the heating which opposes gravity to accelerate the outflow. As the weight of the plasma decreases with the latitude, then the pressure gradient increases along the axis resulting in a larger acceleration close to the base, as explained in Tsinganos & Sauty (1992). This term decreases rapidly as the Alfvén surface is reached. Thus, it is responsible only for the initial acceleration. With the parameter m₁ being negative, the second solution is more accelerated between the base and the Alfvén surface. The velocity of the second solution reaches the velocity of the first one at the Alfvén surface. This effect disappears far from the source.

Fig. 2 Evolution of the radial velocity along the polar axis for solutions in a Schwarzschild metric, with m1 = 0 (blue) and m1 = −0.078 (red). The second case has a smaller terminal velocity.

Fig. 3 Field lines for a solution in a Schwarzschild metric, with parameters λ = 1.0, κ = 0.2, δ = 1.2, ν = 0.8, ℓ = 0, μ = 0.1, e1 = 0 and m1= 0 (left) and m1= −0.078 (right). We note that the case m1 = −0.078 corresponds to a more tightly collimated jet. Lengths are in units of the Schwarzschild radius. The red lines are connected to the magnetosphere of the central object while the green lines are connected to the disk. The separating line is in blue and the light cylinder in black.

5 Solutions in a Kerr metric

In the following, we discuss four different solutions in a Kerr metric to illustrate the present model. A more detailed parametric study is postponed to a following paper. For the purposes of the present paper, we show three cylindrically collimated solutions with high asymptotic Lorentz factor, typical of AGNs and GRBs. Those solutions cross the “light cylinder” and are sorted with increasing magnetic collimation efficiency parameter ϵ. We also exhibit a conical solution crossing the “light cylinder” with high Lorentz factor and strongly negative ϵ, something that was not possible with the previous relativistic meridionally self-similar solutions.

In order to get a Lorentz factor as high as possible in the asymptotic part of the collimated part of the jet, we know from the study of the classical solutions that among all cylindrical solutions, the limiting solutions with the lowest value of Π_⋆ reach the highest terminal velocity. These solutions are the so called limiting solutions in Sauty et al. (2002). As Π_∞ is negative for the limiting solutions, we have to add a positive P₀ value to thepressure. Of course, it is always possible for those cylindrical solutions to have a higher pressure P₀, but by doing so it also increases the effective temperature, in particular in the asymptotic part. For the same set of parameters, it is also possible to get cylindrical solutions by increasing Π_⋆. However, such solutions usually have a strong initial decollimation associated with a peak in both the Lorentz factor and temperature, while the asymptotic jet is decelerated to lower Lorentz factors and smaller radii, a result that we used to interpret the FRI/FRII dichotomy [cf. Meliani et al. (2010)].

Solutions K1 and K2 have been obtained for maximally rotating black holes, that is, a_H close to 1 (a ≃ r_s∕2). In solution K4 the value of a_H has been fixed to 0.73 (a = 0.73r_s∕2). We do not expect all black holes to be maximally rotating. For example, in M 87, the dimensionless spin should be above 0.65 (i.e., a > 0.65r_s∕2), (Li et al. 2009) but not too close to one. Other examples can be found and for K3 we will use the value a_H = 0.6 (a = 0.6r_s∕2) adopted in Mertens et al. (2016) for M 87.

Fig. 4 Three-dimensional (3D) representation of the field lines and streamlines for the thermally collimated solution K1 at the base of the jet and for two flux tubes. The blue lines correspond to streamlines, the red lines to magnetic field lines. The length is in units of the Alfvén radius, that is, ten times the Schwarzschild radius.

5.1 A mildly relativistic collimated solution with oscillations (K1)

The collimated solution K1 corresponds to an over-pressured outflow (κ ≥ 0) in a Kerr metric. As ϵ ≤ 0, the collimation of the jet is not fully magnetic but it has a significant contribution by the gas pressure, at least during the phase of strong acceleration up to ≃30r_s. In this solution, field lines are strongly oscillating compared to the two previous solutions in the Schwarzschild metric. The outflow undergoes a series of strong oscillations connected to the balance between the toroidal magnetic tension and the decollimation forces (centrifugal and electric forces and transverse pressure gradient) of the plasma.

The significant contribution of the transverse pressure gradient also explains these strong oscillations in the flow, as is also shown in the classical solutions. The parameters of this solution are displayed in the first row K1 of Table 1and the output values of m₁ and ϵ in Table 2. We note that m₁ = −0.062 is a relatively small value, which clearly indicates that the Alfvén surface is almost spherically symmetric in this case. As a consequence, the light cylinder is relatively far from the jet axis. Most of the central field lines (see the inner 5 to 7 central red lines in Fig. 5) remain within the “light cylinder”, which means that despite the important role of the magnetic field in the collimation, the jet is pressure or enthalpy driven in the relativistic case. However, unlike in the classical solutions, the electric field is the dominant decollimating force for the lines that cross the “light cylinder”. This decollimation and expansion after the Alfvén surface is associated with a strong pressure gradient yielding a strong acceleration of the jet in the super-Alfvénic regime. More details on this will be given in the following solution. The pressure gradient is the gas pressure gradient close to the axis but assisted by the toroidal magnetic pressure outside the “light cylinder”. This is similar to superfast flows in radially self-similar models for disk winds (see Vlahakis & Königl 2003a and Vlahakis & Königl 2003b). Moreover, in the relativistic regime the inertia increases faster when the flow is accelerated such that the collimation from the magnetic field is delayed to larger distances.

In Fig. 5 we see that the expansion after the Alfvén surface is strong and leads to a late acceleration of the flow. After the large expansion, the jet recollimates smoothly and consequently decelerates slightly because of the compression.

In Fig. 4, we clearly see that there is a strong azimuthal magnetic field, although the scale of the figure tends to exaggerate this phenomenon.

The Lorentz factor γ of this solution reaches a relatively small value around 3.7, typical of less powerful AGN jets like some of the FRI radio-galaxies.

Fig. 5 Poloidal field lines and “light cylinder” for the thermally collimated solution K1, for λ = 1.0, κ = 0.2, δ = 2.3, ν = 0.9, μ = 0.1, ℓ = 0.05, e1 = 0. The length unit is the Schwarzschild radius.

5.2 A highly relativistic collimated solution with oscillations (K2)

The solution K2 is collimated and has an extremely high Lorentz factor, which may be typical of GRBs. This K2 model corresponds to the values of the parameters given in the second line of Table 1, that is, λ = 1.0, κ = 0.2, δ = 1.35, ν = 0.46223, ℓ = 0.5, μ = 0.1, e₁ = 0 and the second line of Table 2 for the output parameters, m₁ = −0.234 and ϵ = −0.04.

For this model, the outflow starts very close to the black hole horizon at r₀ = 1.57 r_s, (see Fig. 7), and thus at the base of the jet the effects of general relativity play an important role. The final velocity is highly relativistic, as is shown in Fig. 6.

The parameter ν of the K2 solution is accurately adjusted (to the fifth digit), so as to obtain a rather high Lorentz factor (larger than 100). This proves the versatility of the model which handles any magnitude of Lorentz factors. In order to obtain such high Lorentz factors, we must carefully tune the parameter directly linked to gravitation, ν, as mentioned above. The same parameter is also responsible for the thermal acceleration in the classical model (Sauty & Tsinganos 1994).

In Figs. 8 and 9 we plot the forces along and perpendicular to a field line defined by α =0.01α_lim where α_lim is the dimensionless magnetic flux between the inner jet and an external accretion disk wind.

The strong decollimation associated with the slow acceleration enhances the electric force as in the previous solution K1. This can be seen in Fig. 9. However due to the higher rotation here, more field lines cross the light cylinder, which is very close to the axis, such that the decollimation from the electric field is much stronger in this solution.

The feedback of this strong electric field is to further increase the decollimation beyond the Alfvén surface at large distances. Again, the large expansion increases the pressure and enthalpy gradient as seen in Fig. 8. Thus, the pressure force increases, resulting in a very long acceleration phase up to 10⁶ r_s. Thus, the thermal acceleration becomes very efficient and the plasma asymptotically reaches an extremely high Lorentz factor. Indeed, the jet radius increase of this solution is also very high with an expansion factor G_∞ ∕G₀ ≃ 1300, as can be seen in Fig. 7.

The value of ϵ is still negative but very close to zero (ϵ = −0.04). As in the classical case, this means the magnetic efficiency to collimate the flow is higher in this model at larges distances. However the decollimating force that ensures the equilibrium is no longer the centrifugal force or the pressure gradient but the electric force on the lines that cross the “light cylinder”. This is a specific feature of relativistic jets.

To further analyze the jet acceleration we may calculate the contribution of the different components of the total energy and the conversion of the magnitude of each component to another form along the streamlines. First, we want to obtain a physically acceptable heating term which goes to zero at infinity along the axis of the flow. By defining the enthalpy in Eqs. (18) and (19) and its analytical expression in our model in Eq. (34), we may fix a streamline limiting the zone where the pressure is positive. Once this field line has been chosen, for example, at α =0.4α_lim, the dimensionless pressure P₀ ∕ρ_⋆c² can be calculated. In the particular case presented here, we have taken P₀∕ρ_⋆c² = 4.3 × 10⁻⁷. The pressure P₀ in Eq. (54) is chosen such that the gas pressure is equal to zero, when Π (R) reaches its minimum value. For some value of α, the term $P_0=P_0(\alpha )=-[B_{\star }^2/8\pi ]{\Pi }_{\text{min}}(1+\kappa \alpha )$ will be too large to ensure that Q goes to zero at infinity along the axis. Indeed, α can take any value below a maximum α_max. We can consider that the solution can be valid only in the region wherein the pressure is positive. When α is fixed and P₀ ∕ρ_⋆c² is deduced, we are able to calculate ξ_⋆. In this particular case, a value ξ_⋆ ≃ 78 is chosen.

Figure 10 shows the normalized total energy on the axis $\mathcal{E}/c^2=h\gamma \xi$, the kinetic component hγξ_K and the external heating hγQ∕c². We find a decrease of the external heating and a related increase of the kinetic part. Thus, the kinetic enthalpy represents the major component up to r = 10²r_s.

Figure 11 shows the same energetic distribution, but on a streamline with α = 0.05α_lim. Out of the polar axis, there are extra energetic components, such as the frame-dragging and the Poynting fluxes. Both these energetic contributions are very small on this field line, as compared to the total energy. Hence, the jet is enthalpy-driven from the axis right up to the limiting line. We also note that at the base of the jet, the frame-dragging energy is of the same order as the Poynting flux. Contrary to the energetic distribution along the axis, the external heating constitutes the larger part of energy at infinity. While along the axis a high value of γ_∞ ≃ 100 is obtained, the Lorentz factor at infinity on this particular line is γ_∞≃ 3.6. Hence, since the acceleration of the plasma and the resulting final flow speed at infinity on this particular field line are small, the external heating is not consumed when accelerating the flow and therefore it is left unused at infinity, contrary to what happens along the axis; Fig. 10.

Fig. 6 Lorentz factor γ for the K1 (blue line) and K2 (red line) solutions. Distances are given in Schwarzschild radius units.

Fig. 7 Poloidal field lines and “light cylinder” for the K2 solution, that is, for λ = 1.0, κ = 0.2, δ = 1.35, ν = 0.46223, μ = 0.1, ℓ = 0.05, e1 = 0. Distances are given in Schwarzschild radius units.

Fig. 8 Longitudinal forces, that is, along the field line, for the K2 solution, along the line α = 0.01αlim. Distances are given in Schwarzschild radius units.

Fig. 9 Transverse forces, that is, perpendicular to the field line, for the K2 solution, along the line α = 0.01αlim. We see that the Lorentz force is collimating and is balanced by the electric force that decollimates. Distances are given in Schwarzschild radius units.

Fig. 10 Relative normalized contribution to the total energy of the kinetic enthalpy hγξk and the external heating distribution hγQ∕c2 along the axis.

Fig. 11 Relative normalized contribution to the total energy of the kinetic enthalpy hγξk and the external heating distribution hγQ∕c2 along the streamline α = 0.05αlim

5.3 A mildly relativistic collimated solution without oscillations (K3)

From VLBI imaging, Mertens et al. (2016) recovered a detailed two-dimensional velocity field in the jet of M 87 at sub-parsec scales. They confirmed the stratification of the flow from the very beginning of the jet and identified a relativistic sheath, that is, an accelerating layer, which is launched from the inner part of the accretion disk at a cylindrical distance of around 5r_s. Mertens et al. (2016) interpret this outer sheath layer as the internal part of an external disk wind. They also interpret the inner spine jet as a component coming from the internal accretion disk. However, it is not clear that the inner spine necessarily has to be one of the disk wind components. Instead we propose that the inner spine jet originates from the black hole corona and has a higher Lorentz factor. The authors note that this fast inner spine jet cannot be detected in their data because of its lower emissivity compared to that of the sheath layer or because its speed is too high.

We propose here, as an alternative scenario, that the spine beam may originate from the magnetosphere of the black hole, either connected to the black hole itself, as in the Blandford–Znajek–Penrose mechanism, or connected to the inner part of the accretion disk. In the first case, the spine jet would be a leptonic plasma, and in the second it would be a hadronic one. In such a case, we can model the spine jet with our meridionally self-similar solutions. We note that unlike Poynting-flux-dominated models, our model is valid on the jet axis.

The Lorentz factor profile inferred by Mertens et al. (2016, see their Fig. 19) supposes that the velocity structure observed in the 43 GHz VLBA maps at a deprojected distance of z ≃ 10 mas ≃ 1300r_s is due to the sheath layer. The Lorentz factor of the sheath has a value γ ≃ 2.4 at R ≃ 1000r_s.

From the curve deduced for the spine jet with an approximated MHD disk wind solution, those authors, following Anderson et al. (2003), model the acceleration and collimation of the flow with a Lorentz factor γ ≈ 7. We propose that it is possible to construct a MHD solution along the lines of the model analysed in this paper, by choosing a positive value of ϵ ≈ 0.5 and a similar velocity profile. The solution K3 presented here has a Lorentz factor γ ≃ 5 at R ≃ 1000r_s. By appropriately tuning the parameters, we could obtain even higher Lorentz factors from 7 to 10, as is observed at the distance of the knot HST-1, where this velocity is observed in the optical band. However, the spine jet is deboosted relative to the sheath and a precise measurement of its Lorentz factor is difficult.

Solution K3 with a maximum Lorentz factor on the axis γ ≃ 5.5 has an asymptotic spine jet radius G_∞≃ 20r_s. At this distance from the axis, the Lorentz factor has dropped to a value γ ≃ 2.4 consistent with the radius and the Lorentz factor at the inner observed distance of the outer sheath jet of Mertens et al. (2016). Thus, our solution may model the initial spine jet inside the sheath layer. However, to confirm that, we need to use this initial solution in simulations similar to those in Hervet et al. (2017). The spine jet/sheath jet interaction will probably produce shocks and rarefaction waves that may further accelerate the jet.

Therefore, although we can obtain such a type of collimated solution for different sets of the parameters, we focus here on the specific solution K3, where λ = 1.2, κ = 0.005, δ = 2.3, ν = 0.409, ℓ = 0.024 and μ = 0.08. Compared to the other Kerr solutions studied in this paper, λ is higher and κ is very small, leading to a solution with a positive value of the magnetic collimation efficiency, ϵ = 0.55. Another advantage of this solution is the fact that the pressure depends only very weakly on the magnetic flux function, that is, on a particular field line.

Interestingly, in Mertens et al. (2016), the radius of the base of this spine jet is equal to r₀ = 2.4r_s, which is consistent with their model of a disk wind solution, by fixing the jet shape and solving the Bernoulli equation. In our self-similar solution, we have a similar radius for the magnetospheric polar cup, where our jet solution starts. Clearly, this is an alternative scenario.

Moreover, the angular velocity of the field lines anchored in this polar cup above the black hole can be calculated from our parameters, $$\begin{array}{lll}{\Omega }_{\star }\hfill & =\hfill & \frac{c\mu }{r_s}\left[\frac{\sqrt{\mu }\lambda }{\nu }\sqrt{1-\frac{\mu }{(1+l^2)}}+\frac{l\mu }{{(1+l^2)}^2}\right]\hfill \\ \hfill & \simeq \hfill & 6.2\times {10}^{-2}\frac{c}{r_s}.\hfill \end{array}$$(98)

By taking the value of the M 87 distance and the black hole mass, as in Mertens et al. (2016), we can calculate Ω_⋆ in the context of this K3 solution. This value may directly be compared to the values they deduced in two jet regions from the conservationof total energy and angular momentum fluxes in the approximation of special relativity.

Hence, we find Ω_⋆ ≃ 1.03 × 10⁻⁶s⁻¹, a value which is for the spine jet almost the same value as the isorotation frequency of a Keplerian speed at the launching location of the sheath layer. It also corresponds to the initial toroidal velocity of the Blandford & Payne (1982) mechanism.

Our model corresponds to an alternative configuration, because the spine jet may either originate from the Keplerian disk, in which case it would be hadronic, or form via the generalized Penrose–Blandford–Znajek mechanism. In this second alternative, the jet would be a leptonic beam with an angular frequency proportional to the spin of the central black hole. The Blandford–Znajek mechanism allows us to extract energy from the black hole when 0 ≤ Ω ≤ ω_BH, with a maximum value for 0.5ω_BH. We note that ω_BH is, by definition, the angular velocity of ZAMO at the location of the outer event horizon and is given by, $$\begin{array}{lll}{\omega }_{\text{BH}}\hfill & =\hfill & \frac{a_{H}c}{r_s\left(1+\sqrt{1-a_H^2}\right)},\hfill \end{array}$$(99)

where a_H is the dimensionless spin of the black hole in units of the gravitational radius r_s ∕2. Indeed simulations of such Poynting-dominated and force-free jets (Tchekhovskoy 2015) have shown that the angular speed of a field line anchored in the magnetosphere is about half the black hole angular speed ω_BH. Our value of Ω_⋆ is one third of 0.5ω_BH ((Nathanail & Contopoulos 2014)). In order to determine if the spine jet of our solutions originates from a Keplerian disk or froma black hole via a generalized Penrose–Blandford–Znajek mechanism, we need to solve the MHD equations up stream up to the black hole horizon. Moreover, in order to model the full jet of M 87, a complete MHD simulation including a disk wind and a spine jet must be developed. Something that should be done in the future.

We notethat as is already known, at the interface of the spine jet and the sheath layer, a re-collimation shock may occur, producing compression and rarefaction waves which may accelerate the flow (Hervet et al. 2017). For those reasons we have chosen here to adjust the value of ν to ν = 0.42, in order to get the solution K3, with a lowerLorentz factor but suppressing completely the oscillations of the field lines. The value calculated above for Ω_⋆ is not changed, but radio emission maps, which are obtained for the M 87 jet, will be produced by re-collimation shocks due to the interaction between the fast spine jet and the sheath layer.

As can be seen in Fig. 12 for the K3 model, the Lorentz factor reaches a nearly constant value at the distance of the B structure observed in the 43 GHz VLBA maps and γ is larger than that deduced for the sheath layer by Mertens et al. (2016). As explained above, the interaction between the sheath layerand the spine jet can induce a bulk flow acceleration. However, as we mention, the Lorentz factor γ is maximum along the axis and decreases with latitude such that at its outer boundary it matches the sheath layer value.

The values of the parameters of the K3 model are given in the third line of Tables 1 and 2. The radius at the base of this spine jet is equal to r₀ = 2.55r_s and the magnetic efficiency to collimate the flow is larger compared to the other two Kerr solutions discussed in this paper, as ϵ is now positive, which means that the jet is fully magnetically collimated. The field lines displayed for K3 in Fig. 13 are nearly cylindrical above the equatorial plane, at distances 10⁴ r_s, with a smooth flaring occurring after the Alfvén surface.

Fig. 12 Lorentz factor for the nonoscillating collimated solution K3.

Fig. 13 Poloidal field lines and “light cylinder” for the nonoscillating collimated solution K3, that is, for λ = 1.2, κ = 0.005, δ = 2.3, ν = 0.42, μ = 0.08, ℓ = 0.024, e1 = 0. Distances are given in Schwarzschild radius units.

Fig. 14 Poloidal field lines and “light cylinder” for the conical solution K4, that is, for λ = 0.0143, κ = 1.451, δ = 3.14, ν = 0.8, μ = 0.41, ℓ = 0.15, e1 = 0. Distances are given in Schwarzschild radius units.

5.4 Conical solution

The parameters of the conical solution K4 are given in Tables 1 and 2. Such radial solutions could be useful to describe relativistic noncollimated outflows, such as those seen in association with radio-quiet galaxies, such as Seyfert galaxies. Conical solutions could also be useful to model GRBs, wherein an unstable noncollimated relativistic wind may fragment into small pieces under some instabilities. In such cases, the apparent collimation of the GRB would be due to fragmentation – see Meliani & Keppens (2010); van Eerten et al. (2011).

A conical solution is a solution in which the limiting value of F at infinity is 0. In this solution, the spherical part of the Mach number diverges, M →∞. The same effect occurs for the cylindrical radius of the flow, G. To obtain this conical solution, we started with the parameters used for modeling the solar wind, as in Sauty et al. (2005). The magnetic collimation parameter must be strongly negative, ϵ∕(2λ²) ≤ 0. We adjust the solution to the relativistic case and increase the velocity by taking a larger value for δ. For the radial solution K4, the parameters λ and κ are adjustedas well in order to obtain a terminal Lorentz factor larger than 8. Thus, we find a conical solution K4 for the following set of parameters: λ = 0.0143, κ = 1.451, δ = 3.14, ν = 0.8, ℓ = 0.15, μ = 0.41, m₁ = −0.004 and e₁ = 0.

For thisset of parameters we get a much more negative value for the magnetic collimation efficiency parameter, ϵ = −5.68. The solution quickly reaches the conical regime (see Fig. 14), and most of the field lines cross the “light cylinder” plotted as a black solid line. The axial radial velocity profile is plotted in Fig. 15. As can be seen, high Lorentz factors are obtained.

5.5 Magnetic collimation efficiency versus black hole spin

As already mentioned, the constant ϵ is a measure of the efficiency of the magnetorotational forces to collimate the flow. We have studied how ϵ relates to the black hole spin a_H = 2l∕μ, keeping unchanged all the other parameters: λ, κ, δ, ν, and μ.

In Fig. 16, we plot ϵ versus the black hole spin a_H, for several cases studied in the context of our model. First we note that for solutions with parameters similar to the K3 solution, referred to here as K3-type solutions, the value of ϵ is positive, it increases with the black hole spin a_H and shows the largest variation in relative magnitude with an absolute total variation equal to 0.05. On the other hand, for K1-type and K4-type solutions, the values of epsilon are negative. The total variation for K1-type solutions is equal to 0.10 and for the K4-type solutions is equal to 0.25. We see that ϵ is increasing with the black hole spin a_H for K1-type and K3-type solutions. It also presents a minimum value for K4 at a black hole spin slightly smaller than a_H = 0. This means that the magnetic collimation efficiency is lower for counter-rotating black holes in relation to their accretion disk for K1-type and K3-type solutions and depends weakly on the spin direction for K4-type ones. However, K4-type solutions are conical, contrary to all other solutions, which are cylindrically collimated. In all cases ϵ does not vary linearly with the black hole spin a_H.

This nonlinear variation of epsilon with the black hole spin a_H is explainedif we try to derive it at the base of the jet where the Alfvénic number is equal to 0, since some terms of the second order in l cannot be neglected in the following equation. $$\begin{array}{ll}ϵ\hfill & =\frac{2{\lambda }^2}{h_z^2}\left(\frac{{\Lambda }^2{N}_B}{D}+\frac{{\overline{\omega }}_z}{\lambda }\right)+{\lambda }^2{\left(\frac{\Lambda N_V}{h_{*}{G}_{0}D}\right)}^2\hfill \\ \hfill & -\frac{{\nu }^2(2e_1-2m_1+\delta -\kappa )R_0}{h_z^2(R_0^2+l^2)}-\frac{{\nu }^2{l}^2{R}_0{G}_0{}^2}{h_z^2{(R_0^2+l^2)}^3}.\hfill \end{array}$$(100)

An increase in ϵ coincides with a decrease in the maximal Lorentz factor for the collimating Kerr solutions, as we will see for K1-type and K3-type solutions. This explains why it was not possible to obtain physical solutions by decreasing l for K2-type solutions. For this second Kerr solution, we performed a fine tuning of the parameters to obtain the largest possible Lorentz factor at large distances from the maximally rotating black hole. Then, decreasing l leads automatically to exceed the value of the speed of light, c, for the polar velocity at some distance in the jet. The acceleration phase does not vary for K1-type and K3-type solutions with l except immediately before R = 500 where the Lorentz factor reaches a plateau. The value of the maximum γ increases when l decreases.

For the two collimated solutions K1 and K3, there is a clear effect of collimation induced by the rotation of the black hole. In Fig. 17 we plot for five values of the black hole spin the evolution of the factor G along the jet, that is„ the ratio of the jet cylindrical radius divided by its value at the Alfvén surface. From this plot it can be seen that the maximum of the radius of the jet is reached at different distances, as the spin varies, and decreases when a_H increases. The same trend is observed for the terminal jet radius but the ratio of G_∞ /G₀ decreases only by a factor of 0.96 between a nonrotating and a maximally rotating black hole (a_H = 0.99). Hence, the faster the black hole rotates, the smaller the maximal jet radius. This result is expected because for a fast black hole rotation, the magnetic collimation efficiency parameter (ϵ) is higher.

The case of K3-type solutions is simple, as the factor G increases with distance until it reaches a constant value. The ratio G_∞ /G₀ gives directly the expansion factor which is decreasing when a_H increases from − 0.99 up to 0.99.

As can be seen from Fig. 18, the Lorentz factor maximum follows the opposite trend to ϵ for the conical K4-type solutions: the minimum value of ϵ is obtained for a_H ≃ 0 and for a nonrotating black hole the Lorentz factor curve reaches a maximum before decreasing up to a plateau at large distances. This type of curve is no longer observed when the absolute value of a_H goes above some threshold. The increase of the plateau value for the Lorentz factor is much more pronounced for a_H > 0 but can be seen also for negative values of the spin. At a distance of r = 1000r_s γ increases from a value of the order of 6 for a_H = 0 up to 14 for a_H = 0.99.

For the conical solutions wherein the asymptotic geometry is fixed, their geometry is less affected. The effect of varying ϵ is different in this case, as compared to this effect for the cylindrically collimated solutions. To analyze the effect of collimation induced by the black hole spin, we plot the ratio of the cylindrical to spherical radius of the field lines for K4-type solutions and five different values of the black hole spin. We note that the ratio of the cylindrical to spherical radius of the field lines gives their opening angle with respect to the axis. The field lines start at the base of the jet with an opening angle which increases with the radial distance and expands away from the polar axis more rapidly for higher black hole rotation. The effect is more pronounced for negative spin parameters. At some distances above the equatorial plane (z ≃ 100r_s) the opening angle becomes constant as the field lines are finally becoming radial. It is found that for a given ratio of α∕α_lim, where α_lim is the last open field line, the asymptotic opening angle is constant for all K4-type solutions, regardless of the spin of the black hole.