© The Authors 2024

1. Introduction

Astrophysical jets are powerful and highly collimated outflows of plasma that emanate from various astrophysical objects, such as active galactic nuclei (AGN), quasars, X-ray binaries, young stellar objects (YSOs), and gamma-ray bursts (GRBs). The formation and collimation of jets has been dealt with analytically and most recently via numerical simulations in different regimes, from hydrodynamics to relativistic and general relativistic magnetohydrodynamics. The YSO jets are nonrelativistic, with typical velocities of ∼100 km s⁻¹ Woitas et al. (2002), whereas AGN and GRB jets are relativistic, with Lorentz factors of γ ∼ 10 and γ ∼ 100, respectively (see e.g., Jorstad et al. 2001; Vlahakis & Königl 2004).

The idea behind magnetohydrodynamic (MHD) models for disk jets or coronal jets is that the gravitational energy extracted from the central object is transferred to the accreting plasma, which produces the jet via a collimation and an acceleration mechanism. The energy released by the accretion increases with the mass of the central object.

With regard to the energy source of relativistic jets, it is believed that they can be powered by a black hole’s gravitational pull or by its rotation. The extraction of energy from a rotating black hole presumes the existence of an ergo-region just outside the black hole’s horizon. When a magnetic field is present, the extraction can occur via either the electromagnetic field or the internal energy of the plasma. In the case in which the electromagnetic flux is dominant, Blandford & Znajek (1977) proposed a model of electromagnetic extraction of a black hole’s rotation energy based on the analogy of the Faraday disk phenomenon, for strongly magnetized jets. This mechanism, which neglects the inertia of the plasma, suggests that the energy is extracted via the Poynting flux of the electromagnetic field threading the black hole. In the second case, in which the plasma inertia dominates, we have the Penrose process Penrose (1969) whereby a particle that falls inward decays into two in the ergo-region. One particle is being absorbed by the black hole and the other one reaches infinity, with energy larger than the initial energy of the infalling one. The energy gain in this (“mechanical”) Penrose process is explained by the negative energy (as seen from infinity) of the ergo-region trapped particle absorbed by the black hole.

In order to explain the acceleration of an astrophysical outflow, many mechanisms have been proposed. The thermally driven models derive from the “de Laval nozzle” analogy of the solar wind (Parker 1963), which requires the presence of hot corona around the central body. The magnetic pressure models that have been examined by Draine (1983) and Contopoulos (1995) assume a toroidal magnetic field, B_ϕ, which is created and amplified by a non-Keplerian rotating disk. Finally, in the magnetocentrifugal-driven models the plasma is flung out even to relativistic velocities from an accretion disk (Blandford & Payne 1982; Pelletier & Pudritz 1992; Cao 1997). In this case, the flow is forced to follow the field lines and rotate with them while being magnetocentrifugally accelerated, as long as the angle between the poloidal magnetic field and the disk is no more than 60°. In the case of relativistic jets, the acceleration may be efficient even at bigger angles, as was shown by Cao & Spruit (1994).

In regard to the collimation of the jet, this is mainly due to the toroidal magnetic field generated by the rotation of the source Bogovalov (1995). Magnetic self-collimation is efficient in a nonrelativistic context Livio (2002), Honda & Honda (2002), Tsinganos & Bogovalov (2002). Moving close to relativistic limit, the electric force as well as the higher inertia of the outflow makes collimation slower, but still remains a possible scenario Vlahakis & Königl (2003, 2004). In the case of relativistic jets, collimation may alternatively be due to the external pressure of a surrounding disk wind Bogovalov & Tsinganos (2005). Obtaining exact solutions with any MHD model, despite several simplifications such as axisymmetry, remains a hard task. Using the hypothesis of self-similarity is a way to overcome this difficulty. Self-similarity assumes a scaling law of one of the variables as functions of one of the coordinates. The choice of the scaling variable depends on the astrophysical problem. Radial self-similarity is usually used in disk wind models, but these solutions cannot describe the flow close the axis. Meridional self-similarity, on the other hand, provides a better way to study the jet close to the symmetry and rotation axis of the central corona.

Tsinganos & Trussoni (1991) developed self-similar models for MHD flows and then Sauty & Tsinganos (1994) proposed two types of solutions and a criterion that allows one to characterize the jets. Meliani et al. (2006) extended these Newtonian models to the case of relativistic jets emerging from a spherical corona surrounding the central part of a Schwarzchild black hole and its inner accretion disk. Then Meliani et al. (2010), using the previous relativistic solutions of Meliani et al. (2006) and the collimation criterion of that model, studied how the FRI, FRII dichotomy may influence the morphology of the inner spine jets. Globus et al. (2014) have extended the model around rotating black holes using the Kerr metric and Chantry et al. (2018) further generalized the model in Kerr metric for non-polytropic solutions and showed that the jet power is of the same order as that determined by numerical solutions (e.g., McKinney 2005; McKinney et al. 2012; Tchekhovskoy et al. 2011). Following a different path than that of the self-similar models, Hou et al. (2024) produce analytical solutions, assuming that the streamline maintains a constant value of polar angle, θ (conical solution). This assumption allows them to derive analytical expressions for the magnetic field, the accretion flow, and the funnel wall jet without the need to solve the transfield equation.

Numerical simulations, both in special relativistic MHD and general relativistic MHD (GRMHD) frames, are an alternative to analytical and semi-analytical models. They are used to model the jet’s formation and study the two-component structure of it. In special relativity, simulations of corona winds have been presented by Bogovalov & Tsinganos (1999) and Tsinganos & Bogovalov (2002), whereas Koide et al. (1999, 2001) worked in the general relativity frame. Bogovalov & Tsinganos (2005) proposed a two-component model in which a relativistic central wind is collimated by a surrounding nonrelativistic disk wind due to the fact that single-component models presented difficulties with regard to collimating the outflow. McKinney & Blandford (2009) studied the stability of jets from rotating, accreting black holes via fully three-dimensional (MHD) simulations. Millas et al. (2017), using the AMR-VAC simulation code, studied the role of toroidal velocity and magnetic fields in the stability of two-component jets. Simulations in the GRMHD frame show that magnetized accretion is able to accumulate poloidal magnetic fields on the black hole. These fields become important enough and the system enters a magnetically arrested state (MAD, Narayan et al. 2003; Igumenshchev 2008) that is able to launch jets whose power exceed the accretion power Tchekhovskoy et al. (2011). A scaling law between the accretion rate and the magnetic flux threading the black hole horizon has been found by McKinney et al. (2012) and Tchekhovskoy et al. (2011). Nakamura et al. (2018) use semianalytical solutions of the steady, axisymmetric force-free electrodynamic jet model (Narayan et al. 2007; Tchekhovskoy et al. 2008) together with the GRMHD simulation code HARM (Gammie et al. 2003) to study the formation of parabolic jets from the spinning black hole in M 87.

Ressler et al. (2021), with the use of 3D GRMHD simulations, study gas accretion onto a rapidly spinning black hole with a Bondi to an event horizon scale separation of 100. They predict that these low-angular-momentum accretion jets will eventually collapse due to kink instabilities. Lalakos et al. (2022) expand this to larger-scale separation. The launching of relativistic jets from Bondi-like accretion flows with either zero or low angular momentum is also studied by Kwan et al. (2023). This study suggests that above a specific threshold the flow can sustain an MAD and eventually launch a powerful jet. Lalakos et al. (2024) considers weakly magnetized zero-angular-momentum gas accretion to study the formation and propagation of jets and specifically the factors that contribute to their collapse.

In this paper, we provide an extension of the nonrelativistic meridional self-similar solutions of Sauty & Tsinganos (1994) in the Kerr metric. As in the case of Chantry et al. (2018), 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 the previous model, this new one assumes a polytropic law for the flow. In Sect. 2, we present the MHD equations in a general relativistic framework under the Kerr metric, in the 3+1 formalism. Section 3 deals with the construction of the model and its equations, and in Sect. 4 we present the solutions. Section 5 is dedicated to summarizing our results and conclusions.

2. Basic equations

2.1. Kerr metric

We begin by defining the metric. The Kerr metric in the Boyer–Lindquist coordinate system is

$$\begin{array}{cc}\hfill d{s}^2& =-(1-\frac{r_{\mathrm{s}}r}{{\rho }^2})c^2\mathrm{d}{t}^2-\frac{2r_{\mathrm{s}}rc\alpha }{{\rho }^2}{sin}^2\theta \mathrm{d}t\mathrm{d}\varphi \hfill \\ \hfill & +\frac{{\rho }^2}{\mathrm{\Delta }}\mathrm{d}{r}^2+{\rho }^2\mathrm{d}{\theta }^2+\frac{{\mathrm{\Sigma }}^2}{{\rho }^2}{sin}^2\theta \mathrm{d}{\varphi }^2,\hfill \end{array}$$(1)

