Abstract

We present a new, approximate method for modelling the acceleration and collimation of relativistic jets in the presence of gravity. This method is self-similar throughout the computational domain where gravitational effects are negligible and, where significant, self-similar within a flux tube. These solutions are applicable to jets launched from a small region (e.g. near the inner edge of an accretion disc). As implied by earlier work, the flow can converge on to the rotation axis, potentially creating a collimation shock.

In this first version of the method, we derive the gravitational contribution to the relativistic equations by analogy with non-relativistic flow.

This approach captures the relativistic kinetic gravitational mass of the flowing plasma, but not that due to internal thermal and magnetic energies. A more sophisticated treatment, derived from the basic general relativistic magnetohydrodynamical equations, is currently being developed.

Here we present an initial exploration of parameter space, describing the effects the model parameters have on flow solutions and the location of the collimation shock. These results provide the groundwork for new, semi-analytic models of relativistic jets which can constrain conditions near the black hole by fitting the jet break seen increasingly in X-ray binaries.

1 INTRODUCTION

Despite decades of study, astrophysical jets are still an enigmatic component of our Universe. Understanding the details of their physics would advance fields as diverse as star formation, gamma-ray bursts, galaxy evolution and cluster dynamics. Jets in all of these objects probably share three basic ingredients: a source of energy and particles such as an accretion flow, magnetic fields created or carried inwards by the flow and rotation. However, the exact details of how these elements combine to launch and collimate jets in the observed systems are not at all clear.

Black holes are ideal test cases to understand jet physics because of their enormous range in mass, allowing us to study jet dynamics over an equally large range of time-scales. Accreting stellar mass black holes in a binary with a companion star (black hole binaries, BHBs) have jets that can undergo an entire launch and quench cycle on time-scales of months, which is related at least in part to the accretion state of the disc (Esin, McClintock & Narayan 1997; Meier 2001; Fender, Belloni & Gallo 2004; McClintock & Remillard 2006). In their supermassive cousins, active galactic nuclei (AGN), such transitions would be expected to occur on time-scales of millions to a billion times longer. To assess the extent to which such a mass scaling operates, we must have a better understanding of which elements of jet physics persist despite a dramatic range of environment and scales.

One key feature of jets regardless of source seems to be their ability to accelerate particles to highly relativistic energy. We observe the outcome of these accelerated particles directly via optically thin synchrotron emission in, e.g. AGN (Marscher & Gear 1985) and in BHBs (Fender 2001). When the jets are compact enough to become self-absorbed, the superposition of many synchrotron-emitting components along the length of the jet leads to a flat or slightly inverted total spectrum (Blandford & Konigl 1979) until the point where the jet emission becomes entirely optically thin. At this point the spectral index will steepen dramatically to α ∼ 0.5- 0.8 typically, where α is defined such that Fν ∝ ν−α. This break in the spectrum can be physically associated with the most compact region in the jets where particle acceleration is present. In low-luminosity AGN jets this break typically occurs in the GHz range (Ho 1999), and if jets are self-similar across the black hole mass range, in BHBs such a break would be predicted to occur in the infrared (IR) bands (Markoff, Falcke & Fender 2001; Heinz & Sunyaev 2003; Markoff et al. 2003). Because the spectral break also indicates the total radiative power of the jets and plays a key role in driving mass-scaling predictions such as the Fundamental Plane of black hole accretion (Merloni, Heinz & di Matteo 2003; Falcke, Körding & Markoff 2004; Plotkin et al. 2012), it is a pivotal parameter linking the inner jets (and assumedly conditions near their launch point) to the outer jets associated with strong radio through IR emission. Deriving this link and overall jet structure in a way consistent with magnetohydrodynamics (MHD) is the focus of this paper.

Constraints on the location of the jet break can come from both observations as well as from spectral fitting. Because the optical/IR bands are often dominated by a stellar companion or the accretion disc, the optically thick-to-thin break has been observed explicitly in only one BHB so far, GX 339-4 [Corbel & Fender 2002; Gandhi et al. 2011; see also Migliari et al. (2006) for a similar break for neutron stars] and has been indirectly constrained for Cyg X-1 based on recent mid-infrared observations (Rahoui et al. 2011). These observations confirm the theoretically predicted break location, and the recent boom in multiwavelength monitoring of many BHBs is leading to new studies of the break and its scaling with luminosity (Russell et al., in preparation). The location in frequency of the break can also be rather tightly constrained by fitting multiple spectra of BHBs and AGN in states associated with compact jets, which has the advantage of also associating a size scale for the region, and distance from the black hole. A series of works have consistently found that the region of the jet where particle acceleration initiates seems to be offset from the black hole at distances of ∼10–1000rg, depending on source luminosity, where rg is the gravitational radius, GM/c2 (Markoff, Falcke & Fender 2001; Markoff et al. 2003, 2008; Markoff, Nowak & Wilms 2005; Gallo et al. 2007; Migliari et al. 2007; Maitra et al. 2009). Accordingly, this zone should also be associated with an offset in the start of the synchrotron core from the black hole in direct imaging. Such a measurement is very difficult to do because of the spatial resolution required. At least for the jet in M87, a relatively nearby AGN (17.0 ± 0.3 Mpc; Tonry et al. 2001) with a large supermassive black hole (6.4 ± 0.5 × 109 M; Gebhardt & Thomas 2009, although note this value is twice the mass found in previous studies), the offset is ∼30–100rg (Junor, Biretta & Livio 1999; Walker et al. 2008; Hada et al. 2011). Interestingly, recent astrometry for M87 has also shown that the location of the radio core with respect to the black hole is extremely stable over several years (Asada et al. 2011). A natural explanation for such a region for the start of particle acceleration would be a shock where continual diffusive acceleration occurs (e.g. Bell 1978; Drury 1983). Because particles need to be continually accelerated, the shock should be a steady feature within the compact jets.

In self-similar, axisymmetric MHD flows, the bulk flows can be accelerated through several singular points in order: the modified slow point (MSP), the Alfvèn point (AP) and the modified fast point (MFP) (Blandford & Payne 1982, hereafter BP82). At the MFP, the velocity of the flow greatly exceeds the fast magnetosonic speed, which is the fastest that signals in the jet can travel. Thus, beyond the MFP the flow is causally disconnected from the flow closer to the central object, allowing disruptions such as shocks to form. A stable location for the MFP in a given flow would therefore provide a natural explanation for a steady feature associated with particle acceleration. This location could, however, be a function of mass and accretion power and vary from source to source within a certain range. In the models mentioned above, the shock location is a fitted free parameter that suggestively falls into the same range for a variety of sources. However, to understand if this is a physically meaningful result, we should ideally be able to show that this location corresponds to a feature that can be derived self-consistently within a theoretical framework and tied to conditions at the jet launch point. The earlier modes are based on a hydrodynamical (HD) velocity profile (Falcke & Biermann 1995; Falcke & Markoff 2000) and thus do not invoke the role of magnetic fields other than as a global parameter. However, if we move to an MHD framework, we can derive the exact location of the MFP (and assumedly the shock) based on boundary conditions defined by the launch point of the jets, and then test if the theory provides the correct scalings as observed, and as seemed to be driving the Fundamental Plane.

In this paper we take a semi-analytical approach to derive new flow solutions (a new jet model) with the explicit aim of exploring the conditions for formation of the MFP by solving the equations of relativistic MHD under the assumptions of self-similarity and axisymmetry. These assumptions greatly simplify the equations describing the jets and are consistent with the results of numerical MHD simulations, which almost always produce a magnetic field line geometry near the launch point that is remarkably self-similar and axisymmetric (see, e.g., figs 2 and 11 in McKinney 2006).

We need to make one further approximation to construct our solutions: in order to link the properties of the MFP to the conditions at the base of the jet, we need to solve for the AP and MSP in a regime which could potentially be very close to the black hole where gravity becomes important. The MSP is a natural point to associate with the launch region of the jets, similar to the sonic nozzle in HD flows, and is thus also the place where conditions could eventually be matched to the accretion inflow, such as a magnetic corona or radiatively inefficient accretion flow (RIAF; Narayan & Yi 1994; Blandford & Begelman 1999; Quataert & Gruzinov 2000; Merloni & Fabian 2002). However, gravity will play a significant role so close to the black hole; therefore, it also needs to be included in order to correctly derive and cross the MSP. Until now no semi-analytical formalism has been developed to describe a relativistic flow passing through all three singular points, because in a relativistic framework gravity is not self-similar. Therefore, in this paper, in order to accomplish this connection, we apply a bridging solution between our previous relativistic self-similar flow model without gravity (Polko, Meier & Markoff 2010, hereafter PMM10, based on the framework for self-similar flow laid out in Vlahakis & Königl 2003, hereafter VK03), which is also valid in the non-relativistic limit (again without gravity), and a non-relativistic flow model that is valid close to the black hole and that does include the effects of gravity self-consistently (Vlahakis et al. 2000, hereafter VTST00). The transition from the non-relativistic flow with gravity to our relativistic flow without gravity from PMM10 occurs sufficiently far from the black hole that gravitational effects are no longer important. The bridging of these two solutions will not describe all possible gravitational fields and velocities, but the range of solutions is wide enough to be applicable to many observed sources, which is the goal of this project. Specifically, within a cylindrical region constituting a ‘flux tube’, that can be associated with the observable ‘sheath’ of compact jets, we demonstrate that the deviations from self-similarity when gravity is present are very small. Within this scenario, we present relativistic MHD solutions that cross all three singular points, allowing us to derive the location of the shock self-consistently from boundary conditions at the base of the jet and to begin to test how well such a model fares against physical sources.

