TY - JOUR
AU - Xia,, Yu-Kai
AB - Abstract We study magnetohydrodynamic (MHD) self-similar collapses and void evolution, with or without shocks, of a general polytropic quasi-spherical magnetofluid permeated by random transverse magnetic fields under the Paczynski–Wiita gravity that captures essential general relativistic effects of a Schwarzschild black hole (BH) with a growing mass. Based on the derived set of non-linear MHD ordinary differential equations, we obtain various asymptotic MHD solutions, the geometric and analytical properties of the magnetosonic critical curve (MSCC) and MHD shock jump conditions. Novel asymptotic MHD solution behaviours near the rim of central expanding voids are derived analytically. By exploring numerical global MHD solutions, we identify allowable boundary conditions at large radii that accommodate a smooth solution and show that a reasonable amount of magnetization significantly increases the mass accretion rate in the expansion-wave-collapse solution scenario. We also construct the counterparts of envelope-expansion-core-collapse solutions that cross the MSCC twice, which are found to be closely paired with a sequence of global smooth solutions satisfying a novel type of central MHD behaviours. MHD shocks with static outer and various inner flow profiles are also examined. Astrophysical applications include dynamic core collapses of magnetized massive stars and compact objects as well as formation of supermassive, hypermassive, dark matter and mixed matter BHs in the Universe, including the early Universe. Such gigantic BHs can be detected in X-ray/gamma-ray sources, quasars, ultraluminous infrared galaxies or extremely luminous infrared galaxies and dark matter overwhelmingly dominated elliptical galaxies as well as massive dark matter halos, etc. Gravitational waves and electromagnetic wave emissions in broad band (including e.g., gamma-ray bursts and fast radio bursts) can result from this type of dynamic collapses of forming BHs involving magnetized media. accretion, accretion discs, black hole physics, magnetic fields, MHD, supernovae: general, quasars: supermassive black holes 1 INTRODUCTION The self-similar dynamic evolution of a spherical gas flow under self-gravity has long been an important subject of extensive theoretical research in astrophysics, primarily because of its relative simplicity in analytical and numerical treatments as well as broad diversity of astrophysical applications. Self-similarity, widely known as a universal manifestation of dynamic evolution sufficiently far from initial and/or boundary conditions, is theoretically backed up by simple yet powerful arguments such as the Buckingham Pi theorem of dimensional analysis (see, e.g. Buckingham 1914; Gibbings 2011). Numerous works have been published for self-similar dynamic solutions of partial differential equations (PDEs) arising from a wide range of hydrodynamic problems (e.g. Taylor 1950; Sedov 1959; Barenblatt 1980; Velikovich, Liberman & Shmalts 1985). More recently, evolutionary convergence to self-similar solutions have been found in numerical simulations for various astrophysical settings both in terms of flow profiles and predictions for relevant physical quantities (e.g. Foster & Chevalier 1993; Harada, Maeda & Semelin 2003; Masuyama, Shigeyama & Tsuboki 2016), adding to more confidence in the reliability and applicability of self-similar dynamic solutions. Meanwhile, we emphasize a completely new perspective when a flow system possesses several distinct classes of different self-similar dynamic solutions. This also bears the relevance to, e.g., self-similar dynamic solutions of Boltzmann equation. In other words, a system may have several possible tracks of self-similar dynamic evolution sufficiently far away from the initial and/or boundary conditions. This is a challenging problem and numerical simulations combined with analytical analyses deem necessary for further progress. The most important point is that a flow system somehow still remembers different types of initial and/or boundary conditions or their various feasible combinations for possible evolution paths. In this study, we have also made a preliminary analysis on this front. One important astrophysical application for the study of spherical self-similar flow under self-gravity is the formation of protostars, brown dwarfs and gaseous planets. It was pioneered by Larson (1969) and Penston (1969) who independently studied a self-similar solution describing the central core of a collapsing isothermal gas sphere under self-gravity. Shu (1977) proposed a competing ‘inside-out collapse’ scenario based on his isothermal self-similar expansion-wave-collapse solution (EWCS). Further analyses were subsequently conducted by Hunter (1977), Whitworth & Summers (1985) and Tsai & Hsu (1995) for an isothermal gas sphere and Suto & Silk (1988) for a conventional polytropic (CP) gas sphere. The phase-matching technique for seeking special smooth solutions across the sonic critical curve (SCC) developed in the above works were fully utilized by Lou & Shen (2004) to construct envelope-expansion-core-collapse (EECC) solutions that may model the asymptotic giant branch (AGB) phase and/or the dynamic phase of planetary nebulae (PNe) in the late evolution of low-mass stars (e.g. Lou & Zhai 2009, 2010). Protostar formation in molecular cloud (MC) cores was pursued further under a more versatile self-similar dynamic model of polytropic gas by Lou & Gao (2006), Fu, Gao & Lou (2011), Lou & Gao (2011), and Gao & Lou (2017) involving self-similar shocks in dynamic gas flow around star-forming MC cores for possibly variable mass accretion rates to avoid the so-called ‘luminosity problem’. Isothermal and general polytropic (GP) self-similar dynamic ‘Champagne flows’ in H ii regions were modelled by Shu et al. (2002) and Hu & Lou (2008), respectively. More recently, GP counterparts of Larson–Penston and Hunter-type solutions were invoked in Lou & Shi (2014b) to study the envelope mass infall rates (EMIRs) of dense cores in MCs in contexts of forming protostars, brown dwarfs and gaseous planets in mini gas globules (MMGs) – smaller versions of Bok globules. Being possibly evolving in isolation, they may form a considerable part of the so-called free floating objects - a fact that challenges the disc formation theory as their sole origin. Besides formation processes of protostars, brown dwarfs and gaseous planets, other astrophysical phenomena addressed by such models include progenitor core collapses prior to supernova explosions and supernova remnants, often involving converging or rebound shocks, interaction of ejecta with surrounding materials and various radiative processes (e.g. Chevalier 1982; Yahil 1983; Chevalier, Blondin & Emmering 1992; Lou 1994; Toptygin 2000; Lou & Wang 2006, 2007; Lou & Cao 2008; Lou & Shi 2014a); formation and dynamic evolution of cosmic voids and galaxy clusters with flows of dark matter (DM) and magnetized intracluster medium (e.g. Fillmore & Goldreich 1984a,b; Bertschinger 1985a,b; Lou, Jiang & Jin 2008); as well as the spherical symmetric dynamic mass accretion of BHs as an important alternative scenario to commonly believed disc accretion (Lou & Jiang 2008; Lou & Lian 2012; Lou & Wu 2012; Lian & Lou 2014; Lou & Hu 2016, Lou 2015a). Specifically, Lou & Jiang (2008) offered a theoretical explanation for the empirical MBH–σ relation between supermassive black holes (SMBHs) and host galactic bulges and elliptical galaxies based on a quasi-static self-similar dynamic evolution scenario, while Lou & Wu (2012) further extended this analysis of MBH–σ relation for globular clusters with possible central intermediate-mass black holes (IMBHs). Such IMBHs and SMBHs may also reside at the centres of dwarf spheroidal galaxies. In parallel, DM black holes (DMBHs) may also form within massive DM haloes and DM overwhelmingly dominated elliptical galaxies as suggested by Lou (2015a), Lou & Hu (2016), and Shi & Lou (2017). Newtonian gravity has been extensively used in most of the foregoing research references. However, when highly compact collapsed objects such as neutron stars or (massive) BHs are involved at the centre, Newtonian model results are regarded as speculative and heuristic from the rigorous general relativistic (GR) perspective. A fully GR treatment would be certainly involved mathematically and has been attempted semi-analytically in Harada & Maeda (2001) and Cai & Shu (2005) for a dynamically collapsing isothermal sphere. To make further GP extensions while avoiding pitfalls of the Newtonian gravity in the strong-field regime, Lian & Lou (2014) adopted the Paczynski–Wiita gravity in their GP hydrodynamic formulation as a next-order correction to the Newtonian gravity for retaining essential GR effects and constructed self-similar dynamic collapse solutions for GP spheres. The Paczynski-Wiita (1980) potential was first introduced as an expedient instrument of keeping GR effects on accretion disc dynamics around non-rotating Schwarzschild BHs and has proven capable of accurately reproducing several key features of the space–time near a Schwarzschild BH (Abramowicz 2009). Using the Paczynski-Wiita gravity, Lian & Lou (2014) obtained generalized EWCSs with a central collapsed Schwarzschild BH within an expanding event horizon. This gives rise to a highly efficient mechanism of dynamically growing SMBHs (with a mass range of ∼106–109 M⊙) and also hypermassive black holes (HMBHs with a mass range of ∼1010–1012M⊙ and even higher mass) in the Universe, including the early Universe (see also Lou 2015a; Lou & Hu 2016 for more information on observational evidence and theoretical aspects). Due to the inferred ubiquitous presence of DM haloes as DM mass reservoirs on various large scales and their gravitational interactions with normal matter (NM), it would be also possible to form DM black holes (DMBHs) and mixed matter black holes (MMBHs) in the Universe, including the early Universe as specifically highlighted in Lou (2015a) via gravitational instability analysis of Bonnor–Ebert sphere. For extremely massive DM haloes, general relativistic instability (GRI) may induce DM gravitational collapses to form SMBHs/HMBHs in the Universe including the early Universe (Shi & Lou 2017). Empirically, HMBHs and SMBHs can be found in X-ray/gamma-ray sources, quasars, ultraluminous infrared galaxies (ULIRGs) and extremely luminous infrared galaxies (ELIRGs), DM overwhelmingly dominated elliptical galaxies, dwarf spheroidals and so forth (Lou & Hu 2016). Naturally, this requires the presence of sufficiently large host (DM) mass reservoirs or haloes with sustained dynamic collapses and/or accretions without interruptions in the Universe, including the early Universe. In particular, unambiguous theoretical predictions of Lian & Lou (2014) for possible HMBHs in the Universe, including the early Universe, were quickly borne out by observations in recent years around the predicted lower mass end of HMBHs (e.g. Wang et al. 2015; Wu et al. 2015; Lou & Hu 2016; Oesch et al. 2016). With the scheduled launch of 6.5 m James Webb Space Telescope (JWST) in October 2018, more exciting discoveries and detections of HMBHs/SMBHs in the early Universe (say, ∼z = 5 to z ≳ 12) are anticipated by JWST superior infrared detection capability (Mather, private communications). Along this line of physical reasoning, it is natural to conceive a scenario in which magnetized gas with a certain fraction, be it large or small, may be involved in such large-scale dynamic collapses of forming HMBHs and SMBHs. At various levels, magnetic field provides dynamic, diagnostic as well as collateral effects including MHD dynamo, reconnection, heating, acceleration of cosmic ray particles (including ultra-high energy cosmic rays - UHECRs), generating electromagnetic radiations in a wide frequency band (including possibly gamma-ray bursts - GRBs and fast radio bursts - FRBs and so forth) and inducing/creating massive MHD turbulence for gravitational wave emissions. The host mass reservoirs would be expected to have higher concentrations of neutral hydrogens (compared to the mean level of H i around the pertinent epoch in the Universe) sufficiently far away from central activities with the characteristic redshifted 21 cm line radiation. From the perspective of a practical observational strategy, the search for such higher concentrations of neutral hydrogens in gigantic mass reservoirs in the redshift range of z ∼ 5 to ≳ 12 would be an important prerequisite and a sensible exercise for detecting diffuse neutral hydrogens of even higher redshifts by radio telescopes. In other words, for the early Universe with the redshift range of z ∼ 5 to ≳ 12, the spatial distribution of neutral hydrogen is most likely already quite clumpy in association with various mass reservoirs, large and small. In this paper, we shall adopt the Paczynski–Wiita gravity to retain key GR effects in a quasi-spherical GP magnetohydrodynamic (MHD) formalism involving random transverse magnetic fields (Wang & Lou 2007, 2008; Lian & Lou 2014) and mainly focus on the dynamic effects of magnetic pressure and tension forces. For galaxies and clusters of galaxies, more massive DM haloes compared to directly observed NM materials have been widely inferred. There are also other astrophysical systems known most likely with much more dominant DM over NM (e.g. Segue 1; see Xiang-Grüss et al. 2009a,b). Logically, there are no physical reasons at all against possible DM haloes large and small to exist by their own in the Universe, including the early Universe, as we have assigned or transferred all gravitational properties of NM to DM (Lou 2015a). It is just a matter of how to detect them or to infer their presence more reliably and to dig them out more effectively. In the absence of direct gravitational interactions with NM, it would be extremely difficult to detect such DM haloes though. Conceptually and under fortuitous situations (still extraordinarily challenging), one may utilize gravitational lensing effects of electromagnetic waves and of gravitational waves from background sources to infer such foreground massive DM haloes. It is also conceivable to have magnetized gas of whatever amount trapped in the gravitational potential of such more massive DM haloes (Lou & Wu 2005; Wu & Lou 2006; Xiang-Grüss et al. 2009a,b; Lou 2015a). For such a situation, we would have an additional valuable diagnostic tool associated with the magnetic field effects for detecting such massive DM haloes. Referring to several very recent publications of van Dokkum et al. (2016 and pertinent references therein), the so-called ultradiffuse galaxy (UDG) Dragonfly 44 (DF 44) detected by the Dragonfly Telephoto Array in New Mexico is the largest (and second-brightest) among 47 UDGs found so far in the Coma cluster of galaxies – a relatively nearby (z = 0.023 and thus a distance of ∼100 Mpc) rich cluster of more than one thousand galaxies. According to the Scientific American summary report by Weitering (2016), what is extremely unusual (or ‘weird’) about UDG DF 44 is that DF 44 contains 99.99 per cent of DM and only one-hundredth of one per cent of this UDG is NM as observed. To us, this would not be surprising (Lou & Wu 2005; Wu & Lou 2006; Xiang-Grüss et al. 2009a,b; Lou 2015a) and, in fact, was more than welcome as we have anticipated this type of dominant DM halo systems all along; as examples, these UDGs form a new class of massive astrophysical halo systems with much more overwhelming DM domination over NM. It is natural to further expect that DM haloes of even less NM fraction in Coma cluster and/or other clusters of galaxies as well as elsewhere of relative isolation in the Universe, including the early Universe. Consequently, we would venture to further predict possible SMBHs and HMBHs at the centres of these massive UDGs and the likes (Lian & Lou 2014; Lou 2015a; Lou & Hu 2016). These are DMBHs and/or MMBHs with very small fractions of NM materials involved in their dynamic formation. In addition to optical signatures, the presence of neutral hydrogens and magnetized gas medium would very much increase the chance of such detections of massive DM and MM haloes by Karl G. Jansky Very Large Array (VLA), the Low-Frequency Array (LOFAR), the Square Kilometre Array (SKA), the Five-hundred-metre Aperture Spherical radio Telescope (FAST) and so forth. The hydrodynamic formalism of Lian & Lou (2014) would be made substantially more versatile by including random transverse magnetic fields, which are ubiquitous in various astrophysical settings (e.g. from stars to clusters of galaxies). Such a magnetic field may be some fossil field, or sustained by convective or turbulent motions in fully or partially ionized plasmas via MHD dynamo mechanism (e.g. Parker 1979; Thompson & Duncan 1993), and possibly enhanced during subsequent collapsing or merging processes (e.g. Braithwaite & Spruit 2004; Ferrario & Wickramasinghe 2006). For instance, magnetic field is known to stabilize MCs against gravitational collapse to some degree and may affect the distribution of stellar masses (e.g. Shu, Adams & Lizano 1987; Lou 1996; Myers 1998) and MHD shocks in supernova remnants serve as the accelerators for high-energy cosmic ray particles (e.g. Blandford & Eichler 1987; Berezhko et al. 1988; Lou 1994; Ackermann et al. 2013). Theoretical works have been envisaging the crucial role of magnetic fields near BHs in the removal of angular momentum of accretion discs through magnetorotational instabilities (MRIs) (e.g. Velikhov 1959; Chandrasekhar 1961; Hawley & Balbus 1991), as well as in extracting the rotational energy of an accretion disc and producing collimated relativistic jets of various scales (e.g. Blandford & Znajek 1977). Recent observations have discovered dynamically important magnetic fields both near the SMBH in our own Galaxy centre and in remote active galactic nuclei (AGNs) (Eatough et al. 2013; Martí-Vidal et al. 2015). On larger scales, the existence of randomly oriented magnetic fields at μG level permeating galaxy clusters is now well established (e.g. Jaffe 1980; Zweibel & Heiles 1997; Govoni & Feretti 2004). In a grossly spherical, self-gravitating and highly conducting gas system, the dynamical behaviour of the magnetic field is dictated by the magnetic flux conservation, commonly referred to as the magnetic field ‘frozen in’ condition. When the magnetic field has random orientation and rotation is sufficiently slow, an overall ‘quasi-spherical’ (e.g. Zeldovich & Novikov 1971) geometry may emerge and maintain itself and is thus suitable for our MHD model treatment. In such a geometry and on large scales, the magnetic field manifests through average magnetic pressure and magnetic tension force in the radial direction. The quasi-spherical MHD problem was considered by Chiueh & Chou (1994) by including only the magnetic pressure gradient term in the radial momentum equation but with a very limited solution exploration. This MHD formalism was further developed by Yu et al. (2006) into the notion of ‘random transverse magnetic fields’, i.e. a ‘ball of thread’ scenario that considers the large-scale effect of magnetic fields on a spatially averaged basis. Both magnetic pressure and tension forces are included. In reference to this scenario, Wang & Lou (2007) and Wang & Lou (2008) investigated CP and GP MHD shocks, whose results were subsequently applied to magnetized supernova explosions (Lou & Wang 2007) and to the formation of magnetars via strongly magnetized core collapses of massive progenitor stars (Hu & Lou 2009). This paper aims at the quasi-spherical self-similar GP MHD problem under the combined effects of the Paczynski–Wiita gravity and random transverse magnetic fields. The magnetofluid equation of state (EoS) is a GP one, i.e. the ratio p/ργ is formally conserved along streamlines but needs not to be a global constant in space and time (see Wang & Lou 2008). We classify the possible asymptotic MHD solution behaviours and systematically present global numerical solutions that do or do not cross the magnetosonic critical curve (MSCC). Sample solutions involving MHD shocks are shown, analysed and discussed. In addition to the MHD counterparts of solutions that have been previously found in the unmagnetized case, our solutions show the influence of magnetic field and exhibit several novel features. Specifically, we can achieve a much higher mass accretion rate towards the central BH when the gas is magnetized in the EWCS scenario and new asymptotic behaviours near the central void emerge for a non-zero magnetic field. Moreover, we report a new type of central asymptotic solution for γ < 4/3 that is unknown before, yet plays an important role in forming the sequence of the discrete ‘pair solutions’ that we find in our numerical MHD solutions involving MSCC crossing or shock waves. The structure of the paper is outlined below: Section 2 describes the basic formulation of our non-linear MHD model and the suitable MHD self-similar transformation leading to two non-linear coupled ordinary differential equations (ODEs). In Section 3, we study the possible analytic asymptotic MHD solutions for large and small x, as well as for x near the event horizon, MSCC and central void rim. Theoretical aspects concerning the construction of MHD shock solutions are the main focus of Section 4, while Section 5 is devoted to geometric and analytical properties of the MSCC. In Section 6, sample numerical solutions are presented and their underlying physical implications are discussed. Finally, we provide a summary of our findings in Section 7. For the convenience of reference, several relevant mathematical details are included in appendices. 2 THEORETICAL FORMULATION OF THE PHYSICAL PROBLEM 2.1 Non-linear MHD partial differential equations We consider a quasi-spherical, self-gravitating, GP ideal magnetofluid with random transverse magnetic fields. Such a GP magnetofluid was formulated and studied by Wang & Lou (2008) with the Newtonian gravity. To account for possible GR effects while keeping the problem formulation readily tractable, we here adopt the Paczynski–Wiita gravity (Paczyński & Wiita 1980; Abramowicz 2009) and keep the rest of the MHD model formulation essentially the same as in Wang & Lou (2008). Paczynski–Wiita gravity has been invoked earlier in Lian & Lou (2014) to study the spherically symmetric GP hydrodynamic collapses of forming Schwarzschild BHs. For such spherical dynamic collapses sustained in sufficiently large mass reservoirs, SMBHS (mass range ∼106–109 M⊙) and HMBHs (mass range ∼1010–1012 M⊙ or higher) can form sufficiently fast in the entire Universe, including the early Universe, with supporting observational evidence (Wang et al. 2015; Wu et al. 2015; Oesch et al. 2016). Several months ago, the Hubble Space Telescope (HST) detected a galaxy of mass ∼109 M⊙ at a record high redshift z = 11.1 corresponding to an early epoch of ∼400 million years since the Big Bang. Recently, fluffy UDG DF44 located in the rich Coma Cluster was inferred to possess a DM halo about ∼10 000 times the visible NM of DF44. It would not be surprising for DF44 DM halo to host a central SMBH or HMBH (Lou & Wu 2005; Wu & Lou 2006; Lian & Lou 2014; Lou 2015a; Lou & Hu 2016). Such phenomena may well happen also in the early Universe within gigantic DM mass reservoirs or haloes. For the Paczynski–Wiita approximation, the gravity experienced by a unit point mass, located at a radius r from the centre of a spherical symmetric gravitating mass M, is \begin{equation} g = -\frac{GM}{\left(r-2GM/c^2\right)^2}, \end{equation} (1) where G = 6.672 59 × 10−8 dyn cm2 g−2 is the universal gravitational constant and c = 3 × 1010 cm s−1 is the speed of light in vacuum. The term 2GM/c2 here is recognized as the Schwarzschild radius of the gravitating sphere mass M and so the modification is easy to remember. Expression (1) for the Paczynski–Wiita gravity is only valid for r > 2GM/c2. If M happens to be a non-rotating central BH, then expression (1) can only be applied outside the event horizon (Lian & Lou 2014). We adopt spherical polar coordinates (r, θ, ϕ), with the origin at the centre of a quasi-spherical magnetofluid. The basic MHD PDEs for the mass conservation are \begin{eqnarray} \frac{\mathrm{\partial} \rho }{\mathrm{\partial} t} + \frac{1}{r^2}\frac{\mathrm{\partial} }{\mathrm{\partial} r}(r^2\rho u) = 0, \end{eqnarray} (2) \begin{eqnarray} \left(\frac{\mathrm{\partial} }{\mathrm{\partial} t}+u\frac{\mathrm{\partial} }{\mathrm{\partial} r}\right)M = 0\ , \qquad \frac{\mathrm{\partial} M}{\mathrm{\partial} r} = 4\pi r^2 \rho , \end{eqnarray} (3) where t, ρ and u denote the time, mass density and radial bulk flow velocity, respectively, and M = M(r, t) is the enclosed mass within a sphere of radius r at time t. The radial momentum equation with the Paczynski–Wiita gravity now reads \begin{eqnarray} \rho \left(\frac{\mathrm{\partial} }{\mathrm{\partial} t} +u\frac{\mathrm{\partial} }{\mathrm{\partial} r}\right)u = \!-\!\frac{\mathrm{\partial} p}{\mathrm{\partial} r} \!-\! \frac{GM\rho }{\left(r-2GM/c^2\right)^2} \!-\! \frac{\mathrm{\partial} }{\mathrm{\partial} r}\frac{\langle B^2_{\rm t}\rangle }{8\pi } \!-\! \frac{\langle B^2_{\rm t} \rangle }{4\pi r}, \nonumber\\ \end{eqnarray} (4) where p denotes the (effective) pressure and |$\langle B^2_{\rm t}\rangle$| is the ensemble mean-square magnitude of random transverse magnetic field |$\boldsymbol {B}_{\rm t}$|⁠. The third term on the right-hand side (RHS) is the negative gradient of a ‘magnetic pressure’ and the fourth term is the mean magnetic ‘tension’ always pointing radially inwards. These two terms together represent the net radial magnetic Lorentz force associated with the random transverse magnetic fields. In the limiting situation where these two magnetic force terms nearly cancel, we would have an approximately force-free magnetic field configuration (e.g. Low & Lou 1990; Wang & Lou 2014) in an MHD flow. Meanwhile, the law of magnetic induction gives \begin{equation} \left(\frac{\mathrm{\partial} }{\mathrm{\partial} t} +u\frac{\mathrm{\partial} }{\mathrm{\partial} r}\right) \left(r^2\langle B^2_{\rm t}\rangle \right) +2r^2 \langle B^2_{\rm t}\rangle \frac{\mathrm{\partial} u}{\mathrm{\partial} r}= 0. \end{equation} (5) Such a magnetofluid is assumed to obey a GP EoS for a formal specific entropy conservation along streamlines \begin{equation} \left(\frac{\mathrm{\partial} }{\mathrm{\partial} t} +u\frac{\mathrm{\partial} }{\mathrm{\partial} r}\right) \ln \left(\frac{p}{\rho ^\gamma }\right)=0. \end{equation} (6) Generally, we would take γ ≥ 1 but not equal to 4/3 in this analysis. For the γ = 4/3 case with Newtonian gravity employed, the interested reader may consult Goldreich & Weber (1980) and Lou & Cao (2008). In several astrophysical contexts, the possible range of 0 < γ < 1 is also explored, e.g. for a turbulent medium. Spaans & Silk (2000) studied a few models with a polytropic EoS of self-gravitating, quiescent interstellar gas clouds. In the regimes of either low or high densities, their models show 0 < γ < 1 (see their figs 1–3). Polytropes with γ < 1 are also explored to account for the observed linewidth–size relation in turbulent MCs (e.g. Maloney 1988). To model turbulent pressure effects in MCs, Lizano & Shu (1989), Gehman, Adams & Watkins (1996) and McLaughlin & Pudritz (1996) investigated polytropes in the γ → 0 regime for the so-called logatropes. To model possible Alfvén wave pressure effects in MCs, McKee & Zweibel (1995) took γ = 1/2 (see also McKee & Holliman 1999). Such γ < 1 models are also referred to as ‘negative-index polytropes’ with γ = 1 + 1/n and an index n < 0 (e.g. Shu et al. 1972; Viala & Horedt 1974a,b). 2.2 Self-similar transformation for quasi-spherical MHD Here, we introduce MHD self-similar transformation so that quasi-spherical non-linear MHD PDEs (1)–(6) can be consistently reduced to two coupled non-linear MHD ODEs in a self-similar form. This MHD self-similar transformation is explicitly given below as, \begin{eqnarray} r &=& k^{1/2}tx\ , \qquad u = k^{1/2}v(x)\ , \qquad M = \frac{k^{3/2}tm(x)}{G}, \nonumber\\ \rho &=& \frac{\alpha (x)}{4\pi G t^2}\ , \qquad p = \frac{k\beta (x)}{4\pi G t^2}\ , \qquad \langle B^2_{\rm t}\rangle = \frac{k w(x)}{G t^2}, \end{eqnarray} (7) where k is a constant coefficient with a dimension of velocity squared, x is the dimensionless self-similar independent variable, while the reduced dependent variables v(x), m(x), α(x), β(x) and w(x) are five dimensionless functions of x only; they stand for the reduced radial bulk flow speed, enclosed mass, mass density, pressure and magnetic energy density (or magnetic pressure), respectively. With MHD self-similar transformation (7), mass conservation PDEs (2) and (3) readily become \begin{equation} (x-v)\left(2\alpha /x + \alpha ^{\prime }\right) =\alpha v^{\prime } \end{equation} (8) and \begin{equation} m = (x-v)m^{\prime } = (x-v)x^2\alpha , \end{equation} (9) where the prime ΄ denotes the differentiation with respect to x and magnetic field induction PDE (5) can be cast into the ODE form of \begin{equation} (x-v)w^{\prime } - 2wv^{\prime } = -2(x-v)w/x. \end{equation} (10) A proper combination of non-linear ODEs (10) and (8) yields the following integral with an integration constant h referred to as the magnetic parameter \begin{equation} w(x) = h x^2\alpha ^2. \end{equation} (11) This integral relation is a straightforward consequence of the well-known ‘frozen-in’ condition of the magnetic flux for ideal MHD processes. Furthermore, GP EoS (6) leads to \begin{equation} \frac{2(\gamma -1)}{(x-v)}=\frac{\beta ^{\prime }}{\beta } - \gamma \frac{\alpha ^{\prime }}{\alpha } \end{equation} (12) by self-similar transformation (7) and can be combined with ODE (9) to yield a valuable algebraic GP EoS relation for reduced pressure β(x), viz. \begin{equation} \beta (x) = C\alpha ^\gamma m^{2\gamma -2} =C\alpha ^{3\gamma -2} (x-v)^{2(\gamma -1)} x^{4(\gamma -1)}, \end{equation} (13) where C is another integration constant. As noted in Lian & Lou (2014), for γ ≠ 4/3, any C can be normalized to unity without loss of generality when an appropriate value of k is chosen. Specifically, if one has a k = kold with C = Cold ≠ 1, then adopting the new choice |$k_{\rm new} = C_{\rm old}^{1/(4-3\gamma )}k_{\rm old}$| for k will do the job. However, this obviously does not work for γ = 4/3 and thus C would remain arbitrary. With two algebraic integral relations (11) and (13), MHD radial momentum ODE (4) takes the form of \begin{eqnarray} &&{ \left[ \alpha (x-v)+\frac{2(\gamma -1)\beta }{(x-v)}\right] v^{\prime } -\left[(3\gamma -2) \frac{\beta }{\alpha } + hx^2\alpha \right]\alpha ^{\prime }}\nonumber\\ && =\frac{2(\gamma -1)\beta }{(x-v)}+4(\gamma -1)\frac{\beta }{x} + \frac{m\alpha }{(x-sm)^2} + 2hx\alpha ^2, \end{eqnarray} (14) where s = 2k/c2 is a dimensionless parameter twice the square of the ratio of sound speed (or effective velocity dispersion) over the speed of light c in vacuum. ODE (14) can be readily combined with ODE (8) to yield two first derivatives α΄ and v΄ separately in the following forms of \begin{eqnarray} X(x,\ \alpha ,\ v)\alpha ^{\prime } =& Y(x,\ \alpha ,\ v), \end{eqnarray} (15) \begin{eqnarray} X(x,\ \alpha ,\ v)v^{\prime } =& Z(x,\ \alpha ,\ v), \end{eqnarray} (16) where the explicit expressions of the three coefficient functionals X, Y and Z are defined below in the order \begin{eqnarray} X(x,\alpha ,v) &\equiv & (x-v)^2 -\gamma \frac{\beta }{\alpha } - hx^2\alpha, \end{eqnarray} (17) \begin{eqnarray} Y(x,\alpha ,v) &\equiv & \frac{2\gamma {-}2}{x{-}v}\beta {-}\frac{2\alpha }{x}(x{-}v)^2 {+} \frac{m\alpha }{(x{-}sm)^2} {+}2hx\alpha ^2, \end{eqnarray} (18) \begin{eqnarray} Z(x,\alpha ,v) &\equiv & \left[ (2\gamma -2) - 2\gamma \frac{x-v}{x} \right] \frac{\beta }{\alpha } + \frac{m(x-v)}{(x-sm)^2}. \end{eqnarray} (19) Definition (19) of Z(x, α, v) does not involve h parameter. Coupled non-linear ODEs (15) and (16) with definitions (17) to (19) are the governing equations for our self-similar MHD model and are mutually consistent with ODEs (15) and (16) of Wang & Lou (2008) (note n = 1 there) by letting s → 0+ and with ODEs (13) and (14) of Lian & Lou (2014) by setting h = 0 there, respectively. ODEs (15) and (16) are also invariant under time-reversal operation (x, v, m, α, β, w)↦( − x, −v, −m, α, β, w). Therefore, we need only consider solutions with x > 0 (i.e. r > 0 and t > 0) and invoke a time reversal operation to obtain solutions with x < 0 (i.e. r > 0 and t < 0). The latter would have quite different interpretations or applications from their x > 0 counterparts (Lou & Shi 2014a). 3 ANALYTIC AND ASYMPTOTIC SOLUTIONS 3.1 Global magnetostatic solution singular as x → 0+ It is easily shown that the following velocity and density profiles \begin{equation} v(x)=0\ ,\qquad \quad \alpha (x)=A/x^2 \quad \qquad (\!\hbox{ for }\ A>0) \end{equation} (20) represent a magnetostatic solution of non-linear ODEs (15) and (16) if and only if parameter A is a positive root of algebraic equation \begin{equation} 2A^{3\gamma -4}(1-sA)^2=1 \end{equation} (21) within the open interval (0 , 1/s). The requirement of A < 1/s is due to the constraint that x − sm should remain positive in all pertinent ODEs for Paczynski–Wiita gravity (1) to be valid. With Paczynski–Wiita gravity (1), solution (20) is novel and represents a singular GP magnetostatic sphere with a force-free magnetic field configuration, as condition (21) is independent of the magnetic parameter h (e.g. Low & Lou 1990; Wang & Lou 2014). For such a force-free magnetic field, the average strength of the random transverse magnetic field, according to integral (11), varies as 1/r and remains independent of time t. While properties of equation (21) have been analysed and discussed by Lian & Lou (2014) in their contexts of two types of EWCSs without magnetic field; the pertinent results are nevertheless displayed here in Fig. 1 for a convenient reference. The possible domain of s and γ is subdivided into three regions. In the region of γ < 4/3, equation (21) has only one root for 0 < A < 1/s. For γ > 4/3, the number of roots for 0 < A < 1/s can be either 0 or 2, depending on the value of s, with a curve serving as the boundary between the two-root region and the no-root region in γ versus s plane. As Lian & Lou (2014) point out, the analysis of equation (21) turns out to be crucial for understanding patterns of the MSCC, a point we would explain more in subsequent sections presently. Figure 1. Open in new tabDownload slide The number of roots A to equation (21) can be 0, 1 or 2 depending on the two parameters γ and s. The condition for one single real root is γ < 4/3, in the ‘type A’ regime. For γ > 4/3, the number of roots can be 2 (type B regime) or 0 (type C regime) depending on s and the curve defined by equation (8/s3γ − 4)(3γ − 4)3γ − 4/(3γ − 2)3γ − 2 = 1 serves as the ‘borderline’, where equation (21) happens to have a double root (see appendix C of Lian & Lou 2014 for details). For all purposes, only the real roots that lie in the range 0 < A < 1/s are counted. In the critical case of γ = 4/3, equation (21) can be explicitly solved to yield |$A=(1\pm 1/\sqrt{2})/s$| and only the minus sign is acceptable whereas the plus sign root is larger than 1/s. Figure 1. Open in new tabDownload slide The number of roots A to equation (21) can be 0, 1 or 2 depending on the two parameters γ and s. The condition for one single real root is γ < 4/3, in the ‘type A’ regime. For γ > 4/3, the number of roots can be 2 (type B regime) or 0 (type C regime) depending on s and the curve defined by equation (8/s3γ − 4)(3γ − 4)3γ − 4/(3γ − 2)3γ − 2 = 1 serves as the ‘borderline’, where equation (21) happens to have a double root (see appendix C of Lian & Lou 2014 for details). For all purposes, only the real roots that lie in the range 0 < A < 1/s are counted. In the critical case of γ = 4/3, equation (21) can be explicitly solved to yield |$A=(1\pm 1/\sqrt{2})/s$| and only the minus sign is acceptable whereas the plus sign root is larger than 1/s. 3.2 Asymptotic MHD solution behaviours at large x 3.2.1 Asymptotic MHD solution approaching a constant velocity Physically, it would be reasonable to assume that the radial flow velocity v approaches a certain finite value V0 at large radii. By this assumption for large x, ODE (8) can be approximated to the zeroth order as \begin{equation} x\alpha ^{\prime }+2\alpha =0\ \qquad ({\rm for }\,\, x \gg 1), \end{equation} (22) suggesting that the leading-order expansion for α(x) bears a form of A0/x2 with an arbitrary constant A0 > 0. That is, a finite asymptotic velocity at large distance sets the profile scaling for the reduced mass density. It is then natural to search for asymptotic expansions in a power series of 1/x, viz. \begin{eqnarray} v(x) = V_0 + \frac{V_1}{x} + \frac{V_2}{x^2} + \frac{V_3}{x^3} + {\cal O}\left( \frac{1}{x^4} \right), \end{eqnarray} (23) \begin{eqnarray} \alpha (x) = \frac{1}{x^2}\left[ A_0 + \frac{A_1}{x} + \frac{A_2}{x^2} + \frac{A_3}{x^3}+\frac{A_4}{x^4} +{\cal O}\left(\frac{1}{x^5}\right)\right]. \end{eqnarray} (24) Substituting two power series (23) and (24) into coupled non-linear MHD ODEs (8) and (16), two sequences of coefficients Ai and Vi (with subscript i = 1, 2, 3, …) are successively determined as \begin{eqnarray} V_1 = 2A_0^{3\gamma -3}-\frac{A_0}{(1-sA_0)^2}, \end{eqnarray} (25) \begin{eqnarray} V_2 = \frac{sA_0^2 V_0}{(1-sA_0)^3}-(3\gamma -4)A_0^{3\gamma -3} V_0, \end{eqnarray} (26) \begin{eqnarray} V_3 &=& \frac{2A_0^{3\gamma -3}}{3}\left[ (2\gamma -5)A_0^{3\gamma -3} - (\gamma -2)(4\gamma -5)V_0^2 -A_0h \right]\nonumber\\ && +\, \frac{A_0^2h-4\gamma A_0^{3\gamma -2}}{3(1-sA_0)^2} + \frac{2A_0^{3\gamma -2}(3-4sA_0)}{3(1-sA_0)^3}\nonumber\\ && +\,\frac{s^2A_0^3V_0^2}{(1-sA_0)^4}-\frac{A_0^2(1-3sA_0)}{6(1-sA_0)^5}, \end{eqnarray} (27) and A1 = 0, \begin{eqnarray} A_2 = A_0^{3\gamma -2}-\frac{A_0^2}{2(1-sA_0)^2}, \end{eqnarray} (28) \begin{eqnarray} A_3 = \frac{2A_0}{3}\left[ (10-6\gamma )A_0^{3\gamma -3}V_0 - \frac{(1-3sA_0)A_0V_0}{(1-sA_0)^3}\right], \end{eqnarray} (29) \begin{eqnarray} A_4 &=& -\frac{A_0^{3\gamma -2}}{2}\left[ (2-2\gamma )A_0^{3\gamma -3}+ (4\gamma ^2-10\gamma +5)V_0^2+A_0h \right]\nonumber\\ && +\, \frac{A_0^3 h}{4(1-sA_0)^2} + \frac{2A_0^{3\gamma -1}[-\gamma +(\gamma -1)sA_0]-A_0^2V_0^2}{4(1-sA_0)^3}\nonumber\\ && +\, \frac{3sA_0^3V_0^2}{4(1-sA_0)^4} + \frac{A_0^3}{4(1-sA_0)^5}. \end{eqnarray} (30) The average (reduced) magnetic field strength can be directly calculated from integral relation (11) as w1/2(x) = h1/2xα(x). In terms of dimensional quantities, this mean field strength has a time-independent main component with a 1/r profile and other remaining time-varying components that decrease no slower than 1/r3 with increasing r. Note that the magnetic parameter h does not enter the leading orders of the asymptotic series expansion at large x. It is only in the third-order term of v(x) and the fourth-order term of α(x) that h appears for the first time. Therefore, dynamic behaviour at very large r (or small t) is not altered very much by the presence of random transverse magnetic fields and the asymptotic expansion solution essentially involves an asymptotic force-free magnetic field configuration (Low & Lou 1990; Wang & Lou 2014). 3.2.2 Absence of MHD thermal expansion solutions as x → +∞ In spite of the limiting behaviour v → V0 at large x just discussed, it may be tempting to assume a ‘Hubble-law’-like behaviour of v(x), i.e. |$v = \tilde{V}x + \tilde{U} + O(1/x)$| as x → +∞. However, this does not lead to physically allowable MHD thermal expansion solutions due to Pazcynski–Wiita gravity (1) and the requirement of sufficiently fast radial falloff of mass density as x → +∞. The specific mathematical analysis parallels that of section 3.2 in Lian & Lou (2014) and is thus not repeated here. For the Newtonian gravity in comparison, such MHD thermal expansion solutions do exist in the regime of large x (see Wang & Lou 2007, 2008 for details). 3.3 MHD solution behaviours approaching the expanding event horizon of an accreting Schwarzschild BH The Paczynski–Wiita gravity brings about one brand new feature to hydrodynamic (Lian & Lou 2014) and MHD formalisms that is missing by the Newtonian gravity model, namely the natural capability to describe the MHD collapse and mass accretion of central Schwarzschild BHs. The defining characteristic of a central Schwarzschild BH is an event horizon at finite and non-zero radius 2GM/c2. That is, there exists a certain x0 > 0 such that x − sm → 0 [or equivalently, (x − v)xα → 1/s] as x → x0 +. If this is indeed the case, it is reasonable to expect v(x) to approach −∞ as x → x0 +. ODE (8) then reduces to the following asymptotic form of \begin{equation} (\alpha v)^{\prime } + 2\alpha v/x = 0. \end{equation} (31) The resulting asymptotic behaviour is −αx2v → m0, where m0 is a constant of integration bearing the physical meaning of (dimensionless) mass accretion rate towards the central BH. Correspondingly, the expansion rate of this event horizon is x0 = sm0. Transforming back to the dimensional physical quantities using MHD self-similar transformation (7), the radius of the event horizon rEH is found to expand linearly with increasing time t in the following manner \begin{equation} r_\mathrm{EH} = \frac{s^{3/2}m_0c}{\sqrt{2}}t=2GM_{\rm s}/c^2, \end{equation} (32) where Ms is the central Schwarzschild BH mass. For our model analysis, it is instructive to write MHD radial momentum ODE (14) in the following equivalent form of \begin{equation} \frac{v}{\alpha }\frac{{\rm d}\beta }{{\rm d}x} +v \frac{{\rm d}(hx^2\alpha )}{{\rm d}x} + (x-v)\frac{{\rm d}}{{\rm d}x} \left( \frac{m}{x-sm} - \frac{v^2}{2} \right) = 0. \end{equation} (33) The factor (x − v) in the last term on the left-hand side (LHS) of ODE (33) can be replaced by −v in the vicinity of the event horizon x0 and we also have dm/dx → 0 as x → x0 from ODE (9). Therefore, m can be replaced by m0 to the first approximation. Using the first equality in integral relation (13) (with C = 1 there) to eliminate β(x), ODE (33) can be integrated with respect to x, resulting in \begin{equation} \frac{v^2}{2} - \frac{m_0}{x-sm_0} + \frac{\gamma }{\gamma -1}\alpha ^{\gamma -1}m_0^{2\gamma -2} + hx^2\alpha = v_d^2/2, \end{equation} (34) where vd is an integration constant.1 Integral relation (34) is the quasi-spherical MHD counterpart of the ‘dynamic Bernoulli relation’ that was first pointed out by Lou (2015b) in the context of GP dynamic cylinders (see also equation 34 of Lou & Hu 2016 without magnetic field and the q parameter there should be the one for spherical geometry; Wang & Lou 2008). The four terms on the LHS are the gas kinetic energy, PW gravitational potential energy, enthalpy and magnetic field energy per unit mass, respectively. As x → x0 +, the first two terms diverge while the last two terms both tend to zero2 for γ ≠ 1, indicating essentially a free-fall gas approaching the event horizon. This reminds us of the free-fall solutions towards the centre in Newtonian gravity models (e.g. Shu 1977; Lou & Shen 2004; Wang & Lou 2008), which may be regarded as a Newtonian analogy to our asymptotically diverging MHD solutions with finite mass accretion rates under PW-gravity towards an expanding event horizon. Outside the event horizon, the internal and magnetic energies will gradually come into play. From MHD Bernoulli relation (34), we immediately find the asymptotic behaviours for reduced velocity v and density α: \begin{eqnarray} v(x) \cong -\left(\frac{2m_0}{x-sm_0}+v_d^2\right)^{1/2}+\cdots , \end{eqnarray} (35) \begin{eqnarray} \alpha (x) \cong \left[ s^2m_0\left(\frac{2m_0}{x-sm_0} +v_d^2\right)^{1/2}\right]^{-1} +\cdots \ \end{eqnarray} (36) to the leading order as x → x0 +. Meanwhile, by integral relation (11), the average magnetic field strength tends towards zero in the same trend as α near the expanding event horizon. The averaged radial component 〈Sr〉 of the Poynting vector |$\boldsymbol {S} = (c/4\pi )\boldsymbol {E}\times \boldsymbol {B}$| represents the electromagnetic energy flux density. For the present MHD model problem, we readily identify \begin{equation} \langle S_{\rm r} \rangle = \frac{u}{4\pi }\langle B_{\rm t}^2 \rangle , \end{equation} (37) which vanishes as x → x0, meaning no electromagnetic energy flux density across the moving event horizon towards the mass growing BH. On the other hand, the conversion rate of magnetic field energy to other forms of energy (i.e. kinetic, internal and gravitational energies) can be found by calculating the averaged power density (PD) of the radial magnetic force, given by \begin{eqnarray} \langle \mathrm{{\rm PD}} \rangle &=& -\frac{u}{4\pi } \left[ \frac{1}{2}\frac{\mathrm{\partial} \langle B_{\rm t}^2 \rangle }{\mathrm{\partial} r} + \frac{\langle B_{\rm t}^2 \rangle }{r} \right] \nonumber\\ &&= {} -\frac{k}{4\pi G t^3}\left\lbrace \frac{v}{2}\frac{{\rm d}w}{{\rm d}x} + \frac{vw}{x} + \frac{s}{2}\left[ v^2 w + \frac{x}{2}\frac{{\rm d}(v^2w)}{{\rm d}x} \right]\right\rbrace ,\nonumber\\ \end{eqnarray} (38) where MHD self-similar transformation (7) has been applied in the second line. Substituting asymptotic free-fall solution (35) and (36) into (38), we arrive at the asymptotic expression below \begin{equation} \langle \mathrm{{\rm PD}} \rangle = \frac{kh}{2^{7/2} \pi G s^2 m_0^{1/2} t^3}\frac{1}{(x-x_0)^{1/2}} + \mathcal {O}(1), \end{equation} (39) for x → x0 +. While |$\langle \mathrm{{\rm PD}} \rangle$| diverges at x0, its spatial integral does not, and the total power remains finite. By expression (39), the electromagnetic energy, instead of being directly accreted towards the central BH, is first converted to kinetic and internal energy of the fluid outside the expanding event horizon. While the Paczynski–Wiita gravity is capable of reproducing the dynamical effects of a curved spacetime around a Schwarzschild BH to a very good accuracy, it does not take into account of special relativistic effects to avoid infinite flow speed though. As a result, the radial flow velocity u = k1/2v(x) tends to −∞ as x → x0 +, rather than being limited by the speed of light c in vacuum. This issue was highlighted by Abramowicz et al. (1996) in the context of accretion disc modelling, and a simple empirical rescaling was suggested, i.e. one should relate the true physical flow speed utrue and the PW-gravity model calculated speed u by the rescaling relation3 \begin{equation} u = u_\mathrm{true}/(1-u_\mathrm{true}^2/c^2)^{1/2}. \end{equation} (40) In our model, the true radial gas inflow speed near the event horizon now approaches −c by this empirical readjustment. Further discussion on this aspect is presented in Section 7. 3.4 MHD solution behaviours near an expanding central void From ODE relation (9), it follows that a non-negative enclosed mass m(x) requires x to be always no smaller than v(x). However, it is physically possible that a solution may have v = x at a particular point x = x1 > 0, implying that no mass resides inside the dimensionless radius x1 and the corresponding physical boundary expands with increasing time t. Under this situation, the self-similar MHD solution should not continue further inwards and we have thus created a central void with its rim dynamically expanding in a self-similar manner. In plausible astrophysical situations like expanding PNe, super bubbles and supernova remnants etc., some sort of central ‘power engine’ such as a hot fast tenuous stellar wind, an energetic photon or neutrino gas sphere, a relativistic electron–positron pair plasma sphere or dark energy (in the case of large-scale cosmic voids), would be needed to drive and maintain the dynamic evolution of such an expanding void. A few years ago, Lou & Zhai (2009, 2010) invoked self-similar void solutions with the void filled by a hot tenuous fast stellar wind and with piecewise isothermal shocks to study the diffuse X-ray emission from PNe (e.g. PN NGC 40 and PN NGC 7662). In the case of PN NGC 7662, its X-ray luminosity was model predicted to vary on time-scales of several decades due to a reverse shock propagating inwards in a hot tenuous fast stellar wind of ∼1500 km s−1 (Lou & Zhai 2009, 2010). For observations over decades, the recent report on SAO 244567, the exciting star of the Stingray Nebula, reveals its rapid evolution from an effective temperature of about 21 kK in 1971 to over 50 kK in the 1990s. Updated data and model analyses of Reindl et al. (2014) further indicate that the effective temperature of SAO 244567 increases from 38 kK in 1988 to a peak value of 60 kK in 2002 involving a terminal wind velocity increase from v∞ = 1800 km s−1 to 2800 km s−1 and since then, the star stopped heating and has cooled down again to 55 kK by 2006. By the PN scenario of Lou & Zhai (2009, 2010), the (proto)PN system SAO 244567 may well involve a similar inward propagating reverse shock in a hot tenuous fast stellar wind that causes its temperature variations over a time scale of several decades. Lou & Cao (2008) discussed the self-similar dynamics of a relativistically degenerate or hot (γ = 4/3) gas sphere, finding that it allows for the emergence of a central void with sharp edge surrounded by an overdense gas shell. The case of γ → 4/3 was investigated and then applied to the dynamic evolution of supernova ejecta in Lou & Wang (2012). Our analysis here closely relates to that of Lian & Lou (2014) (specifically their appendix A), yet with several new solution features emerging in the presence of magnetic field. 3.4.1 Asymptotic MHD void solution with vanishing density To examine the MHD solution behaviour near an expanding central void, we first expand v(x) in terms of δ, the outward displacement from the position of void rim x1 [i.e. x = x1 + δ with v(x1) = x1]: \begin{equation} v = x_1 + K \delta + \mathcal {O}(\delta ^2). \end{equation} (41) We must demand K < 1 in expansion (41) with a positive m(x) for small δ > 0. Assuming K ≠ 0 also for the moment, the asymptotic form of α(x) is easily determined by equation (8) as \begin{equation} \alpha = \alpha _K\delta ^{K/(1-K)} + \cdots , \end{equation} (42) where αK is an integration constant. Clearly, the qualitative behaviour of α near x1 is determined by the sign of K. As a divergent density is not expected in general in the vicinity of a void, we only seek solutions with K > 0. This is achieved by further expanding ODE (16) into a series of δ and balancing the lowest order terms of both sides. Interestingly, the result depends on γ value and is summarized as follows: for 1 < γ < 5/4 (γ = 1 does not fall in this category), there are two possible solutions. The first solution has the radial inertial force balanced by the gas pressure force and we obtain \begin{equation} K = 2(2-\gamma )/(\gamma +1), \end{equation} (43) while the corresponding coefficient αK should satisfy the relation below \begin{equation} K(1-K)^{4-2\gamma } = (\gamma K + 2\gamma -2) \alpha _K^{3\gamma -3}x_1^{4\gamma -4}\ . \end{equation} (44) The second solution is dominated by the magnetic force and the gas pressure. It is determined to be \begin{equation} K = 2(1-\gamma )/(\gamma -2), \end{equation} (45) and the corresponding relation for αK involving magnetic parameter h reads \begin{equation} hK(1-K)^{2-2\gamma } + (\gamma K + 2\gamma - 2)\alpha _K^{3\gamma -4}x_1^{4\gamma -6} = 0\ . \end{equation} (46) In the regime of γ ≥ 5/4, however, the two aforementioned MHD solutions merge into one. For γ = 5/4, the radial inertial force, the gas pressure force and the magnetic force are all of the same order. Yet, for γ > 5/4, the gas pressure drops out as a negligible term. As a result, we have \begin{eqnarray} K = 2/3 \,\,\, {\rm and}\,\,\, hx_1^2\alpha _K = \left\lbrace \begin{array}{@{}l@{\quad }l@{}}\displaystyle 1/9 - 2\alpha _K^{3/4} x_1/\sqrt{3}\ , \ & \gamma = 5/4,\\ \displaystyle 1/9\ , \ & \gamma > 5/4. \end{array}\right.\nonumber\\ \end{eqnarray} (47) The required inequality K > 0 here carries a direct physical interpretation: the outflow velocity of the fluid immediately outside the void rim is larger than the expansion speed of the void itself. Consequently, the area surrounding the void tends to be evacuated and the density α drops to zero there. This is to be contrasted with the void solution obtained with γ = 4/3 by Lou & Cao (2008), where the density diverges at the void rim. A similar divergent behaviour can also be produced in our MHD model if we allow negative K values, which we would not pursue further in this study. 3.4.2 Asymptotic MHD void solution with finite density The foregoing discussion has set aside the important case where the derivative v΄(x1) vanishes, which we must now address by changing the postulated form of velocity profile into \begin{equation} v = x_1 + L \delta ^{P} + \cdots \end{equation} (48) with the power index P > 1. From equation (8), the asymptotic MHD behaviour of the reduced mass density is determined as \begin{eqnarray} \alpha = \alpha _L \times \left\lbrace \begin{array}{@{}l@{\quad }l@{}}\displaystyle { 1 + PL\delta ^{P-1}/(P-1)}+\cdots \ , & \quad P < 2 ; \\ \displaystyle { 1+ 2\left(L - 1/x_1\right)\delta }+\cdots \ , & \quad P=2; \\ \displaystyle { 1 - 2\delta /x_1}+\cdots \ , &\quad P > 2, \end{array}\right.\nonumber\\ \end{eqnarray} (49) with an arbitrary positive integration constant αL. The balance of lowest order terms in ODE (16) again yields γ-dependent results for the two solution parameters L and P. Specifically, for 1 ≤ γ < 2, the dominant forces along the void rim are the magnetic force and the gas pressure force and we arrive at \begin{equation} P = 2 \qquad {\rm and }\qquad L = 1/\left[x_1(hx_1^2\alpha _L+1)\right] \end{equation} (50) for γ = 1, and \begin{equation} P = 2\gamma - 1\ \qquad {\rm and }\qquad L = -\frac{2(\gamma -1)}{h(2\gamma -1)} {\alpha _L^{3\gamma -4}x_1^{4\gamma -6}} \end{equation} (51) for 1 < γ < 2. While for γ ≥ 2, the magnetic force and the gravity are dominant along the void rim (except for the γ = 2 case where the magnetic force, the gravity and the gas pressure are all of the same order) and we derive \begin{equation} P = 3\ {\rm and }\ L = \left\lbrace \begin{array}{@{}l@{\quad }l@{}}\displaystyle -2\alpha _L^2x_1^4/(3h) - 1/(3hx_1^2), &\quad \gamma =2; \\ \displaystyle -1/(3hx_1^2), &\quad \gamma > 2. \end{array}\right. \end{equation} (52) The coefficient L is always negative for γ > 1 by MHD solutions (51) and (52) and its magnitude is inversely proportional to magnetic parameter h. This means that the expansion velocity of the void is always larger than the outflow speed of gas immediately surrounding it, although the difference is gentle (higher than the first order in δ) and can be further smoothed out by a larger h. As a result, the reduced density at the void rim does not approach zero or infinity, but remains at finite positive values. 3.4.3 Role of random transverse magnetic fields in MHD solutions with central void expansion When magnetic parameter h > 0, the mean reduced magnetic field strength w(x)1/2 follows h1/2xα(x) in the vicinity of a void by magnetic flux ‘frozen-in’ condition (11). When h = 0, expressions (46) and (47) for αK as well as expressions (51) and (52) for L all break down, highlighting that these MHD solutions are novel features unique to a magnetofluid, whereas an unmagnetized fluid as studied by Lian & Lou (2014) does not admit these sorts of behaviours. In conceivable astrophysical flow systems, such a sharp transition from unmagnetized to slightly magnetized situations would appear qualitatively similar to the case of magneto-rotational instabilities (MRIs) in the context of generating turbulence in accretion discs. We leave these considerations open for future works and regard this novel behaviour as a theoretical possibility for the moment. Numerical MHD solution examples corresponding to both kinds of asymptotic velocity profiles (41) and (48) will be constructed in Section 6. 3.5 The quasi-magnetostatic asymptotic solution as x → 0 In certain physical situations prior to the self-similar MHD stage, it is entirely possible that neither BHs nor voids emerge around the centre and the self-similar dynamic process is described by an MHD solution that extends to the very centre x = 0. In this section, we investigate quasi-magnetostatic asymptotic solution towards the very centre x → 0. We assume that v(x) ∼ UxQ where the two constant coefficients U and Q can be complex to also describe possible oscillatory behaviours. By ‘quasi-magnetostatic’, we mean Re(Q) > 1 and the leading order asymptotic approximation of ODE (8) reads \begin{equation} \frac{1}{x^2\alpha } \frac{{\rm d}(x^2\alpha )}{{\rm d}x} = UQx^{Q-2}\ . \end{equation} (53) This gives a series expansion of α(x) in the form of \begin{equation} \alpha =\frac{R}{x^2}\left[ 1+\frac{QU}{(Q-1)}x^{Q-1} + \cdots \right], \end{equation} (54) which can be regarded as a global magnetostatic solution (20) perturbed at small x (hence the term ‘quasi-magnetostatic’). For Paczynski–Wiita gravity (1) to be valid, we still require 0 < R < 1/s. The magnetic field has a time-independent dominant term proportional to 1/r as in the global magnetostatic solution and further contains time-decaying terms in the form of (r/t)Q − 1/r + …. Substituting the asymptotic forms of v(x) and α(x) above into ODE (16) and keeping up to the first order of xQ − 1, we obtain \begin{eqnarray} &&{\left(\gamma R^{3\gamma \!-\!3}\!+\!hR\right)QUx^{Q-1} = -2R^{3\gamma -3}\left( 1\!+\!\frac{3\gamma -Q-2}{Q-1}Ux^{Q-1} \right)} \nonumber\\ && +\frac{R}{(1-sR)^2}\left[ 1\!+\!\frac{U}{Q-1} \left(\frac{2sR}{1\!-\!sR}\!-\!Q+2\right) x^{Q\!-\!1}\right]\ . \end{eqnarray} (55) Comparisons of the pertinent coefficients of x0 and xQ − 1 terms on both sides of equation (55) lead to two relations below \begin{eqnarray} 2R^{3\gamma -4}(1-sR)^2 = 1; \end{eqnarray} (56) \begin{eqnarray} &&{\left( \gamma +hR^{4-3\gamma } \right)(Q-1)^2 +\left(\gamma +hR^{4-3\gamma } \right)(Q-1)} \nonumber\\ && \quad +\, 2\left(\frac{1+sR}{1-sR}+3-3\gamma \right)= 0\ . \end{eqnarray} (57) Equation (57) is quadratic in (Q − 1). For a root with Re(Q − 1) > 0 to exist, the two parameters γ and R should satisfy the inequality \begin{equation} R < (3\gamma -4)/[(3\gamma -2)s], \end{equation} (58) and the root of Q is necessarily real. As equation (56) is identical to equation (21), the earlier discussion in Section 3.1 concerning the roots of R remains valid here. For instance, we readily conclude that for the two parameters (γ, s) within the type C regime of Fig. 1, it is impossible to exhibit the behaviour under consideration here. Moreover, the RHS of inequality (58) should be at least positive, requiring γ > 4/3 necessarily.4 This then excludes the type A regime of Fig. 1 for consideration. For (γ, s) in the type B regime of Fig. 1, the upper bound for R as given by inequality (58) coincides with the maximum of the LHS of equation (56) (viewed as a function of R). It follows that only the smaller R solution to equation (56) would satisfy that bound given by inequality (58). In other words, only one of the two global static solutions of the form (20) admits a quasi-magnetostatic asymptotic MHD solution. When s = 0, our results here are consistent with the quasi-magnetostatic solutions of Wang & Lou (2008) for Newtonian gravity, as expected and as an important necessary check. By equa-tions (56) and (57), a finite s > 0 does have an effect on the expansion coefficient R and index Q. Afterall, gravity is one of the dominant forces here. In particular, R increases and Q decreases with increasing s parameter. 3.6 A novel asymptotic MHD solution and the generalized LP-type MHD solution We now study the case of Q = 1, i.e. v(x) ∼ Wx, where W ≠ 0. A positive enclosed mass m(x) > 0 again requires W < 1 and the critical case of W = 1 is excluded. ODE (8) then reduces to \begin{equation} \frac{\alpha ^{\prime }}{\alpha } = \left( -2 + \frac{W}{1-W} \right)x^{-1}, \end{equation} (59) which can be integrated to derive the asymptotic solution \begin{equation} \alpha (x) = T x^{-2 + W/(1-W)} + \cdots , \end{equation} (60) where T is an integration constant (see the dashed MHD solution curve in Fig. 12 and the dashed a΄ MHD shock solution in Fig. 18 as an example of an infinite discrete MHD solution sequence). Imposing the condition x − sm > 0, we must have 0 < W < 1. Sub-stituting the asymptotic forms of v(x) and α(x) for small x > 0 into non-linear ODE (16) and carefully examining every possibility, we identify two different situations depending on the ranges of γ value. For 1 ≤ γ < 4/3, we find a brand new class of MHD solutions that can be further divided into two subclasses.5 The first subclass has W = 2/(3γ) and an arbitrary T > 0, while the second subclass has W = 2/(3γ − 1) and \begin{equation} T = \left[ \frac{(3\gamma -1)}{2}W(1-W)^{4-2\gamma } \right]^{1/(3\gamma -3)}\ . \end{equation} (61) For γ > 1, both subclasses have divergent α towards the centre, whereas for γ = 1 (i.e. isothermal EoS), the second subclass ceases to exist and the first subclass becomes the isothermal Larson–Penston (LP-) type MHD solution (see e.g. Larson 1969; Penston 1969; Lou & Shen 2004; Yu & Lou 2005; Yu et al. 2006), with α approaching a constant T value as x → 0. As will be clear in later sections, in general, the first subclass of solutions, when extended to larger x through numerical integration, will have radial oscillations and cross the MSCC in a manner similar to isothermal LP solutions. For these aspects, these solutions appear similar to LP-type MHD solutions. Nevertheless, accidentally this term has already been given to another type of qualitatively different solutions (Lou & Shi 2014b) that we discuss later. In this paper, we thus refer to the first subclass and second subclass of MHD solutions here as LX-1 and LX-2, respectively. Also, as s parameter does not enter the lowest order expansion of LX-1 and LX-2, we naturally expect these novel solutions to also exist even with Newtonian self-gravity, which are unknown until this research publication. For γ > 4/3, only W = 2/3 provides acceptable results and there is a further relation (1/6 + h)T = 1/9 that fixes T coefficient for a given magnetic parameter h. Here, α approaches a finite positive constant as x → 0, similar to the isothermal LP solution. It is for this reason that this solution was given the name ‘generalized LP solution’ in Lou & Shi (2014b), where it was first discovered for a GP hydrodynamic flow; we now extend this type of solution to GP MHD flows. In this paper, we follow this convention of terminology. The magnetic field of MHD solutions found in this section has a dominant term proportional to (r/t)W/(1 − W)/r. This profile is always shallower than that of the quasi-magnetostatic solution. 4 JUMP CONDITIONS FOR MHD SHOCK WAVES In this section, we discuss MHD shock wave conditions in the framework of self-similar quasi-spherical MHD. MHD shock waves are characterized by drastic changes (or jumps) in pertinent mechanical, electromagnetic and thermodynamic variables of MHD flows across a spherical surface of discontinuity. Such an MHD discontinuity needs to satisfy a set of conservation laws, referred to as the MHD shock conditions. For non-magnetized fluids, these conditions include mass continuity, energy conservation and momentum conservation in the comoving reference framework of a shock. For magnetofluids, there are additional electromagnetic shock conditions to be required. We use subscript ‘−’ for quantities immediately at the front of a shock (where gas flows into the shock, i.e. the upstream side) and by subscript ‘+’ for quantities at the back side of a shock (where gas flows away from the shock, i.e. the downstream side). We adopt a pair of square brackets for a subtraction of whatever physical quantity f on the downstream side from the same quantity on the upstream side, i.e. [f] ≡ f− − f+. The condition for mass continuity across an MHD shock is therefore \begin{equation} [\rho (u-u_{\rm s})] = 0, \end{equation} (62) where us is the radial travel speed of an MHD shock as marked by a subscript s. Meanwhile, the condition for the MHD energy conservation reads \begin{equation} \left[ \frac{\rho }{2}(u-u_{\rm s})^3 + \frac{\gamma p}{(\gamma -1)}(u-u_{\rm s}) + \frac{\langle B^2_{\rm t}\rangle }{4\pi }(u-u_{\rm s})\right] = 0, \end{equation} (63) where the first two terms are kinetic and thermal energy fluxes and the third term is the radial Poynting flux due to the radial transport of magnetic energy density. The radial momentum conservation gives \begin{equation} \left[ p + \frac{\langle B^2_{\rm t}\rangle }{8\pi } + \rho (u-u_{\rm s})^2 \right] = 0\ . \end{equation} (64) Finally, we require the average magnitude of the tangential electric field component to be continuous across such an MHD shock, viz., \begin{equation} \left[ (u-u_{\rm s})^2\langle B^2_{\rm t}\rangle \right] =0\ . \end{equation} (65) The above four MHD shock conditions (62)–(65) can be cast into a self-similar form using MHD transformation (7), with the sound parameter k changing abruptly from k− to k+ as the flow goes across an MHD shock. This is because flowing through an MHD shock wave front always changes the gas entropy (see equations 6 and 13), where k value is chosen to give C = 1 in equation (13). We assume no relationship between the magnetic parameters h− and h+ a priori. These MHD shock conditions in self-similar form are then given by \begin{eqnarray} \left[ k^{1/2}(x-v)\alpha \right] = 0, \end{eqnarray} (66) \begin{eqnarray} \left[ k^{3/2}\left( \frac{\alpha }{2}(x-v)^3 + \frac{\gamma \beta }{\gamma -1} (x-v)+w(x-v)\right)\right]=0, \end{eqnarray} (67) \begin{eqnarray} \left[ k\left(\beta +\frac{w}{2}+\alpha (x-v)^2\right)\right] = 0, \end{eqnarray} (68) \begin{eqnarray} \left[ k^2 w (x-v)^2 \right] = 0. \end{eqnarray} (69) In addition to the above set of four conservation equations for an MHD shock, we also require the radial location of the shock to be consistently represented by the self-similar solutions on both sides, which implies \begin{equation} \left[ k^{1/2}x \right]= 0. \end{equation} (70) The requirement of h− = h+ = h follows immediately from conditions (66), (69) and (70), meaning that the magnetic flux frozen-in condition remains valid across a self-similar MHD shock. General procedures for solving the self-similar MHD shock conditions are described below in order. This procedure is used in the subsequent numerical computations to construct global solutions involving outgoing MHD shocks. Suppose that the variables on one side of the shock, say x+, v+ and α+, are known and a solution for the corresponding variables on the other side is sought. The symbol pair distinguished by subscripts ‘+’ and '−' \begin{equation} \Gamma _{\pm }=(x_\pm - v_\pm )/x_\pm \end{equation} (71) is introduced here for the convenience of MHD shock analysis. We first divide equations (66) through (68) by suitable powers of equation (70) to eliminate k parameter and then use the resulting equations from (67) and (68) to replace the ratio |$\beta _-/x_-^2$| in favour of α− and Γ−. Using equation (66), occurrences of α− can be further replaced in favour of Γ−. Finally, we arrive at a quadratic equation of Γ− in terms of pertinent variables with subscript +: \begin{eqnarray} \frac{\gamma +1}{2\gamma }\Gamma _-^2 \!-\!\left(\frac{\gamma -1}{2\gamma }\Gamma _+ \!+\! \frac{\beta _+}{\Gamma _+\alpha _+x^2_+} \!+\! \frac{h\alpha _+}{2\Gamma _+} \right) \Gamma _- \!-\! \frac{2-\gamma }{2\gamma }h\alpha _+ = 0 .\nonumber\\ \end{eqnarray} (72) In these derivations, we have dropped a factor of Γ+ − Γ− ≠ 0 throughout. The solutions to quadratic equation (72) further yield α− and |$\beta _-/x_-^2$|⁠, which, by the definition of β−, gives the following expression for x−: \begin{eqnarray} x_-^{6\gamma -8} &=& \frac{\Gamma _-^{\gamma -1}}{(\Gamma _+\alpha _+)^{3\gamma -3}} \left[ \frac{(\gamma -1)}{2\gamma }\left( \Gamma _+^2 - \Gamma _-^2 \right) \right.\nonumber\\ && \left. +\frac{(\gamma -1)h\alpha _+}{\gamma } \left(1-\frac{\Gamma _+}{\Gamma _-}\right) + \frac{\beta _+}{\alpha _+x_+^2} \right]\ . \end{eqnarray} (73) Solving equation (73) then determines all the relevant variables at the upstream side (− subscript) of an MHD shock using the downstream variables (+ subscript). 5 MAGNETOSONIC CRITICAL CURVE (MSCC) 5.1 Geometric features of the MSCC Coupled non-linear MHD ODEs (15) and (16) encounter a singularity when X(x, α, v) = 0. Physically, this means that the radial flow speed relative to the local self-similar profile equals the magnetosonic wave speed of radially propagating disturbances. The magnetosonic wave speed is simply the square root of the sum of squares of the sound speed and the well-known Alfvén wave speed. The relation X(x, α, v) = 0 defines a surface in the three-dimensional space spanned by x, α and v, which we call the magnetosonic critical surface (MSCS). Physical solutions cannot smoothly cross this MSCS, unless Y(x, α, v) and Z(x, α, v) vanish simultaneously at the crossing point. The locus of points in the (x, α, v) space along which X, Y and Z all vanish is referred to as the MSCC. By substituting β expression (13) into the equations X = 0 and Y = 0, we have \begin{eqnarray} &(x-v)^2 - \gamma \alpha ^{3\gamma -3} (x-v)^{2\gamma -2} x^{4\gamma -4} - hx^2\alpha = 0, \end{eqnarray} (74) \begin{eqnarray} &&{ \left[ (2\gamma -2) - 2\gamma \frac{x-v}{x} \right] \alpha ^{3\gamma -3} (x-v)^{2\gamma -2} x^{4\gamma -4} } \nonumber\\ && + \frac{(x-v)^2 \alpha }{[1-sx(x-v)\alpha ]^2} = 0\ . \end{eqnarray} (75) The third equation Z = 0 is redundant (i.e. already implied by equations 74 and 75), and we have confirmed that the MSCC is indeed a curve rather than a surface or point in the (x, α, v) space. Besides, the straight line x = v, α = 0 is evidently a very special part of the MSCC, in the sense that physical solution can only reach but not cross it, as the reduced enclosed mass m vanishes along this special line. Naturally, we refer to this line as the zero-mass line (ZML). The LHS of equation (74), with γ > 1 and given x and v, is monotonically decreasing with increasing α. It takes a non-negative value when α = 0 and approaches −∞ as α increases indefinitely. Therefore, a unique non-negative α exists for equation (74) to hold. Symbolically, we denote this as |$\alpha = {\cal F}_{\gamma ,\ h}(x,\ v)$|⁠. Substitution of this α into equation (75) yields an equation in the form of |${\cal G}_{\gamma ,\ h,\ s}(x,\ v)=0$|⁠, where |${\cal G}$| is some function whose analytical expression may not be readily obtainable. However, suitable numerical methods can be employed to construct the solution set of |${\cal G}$|⁠, which is the projection of the MSCC on to the −v versus x plane. To fully understand the shape of the MSCC, some analytical geometric insight is valuable. First, it is instructive to consider the intersections of the MSCC and the half-plane v = 0 with x > 0. Elementary algebraic manipulations show that the x and α values at the intersection(s) (excluding x = 0) are given by \begin{equation} x = (\gamma A^{3\gamma -3} + hA)^{1/2}\ ,\qquad \alpha = A/x^2, \end{equation} (76) where 0 < A < 1/s is any root of algebraic equation (21). The behaviour of the roots of equation (21) has been fully analysed in Section 3.1 and will not be repeated here. The magnitude of magnetic parameter h only influences the position of the intersections but not their number, because A bears no relationship to h. The reason for the emergence of the same algebraic equation from the investigation of global magnetostatic solutions and MSCCs lies in the fact that those global solutions are special solutions lying in the plane v = 0 that cross the MSCC smoothly. Fig. 2 is a plot of MSCCs projected on to the −v versus x plane, with s = 0.2929 and h = 4.0 chosen as fixed and γ varying from 1.05 to 1.8. Fig. 2 here is essentially the ‘magnetized’ version of fig. 2 in Lian & Lou (2014), showing the same morphological patterns. This pattern can be intuitively interpreted with the help of Fig. 1 as follows: as γ increases, the point (γ, s) in the coordinates of Fig. 1 moves along a horizontal line. This line is tangent to the demarcating boundary between regions B and C at its lowest point, by our deliberate choice of the special value s = 0.2929. As shown clearly in Fig. 2, the overall shape of the MSCC is an elongated loop. For γ < 4/3, the point (γ, s) is in the type A regime and the corresponding MSCC loop touches the ZML (not shown) at the point x = v = 0 and makes only one intersection with the v = 0 plane at a positive x value. For γ > 4/3, the point (γ, s) is in the type B regime and the MSCC loop no longer touches the ZML but makes two intersections with the v = 0 plane at positive x values. Interestingly, when γ is close to 1.6 so that the parameter pair (γ, s) falls very close to the type C regime in Fig. 1, the MSCC loop gradually shrinks to very small sizes. One can naturally conceive this limiting process. Figure 2. Open in new tabDownload slide The projection of MSCCs on to the −v versus x plane for different values of γ with fixed s = 0.2929 and h = 4.0, excluding the ZML which is not shown. Those MSCC projections form loops, intersecting the positive x axis twice when γ > 4/3 but only once when γ < 4/3 (excluding x = 0). For γ < 4/3, the MSCC loop also touches the ZML at x = 0. Note that when γ is close to 1.6 so that the parameter pair (γ, s) falls close to the type C regime in Fig. 1, the MSCC loops shrink to very small sizes. Figure 2. Open in new tabDownload slide The projection of MSCCs on to the −v versus x plane for different values of γ with fixed s = 0.2929 and h = 4.0, excluding the ZML which is not shown. Those MSCC projections form loops, intersecting the positive x axis twice when γ > 4/3 but only once when γ < 4/3 (excluding x = 0). For γ < 4/3, the MSCC loop also touches the ZML at x = 0. Note that when γ is close to 1.6 so that the parameter pair (γ, s) falls close to the type C regime in Fig. 1, the MSCC loops shrink to very small sizes. This loop shrinking phenomenon as observed in Fig. 2 strongly suggests that there might be a complete disappearance of the MSCC loop for parameter pair (γ, s) falling inside the type C regime of Fig. 1. This is indeed the case. To better appreciate this fact, we shall further consider, more generally, the intersection of the MSCC loop and the half-plane v = λx for x > 0, with λ being a parameter. The idea is that as λ decreases from 1 to −∞, this half-plane sweeps over the whole physically attainable region in the (x, α, v) spanned space, performing a ‘scan’ of the MSCC. Using the same method above, it can be shown that the x-coordinates of the intersections satisfy \begin{equation} x=[\gamma (1-\lambda )^{2\gamma -2} A^{\prime 3\gamma -3}+hA^{\prime }]^{1/2}/(1-\lambda ), \end{equation} (77) where A΄ is any root of the following algebraic equation \begin{equation} (2-2\gamma \lambda )(1-\lambda )^{2\gamma -4} A^{\prime 3\gamma -4}[1-s(1-\lambda )A^{\prime }]^2 = 1 \end{equation} (78) within the range 0 < A΄ < 1/[(1 − λ)s]. For λ = 0, this equation reduces back to equation (21), as required by consistency. Analysis can be readily carried out concerning the roots A΄ and the results are similar to those obtained for the λ = 0 case. However, there are two important modifications. First, the boundary between the two-roots region B and the no-root region C becomes \begin{equation} s^{3\gamma -4} = \frac{8(3\gamma -4)^{3\gamma -4}}{(3\gamma -2)^{3\gamma -2}} \frac{(1-\gamma \lambda )}{(1-\lambda )^\gamma }\ . \end{equation} (79) Here, the extra multiplication factor, i.e. (1 − γλ)/(1 − λ)γ, is always between 0 and 1 by Bernoulli's inequality.6 In other words, the boundary is always ‘lower’ (in the sense of Fig. 1) compared with the λ = 0 case and any parameter pair (γ, s) within the type C regime of (the original) Fig. 1 will always remain in the ‘no-root zone’ as λ varies. Therefore, the MSCC loop disappears completely, no matter how the magnetic parameter h is chosen. Secondly, for λ ≥ 1/γ, the LHS of equation (78) is always negative and no solution can exist. This immediately leads to the conclusion that, for (γ, s) within the type A regime of Fig. 1, the point x = 0, v = 0 is actually a ‘kink’ on the MSCC loop, whose tangents on the two sides (projected on to the −v versus x plane) are given by x = 0 and v = x/γ, respectively. In fact, it is possible to carry out a more detailed analysis of the asymptotic behaviour of the MSCC near x = 0 and the results are summarized in Appendix A. Fig. 3 shows the variation of MSCC loops with a sequence of different s values, holding both values of γ = 1.4 and h = 4 fixed. Clearly, the loop enlarges as s > 0 decreases and its two intersections with the v = 0 plane move in two opposite directions, with the intersection at larger x approaching infinity as s goes to 0. Thus, the portion of the loop next to the intersection with smaller x is naturally identified as the counterpart of the MSCC in the classical regime of Newtonian gravity (Lou & Shen 2004; Wang & Lou 2007, 2008). We emphasize that this solution behaviour is only characteristic for γ > 4/3. For γ < 4/3, the only intersection moves towards right to some limiting position at finite x as s → 0 and is not shown in Fig. 3 here. Figure 3. Open in new tabDownload slide The projection of MSCCs on to the −v versus x plane for different values of s with fixed values of γ = 1.4 and h = 4. Each MSCC intersects with the half-plane v = 0, x > 0 twice. As s decreases, the MSCC loop enlarges systematically and the intersection with larger x tends to infinity, leaving in effect only one intersection. Portions of the MSCC near that remaining intersection (of smaller x > 0) corresponds to the classical behaviour of Newtonian gravity as parameter s → 0. Figure 3. Open in new tabDownload slide The projection of MSCCs on to the −v versus x plane for different values of s with fixed values of γ = 1.4 and h = 4. Each MSCC intersects with the half-plane v = 0, x > 0 twice. As s decreases, the MSCC loop enlarges systematically and the intersection with larger x tends to infinity, leaving in effect only one intersection. Portions of the MSCC near that remaining intersection (of smaller x > 0) corresponds to the classical behaviour of Newtonian gravity as parameter s → 0. In order to show the influence of random transverse magnetic fields on an MSCC loop, the loop projections on to the −v versus x plane, with various values of h given two fixed values of γ = 1.4 and s = 0.2929 are shown in Fig. 4. As the magnetic field strength increases, the MSCC loop expands and seems to be ‘repelled’ from the origin x = v = 0. This is consistent with equation (77), which indicates that the intersection of the MSCC loop with any half-plane v = λx moves towards larger x as magnetic parameter h increases. The two dash–dotted straight lines emanating from the origin x = v = 0 form envelopes to this family of loops, which correspond to the two special values of λ in equation (78) that allow for a double root in A΄ to emerge. Qualitatively similar behaviour with variation of h is observed for γ < 4/3, despite the fact that the loop always touches the origin x = v = 0. Figure 4. Open in new tabDownload slide The projection of MSCCs on to the −v versus x plane for different values of magnetic parameter h with two fixed values of γ = 1.4 and s = 0.2929. Each MSCC intersects twice with the half-plane v = 0, x > 0. As h increases, the entire MSCC expands and moves away from the origin x = v = 0. Specifically, its two intersections with the v = 0 plane moves towards larger x values according to two expressions in equation (76). The two dash–dotted straight lines emanating from the origin form two envelopes to this family of MSCC loops. Figure 4. Open in new tabDownload slide The projection of MSCCs on to the −v versus x plane for different values of magnetic parameter h with two fixed values of γ = 1.4 and s = 0.2929. Each MSCC intersects twice with the half-plane v = 0, x > 0. As h increases, the entire MSCC expands and moves away from the origin x = v = 0. Specifically, its two intersections with the v = 0 plane moves towards larger x values according to two expressions in equation (76). The two dash–dotted straight lines emanating from the origin form two envelopes to this family of MSCC loops. 5.2 Types of singularities on the MSCC In the foregoing analysis, we mainly focused on the geometrical shape of the MSCC for different values of three pertinent parameters γ, s and h. Of equal importance is the study of the types of singularities along the MSCC, as these different singularity types impose important restrictions on the behaviours that global MHD solutions may take, especially for those solutions that cross the MSCC smoothly. Moreover, extra care is needed when dealing with numerical integrations near singularities of certain types across the MSCC. By applying the L'Hôpital rule to ODE (16) and eliminating α΄ ≡ dα/dx using ODE (8), it is found that for points on the MSCC, the first derivative |$v^{\prime }_{\ast }$| of v with respect to x along the solution satisfies the following quadratic equation: \begin{equation} \mathcal {A}(v^{\prime }_{\ast })^2+(\mathcal {B}-\mathcal {C})v^{\prime }_{\ast } -\mathcal {D}=0, \end{equation} (80) where the lower asterisk highlights the derivative being taken at a point right on the MSCC and the four coefficients |$\mathcal {A},\ \mathcal {B},\ \mathcal {C}$| and |$\mathcal {D}$| are explicitly defined below in order \begin{eqnarray} \mathcal {A} \equiv -\frac{(\gamma -2)}{(x-v)} \frac{\gamma \beta }{\alpha } - 3(x-v)\ ; \end{eqnarray} (81) \begin{eqnarray} \mathcal {B} \equiv 2(x-v) + (2\gamma -2)\left(\frac{1}{x} -\frac{1}{x-v}\right)\frac{\gamma \beta }{\alpha }\ ; \end{eqnarray} (82) \begin{eqnarray} \mathcal {C} \equiv -(2\gamma -2)\left( \frac{1}{x} - \frac{1}{x-v} \right)\frac{\gamma \beta }{\alpha }\ ; \end{eqnarray} (83) \begin{eqnarray} \mathcal {D} &\equiv & \frac{2\gamma \beta }{\alpha (x-v)}\left[ \frac{-4\gamma ^2+5\gamma -2}{\gamma }\left( \frac{x-v}{x} \right) \right.\nonumber\\ && \left. +\,(2\gamma -1)\left(\frac{x-v}{x}\right)^2 + \frac{2(\gamma -1)^2}{\gamma } \right] + \frac{2mv}{(x-sm)^3}\ . \end{eqnarray} (84) Quadratic equation (80) tells us that for any point on the MSCC, at most two directions are allowed for solutions to cross it smoothly. We denote the eigenvalues associated with these two directions by λ1 and λ2, which are easily found to be \begin{equation} \lambda _{1,\ 2} = \mathcal {A}v^{\prime }_{{\ast }1,\ 2} + \mathcal {B}, \end{equation} (85) where |$v^{\prime }_{{\ast }1,\ 2}$| are the two solutions of equation (80). Thus, we have \begin{equation} \lambda _1\lambda _2 = -\mathcal {A}\mathcal {D} + \mathcal {B}\mathcal {C}\ , \qquad \quad \lambda _1 + \lambda _2 = \mathcal {B} + \mathcal {C}\ . \end{equation} (86) The type of singularity is determined by the signs of λ1λ2 and Δ ≡ (λ1 + λ2)2 − 4λ1λ2. If λ1λ2 < 0, the point is a saddle. Meanwhile, if λ1λ2 > 0 and Δ ≥ 0, the point is a node. Otherwise, it is a focus (see Whitworth & Summers 1985 for details). As remarked in Lou & Shen (2004), numerical integration away from a saddle point along any one of the two eigendirections is stable, while numerical integration away from a node is neutrally stable along the eigendirection with larger eigenvalue, and unstable (due to non-uniqueness of solutions) along the eigendirection with smaller eigenvalue. No solution can actually arrive at any focus point on the MSCC. By explicit computation and elementary calculus, one can show that the expression for λ1λ2 contains an overall factor of v in it and that Δ is always positive when v = 0 (see Appendix B for more specific details). This at least tells us that an intersection of the MSCC with the V = 0, x > 0 half-plane is generally a bifurcation point separating nodes and saddle points (see e.g. Whitworth & Summers 1985; Lou & Shen 2004 for similar analyses and discussions in the Newtonian and isothermal hydrodynamic case without magnetic field). Besides these intersections, our extensive numerical explorations indicate that there may be other bifurcation points on the MSCC under certain choices of γ and s parameters, some of which are associated with the emergence of foci for larger values of γ. Fig. 5 demonstrates this point explicitly by plotting the distribution of different types of singularities along the MSCC. Figure 5. Open in new tabDownload slide Different types of singularities along the MSCC loops. The top panel is drawn for s = 0.2929, γ = 1.2 and h = 5.0. The middle and bottom panels both have s = 0.2 and h = 5.0 and the γ values are 1.7 and 1.8, respectively. The solid, dashed and dotted curves stand for the portions of the MSCC that consist of saddle points, nodes and foci, respectively, as classified in the main text. The dash–dotted line is the ZML, drawn for reference. Note that the intersections of the MSCC with the v = 0 plane at positive x > 0 always separates portions of saddles and nodes. Comparison of the middle and bottom panels illustrates the dramatic change of distributions of singularity types as γ value is slightly increased. Our numerical exploration indicates that the occurrence of foci are typical for a combination of large γ values and relatively small s values, i.e. towards the lower right of Fig. 1. However, we did not attempt a rigorous mathematical proof for this claim. Figure 5. Open in new tabDownload slide Different types of singularities along the MSCC loops. The top panel is drawn for s = 0.2929, γ = 1.2 and h = 5.0. The middle and bottom panels both have s = 0.2 and h = 5.0 and the γ values are 1.7 and 1.8, respectively. The solid, dashed and dotted curves stand for the portions of the MSCC that consist of saddle points, nodes and foci, respectively, as classified in the main text. The dash–dotted line is the ZML, drawn for reference. Note that the intersections of the MSCC with the v = 0 plane at positive x > 0 always separates portions of saddles and nodes. Comparison of the middle and bottom panels illustrates the dramatic change of distributions of singularity types as γ value is slightly increased. Our numerical exploration indicates that the occurrence of foci are typical for a combination of large γ values and relatively small s values, i.e. towards the lower right of Fig. 1. However, we did not attempt a rigorous mathematical proof for this claim. In this section, we have focused on the singularity of X = 0, and especially the MSCC. Besides this, there exists another kind of singularity that deserves attention, namely the Paczynski–Wiita singularity of x − sm = 0 for the expanding event horizon. Both functionals Y and Z tend to infinity at this singularity. Fortunately, this important singularity sets the physical limit such that physical MHD solutions can never actually enter it from outside. The proof of this statement is presented in Appendix C. 6 NUMERICAL CONSTRUCTION OF VARIOUS SELF-SIMILAR MHD SOLUTIONS Having examined properties of the MSCC, we now focus on the numerical construction of global self-similar MHD solutions. For this purpose, a fourth-order Runge–Kutta scheme is implemented to integrate two coupled non-linear ODEs (15) and (16) from consistently specified initial values or analytic asymptotic solutions. These initial values or asymptotic solutions are chosen so that the resulting solution falls into certain specific categories that we would investigate and the difficulty of their choice varies considerably among the available solution categories. The pertinent details will be presented as we examine each category separately. Moreover, by analysing bunches of solutions within the same category under various initial conditions and parameter values, a more concrete understanding of MHD solution behaviours could be gained. 6.1 EWCSs and other MHD solutions The numerical solution construction procedure is fairly straightforward based on our past experience, as long as a solution does not hit the MSCS X(x, α, v) = 0. Recall that an asymptotic solution at large x, with v(x) and α(x) approaching V0 and A0/x2, respectively, is analytically represented by asymptotic solutions (23) and (24). By choosing a sufficiently large x, they provide a consistent starting point for inward integration towards small x. Fig. 6 shows several MHD solutions obtained with this numerical procedure, with three parameter values of s = 0.2, γ = 1.7 and h = 5.0 specified (that fall in the type B regime of Fig. 1). One prominent feature of this bunch of MHD solutions is the pair shown with bold solid curves. These are the EWCSs first studied by Shu (1977) for an isothermal gas sphere under the Newtonian gravity. Chiueh & Chou (1994) extended this study to include random transverse magnetic fields but with insufficient solution exploration and no magnetic tension force, while Wang & Lou (2008) further generalized such MHD formalism with both magnetic pressure and tension forces in a quasi-spherical geometry using a GP EoS. The analysis of spherical relativistic EWCSs with the Paczynski–Wiita gravity was pioneered by Lian & Lou (2014), where it was shown that, for each intersection of the SCC with the v = 0 plane, there exists an EWCS whose expanding wavefront sits right at that very intersection. This outgoing wavefront separates the outer static configuration, given by solution (20) in the absence of magnetic field, and the inner infalling region expanding at a constant speed. The inner infalling gas ultimately gets its mass continuously swallowed by a dynamically accreting central Schwarzschild BH whose mass grows with increasing time. Near the growing event horizon, the dynamic solution takes on the ‘free-fall’ asymptotic solutions (35) and (36). The Schwarzschild radius and the mass of the central accreting BH (event horizon) also increase linearly with increasing time t and the dimensionless mass accretion rate m0 can be determined from the numerical solution. The EWCSs are of great theoretical importance, since it captures the basic elements of the spontaneous, quasi-spherically symmetric collapse of an otherwise static GP gas envelope and predicts the mass growth of a central BH given the pertinent parameters. Here in Fig. 6, there are two EWCSs corresponding to the two intersections of the MSCC with the v = 0 plane. The one with a more massive central BH has a higher value for A0 (related to the square of sound speed or random velocity dispersion in proper contexts), such that the outer static spherical envelope is denser in comparison with the other smaller A0 solution in the pair of roots. Figure 6. Open in new tabDownload slide Several numerical self-similar MHD solutions with s = 0.2, γ = 1.7 and h = 5.0. The asymptotic velocity and density parameters V0 and A0, as well as the estimated mass accretion rates m0 near the central BH, are marked along with each solution explicitly. The MSCC is shown as a dash–dotted curve, and the vertical dotted straight lines indicate the estimated positions of event horizons (corresponding to the growing Schwarzschild radii). The EWCSs are shown in two boldface solid curves and solutions with inflows (upper right) or outflows (lower left) at large x are shown by light solid curves. The dashed curves are solutions with zero asymptotic speed at large x but different asymptotic density parameters A0 from those EWCSs (the legend on the lower left is for one dashed curve). The asymptotic behaviour of these solutions near the event horizon satisfies equations (35) and (36). As pointed out in Section 3.3, a simple empirical relation (40) should be invoked for the divergent v near the event horizon to accommodate special relativistic effects so that the true radial flow speed remains less than the vacuum speed of light c. Figure 6. Open in new tabDownload slide Several numerical self-similar MHD solutions with s = 0.2, γ = 1.7 and h = 5.0. The asymptotic velocity and density parameters V0 and A0, as well as the estimated mass accretion rates m0 near the central BH, are marked along with each solution explicitly. The MSCC is shown as a dash–dotted curve, and the vertical dotted straight lines indicate the estimated positions of event horizons (corresponding to the growing Schwarzschild radii). The EWCSs are shown in two boldface solid curves and solutions with inflows (upper right) or outflows (lower left) at large x are shown by light solid curves. The dashed curves are solutions with zero asymptotic speed at large x but different asymptotic density parameters A0 from those EWCSs (the legend on the lower left is for one dashed curve). The asymptotic behaviour of these solutions near the event horizon satisfies equations (35) and (36). As pointed out in Section 3.3, a simple empirical relation (40) should be invoked for the divergent v near the event horizon to accommodate special relativistic effects so that the true radial flow speed remains less than the vacuum speed of light c. The four light solid curves of Fig. 6 are solutions with finite non-zero values of V0 as x → +∞. Two of them correspond to the asymptotic inflows at large x (upper right) and the other two correspond to the asymptotic outflows (lower left). The A0 values of the two inflow solutions are designed to be the same as the ‘dense’ EWCS (i.e. the EWCS with a more massive central BH), while that of the two outflow solutions are equal to that of the ‘dilute’ EWCS. It can be readily seen that as the asymptotic speed |V0| decreases, the inflow solution gets more and more close to the dense EWCS. In other words, the dense EWCS can be regarded as the limiting case for those inflow solutions with diminishing inflow velocity. Similarly, the dilute EWCS can be understood as the limit of a bunch of outflow solutions with vanishing outflow velocity asymptotically. Also shown in Fig. 6 are two dashed-curve solutions. They are constructed with V0 = 0 but different A values from the EWCSs. They have asymptotic behaviours closely mimicking the EWCSs both at large x and near the event horizon, but are smooth and do not contain clear-cut wavefronts of discontinuous first derivative separating static outer envelopes and dynamic collapsing inner regions. Here, it would be appropriate to elucidate a subtle point concerning Fig. 6 and thus to avoid a potential conceptual confusion. It looks as if some of the solutions in Fig. 6 ‘intersect’ the MSCC shown by the loop-like dash–dotted curve. However, this is just an apparent visual effect resulting from the projection of the three-dimensional MSCC in (x, α, v) spanned space on to the −v versus x plane (i.e. the α = 0 plane). In reality, none of these solutions actually ‘cross’ the MSCC, as can be specifically checked by plotting also the projection of everything on to the α versus x plane (i.e. the v = 0 plane), nor do these solutions even touch the MSCS X = 0 in the (x, α, v) spanned space, except for the EWCSs whose expanding wavefronts lie right on that MSCS. Afterall, this is the very basis for our inward integration scheme to be valid. Lian & Lou (2014) was inaccurate on this by saying that some solutions cross the SCC twice in their Fig. 4. The presence of the MSCS X = 0 turns out to restrict the values that A0 and V0 can take for constructing globally smooth MHD solutions. Certain combinations of A0 and V0 will result in a solution that eventually runs towards the MSCS and is therefore unfeasible unless shock waves are properly inserted. Such an (A0, V0) parameter pair will be referred to as ‘forbidden’, otherwise ‘allowed’ in constructing globally smooth MHD solutions. Alternatively, there is, in fact, another possibility that the solution may happen to meet the MSCC. If the meeting point does not lie on the ZML, the solution then could, in principle, be continued to pierce through the MSCS X = 0, but analyticity is not necessarily guaranteed at the crossing point. Fig. 7 shows the result of our attempt to identify the allowed, forbidden and ‘MSCC-destined’ values of (A0, V0) for constructing globally smooth MHD solutions. Such extensive calculations are done with the same parameters γ, s and h as in Fig. 6. The grey-coloured region represents forbidden pairs of A0 and V0, with which the solution terminates on the MSCS X = 0. The light-coloured region, on the other hand, are the allowed values for (A0, V0) pairs. Besides them, there exists a third region, coloured dark and labelled ‘MSCC’ or ‘ZML’. This is the region where the solution heads on to a point on the MSCS X = 0 that lies sufficiently close to the MSCC loop or the ZML.7 The two special (A0, V0) pairs corresponding to the EWCSs lie on the boundary between the light and dark-coloured regions and are specifically marked out by the two plus signs in Fig. 7. From Fig. 7 and in terms of constructing globally smooth MHD solutions, it is clear that a considerable domain of (A0, V0) pairs is ruled out by the existence of the MSCS X = 0. Also it is interesting to note that despite the one-dimensional nature of the MSCC itself, the set of (A0, V0) pairs that lead to a meeting with such a one-dimensional curve is actually two-dimensional. That is, the MSCC has a certain ‘region of attraction’ of positive areas. In this sense, crossing the MSCC is not as rare an event as one might first think. Figure 7. Open in new tabDownload slide A plot of the allowed (light) and forbidden (dark) regions, corresponding to the existence and non-existence of globally smooth MHD solutions with asymptotic parameters of A0 and V0, given s = 0.2, γ = 1.7 and h = 5.0. Besides those two regions, there are other dark-coloured regions labelled ‘MSCC’ and ‘ZML’. They encompass those (A0, V0) parameter pairs whose corresponding solutions hit the MSCS X = 0 at a point sufficiently close to the MSCC loop or the ZML. The two (A0, V0) parameter pairs for the EWCSs lie on the boundary between the light and dark-coloured regions, as shown by the two plus signs on the horizontal line V0 = 0. Figure 7. Open in new tabDownload slide A plot of the allowed (light) and forbidden (dark) regions, corresponding to the existence and non-existence of globally smooth MHD solutions with asymptotic parameters of A0 and V0, given s = 0.2, γ = 1.7 and h = 5.0. Besides those two regions, there are other dark-coloured regions labelled ‘MSCC’ and ‘ZML’. They encompass those (A0, V0) parameter pairs whose corresponding solutions hit the MSCS X = 0 at a point sufficiently close to the MSCC loop or the ZML. The two (A0, V0) parameter pairs for the EWCSs lie on the boundary between the light and dark-coloured regions, as shown by the two plus signs on the horizontal line V0 = 0. As an example of illustration, Fig. 8 highlights the influence of random transverse magnetic field (h parameter) on the dynamic mass accretion rates towards central Schwarzschild BHs. Fig. 8 is prepared by constructing MHD EWCSs for different values of h given two parameters s and γ fixed at 0.2 and 1.7, respectively, and then deriving the m0 parameter for each of these solutions by fitting analytic asymptotic free-fall solutions (35) and (36). From Fig. 8, it is evident that for both of the MHD EWCSs, m0 is a monotonically increasing function of h. As h increases from 0 to 10, the mass accretion rate of the dilute and dense MHD EWCSs increase by a factor of almost 3 and 1/3, respectively. Recall from Section 3.1 that all solutions belonging to the same class (either dilute or dense) have the same α(x) profile in the outer static envelope, yet they differ in the strength of the underlying force-free magnetic field as characterized by the magnetic parameter h. According to magnetic flux conservation integral (11), the magnetic field strength in the outer envelope is therefore directly proportional to h within each class and we conclude that a random transverse magnetic field of greater strength leads to higher dynamic mass accretion towards the central Schwarzschild BH. Figure 8. Open in new tabDownload slide The MHD EWCS mass accretion rate m0 is plotted as a function of the magnetic parameter h, given fixed values s = 0.2 and γ = 1.7. ‘Dense’ (‘dilute’) MHD EWCS corresponds to an expansion wavefront sitting at larger (smaller) x > 0. The curve for dilute MHD EWCS is multiplied by a factor of 20 for an easier visual comparison here. Evidently, the mass accretion rate grows as the magnetic field gets stronger. Figure 8. Open in new tabDownload slide The MHD EWCS mass accretion rate m0 is plotted as a function of the magnetic parameter h, given fixed values s = 0.2 and γ = 1.7. ‘Dense’ (‘dilute’) MHD EWCS corresponds to an expansion wavefront sitting at larger (smaller) x > 0. The curve for dilute MHD EWCS is multiplied by a factor of 20 for an easier visual comparison here. Evidently, the mass accretion rate grows as the magnetic field gets stronger. In Fig. 9, we exhibit several numerical self-similar MHD solutions with a parameter set within the type A regime (i.e. γ < 4/3 regime) in Fig. 1, i.e. s = 0.2929, γ = 1.2 and h = 5.0. Both the −v versus x (i.e. α = 0 projection) and the α versus x (i.e. v = 0 projection) perspectives are shown simultaneously for a convenient comparison. In this case, the MSCC intersects the v = 0 plane only once and the MHD EWCS is unique, as shown by the boldface curve a. Above this curve a, three other solutions (viz., curves b, c and d) are displayed (to the upper right part in the −v versus x panel). Two of them have the same asymptotic mass parameter A0 as the MHD EWCS but with negative (inflow) asymptotic speeds and are shown by two light solid curves (viz., b and c). The other solution, shown as the dashed curve d, has zero asymptotic flow speed at large x and its A0 value is greater than that of the MHD EWCS. Again as in Fig. 6, the MHD EWCS may be regarded as the limiting case of a bunch of solutions, of which b, c and d are representatives. Also shown in Fig. 9 are the three MHD solutions with central expanding voids (i.e. ending at the ZML) instead of central BHs, labelled by e, f and g. From the x versus α perspective, it is clear that they do not cross the MSCC and generally have much smaller mass densities over the entire range of x than their BH counterparts. Their mass density profiles gradually rise from zero as x goes from infinity towards the centre, whereas the radial outflow velocities remain almost unchanged. Below a certain x value, however, the outflow slows down drastically, accompanied by a downturn in mass density. Both the outflow velocity and density continue to decrease until the solution arrives at the void rim. The asymptotic behaviour in the vicinity of the void rim is precisely that described by the solutions (41)–(43). Figure 9. Open in new tabDownload slide Several numerical self-similar MHD solutions with s = 0.2929, γ = 1.2 and h = 5.0. Both the −v versus x (projection on to α = 0 plane) and the α versus x (projection on to v = 0 plane) perspectives are shown. The velocity and mass parameters V0 and A0 for the large-x asymptotic solution, as well as the derived dynamic mass accretion rates m0 towards the central BH or the radii x1 of central voids, are labelled along each solution. The MSCC is shown in dash–dotted curve and the vertical dotted straight lines in the upper panel indicate the estimated positions of event horizons in expansion. The EWCS a is shown in boldface solid curve a and solutions with inflows at large x are shown in light solid curves above the EWCS in the upper panel marked with b and c. The dashed curve d is a solution with asymptotic speed V0 = 0 and asymptotic mass parameter A0 larger than that of the EWCS. The lowest three curves (e, f, g) in the upper panel represent outflow solutions that contain central voids of different radii in expansion. Figure 9. Open in new tabDownload slide Several numerical self-similar MHD solutions with s = 0.2929, γ = 1.2 and h = 5.0. Both the −v versus x (projection on to α = 0 plane) and the α versus x (projection on to v = 0 plane) perspectives are shown. The velocity and mass parameters V0 and A0 for the large-x asymptotic solution, as well as the derived dynamic mass accretion rates m0 towards the central BH or the radii x1 of central voids, are labelled along each solution. The MSCC is shown in dash–dotted curve and the vertical dotted straight lines in the upper panel indicate the estimated positions of event horizons in expansion. The EWCS a is shown in boldface solid curve a and solutions with inflows at large x are shown in light solid curves above the EWCS in the upper panel marked with b and c. The dashed curve d is a solution with asymptotic speed V0 = 0 and asymptotic mass parameter A0 larger than that of the EWCS. The lowest three curves (e, f, g) in the upper panel represent outflow solutions that contain central voids of different radii in expansion. A diagram similar to Fig. 7 can be drawn for the present set of parameters, to show the existence of globally smooth solutions as depending on parameters A0 and V0 that control asymptotic behaviours at large x. The result is displayed in Fig. 10 that is especially interesting in contrast to Fig. 7. Apparently, the major difference between these two figures is that the light-coloured ‘channel’ lying between the two regions labelled ‘ZML’ and ‘MSCC’ in Fig. 7 seems to have been ‘closed’ in Fig. 10. As it turns out, this pattern alteration is closely related to the shape change in the MSCC: the ‘channel’ appears when the MSCC loop and the ZML are detached from each other and closes when they meet. The detached MSCC loop allows for globally smooth solutions with asymptotic radial outflow and small central BHs (e.g. the lowest two solutions shown in Fig. 6), as well as the dilute MHD EWCS that is considered as their limiting case. It is these solutions that form the light-coloured ‘channel’ in Fig. 7. The left ‘bank’ of the channel consists of solutions extending to x = 0, whose asymptotic behaviour there is described by the ‘generalized Larson-Penston solution’ in Section 3.6 (i.e. they have vanishing v and finite α at x = 0). In a loose sense, these solutions can also be regarded as a limiting case of void solutions with vanishing void radius. On the other hand, when the MSCC loop touches the ZML, these aforementioned solutions cease to exist and a lower limit arises for the size of voids that can occur in globally smooth MHD solutions, under any given asymptotic outflow speed V > 0. The MHD solutions with smallest voids are represented in Fig. 10 by the borderline separating the two regions labelled ‘ZML’ and ‘MSCC’, where the ‘channel’ lies originally. Figure 10. Open in new tabDownload slide (Non-)existence of globally smooth MHD solutions whose asymptotic behaviour at large x is characterized by A0 and V0 with s = 0.2929, γ = 1.2 and h = 5.0 fixed. The labels of various regions here bear the same meanings as those in Fig. 7 (see the caption there for details). The only EWCS here is marked out by the plus sign on the horizontal line V0 = 0. Unlike Fig. 7, the regions labelled ‘ZML’ and ‘MSCC’ now border each other; this is caused by the morphology difference in the MSCC. Figure 10. Open in new tabDownload slide (Non-)existence of globally smooth MHD solutions whose asymptotic behaviour at large x is characterized by A0 and V0 with s = 0.2929, γ = 1.2 and h = 5.0 fixed. The labels of various regions here bear the same meanings as those in Fig. 7 (see the caption there for details). The only EWCS here is marked out by the plus sign on the horizontal line V0 = 0. Unlike Fig. 7, the regions labelled ‘ZML’ and ‘MSCC’ now border each other; this is caused by the morphology difference in the MSCC. Fig. 11 displays several numerical MHD solutions with s = 4.0 and γ = 1.4. This set of parameters is selected from the type C regime in Fig. 1 and can host no MHD EWCS. Nevertheless, we fix A0 = 0.18 and plot pairs of solutions by adjusting the value V0 for h = 4.0 and 0, respectively. Comparison of solutions with the same value of h demonstrates that a smaller asymptotic inflow speed or a larger outflow speed will, in general, lead to smaller central BHs and voids occur when the asymptotic outflow speed is increased beyond a certain threshold. By further comparing the solutions of h = 4 with their h = 0 counterparts, it is readily seen that the presence of a random transverse magnetic field increases the radius of central BHs and decreases the size of central voids. That is, magnetization favours the dynamic formation of BHs. In one particular case of V0 = 1.00, a BH is actually ‘created’ by ‘turning on’ random transverse magnetic fields. Figure 11. Open in new tabDownload slide Several numerical self-similar MHD solutions with s = 4.0 and γ = 1.4. The solid-curve and dashed-curve solutions are constructed, given h = 4.0 and 0, respectively. Asymptotic radial velocities V0 are labelled for each pair of solution with the same A0 = 0.18 for all. Here, the estimated dynamic mass accretion rates m0 or void radii x1 are computed for h = 4.0 solutions only. Vertical dotted straight lines indicate the estimated positions of event horizons and are also drawn for h = 4.0 solutions only. The ZML is shown as the straight dash–dotted line in the lower left part. Figure 11. Open in new tabDownload slide Several numerical self-similar MHD solutions with s = 4.0 and γ = 1.4. The solid-curve and dashed-curve solutions are constructed, given h = 4.0 and 0, respectively. Asymptotic radial velocities V0 are labelled for each pair of solution with the same A0 = 0.18 for all. Here, the estimated dynamic mass accretion rates m0 or void radii x1 are computed for h = 4.0 solutions only. Vertical dotted straight lines indicate the estimated positions of event horizons and are also drawn for h = 4.0 solutions only. The ZML is shown as the straight dash–dotted line in the lower left part. 6.2 Smooth MHD solutions across the MSCC We now carry out the general idea for constructing numerical MHD self-similar solutions one step further and examine solutions that cross the MSCC smoothly. The procedure for obtaining local asymptotic smooth solutions near the MSCC has been outlined in Section 5.2 and can be utilized as a starting point for systematic numerical integration. Previously, using an isothermal model with Newtonian gravity, Lou & Shen (2004) constructed EECC hydrodynamic solutions that cross the SCC smoothly twice and proposed their astrophysical applications in describing the behaviour of AGB or post-AGB stars and PNe. Lou & Zhai (2009, 2010) proposed shock models for diffuse X-ray emissions from dynamic PNe NGC 40 and NGC 7662 with variation time-scales of several decades (see also the recent report on proto PN SAO 244567 in Stingray Nebula by Reindl et al. 2014). Wang & Lou (2008) generalized these solutions to the GP formalism with random transverse magnetic fields in a quasi-spherical geometry. In the present context of Paczynski–Wiita gravity, a similar kind of MHD solution can be constructed for 1 ≤ γ < 4/3 (i.e. type A regime in Fig. 1). These smooth GP MHD solutions involve a central mass-accreting BH of extremely small mass or Schwarzschild radius and an MHD envelope that is either contracting (V0 < 0) or expanding (V0 > 0) asymptotically at large x. For a given set of parameters γ, s and h that belong to the type A regime, all smooth MHD solutions that cross the MSCC twice can be arranged into an infinite sequence. For instance, with s = 0.2929, γ = 1.2 and h = 0.1, the first two solutions of the sequence are shown in Fig. 12 (solid curves a and b). These MHD solutions are obtained using the phase-diagram matching scheme (e.g. Lou & Shen 2004). In general, the phase-diagram matching scheme allows one to construct global solutions out of two distinct one-parameter families of local solutions. The two families of local solutions are first integrated to some common meeting plane (also known as phase plane) x = xF to form two phase curves on that plane. Every intersection of the two phase curves then correspond to a legitimate smooth transition between two local solutions, one from each family, leading to a global MHD solution. The phase-diagram matching scheme is numerically stable in the sense that even though numerical integration errors may slightly displace the phase curves off their actual positions, the number and distribution pattern of their intersections will not change in general. Therefore, it is a very reliable method for showing the existence of global MHD solutions of various properties. In our present case, we set xF = 0.7 and take the two families as consisting of local solutions that cross the outer rim of the MSCC smoothly at x*1 and x*2, respectively, along the eigendirection with the larger eigenvalue, where x*1 < xF < x*2. The resulting phase diagram is demonstrated in Fig. 13, with the values of x*1 or x*2 labelled for data points along the phase curves. The phase curve from the first family shows an inward spiral pattern that is common in similar kinds of construction (e.g. Whitworth & Summers 1985; Lou & Shen 2004) and intersects the second phase curve infinitely many times as it converges (under x*1 → 0) to a specific limiting point P0(x = xF, α = 2.041, v = 0) that also belongs to the second phase curve. Correspondingly, the global smooth MHD solutions arising from those phase curve intersections also arrange themselves into a sequence that converges pointwise to the global solution passing through the point P0. It is interesting to observe that the limiting solution is exactly the global static solution described in Section 3.1, which, unlike all solutions in the sequence, lacks a BH at its centre. Figure 12. Open in new tabDownload slide Several MHD solutions that cross the MSCC smoothly, obtained with pertinent parameters s = 0.2929, γ = 1.2 and h = 0.1. The solid curves a and b are the first two in the sequence of MHD solutions that cross the MSCC smoothly twice and harbour a small central BH. The x-coordinates of both MSCC intersections, denoted by x*1 and x*2, are marked in the upper panel. The dimensionless mass accretion rates m0, the event horizon radii x0 and the large-x asymptotic parameters V0 and A0 are shown in the lower panel. The MHD solution a is an EECC solution, while MHD solution b is an envelope-contraction-core-collapse (ECCC) solution. Note that the BH radii are orders of magnitude smaller than those of the EWCSs and a change of variable (i.e. setting ln x, rather than x, to be our independent variable in numerical integration) has to be done to accommodate to this situation. The dashed curve a΄ is the first element in the sequence of solutions that assume LX-1 type central behaviour and cross the MSCC smoothly once. All these MHD solutions are obtained by phase diagram matching procedure. Figure 12. Open in new tabDownload slide Several MHD solutions that cross the MSCC smoothly, obtained with pertinent parameters s = 0.2929, γ = 1.2 and h = 0.1. The solid curves a and b are the first two in the sequence of MHD solutions that cross the MSCC smoothly twice and harbour a small central BH. The x-coordinates of both MSCC intersections, denoted by x*1 and x*2, are marked in the upper panel. The dimensionless mass accretion rates m0, the event horizon radii x0 and the large-x asymptotic parameters V0 and A0 are shown in the lower panel. The MHD solution a is an EECC solution, while MHD solution b is an envelope-contraction-core-collapse (ECCC) solution. Note that the BH radii are orders of magnitude smaller than those of the EWCSs and a change of variable (i.e. setting ln x, rather than x, to be our independent variable in numerical integration) has to be done to accommodate to this situation. The dashed curve a΄ is the first element in the sequence of solutions that assume LX-1 type central behaviour and cross the MSCC smoothly once. All these MHD solutions are obtained by phase diagram matching procedure. Figure 13. Open in new tabDownload slide Phase diagram for the construction of MHD solutions in Fig. 12. The meeting plane is chosen at xF = 0.7. We first select a series of points on the upper limb of MSCC (i.e. the portion with smaller v, see Fig. 12) with x < xF (all these points are saddles) and integrate outwards from them along one of the two eigendirections. The intersections of these solutions with the meeting plane are marked by open circles and these circles are connected to form the spiraling phase curve. Next, we choose a series of points on the upper limb of MSCC with x > xF (these points are either saddles or nodes) and perform an inward integration from them along the eigendirection of larger eigenvalue (thus the numerical integrations remain stable). The intersections of these solutions with the meeting plane are connected to form the second phase curve. Numbers associated with data points are the x values of points on the MSCC from which the integration started. It is not hard to imagine that if the spiral is continued to even smaller initial values of x, it will intersect the other phase curve infinitely many times before tending to a limiting position on the meeting plane. As a result, there is a countably infinite sequence of discrete matched solutions. Fig. 12 shows the first two of them. Note that the limiting position of the phase spiral correspond to the global magnetostatic solution mentioned at the beginning of Section 3. Figure 13. Open in new tabDownload slide Phase diagram for the construction of MHD solutions in Fig. 12. The meeting plane is chosen at xF = 0.7. We first select a series of points on the upper limb of MSCC (i.e. the portion with smaller v, see Fig. 12) with x < xF (all these points are saddles) and integrate outwards from them along one of the two eigendirections. The intersections of these solutions with the meeting plane are marked by open circles and these circles are connected to form the spiraling phase curve. Next, we choose a series of points on the upper limb of MSCC with x > xF (these points are either saddles or nodes) and perform an inward integration from them along the eigendirection of larger eigenvalue (thus the numerical integrations remain stable). The intersections of these solutions with the meeting plane are connected to form the second phase curve. Numbers associated with data points are the x values of points on the MSCC from which the integration started. It is not hard to imagine that if the spiral is continued to even smaller initial values of x, it will intersect the other phase curve infinitely many times before tending to a limiting position on the meeting plane. As a result, there is a countably infinite sequence of discrete matched solutions. Fig. 12 shows the first two of them. Note that the limiting position of the phase spiral correspond to the global magnetostatic solution mentioned at the beginning of Section 3. MHD solutions that cross the MSCC smoothly twice is not the only special sequence of smooth global solutions in the type A regime of Fig. 1. As our numerical exploration reveals, there exists another infinite sequence of special solutions that closely resembles the previous one in every respect, except for central dynamic behaviours. This second sequence of MHD solutions cross the MSCC only once and satisfies LX-1 type asymptotic dynamic behaviour at the centre (see Section 3.6 for this novel asymptotic solution). The first solution (with the same pertinent parameters) of this sequence is illustrated in Fig. 12 by the dashed curve a΄. A phase diagram matching scheme, which replaces the first family of local solutions in our previous scheme by those that have LX-1-type central dynamic behaviour parametrized by the free coefficient T, is used for obtaining this sequence. It would offer more insight to consider the above two sequences of global smooth MHD solutions as just one sequence of solution pairs, by simply pairing the nth solution of both sequences (n = 1, 2, 3, …). One compelling reason for doing so is that the paired MHD solutions lie very close to each other except for very small x where they behave quite differently. In Fig. 12, for example, MHD solutions a and a΄ are almost identical at very large x and both feature an envelope expansion, whereas MHD solution b has an envelope contraction.8 By tracing the phase spiral and counting the intersections in order, we can continue constructing the paired MHD solution sequence. In general, envelope-expansion and envelope-contraction solution pairs occur alternately in the sequence. That is, the third pair again has an envelope expansion, while the fourth pair has an envelope contraction, and so forth. The expansion or contraction velocity diminishes as we go further along the sequence. All MHD solutions in the sequence harbour central BHs, but they are many orders of magnitude smaller in mass and size as compared with the BHs in EWCSs and are surrounded by very dense shells. For example, the second solution in the sequence has a BH of radius x0 = 4.0 × 10−13 and a maximum reduced density αmax = 1.0 × 1012 at x = 1.45x0, although α does abruptly decrease to zero near the event horizon as predicted by the asymptotic solution (36) (this cannot be seen in Fig. 12 due to the small event horizon radius). The third solution (not shown in Fig. 12) has x0 = 4.2 × 10−20 and αmax = 8.7 × 1018, which is about the density ratio of nuclear matter versus ordinary matter. For such unusually high densities, a change in the GP EoS would be expected. While idealized, these model MHD solutions reveal the following physical scenario. A small BH emerges at the centre dynamically accreting from a surrounding magnetized yet extremely dense compact layer that may be capable of bouncing infalling materials to drive an expanding MHD shock. In other words, the demise of a massive progenitor may appear as a supernova yet with a central BH with shredded clumps around (e.g. SN1987A might be a relevant case of interest; Lou & Lian 2012; Lian & Lou 2014). For more information, the dimensionless central BH mass accretion rate m0 as a function of magnetic parameter h is shown in Fig. 14. Unlike those EWCSs, an increasing magnetic field strength suppresses the growth of these BHs and the effect is almost exponential in h (note that this is a semilog plot). Figure 14. Open in new tabDownload slide BH mass accretion rate m0 as a function of magnetic parameter h, drawn for the first two MHD solutions that cross the MSCC twice. Other pertinent parameters are fixed at s = 0.2929 and γ = 1.2. As h increases, mass accretion rates of both MHD solutions decrease exponentially. Figure 14. Open in new tabDownload slide BH mass accretion rate m0 as a function of magnetic parameter h, drawn for the first two MHD solutions that cross the MSCC twice. Other pertinent parameters are fixed at s = 0.2929 and γ = 1.2. As h increases, mass accretion rates of both MHD solutions decrease exponentially. For type B regime parameters of Fig. 1, MHD solutions that cross the MSCC twice were not found in our numerical exploration. The counterparts to the aforementioned solutions in the type A regime of Fig. 1 appear to be those that have a central ‘quasi-magnetostatic’ configuration and cross the outer limb of the MSCC only once. Three such MHD solutions are displayed in Fig. 15. They are obtained by the phase-diagram matching of outward integrations initiated with an analytic ‘quasi-magnetostatic’ core solution (see Section 3.5) and inward integrations from the MSCC along stable (or neutrally stable) directions. The phase curve of the former set of MHD solutions again exhibits a spiral pattern, intersecting the latter set infinitely many times. Therefore, the solutions also form an infinite sequence, which is again found to exhibit alternating envelope expansion and contraction, with the global magnetostatic solution serving as the limiting case as x → 0. Self-similar sub-magnetosonic radial oscillations are found near the centre, where the velocity profiles show an undulating pattern. For each successive MHD solution in the sequence, the number of stagnation points increases by one. The magnitude of the oscillation decreases, as we go further along the MHD solution sequence. Figure 15. Open in new tabDownload slide The first three MHD solutions that have a central ‘quasi-magnetostatic’ configuration v ∼ UxQ and cross the MSCC smoothly once at its outer limb. The pertinent parameters are s = 0.2, γ = 1.7 and h = 0.1. The power index Q in this particular example is found to be 1.543 and different coefficients U are labelled with each solution. The MSCC crossing positions x* as well as the two large-x asymptotic velocity and mass parameters V0 and A0 are also indicated. Envelope expansion and contraction solutions again occur alternately in the present series. Also note that the velocity profiles near the centre show an undulating pattern, with each successive solution possessing one more stagnation point. Physically, these profiles represent radial sub-magnetosonic oscillations in a self-similar manner and the oscillatory pattern travels outwards. Figure 15. Open in new tabDownload slide The first three MHD solutions that have a central ‘quasi-magnetostatic’ configuration v ∼ UxQ and cross the MSCC smoothly once at its outer limb. The pertinent parameters are s = 0.2, γ = 1.7 and h = 0.1. The power index Q in this particular example is found to be 1.543 and different coefficients U are labelled with each solution. The MSCC crossing positions x* as well as the two large-x asymptotic velocity and mass parameters V0 and A0 are also indicated. Envelope expansion and contraction solutions again occur alternately in the present series. Also note that the velocity profiles near the centre show an undulating pattern, with each successive solution possessing one more stagnation point. Physically, these profiles represent radial sub-magnetosonic oscillations in a self-similar manner and the oscillatory pattern travels outwards. The key difference between the two sequences of MHD solutions obtained in this section is that the type B MHD solutions do not possess central BHs, while type A MHD solutions do. This also provides some insight for our proposed removal of central BHs in those type A solutions, due to the expected modifications of GP EoS when the mass density becomes exceedingly high. 6.3 Self-similar MHD shock solutions As shown in Section 6.1, asymptotic profiles at large x in the form of solutions (23) and (24) cannot always be realized for globally smooth MHD solutions by arbitrarily prescribed asymptotic parameters A0 and V0. In realistic astrophysical situations, however, MHD shock waves may be present and accordingly, our global MHD solution space can be substantially enriched by including MHD shocks (note that only fast MHD shocks are involved here). In this section, we construct global MHD solutions with shock waves and examine their properties, by combining our general idea of numerical MHD integration with the specific procedures for implementing shock conditions outlined in Section 4. Much experience has been gained in our earlier studies of constructing shock solutions and their applications. For example, in the context of protostar formation process, Lou & Gao (2006) analysed self-similar polytropic shocks propagating into an expanding or contracting outer envelope. To study the nature of intense magnetic fields in compact stellar remnants produced in aftermaths of supernovae, Lou & Wang (2007) examined a polytropic quasi-spherical rebound MHD shock model with central quasi-magnetostatic MHD configuration. Lou & Zhai (2009, 2010) explored isothermal self-similar void solutions with shocks and modelled X-ray luminosity profiles of dynamic PNe, with further extensions to the polytropic case made by Lou & Wang (2012). Despite our past experience, to enumerate each of the qualitatively distinct shock wave solutions would still be a laborious, if not impossible task, given the vastness of the number of possible combinations. In this study, instead of exhausting all types of MHD shock solutions, we focus mainly on solutions that contain only one shock discontinuity and have an outer MHD EWCS envelope. Taking on a more structural point of view, we aim at obtaining a better understanding of how various types of inner configurations come into existence and how they are linked to each other. We begin by inserting MHD shock waves at arbitrary positions xS on the EWCS9 and then integrate inmagnetostatic EWCS configurationwards from the other side of the shock. The typical central configuration thus obtained has a void, as shown in Fig. 16. The asymptotic behaviour of v(x) and α(x) near these voids satisfies solutions (48) to (51) with a ‘flat’ velocity profile and a finite positive limiting density, in contrast with the other kind of void that appears in globally smooth solutions, such as solutions e, f and g in Fig. 9. Another interesting fact to note in Fig. 16 is that the void radius x1 is not a monotonic function of the shock position xS. For instance, the shock position of solution c lies between those two shocks of solutions b and d, but the void radii of b and d are very close to one another and both are considerably larger than the void radius of solution c. This naturally leads to the question how exactly x1 varies with xS. Our investigation shows that this variation can be highly non-trivial. In particular, there are some special values of xS at which x1 = 0, i.e. the void disappears, giving way to alternative MHD central configurations. Figure 16. Open in new tabDownload slide Void solutions with an MHD shock and an outer EWCS static envelope. These four solutions are constructed with the pertinent parameters s = 1.2, γ = 1.7 and h = 5.0 and are obtained simply by inserting shock waves at four arbitrarily chosen positions xS on the EWCS and integrating inwards. All these four solutions have central voids and the solution behaviour near the void rim x1 is characterized by solutions (48) through (51). The xS and x1 values for each solution are listed in the lower panel. Note that x1 is not monotonic with respect to xS. Figure 16. Open in new tabDownload slide Void solutions with an MHD shock and an outer EWCS static envelope. These four solutions are constructed with the pertinent parameters s = 1.2, γ = 1.7 and h = 5.0 and are obtained simply by inserting shock waves at four arbitrarily chosen positions xS on the EWCS and integrating inwards. All these four solutions have central voids and the solution behaviour near the void rim x1 is characterized by solutions (48) through (51). The xS and x1 values for each solution are listed in the lower panel. Note that x1 is not monotonic with respect to xS. To better understand this issue, consider Fig. 17 that shows void size x1 as a function of xS obtained by direct numerical integration. The solid curve in Fig. 17 is drawn for the pertinent parameters s = 0.2929, γ = 1.2 and h = 0.1 (i.e. a type A regime parameter set), whereas the dashed curve is for s = 0.2, γ = 1.7 and h = 0.1 (i.e. a type B regime parameter set).10 Both of the curves show an ‘arch’-like pattern on the log-linear diagram, which is to be expected from the non-monotonicity. The non-trivial behaviour lies in between those ‘arches’, where x1 decreases to less than 10−4 (as judging purely from the curves themselves). Figure 17. Open in new tabDownload slide Void size x1 as a function of shock position xS, drawn for solutions with a shock and an outer EWCS envelope. The solid curve is for s = 0.2929, γ = 1.2 and h = 0.1, while the dashed curve is for s = 0.2, γ = 1.7 and h = 0.1. The two curves refer to the horizontal axes at the top and the bottom respectively. Both these curves show an ‘arch’-like pattern. The inset is a zoom-in view of the solid curve, in a tiny region between two ‘arches’, where xS is around 1.3645. As can be seen, the solid curve is ‘broken’ between the two ‘arches’. The ‘missing’ part corresponds to solutions that hit the MSCS when being integrated inwards. The hollow circle labelled by ‘BH’ corresponds to the solution that crosses the MSCC smoothly once and connects to the outer EWCS envelope via a shock, while the almost vertical part of the curve labelled by ‘LX-1’ represents the solution that has an LX-1 type behaviour near the centre x = 0 and is connected to the outer EWCS by a shock. These two kinds of special solutions are further presented in Fig. 18. Similar patterns exist between any two ‘arches’ of the solid curve. Meanwhile, the shape of the dashed curve between ‘arches’ is much simpler. It approaches 0 from either side and the corresponding special solution is one that has a central quasi-magnetostatic MHD configuration and an outer EWCS envelope, connected by an MHD shock. Figure 17. Open in new tabDownload slide Void size x1 as a function of shock position xS, drawn for solutions with a shock and an outer EWCS envelope. The solid curve is for s = 0.2929, γ = 1.2 and h = 0.1, while the dashed curve is for s = 0.2, γ = 1.7 and h = 0.1. The two curves refer to the horizontal axes at the top and the bottom respectively. Both these curves show an ‘arch’-like pattern. The inset is a zoom-in view of the solid curve, in a tiny region between two ‘arches’, where xS is around 1.3645. As can be seen, the solid curve is ‘broken’ between the two ‘arches’. The ‘missing’ part corresponds to solutions that hit the MSCS when being integrated inwards. The hollow circle labelled by ‘BH’ corresponds to the solution that crosses the MSCC smoothly once and connects to the outer EWCS envelope via a shock, while the almost vertical part of the curve labelled by ‘LX-1’ represents the solution that has an LX-1 type behaviour near the centre x = 0 and is connected to the outer EWCS by a shock. These two kinds of special solutions are further presented in Fig. 18. Similar patterns exist between any two ‘arches’ of the solid curve. Meanwhile, the shape of the dashed curve between ‘arches’ is much simpler. It approaches 0 from either side and the corresponding special solution is one that has a central quasi-magnetostatic MHD configuration and an outer EWCS envelope, connected by an MHD shock. At the top of Fig. 17, an inset is provided that zooms in one part of the solid curve between its two right-most ‘arches’ at 1.36442 ≤ xS ≤ 1.36448. With the help of this inset, the variation of MHD shock solutions with this range of xS can be described as follows: as xS increases from 1.36442, the void radius x1 gradually decreases, until xS = 1.364 5301 (denoted by a circle in the inset), when the solution meets the MSCS of X(x, α, v) = 0 at a saddle point on the MSCC along one of its two eigendirections. Therefore, this MHD solution continues to cross the MSCC smoothly, resulting in a central small BH configuration (with label ‘BH’). As xS increases further, the solution still meets the critical surface, but the meeting point no longer lies on the MSCC, so the desired global solution does not exist (unless we insert a second shock), resulting in a ‘gap’ in the curve. However, when xS increases to 1.364 5500, a global solution emerges again with a central LX-1 type profile (with label ‘LX-1’ in the inset) and x1 = 0. Finally, as xS increases further beyond this value, a central void with ever increasing radius re-emerges. This whole pattern in its qualitative form, as our investigation demonstrates, not only occurs in the range of xS indicated above, but also appears between all neighbouring ‘arches’ and is very likely to be a generic feature for type A regime parameters. In the above description, we point out the existence of two special kinds of MHD shock solutions with an EWCS envelope for s = 0.2929, γ = 1.2 and h = 0.1. The first kind has a central small BH and the second a central LX-1 type configuration. These two kinds of solutions are again paired, much like the pairing occurring in Section 6.2. Their existence has been further confirmed using the phase-diagram matching method, in light of which it is much easier to see that each of these two kinds of solutions again form an infinite sequence similar to those examined in the previous section.11 Fig. 18 shows one solution of the first kind (a΄) and three successive solutions of the second kind (a, b, c) for illustration purposes, with their detailed data contained in Table 1. No matter whether the shock is on the contracting inner region or the static outer envelope of the EWCS, the shock positions xS (see the second column of Table 1) correspond well to the positions of small x1 on the solid curve of Fig. 17. Moreover, the closer the distance in x between the MHD shock and the EWCS wavefront, the smaller the amplitude of velocity or density discontinuity in that shock. Figure 18. Open in new tabDownload slide Solid curves show the first three successive MHD solutions (in the order a, b, c) that have central LX-1-type behaviour and an MHD shock wave expanding into an outer MHD EWCS type configuration. The dashed curve is MHD solution a΄ that has a central small BH. MHD solutions a and a΄ are paired in the sense that they are almost the same type at reasonably large x. All of these MHD solutions are obtained by the phase-matching scheme with the three pertinent parameters s = 0.2929, γ = 1.2 and h = 0.1, and other relevant data are contained in Table 1. The phase-matching scheme shows that the limiting position of the MHD shock in this series of solutions coincides with the expanding wavefront of the EWCS. Also note that MHD shocks lying closer to the EWCS wavefront tend to have smaller extents of discontinuity in both velocity v(x) and density α(x). The MSCC is shown by the dash–dotted curve in both panels; v(x) versus x are shown in linear-log scales in the top panel and α(x) versus x are shown in log–log scales in the bottom panel. Figure 18. Open in new tabDownload slide Solid curves show the first three successive MHD solutions (in the order a, b, c) that have central LX-1-type behaviour and an MHD shock wave expanding into an outer MHD EWCS type configuration. The dashed curve is MHD solution a΄ that has a central small BH. MHD solutions a and a΄ are paired in the sense that they are almost the same type at reasonably large x. All of these MHD solutions are obtained by the phase-matching scheme with the three pertinent parameters s = 0.2929, γ = 1.2 and h = 0.1, and other relevant data are contained in Table 1. The phase-matching scheme shows that the limiting position of the MHD shock in this series of solutions coincides with the expanding wavefront of the EWCS. Also note that MHD shocks lying closer to the EWCS wavefront tend to have smaller extents of discontinuity in both velocity v(x) and density α(x). The MSCC is shown by the dash–dotted curve in both panels; v(x) versus x are shown in linear-log scales in the top panel and α(x) versus x are shown in log–log scales in the bottom panel. Table 1. Some estimated key parameters for the first three pairs of special MHD shock solutions with s = 0.2929, γ = 1.2 and h = 0.1. The MHD solutions a, b, c and a΄ are shown in Fig. 18, while MHD solutions b΄ and c΄ with central small BHs are paired with MHD solutions b and c, respectively, but are not shown in Fig. 18. The symbol xS is the position of an MHD shock, T is the leading coefficient in asymptotic solution (60), while x*, x0 and m0 denote the position of MSCC-crossing point, event horizon and reduced BH mass, respectively. Central LX-1 . xS . T . . . a 1.364 55 10.332 b 0.894 12 183.99 c 1.160 15 4259.6 Central BH xS x* x0 m0 a’ 1.364 53 5.7 × 10−4 4.02 × 10−5 1.37 × 10−4 b’ 0.894 11 2.8 × 10−9 2.06 × 10−12 7.04 × 10−12 c’ 1.160 15 7.4 × 10−15 3.21 × 10−20 1.10 × 10−19 Central LX-1 . xS . T . . . a 1.364 55 10.332 b 0.894 12 183.99 c 1.160 15 4259.6 Central BH xS x* x0 m0 a’ 1.364 53 5.7 × 10−4 4.02 × 10−5 1.37 × 10−4 b’ 0.894 11 2.8 × 10−9 2.06 × 10−12 7.04 × 10−12 c’ 1.160 15 7.4 × 10−15 3.21 × 10−20 1.10 × 10−19 Open in new tab Table 1. Some estimated key parameters for the first three pairs of special MHD shock solutions with s = 0.2929, γ = 1.2 and h = 0.1. The MHD solutions a, b, c and a΄ are shown in Fig. 18, while MHD solutions b΄ and c΄ with central small BHs are paired with MHD solutions b and c, respectively, but are not shown in Fig. 18. The symbol xS is the position of an MHD shock, T is the leading coefficient in asymptotic solution (60), while x*, x0 and m0 denote the position of MSCC-crossing point, event horizon and reduced BH mass, respectively. Central LX-1 . xS . T . . . a 1.364 55 10.332 b 0.894 12 183.99 c 1.160 15 4259.6 Central BH xS x* x0 m0 a’ 1.364 53 5.7 × 10−4 4.02 × 10−5 1.37 × 10−4 b’ 0.894 11 2.8 × 10−9 2.06 × 10−12 7.04 × 10−12 c’ 1.160 15 7.4 × 10−15 3.21 × 10−20 1.10 × 10−19 Central LX-1 . xS . T . . . a 1.364 55 10.332 b 0.894 12 183.99 c 1.160 15 4259.6 Central BH xS x* x0 m0 a’ 1.364 53 5.7 × 10−4 4.02 × 10−5 1.37 × 10−4 b’ 0.894 11 2.8 × 10−9 2.06 × 10−12 7.04 × 10−12 c’ 1.160 15 7.4 × 10−15 3.21 × 10−20 1.10 × 10−19 Open in new tab The functional relation between x1 and xS shown by the dashed curve of Fig. 17 turns out to be much simpler. Between two neighbouring ‘arches’, the value of x1 simply drops to zero and increases again. The solutions never have the problem of running into the MSCS, hence there are no ‘gaps’ in the curve. The special solutions corresponding to x1 = 0 have a quasi-magnetostatic configuration at the centre and they also form an infinite sequence. The first four solutions in the sequence are displayed in Fig. 19. Our numerical investigations suggest that this pattern is general for the type B regime parameters of γ, s and h. Figure 19. Open in new tabDownload slide The first four solutions that have an MHD shock wave expanding into an outer magnetostatic EWCS configuration are drawn for type B regime parameters s = 0.2, γ = 1.2 and h = 0.1. These MHD shock solutions are close counterparts of those MHD shock solutions shown in Fig. 18, obtained by matching central ‘quasi-magnetostatic’ solutions and the EWCS with MHD shocks inserted at various positions. Despite the disappearance of central BHs, the trend of the positions and strengths of MHD shocks are much the same as those in Fig. 18. Here, U is a key parameter of the ‘quasi-magnetostatic’ solution, as explained at the beginning of Section 3.5. Figure 19. Open in new tabDownload slide The first four solutions that have an MHD shock wave expanding into an outer magnetostatic EWCS configuration are drawn for type B regime parameters s = 0.2, γ = 1.2 and h = 0.1. These MHD shock solutions are close counterparts of those MHD shock solutions shown in Fig. 18, obtained by matching central ‘quasi-magnetostatic’ solutions and the EWCS with MHD shocks inserted at various positions. Despite the disappearance of central BHs, the trend of the positions and strengths of MHD shocks are much the same as those in Fig. 18. Here, U is a key parameter of the ‘quasi-magnetostatic’ solution, as explained at the beginning of Section 3.5. 7 SUMMARY AND DISCUSSION We have formulated here a quasi-spherically symmetric GP MHD model problem with random transverse magnetic fields. The magnetic field manifests on large scales as the net sum of a magnetic pressure plus a magnetic tension force. Paczynski–Wiita gravity (1) is adopted as an approximation to capture GR effects for the appearance of a central Schwarzschild BH, while avoiding the pitfalls of Newtonian gravity at high central mass density (Lian & Lou 2014). With MHD self-similar transformation (7), the PDEs describing our MHD model can be readily reduced to two coupled non-linear ODEs in terms of the reduced radial velocity v(x) and the reduced mass density α(x), along with an integral relation (11) for the magnetic flux ‘frozen-in’ condition. Two important dimensionless parameters s and h arise in our model MHD solutions. As in Lian & Lou (2014), we can readily relate s and h to pertinent physical quantities at small t and/or large r (i.e. initial and far away conditions). In particular, denoting the (effective) sound speed12 and the Keplerian speed by us and uk, respectively, we have |$sA_0^{3\gamma -3}=2u_{\rm s}^2/c^2$|⁠, |$sA_0 = 2u_{\rm k}^2/c^2$| and |$h = [\langle B_{\rm t}^2\rangle /(4\pi p)](u_{\rm s}/u_{\rm k})^2$|⁠, where A0 is the large-x mass parameter in asymptotic solution (24). For pertinent astrophysical applications, here we briefly estimate s and h value ranges for four different settings of astrophysical interests. (1) For forming a protostar in an H2 MC environment or a massive star of first generation in the early Universe, the initial temperature is estimated to be ∼101 K and the sound speed us is ∼0.1 km s−1. The Keplerian speed uk will be on that same order as estimated based on the Jeans mass and length-scales. Therefore, A0 ∼ 1 and s ∼ 10−12. For an MC number density range of 103 ∼ 105 cm−3 and an interstellar magnetic field of ∼10−6G, magnetic parameter h would be around 10−3–10−1. (2) For a post-AGB star or a young PN, the gas envelope may have a temperature around ∼103 K, a sound speed of us ≲ 101 km s− 1 and a Keplerian speed of uk ∼ 102 km s−1. If we take γ = 5/3, then A0 ∼ 0.01 and s ∼ 10−4. With a surface magnetic field strengths of 10−1 ∼ 101 G and a mass density of ∼10−10 g cm−3, magnetic parameter h would be of the order of 10−7 ∼ 10−3. (3) Near an SMBH and/or HMBH (or even a DMBH and an MMBH) at the centre of a galactic bulge, an elliptical galaxy, an ULIRG, an ELIRG, an X-ray/gamma-ray source, a quasar environment and an UDG in the Universe including the early Universe, the (effective) temperature may reach 107 ∼ 108 K or higher and both Keplerian speed uk and (effective) sound speed us would be around ∼102 ∼ 103 km s−1. In these situations, s = 10−7 ∼ 10−5 and A0 ∼ 1. With 1 ∼ 10 mG magnetic field strengths and a mean number density of ∼10 cm−3, magnetic parameter h may fall in the range of 1 ∼ 100. (4) In the collapsing core of a progenitor star prior to a violent Type II supernova explosion, the stellar core temperature can become as high as 1011 ∼ 1012 K or even higher. The gas is extremely relativistic and both us and uk may come close to the speed of light c. From these, we estimate s to be around 0.1 ∼ 1 and A0 ∼ 1. A proto-neutron star or a magnetar can have a typical magnetic field strength of ∼1012–1014 G and a mass density of ∼107 g cm−3, resulting in an h ≲ 1. These sample parameter estimates give a wide range of values our model parameters may take. For more specific model analysis, we mainly pick s and h around 0.1 ∼ 1 in the numerical examples, which happen to be mostly close to the last scenario described above. Nonetheless, the asymptotic MHD solutions and the qualitative behaviours shown by our numerical solutions are fairly general and may be readily adapted for specific combinations of these dimensionless parameters. The class of large-x asymptotic expansion solution of our MHD flows involves an essentially force-free magnetic field and it also includes the global magnetostatic solution as a special case. The central behaviours of this MHD model problem can be broadly classified into three categories: central Schwarzschild BH with growing mass, central void in expansion and allowed central asymptotic MHD solutions regular towards x → 0. First, around a central BH, the gas flow velocity tends to −∞ towards the growing event horizon and we emphasize that this infinity should be properly interpreted, by employing empirical conversion relation (40) that relates the PW model velocity u to the true physical velocity utrue measured by a local (accelerating) observer at rest with respect to the BH (see footnote 3). Clearly, this rescaling only gives approximate results and is not accurate from a strictly GR point of view. Nonetheless, this approximation is found to be a fairly good one in terms of Keplerian velocity profiles (within 5 per cent uncertainly) outside twice the Schwarzschild radius (see Abramowicz et al. 1996 for comparison between PW and GR results and empirical fits by a Lorentz factor). The central Schwarzschild BH in our model accretes mass, but no electromagnetic energy Poynting flux is directly accreted across the event horizon into the central BH. Secondly, for central expanding voids, there are two different types of dynamic behaviours depending on the relative speed of void expansion and gas outflow. If the former is slower than the latter, then the gas density vanishes over the void rim. In contrast, if the former is faster than the latter, then the gas density remains finite over the void rim. Although the detailed density profile depends on polytropic index γ in a somewhat subtle manner, it can be seen some of these behaviours would not occur without random transverse magnetic fields. Thirdly, for solutions that extend to the centre x → 0, we recovered the quasi-magnetostatic solution of Wang & Lou (2008) that results from perturbing one of the two global static solutions for γ > 4/3 at the centre and the generalized LP solution of Lou & Shi (2014b) for γ > 4/3 without magnetic field. More importantly, we discovered a brand new class of MHD asymptotic solutions for γ < 4/3, referred to as LX-1 and LX-2, that plays an indispensable role in constructing all possible numerical MHD solutions. In reference to the procedure of derivation, we expect these novel solutions to also exist even for a GP hydrodynamics under Newtonian gravity in the absence of magnetic field, because the dimensionless parameter s does not explicitly enter the lowest order expansion for such novel asymptotic solutions. In contexts of dynamic formation of BHs in host (DM) mass reservoirs in a wide mass range, we would make clear that magnetic field directly interacts with (partially) ionized plasma of a certain mass fraction. Such magnetized plasma gas medium can be dynamically trapped in a PW gravitational potential of host (DM) mass reservoirs in dynamic collapses. At the centres of such dynamic collapsing mass reservoirs, SMBHs, HMBHs, DMBHs and MMBHs in a wide mass range may emerge. Besides dynamic effects of magnetic field, the collateral effects associated with the magnetic field may also serve as valuable diagnostic tools. The MSCC of our MHD model problem (if present) appears in a loop pattern as projected on to the −v versus x plane. This loop geometry is determined by the number of roots of non-linear algebraic equation (21) and has a major influence on the global MHD solution structure, such as the number of EWCSs and the range of x → +∞ asymptotic parameters leading to various distinct MHD solution behaviours at finite x. The latter aspect of the MHD solution space structure was extensively explored, where we particularly emphasize that the one-dimensional MSCC has a two-dimensional region of attraction on the plane of asymptotic parameters as shown in Figs 7 and 10. Also, a higher central mass accretion rate concurs with a stronger random transverse magnetic field. MHD solutions can cross the MSCS X(x, α, v) = 0 via either MHD shock waves or the MSCC. Our numerical exploration of these MHD solutions show a remarkable feature for γ < 4/3, namely the existence of solution pairs that are almost identical at reasonably large x (i.e. for small t and/or large r) but differ completely at very small x (i.e. for large t and/or small r). Specifically, one solution in this solution pair has a small central BH, while the other solution satisfies the LX-1 asymptotic MHD solution at the centre. When we use the phase-matching technique (for MSCC-crossing smooth solutions) or vary the MHD shock position continuously (for MHD shock solutions), these global MHD solution pairs naturally arrange themselves into an infinite discrete sequence. In realistic astrophysical scenarios, there are conceivable situations where initially similar systems can evolve into completely different end states. For example, supernova explosions of massive progenitor stars are known to produce neutron stars and are also suspected to possibly give rise to stellar mass BHs or simply voids under various conditions. It is therefore particularly revealing that our simple MHD model can also produce several different results under quite similar initial situations and we would pursue astrophysical applications of these MHD solution pairs in a separate paper. We also hope that our new MHD solutions could trigger pertinent numerical simulation studies to determine whether both of these MHD solutions could actually occur in nature by setting some reasonable initial/boundary conditions. It should be made clear that our investigation of solutions with MHD shocks is far from complete, while we have set the stage and paved the way for various further applications. In fact, we have arranged the outer profile to be part of an EWCS in all our MHD shock solutions. To some degree, this may be considered as the opposite extreme of what we have done in the investigation of MSCC-crossing solutions, where the magnitude of the ‘shock’ is, in effect, fixed at zero and no restriction is imposed on the outer envelope profile. We leave it for future works to consider other more general situations. Besides a more thorough exploration of the MHD solution space, we envisage our work to be extended in various other directions as well. For instance, stability analysis will help to determine whether the self-similar MHD solutions can actually occur and how long they will last and well-designed numerical MHD simulation can provide valuable missing links between the initial/boundary condition and the final state (whether it tends to a self-similar flow and which one). One may also explore other magnetic field configurations or other approximation schemes for self-gravity in reference to GR. It is also highly desirable to pursue a complete GR treatment of such a GP MHD model problem with combined analytical and numerical analyses. Acknowledgments This research was supported in parts by the Ministry of Science and Technology (MOST) under the State Key Development Programme for Basic Research grant 2012CB821800, by the National Natural Science Foundation of China (NSFC) grants 10373009, 10533020, 11073014 and 11473018 at Tsinghua University, by MOST grant 2014YQ030975, by Tsinghua University Initiative Scientific Research Programme (20111081008), by Tsinghua Centre for Astrophysics (THCA) and by the Yangtze Endowment, the Special Research Fund for Doctoral Programme (SRFDP) 20050003088, 200800030071 and 20110002110008 as well as 985 grants from the Ministry of Education (MoE) at Tsinghua University. YQL acknowledges the support of the China–Chile Scholarly Exchange Program administered by the Chinese Academy of Sciences South America Center for Astronomy (CASSACA). 1 " For γ = 1, the third term on the LHS of integral (34) becomes ln α. 2 " This argument fails for γ = 1, as do asymptotic solutions (35) and (36), which, for γ = 1, should be changed to |$v(x)\cong -\sqrt{e(x)} -\ln e(x)/[2\sqrt{e(x)}] + \ldots$| and α(x)≅ − 1/[s2m0v(x)], where e(x) = 2m0/(x − sm0). This subtlety was overlooked in the GP hydrodynamic model analysis of Lian & Lou (2014). 3 " In this regard, both Abramowicz (2009) and Lian & Lou (2014) contain a typo by missing the square root in the denominator of conversion formula (40). 4 " Another possible parameter range would be 0 < γ < 2/3 for the RHS of inequality(58) to be positive. In certain polytropic models involving turbulent medium, the range of polytropic index γ may include 0 < γ < 1 as mentioned at the end of Section 2.1. 5 " If one would consider γ ≤ 1, then the first subclass can be generalized down to γ > 2/3, but the second subclass only exists for 1 < γ < 4/3 and cannot be generalized further. 6 " The most common form of Bernoulli's inequality states that (1 + x)γ ≥ 1 + γx for γ > 1 and x > −1. Here, we just take x to be −λ. 7 " In numerical integrations, it is not always possible to have the solution exactly hit the MSCC, as inherent numerical errors cannot be absolutely avoided. Presumably, however, this region gives us some ideas about the range of A0 and V0 that will lead to an MSCC event. 8 " Should we plot b΄, the second solution of the second sequence, it would lie so close to solution b such that the two can no longer be distinguished. 9 " The front and back ends of an MHD shock will, in general, have different x values. To avoid confusion, xS always denotes the front end (the upstream side) of the shock from now on, unless otherwise stated. 10 " Recall that there are two EWCSs for the latter set of parameters and we have only taken the one with dense envelope as the outer configuration in preparing Fig. 17. A similar investigation can be carried out for the dilute EWCS, but our numerical exploration did not find any solution other than those with a central void in that case. Therefore, we leave it out from our main discussion. 11 " As a logical consequence of this, the solid curve in Fig. 17 should be interpreted as to possess an infinite number of ‘arches’. Most of the arches would be packed between xS = 1.12 and 1.17, where the figure fails to present due to the smallness of x1. 12 " In proper astrophysical contexts, this us could represent random velocity dispersions of (DM) particles or even stars. REFERENCES Abramowicz M. A. , 2009 , A&A , 500 , 213 Crossref Search ADS Abramowicz M. A. , Beloborodov A. M., Chen X.-M., Igumenshchev I. V., 1996 , A&A , 313 , 334 Ackermann M. et al. , 2013 , Science , 339 , 807 Crossref Search ADS PubMed Barenblatt G. , 1980 , J. Appl. Math. Mech. , 44 , 267 Crossref Search ADS Berezhko E. G. , Elshin V. K., Krymskii G. F., Petukhov S. I., 1988 , Sov. Phys.-Usp. , 31 , 27 Crossref Search ADS Bertschinger E. , 1985a , ApJS , 58 , 1 Crossref Search ADS Bertschinger E. , 1985b , ApJS , 58 , 39 Crossref Search ADS Blandford R. , Eichler D., 1987 , Phys. Rep. , 154 , 1 Crossref Search ADS Blandford R. D. , Znajek R. L., 1977 , MNRAS , 179 , 433 Crossref Search ADS Braithwaite J. , Spruit H. C., 2004 , Nature , 431 , 819 Crossref Search ADS PubMed Buckingham E. , 1914 , Phys. Rev. , 4 , 345 Crossref Search ADS Cai M. J. , Shu F. H., 2005 , ApJ , 618 , 438 Crossref Search ADS Chandrasekhar S. , 1961 , Hydrodynamic and Hydromagnetic Stability Dover Publications, Inc. , New York Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Chevalier R. A. , 1982 , ApJ , 258 , 790 Crossref Search ADS Chevalier R. A. , Blondin J. M., Emmering R. T., 1992 , ApJ , 392 , 118 Crossref Search ADS Chiueh T. , Chou J.-K., 1994 , ApJ , 431 , 380 Crossref Search ADS Eatough R. P. et al. , 2013 , Nature , 501 , 391 Crossref Search ADS PubMed Ferrario L. , Wickramasinghe D., 2006 , MNRAS , 367 , 1323 Crossref Search ADS Fillmore J. A. , Goldreich P., 1984a , ApJ , 281 , 1 Crossref Search ADS Fillmore J. A. , Goldreich P., 1984b , ApJ , 281 , 9 Crossref Search ADS Foster P. N. , Chevalier R. A., 1993 , ApJ , 416 , 303 Crossref Search ADS Fu T.-M. , Gao Y., Lou Y.-Q., 2011 , MNRAS , 741 , 113 Gao Y. , Lou Y.-Q., 2017 , MNRAS , 466 , L53 Crossref Search ADS Gehman C. S. , Adams F. C., Watkins R., 1996 , ApJ , 472 , 673 Crossref Search ADS Gibbings J. C. , 2011 , Dimensional Analysis Springer Science & Business Media , London, Dordrecht, Heidelberg, New York Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Goldreich P. , Weber S. V., 1980 , ApJ , 238 , 991 Crossref Search ADS Govoni F. , Feretti L., 2004 , Int. J. Mod. Phys. D , 13 , 1549 Crossref Search ADS Harada T. , Maeda H., 2001 , Phys. Rev. D , 63 , 084022 Crossref Search ADS Harada T. , Maeda H., Semelin B., 2003 , Phys. Rev. D , 67 , 084003 Crossref Search ADS Hartman P. , 1982 , Ordinary Differential Equations , 2nd edn Birkhäuser , Boston Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Hawley J. F. , Balbus S. A., 1991 , ApJ , 376 , 223 Crossref Search ADS Hu R.-Y. , Lou Y.-Q., 2008 , MNRAS , 390 , 1619 Hu R.-Y. , Lou Y.-Q., 2009 , MNRAS , 396 , 878 Crossref Search ADS Hunter C. , 1977 , ApJ , 218 , 834 Crossref Search ADS Jaffe W. , 1980 , ApJ , 241 , 925 Crossref Search ADS Larson R. B. , 1969 , MNRAS , 145 , 271 Crossref Search ADS Lian B. , Lou Y.-Q., 2014 , MNRAS , 438 , 1242 Crossref Search ADS Lizano S. , Shu F. H., 1989 , ApJ , 342 , 834 Crossref Search ADS Lou Y.-Q. , 1994 , ApJ , 428 , L21 Crossref Search ADS Lou Y.-Q. , 1996 , MNRAS , 279 , 67 Crossref Search ADS Lou Y.-Q. , 2015a , Proc. SALT Science Conf. Stellenbosch Institute of Advanced Study, South Africa Lou Y.-Q. , 2015b , MNRAS , 454 , 2815 Crossref Search ADS Lou Y.-Q. , Cao Y., 2008 , MNRAS , 384 , 611 Crossref Search ADS Lou Y.-Q. , Gao Y., 2006 , MNRAS , 373 , 1610 Crossref Search ADS Lou Y.-Q. , Gao Y., 2011 , MNRAS , 412 , 1755 Crossref Search ADS Lou Y.-Q. , Hu X.-Y., 2016 , MNRAS , 459 , 2014 Lou Y.-Q. , Jiang Y.-F., 2008 , MNRAS , 391 , L44 Crossref Search ADS Lou Y.-Q. , Lian B., 2012 , MNRAS , 420 , 2147 Crossref Search ADS Lou Y.-Q. , Shen Y., 2004 , MNRAS , 348 , 717 Crossref Search ADS Lou Y.-Q. , Shi C.-H., 2014a , MNRAS , 442 , 741 Crossref Search ADS Lou Y.-Q. , Shi C.-H., 2014b , MNRAS , 445 , 1186 Crossref Search ADS Lou Y.-Q. , Wang L., 2012 , MNRAS , 420 , 1897 Crossref Search ADS Lou Y.-Q. , Wang W.-G., 2006 , MNRAS , 372 , 885 Crossref Search ADS Lou Y.-Q. , Wang W.-G., 2007 , MNRAS , 378 , L54 Crossref Search ADS Lou Y.-Q. , Wu Y., 2005 , MNRAS , 364 , 475 Crossref Search ADS Lou Y.-Q. , Wu Y., 2012 , MNRAS , 422 , L28 Crossref Search ADS Lou Y.-Q. , Zhai X., 2009 , Ap&SS , 323 , 17 Crossref Search ADS Lou Y.-Q. , Zhai X., 2010 , MNRAS , 408 , 436 Crossref Search ADS Lou Y.-Q. , Jiang Y.-F., Jin C.-C., 2008 , MNRAS , 386 , 835 Crossref Search ADS Low B. C. , Lou Y. Q., 1990 , ApJ , 352 , 343 Crossref Search ADS Maloney P. , 1988 , ApJ , 334 , 761 Crossref Search ADS Martí-Vidal I. , Muller S., Vlemmings W., Horellou C., Aalto S., 2015 , Science , 348 , 311 Crossref Search ADS PubMed Masuyama M. , Shigeyama T., Tsuboki Y., 2016 , PASJ , 68 , 22 McKee C. F. , Holliman J. H. II, 1999 , ApJ , 522 , 313 Crossref Search ADS McKee C. F. , Zweibel E. G., 1995 , ApJ , 440 , 686 Crossref Search ADS McLaughlin D. E. , Pudritz R. E., 1996 , ApJ , 469 , 194 Crossref Search ADS Myers P. C. , 1998 , ApJ , 496 , L109 Crossref Search ADS Oesch P. A. et al. , 2016 , ApJ , 819 , 129 Crossref Search ADS Paczyński B. , Wiita P. J., 1980 , A&A , 88 , 23 Parker E. N. , 1979 , Cosmical Magnetic Fields: Their Origin and Their Activity Oxford Univ. Press , Oxford Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Penston M. V. , 1969 , MNRAS , 144 , 425 Crossref Search ADS Reindl N. , Rauch T., Parthasarathy M., Werner K., Kruk J. W., Hamann W.-R., Sander A., Todt H., 2014 , A&A , 565 , A40 Crossref Search ADS Sedov L. I. , 1959 , Similarity and Dimensional Methods in Mechanics Academic Press , New York Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Shi C. H. , Lou Y.-Q., 2017 , MNRAS , preprint (arXiv:e-print) Shu F. H. , 1977 , ApJ , 214 , 488 Crossref Search ADS Shu F. H. , Milione V., Gebel W., Yuan C., Goldsmith D. W., Roberts W. W., 1972 , ApJ , 173 , 557 Crossref Search ADS Shu F. H. , Adams F. C., Lizano S., 1987 , ARA&A , 25 , 23 Crossref Search ADS Shu F. H. , Lizano S., Galli D., Cantó J., Laughlin G., 2002 , ApJ , 580 , 969 Crossref Search ADS Spaans M. , Silk J., 2000 , ApJ , 538 , 115 Crossref Search ADS Suto Y. , Silk J., 1988 , ApJ , 326 , 527 Crossref Search ADS Taylor G. , 1950 , Proc. R. Soc. A , 201 , 175 Crossref Search ADS Thompson C. , Duncan R. C., 1993 , ApJ , 408 , 194 Crossref Search ADS Toptygin I. N. , 2000 , Astron. Lett. , 26 , 356 Crossref Search ADS Tsai J. C. , Hsu J. J. L., 1995 , ApJ , 448 , 774 Crossref Search ADS van Dokkum P. et al. , 2016 , ApJ , 828 , L6 Crossref Search ADS Velikhov E. P. , 1959 , J. Exp. Theor. Phys. , 36 , 1398 Velikovich A. L. , Liberman M. A., Shmalts R. F., 1985 , Sov. J. Exp. Theor. Phys. Lett. , 41 , 97 Viala Y. , Horedt G. P., 1974a , A&AS , 16 , 173 Viala Y. P. , Horedt G., 1974b , A&A , 33 , 195 Wang L. , Lou Y.-Q., 2014 , MNRAS , 439 , 2323 Crossref Search ADS Wang F. et al. , 2015 , ApJ , 807 , L9 Crossref Search ADS Wang W.-G. , Lou Y.-Q., 2007 , Ap&SS , 311 , 363 Crossref Search ADS Wang W.-G. , Lou Y.-Q., 2008 , Ap&SS , 315 , 135 Crossref Search ADS Weitering H. , 2016 , Scientific American , SPACE.com Whitworth A. , Summers D., 1985 , MNRAS , 214 , 1 Crossref Search ADS Wu Y. , Lou Y.-Q., 2006 , MNRAS , 372 , 992 Crossref Search ADS Wu X.-B. et al. , 2015 , Nature , 518 , 512 Crossref Search ADS PubMed Xiang-Grüss M. , Lou Y.-Q., Duschl W., 2009a , MNRAS , 397 , 815 Crossref Search ADS Xiang-Grüss M. , Lou Y.-Q., Duschl W., 2009b , MNRAS , 400 , L52 Crossref Search ADS Yahil A. , 1983 , ApJ , 265 , 1047 Crossref Search ADS Yu C. , Lou Y.-Q., 2005 , MNRAS , 364 , 1168 Crossref Search ADS Yu C. , Lou Y.-Q., Bian F.-Y., Wu Y., 2006 , MNRAS , 370 , 121 Crossref Search ADS Zeldovich Y. B. , Novikov I. D., 1971 , Relativistic Astrophysics. Vol. 1: Stars and Relativity Univ. Chicago Press , Chicago Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Zweibel E. G. , Heiles C., 1997 , Nature , 385 , 131 Crossref Search ADS APPENDIX A: ASYMPTOTIC BEHAVIOUR OF MSCC NEAR THE CENTRE In this appendix, we perform a detailed analysis on the asymptotic behaviour of MSCC near x = 0. The results of the analysis are valuable for constructing numerical solutions, because we, at times, need positions on the MSCC with extremely small x, as the starting points for numerical integration (see Section 6.2). Recall that the MSCC loop (excluding the ZML) touches x = 0 only when 1 < γ < 4/3, and its projection on to the −v versus x (i.e. α = 0) plane displays a ‘kink’-like feature near the touching point. We have already indicated that the tangent at one side of the ‘kink’ is v = x/γ. Assuming α = ΛxN for x ≪ 1 and substituting it into non-linear algebraic equation (74), we immediately determine N and Λ as \begin{equation} N = -2+2/[3(\gamma -1)] > 0, \end{equation} (A1) and \begin{equation} \Lambda =\left[ \frac{1}{\gamma }\left( 1-\frac{1}{\gamma } \right)^{4-2\gamma }\right]^{1/(3\gamma -3)}\ . \end{equation} (A2) On the other side (i.e. the v < 0 side) of the kink, the MSCC is singular in x and it would be better to regard v as an independent variable. We then assume x = Σ(−v)n and α = Θ(−v)−m for v ≪ 1, where the two coefficients, Σ and Θ, and the two exponents n and m are required to satisfy \begin{equation} \Sigma >0\ ,\qquad \Theta >0\ ,\qquad n>1\ ,\qquad m>0\ . \end{equation} (A3) Balance of the dominant terms in two algebraic (74) and (75) leads to n = 2/(γ − 1) − 5 and m = n − 1, and the vanishing of their coefficients yields Σ = 1/[8γ1/(γ − 1)] and Θ = 2/Σ. One can readily verify that all inequalities in (A3) are satisfied. Switching back to x as the independent variable, on the v < 0 side, we have equivalently \begin{equation} -v = \left( \frac{x}{\Sigma } \right)^{1/n}\ ,\qquad \alpha = \frac{2}{\Sigma }\left( \frac{x}{\Sigma } \right)^{-(n-1)/n}\ . \end{equation} (A4) To better appreciate the three-dimensional features of the MSCC, the interested reader may consult the dash–dotted curves in Fig. 9 and Fig. 12 and relate the MSCC behaviour near x = 0 to expressions (A4) derived above. APPENDIX B: ON THE BIFURCATION PROPERTY OF THE POINT WITH v = 0 ALONG MSCC This appendix is to determine the point with v = 0 on the MSCC as a saddle-node bifurcation point. We are not interested in the origin x = v = 0 and only consider points with v = 0 and x > 0. Explicit calculations using expressions (81) through (86) lead to the following expression for the product λ1λ2 as \begin{eqnarray} \lambda _1\lambda _2 &=&{ } \frac{2v[(-3\gamma ^2+8\gamma -4)x-\gamma ^2v]}{\gamma x^2 (x-v)^2}\mathcal {M}^2 \nonumber\\ &&+ \frac{2v[(2\gamma ^2-11\gamma +6)x + (6\gamma ^2 -3\gamma )v]}{\gamma x^2}\mathcal {M} \nonumber\\ &&+ \left[\frac{\gamma -2}{x-v} \mathcal {M} + 3(x-v)\right]\frac{2mv}{(x-sm)^3}, \end{eqnarray} (B1) where |$\mathcal {M} \equiv (x-v)^2-hx^2\alpha$|⁠. Given that there is an overall factor of v on the RHS of equation (B1), it is evident that λ1λ2 must change its sign at v = 0 under generic situations, i.e. when limv → 0λ1λ2/v ≠ 0, which is almost always the case. Therefore, we see that v = 0 is a bifurcation point and we always have saddle points on one side of the bifurcation. Meanwhile, the expression for Δ at v = 0 reads \begin{equation} \Delta |_{v=0} = 2x^2 + 4hx\alpha (\gamma -2)\ . \end{equation} (B2) Recall that at the point on MSCC with v = 0, x and α satisfy relation (76), leading to the conclusion that expression (B2) and \begin{equation} \mathcal {P}(A) \equiv (\gamma A^{3\gamma -3} + hA)^{3/2} - hA(2-\gamma ) \end{equation} (B3) bear the same sign, where A > 0 is the corresponding root of equation (21) that lies within (0, 1/s). Basic calculus reveals that |$\mathcal {P}(0) = 0$| and |$\mathcal {P}(A)$| is monotonically non-decreasing in the open interval (0, 1/s). Therefore, Δ|v = 0 always stays positive and the bifurcation is always between nodes and saddle points. APPENDIX C: A MATHEMATICAL PROOF OF THE IMPOSSIBILITY OF CROSSING THE PACZYNSKI–WIITA SINGULARITY The Paczynski–Wiita gravity is physically valid only when r > 2GM/c2 or x − sm > 0 after MHD self-similar transformation (7). Should this physical restriction be relaxed, then x − sm = 0 would become a singularity in our model. However, this singularity would not cause essential difficulty since any solution, as long as it satisfies x − sm > 0 initially, will always stay x − sm > 0 and here is a mathematical proof. We first introduce the set Ω defined by \begin{equation} \Omega = \lbrace (x,\ \alpha ,\ v)|x>0, \alpha > 0, x-v>0\rbrace \end{equation} (C1) as a region of the three-dimensional Euclidean space spanned by x, α and v. In this region Ω, we can define a vector field |$\boldsymbol {V}(x, \alpha , v) = (x-sm)^2(X\ , Y\ , Z)$|⁠, where X, Y and Z are defined by expressions (17)–(19), i.e., they are the same three functionals that appear in our coupled non-linear MHD ODEs (15) and (16). Therefore, the integral curves of |$\boldsymbol {V}$|⁠, which we denote by [x(τ) , α(τ), v(τ)] where τ is a real parameter, satisfy \begin{equation} \frac{{\rm d}}{{\rm d}\tau }(x,\alpha ,v)=\boldsymbol {V}(x\ ,\alpha \ ,v)\ . \end{equation} (C2) Clearly, this equation is equivalent to ODEs (15) and (16) when (x − sm)2X ≠ 0. Therefore, solutions to ODEs (15) and (16) are integral curves of |$\boldsymbol {V}$|⁠. Consider now the surface Σ defined by x − sm = 0, i.e. α = 1/[sx(x − v)]. It separates Ω into two regions, one with x − sm > 0 and the other with x − sm < 0. Equivalently, we define the surface Σ as |${\cal F}(x,\ \alpha ,\ v)\equiv x(x-v)\alpha -1/s=0$|⁠. The normal at a point (x, α, v) on this Σ surface is then given by the gradient operation |$\nabla {\cal F}(x,\ \alpha ,\ v)=[\alpha (2x-v),\ x(x-v),\ -x\alpha ]$|⁠. We aim at proving that no integral curve of |$\boldsymbol {V}$| can pass from one region to the other and the key is to notice that |$\boldsymbol {V}$| becomes \begin{equation} \boldsymbol {V}|_{\Sigma } = (x-v)x^2\alpha [0,\ \alpha ,\ (x-v)] \end{equation} (C3) on the Σ surface |${\cal F}(x,\ \alpha ,\ v) \equiv x(x-v)\alpha -1/s=0$|⁠, which is tangent to Σ surface itself as the scalar product of |$\boldsymbol {V}|_{\Sigma }\cdot \nabla {\cal F}$| vanishes. Therefore, Σ consists of the solution integral curves of |$\boldsymbol {V}|_{\Sigma }$|⁠. Also it is not hard to see that every component of |$\boldsymbol {V}$| has continuous derivative near Σ, so the Picard–Lindelöf theorem of ODEs (Hartman 1982) can be applied to conclude that through any point on Σ, there is one unique integral curve of |$\boldsymbol {V}$|⁠. If there had been one integral curve of |$\boldsymbol {V}$| that passes from x − sm > 0 to x − sm < 0, it will cross Σ at some point |$\mathcal {Q}$| that already belongs to an integral curve that lies entirely on Σ. This contradicts the uniqueness and our proof is complete. © 2016 The Authors Published by Oxford University Press on behalf of the Royal Astronomical Society
TI - Relativistic self-similar dynamic gravitational collapses of a quasi-spherical general polytropic magnetofluid
JO - Monthly Notices of the Royal Astronomical Society
DO - 10.1093/mnras/stw3119
DA - 2017-05-11
UR - https://www.deepdyve.com/lp/oxford-university-press/relativistic-self-similar-dynamic-gravitational-collapses-of-a-quasi-8xRT5eqkSm
SP - 2
VL - 467
IS - 1
DP - DeepDyve
ER -