using the following notations for the elements:

$$\begin{array}{cc}& \mathrm{\Delta }=r^2+{\alpha }^2-r_{\mathrm{s}}r,\hfill \end{array}$$(2)

$$\begin{array}{cc}& {\rho }^2=r^2+{\alpha }^2{cos}^2\theta ,\hfill \end{array}$$(3)

$$\begin{array}{cc}& {\mathrm{\Sigma }}^2={(r^2+{\alpha }^2)}^2-{\alpha }^2\mathrm{\Delta }{sin}^2\theta ,\hfill \end{array}$$(4)

$\alpha =\frac{J}{\mathit{Mc}}$, and $r_{\mathrm{s}}=\frac{2GM}{c^2}$. J is the angular momentum of the massive black hole and M is its mass. Below, we write the lapse function, h, the angular velocity, ω, of the zero-angular-momentum observers (ZAMOs), the shift vector’s coordinates, β_ϕ,  β^ϕ, and the line elements, h_r,  h_θ,  h_ϕ, respectively.

$$\begin{array}{cc}& h={(1-\frac{r_{\mathrm{s}}r}{{\rho }^2}+{\beta }^{\varphi }{\beta }_{\varphi })}^{1/2}=\frac{\rho }{\mathrm{\Sigma }}\sqrt{\mathrm{\Delta }},\hfill \end{array}$$(5)

$$\begin{array}{cc}& \omega =\frac{\alpha c{r}_{\mathrm{s}}r}{{\mathrm{\Sigma }}^2},{\beta }^{\varphi }=-\frac{\omega }{c},{\beta }_{\varphi }=-\frac{\omega }{c}{\varpi }^2,\hfill \end{array}$$(6)

$$\begin{array}{cc}& h_{\mathrm{r}}=\frac{\rho }{\mathrm{\Delta }},h_{\theta }=\rho ,h_{\varphi }=\varpi =\frac{\mathrm{\Sigma }}{\rho }sin\theta .\hfill \end{array}$$(7)

2.2. Maxwell’s equations

Starting with Maxwell’s equations in curved spacetime (Thorne & Macdonald 1982; Gourgoulhon 2012), we applied the axisymetric and stationary conditions, leading to

$$\begin{array}{cc}& \mathbf{\nabla }·\mathit{E}=4\pi {\rho }_{\mathrm{e}},\hfill \end{array}$$(8)

$$\begin{array}{cc}& \mathbf{\nabla }·\mathit{B}=0,\hfill \end{array}$$(9)

$$\begin{array}{cc}& \mathbf{\nabla }\times \left(h\mathit{E}\right)=(\mathit{B}·\mathrm{\nabla }\frac{\omega }{c})\varpi \widehat{\varphi },\hfill \end{array}$$(10)

$$\begin{array}{cc}& \mathbf{\nabla }\times \left(h\mathit{B}\right)=\frac{4\pi h}{c}\mathit{J}-(\mathit{E}·\mathbf{\nabla }\frac{\omega }{c})\varpi \widehat{\varphi }.\hfill \end{array}$$(11)

The quantities above are all given in the ZAMO frame. We split all vector fields into a poloidal and toroidal component. The poloidal component of the magnetic field, B_p, can be expressed in terms of the magnetic flux function, A:

$$\begin{array}{c}\hfill {\mathit{B}}_{\mathrm{p}}=\mathbf{\nabla }\times \left(\frac{A}{\varpi }\widehat{\varphi }\right).\end{array}$$(12)

Faraday’s law gives us the electric field:

$$\begin{array}{c}\hfill \mathbf{\nabla }\times (h\mathit{E}-\frac{\omega }{c}\mathrm{\nabla }A)=0.\end{array}$$(13)

If we introduce the electric potential, Φ, we can write the electric field, E, as a function of the magnetic flux and the electric potential as

$$\begin{array}{c}\hfill h\mathit{E}=\frac{\omega }{c}\mathbf{\nabla }A-\mathbf{\nabla }\mathrm{\Phi }.\end{array}$$(14)

Finally, we write the condition of ideal MHD for infinite conductivity:

$$\begin{array}{c}\hfill \mathit{E}+\frac{\mathit{V}\times \mathit{B}}{c}=0.\end{array}$$(15)

V is the fluid velocity relative to the Eulerian observer. It should be noted that the symbol ∇ stands for the covariant derivative on a space hypersurface. The expressions of the components can be found in Appendix A.

2.3. Equations of motion

Using 3+1 formalism in the Kerr metric, in a steady state, we derived the equation of mass conservation, the Euler equation, and the energy conservation, respectively

$$\begin{array}{cc}& \mathbf{\nabla }·\left({\rho }_0\gamma h\mathit{V}\right)=0,\hfill \end{array}$$(16)

$$\begin{array}{cc}& {\rho }_0\gamma (\mathit{V}·\mathbf{\nabla })\left(\gamma \xi \mathit{V}\right)+{\rho }_0\xi {\gamma }^2[c^2\mathbf{\nabla }lnh+\frac{\varpi \omega V^{\widehat{\varphi }}}{h}\mathbf{\nabla }ln\omega ]\hfill \\ \hfill & +\mathbf{\nabla }P={\rho }_{\mathrm{e}}\mathit{E}+\frac{\mathit{J}\times \mathit{B}}{c},\hfill \end{array}$$(17)

$$\begin{array}{cc}& {\gamma }^2{\rho }_0\xi c[\mathit{V}·\mathbf{\nabla }ln(\gamma \xi h)+\frac{\varpi \omega V^{\widehat{\varphi }}}{h{c}^2}\mathit{V}·\mathbf{\nabla }ln\omega ]=\frac{\mathit{J}·\mathit{E}}{c},\hfill \end{array}$$(18)

where

$$\begin{array}{c}\hfill \gamma \equiv {(1-\frac{{\mathit{V}}^2}{c^2})}^{-1/2},\end{array}$$(19)

is the Lorenz factor, ρ₀ is the rest mass density, and ξc² is the specific enthalpy measured in the co-moving frame of the flow. Furthermore, $V^{\widehat{\varphi }}$ is the toroidal component of the speed of the flow as measured by the ZAMOs. The specific enthalpy is given by

$$\begin{array}{c}\hfill \xi =1+\frac{\mathrm{\Gamma }}{\mathrm{\Gamma }-1}\frac{P}{{\rho }_0{c}^2},\end{array}$$(20)

with Γ being the polytropic index. In this work, since the flow is extremely relativistic, the polytropic index is equal to 4/3.

The first law of thermodynamics can be obtained by projecting the conservation of energy – momentum tensor along the fluid four-velocity u_β∇_αT^αβ = 0. The contribution of the electromagnetic field is zero, since we have ideal MHD. Using the conditions of steady state and axisymmetry we get

$$\begin{array}{c}\hfill {\rho }_0{c}^2{\mathit{V}}_{\mathrm{p}}·\mathbf{\nabla }\xi ={\mathit{V}}_{\mathrm{p}}·\mathbf{\nabla }P.\end{array}$$(21)

The last equation of motion will be the wind equation, which is derived in a separate section further below.

2.4. Constants of motion

As was stated above, the flow is axisymetric and steady; that is, $\frac{\partial }{\partial \varphi }=0$ and $\frac{\partial }{\partial t}=0$. In this case, the MHD equations in a general relativity frame can be partially integrated so that we can find several constants of motion (Beskin 2009). These constants, including the magnetic field, have already been deduced in previous work such as Cayatte et al. (2014). This type of flow is governed by a function, A, that defines the geometry of the magnetic flux surfaces, so field lines and first integrals are lines of constant magnetic flux, A. Using Eq. (16) and the axisymetric condition we can write

$$\begin{array}{c}\hfill 4\pi {\rho }_0\gamma h{\mathit{V}}_{\mathrm{p}}=\mathbf{\nabla }\times \left(\frac{\mathrm{\Psi }}{\varpi }\widehat{\varphi }\right),\end{array}$$(22)

where Ψ is the mass flux.

Using Eqs. (12), (15), and (22) we get that Ψ is constant on surfaces of constant A on which the corresponding streamlines and field lines are roped, V_p‖B_p. If we define ${\mathrm{\Psi }}_A=\frac{\mathrm{d}\mathrm{\Psi }}{\mathrm{d}A}$, which is a function of A, then we can write

$$\begin{array}{c}\hfill {\mathit{V}}_{\mathrm{p}}=\frac{{\mathrm{\Psi }}_A}{4\pi {\rho }_0\gamma h}{\mathit{B}}_{\mathrm{p}}.\end{array}$$(23)

