TY - JOUR
AU - Tafoya,, D
AB - ABSTRACT Outflows, spanning a wide range of dynamical properties and spatial extensions, have now been associated with a variety of accreting astrophysical objects, from supermassive black holes at the core of active galaxies to young stellar objects. The role of such outflows is key to the evolution of the system that generates them, for they extract a fraction of the orbiting material and angular momentum from the region close to the central object and release them in the surroundings. The details of the launching mechanism and their impact on the environment are fundamental to understand the evolution of individual sources and the similarities between different types of outflow-launching systems. We solve semi-analytically the non-relativistic, ideal, magnetohydrodynamics equations describing outflows launched from a rotating disc threaded with magnetic fields using our new numerical scheme. We present here a parameter study of a large sample of new solutions. We study the different combinations of forces that lead to a successfully launched jet and discuss their global properties. We show how these solutions can be applied to the outflow of the water fountain W43A for which we have observational constraints on magnetic field, density and velocity of the flow at the location of two symmetrical water maser emitting regions. MHD, stars: AGB and post-AGB, stars: winds, outflows, stars: jets 1 INTRODUCTION Jets, and more generally speaking outflows, are a widespread phenomena in many different systems, from protostars to supermassive black holes. Until only recently, outflows were believed to be launched by extraction of rotational energy either from a magnetically threaded disc (Blandford & Payne 1982) or from a rotating black hole (Blandford & Znajek 1977). In fact, it was common to think that there were at least two separated kinds of jets, the magnetically dominated jets from black holes and the pressure-dominated jets from any other jetted source. However, with the advent of cutting-edge GRMHD simulations and the first post-processed emission spectrum associated with it (Mościbrodzka, Falcke & Shiokawa 2016; Liska et al. 2017; Davelaar et al. 2019), it is becoming clear that there is not such a dichotomy and, most likely, at least in black hole systems, the two mechanisms can coexist. Furthermore, the emission is likely dominated by the outer, more mass loaded, jet sheath rooted on to the accretion disc, whereas the inner core of the jet is lighter and magnetically dominated (Mościbrodzka et al. 2016). Similarly, when the jet-launching object is a protostar or a (non-BH) compact object, the outflow is likely to be a composition of a stellar wind (e.g. Shu et al. 1994) or an equivalent Blandford–Znajek process for highly magnetized neutron stars (Parfrey, Spitkovsky & Beloborodov 2016) and a disc-driven outflow (e.g. Pudritz & Norman 1983; Contopoulos & Lovelace 1994; Ferreira 1997; Vlahakis et al. 2000). Semi-analytical models describing the various launching mechanisms listed above have been continuously developed in parallel with simulations because they capture the underlying physics while allowing a time-efficient exploration of the parameter space and fitting of astrophysical sources. However, in order to make the equations treatable with a semi-analytical approach, the dimensionality of the problem is reduced by assuming symmetries in the system and a non-linear separation of variables is performed. The separation of variables is commonly referred to as the self-similarity assumption. There are two distinct classes of self-similar models depending on how the separation of variables is carried out (Vlahakis & Tsinganos 1998): the meridional self-similar models (e.g. Sauty & Tsinganos 1994; Trussoni, Tsinganos & Sauty 1997; Sauty, Tsinganos & Trussoni 1999; Chantry et al. 2018), where the dependent variables are functions of r and radial self-similar models (e.g. Contopoulos & Lovelace 1994; Ferreira 1997; Vlahakis et al. 2000; Vlahakis & Königl 2003; Polko, Meier & Markoff 2010, 2013, 2014; Ceccobello et al. 2018), where the independent variable is θ. In both classes of self-similar models, we are left with a mixed system of differential and algebraic equations describing the accelerating flow along a magnetic field line threading a rotating disc. To determine the motion of the fluid element, one needs to solve simultaneously the forces acting along and perpendicular to a streamline. This is a notoriously cumbersome problem, which can be tackled by introducing further simplifications, such as assuming a fixed structure of the magnetic field and/or neglecting the gas pressure force and/or using asymptotic extensions of the models to replace the region of the solutions at large distance from the disc. In Ceccobello et al. (2018, hereafter Paper I), we presented our newly developed algorithm to self-consistently solve the poloidal and transverse forces, given by the Bernoulli and Grad–Shafranov equations, respectively, for a relativistic fluid in the presence of gravity, under the assumption of radial self-similarity. We showed that with our numerical algorithm it is possible to obtain solutions with a broad variety of jet structures and dynamical properties and work is ongoing to couple these solutions with a radiative code and apply those to black hole systems (Lucchini et al., in preparation). In this paper, we adopt the equations presented in Vlahakis et al. (2000, hereafter VTST00) and we adapt our algorithm, described in Paper I, to perform a parameter study to model astrophysical sources with more moderate speeds, such as young stellar objects and evolved stars outflows. In Section 2, we summarize the basic equations and give a short description of the algorithm. In Section 3, we show the results of our parameter space exploration and discuss the solution properties as they transition from cold jets to hot ones. In Section 4, we show an example of an application to the post-AGB star W43A and we give the selection criteria we used to isolate the solutions that better resemble the jet of W43A and discuss the characteristics of the selected jet configuration in relation to the source. Finally, in Section 5, we summarize the study presented in this paper. 2 EQUATIONS AND NUMERICAL METHOD 2.1 Problem description The equations that we are going to solve with our numerical algorithm are the ones describing an axisymmetric, radial self-similar, non-relativistic, disc-driven outflow with non-negligible enthalpy (Fig. 1). Since we adopted the prescription given in Vlahakis et al. (2000),1 we present here just a brief summary. In Appendix B, we report the conversion from dimensionless to physical quantities as a function of the input parameters and the scaling relations. The dependent variables of the equations described in VTST00 are the poloidal Mach number M, the dimensionless cylindrical radius, G, and the angle describing the inclination of the streamline with respect to the disc plane, ψ. These are all functions of θ once their radial dependence has been defined as power laws of the function |$\alpha \equiv \varpi _{\rm A}^2/\varpi _*^2$|⁠, where ϖA is the cylindrical radius at the Alfvén point (AP) and ϖ* is the chosen scaling length of the problem and effectively is the cylindrical radius of the AP on the streamline with α = 1 (see Fig. 1). Figure 1. Open in new tabDownload slide System of coordinates we adopt to describe a solution of the MHD system of equations presented in VTST. We identify a solution with the ‘reference’ streamline identified with the label |$\alpha \equiv \varpi _{\rm A}^2/\varpi _*^2 = 1$|⁠. The two dependent variables (M, G), together with all the other quantities describing the system, are functions of θ (the independent variable in radial self-similarity), which is the angle between a point on the streamline and the z-axis. The angle ψ is defined by the tangent to the streamline and the horizontal axis, while the distance from a point on the streamline to the z-axis is defined by its cylindrical radius ϖ in units of the Alfvén cylindrical radius ϖA. Figure 1. Open in new tabDownload slide System of coordinates we adopt to describe a solution of the MHD system of equations presented in VTST. We identify a solution with the ‘reference’ streamline identified with the label |$\alpha \equiv \varpi _{\rm A}^2/\varpi _*^2 = 1$|⁠. The two dependent variables (M, G), together with all the other quantities describing the system, are functions of θ (the independent variable in radial self-similarity), which is the angle between a point on the streamline and the z-axis. The angle ψ is defined by the tangent to the streamline and the horizontal axis, while the distance from a point on the streamline to the z-axis is defined by its cylindrical radius ϖ in units of the Alfvén cylindrical radius ϖA. To obtain a full solution, i.e. a streamline rooted at the disc mid-plane and terminating infinity, the adopted numerical scheme must handle three singular points that are present in the Bernoulli and Grad–Shafranov equations when solved simultaneously: the AP and the magnetosonic fast/slow points (MFP/MSP).2 At each singular point the equations can be regularized either analytically, in the case of the AP, or numerically, for the MSP and the MFP, analogously to the simpler case of the sonic point in the Parker wind model (Parker 1958). The AP has been studied extensively, due to the possibility of manipulating the equations analytically there. The other two singular points present a more complex case. On the one hand, the position of both the MSP and the MFP is not known before the full solution for a given set of initial parameters is calculated, on the other hand, a full solution cannot be computed without knowing the position of these two singular points and the AP. Due to this intrinsic difficulty, the MSP and MFP are often neglected by assuming cold flows, i.e. thermal pressure plays no role in accelerating the flow (no MSP), and/or by adopting a given asymptotic behaviour of the streamline once the flow has become superalfvenic, which effectively pushes the MFP at infinity. Typically, either one or both of the above assumptions are made to avoid dealing with the complexity of determining these singular points. Moreover, when the MSP and/or the MFP are not removed from the equations, finding solutions across large volumes of the parameter space is a difficult task that requires a solid numerical algorithm capable of recovering the unknown positions of the singular points and properly handling the equations at these locations for wide ranges of the input parameters. However, the role of the MFP in self-similar theories is fundamental when solving the Bernoulli and Grad–Shafranov equations combined, because it is the singular point where the flow loses causal contact with the source (Li, Chiueh & Begelman 1992; Bogovalov & Tsinganos 1999; Meier 2012). Downstream of the MFP the flow starts to focus rapidly towards the polar axis up until the last recollimation point (LRP). We identify the LRP with the region where the jet terminates in our solutions (see Paper I). This region has been connected in relativistic jets with the standing shock/particle acceleration regions in active galactic nuclei and in stellar-mass black hole systems (e.g. Markoff, Falcke & Fender 2001; Markoff, Nowak & Wilms 2005; Markoff 2010; Polko et al. 2010; Meier 2012; Cohen et al. 2014; Ceccobello et al. 2018). Weber & Davis (1967) showed that there can be multiple families of solutions with different velocity profiles, crossing either none or one/two/three singular points. We are looking at those that cross all three points, which are characterized by an increasing poloidal Mach number. VTST00 were the first to calculate complete solutions with all these characteristics for the non-relativistic case. 2.2 Non-relativistic MHD system of equations The Bernoulli and Grad–Shafranov equations for a steady-state axisymmetric system describe the energy flux balance along the poloidal direction and the equilibrium configuration of the magnetic field lines. Both can be derived from the conservation of momentum equation, which describes the forces acting on a streamline: $$\begin{eqnarray*} \rho (\boldsymbol{V}\cdot \nabla) \boldsymbol{V} -\frac{1}{4\pi }\boldsymbol{B} \times (\nabla \times \boldsymbol{B}) + \nabla P -\rho \nabla \frac{\mathcal {G}M}{r} =0 , \end{eqnarray*}$$(1) where |$\rho , P, \boldsymbol{V}, \boldsymbol{B}$| are the density, pressure, velocity, and magnetic field of the flow. |$\mathcal {G}$| and |$\mathcal {M}$| are the gravitational constant and the mass of the central object, respectively. If we adopt either cylindrical (⁠|$\boldsymbol{\hat{z}}, \boldsymbol{\hat{\varpi }}, \boldsymbol{\hat{\phi}}$|⁠) or spherical coordinates (⁠|$\boldsymbol{\hat{r}}, \boldsymbol{\hat{\theta }}, \boldsymbol{\hat{\phi}}$|⁠), the poloidal and perpendicular unit vectors (⁠|$\hat{\boldsymbol{b}}$|⁠, |$\hat{\boldsymbol{n}}$|⁠) can be written as follows: $$\begin{eqnarray*} &\hat{\boldsymbol{n}} = \cos (\psi)\boldsymbol{\hat{z}} - \sin (\psi) \boldsymbol{\hat{\varpi }} = \cos (\theta +\psi)\hat{\boldsymbol{r}} - \sin (\theta +\psi)\hat{\boldsymbol{\theta }} , \end{eqnarray*}$$(2) $$\begin{eqnarray*} &\hat{\boldsymbol{b}} = \sin (\psi)\boldsymbol{\hat{z}} + \cos (\psi) \boldsymbol{\hat{\varpi }} = \sin (\theta +\psi)\hat{\boldsymbol{r}}+ \cos (\theta +\psi)\hat{\boldsymbol{\theta }} , \end{eqnarray*}$$(3) $$\begin{eqnarray*} &\hat{\boldsymbol{\phi }} = \hat{\boldsymbol{\phi }}. \end{eqnarray*}$$(4) The projection of equation (1) along |$\hat{\boldsymbol{b}}$|⁠, the Bernoulli equation, describes how the different types of energies can be converted to one another. The projection of equation (1) along |$\hat{\boldsymbol{n}}$|⁠, the Grad–Shafranov equation or transfield equation, provides the shape of the magnetic field lines. The projections of the Bernoulli and the transfield equation can be rewritten using the scaling equations given in Appendix B and then rearranged in the following form: $$\begin{eqnarray*} A_i\frac{\mathrm{ d}M^2}{\mathrm{ d}\theta } + B_i\frac{\mathrm{ d}\psi }{\mathrm{ d}\theta } = C_i, \end{eqnarray*}$$(5) with i = 1 representing the coefficients of the Bernoulli equation and i = 2 the coefficients of the transfield equation. The Bernoulli and transfield equations arranged in the way described above can further be recast into a system of two first-order differential equations for the evolution of the poloidal Mach number |$M(\theta) = \sqrt{4\pi \rho } V_p/B_p$| and the angle ψ describing the inclination of the streamline with respect to the horizontal axis: $$\begin{eqnarray*} \frac{\mathrm{ d}M^2}{\mathrm{ d}\theta } = \frac{\mathcal {N_{\rm 1}}}{\mathcal {D}} = \frac{B_2 C_1 - B_1 C_2}{ A_1B_2- A_2 B_1} , \end{eqnarray*}$$(6) $$\begin{eqnarray*} \frac{\mathrm{ d}\psi }{\mathrm{ d}\theta } = \frac{\mathcal {N_{\rm 2}}}{\mathcal {D}} = \frac{ A_1 C_2 - A_2 C_1}{ A_1B_2- A_2 B_1}, \end{eqnarray*}$$(7) with the numerators |$\mathcal {N}_i$| (i = 1, 2) and the denominator |$\mathcal {D}$| being functions of the coefficients Ai, Bi, Ci (i = 1, 2) which are given in Appendix A. As described in Paper I, in order to minimize the intrinsic errors we chose not to solve equation (7), but instead derive ψ(θ) from the Bernoulli integral equation (12) (see also VTST00, and Appendix A) from the MSP to the LRP. Upstream of the MSP, the streamlines can undergo oscillations, depending on the given set of input parameters, so the sign of cos (ψ + θ) can change. Hence, equation (7) must be integrated with care in this region to ensure the correct radial profile of the solutions from the disc to the MSP. Additionally, we solve a differential equation for the unknown function G(θ), which is defined as the cylindrical radius to the polar axis of a streamline labelled by α, normalized to its cylindrical radius at the AP. The equation for G is the following: $$\begin{eqnarray*} G(\theta) = \frac{\varpi }{\varpi _\alpha } = \frac{\varpi }{\varpi _\star } \alpha ^{-1/2}, \end{eqnarray*}$$(8) $$\begin{eqnarray*} \frac{\mathrm{ d}G^2}{\mathrm{ d}\theta } = \frac{2 G^2\cos (\psi)}{\sin (\theta)\cos (\psi +\theta)}. \end{eqnarray*}$$(9) The solution of these equations depends on six parameters: Γ, F, kVTST, λVTST, μVTST, and ϵVTST (see VTST00). The first parameter Γ is the polytropic index in the equation of state q = P/ρΓ, where q is the specific gas entropy and a constant of motion of the problem. The parameter F determines the initial current distribution in the radial direction, |$-\varpi B_\phi = \mathcal {C}_2 \sin (\theta) r^{F-1}$|⁠, which is an increasing or decreasing function of r depending on the value of F. This parameter also determines the radial dependence of the magnetic field lines through B ∼ ϖF − 2. kVTST is proportional to the ratio between the Keplerian speed and the poloidal flow speed at the Alfvén radial distance, and often is referred to as the mass-loss parameter (see e.g. Ferreira 1997, but also VTST00). λVTST is the specific angular momentum in units of V⋆ϖ⋆ and μVTST is proportional to the gas entropy. The parameters kVTST, μVTST, and λVTST are defined by the following relations: $$\begin{eqnarray*} k_{\rm VTST} = \sqrt{\frac{\mathcal {GM}}{\varpi _\star V_\star ^2}}; \quad \mu _{\rm VTST} = \frac{8 \pi P_\star }{B_\star ^2}; \quad \lambda _{\rm VTST} = \frac{L}{V_\star \varpi _\star }. \end{eqnarray*}$$(10) It is worth noticing that the starred quantities found across the paper are scaling factors and can be related to the quantities calculated at the AP on the reference streamline (α = 1), namely ρ⋆ = ρA, ϖ⋆ = ϖA, and $$\begin{eqnarray*} \left(B_*,V_*\right) = -\frac{\cos (\theta _{\rm A} + \psi _{\rm A})}{\sin (\theta _{\rm A})}\left(B_{p,\rm A}, V_{p,\rm A}\right), \quad {\rm with}\ B_* = \sqrt{4\pi \rho _*} V_* \nonumber \\ \end{eqnarray*}$$(11) Finally, ϵVTST is the sum of kinetic, enthalpy, gravitational, and Poynting energy flux densities per unit of mass flux density, rescaled by |$\alpha ^{-1/2}V^{2}_\star$|⁠, i.e. $$\begin{eqnarray*} \epsilon _{\rm VTST} &=& \frac{\alpha ^{1/2}}{V_\star ^2} E = \left[\epsilon _{\rm {K}, p} + \epsilon _{\rm {K}, \phi } +\epsilon _{\rm {T}} +\epsilon _{\rm {M} }+ \epsilon _{\rm G} \right] \nonumber \\ &=& \left[\frac{1}{2} \left(\frac{M^2}{G^2}\frac{\sin (\theta)}{\cos (\theta +\psi)}\right)^2 + \frac{1}{2}\left(\frac{\lambda _{\rm VTST}}{G^2}\frac{G^2 - M^2}{1-M^2}\right)^2 \right. \nonumber \\ && \left. +\, \frac{\mu _{\rm VTST}}{2} \frac{\Gamma }{\Gamma -1}M^{2(1-\Gamma)} + \lambda _{\rm VTST}^2\frac{1-G^2}{1-M^2}- k_{\rm VTST}^2 \frac{\sin (\theta)}{G}\right]. \nonumber \\ \end{eqnarray*}$$(12) The total energy flux per unit mass can be rescaled with the Alfvén poloidal velocity as |$2 E/V_{\mathrm{ A},p}^2 = 2\epsilon _{\rm VTST} V_*^2\alpha ^{-1/2}/V_{\mathrm{ A},p}^2$|⁠, which becomes $$\begin{eqnarray*} \tilde{\epsilon }= 2\epsilon _{\rm VTST} \frac{\cos ^2(\theta _{\rm A}+\psi _{\rm A})}{\sin ^2(\theta _{\rm A})} \end{eqnarray*}$$(13) and with the use of the De L’Hôpital rule to regularize the indefinite terms (see equation 18), we can write it at the AP and obtain the Alfvén Regularity Condition (ARC, see VTST00) in the compact form $$\begin{eqnarray*} \tilde{\epsilon } &=& 1 + \frac{\cos ^2(\theta _{\rm A}+\psi _{\rm A})}{\sin ^2(\theta _{\rm A})}\left[- 2k_{\rm VTST}^2 \sin (\theta _{\rm A}) + \mu _{\rm VTST}\frac{\Gamma }{\Gamma -1} \right.\nonumber \\ && \left. +\lambda _{\rm VTST}^2\left(1+g_{\rm A}^2\right)\right]. \end{eqnarray*}$$(14) The function gA is the fastness parameter calculated at the AP. A general definition of the fastness parameter given by Pelletier & Pudritz (1992) is $$\begin{eqnarray*} \frac{V_\phi }{\varpi } = \Omega (1- g) , \end{eqnarray*}$$(15) where $$\begin{eqnarray*} \Omega = \frac{1}{\varpi } \left(V_\phi - V_p \frac{B_\phi }{B_p}\right) \end{eqnarray*}$$(16) is the angular frequency of the streamline, which is a constant of motion of the problem. The fastness parameter gives a measure of how large the angular velocity of the gas is in relation to the angular velocity of the magnetic surface on which it moves. We can derive gA from the application of the De L’Hôpital rule to the indefinite forms $$\begin{eqnarray*} \left.\frac{1-G^2}{1-M^2}\right|_{\rm A} \equiv g_{\rm A} = \frac{2 \cos (\psi _{\rm A})}{p_{\rm A} \sin (\theta _{\rm A})\cos (\theta _{\rm A}+\psi _{\rm A})} , \end{eqnarray*}$$(17) $$\begin{eqnarray*} \left.\frac{G^2-M^2}{1-M^2}\right|_{\rm A}& = 1 - g_{\rm A} , \end{eqnarray*}$$(18) where pA = dM2/dθ|A. In the following section, we summarize the method we developed in Paper I that we now adapt to solve the non-relativistic equations. For the details of the algorithm, we address the interested reader to Paper I. Indeed, there is no substantial difference in the mechanics of the algorithm, although the non-relativistic equations are noticeably easier to handle. 2.3 Method In Paper I, we described a new numerical method to find solutions to the relativistic radial self-similar magnetohydrodynamics (MHD) equations for a disc-launched jet in the presence of gravity (Vlahakis & Königl 2003; Polko et al. 2014). As discussed in Section 2.1, even under the simplifying assumption of self-similarity, solving self-consistently and simultaneously the Bernoulli and Grad–Shafranov equations is known to be a rather difficult task because of the singular surfaces. At the location of the singular points, the equations (6) and (7) are indeterminate but finite, e.g. $$\begin{eqnarray*} \frac{\mathrm{ d}M^2}{\mathrm{ d}\theta } = \frac{\mathcal {N}_1}{\mathcal {D}} = \frac{0}{0} = {\rm finite}. \end{eqnarray*}$$(19) However, only at the AP one can derive an analytical expression that gives the finite value of the derivative of the poloidal Mach number (ARC). The location of the AP and |$G^2_{\rm A}, M^2_{\rm A}, \mathrm{ d}G^2/\mathrm{ d}\theta |_{\rm A}$| and dM2/dθ|A can be determined from the values of the input parameters and the ARC (equation 14). The regularity conditions at the MFP and MSP can exclusively be derived numerically together with their position on the streamline. As a result, the most frequent approach is to determine all the unknown functions and parameters at AP and then integrate the system with a shooting method towards the other two singular points. However, given the high accuracy needed to determine the values of the parameters and the intrinsic numerical difficulties of treating, under these conditions, the form 0/0, this method presents serious drawbacks and does not allow to easily find and convincingly identify solutions to the required accuracy threshold. Therefore, it impedes a full exploration of the parameter space. The structure of our numerical method is the following: We guess the locations of the critical points, θMSP and θMFP, and derive values for M2, G2, and their derivatives given by the condition that the numerators and the denominator of equation (19), and of the similar equation for ψ, i.e. |$\mathrm{ d}\psi /\mathrm{ d}\theta = \mathcal {N}_2/\mathcal {D} = 0/0$|⁠, are zero at the MSP/MFP of choice. We integrate away from AP, MSP, and MFP towards the mid-points θmid,MSP = (θA + θMSP)/2, and θmid,MFP = (θA + θMFP)/2. We determine the parameters that give a match at the mid-points using the Bayesian open-source code multinest (Feroz & Hobson 2008; Feroz, Hobson & Bridges 2009; Feroz et al. 2013). The specific choice of input parameters and fitted parameters is given in Table 1. Once a particular family of solutions is specified through the choice of F, Γ, θA, ψA, and kVTST, we identify the location of the MSP and MFP and the best-fitting values of the remaining parameters μVTST and λVTST and we extend the solutions upstream of the MSP towards the disc mid-plane and downstream of the MFP towards the LRP. In Paper I, we defined this point as the last point we were able to calculate with our algorithm. The last few integration points before LRP seem to indicate the onset of a recollimation shock where the fluid is compressed in a small section around the polar axis. Indeed, we noticed that the denominator is approaching zero again in equations (6) and (7), while the numerator is not. This means that both the derivative of M2 and ψ become infinite close to LRP, making the integration towards this (singular) point impossible. Table 1. Model parameters. The parameters are equivalent to VTST00, but we changed the notation of some of them to avoid confusion with the relativistic parameters and physical quantities described in Paper I. Input parameters . . F Exponent of the radial scaling of the current Γ Polytropic index of the gas θA Angular distance of the AP from the jet axis ψA Inclination of the streamline with respect to the horizontal axis at the AP kVTST Mass-loss parameter Fitted parameters θMFP Angular distance of the MFP from the jet axis θMSP Angular distance of the MSP from the jet axis μVTST Scaling of the gas-to-magnetic pressure ratio λVTST Specific angular momentum in units of V⋆ϖ⋆ Input parameters . . F Exponent of the radial scaling of the current Γ Polytropic index of the gas θA Angular distance of the AP from the jet axis ψA Inclination of the streamline with respect to the horizontal axis at the AP kVTST Mass-loss parameter Fitted parameters θMFP Angular distance of the MFP from the jet axis θMSP Angular distance of the MSP from the jet axis μVTST Scaling of the gas-to-magnetic pressure ratio λVTST Specific angular momentum in units of V⋆ϖ⋆ Open in new tab Table 1. Model parameters. The parameters are equivalent to VTST00, but we changed the notation of some of them to avoid confusion with the relativistic parameters and physical quantities described in Paper I. Input parameters . . F Exponent of the radial scaling of the current Γ Polytropic index of the gas θA Angular distance of the AP from the jet axis ψA Inclination of the streamline with respect to the horizontal axis at the AP kVTST Mass-loss parameter Fitted parameters θMFP Angular distance of the MFP from the jet axis θMSP Angular distance of the MSP from the jet axis μVTST Scaling of the gas-to-magnetic pressure ratio λVTST Specific angular momentum in units of V⋆ϖ⋆ Input parameters . . F Exponent of the radial scaling of the current Γ Polytropic index of the gas θA Angular distance of the AP from the jet axis ψA Inclination of the streamline with respect to the horizontal axis at the AP kVTST Mass-loss parameter Fitted parameters θMFP Angular distance of the MFP from the jet axis θMSP Angular distance of the MSP from the jet axis μVTST Scaling of the gas-to-magnetic pressure ratio λVTST Specific angular momentum in units of V⋆ϖ⋆ Open in new tab 3 PARAMETER STUDY Given the wealth of solutions that we are able to retrieve using this algorithm, we focus on a grid of solutions obtained by fixing the adiabatic index Γ to 5/3, the exponent of the radial scaling of the current F to 0.75, as in Blandford & Payne (1982, hereafter BP), and the mass-loss parameter kVTST between 1.5 and 5.0 in steps of 0.5. We note that the resulting solutions will be generally different from the stereotypical BP-like solution because we include gas pressure and the crossing of all the three singular points. For each kVTST, we seek solutions with all the allowed combinations of θA and ψA, which are the angles determining the position and the collimation of the streamline at the AP, respectively. In Fig. 2, we show the distribution of these solutions in the plane of dimensionless angular momentum and entropy, i.e. the (λVTST, μVTST)-plane. Each line represents solutions for a constant kVTST and θA, while only ψA varies. Figure 2. Open in new tabDownload slide Upper panel: Grid of solutions presented in the angular momentum and entropy plane, i.e. (λVTST, μVTST)-plane for Γ = 5/3, F = 0.75. The location (θA) of the AP, collimation (ψA) at the AP, and the mass-loss parameter kVTST are allowed to vary within a chosen grid. In each kVTST = constant subset, the lines connect solutions with constant θA and variable ψA. Lower panel: Same plot for a subset of solutions for kVTST = 3.0. Neighbouring lines differ in θA by 5°. Along each line we have indicated the value of ψA for some of them to illustrate how sensitive the equations are for a small change of this angle. Particularly, a tiny change in ψA translates into a large step in μVTST, when λVTST is large. Figure 2. Open in new tabDownload slide Upper panel: Grid of solutions presented in the angular momentum and entropy plane, i.e. (λVTST, μVTST)-plane for Γ = 5/3, F = 0.75. The location (θA) of the AP, collimation (ψA) at the AP, and the mass-loss parameter kVTST are allowed to vary within a chosen grid. In each kVTST = constant subset, the lines connect solutions with constant θA and variable ψA. Lower panel: Same plot for a subset of solutions for kVTST = 3.0. Neighbouring lines differ in θA by 5°. Along each line we have indicated the value of ψA for some of them to illustrate how sensitive the equations are for a small change of this angle. Particularly, a tiny change in ψA translates into a large step in μVTST, when λVTST is large. As the upper panel of Fig. 2 shows, although our solutions cover a good extent of this region of the parameter space, a few series could not be completed because of the disappearing of the MSP below the disc mid-plane, e.g. for kVTST = 2.5 (grey line), which are physically not meaningful. In the lower panel of Fig. 2, we show how the collimation angle at the AP, ψA, is changing for a few lines on which the position of the AP, θA, is constant and the mass-loss parameter, kVTST, has been set to 3.0 for all the lines. We see that the parameters of a solution change significantly with only a small change in ψA. This is particularly true in the top part of the figure where the dimensionless angular momentum, λVTST, is large. We only find a solution when the sum of the angles θA + ψA is roughly within the interval 93–111°. This range varies depending on the value of kVTST (see Table 2). In general, the allowed range of this sum is between 90° and 180° to ensure that the derivative of the poloidal Mach number is negative, i.e. the fluid is accelerating at Alfvén. For a constant location of the AP, θA, the collimation angle ψA is small when the entropy μVTST approaches zero and is large when the angular momentum λVTST approaches zero. As we will discuss later, the combination of these two angles ultimately determines the dynamics and the geometry of the jet and the narrower range of their sum that we find is likely due to the minimum and maximum energy fluxes allowed in this region of the parameter space (see Fig. 4). Table 2. The maximum value of θA and the sum θA + ψA with a constant kVTST. The minimum of θA = 10° and of the sum θA + ψA = 93° are the same for all kVTST values. kVTST . θA,max . (θA + ψA)max . 1.5 30° 98° 2.0 45° 102° 2.5 65° 106° 3.0 65° 107° 3.5 70° 109° 4.0 75° 110° 4.5 75° 110° 5.0 80° 111° kVTST . θA,max . (θA + ψA)max . 1.5 30° 98° 2.0 45° 102° 2.5 65° 106° 3.0 65° 107° 3.5 70° 109° 4.0 75° 110° 4.5 75° 110° 5.0 80° 111° Open in new tab Table 2. The maximum value of θA and the sum θA + ψA with a constant kVTST. The minimum of θA = 10° and of the sum θA + ψA = 93° are the same for all kVTST values. kVTST . θA,max . (θA + ψA)max . 1.5 30° 98° 2.0 45° 102° 2.5 65° 106° 3.0 65° 107° 3.5 70° 109° 4.0 75° 110° 4.5 75° 110° 5.0 80° 111° kVTST . θA,max . (θA + ψA)max . 1.5 30° 98° 2.0 45° 102° 2.5 65° 106° 3.0 65° 107° 3.5 70° 109° 4.0 75° 110° 4.5 75° 110° 5.0 80° 111° Open in new tab The lowest value of the sum, i.e. 93 (small θA, large ψA), coincides with the jet configurations with lowest total energy-to-mass flux ratios which is around |$\sim V_{\rm A,p}^2/2$| at the jet base (z = 0) for kVTST of the order of unity (equation 14). These solutions have little-to-none magnetic field (λVTST → 0), and represent a tenuous jet (low total energy flux, see Sections 3.1 and 3.2) supported by some (μVTST small) gas pressure, which provides the balance to gravity. By varying the two angles within the allowed range, we recover a large collection of solutions where we see low-energy hot jets transform into cold and fast jets with a large angular momentum (λVTST ≳ 20) and a relatively small contribution of the gas pressure (small μVTST) to the total energy. The large variety of physical properties within this sample of solutions provides an ideal framework to study the different jet configurations and to devise a method for the comparison of such solutions to astrophysical sources. 3.1 General trends In this section, we discuss some general properties and trends observed while inspecting the whole ensemble of solutions. In Fig. 3, we show how the total energy is divided up between rotational energy and generalized pressure (Ferreira 1997). The rotational energy is the difference between the total energy in an inertial frame and the total energy in a frame rotating with a frequency Ω (equation 16), i.e. $$\begin{eqnarray*} E_{\rm rot} = L\Omega = \varpi _{\rm A}^2\Omega ^2. \end{eqnarray*}$$(20) In the above equation, L is the angular momentum, defined as $$\begin{eqnarray*} L = \varpi \left(V_\phi - \frac{B_\phi }{4\pi \rho }\frac{B_p}{V_p} \right), \end{eqnarray*}$$(21) which is also a constant of motion along the streamline. Fig. 3 shows the total energy in the rotating frame rescaled by the poloidal kinetic energy at the AP (blue dots): $$\begin{eqnarray*} 2\frac{E - L\Omega }{V_{\rm A,p}^2} = \tilde{\epsilon }- \tilde{\epsilon }_{\rm rot} \end{eqnarray*}$$(22) and the rotational energy rescaled, |$2L\Omega /V_{\rm A,p}^2= \tilde{\epsilon }_{\rm rot}$| (yellow dots) versus the rescaled total energy in the inertial frame |$\tilde{\epsilon }= 2E/V_{{\rm A},p}^2$|⁠. Figure 3. Open in new tabDownload slide The rotational energy as a function of the total energy scaled with the poloidal velocity at Alfvén for all solutions (yellow dots). In blue is shown the total energy in a frame rotating with the constant angular frequency Ω versus the total energy in the inertial frame, also scaled with the poloidal velocity at Alfvén. The red crosses are solutions with kVTST = 3.0 and θA = 60°. They will be used to discuss other trends later on. Figure 3. Open in new tabDownload slide The rotational energy as a function of the total energy scaled with the poloidal velocity at Alfvén for all solutions (yellow dots). In blue is shown the total energy in a frame rotating with the constant angular frequency Ω versus the total energy in the inertial frame, also scaled with the poloidal velocity at Alfvén. The red crosses are solutions with kVTST = 3.0 and θA = 60°. They will be used to discuss other trends later on. All the points lie on two narrow curves. The solutions highlighted in the bottom panel of Fig. 2 are marked as red crosses in Fig. 3. The total energy in the rotational frame, |$\tilde{\epsilon }{\!-\!}\tilde{\epsilon }_{\rm rot}$| (blue dots), otherwise called the generalized pressure (Pelletier & Pudritz 1992; Ferreira 1997), achieves a maximum when. the rotational energy, |$\tilde{\epsilon }_{\rm rot}$| (yellow dots), is negligible. Since the total energy flux sustaining a jet, i.e. the Bernoulli constant, is positive, the generalized pressure can change sign depending on the relative contribution of the rotational energy, |$\tilde{\epsilon }_{\rm rot}$|⁠, to the total energy. As the rotational energy, |$\tilde{\epsilon }_{\rm rot}$|⁠, increases, it approaches equipartition with the generalized pressure which occurs in the regime where the latter is still positive. When the sign flip occurs, we start to see a dominant contribution of the magnetic energy in the Bernoulli equation (equation 12). Based on the ratio between the thermal energy and the magnetic energy flux we distinguish three categories of solutions: thermally dominated hot, equipartition/centrifugal and magnetically dominated cold jets (see Fig. 4). We show this ratio at the disc mid-plane (blue squares) and at the MSP (pink crosses) for the full sample of solutions versus the total energy rescaled with the poloidal kinetic energy. The vertical lines are drawn to guide the eye. We see that the distribution of the jet models in this plane is very similar between z = 0 and the MSP. The hot jets are low-energy solutions and as the ETH/EM increases the total energy flux, |$\tilde{\epsilon }$|⁠, remains constant and at its minimum value. When the thermal and magnetic energy flux are roughly at equipartition, the total energy is increasing steadily as the solutions become more magnetically dominated. As we enter the cold regime, the magnetic energy grows more rapidly for a small variation of the input parameters (see bottom panel of Fig. 2), but the jet configurations do not increase so much in total energy anymore, approaching its maximum. Since there is this correspondence between total energy and hot/cold regime, we will use it interchangeably across the paper. Figure 4. Open in new tabDownload slide Ratio between the thermal energy and magnetic energy as a function of the total energy with respect to the poloidal kinetic energy at Alfvén for all solutions. Figure 4. Open in new tabDownload slide Ratio between the thermal energy and magnetic energy as a function of the total energy with respect to the poloidal kinetic energy at Alfvén for all solutions. In Fig. 5, we show the different contributions to the total energy at the base (top panel) and at the MSP (bottom panel) for a series of solutions with kVTST = 3.0 and θA = 60°. The trends discussed here are also observed in other series. We start by noticing that when the magnetic energy is larger, the total energy is larger too. When the total energy is low, the gravitational energy and the thermal energy dominate with almost equal magnitude, cancelling each other. Only at higher total energies the thermal energy becomes negligible. Apart from the most energetic solutions, the kinetic energy consists mainly of the poloidal component. At higher energies the poloidal component of the velocity of the gas leaving the mid-plane is relatively low while the toroidal speed gives the largest contribution to the total kinetic energy. Figure 5. Open in new tabDownload slide The total energy at the base (top panel) and at the MSP (bottom panel) split up in its components for a series of solutions for kVTST = 3.0 and θA = 60°. All energies have been scaled with the poloidal kinetic energy at Alfvén. Figure 5. Open in new tabDownload slide The total energy at the base (top panel) and at the MSP (bottom panel) split up in its components for a series of solutions for kVTST = 3.0 and θA = 60°. All energies have been scaled with the poloidal kinetic energy at Alfvén. In Fig. 6, we present the components of the velocity and of the angular frequency Ω of the streamlines (equation 16) for the same series of solutions presented in Fig. 5. At lower total energies the poloidal velocity is relatively large with respect to the toroidal velocity. As a consequence, even when the ratio of the magnetic field components (grey line with dots) is low, i.e. the magnetic field is almost not twisted at all, the second term on the RHS of the equation describing Ω (equation 16, magenta line with stars) is dominant, while the toroidal velocity (brown line with pentagons) is negative and smaller. This means that the gas is lagging behind the rotation of the disc and the magnetic field is weak, while Ω is at its minimum. Figure 6. Open in new tabDownload slide Velocities at the base (top panel) and at the MSP (bottom panel) scaled by the poloidal velocity at Alfvén as a function of the scaled total energy for a series of solutions for kVTST = 3.0 and θA = 60°. Figure 6. Open in new tabDownload slide Velocities at the base (top panel) and at the MSP (bottom panel) scaled by the poloidal velocity at Alfvén as a function of the scaled total energy for a series of solutions for kVTST = 3.0 and θA = 60°. Only the very last solution with the highest energy of this series is rotating at Keplerian speed, which can be seen by noticing that the last dot of the pink line with crosses (Vk) coincides with the last point of the brown line with pentagons (Vϕ) in the top panel of Fig. 6. As expected by the non-negligible contribution of the enthalpy, the overwhelming majority of the solutions in this ensemble is sub-Keplerian at the disc mid-plane, with a deviation increasingly larger as the solutions become warmer and warmer. This means that the typical approximation Ω ∼ Ωgas = Vϕ/ϖ ∼ Ωk cannot be taken as a general property of this sample of solutions. Only a small fraction of the solutions presented in this paper can be considered corotating with the disc, like for instance the last three high-energy solutions in Fig. 6, where we see that Ωϖ (green line with crosses) matches Vϕ (brown line with pentagons), while −VpBϕ/Bp (magenta line with stars) is close to zero. From a geometrical point of view, the radial profile of the streamlines varies depending on how hot the jet is, typically with highly oscillating jet bases for cold jets while no oscillations are present for warm and hot jets (see Fig. 7). This is a consequence of the oscillatory nature of the transverse component of the forces that define the collimation of the streamline. We will discuss this topic in detail in Section 3.2. Since it is very likely that such oscillations may be unstable and considering that the MSP is a more robust point in our solutions, we identify the MSP with the jet base from now on. Figure 7. Open in new tabDownload slide Examples of streamlines of a series of solutions for increasing μVTST. Figure 7. Open in new tabDownload slide Examples of streamlines of a series of solutions for increasing μVTST. Different jet configurations can be also classified based on the amount of acceleration that the gas experiences from the MSP to the MFP, being that the point where the flow loses causal contact with the source and the flow upstream. In Fig. 8, we plot all the solutions divided in subgroups with constant kVTST in the plane defined by the increase in the poloidal velocity experienced by the matter from the MSP to the MFP and the rescaled total energy flux. The low-energy flux, pressure-driven, solutions have also low ΔV/VMFP since they are characterized by large poloidal velocities at the MSP which do not increases much approaching the MFP. As the energy flux increases, the poloidal velocity decreases (see bottom panel of Fig. 6) and the increment of the velocity ΔV/Vp,MFP approaches 1. Figure 8. Open in new tabDownload slide Acceleration fraction versus the total energy rescaled. Figure 8. Open in new tabDownload slide Acceleration fraction versus the total energy rescaled. Similarly, the acceleration of the flow is also traced by the increase in the poloidal kinetic energy. In Fig. 9, we show the relative increment/decrement of the energy fluxes between the MSP and the AP (first phase of the acceleration, top three panels) and the AP and the MFP (second phase of acceleration, bottom three panels) in a transition from hot to cold solutions (low- to high-energy flux). In the first phase of the acceleration, hot solutions are driven by the thermal energy which suffers the largest decrement. However, as highlighted by the zoom around zero, a fraction of the thermal energy is transferred to the magnetic energy, which is increasing for hot solutions with energy fluxes <5. This behaviour is followed closely by the relative increment/decrement of the components of the angular momentum. For these hot solutions the hydrodynamical component of L decreases, while the magnetic component increases, showing that the angular momentum of the gas is transferred to the angular momentum of the magnetic field. Such additional channel of energy transfer has been seen in simulations such as e.g. Komissarov et al. (2009), Cayatte et al. (2014), and in Paper I. This effect is seen as well in the bottom panel of Fig. 13 as a small rise in the magnetic energy around the AP. As the jet models move to higher energy configurations, the magnetic energy increases while the thermal energy is still important, leading to an increasing poloidal kinetic energy. The peak of the poloidal kinetic energy occurs in correspondence to ΔEM/Etot ≃ ΔETH/Etot. Then, it decreases again due to a decrease in toroidal kinetic energy. In the second phase of the acceleration, the thermal energy still dominates for hot low-energy solutions. Equipartition/MC and cold solutions instead are accelerated all the way from the MSP to the MFP by the magnetic field. In the upper part of the jet, the relative increment of the components of the angular momentum do not change sign and the magnetic angular momentum is always transferred to the gas component. Figure 9. Open in new tabDownload slide Upper panels: Relative increment of the energy fluxes with respect to the total energy flux, |$(E_{\cdot ,\rm AP} - E_{\cdot ,\rm MSP})/E_{\rm tot}$|⁠, versus the total energy rescaled from the MSP to the AP. The small left panel is a zoom around zero. The small right panel shows the relative increment of the components of the angular momentum with respect to the total angular momentum between the MSP and the AP versus the total energy rescaled. The solutions are obtained for θA = 65° and kVTST = 4.0. Lower panels: Same as the upper panel but between the AP and the MFP. Figure 9. Open in new tabDownload slide Upper panels: Relative increment of the energy fluxes with respect to the total energy flux, |$(E_{\cdot ,\rm AP} - E_{\cdot ,\rm MSP})/E_{\rm tot}$|⁠, versus the total energy rescaled from the MSP to the AP. The small left panel is a zoom around zero. The small right panel shows the relative increment of the components of the angular momentum with respect to the total angular momentum between the MSP and the AP versus the total energy rescaled. The solutions are obtained for θA = 65° and kVTST = 4.0. Lower panels: Same as the upper panel but between the AP and the MFP. Moreover, since the downstream portion of the MFP might be already affected by a shock given by the loss of causal contact with the flow upstream, we take as a proxy the total jet length the distance between the MSP and the MFP. The top panel of Fig. 10 shows that low-energy solutions can be as short as 102/ϖ* and as long as 106/ϖ*. As the total energy increases this interval narrows by ∼2 order of magnitude (103−7 × 104). We note that when the streamlines become more vertical (increasing ψA), this leads to a decrease in total energy (the lines in the plot are drawn for constant angular position of the AP, θA), while increasing θA (from top to bottom) makes the Δz decrease. If we were to focus on one of the most extended lines across the energy range, for instance the red dashed line that is for an intermediate constant value of the angular position of the AP (θA = 50°) and a fixed mass-loss parameter (kVTST = 4.0), we would see a correlation between the distance between the MSP and the MFP and the total energy: the higher the energy the larger is the distance, until it reaches an almost constant length (∼20 000/ϖ*). We note that the maximum of Δz does not coincide though with the highest energy in the line. Therefore, beyond a certain total energy, the jets do not grow taller, but their ΔV/Vp,MSP increases as shown in the middle panel of Fig. 10 and in Fig. 8. Low-energy hot solutions increase in length by a factor of 10 as the collimation angle, ψA, increases, maintaining their velocity increment roughly constant (bottom panel of Fig. 10). Figure 10. Open in new tabDownload slide Upper panel: Distance between the MFP and MSP versus total energy rescaled |$\tilde{\epsilon }= 2 E_{\rm tot}/V_{p,\rm A}^2$|⁠. The lines are connecting solutions with constant angular position of the AP (θA) and the colours mark different values of the mass-loss parameter kVTST. The red dashed line is a series of solutions obtained with θA = 50° and kVTST = 4.0. Middle panel: Poloidal velocity scaled with the Alfvén poloidal velocity (⁠|$\tilde{V}_p = V_p/V_{p,\rm A}$|⁠) along the streamline for a set of solutions with θA = 50° and kVTST = 4.0. The solutions in this plot have a total energy flux rescaled (⁠|$\tilde{\epsilon }$|⁠) going from 1.5 to 10.5 and distances between the MSP and MFP changing by a factor of ∼50. Bottom panel: Same plot as the middle one, but for θA = 30°. The solutions in this plot have roughly a constant total energy flux rescaled (⁠|$\tilde{\epsilon }\sim 1.3$|⁠) and distances between the MSP and MFP changing by one order of magnitude, going from 7900 to 72 500. The second solutions in the middle and bottom panel have the largest Δz, but not the largest |$\tilde{\epsilon }$|⁠. Figure 10. Open in new tabDownload slide Upper panel: Distance between the MFP and MSP versus total energy rescaled |$\tilde{\epsilon }= 2 E_{\rm tot}/V_{p,\rm A}^2$|⁠. The lines are connecting solutions with constant angular position of the AP (θA) and the colours mark different values of the mass-loss parameter kVTST. The red dashed line is a series of solutions obtained with θA = 50° and kVTST = 4.0. Middle panel: Poloidal velocity scaled with the Alfvén poloidal velocity (⁠|$\tilde{V}_p = V_p/V_{p,\rm A}$|⁠) along the streamline for a set of solutions with θA = 50° and kVTST = 4.0. The solutions in this plot have a total energy flux rescaled (⁠|$\tilde{\epsilon }$|⁠) going from 1.5 to 10.5 and distances between the MSP and MFP changing by a factor of ∼50. Bottom panel: Same plot as the middle one, but for θA = 30°. The solutions in this plot have roughly a constant total energy flux rescaled (⁠|$\tilde{\epsilon }\sim 1.3$|⁠) and distances between the MSP and MFP changing by one order of magnitude, going from 7900 to 72 500. The second solutions in the middle and bottom panel have the largest Δz, but not the largest |$\tilde{\epsilon }$|⁠. Lastly, we discuss the variation of the plasma-β and mass load η at the MSP, which we identify with the jet base as discussed above. These two quantities are given by $$\begin{eqnarray*} \beta = \frac{P}{B^2/8\pi } \qquad {\rm and} \qquad \eta = \frac{4\pi \rho V_p\varpi \Omega }{B_p^2}, \end{eqnarray*}$$(23) following the definitions of e.g. Anderson et al. (2005) and Spruit (1996). Since the general trend is the same within subsets of solutions with constant kVTST, we present here the series of solutions obtained for kVTST = 4.5 and for Alfvén position angle (PA), θA, going from 10° to 75° roughly from the bottom up (Fig. 11) and then discuss how they change for increasing kVTST at constant θA (Fig. 12). Solutions in both figures have increasing collimation angle, ψA, along each line from left to right (from 28° to 83° in Fig. 11 and from 48° to 57° in Fig. 12). In Fig. 11, we see that thermally dominated, low-energy flux solutions have the largest plasma-β (∼1). Then as the collimation angle, ψA, decreases, the energy flux increases and the plasma-β experiences a first decrease. For the Alfvén angular positions for which more ψA values are allowed, we see the plasma-β remaining constant for many consecutive solutions of increasing energy flux. However, when the solutions become magnetically dominated, the plasma-β has a drop. The mass load has a minimum which coincides with the beginning of the plateaux of the plasma-β, to again rise to higher energy fluxes. The relatively large mass load of the low-energy flux solutions is due to high density of the gas, while a similar value is reached for the high-energy flux solutions because the magnetic field is more tightly wound up (|Bϕ/Bp| > 1, see e.g. Spruit 1996; Anderson et al. 2005). Figure 11. Open in new tabDownload slide Plasma-β (black lines with dots, top panel) and mass load η (red lines with crosses, bottom panel) versus total energy flux rescaled for all the solutions with kVTST = 4.5 at the MSP. Each line connects solutions with constant θA and increasing ψA. The arrows show the approximate direction of increasing θA and ψA. Figure 11. Open in new tabDownload slide Plasma-β (black lines with dots, top panel) and mass load η (red lines with crosses, bottom panel) versus total energy flux rescaled for all the solutions with kVTST = 4.5 at the MSP. Each line connects solutions with constant θA and increasing ψA. The arrows show the approximate direction of increasing θA and ψA. Figure 12. Open in new tabDownload slide Plasma-β (lines with dots, top panel) and mass load η (lines with crosses, bottom panel) versus total energy flux rescaled for all the solutions with θA = 45° and kVTST = 3.0 (blue), 3.5 (pink), 4.0 (grey), 4.5 (green), 5.0 (yellow) at the MSP. Each line connects solutions with constant θA = 45° and varying ψA. Figure 12. Open in new tabDownload slide Plasma-β (lines with dots, top panel) and mass load η (lines with crosses, bottom panel) versus total energy flux rescaled for all the solutions with θA = 45° and kVTST = 3.0 (blue), 3.5 (pink), 4.0 (grey), 4.5 (green), 5.0 (yellow) at the MSP. Each line connects solutions with constant θA = 45° and varying ψA. In Fig. 12, we show how the same quantities vary in relation to an increase in kVTST. The plasma-β and the mass load, η, show a similar behaviour with respect to the mass-loss parameter, kVTST: the larger is kVTST the larger is the plasma-β and the mass load. However, we notice that η has a weaker dependence on kVTST both at low- and high-energy fluxes, while the plasma-β responds to a change in kVTST more homogeneously across the energy flux interval. 3.2 Hot and cold jets To illustrate the qualitative changes of the outflow properties along a series of solutions for increasing collimation angle ψA, we describe the transition looking at the two extreme solutions plots of the components of the Bernoulli equation (equation 12) and an intermediate one which resembles a more classical magneto-centrifugally launched jet. We will refer to these solutions as cold, magneto-centrifugal (MC), and hot jet models and list their parameters in Table 3. As shown in Fig. 13, the energy fluxes along the poloidal direction are substantially different going from the cold (upper panel) to the hot (lower panel) jet solution. The cold jet has a high Poynting-to-enthalpy flux ratio. The magnetic energy is then converted into kinetic energy downstream of the AP. Upstream of the MSP, all the energy fluxes are oscillating, following the oscillations of the radial profile of the streamline (see Fig. 7). The intermediate MC jet solution has qualitatively the same characteristics of the cold one, but the oscillations are gone. The hot jet has an uneventful behaviour of the energy fluxes along the streamline. The enthalpy is dominant and roughly equal to gravity in absolute value and opposite in sign. Right after the AP, initially the thermal energy flux is the main source of energy being transformed into kinetic energy and into magnetic energy, which shows a small increase, as discussed in Section 3.1. Then, the magnetic energy flux takes over the final acceleration. For constant mass-loss parameter, kVTST, the total energy flux is ∼2 orders of magnitude larger for the cold jet. This larger energy reservoir allows the cold jet to extend in length a factor of ∼100 more than the hot jet, when the same reference scale length, ϖ* is applied. Figure 13. Open in new tabDownload slide Energy fluxes along the streamline. The green line is Poynting flux, the pink line is the enthalpy energy flux, the brown line is the kinetic energy flux, and the dashed purple line is the gravitational energy with reversed sign. The total energy flux is shown as a solid black line. Note that the y-axis scale is different in the three plots, while the x-axis scale is the same. From top to bottom the solutions go from cold to hot. The parameters are listed in Table 3. Figure 13. Open in new tabDownload slide Energy fluxes along the streamline. The green line is Poynting flux, the pink line is the enthalpy energy flux, the brown line is the kinetic energy flux, and the dashed purple line is the gravitational energy with reversed sign. The total energy flux is shown as a solid black line. Note that the y-axis scale is different in the three plots, while the x-axis scale is the same. From top to bottom the solutions go from cold to hot. The parameters are listed in Table 3. Table 3. Parameters of the solutions used in Section 3.2. The solutions have the following common parameters kVTST = 3.0, θA = 60°, Γ = 5/3, and F = 0.75. The cold jet and the hot jet models are the extremes of the series, while the MC jet is an intermediate one which is closest to the classical magneto-centrifugal jets encountered in the literature. Model . μVTST . λVTST . θMFP . θMSP . ψA . Cold jet 7.58 × 10−2 17.3853 0.11803 1.2462 37.07 MC jet 0.5396 16.8723 0.11778 1.2578 37.09 Hot jet 6.5510 1.6861 0.12347 1.3852 46.00 Model . μVTST . λVTST . θMFP . θMSP . ψA . Cold jet 7.58 × 10−2 17.3853 0.11803 1.2462 37.07 MC jet 0.5396 16.8723 0.11778 1.2578 37.09 Hot jet 6.5510 1.6861 0.12347 1.3852 46.00 Open in new tab Table 3. Parameters of the solutions used in Section 3.2. The solutions have the following common parameters kVTST = 3.0, θA = 60°, Γ = 5/3, and F = 0.75. The cold jet and the hot jet models are the extremes of the series, while the MC jet is an intermediate one which is closest to the classical magneto-centrifugal jets encountered in the literature. Model . μVTST . λVTST . θMFP . θMSP . ψA . Cold jet 7.58 × 10−2 17.3853 0.11803 1.2462 37.07 MC jet 0.5396 16.8723 0.11778 1.2578 37.09 Hot jet 6.5510 1.6861 0.12347 1.3852 46.00 Model . μVTST . λVTST . θMFP . θMSP . ψA . Cold jet 7.58 × 10−2 17.3853 0.11803 1.2462 37.07 MC jet 0.5396 16.8723 0.11778 1.2578 37.09 Hot jet 6.5510 1.6861 0.12347 1.3852 46.00 Open in new tab The forces acting along (⁠|$\hat{b}$|⁠) and perpendicular (⁠|$\hat{n}$|⁠) to the streamline highlight the transition from cold to hot jet configurations. Here, we give the compact form of the forces in both direction, while we provide the full derivation in Appendix C. $$\begin{eqnarray*} \hat{b}: \frac{\rho }{2} \frac{\partial V_p^2}{\partial l} &=& \rho V_\phi ^2 \frac{\cos (\psi)}{\varpi } - \frac{\partial P}{\partial l} + \rho \frac{\partial }{\partial l}\left(\frac{\mathcal {GM}}{r}\right) \nonumber \\ &&- \frac{1}{8\pi }\frac{\partial B_\phi ^2}{\partial l} - B_\phi ^2\frac{\cos (\psi)}{4\pi \varpi } , \end{eqnarray*}$$(24) $$\begin{eqnarray*} \hat{n}: \left(\rho V_p^2 - \frac{B_p^2}{4\pi }\right)\frac{\partial \psi }{\partial l} &=& + \rho V_\phi ^2 \frac{\sin (\psi)}{\varpi } - \frac{\partial P}{\partial n} + \rho \frac{\partial }{\partial n}\left(\frac{\mathcal {GM}}{r}\right) \nonumber \\ && - \frac{1}{8\pi }\frac{\partial }{\partial n}\left(B_p^2 + B_\phi ^2\right) + B_\phi ^2\frac{\sin (\psi)}{4\pi \varpi }. \end{eqnarray*}$$(25) The term on the LHS of the equation (24) is the acceleration along the streamline, the first term on the RHS is the centrifugal force, the second term is the gas pressure force, the third term is the gravitational force, and the last two terms are the magnetic pressure gradient and the magnetic tension. On the LHS of equation (25) there is the derivative of the angle ψ along the streamline. The inverse of this derivative is also called the collimation radius, Rc = (∂ψ/∂l)−1. On the RHS there are the centrifugal force, the gas pressure force, the gravitational force, and the magnetic pressure gradient and the magnetic tension. In the following discussion, we refer to accelerating/collimating forces when such terms are positive, and to decelerating/decollimating forces when they are negative. In Fig. 14, we show the forces perpendicular to the streamline and in Fig. 15 the forces along the streamline for the same three solutions. Figure 14. Open in new tabDownload slide Forces perpendicular to the streamline. The black line shows the derivative of the collimation angle along the streamline. The lines are solid when the force is providing collimation (positive) and dashed when it is instead decollimating (negative). From top to bottom the solutions go from cold to hot. The parameters are listed in Table 3. Figure 14. Open in new tabDownload slide Forces perpendicular to the streamline. The black line shows the derivative of the collimation angle along the streamline. The lines are solid when the force is providing collimation (positive) and dashed when it is instead decollimating (negative). From top to bottom the solutions go from cold to hot. The parameters are listed in Table 3. Figure 15. Open in new tabDownload slide Forces along the streamline. The black line is the poloidal acceleration along the streamline. The lines are solid when the force is providing collimation (positive) and dashed when it is instead decollimating (negative). From top to bottom the solutions go from cold to magneto-centrifugal to hot. The parameters are listed in Table 3. Figure 15. Open in new tabDownload slide Forces along the streamline. The black line is the poloidal acceleration along the streamline. The lines are solid when the force is providing collimation (positive) and dashed when it is instead decollimating (negative). From top to bottom the solutions go from cold to magneto-centrifugal to hot. The parameters are listed in Table 3. The cold jet has a troublesome start, since it lacks a vertical velocity component that allows for a straightforward launching (top panel in Fig. 16). At the very beginning, the jet is decollimating (∂ψ/∂l < 0, black thin line) under the action of the gas pressure force (pink line). Soon, the gas pressure gradient changes sign and together with the other positive forces, i.e. gravity (purple line), centrifugal (brown line), and magnetic tension (teal line), is collimating the jet against magnetic pressure gradients. Around the peak of gravity and the centrifugal force, the pressure gradient becomes negative but smaller in modulus, resulting in a converging streamline (thick solid black line). After that the previous configuration of the forces is mirrored to the right-hand side of the peak, until shortly before the MSP, the streamline starts to decollimate again. However, downstream of the MSP the sign switches again when the magnetic tension becomes dominant, keeping the jet collimated up until also the magnetic pressure gradient becomes positive, about half way between the AP and the MFP (top panel of Fig. 14). The MC jet model shows the same behaviour downstream of the MSP, while it presents no oscillations in the region between the disc and the MSP. Figure 16. Open in new tabDownload slide Zoom of the forces along the streamline (bottom) and perpendicular (top) to it for the cold jet model (Table 3). The thick solid black line is the streamline. Figure 16. Open in new tabDownload slide Zoom of the forces along the streamline (bottom) and perpendicular (top) to it for the cold jet model (Table 3). The thick solid black line is the streamline. The hot jet is always collimating. Until the AP, gravity is the main force driving the collimation against the gas pressure gradient that remains negative until past the AP. Beyond this point, the magnetic forces become dominant in keeping the jet focused. The poloidal forces also presents oscillations in the cold jet model, while they do not in the MC and hot jet models. In the top panel of Fig. 15, we see some more moderate oscillations for the initial segment of the cold jet. In the upstream region of the MSP, the jet is initially slowly accelerating (⁠|$\partial V_p^2/\partial l\gt 0$|⁠). Then the pressure force (pink line) becomes negative and gravity (purple line) is attracting the fluid back to the centre (the thick solid black line in Fig, 16 shows the streamline focussing towards the axis), now increasing its speed, providing acceleration while the poloidal motion of the flow (think black line) is actually decelerating (bottom panel in Fig. 16). At the minimum radius of the streamline, all the forces change sign and the jet starts to accelerate driven by a combination of centrifugal (brown) and magnetic force (teal line). Half way between the MSP and the AP, the magnetic force takes over and it will sustain the acceleration for the remaining (and larger) fraction of the jet extent. In the intermediate jet solution the flow is accelerated by the gas pressure force (and for a small segment just downstream of the MSP by the centrifugal and the magnetic forces) until halfway the MSP and the AP, when the magnetic force drives again the acceleration of the jet until the LRP. The hot jet is instead decelerating until past the MSP, then the pressure provides acceleration working against the gravitational pull. Finally, halfway between the AP and the MFP, the magnetic force becomes the dominant accelerating force for the rest of the jet length. 4 PROOF OF CONCEPT: APPLICATION TO W43A In this section, we describe how to compare our solutions to an astrophysical source, the water fountain W43A. W43A is a pre-planetary nebula (PPN; plural, PPNe), located at a distance of 2.2 kpc from the sun (Tafoya et al. 2020), that is thought to be hosting an asymptotic giant branch (AGB) star (Imai & Diamond 2005; Tafoya et al. 2020). It has been observed that during the transition from the AGB to planetary nebula (PN) phases, the star’s ejecta change from a roughly spherical symmetric wind to an envelope with a highly non-spherical configuration (Balick & Frank 2002). These non-spherical post-AGB or PPNe envelopes often exhibit (collimated) bipolar outflows and/or jets, which are most likely formed at the time the star leaves the AGB (e.g. Sahai & Trauger 1998). The origin of the non-spherical outflows around W43A and other PPNe is a matter of debate, and is typically thought to include a common envelope evolution phase (e.g. Nordhaus & Blackman 2006). It is suggested that W43A also hosts a close companion embedded in the circumstellar envelope of the AGB star, likely a main-sequence star or a white dwarf, although such companion has not been directly observed (e.g. Imai et al. 2002; Imai & Diamond 2005; Tafoya et al. 2020). The binary interaction between the two stars is expected to lead to the ejection of the envelope. During this phase, both a circumbinary disc and an accretion disc around the companion can form. It has been proposed that fast outflows, either collimated or wide, can be launched before, during and/or after the common envelope phase, contributing to the evolution of the system by heating and mechanically re-disturbing the material of the envelope, possibly leading to its ejection (Chamandy et al. 2018; Soker 2020). A scenario in which jets are launched at the onset of the short-lived water fountain phase of W43A life cycle seems plausible considering the current properties of the source (Tafoya et al. 2020). Following the argument that Sahai et al. (2017) used for the water fountain IRAS 16342−3814, if we were to assume that the radiation pressure is the main force responsible for the launching and acceleration the jets of W43A, we could estimate the time-scale for ejecting such radiation-driven jets as $$\begin{eqnarray*} \tau _{\rm rad} = \frac{Pc}{L} , \end{eqnarray*}$$(26) where P is the total momentum, L is the luminosity of the source, and c is the speed of light. The momentum derived from observational constraints is P ∼ 3.06 × 1037 g cm s−1. Adopting a luminosity of 6000 L⊙ given by Duran-Rojas et al. (2014), we obtain a time-scale of ∼1268 yr which is almost 20 times larger than the dynamical time-scale (tdyn ∼ 65 yr) estimated by Tafoya et al. (2020). Thus, radiation can be ruled out as the mechanism responsible for launching and accelerating the jet. Several mechanisms have been proposed to produce collimated jets, many of which make use of magnetic fields to drive, or at least to strongly contribute to, the acceleration and collimation of the material from a rotating object, i.e. a star, a compact object, or a disc (e.g. Blandford & Znajek 1977; Blandford & Payne 1982; Ferreira 1997; Shu et al. 2000; Parfrey et al. 2016) and a similar contribution has been proposed for PPNe as well (e.g. García-Segura, López & Franco 2005). 4.1 Observational constraints Recent observations by Tafoya et al. (2020) show that W43A possesses a dense (n ∼ 2 × 107 cm−3), collimated (z/ϖ ∼ 20, where ϖ is the radius of the jet and z is its height) molecular jet. The molecular jet inclination angle with respect to the plane of the sky is 35°, and its PA (with respect to the north) is 68°. The jet extends with constant collimation angle out to a distance from the central source of ≈1600 au, and it is surrounded by two lobes of shocked material with a lower density (n ∼ 3 × 106 cm−3) (see Fig. 17). Figure 17. Open in new tabDownload slide Sketch of W43A masers emission regions. Figure 17. Open in new tabDownload slide Sketch of W43A masers emission regions. W43A is known to host maser emission from different chemical species, such as OH, H2O, and SiO. The OH masers are located on an expanding torus of radius ∼500 au with an expansion velocity of ∼18 km s−1 and a velocity separation of ∼16 km s−1. The density required for the excitation of the OH masers at that distance is ∼104−106 cm−3 (Elitzur, Hollenbach & McKee 1992). The H2O maser emission is observed at the two regions where the jet seems to be interacting with the lobes. The H2O maser spots have velocities ∼150 km s−1 and hydrogen densities ∼108−1010 cm−3 (Imai et al. 2002; Vlemmings & Diamond 2006; Vlemmings, Diamond & Imai 2006). SiO masers have also been observed, at ∼70 au from the star Imai & Diamond (2005), and were modelled as an expanding shell of shocked material surrounding a high-velocity outflow. The magnetic field in the material surrounding W43A has been measured using observations of the Zeeman splitting of H2O and OH masers (Vlemmings et al. 2006; Amiri, Vlemmings & van Langevelde 2010). The magnetic field strength measured in the H2O maser regions is ∼200 mG (Vlemmings et al. 2006). The magnetic field in the maser regions is likely enhanced, due to compression of the field lines in the shocked interaction region between the jet and the surrounding medium. The H2 number density in the lobes around the jet is estimated to be 3 × 106 cm−3 and that in the surrounding shell is 5 × 108 cm−3 (Tafoya et al. 2020, Fig. 17). Using these densities to update the uncompressed magnetic field estimates from Vlemmings et al. (2006) and Amiri et al. (2010) and assuming a typical H2O maser region number density of 109 cm−3 and a magnetic C-shock, we find a magnetic field strength in the range of ∼0.6 mG, when the shock occurs in the lower density material of the lobes, to ∼100 mG, if the shock occurs in the denser shell surrounding the lobes. Since the exact maser density is unknown, the uncertainty on these values is large. Although it is unclear exactly which component of the magnetic field is traced by the H2O maser measurements, the linear polarization direction and evidence of change in sign of the measured magnetic field across the jet indicate that the masers likely probe the toroidal magnetic field (Amiri et al. 2010). We refer to the bipolar high-velocity outflow traced by the H2O masers as the molecular jet of W43A. 4.2 Modelling assumptions How such molecular jet is launched and how it maintains its collimation throughout its length has not yet been established. We hypothesize that what is shaping the molecular jet of W43A is a disc-driven MHD jet. More specifically, we assume that an MHD jet is launched by an accretion disc formed around a white dwarf companion (⁠|$\mathcal {M} = 0.6\, \mathrm{ M}_\odot$|⁠, |$\mathcal {R} = 0.01\, \mathrm{ R}_\odot$|⁠) orbiting around the AGB star (see Fig. 18). In this scenario, the MHD jet is accelerated outwards and entrains material from the surroundings, building up a more mass loaded, slower cocoon which is observed as a molecular jet (e.g. Hardee 1996; Rosen et al. 1999). In this paper, we show that the properties of disc-driven MHD jets can be very diverse. In order to reduce the allowed range of such properties, we compare our solutions to the observational constraints of the observed molecular jet of W43A. If the MHD jet is driving the molecular jet, its momentum has to be at least equal, or larger, than the momentum carried by the molecules (⁠|$P_j \ge P_{j,\rm obs}$|⁠). Since the way in which the composition of the MHD jet relates to the molecular content is unknown, we assume that its hydrogen number density is at most equal to the hydrogen density estimated from the CO mass (n ≥ nobs). With a lower density, the MHD jet is also likely to travel at a faster speed than the water maser spots in that region (⁠|$V_{\rm H_{2}O} \ge V_{\rm H_{2}O, obs}$|⁠). Finally, we adopt the full range of the toroidal magnetic field strength (Bϕ = 0.6−100 mG) to look for solutions which have the requirements listed above and extend for 2000 au. Given these uncertainties, we choose to compare the momentum rate carried by our solutions with the one estimated with the observational constraints to avoid a direct comparison between densities and velocities where we should instead make more assumptions such as on the ionization fraction or on the intrinsic speed of the jet. This allow us to predict the general properties of the MHD jet to drive the molecular jet sheath surrounding it. Figure 18. Open in new tabDownload slide Sketch of the central region of W43A. Figure 18. Open in new tabDownload slide Sketch of the central region of W43A. To estimate the momentum rate from observations, we approximate the molecular jet of W43A as a full cylinder with a length of 2000 au and a radius of 45 au. The momentum rate can be estimated as |$\dot{P} = M_j V_{\rm H_{2}O}/t_{\rm dyn}$|⁠, which for a jet with a total mass of ∼10−3 M⊙, a velocity of 150 km s−1 and a dynamical time-scale of 65 yr is ∼0.0023 M⊙ yr−1 km s−1. This is equivalent of a total momentum over 65 yr of P ∼ 3.06 × 1037 g cm s−1. We note that this estimated value of the momentum of W43A lies within the range (1035−1040 g cm s−1) reported by Blackman & Lucchini (2014) for a sample of PPNe showing high-velocity and extreme high-velocity outflows. Finally, we derive the mass-loss rate of |$\dot{M}_j = \pi \varpi _0^2\rho V_{\rm H_{2}O} \sim 1.58\times 10^{-05}$| M⊙ yr−1. For these calculations, we have considered a constant density and velocity along the jet axis and along the jet radius as well. In the next section, we will discuss how to derive averaged quantities in physical units from a single scale-invariant streamline, which is what we call a solution. 4.3 Scaling of the solutions We use a sample of roughly 1500 solutions mapping the parameter space and we start by selecting the ones that satisfy the criterion: $$\begin{eqnarray*} \left(\frac{\varpi }{z}\right)_{\rm H_{2}O} = \tan (\theta _{\rm H_{2}O}), \end{eqnarray*}$$(27) which is equivalent to determining which solutions terminate beyond the H2O spots (⁠|$\theta _{\rm LRP} \lt \theta _{\rm H_{2}O}$|⁠). The quantity that determines |$\theta _{\rm H_{2}O}$| is the jet cylindrical radius at |$z_{\rm H_{2}O} = 1000$| au, while we keep the jet height at the water maser spots constant. Observational constraints on the cylindrical radius at H2O vary from ∼10 au (Imai et al. 2002) to ∼45 au (Tafoya et al. 2020). We treat the cylindrical radius at H2O as a free parameter within the above interval. Since our solutions are calculated in dimensionless units, e.g. (ϖ/ϖ*, V/V*, B/B*, ρ/ρ*, …), the first step to compare them with a physical system is to introduce a characteristic length to scale them. We use the cylindrical radius of the jet at the H2O masers spots as reference length to scale the cylindrical radius of each solution, ϖ, through the relation $$\begin{eqnarray*} \varpi _* = \left(\frac{\varpi }{G}\right)_{\rm H_{2}O} \end{eqnarray*}$$(28) for the reference streamline (α = 1, see Section 2.2). Once the reference length is fixed, the scaling of the velocity is also defined as $$\begin{eqnarray*} V_* = \sqrt{\frac{\mathcal {G M}}{\kappa _{\rm VTST}^2\varpi _*}} \end{eqnarray*}$$(29) for a given central object with mass (for a white dwarf, |$\mathcal {M} = 0.6\, \mathrm{ M}_\odot$|⁠). We use the observational constraints on the toroidal magnetic field component to determine the maximum and minimum B* as $$\begin{eqnarray*} B_* = \frac{B_{\phi ,\rm {obs}}}{B_{\phi ,\rm {H_{2}O}}}, \end{eqnarray*}$$(30) where |$B_{\phi ,\rm {obs}}$| are the values of the updated magnetic field limits, 0.6 and 100 mG, discussed in Section 4.1 and |$B_{\phi ,\rm {H_{2}O}}$| is the value of the toroidal component of the magnetic field for the reference line of our solutions at the position of the H2O maser spot. Then we introduce a third value of |$B_{\phi , \rm obs}$| that matches the momentum rate deduced from observational constraints. Finally, we derive the scaling for the mass density and the pressure from the above as follows: $$\begin{eqnarray*} \rho _* = \frac{B_*^2}{4\pi V_*^2}\qquad {\rm and} \qquad P_* = \frac{\mu _{\rm VTST} B_*^2 }{8\pi }. \end{eqnarray*}$$(31) 4.4 Integrated quantities As is generally the case for self-similar models, the properties of the jet at a given radius are derived from the reference streamline and extended with the appropriate radial dependence in the form of power law of the parameter α (see Section 2; Ferreira & Pelletier 1993; Vlahakis & Tsinganos 1998). This parameter is defined as $$\begin{eqnarray*} \alpha = \left(\frac{\varpi _\alpha }{\varpi _\star }\right)^2 = \left(\frac{\varpi (\theta)}{\varpi _\star G(\theta)}\right)^2 , \end{eqnarray*}$$(32) where ϖ(θ) is the radial profile of the reference streamline, G(θ) = ϖ(θ)/ϖα, ϖα is the cylindrical radius at the AP for the streamline with a given α and ϖ⋆ is the scaling length defined in equation (28) using the criterion described in equation (27). We give the radial scalings for all the relevant quantities in Appendix B. Using these relations, a given solution can be extended to infinity and towards the polar axis. Expanding a solution over the radial direction is necessary to calculate quantities such as the jet mass-loss and the momentum rate which require an integration over a surface perpendicular to the jet axis. Since the geometry of the equations that we adopted has a singularity on the polar axis, for the following calculations of integrated quantities we will consider a flux tube defined by inner and outer cylindrical radii, ϖin and ϖout, or, equivalently, αin and αout. Once that the scaling length ϖ* is defined, the inner and outer radii are determined and so are also the streamline labels αin and αout, through the equation (32). First, we evaluate a density-weighted average velocity for each jet solution at the height of the H2O maser spots over the flux tube area as follows: $$\begin{eqnarray*} \langle V_{\rm H_2O}\rangle = \frac{ 2\pi z_{\rm H_{2}O}^2\int _{\theta _{\rm LRP}}^{\theta _{\rm H_{2}O}} \rho V\frac{\sin (\theta)}{\cos ^3(\theta)}\mathrm{ d}\theta }{ 2\pi z_{\rm H_{2}O}^2 \int _{\theta _{\rm LRP}}^{\theta _{\rm H_{2}O}} \rho \frac{\sin (\theta)}{\cos ^3(\theta)}\mathrm{ d}\theta }, \end{eqnarray*}$$(33) where the relation between θ and α (or ϖ) is defined as $$\begin{eqnarray*} \alpha (\theta) = \left(\frac{z_{\rm H_{2}O}\tan (\theta)}{\varpi (\theta)}\right)^2, \end{eqnarray*}$$(34) where ϖ(θ) is calculated on the reference streamline with α = 1. We note that, since the velocity decreases with increasing α, we expect this average to be dominated by the inner streamlines in the flux tube, while the streamlines close to the outer edge of the flux tube will be slower. For this reason, it is more meaningful to compare the density-averaged velocity with the observed (almost constant) velocity. The mass-loss rate of the jet can be derived as the mass flux flowing from the z = 0 surface of the flux tube as follows: $$\begin{eqnarray*} \dot{M_j} = \int _{\varpi _{\rm in}}^{\varpi _{\rm out}} \rho V \cdot \mathrm{ d}A. \end{eqnarray*}$$(35) The mass-loss rate is dependent on the value of the toroidal magnetic field we are considering. Since there is still considerable uncertainty on the strength of the toroidal component of the magnetic fields, we can associate to each solution three mass-loss rates corresponding to the minimum and maximum Bϕ in Section 4.1 and the minimum value of Bϕ for which we find matching solutions (see Fig. 19). Figure 19. Open in new tabDownload slide Total jet height versus momentum rate for a jet width of 20 au at H2O for all the solutions in the sample. The blue squares and pink triangles represent the solutions |$\dot{P}$| for a toroidal magnetic field of 0.6 and 100 mG, respectively. The yellow dots are the momentum rates of the solutions with the minimum Bϕ (14 mG) that matches the observed |$\dot{P}$|⁠. The light yellow areas highlights the allowed momentum rates within this interval of Bϕ. The shaded grey areas show the acceptance intervals for zLRP and |$\dot{P}$|⁠. Figure 19. Open in new tabDownload slide Total jet height versus momentum rate for a jet width of 20 au at H2O for all the solutions in the sample. The blue squares and pink triangles represent the solutions |$\dot{P}$| for a toroidal magnetic field of 0.6 and 100 mG, respectively. The yellow dots are the momentum rates of the solutions with the minimum Bϕ (14 mG) that matches the observed |$\dot{P}$|⁠. The light yellow areas highlights the allowed momentum rates within this interval of Bϕ. The shaded grey areas show the acceptance intervals for zLRP and |$\dot{P}$|⁠. Similarly, we will give three values for the momentum rate of each jet configurations. The momentum rate of a jet model is $$\begin{eqnarray*} \dot{P} = 2\pi z_{\rm H_{2}O}^2 \int _{\theta _{\rm LRP}}^{\theta _{\rm H_{2}O}} \rho V^2 \frac{\sin (\theta)}{\cos ^3(\theta)}\mathrm{ d}\theta , \end{eqnarray*}$$(36) where the integration is done over all the streamlines contributing to the flux tube above the H2O maser spot. 4.5 Comparison results In Fig. 19, we present the total jet height (zLRP) versus the momentum rate of all the solutions in our sample. The black diamond marks the observed |$\dot{P}$| at the observed total jet height. The shaded grey horizontal and vertical areas show the intervals for |$\dot{P}$| and zLRP we use to define a solution as a good match. The shaded light yellow area between the blue squares and the magenta triangles define the values of the momentum rates that are allowed within the range of the toroidal magnetic field derived from observations. We see that for any value of Bϕ the solutions fall on a curve with little-to-none scattering introduced by the variation of the other jet properties. We produce this plot once we have set the half-width of the jet, but before introducing the other constraints on velocity and density and we find that a toroidal magnetic field at the H2O maser spots of at least 14 mG is required for the jet solutions to have a comparable or higher momentum rate than the observed one. Among the solutions found at the interception of the shaded grey areas in Fig. 19 for the given choice of the jet radius (⁠|$\varpi_{\rm H_2O}$| = 20 au) and toroidal magnetic field (Bϕ = 14 mG), we present a sample of eight jet configurations which satisfy all the constraints on density, velocity, total jet height and momentum rate. We report the parameters and the relevant scaled quantities in Table 4. Given the observational constraints (Section 4.1) and the tight correlation that exist between the jet total extent and the momentum rate, we are left with solutions having the same angular position and collimation angle at the AP, θA = 14° and ψA = 79° (and ψ0 = 35°), respectively, for the given choice of the parameter F (F = 0.75) and the polytropic index of the gas (Γ = 5/3). As a reference, we give the typical BP solution parameters (kBP = 0.03, λBP = 30, ψ0 = 32°) and we report in Table 4 our parameters in BP units. We remind the reader that the equations that we adopted differ from the classical BP because we do not neglect the enthalpy of the gas. Table 4. Parameters of the selected solutions as a good match to W43A. All solutions have models parameters F = 0.75, Γ = 5/3, θA = 14°, ψA = 79° and scaling parameters |$\varpi_{\rm H_2O}$| = 20 au, |$z_{\rm H_2O}$| = 1000 au, and Bϕ = 14 mG. Model . S1 . S2 . S3 . S4 . S5 . S6 . S7 . S8 . kVTST 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 μVTST 0.7957 1.3122 1.9295 2.6433 3.4504 4.3488 5.3373 6.4149 λVTST 0.3155 0.3185 0.3214 0.3239 0.3261 0.3280 0.3295 0.3307 ϵVTST 11.2599 11.4833 11.7117 11.9394 12.1626 12.3788 12.5862 12.7840 θMFP 2.1866 × 10−2 2.1753 × 10−2 2.1623 × 10−2 2.1480 × 10−2 2.1327 × 10−2 2.1166 × 10−2 2.1000 × 10−2 2.0830 × 10−2 θMSP 0.7494 0.6046 0.5194 0.4631 0.4232 0.3938 0.3712 0.3535 |$\langle V_{\rm H_{2}O} \rangle$| (km s−1) 1982 1510 1238 1059 934 843 775 722 〈Vp,MSP〉 (km s−1) 1405 1057 859 731 642 578 530 494 |$\langle n_{\rm H_{2}O} \rangle$| (cm−3) 8.05 × 105 1.39 × 106 2.10 × 106 2.87 × 106 3.72 × 106 4.60 × 106 5.49 × 106 6.37 × 106 |$\dot{M}_j$| (M⊙ yr−1) 1.24 × 10−6 1.58 × 10−6 1.95 × 10−6 2.24 × 10−6 2.64 × 10−6 2.87 × 10−6 3.08 × 10−6 3.26 × 10−6 |$\dot{P}$| (M⊙ yr−1 km s−1) 2.38 × 10−3 2.41 × 10−3 2.41 × 10−3 2.42 × 10−3 2.43 × 10−3 2.45 × 10−3 2.48 × 10−3 2.51 × 10−3 P (g cm s−1) 3.07 × 1037 3.12 × 1037 3.12 × 1037 3.13 × 1037 3.14 × 1037 3.17 × 1037 3.20 × 1037 3.24 × 1037 kBP 0.0032 0.0042 0.0052 0.0062 0.0071 0.0080 0.0088 0.0097 μBP 13.1344 8.3814 5.9173 4.4605 3.5194 2.8716 2.4038 2.0532 λBP 1.6307 1.2404 1.0059 0.8491 0.7364 0.6513 0.5844 0.5303 ϵBP 8.3257 × 10−2 4.7334 × 10−2 3.0601 × 10−2 2.1451 × 10−2 1.5895 × 10−2 1.2265 × 10−2 9.7593 × 10−3 7.9553 × 10−3 Model . S1 . S2 . S3 . S4 . S5 . S6 . S7 . S8 . kVTST 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 μVTST 0.7957 1.3122 1.9295 2.6433 3.4504 4.3488 5.3373 6.4149 λVTST 0.3155 0.3185 0.3214 0.3239 0.3261 0.3280 0.3295 0.3307 ϵVTST 11.2599 11.4833 11.7117 11.9394 12.1626 12.3788 12.5862 12.7840 θMFP 2.1866 × 10−2 2.1753 × 10−2 2.1623 × 10−2 2.1480 × 10−2 2.1327 × 10−2 2.1166 × 10−2 2.1000 × 10−2 2.0830 × 10−2 θMSP 0.7494 0.6046 0.5194 0.4631 0.4232 0.3938 0.3712 0.3535 |$\langle V_{\rm H_{2}O} \rangle$| (km s−1) 1982 1510 1238 1059 934 843 775 722 〈Vp,MSP〉 (km s−1) 1405 1057 859 731 642 578 530 494 |$\langle n_{\rm H_{2}O} \rangle$| (cm−3) 8.05 × 105 1.39 × 106 2.10 × 106 2.87 × 106 3.72 × 106 4.60 × 106 5.49 × 106 6.37 × 106 |$\dot{M}_j$| (M⊙ yr−1) 1.24 × 10−6 1.58 × 10−6 1.95 × 10−6 2.24 × 10−6 2.64 × 10−6 2.87 × 10−6 3.08 × 10−6 3.26 × 10−6 |$\dot{P}$| (M⊙ yr−1 km s−1) 2.38 × 10−3 2.41 × 10−3 2.41 × 10−3 2.42 × 10−3 2.43 × 10−3 2.45 × 10−3 2.48 × 10−3 2.51 × 10−3 P (g cm s−1) 3.07 × 1037 3.12 × 1037 3.12 × 1037 3.13 × 1037 3.14 × 1037 3.17 × 1037 3.20 × 1037 3.24 × 1037 kBP 0.0032 0.0042 0.0052 0.0062 0.0071 0.0080 0.0088 0.0097 μBP 13.1344 8.3814 5.9173 4.4605 3.5194 2.8716 2.4038 2.0532 λBP 1.6307 1.2404 1.0059 0.8491 0.7364 0.6513 0.5844 0.5303 ϵBP 8.3257 × 10−2 4.7334 × 10−2 3.0601 × 10−2 2.1451 × 10−2 1.5895 × 10−2 1.2265 × 10−2 9.7593 × 10−3 7.9553 × 10−3 Open in new tab Table 4. Parameters of the selected solutions as a good match to W43A. All solutions have models parameters F = 0.75, Γ = 5/3, θA = 14°, ψA = 79° and scaling parameters |$\varpi_{\rm H_2O}$| = 20 au, |$z_{\rm H_2O}$| = 1000 au, and Bϕ = 14 mG. Model . S1 . S2 . S3 . S4 . S5 . S6 . S7 . S8 . kVTST 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 μVTST 0.7957 1.3122 1.9295 2.6433 3.4504 4.3488 5.3373 6.4149 λVTST 0.3155 0.3185 0.3214 0.3239 0.3261 0.3280 0.3295 0.3307 ϵVTST 11.2599 11.4833 11.7117 11.9394 12.1626 12.3788 12.5862 12.7840 θMFP 2.1866 × 10−2 2.1753 × 10−2 2.1623 × 10−2 2.1480 × 10−2 2.1327 × 10−2 2.1166 × 10−2 2.1000 × 10−2 2.0830 × 10−2 θMSP 0.7494 0.6046 0.5194 0.4631 0.4232 0.3938 0.3712 0.3535 |$\langle V_{\rm H_{2}O} \rangle$| (km s−1) 1982 1510 1238 1059 934 843 775 722 〈Vp,MSP〉 (km s−1) 1405 1057 859 731 642 578 530 494 |$\langle n_{\rm H_{2}O} \rangle$| (cm−3) 8.05 × 105 1.39 × 106 2.10 × 106 2.87 × 106 3.72 × 106 4.60 × 106 5.49 × 106 6.37 × 106 |$\dot{M}_j$| (M⊙ yr−1) 1.24 × 10−6 1.58 × 10−6 1.95 × 10−6 2.24 × 10−6 2.64 × 10−6 2.87 × 10−6 3.08 × 10−6 3.26 × 10−6 |$\dot{P}$| (M⊙ yr−1 km s−1) 2.38 × 10−3 2.41 × 10−3 2.41 × 10−3 2.42 × 10−3 2.43 × 10−3 2.45 × 10−3 2.48 × 10−3 2.51 × 10−3 P (g cm s−1) 3.07 × 1037 3.12 × 1037 3.12 × 1037 3.13 × 1037 3.14 × 1037 3.17 × 1037 3.20 × 1037 3.24 × 1037 kBP 0.0032 0.0042 0.0052 0.0062 0.0071 0.0080 0.0088 0.0097 μBP 13.1344 8.3814 5.9173 4.4605 3.5194 2.8716 2.4038 2.0532 λBP 1.6307 1.2404 1.0059 0.8491 0.7364 0.6513 0.5844 0.5303 ϵBP 8.3257 × 10−2 4.7334 × 10−2 3.0601 × 10−2 2.1451 × 10−2 1.5895 × 10−2 1.2265 × 10−2 9.7593 × 10−3 7.9553 × 10−3 Model . S1 . S2 . S3 . S4 . S5 . S6 . S7 . S8 . kVTST 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 μVTST 0.7957 1.3122 1.9295 2.6433 3.4504 4.3488 5.3373 6.4149 λVTST 0.3155 0.3185 0.3214 0.3239 0.3261 0.3280 0.3295 0.3307 ϵVTST 11.2599 11.4833 11.7117 11.9394 12.1626 12.3788 12.5862 12.7840 θMFP 2.1866 × 10−2 2.1753 × 10−2 2.1623 × 10−2 2.1480 × 10−2 2.1327 × 10−2 2.1166 × 10−2 2.1000 × 10−2 2.0830 × 10−2 θMSP 0.7494 0.6046 0.5194 0.4631 0.4232 0.3938 0.3712 0.3535 |$\langle V_{\rm H_{2}O} \rangle$| (km s−1) 1982 1510 1238 1059 934 843 775 722 〈Vp,MSP〉 (km s−1) 1405 1057 859 731 642 578 530 494 |$\langle n_{\rm H_{2}O} \rangle$| (cm−3) 8.05 × 105 1.39 × 106 2.10 × 106 2.87 × 106 3.72 × 106 4.60 × 106 5.49 × 106 6.37 × 106 |$\dot{M}_j$| (M⊙ yr−1) 1.24 × 10−6 1.58 × 10−6 1.95 × 10−6 2.24 × 10−6 2.64 × 10−6 2.87 × 10−6 3.08 × 10−6 3.26 × 10−6 |$\dot{P}$| (M⊙ yr−1 km s−1) 2.38 × 10−3 2.41 × 10−3 2.41 × 10−3 2.42 × 10−3 2.43 × 10−3 2.45 × 10−3 2.48 × 10−3 2.51 × 10−3 P (g cm s−1) 3.07 × 1037 3.12 × 1037 3.12 × 1037 3.13 × 1037 3.14 × 1037 3.17 × 1037 3.20 × 1037 3.24 × 1037 kBP 0.0032 0.0042 0.0052 0.0062 0.0071 0.0080 0.0088 0.0097 μBP 13.1344 8.3814 5.9173 4.4605 3.5194 2.8716 2.4038 2.0532 λBP 1.6307 1.2404 1.0059 0.8491 0.7364 0.6513 0.5844 0.5303 ϵBP 8.3257 × 10−2 4.7334 × 10−2 3.0601 × 10−2 2.1451 × 10−2 1.5895 × 10−2 1.2265 × 10−2 9.7593 × 10−3 7.9553 × 10−3 Open in new tab We notice that the only other parameter that has not been constrained is the mass-loss parameter kVTST. This parameter is proportional to the mass-to-magnetic flux ratio, which leads to an uncertainty on the dimensionless angular momentum, λVTST, and the parameter regulating the entropy of the gas, μVTST. Such spread in the values of the above parameters is reflected in the uncertainties on the average velocities and densities at the maser spot. However, the resulting momentum rate among the selected models is roughly constant and very close to the momentum rate estimated from the observations (⁠|$\dot{P} \sim 0.0024{\!-\!}0.0025$| M⊙ yr−1 km s−1). We note that substantially larger |$\dot{P}$| could be achieved by adopting a larger toroidal magnetic field. A shorter jet (<2000 au) would still require a larger Bϕ, while a taller jet could yield a larger momentum rate for the same choice of Bϕ. In Fig. 20, we show the velocity profiles of the jet models given in Table 4. Given the model assumptions given in Section 4.2, we selected solutions that have velocity profiles of the reference (α = 1) outermost line of the flux tube largely exceeding the average observed velocity of 150 km/s (horizontal black line), so that the MHD jet would be able to transfer momentum to and accelerate the molecular cocoon. The inner streamlines (α < 1) have the same acceleration profile, but higher speeds due to the relations (B5) and (B7) given in Appendix B. We notice that the smaller the mass load parameter, kVTST, the larger is the speed. The uncertainty left on kVTST, and therefore on the velocity of the MHD jet layer, can only be removed with further observations of the core of the molecular jet of W43A. Figure 20. Open in new tabDownload slide Velocity versus jet height of the scaled reference streamlines for α = 1, which are the outermost streamlines of our selected jet configurations. The solid horizontal line is the observed velocity at the H2O masers. The vertical thin solid line is the ‘base’ of the jet as seen in the CO observations, and the vertical dashed solid line is the location of the H2O maser spot. The thick solid vertical line is at z = 2000 au. Figure 20. Open in new tabDownload slide Velocity versus jet height of the scaled reference streamlines for α = 1, which are the outermost streamlines of our selected jet configurations. The solid horizontal line is the observed velocity at the H2O masers. The vertical thin solid line is the ‘base’ of the jet as seen in the CO observations, and the vertical dashed solid line is the location of the H2O maser spot. The thick solid vertical line is at z = 2000 au. While there is a moderate acceleration taking place from the MSP to shortly downstream of the AP, in the portion of the jet observed through the emission of CO, i.e. from 45 to 2000 au (the region between the thin vertical solid line and the thick vertical solid line in Fig. 20), the velocity has already reached its maximum and it stays constant up until the jet tip. The total velocity is entirely poloidal, while the toroidal component is close to zero along the entire jet extent. Under these circumstances, the magnetic field and the gas are not corotating even upstream of the Alfvén surface, which is an indication of a jet driven by thermal pressure (see bottom panels of Figs 13–15) as opposed to a magnetically driven jet (top and middle panels of Figs 13–15). Typically these winds are less powerful and they can only achieve higher speeds if a large injection speed (Vp,MSP) is provided (see Table 4). Such high injections speeds are consistent or higher than the initial speeds considered in recent MHD simulations by Balick, Frank & Liu (2020), which are successful in reproducing the qualitative shapes of a sample of PPNe. In order to make this comparison as complete as possible, we investigate the effect of varying the radius of the jet and of having a main-sequence star as the accreting object. Increasing or decreasing the jet radius within the observed range 10–45 au has the effect of decreasing/increasing the angular position of the AP and increasing/decreasing its collimation angle. The thinnest jet (⁠|$\varpi _{\rm H_{2}O} = 10$| au, θA = 10°, ψA = 83°) allows two values of the mass-loss parameter (1.5 and 2.0) instead of eight, limiting the selection to just two models. The thickest jet (⁠|$\varpi _{\rm H_{2}O} = 45$| au) leaves us solutions with θA (∼25°) and ψA (68°) and excludes kVTST = 1.5. We also considered the possibility that the central object may be a main-sequence star (M ∼ M⊙ and R⋆ ∼ R⊙), finding our conclusions unaltered, as expected by the mild dependence that our scaling scheme has on the mass of the central star. 5 SUMMARY In this paper, we discussed the adaptation of the numerical algorithm we presented in Paper I to solve the non-relativistic, radial self-similar MHD equations describing a disc-driven outflow. We focused on the study of a large sample of solutions defined by constant Blandford–Payne-like parameter F (F = 0.75) and polytropic index Γ = 5/3. We recognized similar patterns within the collection of jet configurations that are ultimately ascribed to the cold-to-hot transition that we find recurrently for similar values of the angular position of the AP and the collimation angle at the same position. We analysed the behaviour of all the relevant jet quantities undergoing this transition and found that Cold jets have the largest (dimensionless) angular momentum and they have the lowest enthalpy and plasma-β much lower than unity. They are therefore magnetically dominated jets. They have little-to-none vertical speed upstream of the MSP, but have |Bϕ/Bp| ratios larger than unity. This combination produces twisted streamlines with variable radius, due to the oscillatory behaviour of the transverse forces. The highly wound up magnetic field is also responsible for the relatively high-mass load (η ≲ 1) of these solutions. At approximately half-way between the AP and the MFP, a large fraction of the magnetic energy has turned into kinetic energy and the jet becomes kinetically dominated until the LRP. MC jets are similar to cold jets; however, the enthalpy is slightly larger, and it plays a role in lifting the gas. These jet models do not suffer oscillations upstream of the MSP. From this point on, these models resemble closely the cold jets. They are, however, the most efficient at accelerating the flow, even though their total energy flux is lower than the purely cold jets. Hot jets are thermally dominated jet configurations (plasma-β = P/(B2/8π) ≲ 1), where the magnetic field is contributing significantly to the acceleration and collimation of the jet only downstream of the AP. These solutions start off with a large poloidal speed, negative toroidal velocity and |Bϕ/Bp| < <1. The acceleration is only mild and they have low-energy flux densities. Within this regime, we see two types of energy transfer channels that lead to an increase in kinetic energy. In hot jets, the gas pressure is responsible for the acceleration in the initial jet segment, which can extend even just downstream of the AP. A fraction of the thermal energy is transferred to the magnetic energy, which then is used for the last acceleration until the tip of the jet. We then describe a procedure for the identification of specific jet solutions to be compared to an astrophysical source, in our case the water fountain W43A. W43A is believed to be an AGB star in the process of becoming a PN. During this current, short-lived phase, the source is launching collimated molecular jets, the nature of which is debated. We assume that the jets of W43A is launched by a disc-driven ionized inner shell, which is surrounded by a molecular jet sheath. Since the true nature or even the existence of such a jet core is unknown, we adopted the constraints on the molecular gas as upper/lower limits to the corresponding quantities of the atomic jet, namely the size, hydrogen number density, velocity, and magnetic field strength to identify possible jet configuration. We conducted an exhaustive examination of our sample and we established that, given the observed molecular properties within the jets of W43A and our (large, but finite) collection of solutions, the most suitable jet model for the jets of W43A is a thermally dominated jet configuration with a high injection speed, but not efficiently accelerating for most of its extent. We found that the strength of the toroidal component of the magnetic field is the parameter that affects the most this comparison. This procedure can be applied to other sources, for which the magnetic field has been difficult to measure, in order to determine a range of plausible magnetic field strengths given observational constraints on density, velocity, and jet size. In future works, we will expand our grid of solutions to the additional two dimensions, namely the radial scaling of the current, F, and the polytropic index, Γ, and compare the full sample to other astrophysical sources. ACKNOWLEDGEMENTS Chiara Ceccobello and Wouter Vlemmings acknowledge support from the Swedish Research Council (VR). DATA AVAILABILITY A catalogue of all the solutions that have been found is available on request to the main author. Footnotes 1 We will use Γ for the polytropic index and F for the power-law exponent, instead of the symbols γ and x as was done in VTST00 to maintain the same convention we had for the relativistic equations in Paper I. 2 Note that after the separation of variables, they are points (not surfaces) on a single streamline and they are modified because their position and definition of the phase speeds of the slow and fast magnetosonic waves are affected by the geometry of the magnetic field (e.g. Sauty & Tsinganos 1994; Ferreira & Pelletier 1995). REFERENCES Amiri N. , Vlemmings W., van Langevelde H. J., 2010 , A&A , 509 , A26 10.1051/0004-6361/200913194 Crossref Search ADS Anderson J. M. , Li Z.-Y., Krasnopolsky R., Blandford R. D., 2005 , ApJ , 630 , 945 10.1086/432040 Crossref Search ADS Balick B. , Frank A., 2002 , ARA&A , 40 , 439 10.1146/annurev.astro.40.060401.093849 Crossref Search ADS Balick B. , Frank A., Liu B., 2020 , ApJ , 889 , 13 10.3847/1538-4357/ab5651 Crossref Search ADS Blackman E. G. , Lucchini S., 2014 , MNRAS , 440 , L16 10.1093/mnrasl/slu001 Crossref Search ADS Blandford R. D. , Payne D. G., 1982 , MNRAS , 199 , 883 (BP) 10.1093/mnras/199.4.883 Crossref Search ADS Blandford R. D. , Znajek R. L., 1977 , MNRAS , 179 , 433 10.1093/mnras/179.3.433 Crossref Search ADS Bogovalov S. , Tsinganos K., 1999 , MNRAS , 305 , 211 10.1046/j.1365-8711.1999.02413.x Crossref Search ADS Cayatte V. , Vlahakis N., Matsakos T., Lima J. J. G., Tsinganos K., Sauty C., 2014 , ApJ , 788 , L19 10.1088/2041-8205/788/1/L19 Crossref Search ADS Ceccobello C. , Cavecchi Y., Heemskerk M. H. M., Markoff S., Polko P., Meier D., 2018 , MNRAS , 473 , 4417 (Paper I) 10.1093/mnras/stx2567 Crossref Search ADS Chamandy L. et al. , 2018 , MNRAS , 480 , 1898 10.1093/mnras/sty1950 Crossref Search ADS Chantry L. , Cayatte V., Sauty C., Vlahakis N., Tsinganos K., 2018 , A&A , 612 , A63 10.1051/0004-6361/201731793 Crossref Search ADS Cohen M. H. et al. , 2014 , ApJ , 787 , 151 10.1088/0004-637X/787/2/151 Crossref Search ADS Contopoulos J. , Lovelace R. V. E., 1994 , ApJ , 429 , 139 10.1086/174307 Crossref Search ADS Davelaar J. et al. , 2019 , A&A , 632 , A2 10.1051/0004-6361/201936150 Crossref Search ADS Duran-Rojas M. et al. , 2014 , in Morisset C., Delgado-Inglada G., Torres-Peimbert S., eds, Asymmetrical Planetary Nebulae VI Conference . p. 19 Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Elitzur M. , Hollenbach D. J., McKee C. F., 1992 , ApJ , 394 , 221 10.1086/171574 Crossref Search ADS Feroz F. , Hobson M. P., 2008 , MNRAS , 384 , 449 10.1111/j.1365-2966.2007.12353.x Crossref Search ADS Feroz F. , Hobson M. P., Bridges M., 2009 , MNRAS , 398 , 1601 10.1111/j.1365-2966.2009.14548.x Crossref Search ADS Feroz F. , Hobson M. P., Cameron E., Pettitt A. N., 2019 , The Open J. Astrophys. , 2 , 10 10.21105/astro.1306.2144 Crossref Search ADS Ferreira J. , 1997 , A&A , 319 , 340 Ferreira J. , Pelletier G., 1993 , A&A , 276 , 625 Ferreira J. , Pelletier G., 1995 , A&A , 295 , 807 García-Segura G. , López J. A., Franco J., 2005 , ApJ , 618 , 919 10.1086/426110 Crossref Search ADS Hardee P. E. , 1996 , in Hardee P. E., Bridle A. H., Zensus J. A., eds, ASP Conf. Ser. Vol. 100, Energy Transport in Radio Galaxies and Quasars . Astron. Soc. Pac , San Francisco , p. 273 Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Imai H. , Diamond P. J., 2005 , in Romney J., Reid M., eds, ASP Conf. Ser. Vol. 340, Future Directions in High Resolution Astronomy . Astron. Soc. Pac , San Francisco , p. 399 Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Imai H. , Obara K., Diamond P. J., Omodaka T., Sasao T., 2002 , Nature , 417 , 829 10.1038/nature00788 Crossref Search ADS PubMed Komissarov S. S. , Vlahakis N., Königl A., Barkov M. V., 2009 , MNRAS , 394 , 1182 10.1111/j.1365-2966.2009.14410.x Crossref Search ADS Li Z.-Y. , Chiueh T., Begelman M. C., 1992 , ApJ , 394 , 459 10.1086/171597 Crossref Search ADS Liska M. , Hesp C., Tchekhovskoy A., Ingram A., van der Klis M., Markoff S., 2018 , MNRAS , 474 , L81 10.1093/mnrasl/slx174 Crossref Search ADS Markoff S. , 2010 , in Belloni T., ed., Lecture Notes in Physics, Vol. 794, The Jet Paradigm – From Microquasars to Quasars . Springer-Verlag , Berlin , p. 143 10.1007/978-3-540-76937-8_6 Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Crossref Markoff S. , Falcke H., Fender R., 2001 , A&A , 372 , L25 10.1051/0004-6361:20010420 Crossref Search ADS Markoff S. , Nowak M. A., Wilms J., 2005 , ApJ , 635 , 1203 10.1086/497628 Crossref Search ADS Meier D. L. , 2012 , Black Hole Astrophysics: The Engine Paradigm , Springer-Verlag , Berlin Heidelberg Mościbrodzka M. , Falcke H., Shiokawa H., 2016 , A&A , 586 , A38 10.1051/0004-6361/201526630 Crossref Search ADS Nordhaus J. , Blackman E. G., 2006 , MNRAS , 370 , 2004 10.1111/j.1365-2966.2006.10625.x Crossref Search ADS Parfrey K. , Spitkovsky A., Beloborodov A. M., 2016 , ApJ , 822 , 33 10.3847/0004-637X/822/1/33 Crossref Search ADS Parker E. N. , 1958 , ApJ , 128 , 664 10.1086/146579 Crossref Search ADS Pelletier G. , Pudritz R. E., 1992 , ApJ , 394 , 117 10.1086/171565 Crossref Search ADS Polko P. , Meier D. L., Markoff S., 2010 , ApJ , 723 , 1343 10.1088/0004-637X/723/2/1343 Crossref Search ADS Polko P. , Meier D. L., Markoff S., 2013 , MNRAS , 428 , 587 10.1093/mnras/sts052 Crossref Search ADS Polko P. , Meier D. L., Markoff S., 2014 , MNRAS , 438 , 959 10.1093/mnras/stt2155 Crossref Search ADS Pudritz R. E. , Norman C. A., 1983 , ApJ , 274 , 677 10.1086/161481 Crossref Search ADS Rosen A. , Hardee P. E., Clarke D. A., Johnson A., 1999 , ApJ , 510 , 136 10.1086/306565 Crossref Search ADS Sahai R. , Trauger J. T., 1998 , AJ , 116 , 1357 10.1086/300504 Crossref Search ADS Sahai R. , Vlemmings W. H. T., Gledhill T., Sánchez Contreras C., Lagadec E., Nyman L. Å., Quintana-Lacaci G., 2017 , ApJ , 835 , L13 10.3847/2041-8213/835/1/L13 Crossref Search ADS Sauty C. , Tsinganos K., 1994 , A&A , 287 , 893 Sauty C. , Tsinganos K., Trussoni E., 1999 , A&A , 348 , 327 Shu F. , Najita J., Ostriker E., Wilkin F., Ruden S., Lizano S., 1994 , ApJ , 429 , 781 10.1086/174363 Crossref Search ADS Shu F. H. , Najita J. R., Shang H., Li Z. Y., 2000 , in Mannings V., Boss A. P., Russell S. S., eds, Protostars and Planets IV . Univ. Arizona Press , Tucson, AZ , p. 789 Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Soker N. , 2020 , Galaxies , 8 , 26 10.3390/galaxies8010026 Crossref Search ADS Spruit H. C. , 1996 , in Wijers R. A. M. J., Davies M. B., Tout C. A., eds, NATO Advanced Science Institutes (ASI) Series C, Vol. 477 . Kluwer Academic Publ , p. 249 Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Tafoya D. , Imai H., Gómez J. F., Nakashima J.-i., Orosz G., Yung B. H. K., 2020 , ApJ , 890 , L14 10.3847/2041-8213/ab70b8 Crossref Search ADS Trussoni E. , Tsinganos K., Sauty C., 1997 , A&A , 325 , 1099 Vlahakis N. , Königl A., 2003 , ApJ , 596 , 1080 10.1086/378226 Crossref Search ADS Vlahakis N. , Tsinganos K., 1998 , MNRAS , 298 , 777 10.1046/j.1365-8711.1998.01660.x Crossref Search ADS Vlahakis N. , Tsinganos K., Sauty C., Trussoni E., 2000 , MNRAS , 318 , 417 (VTST00) 10.1046/j.1365-8711.2000.03703.x Crossref Search ADS Vlemmings W. H. T. , Diamond P. J., 2006 , ApJ , 648 , L59 10.1086/507679 Crossref Search ADS Vlemmings W. H. T. , Diamond P. J., Imai H., 2006 , Nature , 440 , 58 10.1038/nature04466 Crossref Search ADS PubMed Weber E. J. , Davis L. Jr, 1967 , ApJ , 148 , 217 10.1086/149138 Crossref Search ADS APPENDIX A: COEFFICIENTS OF THE BERNOULLI AND TRANSVERSE EQUATIONS In equation (5), we gave a general form to which both the Bernoulli and the transfield equations can be reduced. Here, we provide the explicit form of the coefficients Ai, Bi, Ci with i = 1, 2 for both equations as follows: $$\begin{eqnarray*} A_1 = \frac{\cos ^2(\psi +\theta)}{\sin ^2(\theta)}\left[-2\lambda _{\rm VTST}^2 \frac{M^2}{G^2} \frac{(1-G^2)^2}{(1-M^2)^3}+\Gamma \frac{\mu _{\rm VTST}}{M^{2\Gamma }}\right] - 2\frac{M^2}{G^4}, \end{eqnarray*}$$(A1) $$\begin{eqnarray*} B_1 = -2\frac{M^4}{G^4}\tan (\psi +\theta), \end{eqnarray*}$$(A2) $$\begin{eqnarray*} C_1 = 2k_{\rm VTST}^2\frac{\sin (\theta +\psi)\cos (\theta +\psi)}{G\sin (\theta)} - 2\frac{M^4}{G^4}\frac{\cos (\psi)}{\sin (\theta)\cos (\theta +\psi)} + 2\lambda _{\rm VTST}^2\frac{G^2(2M^2-1)-M^4}{(1-M^2)^2}\frac{\cos (\psi)\cos (\psi +\theta)}{\sin ^3(\theta)} , \end{eqnarray*}$$(A3) $$\begin{eqnarray*} A_2 = \sin (\theta +\psi)\cos (\theta +\psi)\left[\lambda _{\rm VTST}^2G^2\frac{(1-G^2)^2}{(1-M^2)^3} -\Gamma \mu _{\rm VTST}\frac{G^4}{2M^{2(\Gamma +1)}} \right] , \end{eqnarray*}$$(A4) $$\begin{eqnarray*} B_2 = \sin ^2(\theta)\left[\frac{1}{\cos ^2(\theta +\psi)} - M^2\right], \end{eqnarray*}$$(A5) $$\begin{eqnarray*} C_2 &=& \lambda _{\rm VTST}^2\frac{G^2}{M^2}\left(\frac{G^2-M^2}{1-M^2}\right)^2\frac{\sin (\psi)\cos (\psi +\theta)}{\sin (\theta)} +\frac{\sin ^2(\theta)}{\cos ^2(\psi +\theta)}\left[F-2+\frac{\cos (\psi)\sin (\theta +\psi)}{\sin (\theta)}\right] \nonumber \\ && +\lambda _{\rm VTST}^2G^2\frac{(1-G^2)^2}{(1-M^2)^2}\left[F-1+2\frac{G^2}{(1-G^2)}\frac{\cos (\psi)\sin (\psi +\theta)}{\sin (\theta)}\right] +\mu _{\rm VTST}(F-2)\frac{G^4}{M^{2\Gamma }} +k_{\rm VTST}^2\frac{G^3}{M^2}\sin (\theta)\cos ^2(\psi +\theta) . \end{eqnarray*}$$(A6) APPENDIX B: PHYSICAL SCALING We report here the definitions of the density and the velocity and magnetic field components from VTST00. $$\begin{eqnarray*} \rho = \frac{\rho _\star }{M^2}\alpha ^{F-3/2}, \quad P = \frac{P_\star }{M^{2\Gamma }}\alpha ^{F-2}, \quad {\rm with}\quad \alpha = \left(\frac{\varpi _\alpha }{\varpi _\star }\right)^2 = \left(\frac{\varpi }{\varpi _\star G}\right)^2 , \end{eqnarray*}$$(B1) $$\begin{eqnarray*} \boldsymbol{B_p} = - \frac{B_\star \alpha ^{(F-2)/2}}{G^2}\frac{\sin (\theta)}{\cos (\theta +\psi)}\hat{\boldsymbol{b}} = B_r \hat{\boldsymbol{r}} + B_\theta \hat{\boldsymbol{\theta }} = B_p\sin (\theta +\psi)\hat{\boldsymbol{r}}+B_p\cos (\theta +\psi)\hat{\boldsymbol{\theta }}, \end{eqnarray*}$$(B2) $$\begin{eqnarray*} B_p = - \frac{B_\star \alpha ^{(F-2)/2}}{G^2}\frac{\sin (\theta)}{\cos (\theta +\psi)}, \end{eqnarray*}$$(B3) $$\begin{eqnarray*} \boldsymbol{B_\phi } = - \lambda _{\rm VTST}\frac{B_\star \alpha ^{(F-2)/2}}{G}\frac{(1-G^2)}{(1-M^2)} \hat{\boldsymbol{\phi }}= B_\phi \hat{\boldsymbol{\phi }} , \end{eqnarray*}$$(B4) $$\begin{eqnarray*} \boldsymbol{V_p} = - \frac{V_\star \alpha ^{-1/4}M^2}{G^2}\frac{\sin (\theta)}{\cos (\theta +\psi)}\hat{\boldsymbol{b}} = V_r \hat{\boldsymbol{r}} + V_\theta \hat{\boldsymbol{\theta }} = V_p\sin (\theta +\psi)\hat{\boldsymbol{r}}+V_p\cos (\theta +\psi)\hat{\boldsymbol{\theta }}, \end{eqnarray*}$$(B5) $$\begin{eqnarray*} V_p = - \frac{V_\star \alpha ^{-1/4}M^2}{G^2}\frac{\sin (\theta)}{\cos (\theta +\psi)}, \end{eqnarray*}$$(B6) $$\begin{eqnarray*} \boldsymbol{V_\phi } = \lambda _{\rm VTST}\frac{V_\star \alpha ^{-1/4}}{G}\frac{(G^2-M^2)}{(1-M^2)} \hat{\boldsymbol{\phi }} = V_\phi \hat{\boldsymbol{\phi }} , \end{eqnarray*}$$(B7) $$\begin{eqnarray*} B_\star = \sqrt{4\pi \rho _\star }V_\star , \quad \quad P_\star = \mu _{\rm VTST} \frac{B_\star ^2}{8\pi }, \quad V_\star ^2 = \frac{\mathcal {GM}}{\varpi _\star k_{\rm VTST}^2} . \end{eqnarray*}$$(B8) If we want to compare our solutions to a real jets, we need to determine the scaling of the physical quantities for a given streamline. Following VTST00, we give here general scaling relations similar to the ones given in their section 5.1, which are tailored for a subclass of solutions (⁠|$M_0 \simeq 0.01, \, G_0 \simeq 0.1$|⁠, and F = 0.75). Using the definition of α in (B1) and of V⋆ in (B8), the Keplerian velocity at the footpoint of a streamline (z = 0) can be written as $$\begin{eqnarray*} V_{\rm kep,0} = \sqrt{\frac{\mathcal {GM}}{\varpi _0}} = k V_\star G_0^{-1/2} \alpha ^{-1/4} . \end{eqnarray*}$$(B9) The sound speed for a polytropic gas is defined as $$\begin{eqnarray*} C_{\rm s}^2 = \frac{\mathrm{ d} P}{\mathrm{ d} \rho } = \frac{\Gamma P}{\rho }, \end{eqnarray*}$$(B10) which can be calculated at z = 0 and recast by making use of the equations (B1) and (B8) as follows: $$\begin{eqnarray*} C_{\rm s,0}^2 = \frac{\mu _{\rm VTST}}{2} \Gamma V_\star ^2\alpha ^{-1/2} M_0^{2(1-\Gamma)}. \end{eqnarray*}$$(B11) The vertical component of the velocity at z = 0 is V0 ≡ Vz(z = 0) = Vp,0sin ψ0, or $$\begin{eqnarray*} V_0 \equiv V_z (z=0) = \alpha ^{-1/4}M_0^2 G_0^{-2} V_\star . \end{eqnarray*}$$(B12) With these quantities, we can calculate the following ratios: $$\begin{eqnarray*} \left(\frac{C_{\rm s}}{V_0}\right)_0 = \sqrt{\frac{\mu _{\rm VTST} \Gamma }{2}} M_0^{-(1+\Gamma)}G_0^{2}, \end{eqnarray*}$$(B13) $$\begin{eqnarray*} \left(\frac{V_{\rm kep}}{V_0}\right)_0 & = k_{\rm VTST}\, G_0^{3/2} M_0^{-2}, \end{eqnarray*}$$(B14) $$\begin{eqnarray*} \left(\frac{C_{\rm s}}{V_{\rm kep}}\right)_0 = \sqrt{\frac{\mu _{\rm VTST} \Gamma }{2 k_{\rm VTST}^2}} M_0^{(1-\Gamma)} G_0^{1/2}, \end{eqnarray*}$$(B15) and, as well, the reference length to scale the solutions to real objects $$\begin{eqnarray*} \varpi _0=\frac{\mathcal {GM}}{V_{\rm kep}^2} = \frac{\mathcal {GM}}{\left(\frac{V_{\rm kep}}{V_0}\right)^2 V_0^2} \simeq 1.91\times 10^5 k_{\rm VTST}^{-2} M_0^4 G_0^{-3} \left(\frac{\mathcal {M}}{\mathcal {M}_\odot }\right) \left(\frac{V_0}{\rm {km\, s^{-1}}}\right)^{-2} R_\odot . \end{eqnarray*}$$(B16) APPENDIX C: POLOIDAL AND TRANSFIELD FORCES ACTING ON A FIELDLINE In the non-relativistic case, it is customary to divide the force acting on the field line in four contributions: the kinetic force, the thermal pressure force, the electromagnetic force, and the gravitational force. The net force is zero, so we have the following equation: $$\begin{eqnarray*} \mathcal {F}_{\rm K} + \mathcal {F}_{\rm T} + \mathcal {F}_{\rm EM} + \mathcal {F}_{\rm G} \equiv \rho (\boldsymbol{V} \cdot \nabla)\boldsymbol{V} + \nabla P -\frac{(\nabla \times \boldsymbol{B})\times \boldsymbol{B}}{4\pi } - \rho \nabla \frac{\mathcal {GM}}{r} = 0. \end{eqnarray*}$$(C1) By taking the inner product with |$\hat{\boldsymbol{n}}$|⁠, we obtain the projection of the forces in the direction perpendicular to the field line, whereas with the inner product with |$\hat{\boldsymbol{b}}$|⁠, we obtain the forces along the poloidal direction. C1 Transfield forces The transfield kinetic force: $$\begin{eqnarray*} \mathcal {F}_{\rm K, \perp } &\equiv & \rho (\boldsymbol{V} \cdot \nabla)\boldsymbol{V} \cdot \hat{\boldsymbol{n}} = \, \rho \cos (\theta +\psi)\left[V_r\frac{\partial V_r}{\partial r} + \frac{V_\theta }{r}\frac{\partial V_r}{\partial \theta } - \frac{V_\theta ^2+V_\phi ^2}{r}\right] -\rho \sin (\theta +\psi)\left[V_r\frac{\partial V_\theta }{\partial r} + \frac{V_\theta }{r}\frac{\partial V_\theta }{\partial \theta } + \frac{V_\theta V_r}{r} - \frac{V_\phi ^2}{r\tan (\theta)}\right] \nonumber \\ &= &\, \rho \frac{V_p^2}{\varpi }\sin (\theta)\cos (\theta +\psi)\frac{\mathrm{ d}\psi }{\mathrm{ d}\theta } +\rho \frac{V_\phi ^2}{\varpi }\sin (\psi)\nonumber \\ &= &\, \rho \left(- \frac{V_\star \alpha ^{-1/4}M^2}{G^2}\frac{\sin (\theta)}{\cos (\theta +\psi)}\right)^2\sin (\theta)\frac{\cos (\theta +\psi)}{\varpi }\frac{\mathrm{ d}\psi }{\mathrm{ d}\theta } + \rho \left(\lambda _{\rm VTST}\frac{V_\star \alpha ^{-1/4}}{G}\frac{(G^2-M^2)}{(1-M^2)}\right)^2\frac{\sin (\psi)}{\varpi }\nonumber \\ &= & \left[\frac{B_\star ^2\alpha ^{F-2}}{4\pi \varpi G^4}\frac{\sin (\theta)}{\cos (\theta +\psi)}\right] \left\lbrace M^2\sin ^2(\theta)\frac{\mathrm{ d}\psi }{\mathrm{ d}\theta }+\lambda _{\rm VTST}^2\frac{G^2}{M^2}\frac{(G^2-M^2)^2}{(1-M^2)^2}\frac{\sin (\psi)}{\sin (\theta)}\cos (\theta +\psi)\right\rbrace . \end{eqnarray*}$$(C2) The transfield thermal pressure force: $$\begin{eqnarray*} \mathcal {F}_{\rm T, \perp } &\equiv & \nabla P \cdot \hat{\boldsymbol{n}} = \left\lbrace \frac{\partial P}{\partial r} \hat{\boldsymbol{r}} + \frac{1}{r} \frac{\partial P}{\partial \theta } \hat{\boldsymbol{\theta }}\right\rbrace \cdot \hat{\boldsymbol{n}} \nonumber \\ &=&\, 2P \frac{\sin (\theta)}{\varpi \cos (\theta +\psi)}\left((F-2)(\cos ^2(\theta +\psi)+\sin ^2(\theta +\psi))+\frac{\Gamma }{2M^2}\sin (\theta +\psi) \cos (\theta +\psi)\frac{\mathrm{ d}M^2}{\mathrm{ d}\theta }\right)\nonumber \\ &=& \left[\frac{B_\star ^2\alpha ^{F-2}}{4\pi \varpi G^4}\frac{\sin (\theta)}{\cos (\theta +\psi)}\right] \left(\frac{G^4\mu _{\rm VTST}}{M^{2\Gamma }}\right)\left\lbrace (F-2)+\frac{\Gamma }{2M^2}\sin (\theta +\psi) \cos (\theta +\psi)\frac{\mathrm{ d}M^2}{\mathrm{ d}\theta }\right\rbrace . \end{eqnarray*}$$(C3) The transfield electromagnetic force: $$\begin{eqnarray*} \mathcal {F}_{\rm EM, \perp } &\equiv &- \frac{(\nabla \times \boldsymbol{B})\times \boldsymbol{B}}{4\pi } \cdot \hat{\boldsymbol{n}} = - \frac{1}{4\pi } \left[\frac{B_\phi }{r}\left(\frac{1}{\sin (\theta)}\frac{\partial B_r}{\partial \phi }-\frac{\partial (rB_\phi)}{\partial r}\right) - \frac{B_\theta }{r}\left(\frac{\partial (rB_\theta)}{\partial r} -\frac{\partial B_r}{\partial \theta }\right)\right] \hat{\boldsymbol{r}}\cdot \hat{\boldsymbol{n}}\nonumber \\ && -\frac{1}{4\pi } \left[ \frac{B_r}{r}\left(\frac{\partial (rB_\theta)}{\partial r} -\frac{\partial B_r}{\partial \theta }\right) - \frac{B_\phi }{r\sin (\theta)}\left(\frac{\partial (\sin (\theta)B_\phi)}{\partial \theta } - \frac{\partial B_\theta }{\partial \phi } \right)\right] \hat{\boldsymbol{\theta }}\cdot \hat{\boldsymbol{n}}\nonumber \\ &= & \frac{B_\phi ^2}{4\pi \varpi }\frac{\sin (\theta)}{\cos (\theta +\psi)}\left\lbrace (F-1)+ \frac{2\cos (\psi)\sin (\theta +\psi)}{\sin (\theta)}\frac{G^2}{1-G^2} - \frac{\sin (\theta +\psi)\cos (\theta +\psi) }{(1-M^2)}\frac{\mathrm{ d}M^2}{\mathrm{ d}\theta }\right\rbrace \nonumber \\ && +\frac{B_p^2}{4\pi \varpi }\frac{\sin (\theta)}{\cos (\theta +\psi)}\left\lbrace (F-2)+\frac{\cos (\psi)\sin (\theta +\psi)}{\sin (\theta)}-\frac{\mathrm{ d}\psi }{\mathrm{ d}\theta } \right\rbrace \nonumber \\ &= & \left[\frac{B_\star ^2\alpha ^{F-2}}{4\pi \varpi G^4}\frac{\sin (\theta)}{\cos (\theta +\psi)}\right]\lambda _{\rm VTST}^2\left(\frac{G(1-G^2)}{(1-M^2)}\right)^2\left\lbrace (F-1)+ \frac{2\cos (\psi)\sin (\theta +\psi)}{\sin (\theta)}\frac{G^2}{1-G^2} - \frac{\sin (\theta +\psi)\cos (\theta +\psi) }{(1-M^2)}\frac{\mathrm{ d}M^2}{\mathrm{ d}\theta }\right\rbrace \nonumber \\ && +\left[\frac{B_\star ^2\alpha ^{F-2}}{4\pi \varpi G^4}\frac{\sin (\theta)}{\cos (\theta +\psi)}\right] \frac{\sin ^2(\theta)}{\cos ^2(\theta +\psi)}\left\lbrace (F-2)+\frac{\cos (\psi)\sin (\theta +\psi)}{\sin (\theta)}-\frac{\mathrm{ d}\psi }{\mathrm{ d}\theta } \right\rbrace . \end{eqnarray*}$$(C4) The transfield gravitational force: $$\begin{eqnarray*} \mathcal {F}_{\rm G, \perp } &\equiv & -\rho \nabla \frac{\mathcal {GM}}{r} \cdot \hat{\boldsymbol{n}} = - \frac{\rho _\star }{M^2}\alpha ^{F-3/2} \mathcal {GM}\frac{\partial }{\partial r}\left(\frac{1}{r}\right) \hat{\boldsymbol{r}}\cdot \hat{\boldsymbol{n}}\nonumber \\ &= & \frac{\rho _\star }{M^2}\alpha ^{F-3/2} \mathcal {GM}\frac{\sin ^2(\theta)}{\varpi ^2}\cos (\theta +\psi)\nonumber \\ &= & \left[\frac{B_\star ^2\alpha ^{F-2}}{4\pi \varpi G^4}\frac{\sin (\theta)}{\cos (\theta +\psi)}\right] \left(\frac{k_{\rm VTST}^2\sin (\theta)}{G}\right)\frac{G^4}{M^2}\cos ^2(\theta +\psi) . \end{eqnarray*}$$(C5) C2 Poloidal forces In this section, we give the forces along the poloidal directions. The poloidal kinetic force: $$\begin{eqnarray*} \mathcal {F}_{\rm K, \parallel } &\equiv &\rho (\boldsymbol{V} \cdot \nabla)\boldsymbol{V} \cdot \hat{\boldsymbol{b}} = \rho \sin (\theta +\psi)\left[V_r\frac{\partial V_r}{\partial r} + \frac{V_\theta }{r}\frac{\partial V_r}{\partial \theta } - \frac{V_\theta ^2+V_\phi ^2}{r}\right] +\rho \cos (\theta +\psi)\left[V_r\frac{\partial V_\theta }{\partial r} + \frac{V_\theta }{r}\frac{\partial V_\theta }{\partial \theta } + \frac{V_\theta V_r}{r} - \frac{V_\phi ^2}{r\tan (\theta)}\right] \nonumber \\ &=& \rho \frac{V_p^2}{r}\frac{\cos (\theta +\psi)}{M^2} \frac{\mathrm{ d}M^2}{\mathrm{ d}\theta } +\rho \frac{V_p^2}{r}\sin (\theta +\psi)\frac{\mathrm{ d}\psi }{\mathrm{ d}\theta } -\rho \frac{V_p^2 + V_\phi ^2}{r} \frac{\cos (\psi)}{\sin (\theta)}\nonumber \\ &=& \left(\frac{B_\star ^2}{4\pi V_\star ^2}\frac{\alpha ^{F-3/2}}{M^2}\right)\frac{\sin (\theta)}{\varpi } \left\lbrace \left(- \frac{V_\star \alpha ^{-1/4}M^2}{G^2}\frac{\sin (\theta)}{\cos (\theta +\psi)}\right)^2 \left[\frac{\cos (\theta +\psi)}{M^2} \frac{\mathrm{ d}M^2}{\mathrm{ d}\theta } +\sin (\theta +\psi)\frac{\mathrm{ d}\psi }{\mathrm{ d}\theta } -\frac{\cos (\psi)}{\sin (\theta)}\right] \right.\nonumber \\ && \left. -\left(\lambda _{\rm VTST}\frac{V_\star \alpha ^{-1/4}}{G}\frac{(G^2-M^2)}{(1-M^2)}\right)^2 \frac{\cos (\psi)}{\sin (\theta)}\right\rbrace \nonumber \\ &=& \left[\frac{B_\star ^2\alpha ^{F-2}}{4\pi \varpi G^4}\frac{\sin (\theta)}{\cos (\theta +\psi)}\right]\left\lbrace \sin ^2(\theta)\frac{\mathrm{ d}M^2}{\mathrm{ d}\theta } + M^2\sin ^2(\theta)\tan (\theta +\psi)\frac{\mathrm{ d}\psi }{\mathrm{ d}\theta } - M^2\frac{\sin (\theta)\cos (\psi)}{\cos (\theta +\psi)} \right.\nonumber \\ && \left.- \lambda _{\rm VTST}^2\frac{G^2}{M^2}\frac{(G^2 -M^2)^2}{(1-M^2)^2}\frac{\cos (\psi)\cos (\theta +\psi)}{\sin (\theta)}\right\rbrace . \end{eqnarray*}$$(C6) The thermal pressure force: $$\begin{eqnarray*} \mathcal {F}_{\rm T, \parallel } & \equiv & \nabla P \cdot \hat{\boldsymbol{b}} = \left\lbrace \frac{\partial P}{\partial r} \hat{\boldsymbol{r}} + \frac{1}{r} \frac{\partial P}{\partial \theta } \hat{\boldsymbol{\theta }}\right\rbrace \cdot \hat{\boldsymbol{b}} \nonumber \\ &=& - P \frac{\sin (\theta)}{\varpi \cos (\theta +\psi)}\left(\cos ^2(\theta +\psi)\frac{\Gamma }{M^2}\frac{\mathrm{ d}M^2}{\mathrm{ d}\theta }\right)\nonumber \\ &=& - \left[\frac{B_\star ^2\alpha ^{F-2}}{4\pi \varpi G^4}\frac{\sin (\theta)}{\cos (\theta +\psi)}\right] \left(\frac{G^4\mu _{\rm VTST}}{M^{2\Gamma }}\right)\cos ^2(\theta +\psi)\frac{\Gamma }{2M^2}\frac{\mathrm{ d}M^2}{\mathrm{ d}\theta }. \end{eqnarray*}$$(C7) The electromagnetic force: $$\begin{eqnarray*} \mathcal {F}_{\rm EM, \parallel } & \equiv & - \frac{(\nabla \times \boldsymbol{B})\times \boldsymbol{B}}{4\pi } \cdot \hat{\boldsymbol{b}} = - \frac{1}{4\pi } \left[\frac{B_\phi }{r}\left(\frac{1}{\sin (\theta)}\frac{\partial B_r}{\partial \phi }-\frac{\partial (rB_\phi)}{\partial r}\right) - \frac{B_\theta }{r}\left(\frac{\partial (rB_\theta)}{\partial r} -\frac{\partial B_r}{\partial \theta }\right)\right] \hat{\boldsymbol{r}}\cdot \hat{\boldsymbol{b}}\nonumber \\ && - \frac{1}{4\pi } \left[ \frac{B_r}{r}\left(\frac{\partial (rB_\theta)}{\partial r} -\frac{\partial B_r}{\partial \theta }\right) - \frac{B_\phi }{r\sin (\theta)}\left(\frac{\partial (\sin (\theta)B_\phi)}{\partial \theta } - \frac{\partial B_\theta }{\partial \phi } \right)\right] \hat{\boldsymbol{\theta }}\cdot \hat{\boldsymbol{b}}\nonumber \\ &=& - \frac{1}{4\pi }\frac{B_\phi ^2}{r}\left[\frac{2G^2}{1-G^2}\frac{\cos (\psi)}{\sin (\theta)} - \frac{\cos (\theta + \psi)}{1-M^2}\frac{\mathrm{ d}M^2}{\mathrm{ d}\theta }\right]\nonumber \\ &=& - \left[\frac{B_\star ^2\alpha ^{F-2}}{4\pi \varpi G^4}\frac{\sin (\theta)}{\cos (\theta +\psi)}\right] \lambda _{\rm VTST}^2\frac{G^2(1-G^2)^2}{(1-M^2)^2}\left[\frac{2G^2}{1-G^2}\frac{\cos (\psi)\cos (\theta +\psi)}{\sin (\theta)} - \frac{\cos ^2(\theta + \psi)}{1-M^2}\frac{\mathrm{ d}M^2}{\mathrm{ d}\theta }\right] . \end{eqnarray*}$$(C8) The gravitational force: $$\begin{eqnarray*} \mathcal {F}_{\rm G, \parallel } & \equiv & - \rho \nabla \frac{\mathcal {GM}}{r} \cdot \hat{\boldsymbol{b}} = - \frac{\rho _\star }{M^2}\alpha ^{F-3/2} \mathcal {GM}\frac{\partial }{\partial r}\left(\frac{1}{r}\right) \hat{\boldsymbol{r}}\cdot \hat{\boldsymbol{b}}\nonumber \\ &= & \frac{\rho _\star }{M^2}\alpha ^{F-3/2} \mathcal {GM}\frac{\sin ^2(\theta)}{\varpi ^2}\sin (\theta +\psi)\nonumber \\ &= & \left[\frac{B_\star ^2\alpha ^{F-2}}{4\pi \varpi G^4}\frac{\sin (\theta)}{\cos (\theta +\psi)}\right] \left(\frac{k_{\rm VTST}^2\sin (\theta)}{G}\right)\frac{G^4}{M^2}\sin (\theta +\psi)\cos (\theta +\psi). \end{eqnarray*}$$(C9) © 2020 The Author(s) Published by Oxford University Press on behalf of Royal Astronomical Society This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://academic.oup.com/journals/pages/open_access/funder_policies/chorus/standard_publication_model)
TI - A study of radial self-similar non-relativistic MHD outflow models: parameter space exploration and application to the water fountain W43A
JO - Monthly Notices of the Royal Astronomical Society
DO - 10.1093/mnras/staa3660
DA - 2021-01-05
UR - https://www.deepdyve.com/lp/oxford-university-press/a-study-of-radial-self-similar-non-relativistic-mhd-outflow-models-kQpCa0c8PS
SP - 2071
EP - 2090
VL - 501
IS - 2
DP - DeepDyve
ER -