Gravity is treated in the weak (Newtonian) limit, and in this paper we have opted for a Paczyńsky–Wiita (pseudo-Newtonian) augmentation since some of our solutions approach the Schwarzschild radius, where general relativistic effects begin to be important. No other general relativistic correction terms, nor black hole spin, are treated at this time.

In Section 2, we present the details of the bridging formalism. In Section 3 we use the new formalism to obtain solutions crossing the MSP, AP and MFP, and describe their properties. We also show the effects that changing the model parameters has on the location of the MSP and MFP. In Section 4, we present an in-depth discussion of our model and results, and we draw our conclusions in Section 5.

2 METHOD

In order to tie the conditions at the start of the particle acceleration region to the conditions at the base of the jet, we need to have a full description of a hot, relativistic flow. Because close to the central object gravity is as important as temperature and the magnetic field strength, the gravitational potential must be taken into account in order to be able to describe the jet in this region. Because this region is generally subrelativistic or mildly relativistic, we start with a solution that includes gravity in a non-relativistic flow. Far from the black hole the flow will be relativistic, but gravity will no longer be important in the solution, so we can use the equations derived in PMM10 to describe this flow. We then seek solutions that satisfy both the relativistic flow without gravity and the non-relativistic flow with gravity at opposite ends of the jet.

2.1 A physical description of the flow

In the framework of self-similar relativistic MHD, we find that jets in our solutions are accelerated in the following way: starting from the base of the flow, the initial bulk acceleration is provided through thermal pressure gradients by means of sound waves by the very hot particles surrounding the central object. When the flow exceeds the sound speed, this mechanism becomes inefficient, and acceleration due to the centrifugal force of the rotating magnetic field by means of Alfvén waves takes over. After the flow surpasses the Alfvén speed, magnetic pressure becomes the dominant mechanism. When the flow velocity increases beyond the fast magnetosonic speed, the final boost is given by the pinching of the magnetic field. If at this point the flow overcollimates, a shock may form, which can accelerate the particles into the observed power-law distributions, thus constituting the start of the particle acceleration region.

In prior work (PMM10) we were unable to probe the region where the initial bulk acceleration occurs due to the absence of gravity in our formalism. However, the singular MSP can be produced only by the inward gravitational force balancing the outward magnetic and thermal pressure forces. Thus, by including gravity we can describe the jet from very close to the central object to the first instance of overcollimation, providing the connection between the conditions at these two points. We will discuss how we deal with the issues regarding self-similarity in Section 4.

2.2 A new solution technique: the $$\boldsymbol {C^\infty }$$-continuous bridging method

The matching of two flow solutions valid in different regimes of parameter or physical space is a common technique in physics. Some techniques simply solve the two different equations in the different regimes and then match the solutions at a single (and arbitrary) point. Such bridging solutions are only C0-continuous, meaning their derivatives are not continuous functions. Other methods may involve spline fitting the two solutions together in a finite bridging region. In this case, the two solutions are valid in each of their respective regimes, but in the bridging region neither is strictly valid. Instead, the spline fit is an approximation of the transition between the two solutions, and it is C1-continuous or better.

Our method for bridging the non-relativistic VTST00 MHD wind equation with gravity and the relativistic VK03 one without gravity is similar to the above approaches. That is, when gravitational forces are important and the flow is non-relativistic, we obtain the VTST00 MHD solution. However, when gravity is no longer important, the flow has the character of the VK03 relativistic solutions, whether it is relativistic or not. The difference here is that because of the nearly identical structures of the VTST00 and VK03 equations, we can perform this task with the construction of a single wind equation. The method, then, will be C-continuous, with all derivatives being continuous functions – a distinct advantage when one is dealing with equations which have singular points through which the solution must pass.

2.3 The basic $$\boldsymbol {C^\infty }$$-continuous method

We now give a brief overview of our previous work in PMM10 and describe how we here extend it by bridging it with a non-relativistic formulation including gravity valid close to the black hole.

In PMM10 we combined the Bernoulli equation (describing the forces along a field line) and the transfield equation (describing the forces perpendicular to a field line) to construct a wind equation, which describes the bulk acceleration of the flow and fully specifies the jet geometry. Using this single differential equation, it was possible to obtain jet solutions that crossed an MFP at a finite distance from the origin. The terms in the wind equation thus obtained using the formalism of VK03 turned out to be almost identical to the terms in the wind equation given in VTST00. Apart from relativistic corrections to all the corresponding terms, the only difference is the gravity term in the VTST00 wind equation, which is non-relativistic. If we modify the VK03 wind equation to include gravity, our combined wind equation will reduce exactly to the formalism of VTST00 when relativistic effects apart from gravity are negligible as is the case close to the black hole, while it also reduces exactly to the formalism of VK03 when gravity becomes negligible as is the case far away from the black hole. Between these two regimes the wind equation describes the bridging regime where relativistic effects may begin to dominate the gravitational effects. Since all these regimes are described by one continuous equation, the resulting solution will be a smooth transition between the two regimes. We prefer this approach over matching the solutions of the two regimes at an arbitrary location in the jet while trying to satisfy any number of continuity conditions.

Since the notation in VTST00 differs from VK03, we use equations (8)–(12) in VTST00 and equations (24) in VK03 to translate the notation to that used in the latter. This conversion also provides the relativistic corrections to the gravitational mass. Using this relativistic velocity in the definition of the gravity term in VTST00 yields a relativistic form of the gravity term, allowing us to include it into the VK03 wind equation. The term-by-term comparison of the two wind equations ensures that we use the proper scaling. Since some of our solutions get very close to the central black hole, we adapted the gravity term to include a pseudo-Newtonian potential (Paczyńsky & Wiita 1980).

After these steps (see Appendix A) the gravity term has the form  
\begin{equation} - \frac{\mu ^2 x_\mathrm{A}^4}{F^2 \sigma _\mathrm{M}^2} \frac{\left(1 - M^2 - x_\mathrm{A}^2\right)^2}{(1 - M^2 - x^2)^2} \left[ \frac{c^2 \varpi _\mathrm{A}G}{\mathcal {G M} \sin (\theta )} - 2 \right]^{-1}, \end{equation}
(1)
where the subscript A denotes values at the AP and μc2 is the total energy-to-mass flux ratio, xA is the cylindrical radius distance of the AP scaled to the light cylinder radius, F controls the current distribution, σM is the magnetization parameter, M is the Alfvénic Mach number, x is the cylindrical radius distance scaled to the light cylinder radius, c is the speed of light, ϖA is the cylindrical radius distance to the AP, G is the cylindrical radius distance scaled to the AP, $$\mathcal {G}$$ is the gravitational constant, $$\mathcal {M}$$ is the mass of the black hole and θ is the spherical polar angle (see Table 1 for an overview of all model parameters).

The first two fractions consist purely of constants along a given field line, providing the overall scaling. The third fraction is proportional to (γξ)2 (see equation A1), where γ is the Lorentz factor and ξc2 is the specific relativistic enthalpy (per baryon mass). The last fraction corresponds to (r/rg − 2)−1, with r the spherical radius and rg the gravitational radius, and it is this term that represents the pseudo-Newtonian potential.

Several equations besides the wind equation are needed to fully determine a solution, so we have made sure that we take into account all equations in which the gravity term appears. As explained above, gravity, as described by equation (1), shows up in the numerator of the wind equation. By evaluating the wind equation at the AP, we obtain the Alfvén regularity condition (ARC; see equation B1), which also depends on gravity. The ARC allows us to calculate the slope of M2 at the AP, pA, one of the initial parameters of the integration. The gravity term, as it appears on the right-hand side of the ARC in VK03, is given by  
\begin{equation} - \frac{x_\mathrm{A}^2}{1 - x_\mathrm{A}^2} \left[ \frac{c^2 \varpi _\mathrm{A}}{\mathcal {G M} \sin (\theta _\mathrm{A})} - 2 \right]^{-1}. \end{equation}
(2)

Another location in which the gravity term may appear is in the Bernoulli equation evaluated at the AP. This equation is used to determine μ, another initial parameter of the integration (see equation B2). However, by evaluating the Bernoulli equation of VTST00 at the AP, the gravity term in that equation vanishes there. So equation (B5) in VK03 is still valid (see Appendix C). The transfield equation is incorporated into the wind equation but not used independently in our calculations, so its dependence on gravity is irrelevant here.

By making these changes we have been able to find solutions that can pass through all three singular points of a hot, relativistic MHD flow, extending all models so far published (for a qualitative overview see Table 2).

Table 2.

Overview of self-similar, axisymmetric MHD models, classified as whether they include gravity, allow for a warm flow, allow for relativistic velocities and cross the MSP, AP and MFP, respectively.