Starting with Eq. (14), we can show that surfaces of constant electric potential are also surfaces of constant magnetic flux, Φ = Φ(A). So using the poloidal components of Eq. (15) we get the law of iso-rotation,

$$\begin{array}{c}\hfill \mathrm{\Omega }=\omega +\frac{h{V}^{\widehat{\varphi }}}{\varpi }-\frac{{\mathrm{\Psi }}_A{B}^{\widehat{\varphi }}}{4\pi {\rho }_0\gamma \varpi }.\end{array}$$(24)

$\mathrm{\Omega }(A)\equiv c\frac{\mathrm{d}\mathrm{\Phi }}{\mathrm{d}A}$ is the iso-rotation frequency, which is constant along each magnetic flux line.

Moving on, we can get the conservation of angular momentum flux, L(A), which is derived by integrating the Euler equation in the toroidal direction:

$$\begin{array}{c}\hfill L=\varpi (\gamma \xi V^{\widehat{\varphi }}-\frac{h{B}^{\widehat{\varphi }}}{{\mathrm{\Psi }}_A}).\end{array}$$(25)

Then we have the energy conservation, which derives directly from Eq. (18), substituting J from Eq. (11) using Eq. (14) and the definition of Ω,

$$\begin{array}{c}\hfill \mathcal{E}-L\omega =\gamma \xi h{c}^2-\frac{h\varpi (\mathrm{\Omega }-\omega )}{{\mathrm{\Psi }}_A}{B}^{\widehat{\varphi }}.\end{array}$$(26)

Finally, starting from the first law of thermodynamics Eq. (21) and assuming the equation of state Eq. (20), we get the polytropic law

$$\begin{array}{c}\hfill P=Q(A){\rho }_0^{\mathrm{\Gamma }},\end{array}$$(27)

where Q is a function of the magnetic flux, A, and remains constant along each flux line. This equation actually replaces Eq. (18) in our equations of motion.

2.5. Fields and velocities

Using the integrals of motion, we can express the Lorentz factor, and the toroidal components of the velocity and the magnetic field as functions of these integrals and the poloidal components. So,

$$\begin{array}{cc}& \gamma h\xi =\frac{M^2(\mathcal{E}-L\omega )-h^2(\mathcal{E}-L\mathrm{\Omega })}{(M^2-h^2)c^2+{\varpi }^2{(\mathrm{\Omega }-\omega )}^2},\hfill \end{array}$$(28)

$$\begin{array}{cc}& B^{\widehat{\varphi }}=\frac{{\mathrm{\Psi }}_A}{h\varpi }\frac{L{h}^2{c}^2-{\varpi }^2(\mathrm{\Omega }-\omega )(\mathcal{E}-L\omega )}{(M^2-h^2)c^2+{\varpi }^2{(\mathrm{\Omega }-\omega )}^2},\hfill \end{array}$$(29)

$$\begin{array}{cc}& V^{\widehat{\varphi }}=\frac{h}{\varpi }\frac{M^{2}L{c}^2-(\mathcal{E}-L\mathrm{\Omega }){\varpi }^2(\mathrm{\Omega }-\omega )}{M^2(\mathcal{E}-L\omega )-h^2(\mathcal{E}-L\mathrm{\Omega })},\hfill \end{array}$$(30)

where M is the Alfvén Mach number, defined as

$$\begin{array}{c}\hfill M^2=\frac{4\pi h^2{\rho }_0\xi {\gamma }^2{V}_{\mathrm{p}}^2}{B_{\mathrm{p}}^2}=\frac{\xi {\mathrm{\Psi }}_A^2}{4\pi {\rho }_0}.\end{array}$$(31)

It is more convenient for later use to study quantities with the same denominator. Thus, we use $H=h\gamma \xi \frac{V^{\widehat{\varphi }}}{c}$

$$\begin{array}{c}\hfill H=h\gamma \xi \frac{V^{\widehat{\varphi }}}{c}=\frac{h}{c\varpi }\frac{M^{2}L{c}^2-(\mathcal{E}-L\mathrm{\Omega }){\varpi }^2(\mathrm{\Omega }-\omega )}{(M^2-h^2)c^2+{\varpi }^2{(\mathrm{\Omega }-\omega )}^2}.\end{array}$$(32)

We also use the symbol Z = γhξc² in order to make the notation a bit “lighter.” At the Alfvén transition surface, both the numerators and the denominators of the above Eqs. (28), (29), (32) must be zero. So, we get

$$\begin{array}{c}\hfill M_a^2=h_a^2[1-\frac{{\varpi }_a^2{(\mathrm{\Omega }-{\omega }_a)}^2}{h_a^2{c}^2}].\end{array}$$(33)

Setting the numerator of Eq. (28) to zero, we get

$$\begin{array}{c}\hfill \frac{{\varpi }_a^2{(\mathrm{\Omega }-{\omega }_a)}^2}{h_a^2{c}^2}=\frac{L(\mathrm{\Omega }-{\omega }_a)}{\mathcal{E}-L{\omega }_a}.\end{array}$$(34)

Both of these equations must be satisfied at the Alfvén transition surface. The numerators of Eqs. (29) and (32) are also zero when Eq. (34) is satisfied, since they are just a linear combination of Eqs. (33) and (34). The a index refers to the quantities measured in the Alfvén surface.

2.6. Forces

In this section, we present the force densities in general and then their components as we project them perpendicular to the magnetic field onto the ∇A/|∇A| vector. By doing so, we may calculate the contribution of each term and study the significance of each force to the collimation of the outflow.

The centrifugal force density is defined as

$$\begin{array}{c}\hfill {\mathit{f}}_{\mathrm{cen}}={\rho }_0\xi {\gamma }^2{\left(V^{\widehat{\varphi }}\right)}^2\mathrm{\nabla }ln\varpi .\end{array}$$(35)

The inertia force density, f_in, is defined as

$$\begin{array}{c}\hfill {\mathit{f}}_{\mathrm{in}}=-{\rho }_0\gamma (\mathit{V}·\mathbf{\nabla })\left(\gamma \xi \mathit{V}\right)-{\mathit{f}}_{\mathrm{cen}}.\end{array}$$(36)

The electromagnetic force density is given by

$$\begin{array}{c}\hfill {\mathit{f}}_{\mathrm{electromagnetic}}={\rho }_{\mathrm{e}}\mathit{E}+\frac{\mathit{J}\times \mathit{B}}{c}.\end{array}$$(37)

We can split the terms even more and separate the electric field and the magnetic field parts. So, the electric field is

$$\begin{array}{c}\hfill {\mathit{f}}_{\mathrm{elec}}={\rho }_{\mathrm{e}}\mathit{E},\end{array}$$(38)

where the electron density, ρ_e, can be replaced by Eq. (8). The magnetic field force density can be written as

$$\begin{array}{c}\hfill {\mathit{f}}_{\mathrm{mag}}=\frac{\mathrm{\nabla }\times \left(h\mathit{B}\right)\times \mathit{B}}{4\pi h}+\frac{\mathit{E}·\mathbf{\nabla }\omega }{4\pi ch}\mathbf{\nabla }A.\end{array}$$(39)

We can furthermore separate the magnetic field terms in f_pol, which is the contribution of the poloidal magnetic field term, B_p, and f_tor, the contribution of the toroidal magnetic field term, $B^{\widehat{\varphi }}$.

$$\begin{array}{c}\hfill {\mathit{f}}_{\mathrm{pol}}=\frac{\mathrm{\nabla }\times \left(h{\mathit{B}}_{\mathbf{p}}\right)\times {\mathit{B}}_{\mathbf{p}}}{4\pi h}\end{array}$$(40)

and the toroidal magnetic force density,

$$\begin{array}{c}\hfill {\mathit{f}}_{\mathrm{tor}}=\frac{\mathrm{\nabla }\times \left(h{B}^{\widehat{\varphi }}\widehat{\varphi }\right)\times B^{\widehat{\varphi }}\widehat{\varphi }}{4\pi h}.\end{array}$$(41)

The gravitational force density as well as the Lense-Thirring one are written as

$$\begin{array}{c}\hfill {\mathit{f}}_{\mathrm{gr}}=-{\rho }_0\xi {\gamma }^2{c}^2\mathbf{\nabla }lnh\end{array}$$(42)

and

$$\begin{array}{c}\hfill {\mathit{f}}_{\mathrm{lt}}=-{\rho }_0\xi {\gamma }^2\frac{\varpi \omega V^{\widehat{\varphi }}}{h}\mathbf{\nabla }ln\omega ,\end{array}$$(43)

and finally the pressure force is

$$\begin{array}{c}\hfill {\mathit{f}}_{\mathrm{p}}=-\mathbf{\nabla }P.\end{array}$$(44)

Projecting the forces onto the ∇A/|∇A| vector, and after some algebra, we get

$$\begin{array}{cc}& f_{\mathrm{in}\perp }=\frac{{\mathrm{\Psi }}_A{(\mathrm{\nabla }A)}^2}{4\pi h}\mathrm{\nabla }\left(\frac{M^2}{{\mathrm{\Psi }}_A{\varpi }^{2}h}\right)·\frac{\mathrm{\nabla }A}{|\mathrm{\nabla }A|}-\frac{{\mathrm{\Psi }}_A^{2}Z}{4\pi M^{2}h{c}^2}\mathrm{\nabla }\left(\frac{Z}{h}\right)·\frac{\mathrm{\nabla }A}{|\mathrm{\nabla }A|}\hfill \\ \hfill & +\frac{{\mathrm{\Psi }}_A^2{c}^2\xi }{4\pi M^2}\mathrm{\nabla }\xi ·\frac{\mathrm{\nabla }A}{|\mathrm{\nabla }A|}+\frac{{\mathrm{\Psi }}_A^{2}H}{4\pi M^2{h}^2}\mathrm{\nabla }H·\frac{\mathrm{\nabla }A}{|\mathrm{\nabla }A|}\hfill \\ \hfill & -\frac{{\mathrm{\Psi }}_A^2{H}^2}{4\pi M^2{h}^2}\mathrm{\nabla }lnh·\frac{\mathrm{\nabla }A}{|\mathrm{\nabla }A|}+\frac{M^2{(\mathrm{\nabla }A)}^2}{4\pi {\varpi }^2{h}^2|\mathrm{\nabla }A|}{\mathrm{\nabla }}^{2}A,\hfill \end{array}$$(45)

$$\begin{array}{cc}& f_{\mathrm{cen}\perp }=\frac{{\mathrm{\Psi }}_A^2{H}^2}{4\pi M^2{h}^2}\mathrm{\nabla }ln\varpi ·\frac{\mathrm{\nabla }A}{|\mathrm{\nabla }A|}\hfill \end{array}$$(46)

$$\begin{array}{cc}& f_{\mathrm{elec}\perp }=\frac{(\mathrm{\Omega }-\omega ){(\mathrm{\nabla }A)}^2}{ch|\mathrm{\nabla }A|}{\rho }_{\mathrm{e}}\hfill \end{array}$$(47)

$$\begin{array}{cc}& f_{\mathrm{pol}\perp }=-\frac{{(\mathrm{\nabla }A)}^2}{4\pi h|\mathrm{\nabla }A|}[\mathrm{\nabla }\left(\frac{h}{\varpi }\right)·\mathrm{\nabla }A+\frac{h}{{\varpi }^2}{\mathrm{\nabla }}^{2}A]\hfill \end{array}$$(48)

$$\begin{array}{cc}& f_{\mathrm{tor}\perp }=-\frac{B^{\widehat{\varphi }}}{4\pi \varpi h|\mathrm{\nabla }A|}\mathrm{\nabla }(\varpi h{B}^{\widehat{\varphi }})·\mathrm{\nabla }A\hfill \end{array}$$(49)

$$\begin{array}{cc}& f_{\mathrm{p}\perp }=-\frac{\mathrm{\nabla }P·\mathrm{\nabla }A}{|\mathrm{\nabla }A|}\hfill \end{array}$$(50)

$$\begin{array}{cc}& f_{\mathrm{gr}\perp }=-\frac{{\mathrm{\Psi }}_A^2{Z}^2}{4\pi M^2{c}^2{h}^2}\mathrm{\nabla }lnh·\frac{\mathrm{\nabla }A}{|\mathrm{\nabla }A|}\hfill \end{array}$$(51)

$$\begin{array}{cc}& f_{\mathrm{lt}\perp }=-\frac{{\mathrm{\Psi }}_A^2\omega \varpi }{4\pi M^{2}c{h}^3}ZH\mathrm{\nabla }ln\omega ·\frac{\mathrm{\nabla }A}{|\mathrm{\nabla }A|}\hfill \end{array}$$(52)

3. Model details and equations

3.1. Angular expansion

As was mentioned in the previous section, the constants of motions are functions of the magnetic flux, A. Moreover, both MHD equations and the Kerr metric constitute a system of nonlinear equations that cannot be solved analytically. The typical method of separation of variables that is used in Newtonian flows cannot be applied in this particular case, the reason being the complexity of the Kerr metric. Having said that, we approach this problem, describing the jet close to its symmetry axis, by expanding all variables with respect to the polar angle around the rotation axis in the form w(r, θ) = w₀(r)+w₁(r)sin²θ. The quantities with the zero index, imply the value of the function exactly on the axis, and the ones with an index of one will give the value of the function away from the axis at an angle, θ. It is also equivalent to expand all quantities over different flux lines, since the magnetic flux will be written as a function of both the distance, r, and the polar angle, θ, such as A = f(r)sin²θ, with f(r) being a function of r. So, equivalently, each quantity, w, can be written as w(r, A) = w₀(r)+w₁′(r)A. The same as above, the quantities with the zero index are over the axis and an index of 1 gives the value of the function over a magnetic flux line, A.

We begin by writing all physical quantities in a dimensionless form. So let us define the dimensionless spherical radius, R, as

$$\begin{array}{c}\hfill R=\frac{r}{{\varpi }_0}.\end{array}$$(53)

In this model, we arbitrarily selected ϖ₀ to be equal to 10³ Schwarzschild radii, r_s. We used the following dimensionless parameters,

$$\begin{array}{c}\hfill \mu =\frac{r_{\mathrm{s}}}{{\varpi }_0},l=\frac{a}{{\varpi }_0},\end{array}$$(54)

and we expanded in the first order of α all the physical quantities, where α is the dimensionless flux function,

$$\begin{array}{c}\hfill \alpha =\frac{A}{B_0{\varpi }_0^2}.\end{array}$$(55)

We can now define the radial part of the magnetic flux, the function, f, as a function of the dimensionless radius, R,

$$\begin{array}{c}\hfill f(R)=\frac{B_0{\varpi }_0^2}{2}\frac{R^2+l^2}{G{(R)}^2},\end{array}$$(56)

with G(R) the dimensionless cylindrical radius. So, the dimensionless magnetic flux can now be written as

$$\begin{array}{c}\hfill \alpha =\frac{R^2+l^2}{2G{(R)}^2}{sin}^2\theta .\end{array}$$(57)

We can now expand all the physical quantities for small angles, θ, or for small values of α, respectively. We kept only terms up to the first order in α so we can express the quantities as ξ(R, α) = ξ₀(R)+ξ₁(R)α. Now, we can write the quantities of the metric such as the line elements, the lapse function, and the angular velocity in that form as follows:

$$\begin{array}{c}\hfill h=\sqrt{1-\frac{\mu R}{R^2+l^2}}(1-\frac{\mu R{l}^2{G}^2}{{(R^2+l^2)}^3}\alpha ).\end{array}$$(58)

The terms on the axis, where θ = 0, are denoted with a zero index and the first-order terms on α with an index of one. So, the lapse function on the z axis, h(R, α = 0), will be simply written as h₀. As is obvious, the terms on the z axis are independent of α. In the same manner, the angular velocity is ω = ω₀ + ω₁ α, with

$$\begin{array}{c}\hfill {\omega }_0=\frac{clR\mu }{{\varpi }_0{(R^2+l^2)}^2},\\ \hfill {\omega }_1={\omega }_0\frac{2l^2{G}^2{h}_0^2}{{(R^2+l^2)}^2}.\end{array}$$(59)

The line elements of the metric take the form

$$\begin{array}{cc}\hfill h_{\mathrm{r}}& =\frac{1}{h_0}(1-\frac{l^2{G}^2}{{(R^2+l^2)}^2}\alpha ),\hfill \\ \hfill h_{\theta }& ={\varpi }_0\sqrt{R^2+l^2}(1-\frac{l^2{G}^2}{{(R^2+l^2)}^2}\alpha ),\hfill \\ \hfill h_{\varphi }& =\varpi ={\varpi }_0\sqrt{2}G\sqrt{\alpha }.\hfill \end{array}$$(60)

We also introduce the expansion factor, F, that expresses the geometry of the flux lines:

$$\begin{array}{c}\hfill F=\frac{\mathrm{d}lnf}{\mathrm{d}lnR}=2(\frac{R^2}{R^2+l^2}-\frac{\mathrm{d}lnG}{\mathrm{d}lnR}).\end{array}$$(61)

3.2. Free integrals

We now expand to the first order in the magnetic flux the integrals of motion, the mass to magnetic flux ratio, ${\mathrm{\Psi }}_A^2$, the isorotation frequency, Ω, the total angular momentum to mass flux ratio, L, and Q,