Paper Gravity Warm Relativistic MSP AP MFP 
BP82 ✓    ✓  
Li, Chiueh & Begelman (1992  ✓  ✓  
VTST00 ✓ ✓  ✓ ✓ ✓ 
VK03  ✓ ✓  ✓  
PMM10  ✓ ✓  ✓ ✓ 
Polko et al. (this work) ✓ ✓ ✓ ✓ ✓ ✓ 
Paper Gravity Warm Relativistic MSP AP MFP 
BP82 ✓    ✓  
Li, Chiueh & Begelman (1992  ✓  ✓  
VTST00 ✓ ✓  ✓ ✓ ✓ 
VK03  ✓ ✓  ✓  
PMM10  ✓ ✓  ✓ ✓ 
Polko et al. (this work) ✓ ✓ ✓ ✓ ✓ ✓ 
Table 2.

Overview of self-similar, axisymmetric MHD models, classified as whether they include gravity, allow for a warm flow, allow for relativistic velocities and cross the MSP, AP and MFP, respectively.

Paper Gravity Warm Relativistic MSP AP MFP 
BP82 ✓    ✓  
Li, Chiueh & Begelman (1992  ✓  ✓  
VTST00 ✓ ✓  ✓ ✓ ✓ 
VK03  ✓ ✓  ✓  
PMM10  ✓ ✓  ✓ ✓ 
Polko et al. (this work) ✓ ✓ ✓ ✓ ✓ ✓ 
Paper Gravity Warm Relativistic MSP AP MFP 
BP82 ✓    ✓  
Li, Chiueh & Begelman (1992  ✓  ✓  
VTST00 ✓ ✓  ✓ ✓ ✓ 
VK03  ✓ ✓  ✓  
PMM10  ✓ ✓  ✓ ✓ 
Polko et al. (this work) ✓ ✓ ✓ ✓ ✓ ✓ 

2.4 The effects of including gravity on the solutions

In every instance of gravity the ratio $$\varpi _\mathrm{A}/ \mathcal {M}$$ appears. Thus, for a central object that is twice as massive, if the AP is moved twice as far from the axis of symmetry, the solution remains unchanged. This result is a direct consequence of mass scaling, since all properties are expressed in gravitational radii. Regulating the gravitational force exerted by the compact object can therefore be achieved by changing either $$\mathcal {M}$$ or ϖA and, apart from the overall physical size of the system, these two actions are equivalent. Increasing the effect of gravity can thus be thought of either as increasing the mass of the central object (increasing the reach of the gravitational well) or as decreasing the radius of the AP (moving the AP deeper into the gravitational well).

One way to show these effects would be by starting with a singular solution containing both an MFP and an MSP, slowly increasing ϖA to decrease the effect of gravity, while keeping solutions with an MFP and MSP. Unfortunately, these solutions do not smoothly transition into solutions that only have an MFP. It is therefore impossible to directly compare these new solutions crossing all three singular points to solutions without an MSP and we are left noting the effects within these new solutions.

The significance of having a solution with an MSP is that only then it is possible to connect conditions at the MFP to conditions very close to the central object. Even though the extension of a field line in linear distance may not be much, we can be sure that the results in the sub-Alfvénic regime have meaning by satisfying a strict boundary condition. We equate the MSP with the launch point of the jet and fit the conditions at the MSP to the accretion flow. Hence, extending the solutions to the MSP allows a link between the accretion flow and the start of the particle acceleration region, which we will use in later work.

Compared with the Paczyńsky–Wiita potential, the Newtonian potential has a weaker gravitational force at the same radius. As it is the gravitational force that balances all other forces at the MSP, using a Newtonian potential would move the MSP closer to the black hole, possibly beyond the event horizon for certain solutions. To avoid this, we adopt the Paczyńsky–Wiita potential to approximate the gravitational potential close to the black hole. The adaptation from a Newtonian potential is straightforward (see Appendix A).

2.5 Model parameters

Table 1.

List of model parameters.

Fundamental parameters 
F Sets the radial dependence of the magnetic field strength 
Γ Adiabatic index 
xA Cylindrical radius in terms of the light cylinder radius 
σM Magnetization parameter 
q Dimensionless adiabatic coefficient 
ϖA Cylindrical radius 
$$\mathcal {M}$$ Mass of the central object 
Derived parameters 
θA Angle with respect to the axis of symmetry 
ψA Angle of the field line with respect to the disc 
pA Derivative of the Alfvénic Mach number squared with respect to θ 
σA Magnetization function (Poynting-to-matter energy flux ratio) 
μc2 Total energy-to-mass flux ratio 
ξc2 Specific (per baryon mass) relativistic enthalpy 
Fundamental parameters 
F Sets the radial dependence of the magnetic field strength 
Γ Adiabatic index 
xA Cylindrical radius in terms of the light cylinder radius 
σM Magnetization parameter 
q Dimensionless adiabatic coefficient 
ϖA Cylindrical radius 
$$\mathcal {M}$$ Mass of the central object 
Derived parameters 
θA Angle with respect to the axis of symmetry 
ψA Angle of the field line with respect to the disc 
pA Derivative of the Alfvénic Mach number squared with respect to θ 
σA Magnetization function (Poynting-to-matter energy flux ratio) 
μc2 Total energy-to-mass flux ratio 
ξc2 Specific (per baryon mass) relativistic enthalpy 

Note. The subscript A means the value of the variable at the AP. For a complete description of the parameters, please see PMM10.

Table 1.

List of model parameters.

Fundamental parameters 
F Sets the radial dependence of the magnetic field strength 
Γ Adiabatic index 
xA Cylindrical radius in terms of the light cylinder radius 
σM Magnetization parameter 
q Dimensionless adiabatic coefficient 
ϖA Cylindrical radius 
$$\mathcal {M}$$ Mass of the central object 
Derived parameters 
θA Angle with respect to the axis of symmetry 
ψA Angle of the field line with respect to the disc 
pA Derivative of the Alfvénic Mach number squared with respect to θ 
σA Magnetization function (Poynting-to-matter energy flux ratio) 
μc2 Total energy-to-mass flux ratio 
ξc2 Specific (per baryon mass) relativistic enthalpy 
Fundamental parameters 
F Sets the radial dependence of the magnetic field strength 
Γ Adiabatic index 
xA Cylindrical radius in terms of the light cylinder radius 
σM Magnetization parameter 
q Dimensionless adiabatic coefficient 
ϖA Cylindrical radius 
$$\mathcal {M}$$ Mass of the central object 
Derived parameters 
θA Angle with respect to the axis of symmetry 
ψA Angle of the field line with respect to the disc 
pA Derivative of the Alfvénic Mach number squared with respect to θ 
σA Magnetization function (Poynting-to-matter energy flux ratio) 
μc2 Total energy-to-mass flux ratio 
ξc2 Specific (per baryon mass) relativistic enthalpy 

Note. The subscript A means the value of the variable at the AP. For a complete description of the parameters, please see PMM10.

By including gravity in the equations, two additional parameters need to be specified before an integration. These are the mass of the central object, $$\mathcal {M}$$, and the cylindrical radius of the AP, ϖA. As mentioned above, only the ratio of these parameters appears in the equations, so these two parameters are linearly dependent and solutions with the same ratio are indistinguishable if all other parameters are fixed. Since a solution is given in gravitational radii, this means that any solution can be scaled up or down to any physical size. For our initial solution we have chosen to set $$\mathcal {M}$$ to the mass of a typical stellar mass black hole of 10 M (although this value will eventually be set to the deduced mass of an observed black hole system) and to set ϖA to 5rg. We will look at the effects of mass scaling at a later time.

The parameters F and Γ are set to a single value for now, but can be changed if necessary. For example, since F influences the geometry of the jet, it also affects the collimation. By allowing it to vary, we will have the freedom to adjust the degree of collimation to fit the broad-band data on our sources. Due to the self-similar nature of the equations, the height of the MFP will depend on the initial radius at the disc. Because we want a representative location, we choose to anchor the field line in the most active part of the accretion disc close to the central object. Within the innermost stable circular orbit (ISCO) the accreting matter can no longer form a disc and will move towards the central object as a radiatively inefficient flow. Therefore, we set F = 0.75, which corresponds to a wide range of accretion disc models, including BP82, Shakura–Sunyaev (Shakura & Sunyaev 1973) and RIAF models. We also set Γ = 5/3, appropriate for non-radiation-dominated jets in the hard state, where the radiation and the particles do not behave as a single fluid.

The choice of parameters to fit the singular points will affect the solution we find. Were we to use θA and ψA to fit for the MFP and MSP positions, the angle of the AP and the slope of the field line there would be varied to find a singular solution. If, on the other hand, x2A and q were used to fit for the MFP and MSP, the light cylinder radius (and hence angular velocity) and dimensionless adiabatic coefficient would be varied instead. This last combination most closely approaches self-similarity out of all the pairs of parameters that we have studied. Therefore, in this work we will use x2A and q as fitting parameters. It is important to note that by choosing different fitting parameters only the way we move through parameter space changes, and not the set of solutions themselves.

3 RESULTS

3.1 First solution

Table 3.

Parameters of solutions.

   θA ψA  ϖA $$\mathcal {M}$$      
 F Γ (deg) (deg) σM (rgM x2A q pA σA μ 
PMM10a 0.75 5/3 50 55 2.539 81 − − 0.75 0.12 −1.666 51 2.603 90 9.851 17 
This work 0.75 5/3 50 60 2.5 10 0.755 500 0.271 221 −1.277 54 2.257 11 12.4036 
   θA ψA  ϖA $$\mathcal {M}$$      
 F Γ (deg) (deg) σM (rgM x2A q pA σA μ 
PMM10a 0.75 5/3 50 55 2.539 81 − − 0.75 0.12 −1.666 51 2.603 90 9.851 17 
This work 0.75 5/3 50 60 2.5 10 0.755 500 0.271 221 −1.277 54 2.257 11 12.4036 

aThe parameters given are for the solution c in PMM10. This is the solution we will compare our new result with. The values for the first seven parameters (F through $$\mathcal {M}$$) of the solution in this work are exact, and for the last five (x2A through μ) they are rounded off. Because singular solutions require high precision, the rounded-off numbers are given with six significant digits.

Table 3.

Parameters of solutions.

   θA ψA  ϖA $$\mathcal {M}$$      
 F Γ (deg) (deg) σM (rgM x2A q pA σA μ 
PMM10a 0.75 5/3 50 55 2.539 81 − − 0.75 0.12 −1.666 51 2.603 90 9.851 17 
This work 0.75 5/3 50 60 2.5 10 0.755 500 0.271 221 −1.277 54 2.257 11 12.4036 
   θA ψA  ϖA $$\mathcal {M}$$      
 F Γ (deg) (deg) σM (rgM x2A q pA σA μ 
PMM10a 0.75 5/3 50 55 2.539 81 − − 0.75 0.12 −1.666 51 2.603 90 9.851 17 
This work 0.75 5/3 50 60 2.5 10 0.755 500 0.271 221 −1.277 54 2.257 11 12.4036 

aThe parameters given are for the solution c in PMM10. This is the solution we will compare our new result with. The values for the first seven parameters (F through $$\mathcal {M}$$) of the solution in this work are exact, and for the last five (x2A through μ) they are rounded off. Because singular solutions require high precision, the rounded-off numbers are given with six significant digits.

Starting with the warm (initial ξ ≈ 4.7) and fast (final γ ≈ 10) solution c in PMM10, we use our new method to explore parameter space to determine a solution with an MSP. The parameter values of this new solution are very close to those of solution c (see Table 3), establishing the relative ease with which solutions can be found. Fig. 1 shows the numerator and denominator of the wind equation, and their ratio corresponding to the derivative of the square of the Alfvénic Mach number with respect to the poloidal spherical angle θ. From this ratio, it is clear that we have indeed found a smooth crossing. In comparison with the solution c, the AP is still located at θ = 0.87 rad or 50° from the axis of symmetry, but the MFP has moved outwards to θ = 0.029 rad or 1$ ${.\!\!\!\!\!\!^{\prime}}$$65, while the MSP (only ostensible in the solution c) has moved inwards to θ = 0.92 rad or 52$ ${.\!\!\!\!\!\!^{\prime}}$$6. This angle corresponds to a spherical radius of 5.96rg, which is just within the ISCO.

Figure 1.

The reference solution showing both an MFP and an MSP. The dotted/red and dashed/blue lines show the numerator and denominator of the wind equation, respectively, while the solid/black line shows their ratio. The latter determines the total bulk acceleration of M2 with respect to polar angle θ. The vertical line shows the location of the AP. The vertical axis is a ‘scaled logarithm’ of the plotted parameters, i.e. y = sign(x)log 10[1 + abs(x)/10− 12] to clearly show the variables over many orders of magnitude. At low θ, a small change in angle corresponds to a vast change in height, also contributing to the near-vertical change of sign of the bulk acceleration.

Figure 1.

The reference solution showing both an MFP and an MSP. The dotted/red and dashed/blue lines show the numerator and denominator of the wind equation, respectively, while the solid/black line shows their ratio. The latter determines the total bulk acceleration of M2 with respect to polar angle θ. The vertical line shows the location of the AP. The vertical axis is a ‘scaled logarithm’ of the plotted parameters, i.e. y = sign(x)log 10[1 + abs(x)/10− 12] to clearly show the variables over many orders of magnitude. At low θ, a small change in angle corresponds to a vast change in height, also contributing to the near-vertical change of sign of the bulk acceleration.

The values of the velocity, the Lorentz factor, the magnetic and electric field strength, the density and the pressure along a reference field line are given in Fig. 2. The poloidal velocity monotonically increases from 0.067c, through 0.11c at the MSP, to very close to the speed of light at the MFP. The toroidal velocity starts off positive, turns negative after the AP and returns to positive before the MFP. The Lorentz factor slowly increases from 1.17 at the MSP to a final value of 12.3.

Figure 2.

Various physical quantities that result from integrating the wind equation depicted in Fig. 1, plotted against polar angle θ. The MFP, AP and MSP are marked. Panel (a) shows the Lorentz factor (γ), the poloidal velocity (Vp) and toroidal velocity (Vϕ) of our canonical solution. Please note that the toroidal velocity starts out positive, turns negative just beyond the AP (the dotted line) and then becomes positive again later in the outflow. Panel (b) shows the electric field energy (E2/8π), the toroidal magnetic field energy (B2ϕ/8π), the poloidal magnetic field energy (B2p/8π), the density (ρ0c2) and the pressure (P). These have all been divided by the scaling factor B20αF − 2, where B0 is a reference magnetic field strength and α is the square of the cylindrical distance of the AP in terms of a reference length (α = ϖ2A02), as defined in VK03. Although the square is plotted, Bϕ is negative everywhere in this model.

Figure 2.

Various physical quantities that result from integrating the wind equation depicted in Fig. 1, plotted against polar angle θ. The MFP, AP and MSP are marked. Panel (a) shows the Lorentz factor (γ), the poloidal velocity (Vp) and toroidal velocity (Vϕ) of our canonical solution. Please note that the toroidal velocity starts out positive, turns negative just beyond the AP (the dotted line) and then becomes positive again later in the outflow. Panel (b) shows the electric field energy (E2/8π), the toroidal magnetic field energy (B2ϕ/8π), the poloidal magnetic field energy (B2p/8π), the density (ρ0c2) and the pressure (P). These have all been divided by the scaling factor B20αF − 2, where B0 is a reference magnetic field strength and α is the square of the cylindrical distance of the AP in terms of a reference length (α = ϖ2A02), as defined in VK03. Although the square is plotted, Bϕ is negative everywhere in this model.

The poloidal and (negative) toroidal magnetic field strengths first increase in magnitude and then decrease, both peaking around the AP. The electric field has the same behaviour. The density and pressure both drop monotonically as is expected for an expanding jet. Just beyond the MFP the jet overcollimates, causing both density and pressure to increase again.

Fig. 3 shows the geometry, the energetics and the Alfvénic Mach number of the jet. The reference field line, plotted in the upper-left panel, shows that for most of the expansion the jet is parabolic with the height of the jet being approximately the radius to the power of 4/3. At its highest point beyond the MFP, the jet has overcollimated to only 2 per cent of the maximum width, attained just before the MFP. The main difference between this solution and the solution c from PMM10 is that the height-to-width ratio has increased by a factor of 2.5, which means it has comparatively a narrower opening angle.

Figure 3.

Various physical quantities that result from integrating the wind equation depicted in Fig. 1, now plotted against cylindrical radius ϖ instead of θ. These plots are similar to those in fig. 4 of PMM10. Panel (a) shows the geometry of the field line. Panel (b) shows the magnetic energy (S ≡ −ϖΩBϕAc2), the thermal energy including the rest mass (ξ) and the kinetic energy ([γ − 1]ξ). Panel (c) shows the opening half-angle of the outflow (π/2 − ψ) and the causal connection opening angle ($$\arcsin [1 / \gamma ]$$). Panel (d) shows the cylindrical radius in units of the ‘light cylinder’ radius (x) and the Alfvénic Mach number (M). The MSP, AP and MFP (from left to right) are indicated for each quantity.

Figure 3.

Various physical quantities that result from integrating the wind equation depicted in Fig. 1, now plotted against cylindrical radius ϖ instead of θ. These plots are similar to those in fig. 4 of PMM10. Panel (a) shows the geometry of the field line. Panel (b) shows the magnetic energy (S ≡ −ϖΩBϕAc2), the thermal energy including the rest mass (ξ) and the kinetic energy ([γ − 1]ξ). Panel (c) shows the opening half-angle of the outflow (π/2 − ψ) and the causal connection opening angle ($$\arcsin [1 / \gamma ]$$). Panel (d) shows the cylindrical radius in units of the ‘light cylinder’ radius (x) and the Alfvénic Mach number (M). The MSP, AP and MFP (from left to right) are indicated for each quantity.

The upper-right panel of Fig. 3 shows the partition of energy of the jet. After a small drop, the kinetic energy first increases at the expense of the thermal energy, signifying initial bulk acceleration due to a gas pressure gradient. After the flow has cooled, magnetic acceleration takes over and continues to accelerate the flow also after overcollimation. Only just before the flow hits the axis does the Lorentz factor decrease again with a corresponding increase in the magnetic field strength.

The lower-left panel of Fig. 3 shows the opening half-angle (π/2 − ψ) of the flow (also compare with Fig. 5) and the causal connection opening angle. Both angles are measured from the axis of symmetry. The opening half-angle shows that after a relatively cylindrical start, the flow widens before slowly collimating again. The causal connection opening angle is the equivalent of a sonic Mach cone. The flow cannot influence anything outside of the forward pointing cone with this angle. As the velocity increases, this cone becomes more narrow. This plot is almost identical to that for the solution c.

The lower-right panel of Fig. 3 shows the Alfvénic Mach number, the velocity in units of the Alfvén speed. After a long initial rise, shortly after overcollimation it decreases again. This decrease relates to the leftmost part of Fig. 1 where the numerator crosses zero again, without a corresponding crossing of the denominator, causing a decrease of the Mach number. Because the Mach number is defined in terms of the Alfvén speed, it is not the overall velocity that is decreasing, as there is no corresponding decrease in the Lorentz factor. It is rather the magnetic field strength that increases due to the overcollimation, leading to a higher Alfvén speed. The plot of the Mach number is very similar to that for the solution c, but there is a 20-fold increase of M. Note that in PMM10 we plotted the square of the Mach number.

3.2 An initial exploration of parameter space

In this section we describe our exploration of parameter space. As a starting point we pick our first solution. We change the value of a single parameter (θA, ψA, σM or $$\mathcal {M}$$) until it is no longer possible to obtain a singular solution by fitting x2A, q and pA. The effect of changing these four parameters on the height of the MFP, AP and MSP is shown in Fig. 4.

Figure 4.

Examples of how the heights (in gravitational radii) of the solutions’ three singular points change as each of the four principal free parameters is changed – solid/red line: MFP; dashed/blue line: AP; dotted/black line: MSP. The initial solution (see Table 3) is indicated by the vertical black line. The parameters varied are as follows: panel (a) – poloidal spherical angle of the AP (θA); panel (b) – the angle the field line makes with the disc at the AP (ψA); panel (c) – magnetization (σM); and panel (d) – black hole mass ($$\mathcal {M}$$) in units of M. The horizontal axis in each of these plots is linear, not logarithmic. Note the dramatic excursions in the height of the MFP (three orders of magnitude) with only modest changes in the free parameters. Note also that the MSP and AP heights are still outside the black hole horizon (see Fig. 5), i.e. even though the height of the MSP is <2rg, the MSP spherical radius for these solutions remains outside the black hole horizon (r > 2rg).

Figure 4.

Examples of how the heights (in gravitational radii) of the solutions’ three singular points change as each of the four principal free parameters is changed – solid/red line: MFP; dashed/blue line: AP; dotted/black line: MSP. The initial solution (see Table 3) is indicated by the vertical black line. The parameters varied are as follows: panel (a) – poloidal spherical angle of the AP (θA); panel (b) – the angle the field line makes with the disc at the AP (ψA); panel (c) – magnetization (σM); and panel (d) – black hole mass ($$\mathcal {M}$$) in units of M. The horizontal axis in each of these plots is linear, not logarithmic. Note the dramatic excursions in the height of the MFP (three orders of magnitude) with only modest changes in the free parameters. Note also that the MSP and AP heights are still outside the black hole horizon (see Fig. 5), i.e. even though the height of the MSP is <2rg, the MSP spherical radius for these solutions remains outside the black hole horizon (r > 2rg).

This plot also shows a physical reason why the solutions do not fill up the whole of parameter space. For high values of ψA and for low values of σM and $$\mathcal {M}$$, the MSP approaches the AP. Beyond the location where they coincide, no further solution is possible.

If θA and ψA are changed together so that their sum is roughly constant, the range of solutions spans approximately 45°. If they are changed separately, the range is significantly smaller. σM has the largest spread in this case ranging from 0.4 to above 60, monotonically increasing the height of the MFP. $$\mathcal {M}$$ has the range 4–12 M, spanning most of the mass range of astrophysical black holes for an AP distance of ϖA ≈ 74 km.

Fig. 4 gives the bounds for our solutions. Since we connect the start of the particle acceleration region with the MFP, we are mainly interested in its location within the jet. In the solutions found so far, the height of the MFP ranges from 104rg to 108rg. Because this work is by no means an exhaustive exploration of parameter space, this range provides only a preliminary indication. Nevertheless, it already allows us to model a variety of sources. The best way to change the height of the MFP seems to be adjusting θA. The height monotonically decreases while increasing this parameter. Another possibility is decreasing σM. ψA and $$\mathcal {M}$$ seem to have little effect over a wide range of values. The MSP should lie in or near the corona; thus, locations too far away and too close to the black hole should be avoided. The height of the MSP ranges from 2rg to 10rg. This lower value is very close to the Schwarzschild radius and justifies the inclusion of the pseudo-Newtonian potential. The parameter $$\mathcal {M}$$ has the largest effect on the location of the MSP.

The effect of changing $$\mathcal {M}$$ on the geometry of the field lines is shown in Fig. 5. This plot shows only the single reference field line of a solution. It is clear that the height of the MSP can be smaller than the Schwarzschild radius, but the spherical radius cannot. This geometry explains why the height of the MSP can be smaller than 2rg in Fig. 4. As shown by the lower-right panel of this figure, the solutions with a higher black hole mass collimate first and reach a lower total height (in rg).

Figure 5.

The effect of changing $$\mathcal {M}$$ on the field line geometry. The solid black/line is our canonical solution in Table 3. $$\mathcal {M}$$ ranges from 7 to 10 M in steps of 0.5 M. In panel (a), the MSP and AP are indicated by the black dots. The arrows point to the AP. In panel (b) the MFP is indicated by black dots. Note that the height of the MSP (indicated at the lower end of the line) can be smaller than the Schwarzschild radius, but its spherical radius always remains outside the horizon.

Figure 5.

The effect of changing $$\mathcal {M}$$ on the field line geometry. The solid black/line is our canonical solution in Table 3. $$\mathcal {M}$$ ranges from 7 to 10 M in steps of 0.5 M. In panel (a), the MSP and AP are indicated by the black dots. The arrows point to the AP. In panel (b) the MFP is indicated by black dots. Note that the height of the MSP (indicated at the lower end of the line) can be smaller than the Schwarzschild radius, but its spherical radius always remains outside the horizon.

3.3 Self-similarity

By including gravity in the relativistic equations of MHD, the assumption of self-similarity along the field line is broken. Near the black hole, the flow is non-relativistic and the wind equation reduces to the form of VTST00, which is self-similar. Far away from the centre, the gravitational energy is negligible and the wind equation assumes the form of VK03, which is also self-similar. This combination means that along the integrated field line, the flow transitions from one form of self-similarity to another. It should be noted that the reference field line we integrate will always be continuous.

Since one field line does not constitute a jet, we would like to be able to describe the flow in a region around our reference field line. We move from one field line to another by varying the parameter α = ϖ2A02, where ϖA is the cylindrical radius of the AP for a specific field line and ϖ0 is a reference length. By choosing this reference length to be the cylindrical radius of the AP for our reference field line, without loss of generality, that field line has α = 1.

Changing α to select different field lines then becomes equal to changing ϖA or $$\mathcal {M}$$. Since this parameter is one of our model parameters, the other model parameters have to be fitted in order to obtain solutions passing through all three singular points. Although the best result would be obtained by changing all parameters simultaneously, for simplicity, we limit ourselves to two. We are free to choose any two parameters, and when using x2A and q the field lines do not cross, which is required for self-similarity. Because θA and ψA are therefore constant for all field lines in a region around our reference field line, at the AP self-similarity is exact. We will define this region, where the field lines do not cross and satisfy self-similarity to within a specified error, as a flux tube.

To show that the non-crossing of field lines holds for a large part of parameter space, we have performed this test at different parameter values. Fig. 6 shows the deviations from self-similarity of this solution by dividing several field lines around a reference field line by this central field line. For perfect self-similarity, this ratio should be a constant. In our case, self-similarity is maintained for the majority of the jet. Only very close to the MSP do the deviations become more pronounced due to gravity. These deviations are also demonstrated by the MSP occurring at slightly different angles for different values of the black hole mass $$\mathcal {M}$$ (see Fig. 5). If we allow a maximum deviation of 10 per cent, which in this case will occur at the MSP since we will terminate our solutions there, the width of a flux tube would be about 0.30 of the radius of the central field line, which is a significant fraction.

Figure 6.

The ratios of the radial size of several field lines with regularly increasing ϖA fitted with x2A and q. The bottom line has α = 0.64 and the top line has α = 1.6 (where α ≡ ϖ2A02). The horizontal dashed lines show the values for exact self-similarity. The deviations from self-similarity only get pronounced near the MSP. The shaded region shows the width of the flux tube for which the combined deviation is smaller than 10 per cent beyond the MSP. The parameters of the reference solution are x2A = 0.01, σM = 0.01, q = 0.01, ϖA = 5.08rg, $$\mathcal {M} = 10\,\mathrm{M}_{\odot }$$, θA = 0.952 915 and ψA = 89.1119. The Lorentz factor at the MFP is just above 11, showing that this solution is relativistic.

Figure 6.

The ratios of the radial size of several field lines with regularly increasing ϖA fitted with x2A and q. The bottom line has α = 0.64 and the top line has α = 1.6 (where α ≡ ϖ2A02). The horizontal dashed lines show the values for exact self-similarity. The deviations from self-similarity only get pronounced near the MSP. The shaded region shows the width of the flux tube for which the combined deviation is smaller than 10 per cent beyond the MSP. The parameters of the reference solution are x2A = 0.01, σM = 0.01, q = 0.01, ϖA = 5.08rg, $$\mathcal {M} = 10\,\mathrm{M}_{\odot }$$, θA = 0.952 915 and ψA = 89.1119. The Lorentz factor at the MFP is just above 11, showing that this solution is relativistic.

Beyond the AP, deviations can be caused by changes in the fitting parameters to obtain singular solutions for the different field lines. The bigger the differences in these parameters, the greater the deviations and the narrower the resultant flux tube. The required changes in the parameters depend on the location in parameter space, and so the allowed width depends on the parameter values chosen. For the parameters of Fig. 6, the deviations beyond the AP are very small. While our equations do not strictly satisfy self-similarity, we have shown that it holds very well and that determining the width of the flux tube where it does is straightforward.

4 DISCUSSION

By including gravity in the MHD equations, we are able to extend earlier solutions that crossed the MFP and AP to also go through the MSP, allowing us to describe a jet that crosses all three singular points. Because of this lower boundary condition, we now have a reliable description of the jet below the AP. This description gives a smooth solution from very near the central object out to the point of overcollimation, using a single partial differential equation to describe all regimes. This approach produces a workable physical model, allowing us to tie the jet properties to the conditions at the base and providing a self-consistent determination of the start of the particle acceleration region.

At the same time, however, the addition of gravity also violates the conditions for self-similarity. The reason why a single self-similar relativistic flow equation, with gravity, has not been derived is that relativistic flow has one scaling with radius, while gravity has another. Our C-continuous bridging method has not changed that situation: while describable with a single, continuous equation, our solutions in the relativistic part of the flow without gravity will have different dependences on the radius parameter α than in the non-relativistic part near the black hole that includes gravity. That is, at least one term in our wind equation will have a dependence on the radius parameter α.

In effect, our modified wind equation creates two regions with self-similar geometry, but with different self-similar dependences in each region. Because our focus is on the bulk acceleration and collimation of jets in relativistic sources, our approach to this issue is to choose the self-similarity of the relativistic VK03 equations and restrict the different dependence on α to the gravity term only. That is, our solutions will not be strictly self-similar in the low-speed part of the flow with gravity (the VTST00 regime), but they will be self-similar far from the black hole in the VK03 regime (see Fig. 6). This highlights another advantage of this approach: since the gravity term is an algebraic one only (not involving any derivatives of the flow parameters with respect to either radius or polar angle), all the physical, radial and angular dependences of the original equations will be preserved.

The different self-similarities within a single solution are not unique to our approach; it will be true of any method that attempts to bridge self-similar solutions with and without gravity. In this regard, interpretation of the solutions will be an important aspect of this study. There are two main issues to consider.

First, while accretion discs and jets may have a self-similar character, both simulations and observations show that their activity is concentrated in specific regions. For example, much of the accretion, and therefore jet, power is generated near the disc inner edge. And features like the one we seek (a strong collimation shock in the flow) will occur primarily at a specific point in the jet and be observed as such, rather than being spread over a large range in radius. Therefore, we shall concentrate on only one bundle (‘flux tube’) of magnetic field and stream lines, and assume that this flux tube is anchored near the disc inner edge and that it passes through the most important features of the jet. This will limit the range in α in which we are interested.

Secondly, in order that the solutions remain reasonably self-similar-like, we shall choose solutions in which the field lines do not cross within that flux tube. By choosing a specific combination of the free parameters to fit a solution with three singular points, the field lines around our reference field line behave well. We have found that selecting x2A and q as fitting parameters for the singular points produces results that satisfy self-similarity well within a flux tube of finite width.

With our method, therefore, we will restrict our solutions to a limited, and conservative, region of parameter space. The flow in our solutions must remain reasonably non-relativistic in the VTST00 part of the flow (when gravity is important), and we must consider only flux tubes of narrow enough width that field lines satisfy self-similarity to within a specified error in this region with gravity.

Fig. 7 can be used to see whether the assumptions of the two regimes are not violated for a particular solution. The kinetic energy (γ − 1) should be negligible near the black hole, while the gravitational potential energy $$\left( \frac{\mathcal {G M}}{c^2 r} \right)$$ should be significant, satisfying the conditions of the non-relativistic regime and demonstrating the importance of gravity. Around the AP the kinetic energy overtakes, without becoming relativistic. This is the bridging region where gravity becomes unimportant. Only beyond the AP does the flow become truly relativistic, while gravity becomes negligible, satisfying the conditions of the relativistic regime.

Figure 7.

Comparison of the gravitational potential energy and the kinetic energy along the jet of the bottom solution detailed in Table 3. With this plot we can select solutions where the gravitational energy $$\mathcal {G M}/\left( c^2 r \right)$$ is important close to the black hole, whereas the kinetic energy (γ − 1) only becomes relativistic further along the jet.

Figure 7.

Comparison of the gravitational potential energy and the kinetic energy along the jet of the bottom solution detailed in Table 3. With this plot we can select solutions where the gravitational energy $$\mathcal {G M}/\left( c^2 r \right)$$ is important close to the black hole, whereas the kinetic energy (γ − 1) only becomes relativistic further along the jet.

In this work we have only made a cursory exploration of the full parameter space, and it seems likely that part of that space will allow for solutions with smaller MFP values. The location of the MSP falls within the range of an RIAF disc model. Also, the velocity at the MSP, around c/3, is very reasonable. These features make it plausible that we will be able to match the conditions to an accretion flow as well as to the size of the jet base as indicated by fits to the spectrum (Markoff et al. 2005), allowing us to determine the location of the shock region from the conditions at the base of the jet.

It is important to note that since our assumption of non-relativistic flow (when gravity is important) eliminates the part of parameter space that produces relativistic flows close to the black hole, we exclude solutions that may relate to physical sources. Future work would include a more general model that allows relativistic flow in the gravitational field. However, the model in this paper provides a good description of the jet structure near the MSP when the flow is non-relativistic.

We should also note that the gravitational field we chose to include here is, at best, a pseudo-Newtonian (or pseudo-Schwarzschild) one. Our solutions, therefore, do not take into account any metric rotation that would occur near a Kerr black hole, for example. However, our method might be able to hint at the presence of a Kerr hole, if, for example, after fitting to broad-band data, we find that the footpoint of the field line may be significantly less than the ISCO of a Schwarzschild black hole (≪6 rg) or that the required magnetization is significantly greater than would be typical at the ISCO of a normal Keplerian disc.

Another issue neglected by our computations (and, indeed, by all steady-state analyses) is the effect of instabilities on the accelerating jet flow. We will briefly discuss fluid and MHD instabilities in the weak-field limit and then MHD instabilities in the strong-field limit as well. Fig. 2(b) compares the relative strengths of fluid dynamical forces (i.e. pressure P) with electromagnetic forces (e.g., (B2 + E2)/8π). We see that the former dominate only below the AP, in the slow magnetosonic region, so it is only there where we expect Kelvin–Helmholtz or magnetorotational instabilities possibly to be important. Indeed, the MSP itself may lie in the atmosphere of a turbulent accretion disc, giving rise to a somewhat more complex structure in the sub-Alfvénic region than that assumed here. Nevertheless, recent two- and three-dimensional simulations of relativistic jets launched from turbulent accretion discs (McKinney 2006; McKinney & Blandford 2009) show that MHD jets similar to those investigated here (including a well-defined classical slow point) do indeed form and are not disrupted by weak-field instabilities near the accretion disc. Above this point, like other outflowing magnetospheres (solar, pulsar, etc.), the dynamical forces in the flow are so dominated by electromagnetic forces (eventually by many orders of magnitude) that any weak-field instabilities will be strongly suppressed or at least unimportant in affecting the kinematics of the accelerating jet.

However, strong-field instabilities [in particular, the current-driven helical kink (CDI)] need to be looked at more closely. Nakamura & Meier (2004), for example, showed that in non-relativistic jet flow, rotation velocities well above the Alfvén speed can suppress the helical kink, as can a steep external pressure gradient. They speculated that relativistic flow would further reduce the development of kinks, but a similar detailed relativistic study has not yet been performed. Again, like one of the few relativistic three-dimensional MHD simulations, McKinney & Blandford (2009) can give some insight into the behaviour of real MHD jets in the strong-field region. While the accelerating jet appears to be affected somewhat by helical kinking in that study, it does not appear to be disrupted by the CDI. It remains a viable jet well beyond the classical fast point, where the rotational speed should be significantly super-Alfvénic. However, no three-dimensional simulation so far has followed the acceleration out to the MFP region, let alone for long model times to achieve a quasi-steady state. And detailed numerical parameter studies are even further in the future. So, the question of whether MHD jet acceleration can be largely stable to strong-field MHD instabilities is still an open question.

5 CONCLUSIONS

We have shown that it is possible to extend the time-independent, semi-analytic solutions of PMM10 to solutions that include gravity using a C-continuous bridging method, connecting two valid regimes with a smooth transition region. We also extend the MHD wind/jet equation to include a pseudo-Newtonian potential for gravity. For the first time these solutions cross all three singular points and describe a relativistic jet from the vicinity of a central black hole to the point where the flow hits the axis of symmetry. Although these solutions do not satisfy strict self-similarity, the deviations within a flux tube of a certain width can be estimated and controlled, and thus a physically relevant model can be constructed.

We have discovered a solution crossing all three singular points with parameters very close to the one from our previous paper. This new solution shows the relative ease with which solutions can be found and allows for a comparison. The main differences are the higher distance of the MFP and a significantly higher velocity, both owing to a higher bulk acceleration at the AP. Otherwise, the two solutions look very similar, showing that gravity has little influence beyond the region close to the central object. Allowing the creation of an MSP, including gravity, enables us to have a reliable description of the flow below the AP.

We have made a cursory examination of parameter space by changing one parameter of a single solution at a time. This exploration gives a preliminary indication of the extent of the solution space, while at the same time showing the effects it has on the solutions. In particular, we have found that the poloidal spherical angle of the AP has the most significant effect on the height of the MFP. On the other hand, the ratio of the AP to the light cylinder radius has the greatest effect on the height of the MSP.

The multidimensional solution space is constrained by physical and mathematical conditions. Changing only one parameter at a time, when the light cylinder radius approaches the AP, the AP moves outwards, the temperature is increased, or the magnetization decreased, and the jets become wider and shorter, with the MSP approaching the AP. By increasing the cylindrical radius of the AP, the effects of gravity decrease. Eventually this decrease leads to the regime where an MSP cannot be created anymore, approaching the situation without gravity. Despite all these constraints, the solution space allows a wide range of parameter values to be chosen, translating into a wide range of properties, like the location of the MSP and MFP, the velocity, magnetic energy, density and pressure of the flow.

By matching conditions at the MSP to an accretion flow model, the mathematical parameters considered free in this work will be tied to the conditions at the base of the jet and subsumed into the model, with the wide range of properties in the solution space assuring a good fit. Consequently, the height of the MFP, and the start of the particle acceleration region, will be uniquely determined by the conditions close to the central object, providing a self-consistent connection. After integration into a model that can determine the spectral emission of a jet solution, we will be able to ascertain the conditions that best describe the overall spectrum of a given black hole system. In future work we will explore how well this model succeeds in predicting the correct location of the optically thick-to-thin break observed in the broad-band spectra of compact jets, hopefully with new insights into the physics of jet launching conditions.

PP and SM gratefully acknowledge support from a Netherlands Organization for Scientific Research (NWO) Vidi Fellowship. In addition, SM is grateful for support from the European Community's Seventh Framework Programme (FP7/2007-2013) under grant agreement number ITN 215212 ‘Black Hole Universe’. Part of the research described in this paper was carried out at the Jet Propulsion Laboratory, California Institute of Technology, under a contract with the National Aeronautics and Space Administration.

REFERENCES

Asada
K.
Nakamura
M.
Doi
A.
Nagai
H.
Inoue
M.
Romero
G. E.
Sunyaev
R. A.
Belloni
T.
Proc. IAU Symp. 275, EVN Monitoring Observation of M87 Jet
 , 
2011
Cambridge
Cambridge Univ. Press
pg. 
198
 
Bell
A. R.
MNRAS
 , 
1978
, vol. 
182
 pg. 
147
 
Blandford
R. D.
Begelman
M. C.
MNRAS
 , 
1999
, vol. 
303
 pg. 
L1
 
Blandford
R. D.
Konigl
A.
ApJ
 , 
1979
, vol. 
232
 pg. 
34
 
Blandford
R. D.
Payne
D. G.
MNRAS
 , 
1982
, vol. 
199
 pg. 
883
  
BP82
Corbel
S.
Fender
R. P.
ApJ
 , 
2002
, vol. 
573
 pg. 
L35
 
Drury
L. O.
Rep. Prog. Phys.
 , 
1983
, vol. 
46
 pg. 
973
 
Esin
A. A.
McClintock
J. E.
Narayan
R.
ApJ
 , 
1997
, vol. 
489
 pg. 
865
 
Falcke
H.
Biermann
P. L.
A&A
 , 
1995
, vol. 
293
 pg. 
665
 
Falcke
H.
Markoff
S.
A&A
 , 
2000
, vol. 
362
 pg. 
113
 
Falcke
H.
Körding
E.
Markoff
S.
A&A
 , 
2004
, vol. 
414
 pg. 
895
 
Fender
R. P.
MNRAS
 , 
2001
, vol. 
322
 pg. 
31
 
Fender
R. P.
Belloni
T. M.
Gallo
E.
MNRAS
 , 
2004
, vol. 
355
 pg. 
1105
 
Gallo
E.
Migliari
S.
Markoff
S.
Tomsick
J. A.
Bailyn
C. D.
Berta
S.
Fender
R.
Miller-Jones
J. C. A.
ApJ
 , 
2007
, vol. 
670
 pg. 
600
 
Gandhi
P.
, et al. 
ApJ
 , 
2011
, vol. 
740
 pg. 
L13
 
Gebhardt
K.
Thomas
J.
ApJ
 , 
2009
, vol. 
700
 pg. 
1690
 
Hada
K.
Doi
A.
Kino
M.
Nagai
H.
Hagiwara
Y.
Kawaguchi
N.
Nat
 , 
2011
, vol. 
477
 pg. 
185
 
Heinz
S.
Sunyaev
R. A.
MNRAS
 , 
2003
, vol. 
343
 pg. 
L59
 
Ho
L. C.
ApJ
 , 
1999
, vol. 
516
 pg. 
672
 
Junor
W.
Biretta
J. A.
Livio
M.
Nat
 , 
1999
, vol. 
401
 pg. 
891
 
Li
Z.-Y.
Chiueh
T.
Begelman
M. C.
ApJ
 , 
1992
, vol. 
394
 pg. 
459
 
Maitra
D.
Markoff
S.
Brocksopp
C.
Noble
M.
Nowak
M.
Wilms
J.
MNRAS
 , 
2009
, vol. 
398
 pg. 
1638
 
Markoff
S.
Falcke
H.
Fender
R.
A&A
 , 
2001
, vol. 
372
 pg. 
L25
 
Markoff
S.
Nowak
M.
Corbel
S.
Fender
R.
Falcke
H.
A&A
 , 
2003
, vol. 
397
 pg. 
645
 
Markoff
S.
Nowak
M. A.
Wilms
J.
ApJ
 , 
2005
, vol. 
635
 pg. 
1203
 
Markoff
S.
, et al. 
ApJ
 , 
2008
, vol. 
681
 pg. 
905
 
Marscher
A. P.
Gear
W. K.
ApJ
 , 
1985
, vol. 
298
 pg. 
114
 
McClintock
J. E.
Remillard
R. A.
Black Hole Binaries
 , 
2006
Cambridge
Cambridge Univ. Press
McKinney
J. C.
MNRAS
 , 
2006
, vol. 
368
 pg. 
1561
 
McKinney
J. C.
Blandford
R. D.
MNRAS
 , 
2009
, vol. 
394
 pg. 
L126
 
Meier
D. L.
ApJ
 , 
2001
, vol. 
548
 pg. 
L9
 
Merloni
A.
Fabian
A. C.
MNRAS
 , 
2002
, vol. 
332
 pg. 
165
 
Merloni
A.
Heinz
S.
di Matteo
T.
MNRAS
 , 
2003
, vol. 
345
 pg. 
1057
 
Migliari
S.
Tomsick
J. A.
Maccarone
T. J.
Gallo
E.
Fender
R. P.
Nelemans
G.
Russell
D. M.
ApJ
 , 
2006
, vol. 
643
 pg. 
L41
 
Migliari
S.
, et al. 
ApJ
 , 
2007
, vol. 
670
 pg. 
610
 
Nakamura
M.
Meier
D. L.
ApJ
 , 
2004
, vol. 
617
 pg. 
123
 
Narayan
R.
Yi
I.
ApJ
 , 
1994
, vol. 
428
 pg. 
L13
 
Paczyńsky
B.
Wiita
P. J.
A&A
 , 
1980
, vol. 
88
 pg. 
23
 
Plotkin
R. M.
Markoff
S.
Kelly
B. C.
Körding
E.
Anderson
S. F.
MNRAS
 , 
2012
, vol. 
419
 pg. 
267
 
Polko
P.
Meier
D. L.
Markoff
S.
ApJ
 , 
2010
, vol. 
723
 pg. 
1343
  
PMM10; Paper I
Quataert
E.
Gruzinov
A.
ApJ
 , 
2000
, vol. 
539
 pg. 
809
 
Rahoui
F.
Lee
J. C.
Heinz
S.
Hines
D. C.
Pottschmidt
K.
Wilms
J.
Grinberg
V.
ApJ
 , 
2011
, vol. 
736
 pg. 
63
 
Shakura
N. I.
Sunyaev
R. A.
A&A
 , 
1973
, vol. 
24
 pg. 
337
 
Tonry
J. L.
Dressler
A.
Blakeslee
J. P.
Ajhar
E. A.
Fletcher
A. B.
Luppino
G. A.
Metzger
M. R.
Moore
C. B.
ApJ
 , 
2001
, vol. 
546
 pg. 
681
 
Vlahakis
N.
Königl
A.
ApJ
 , 
2003
, vol. 
596
 pg. 
1080
  
VK03
Vlahakis
N.
Tsinganos
K.
Sauty
C.
Trussoni
E.
MNRAS
 , 
2000
, vol. 
318
 pg. 
417
  
VTST00
Walker
R. C.
Ly
C.
Junor
W.
Hardee
P. J.
J. Phys. Conf. Ser.
 , 
2008
, vol. 
131
 pg. 
012053
 

APPENDIX A: DERIVATION OF THE GRAVITY TERM

To obtain a relativistic form of gravity, we need the following conversions: ϖ* in VTST00 is equal to ϖ0 in VK03 and V* similarly corresponds to $$\frac{\alpha ^{1/4} F \sigma _\mathrm{M}c}{\gamma \xi x_\mathrm{A}^2}$$.

Substituting these into equation (17) of VTST00 produces  
\begin{eqnarray} \kappa &=& \displaystyle\sqrt{\frac{\mathcal {G M}}{\varpi _* V_*^2}} = \sqrt{\frac{\mathcal {G M} \gamma ^2 \xi ^2 x_\mathrm{A}^4}{\varpi _0 \alpha ^{1/2} F^2 \sigma _\mathrm{M}^2 c^2}} \nonumber \\ &=& \displaystyle\sqrt{\frac{\mathcal {G M}}{c^2 \varpi _\mathrm{A}} \frac{\mu ^2 x_\mathrm{A}^4}{F^2 \sigma _\mathrm{M}^2} \frac{(1 - M^2 - x_\mathrm{A}^2)^2}{(1 - M^2 - x^2)^2}}. \end{eqnarray}
(A1)
The full gravity term then becomes  
\begin{equation} - \frac{\kappa ^2 \sin (\theta )}{G} = - \frac{\mathcal {G M}}{c^2} \frac{\mu ^2 x_\mathrm{A}^4}{F^2 \sigma _\mathrm{M}^2} \frac{\left(1 - M^2 - x_\mathrm{A}^2\right)^2}{(1 - M^2 - x^2)^2} \frac{\sin (\theta )}{\varpi _\mathrm{A}G}. \end{equation}
(A2)
In this expression, (ϖAG)/sin (θ) is equal to the spherical radius, r. To include a pseudo-Newtonian potential, we replace 1/r by 1/(rrS), where rS is the Schwarzschild radius ($$= 2 \mathcal {G M}/c^2$$), and simplify  
\begin{eqnarray} - \displaystyle\frac{\kappa ^2 \sin (\theta )}{G} = - \frac{\mu ^2 x_\mathrm{A}^4}{F^2 \sigma _\mathrm{M}^2} \frac{\left(1 - M^2 - x_\mathrm{A}^2\right)^2}{(1 - M^2 - x^2)^2} \left[ \frac{c^2 \varpi _\mathrm{A}G}{\mathcal {G M} \sin (\theta )} - 2 \right]^{-1}. \nonumber \\ \end{eqnarray}
(A3)

APPENDIX B: EQUATIONS FOR THE INITIAL PARAMETER VALUES

The ARC is given by evaluating the wind equation at the AP,  
\begin{eqnarray} &&\frac{F^2 \sigma _\mathrm{M}^2 (1 - x_\mathrm{A}^2) (\sigma _\mathrm{A}+ 1)^2 \sin (\theta _\mathrm{A})}{\mu ^2 \cos ^2(\theta _\mathrm{A}+ \psi _\mathrm{A})} \bigg \lbrace \nonumber \\ && -\, 2 \frac{\Gamma - 1}{\Gamma } \frac{(F - 2)(\xi _\mathrm{A}- 1) (1 - x_\mathrm{A}^2)}{\xi _\mathrm{A}x_\mathrm{A}^2} \sin (\theta _\mathrm{A}) \nonumber \\ && +\, 2 \cos (\psi _\mathrm{A}) \sin (\theta _\mathrm{A}+ \psi _\mathrm{A}) \frac{\sigma _\mathrm{A}+ 1}{\sigma _\mathrm{A}} \nonumber \\ && +\, \frac{\sin (\theta _\mathrm{A})}{x_\mathrm{A}^2} \left[ (F - 1) (1 - x_\mathrm{A}^2) - 1\right] \bigg \rbrace \nonumber \\ &=& \left[ x_\mathrm{A}^2 - \sigma _\mathrm{A}(1 - x_\mathrm{A}^2) \right]^2 - (F - 1) \sigma _\mathrm{A}^2 (1 - x_\mathrm{A}^2) \nonumber \\ && -\, 2 \frac{\Gamma - 1}{\Gamma } (F - 2) \frac{\xi _\mathrm{A}- 1}{\xi _\mathrm{A}} \left\lbrace x_\mathrm{A}^2 - \left[ x_\mathrm{A}^2 - \sigma _\mathrm{A}(1 - x_\mathrm{A}^2) \right]^2 \right\rbrace \nonumber \\ && -\, \frac{x_\mathrm{A}^2}{1 - x_\mathrm{A}^2} \left[ \frac{c^2 \varpi _\mathrm{A}}{\mathcal {G M} \sin (\theta _\mathrm{A})} - 2 \right]^{-1}. \end{eqnarray}
(B1)
We obtain μ2 by evaluating the Bernoulli equation at the AP  
\begin{eqnarray} \mu ^2 &=& \displaystyle\frac{(\sigma _\mathrm{A}+1)^2}{x_\mathrm{A}^2 - \left[ x_\mathrm{A}^2 - \sigma _\mathrm{A}\left( 1 - x_\mathrm{A}^2 \right) \right]^2} \bigg [ x_\mathrm{A}^2 \xi _\mathrm{A}^2 \nonumber \\ && +\, \displaystyle\frac{F^2 \sigma _\mathrm{M}^2 \left( 1 - x_\mathrm{A}^2 \right)^2 \sin ^2(\theta _\mathrm{A})}{x_\mathrm{A}^2 \cos ^2(\theta _\mathrm{A}+ \psi _\mathrm{A})} \bigg ]. \end{eqnarray}
(B2)

APPENDIX C: GRAVITY IN THE BERNOULLI EQUATION

The only terms that include the effects of gravity in the Bernoulli equation given by (A3) in VTST00 are the Bernoulli constant and the gravity term. Dividing out a common factor of 2, we have  
\begin{equation} \varepsilon + \frac{\kappa ^2 \sin (\theta )}{G}. \end{equation}
(C1)
At the AP, G = 1 and θ = θA, reducing these terms to  
\begin{equation} \varepsilon + \kappa ^2 \sin (\theta _\mathrm{A}). \end{equation}
(C2)
If we write ε using equation (19) of VTST00, G ≡ ϖ/ϖA, and equation (2.7a) of BP82, we obtain (please note that to keep to the notation of the preceding two papers, here the subscript 0 means the value at the base of the outflow, not the reference values as used in VK03; unfortunately, r0 as used in BP82 is ϖ0 as used in VTST00)  
\begin{equation} \varepsilon = \varepsilon \frac{\kappa ^2}{G_0} = \frac{e}{\mathcal {G M} / \varpi _0} \frac{\kappa ^2}{\varpi _0 / \varpi _\mathrm{A}} = \kappa ^2 \frac{e}{\mathcal {G M} / \varpi _\mathrm{A}}. \end{equation}
(C3)
Since e, the specific energy, is a constant of motion, we can evaluate it at the AP, just like the Bernoulli equation in VK03 (equation 2.2 of BP82 with the gravitational potential expanded):  
\begin{equation} e = e_\mathrm{A}= \frac{V_\mathrm{A}^2}{2} + h_\mathrm{A}- \frac{\mathcal {G M}}{r_\mathrm{A}} - \frac{\Omega \varpi _\mathrm{A}B_{\phi ,\mathrm{A}}}{\Psi _{A,\mathrm{A}}}. \end{equation}
(C4)
Here, rA is the spherical radius of the AP. Substituting (C4) and (C3) into (C2), we obtain  
\begin{equation} \kappa ^2 \left[ \frac{\varpi _\mathrm{A}}{\mathcal {G M}} \left( \frac{V_\mathrm{A}^2}{2} + h_\mathrm{A}- \frac{\Omega \varpi _\mathrm{A}B_{\phi ,\mathrm{A}}}{\Psi _{A,\mathrm{A}}} \right) - \frac{\varpi _\mathrm{A}}{r_\mathrm{A}} + \sin (\theta _\mathrm{A}) \right]. \end{equation}
(C5)
However, ϖA/rA = sin (θ), so the last two terms cancel. The κ2 term cancels the $$\mathcal {G M}$$ factor (equation A1), leaving no terms that depend on gravity. Thus, there is no dependence on gravity in the Bernoulli equation in the non-relativistic regime with gravity at the AP, and consequently there is no gravity term in the relativistic Bernoulli equation at the AP (equation B2).