$$\begin{array}{c}\hfill {\mathrm{\Psi }}_A^2={\mathrm{\Psi }}_{A_0}^2+{\mathrm{\Psi }}_{A_1}^{2}A.\end{array}$$(62)

Alternatively, we can express Ψ_A through dimensionless quantities, ψ₀ and ψ₁, as is shown below:

$$\begin{array}{c}\hfill {\mathrm{\Psi }}_A^2=\frac{4\pi B_0^2}{c^2}({\psi }_0+{\psi }_1\alpha ).\end{array}$$(63)

The total energy is given as

$$\begin{array}{c}\hfill \frac{\mathcal{E}}{c^2}={\epsilon }_0+{\epsilon }_1\alpha ,\end{array}$$(64)

where ε₀ = γ₀ξ₀h₀. Although the isorotation frequency will be expanded, the first-order term vanishes. The total angular momentum to mass flux ratio, L, is written as

$$\begin{array}{c}\hfill L=L_1\alpha ,\end{array}$$(65)

and finally Q is written as

$$\begin{array}{c}\hfill Q=Q_0+Q_{1}A.\end{array}$$(66)

Expressing Q through dimensionless quantities, q₀ and q₁, we can type

$$\begin{array}{c}\hfill Q=\frac{\mathrm{\Gamma }-1}{\mathrm{\Gamma }}\frac{c^{2\mathrm{\Gamma }}}{{\left(B_0^2{\psi }_0\right)}^{\mathrm{\Gamma }-1}}(q_0+q_1\alpha ).\end{array}$$(67)

3.3. Wind equation

We started with the definition of the Lorenz factor and analyzed the speed of the flow in the poloidal and toroidal component, $V^2=V_{\mathrm{p}}^2+{\left(V^{\widehat{\varphi }}\right)}^2$. We substituted the poloidal and toroidal terms by using Eqs. (23) and (32), respectively, and using the expression of the poloidal magnetic field Eq. (12) we get

$$\begin{array}{c}\hfill \frac{Z^2}{c^2}={(c\xi h)}^2+{(cH)}^2+{\left(\frac{M^2\mathrm{\nabla }A}{\varpi {\mathrm{\Psi }}_A}\right)}^2.\end{array}$$(68)

By expanding the above equation in the first order on α, we can calculate the dimensionless cylindrical radius, G, and also the first-order term of the specific enthalpy, ξ₁,

$$\begin{array}{cc}& G=\frac{M_0}{\sqrt[4]{4\pi {\psi }_0({\epsilon }_0^2-h_0^2{\xi }_0^2)}},\hfill \end{array}$$(69)

$$\begin{array}{cc}& {\xi }_1=\frac{N_{{\xi }_1}}{D_{{\xi }_1}},\hfill \end{array}$$(70)

$$\begin{array}{cc}& N_{{\xi }_1}=-({\xi }_0-1){\xi }_0[2R^2(2q_1{\psi }_0+q_0(\mathrm{\Gamma }-1){\psi }_1)M_0^4\hfill \\ \hfill & +q_0(\mathrm{\Gamma }-1){\psi }_0{F}^2{G}^2{h}_0^2{M}_0^4\hfill \\ \hfill & +8\pi q_0{R}^2(\mathrm{\Gamma }-1){\psi }_0^2{G}^4(H_1^2-2{\epsilon }_0{\zeta }_1+2h_0{h}_1{\xi }_0^2)],\hfill \end{array}$$(71)

$$\begin{array}{cc}& D_{{\xi }_1}=4q_0{R}^2{\psi }_0[4\pi (\mathrm{\Gamma }-1){\psi }_0{G}^4{h}_0^2({\xi }_0-1){\xi }_0^2\hfill \\ \hfill & +M_0^4(1-\mathrm{\Gamma }+(\mathrm{\Gamma }-2){\xi }_0)].\hfill \end{array}$$(72)

3.4. Transfield equation

Starting from the Euler Eq. (17), we expanded each term over the first order of the magnetic flux and by projecting the result over the gradient of the flux, ∇A, we derived the transfield equation. It is a second-order differential equation, and by solving it numerically we can calculate all the quantities and forces of the outflow with respect to the dimensionless radius, R. With a little algebra instead of the second-order differential equation, we can equivalently get a system of first-order differential equations, as is shown in Appendix B. Although we can only get numerical solutions of the transfield equation, is it very important to do so, because this way we can self-consistently determine the geometry of the flux lines, without needing to make arbitrary assumptions about the shape of the flow.

3.5. Velocities and magnetic fields

Furthermore, we can express the velocity and the magnetic field components over the magnetic flux. So, after expanding and by applying the stationary and axissymetry criterion starting from Eq. (12), we get the radial and polar terms of the poloidal magnetic field, respectively,

$$\begin{array}{c}\hfill B^{\widehat{r}}=\frac{B_0}{G^2}[1-2\frac{R^4+l^2{R}^2+l^{2}R\mu }{{(R^2+l^2)}^3}{G}^2\alpha ]\end{array}$$(73)

and

$$\begin{array}{c}\hfill B^{\widehat{\theta }}=-\frac{B_{0}F{h}_0}{\sqrt{2}GR}\sqrt{\alpha }.\end{array}$$(74)

Using Eqs. (29) and (30), we derived the expressions for the toroidal magnetic field and the toroidal speed, respectively,

$$\begin{array}{c}\hfill B^{\widehat{\varphi }}=-\frac{2B_0\sqrt{\pi {\psi }_0}[L_1{h}_0^2-{\epsilon }_0{\varpi }_1^2({\mathrm{\Omega }}_0-{\omega }_0)]}{h_0{\varpi }_1(h_0^2-M_0^2)}\sqrt{\alpha }\end{array}$$(75)

and

$$\begin{array}{c}\hfill V^{\widehat{\varphi }}=\frac{c{h}_0[{\epsilon }_0{\varpi }_1^2({\mathrm{\Omega }}_0-{\omega }_0)-L_1{M}_0^2]}{{\epsilon }_0{\varpi }_1(h_0^2-M_0^2)}\sqrt{\alpha }.\end{array}$$(76)

Using the relation connecting the magnetic fields with the Mach number and the velocities Eq. (31), we can write down the equations of the radial and polar speed,

$$\begin{array}{c}\hfill V^{\widehat{r}}=\frac{c{M}_0^2}{2\sqrt{\pi {\psi }_0}{\epsilon }_0{G}^2}(1+\frac{V^{\widehat{r_1}}}{V^{\widehat{r_0}}}\alpha )\end{array}$$(77)

and

$$\begin{array}{c}\hfill V^{\widehat{\theta }}=-\frac{cF{h}_0{M}_0^2}{2\sqrt{2\pi {\psi }_0}{\epsilon }_{0}RG}\sqrt{\alpha }.\end{array}$$(78)

The expressions of $V^{\widehat{r_0}}$ and $V^{\widehat{r_1}}$ are given in Appendix B.

We also present the radial and polar components of the electric field, $E^{\widehat{r}}$ and $E^{\widehat{\theta }}$:

$$\begin{array}{cc}& E^{\widehat{r}}=-\frac{F({\mathrm{\Omega }}_0-{\omega }_0)}{\mathit{cR}}\alpha ,\hfill \end{array}$$(79)

$$\begin{array}{cc}& E^{\widehat{\theta }}=-\frac{\sqrt{2}({\mathrm{\Omega }}_0-{\omega }_0)}{c{\varpi }_0{h}_{0}G}\sqrt{\alpha }.\hfill \end{array}$$(80)

Using Maxwell’s Eq. (8), we can calculate the charge density. We will write down the zero-order term only, since the first order vanishes after the expansion takes place:

$$\begin{array}{c}\hfill {\rho }_{\mathrm{e}}=-\frac{{\mathrm{\Omega }}_0-{\omega }_0}{2\pi h_0{G}^2}.\end{array}$$(81)

Finally, using the specific enthalpy Eq. (20) and the polytropic law Eq. (27), we may express the zero- and first-order terms of the rest mass density and the pressure:

$$\begin{array}{c}\hfill {\rho }_0={\rho }_{0_0}+{\rho }_{0_1}\alpha ,\end{array}$$(82)

and

$$\begin{array}{c}\hfill P=P_0+P_1\alpha .\end{array}$$(83)

The expressions of ρ₀₀, ρ₀₁ as well as P₀, P₁ are given in Appendix B.

4. Solutions

Although the model is meant to cover a wide variety of solutions, in this paper we present solutions that begin from a stagnation surface and cross the sound-critical point close to the event horizon, here at a radius of r = 2r_s. The system consists of three nonlinear differential equations as functions of the radius, R, the expansion function, F, and a function Δ that depend on specific enthalpy ξ which will be solved numerically. The system depends on eleven more parameters: q₀,  q₁, which define the pressure on the axis and its first-order term; ψ₀,  ψ₁, the dimensionless mass to magnetic flux ratio; ε₀,  ε₁, energy on the axis and on each field line, respectively; L₁, the angular momentum flux; the scaling factor, μ; the dimensionless spin of the black hole, l; the polytropic index, Γ; and the isorotation frequency on the axis, Ω₀. In this paper, the polytropic index is set to 4/3 and we assume a maximally rotating black hole. What is more, two of the above parameters can be evaluated from the regularity conditions as the flow crosses the critical point. We have chosen to evaluate q₁, ψ₁ since this way we have first-order equations to solve (see Appendix C for more details). The solution below has an asymptotic Lorentz factor, γ ≃ 70. We have chosen parameters in such a manner in order to obtain a solution with a super Alfvénic speed right from the start. The flow starts at subsonic speed and becomes supersonic as it crosses the fast magnetosonic critical point at a distance of 2 Schwarzschild radii. We note that on the rotational axis the fast magnetosonic speed equals either the sound speed or the Alfvén speed. In this particular case it equals the sound speed. For this solution, we have chosen a value for the isorotation frequency ${\mathrm{\Omega }}_0=\frac{{\omega }_{0_{\mathrm{EH}}}}{2}$, where ω_0EH is the angular velocity of the ZAMOs at the event horizon. Furthermore, we assume that the black hole is maximally rotating, α = r_s/2, or equivalently l = μ/2.

In Fig. 1 we compare the velocity of the flow, γ₀u₀ (red line), with the sound speed, γ_su_s (green line), as well as the Alfvén speed γ_au_a (blue line), on a logarithmic scale. The flow (red line) starts from a stagnation surface at a distance of approximately 1.5r_s Schwarzschild radii, moving subsonically, as is shown in Fig. 1, until it intersects the fast magnetosonic point at a distance of two Schwarzschild radii, r = 2r_s, and then becomes superfast. The flow is well into the super Alfvénic regime from the start of the acceleration, and quickly becomes highly relativistic. Figure 2 zooms close to the stagnation surface up until the critical point crossing.

Fig. 1. Velocity of the flow, γ0u0, in red, the sound speed, γsus, in green, and the Alfvén speed, γaua, in blue color, on a logarithmic scale. The stagnation surface is close to the black hole and the critical point crossing is at a distance of two Schwarzschild radii, r = 2rs.

Fig. 2. Zoom-in on the stagnation surface in the subsonic area of the flow until the critical point. Same as in Fig. 1, the velocity of the flow, γ0u0, in red, the sound speed, γsus, in green, and the Alfvén speed, γaua, in blue.

Similarly, in Fig. 3 we can see that the specific enthalpy is rapidly decreasing even close to the starting point of the flow as the velocity and, respectively, the kinetic energy, rapidly increase, making the flow extremely relativistic, as is seen in Fig. 1. At a large distance where the flow has become cold, the kinetic energy has reached its maximum value. The flow is thermally driven, despite the fact that the magnetic field is present. This result agrees with previous works by Meliani et al. (2006) and Chantry et al. (2018). The radial velocity, V_r0, on the axis shown in Fig. 4 rapidly reaches a maximum value close to the speed of light. This was quite expected for a highly relativistic flow with a Lorentz factor close to 70. The analytical form of the velocity can be seen in Eq. (B.31).

Fig. 3. Specific enthalpy of the flow, ξ0, on a logarithmic plot. The blue and red lines are the parts of the specific enthalpy before and after the critical point, respectively. The speed of light is considered to be c = 1.

Fig. 4. Zeroth term of the radial speed, Vr0. The blue and red lines represent the part of the velocity before and after it crosses the critical point, respectively.

The toroidal component of the speed, V_ϕ1, plotted in Fig. 5 as seen from the ZAMO frame, is just the toroidal speed described in Eq. (79) over the flux function, $V_{{\varphi }_1}=\frac{V^{\widehat{\varphi }}}{\sqrt{\alpha }}$. The speed is negative very close to the stagnation surface and this has to do to with the chosen value of the integral, Ω₀. As was said at the beginning of this section, we assume that ${\mathrm{\Omega }}_0=\frac{{\omega }_{0_{\mathrm{EH}}}}{2}$, corresponding to the maximum energy extraction predicted for the Blandford & Znajek (1977) mechanism.

Fig. 5. Toroidal component, Vϕ1, of the speed flow. Due to the chosen value of the integral, Ω0, this leads to a negative toroidal speed close to the stagnation surface.

In Figs. 6 and 7 we present the radial and toroidal components of the magnetic field as they are written in Eqs. (B.29) and (B.28), respectively. At first glance, it might seem that the value of the toroidal field, arithmetically, is much higher than the value of the radial field, but we need to have in mind that the contribution of the toroidal field to each flux line is dependent on $\sqrt{\alpha }$. Due to the short region where the model is valid, something that we discuss in the next section, the value of $\sqrt{\alpha }$ significantly reduces the contribution of the toroidal field. Nevertheless, the high value of it close to the stagnation surface suggests that the field lines are over-collimated in this region, as is seen in Fig. 8.

Fig. 6. Radial component of the magnetic field on the axis.

Fig. 7. Toroidal component of the magnetic field.

Fig. 8. Field lines close to the stagnation surface. The lines are over-collimated at the beginning, which is in agreement with the high value of the toroidal magnetic field. All the distances are on a scale of a thousand Schwarzschild radii (1000 rs). Blue lines refer to the flow before the critical point, whereas the red ones refer to after it has crossed the critical point. $\varpi =\sqrt{R^2+l^2}sin\theta$ and $z=\sqrt{R^2+l^2}cos\theta$.

In Fig. 9 we plot, on a logarithmic scale, the force densities of the flow as we deduce them from the Euler equation in the transfield direction. As we can observe close to the critical point, and thus close to the black hole, the inertia including the gravitational force and the gradient of pressure force are the most significant ones. This changes quite fast, and a little bit further away from the black hole these forces decline rapidly. The most significant forces of the flow become the electric force and the force of the toroidal component of the magnetic field, f_tor. The flow is thermally driven but the electromagnetic force is the one responsible for the collimation of the jet. The poloidal magnetic force, f_pol, together with the centrifugal force, f_cen, remain the least significant forces throughout the whole range of the flow. The inertia force, f_in + f_gr, and the gradient of pressure, f_p, as well as the f_tor and f_elect remain competitive forces throughout the whole region of the flow. In Fig. 10 we compare these competitive groups. We did so in order to confirm our last claim that the centrifugal force and the force of the poloidal magnetic field contribute the least to the acceleration of the flow. As is indeed clear, those are the least important forces, even compared to the (f_in+f_p+f_gr) and (f_elect+f_tor) terms.

Fig. 9. Transfield components of force densities on a logarithmic scale. ftor and felect represent the forces of the toroidal magnetic field and the electric field, respectively. The contribution of the inertia plus the gravitational force (fin+fgr) as well as the gradient of pressure, fp, play a significant role closer to the event horizon. The poloidal and centrifugal components, fpol and fcen, respectively, do not play a significant role.

Fig. 10. Comparison of the poloidal magnetic force with the difference of the competitive forces on a logarithmic scale.

In Fig. 11 we compare the energy of the flow minus the energy on the axis over the magnetic flux. Although the flow is thermally driven, the magnetic acceleration (red line) plays a significant role in small angles near the axis. The green line shows the contribution of spacetime rotation. The energy contribution of the spacetime rotation, although present, does not play a significant role in accelerating the flow, which is something to be expected since the stagnation surface is already outside the ergosphere and the jet travels to great distances. As we can see, the magnetic force contributes significantly to the acceleration of the flow, but the thermal acceleration on the axis still outweighs the contribution of the magnetic one, overall. This is an inherent property of the model, which we will discuss in the last section. The geometry of the lines of the flow, starting from the stagnation surface and onward to large distances, are shown in Fig. 12.

Fig. 11. Energy of the flow minus the energy on the axis over the magnetic flux. The corresponding term of the magnetic field is shown in red. The green line shows the contribution of spacetime rotation. The kinetic and thermal terms are shown in blue.

Fig. 12. Lines of the flow. $\varpi =\sqrt{R^2+l^2}sin\theta$, $z=\sqrt{R^2+l^2}cos\theta$. The scale of the distance is a thousand Schwarzschild radii (1000 rs).

4.1. Region where the model is valid

Since this model relies on expanding all physical quantities with respect to the polar angle around the axis, it is quite reasonable to wonder what region this expansion is valid in. We assume a percentage of relative error, the same for all physical quantities, that the terms of the expansion are “allowed” to have and find the angle that corresponds to this condition. We do this for all the physical quantities and keep the smallest angle. We calculate the angle assuming the relative error at a distance of r = 2r_s, where the flow becomes supersonic. Having said that, we can express each physical quantity, such as special enthalpy, as ξ = ξ₀ + Δξ with Δξ = ξ₁α. The dimensionless magnetic flux, α, in Eq. (57) “hides” the dependence of the polar angle, θ, and thus we can calculate it. The relative error in this solution has been chosen as 5%, meaning $\frac{\mathrm{\Delta }\xi }{{\xi }_0}=5\%$, and is similar for the other physical quantities.

As is obvious from the Table 1, for a relative error of 5% the pressure restricts the expansion to being valid for an angle of 11.6°. Since this is the smallest of all other angles, this is considered to be the maximum angle where the expansion is valid. The angle that corresponds to the pressure gets bigger with distance, as we can see in Fig. 13, but still remains smaller from the rest, and thus remains the determinant factor that restricts the region of validity.

Fig. 13. Lines of the flow alongside the terminal lines defined by the maximum allowed field line. The terminal line is defined by an arbitrary percentage of error. While the percentage is the same for all physical quantities, in this case 5%, the angle where the model is valid takes different values depending on the physical quantity, and therefore corresponds to a different maximum field line. The lines, starting with the closest to the axis, are assigned to pressure, density, Mach number, radial velocity, and the radial magnetic field. The errors of speed, γu, and specific enthalpy remain significantly smaller than 5% in this region.

Since pressure allows for the smallest angle away from the axis, it is this physical quantity that defines the overall maximum terminal line of the flow. This line defines the boundary inside of which the approximation of small angles is valid. Although we can assume a different relative error for the expansion to be valid, the pressure always remains the physical quantity that defines the limit of the area of validity. Since the lines of the flow have a complex shape, the angle that defines the validity area does not remain constant as we move further away from the black hole. In any case, the pressure line (the magenta line in Fig. 13) always defines the maximum terminal line. All physical quantities inside this area are always below the 5% threshold that we have assumed. Comparing the validity region to a nonrotating model, the rotating one is less wide closer to the black hole, while the overall shapes of the lines are similar.

5. Conclusions

In this paper we have studied the formation and collimation of relativistic jets. This was done via a semi-analytical meridional self-similar model. This model is based on the classical work developed by Sauty & Tsinganos (1994) and the generalization by Meliani et al. (2006) in the Schwarzschild metric and Globus et al. (2014), Chantry et al. (2018) in the Kerr metric.

Specifically, we presented MHD solutions in the Kerr metric and we expanded the equations in the first order of the magnetic flux function. The system depends on a large number of parameters. Two of them can be calculated by applying the regularity condition at the critical point, as was mentioned in Sect. 4.

The numerical solution presented here starts from a stagnation surface very close to the event horizon and crosses the fast magnetosonic surface, becoming thus supersonic. The radial magnetic field of the flow was chosen to be relatively small, in order to begin with a super Alfvénic outflow. Although the energy extraction is due to both electromagnetic and thermal mechanisms, the electromagnetic field does not seem to be the dominant acceleration factor of the outflow. Indeed, the flow is mostly thermally driven, with a regime allowed for the solutions to be valid, just a few degrees away from the rotating axis. On the other hand, the magnetic field, and to be more specific the toroidal component, is the dominant factor responsible for the collimation of the jet. These results are in agreement with the previous works of Meliani et al. (2006), Globus et al. (2014), and Chantry et al. (2018). The last two works were also done in Kerr metric. That being said, as we can see from Fig. 11, along a flux line the magnetic acceleration contributes to the jet driving, but overall, comparing it with the total thermal acceleration on the axis, it is not significant. Since the model was based on expanding all physical quantities over small angles away from the rotational axis, and due to the fact that the magnetic terms, other than the radial component, B_r, do not possess a component on the axis, it seems quite inevitable that the thermal acceleration that is present on the axis is always the dominant term. It is possible to increase the contribution of the magnetic field to the acceleration of the flow by simply allowing for larger angles in the expansion. By doing so, we expand the region where the model is valid, but we decided to choose the more conservative approach by setting the relative error in all quantities to be less than 5%. We have also produced solutions for a nonrotating black hole. Comparing the results of the rotating case with the behavior of the model of a nonrotating black hole, these are very close, without significant differences and that is why we do not display them separately. The magnetic fields, the velocity, and the density exhibit similar behaviors in both the rotating and nonrotating models. This is something to be expected, since the region where the acceleration and in general the phenomena take place is quite far from the event horizon where general relativity is the dominant factor. Even the start of the flow, at the stagnation surface, is already far away from the event horizon, so the frame dragging effects do not alter the overall behavior of the model.

On the other hand, and contrary to the model of Chantry et al. (2018) in which there is external heating, in this model the pressure is governed by a polytropic law without the need to assume any heating. This model can be used as an initial condition on numerical simulations.

A wider variety of solutions crossing both Alfvén and sound critical points will be the subject of future work.

Acknowledgments

We would like to thank the referee Zakaria Meliani for his helpful comments.

References

Appendix A: Operators of vectors in Boyer-Lindquist coordinates in the Kerr Metric

The vector components are expressed in the orthonormal base, ϵ_{i = 1, 2, 3}. The natural base, e_{i = 1, 2, 3}, is associated with the Boyer-Lindquist coordinates. The transformation relation between e_i and ϵ_i can be given by the formula e_i = h_iϵ_i, ∀i = 1, 2, 3. A vector component in the orthonormal base can be written as $\mathit{V}=V^{\widehat{i}}{\mathit{ϵ}}_i$, with $V^{\widehat{i}}=V_{\widehat{i}}$. We note V_i, Vⁱ, the components in the natural base, and $V^{\widehat{i}},V_{\widehat{i}}$, the components if the ZAMO orthonormal base.

Next, we summarize the results:

$$\begin{array}{cc}\hfill {\mathit{e}}_i& =h_i{\mathit{ϵ}}_i\hfill \\ \hfill V_i& =h_i^2{V}^i=h_i{V}^{\widehat{i}}.\hfill \end{array}$$(A.1)

We note that for ϵ₃ = ϵ_ϕ in the main body of this paper we use $\widehat{\varphi }.$

Below, we write down the vectorial operators under the assumption of axisymmetry, ∂_ϕ = 0.

• The gradient vector, ∇, in the ZAMO orthonormal bases:

  $$\begin{array}{c}\hfill \mathbf{\nabla }={\displaystyle \sum _{i=1}^3}\frac{{\mathit{ϵ}}_i}{h_i}.\end{array}$$(A.2)

• The divergence of the vector, V:

  $$\begin{array}{c}\hfill \mathbf{\nabla }·\mathit{V}=\frac{1}{h_{\mathrm{r}}{h}_{\theta }{h}_{\varphi }}[{\partial }_{\mathrm{r}}\left(h_{\theta }{h}_{\varphi }{V}^{\widehat{r}}\right)+{\partial }_{\theta }\left(h_{\mathrm{r}}{h}_{\varphi }{V}^{\widehat{\theta }}\right)].\end{array}$$(A.3)

• The Laplace operator:

  $$\begin{array}{c}\hfill {\mathrm{\nabla }}^{2}A=\frac{1}{h_{\mathrm{r}}{h}_{\theta }{h}_{\varphi }}[\frac{1}{{\varpi }_0^2}{\partial }_R\left(\frac{h_{\theta }{h}_{\varphi }}{h_{\mathrm{r}}}{\partial }_{R}A\right)+{\partial }_{\theta }\left(\frac{h_{\mathrm{r}}{h}_{\varphi }}{h_{\theta }}{\partial }_{\theta }A\right)].\end{array}$$(A.4)

• The curl operator on a vector, V:

  $$\begin{array}{cc}\hfill \mathbf{\nabla }\times \mathit{V}=& \frac{h_k}{h_{\mathrm{r}}{h}_{\theta }{h}_{\varphi }}{ϵ}^{\mathit{ijk}}{\partial }_i\left(h_j{V}^{\widehat{j}}\right){\mathit{ϵ}}_k\hfill \\ \hfill \\ \hfill \mathbf{\nabla }\times \mathit{V}=& \left[\begin{array}{c}\frac{1}{h_{\theta }{h}_{\varphi }}{\partial }_{\theta }\left(\varpi V^{\widehat{\varphi }}\right)\\ -\frac{1}{h_{\mathrm{r}}{h}_{\varphi }{\varpi }_0}{\partial }_R\left(\varpi V^{\widehat{\varphi }}\right)\\ \frac{1}{h_{\mathrm{r}}{h}_{\theta }}[\frac{1}{{\varpi }_0}{\partial }_R\left(h_{\theta }{V}^{\widehat{\theta }}\right)-{\partial }_{\theta }\left(h_{\mathrm{r}}{V}^{\widehat{r}}\right)]\end{array}\right].\hfill \end{array}$$(A.5)

• The advection term:

  $$\begin{array}{c}\hfill (\mathit{V}·\mathbf{\nabla })\mathit{V}=V^{\alpha }{D}_{\alpha }{V}^{\beta }{\mathit{e}}_{\mathit{\beta }}.\end{array}$$(A.6)

The poloidal component of the above gives

$$\begin{array}{c}\hfill {[(\mathit{V}·\mathbf{\nabla })\mathit{V}]}_{\mathrm{p}}=\left[\begin{array}{c}\frac{V^{\widehat{r}}}{{\varpi }_0{h}_{\mathrm{r}}}{\partial }_R{V}^{\widehat{r}}+\frac{V^{\widehat{\theta }}}{h_{\theta }}{\partial }_{\theta }{V}^{\widehat{r}}+\frac{V^{\widehat{r}}{V}^{\widehat{\theta }}}{h_{\theta }}{\partial }_{\theta }ln(h_{\mathrm{r}})-\frac{{(V^{\widehat{\theta }})}^2}{h_{\mathrm{r}}}{\partial }_{\mathrm{r}}ln(h_{\theta })-\frac{{(V^{\widehat{\varphi }})}^2}{h_{\mathrm{r}}}{\partial }_{\mathrm{r}}ln(h_{\varphi })\\ \frac{V^{\widehat{r}}}{{\varpi }_0{h}_{\mathrm{r}}}{\partial }_R{V}^{\widehat{\theta }}+\frac{V^{\widehat{\theta }}}{h_{\theta }}{\partial }_{\theta }{V}^{\widehat{\theta }}+\frac{V^{\widehat{r}}{V}^{\widehat{\theta }}}{{\varpi }_0{h}_{\mathrm{r}}}{\partial }_{R}ln(h_{\theta })-\frac{{(V^{\widehat{r}})}^2}{h_{\theta }}{\partial }_{\theta }ln(h_{\mathrm{r}})-\frac{{(V^{\widehat{\varphi }})}^2}{h_{\theta }}{\partial }_{\theta }ln(h_{\varphi })\end{array}\right],\end{array}$$(A.7)

and finally the nonsymmetric term

$$\begin{array}{c}\hfill [(\mathit{B}·\mathbf{\nabla })\mathit{C}]=\left[\begin{array}{c}\frac{B^{\widehat{r}}}{{\varpi }_0{h}_{\mathrm{r}}}{\partial }_R{C}^{\widehat{r}}+\frac{B^{\widehat{\theta }}}{h_{\theta }}{\partial }_{\theta }{C}^{\widehat{r}}+\frac{B^{\widehat{r}}{C}^{\widehat{\theta }}}{h_{\theta }}{\partial }_{\theta }ln(h_{\mathrm{r}})-\frac{B^{\widehat{\theta }}{C}^{\widehat{\theta }}}{h_{\mathrm{r}}}{\partial }_{\mathrm{r}}ln(h_{\theta })-\frac{B^{\widehat{\varphi }}{C}^{\widehat{\varphi }}}{h_{\mathrm{r}}}{\partial }_{\mathrm{r}}ln(h_{\varphi })\\ \frac{B^{\widehat{r}}}{{\varpi }_0{h}_{\mathrm{r}}}{\partial }_R{C}^{\widehat{\theta }}+\frac{B^{\widehat{\theta }}}{h_{\theta }}{\partial }_{\theta }{C}^{\widehat{\theta }}+\frac{B^{\widehat{\theta }}{C}^{\widehat{r}}}{{\varpi }_0{h}_{\mathrm{r}}}{\partial }_{R}ln(h_{\theta })-\frac{B^{\widehat{r}}{C}^{\widehat{r}}}{h_{\theta }}{\partial }_{\theta }ln(h_{\mathrm{r}})-\frac{B^{\widehat{\varphi }}{C}^{\widehat{\varphi }}}{h_{\theta }}{\partial }_{\theta }ln(h_{\varphi })\\ \frac{B^{\widehat{r}}}{{\varpi }_0{h}_{\mathrm{r}}}{\partial }_R{C}^{\widehat{\varphi }}+\frac{B^{\widehat{\theta }}}{h_{\theta }}{\partial }_{\theta }{C}^{\widehat{\varphi }}+\frac{B^{\widehat{\varphi }}{C}^{\widehat{r}}}{{\varpi }_0{h}_{\mathrm{r}}}{\partial }_{R}ln(h_{\varphi })+\frac{B^{\widehat{\varphi }}{C}^{\widehat{\theta }}}{h_{\theta }}{\partial }_{R}ln(h_{\varphi })\end{array}\right].\end{array}$$(A.8)

Appendix B: Differential equations

The system of the two first-order ordinary differential equations is shown below:

$$\begin{array}{c}\hfill \{\begin{array}{cc}& \frac{d\mathrm{\Delta }}{\mathit{dR}}=\frac{N_1(R,\mathrm{\Delta },F)}{D(R,\mathrm{\Delta })}\hfill \\ \hfill & \frac{\mathit{dF}}{\mathit{dR}}=\frac{N_2(R,\mathrm{\Delta },F)}{D(R,\mathrm{\Delta })}\hfill \end{array}.\end{array}$$(B.1)

D(R, Δ) is the common denominator that becomes zero at some specific set of values (R_c, Δ_c, F_c).

We define $\frac{\mathit{dR}}{\mathit{ds}}=D(R,\mathrm{\Delta })$. So, the system of two first-order differential equations becomes a system of three equations using (B.1)

$$\begin{array}{c}\hfill \{\begin{array}{cc}& \frac{d\mathrm{\Delta }}{\mathit{ds}}=N_1(R(s),\mathrm{\Delta }(s),F(s))\hfill \\ \hfill & \frac{\mathit{dF}}{\mathit{ds}}=N_2(R(s),\mathrm{\Delta }(s),F(s))\hfill \\ \hfill & \frac{\mathit{dR}}{\mathit{ds}}=D(R(s),\mathrm{\Delta }(s))\hfill \end{array}\end{array}$$(B.2)

$$\begin{array}{c}\hfill N_1=\frac{(q_0{\mathrm{\Delta }}_0^{1-\mathrm{\Gamma }}+1)\{{\epsilon }_0^2[F(R^2+l^2)-2R^2]{\mathrm{\Delta }}_0^{2\mathrm{\Gamma }}-[F(R^2+l^2)-2R^2]h_0^2{(q_0{\mathrm{\Delta }}_0+{\mathrm{\Delta }}_0^{\mathrm{\Gamma }})}^2+h_{0}R(R^2+l^2){(q_0{\mathrm{\Delta }}_0+{\mathrm{\Delta }}_0^{\mathrm{\Gamma }})}^2{\partial }_{r_0}h\}}{{\mathrm{\Delta }}_0^{2\mathrm{\Gamma }-1}R(R^2+l^2)}\end{array}$$(B.3)

$$\begin{array}{c}\hfill N_2=\frac{4\pi R(R^2+l^2)h_0^2{G}^4}{(h_0^2-M_0^2)(l^2+R^2-\mu R)}({\mathrm{\Lambda }}_1{N}_1+{\mathrm{\Lambda }}_2{N}_{{\xi }_1}+{\mathrm{\Lambda }}_{3}D),\end{array}$$(B.4)

where we give the expressions of Λ₁, Λ₂, Λ₃ such as

$$\begin{array}{c}\hfill {\mathrm{\Lambda }}_1=\frac{(\mathrm{\Gamma }-1)q_0{\psi }_0{f}_{r_1}}{{\mathrm{\Delta }}_0^{\mathrm{\Gamma }+1}}+\frac{f_{r_1}{f}_{{\theta }_1}^2[1+q_0(2-\mathrm{\Gamma }){\mathrm{\Delta }}_0^{1-\mathrm{\Gamma }}]}{4\pi h_0^2{h}_{r_0}{\varpi }_1^2}\end{array}$$(B.5)

$$\begin{array}{c}\hfill {\mathrm{\Lambda }}_2=\frac{M_{1_b}}{2\pi h_0^2{G}^6}+\frac{{\psi }_0{f}_{{\theta }_1}^2(1+q_0{\mathrm{\Delta }}_0^{1-\mathrm{\Gamma }})}{M_0^2}-\frac{{\psi }_0{f}_{{\theta }_1}^2}{{\mathrm{\Delta }}_0}-\frac{f_{{\theta }_1}^3{M}_{1_b}}{2\pi h_0^2{\varpi }_1^3}+\frac{f_{{\theta }_1}^4{M}_{1_b}}{4\pi h_0^2{\varpi }_1^2}\end{array}$$(B.6)