A&A 491, 671-680 (2008)
DOI: 10.1051/0004-6361:200809848

A three-dimensional magnetohydrodynamic model of planetary nebula jets, knots, and filaments

K. H. Tsui

Instituto de Física - Universidade Federal Fluminense, Campus da Praia Vermelha, Av. General Milton Tavares de Souza s/n Gragoatá, 24.210-346, Niterói, Rio de Janeiro, Brasil

Received 26 March 2008 / Accepted 4 September 2008

Abstract
Context. Aims. A self-organizational view of planetary nebulae, driven by global conservation properties of the magnetohydrodynamic (MHD) plasma, is presented. Self-similarity with a self-consistent evolution function is used, as a method and as a model for self-organized states, to solve the time-dependent MHD equations in the gravitational field of a central star.
Methods. The self-similar configurations are constructed on the basis of a limited radially expanding MHD plasma sphere, which could be associated with the ejected hydrogen layer of the AGB star due to its strong pulsations.
Results. Although the plasma expansion velocity is isotropically radial, driven by self-organization, the interactions between the conducting plasma and the magnetic fields steer the MHD system to a highly inhomogeneous spatial distribution. This spatial redistribution involves the focusing of magnetic energy, plasma thermal energy, and plasma momentum, generating collimated jets, point-symmetric knots, and circulating filaments.

Key words: magnetohydrodynamics (MHD) - methods: analytical - ISM: planetary nebulae: general

1 Introduction

A star with a mass of up to approximately eight solar masses evolves along the main sequence of the HR (temperature-luminosity) diagram to a stage where the accelerated fusion of the hydrogen layer, immediately above the helium core, causes the star to expand. The expansion is faster than the energy production rate, such that the expansion cools the star, moving it to the red giant branch and launching a slow wind. The expanding outer layers of the star are convective, and they bring the material in the fusing region below to the surface by turbulence. As the helium core ignites and expands, it pushes the hydrogen layer up and slows down the hydrogen fusion rate. The total energy generation decreases, and the star contracts in size and heats up in surface temperature, migrating it to the horizontal branch. As nuclear fusions continue, carbon and oxygen become the core material, surrounded by an inner helium layer and an outer hydrogen layer, plus a highly convective envelope, bringing the star to the asymptotic giant branch (AGB). The energy output is shifted to long wavelength emission towards the cool side of the HR diagram. Helium burning is quite sensitive to temperature, which causes burning instabilities and large star pulsations. The pulsations can become so strong that they give the outer layers of the star enough energy to eject matter as a fast wind, forming a planetary nebula. The remaining core of the star cools down and becomes a small but dense white drawf. Due to the rotational velocity of the star, we believe that strong magnetic fields are present in the star, especially in the helium and hydrogen layers, and in the outer convective envelope, and such fields have been actually detected in central stars (Jordan et al. 2005). Considering a rigidly rotating star and self-gravity of the distributed stellar mass, our preliminary research indicates that magnetic field lines in the star are wound in a torus, in layers. These magnetic field lines could be exposed during the pulsation phase, and released to the fast wind. The initial configuration of the field lines need not stay the same as the fast wind expands in space, because of its interaction with the conducting plasma and its transformation through dissipative processes.

As observational techniques and instruments become more sophisticated, the morphologies of planetary nebulae have become increasingly complex. As a result, the early pioneering model of interacting stellar winds for spherical planetary nebulae (Dyson & de Vries 1972; Kwok et al. 1978) has evolved to generalized interacting stellar winds to account for axisymmetric bipolar planetary nebulae (Morris 1981; Kahn & West 1985; Mellema 1997; Zhang & Kwok 1998). The bipolar features could be reproduced by the presence of a dense equatorial cloud, supposedly generated by the slow wind (Kahn & West 1985; Mellema et al. 1991; Mellema & Frank 1997; Dijkstra & Speck 2006), which has been imaged by high resolution instruments. Besides the bipolar features, images have also revealed point-symmetric knots and highly collimated jets (Miranda & Solf 1992; Lopez et al 1993; Balick et al. 1993). Under the generalized interacting stellar wind model, the possibility of the AGB star having a secondary close companion makes it possible to have an accretion disk around the secondary star to launch the jets (Morris 1981; Soker & Livio 1994; Mastrodemos & Morris 1999).

Nevertheless, these collimated jets and point-symmetric knots have led to the alternative model of magnetized wind-blown bubbles where an embedded azimuthal magnetic field in the fast stellar wind is responsable for collimated jets (Pascoli 1993; Chevalier & Luo 1994; Garcia-Segura 1997; Matt et al. 2000; Gardiner & Frank 2001). Point-symmetric knots can be reproduced with magnetized wind-blown bubbles by considering jets that precess about the polar axis generated by a tilted rotating star or a binary system at the center (Garcia-Segura & Lopez 2000). In particular, using the thin shell approximation of Giuliani (1982), Chevalier & Luo (1994) developed a self-similar model of planetary nebulae with a pre-selected constant radially expanding velocity. Observations over time have indicated that the winds have such extraordinarily high momentum and energy that they cannot be accounted for by the radiation pressure of the AGB star (Bujarrabal et al. 2001). These observational results appear to favor the magnetized wind-blown bubbles scenario where jets could be locally magnetically driven (Garcia-Segura et al. 2005) provided that there is a strong enough toroidal magnetic field. To explain the morphologies of planetary nebulae, current models attribute dynamic structures to the central star, or binary, to reproduce nebula morphologies in a deterministic way.

2 Self-similar methods and self-organization

To account for the highly collimating jets, point-symmetric knots, and spatially resolved filaments, we consider time-dependent full MHD equations and take an indirect form to describe the morphologies, through self-organizational processes via turbulent dissipation. From an analytic approach, the end configurations are disconnected from the initial conditions of the central engine, that are lost in the chaotic processes, although numerically we can always simulate over the turbulence, given enough computing power and time. In order to solve for the end configurations of the time-dependent MHD equations, we look for solutions of reduced dimensions in a self-similar form and with self-consistent evolution in time. Self-similarity between time and one spatial coordinate lowers the dimensions of the solution, and takes the MHD system through a specific time trajectory. We believe that this kind of self-similar solution, with self-consistent time evolution, is able to represent stable self-organized structures (Low 1982a,b, 1984a,b). This self-similar approach relies on the global conservation properties of MHD plasma to reach self-organized configurations, regardless of the initial conditions. This approach is hypothetical only, however we recall that Earth's magnetosphere is a stationary object in space with respect to the moving frame of Earth about the Sun. The magnetic field lines are stationary in such a frame. Nevertheless, near the terrestrial surface, the field lines co-rotate with Earth. In between, there is a transition layer that joins the stationary outer part to the co-rotating inner part of the magnetosphere. This is an example of dissipative processes in driving the end configuration, and where the initial conditions are lost in the process. The same can be said for the heliosphere. Therefore, on the scale of planetary nebulae, it is likely that the details of the AGB star initial conditions could be lost through dissipation.

Self-similarity, as a method, can be applied in different ways. In the case of disk winds (Blandford & Payne 1982; Pelletier & Pudritz 1992), steady state MHD equations are solved for self-similar solutions in space between two spatial labels in cylindrical coordinates to describe the transport of plasma from the equatorial disk to the polar axes. This amounts to a deterministic cause and effect approach through self-similar solutions in space. Under a specific magnetic field configuration of the disk, plasma transport with a continuous velocity field generates spatial singularities in the governing equation. In the case of the stellar envelope, Lou and his collaborators have treated an aggregating fluid under its self-gravitational field with a similarity variable x=r(t)/at and a pre-selected evolution function at, where a is the sound speed (Lou & Shen 2004; Bian & Lou 2005; Lou & Wang 2006; Lou & Gao 2006). Because of the linear dependence on time, this similarity variable x refers to a reference frame moving at speed a, which is different from the radial flow velocity. For this reason, the convective derivative remains explicit in the x representation. As a result, this approach amounts to finding the plasma structures in an adequate moving frame in the Eulerian x fluid description. Although an evolution function with a constant radial velocity for a spherically symmetric expanding plasma appears reasonable from a kinematic point of view, it is not true in most cases for plasmas. This is because the energy density is shared by the magnetic field, the kinetic energy of the plasma, and the thermal energy of the plasma. Consequently, the radial plasma velocity could slow down or speed up, depending on the interactions between the magnetic field and the conducting plasma. Without self-consistent time evolution of the system, this particular approach amounts to an application of self-similarity as a method.

The driving force for self-organization relies on the MHD global conservation properties through the quadratic invariants in the absence of dissipation (Hasegawa 1985; Biskamp 1993; Zhu et al. 1995; Yoshida & Mahajan 2002; Kondoh et al. 2004). These invariants are the total energy density, which includes magnetic energy, plasma kinetic energy, and plasma thermal energy, and the magnetic helicity, which is the scalar product of the vector potential and the magnetic field. In the case of incompressible fluids, there is also the cross-helicity, which is the scalar product of the plasma velocity and the magnetic field. In the presence of dissipation, these invariants could undergo constrained changes. These constrained variations could change one MHD configuration to another topologically accessible configuration of lower energy. These configurations are isolated configurations in the ideal case. The continuous topological transformation from one configuration to another is only possible through dissipative paths, such as magnetic reconnection. Mathematically speaking, this evolution amounts to the application of the variational principle on the total energy density under the constraint of global magnetic helicity conservation (Hasegawa 1985). The exact processes of dissipation in taking the system to self-organization need not be specified under the variational description.

Tsui (2008) has described the equatorial plasma torus of planetary nebulae from this self-organizational view, by using temporally self-similar MHD solutions with self-consistent evolution functions in spherical coordinates with radially expanding plasma flow, as a method and as a model. For the present case, we follow the same apporach but remove the axisymmetric assumption to seek three-dimensional self-similar solutions to model jets, knots, and filaments. Earlier works on temporally self-similar MHD with self-consistent evolution functions include one-dimensional (Osherovich et al. 1993, 1995) and two-dimensional (Tsui & Tavares 2005) cylindrical models of magnetic ropes in the solar-terrestrial environment, two-dimensional axisymmetric model of atmospheric ball lightning in spherical coordinates (Tsui 2006; Tsui et al. 2006), active galactic nucleus jets in astrophysical plasmas (Tsui & Serbeto 2007), traditionally treated as a steady state accretion-ejection MHD transport phenomenon from the accretion disk to the polar axis (Blandford & Payne 1982; Pelletier & Pudritz 1992). As for three-dimensional self-similar MHD solutions, there is one earlier case of interplanetary solar magnetic ropes (Gibson & Low 1998), which is not appropriate for astrophysical applications.

3 Self-similar formulation

The basic MHD equations in a Eulerian fluid description are given by

$${\partial\rho\over\partial t}+\nabla\cdot(\rho\vec v)~ =~0,$$ (1)

$$\rho\left\{{\partial\vec v\over\partial t} +(\vec v\cdot\nab... ... =~\vec J\times\vec B-\nabla p -\rho{GM\over r^{3}}\vec r ,$$ (2)

$${\partial\vec B\over\partial t}~ =~-\nabla\times\vec E~ =~\nabla\times(\vec v\times\vec B) ,$$ (3)

$$\nabla\times\vec B = \mu\vec J ,$$ (4)

$$\nabla\cdot\vec B = 0 ,$$ (5)

$${\partial\over\partial t}\left({p\over\rho^{\gamma}}\right) ... ...ec v\cdot \nabla\right)\left({p\over\rho^{\gamma}}\right)=0 .$$ (6)

Here, $\rho$ is the mass density, $\vec v$ is the bulk velocity, $\vec J$ is the current density, $\vec B$ is the magnetic field, p is the plasma pressure, $\mu$ is the free space permeability, $\gamma$ is the polytropic index, and M is the central mass that provides the gravitational field. This set of MHD equations describes the self-consistent interactions of a conducting fluid with the magnetic field. The fluid is governed by the Navia-Stoke equations coupled to the magnetic field through the $\vec J\times\vec B$ magnetic force. The magnetic field is governed by the slowly time varying Maxwell equations coupled to the fluid through the $\vec E+\vec v\times\vec B=0$ infinite conductivity condition. We have modelled the planetary nebula gas as a fully ionized plasma, although there is neutral material in the gas. If the gas is partially ionized, the neutral components could be collisionally coupled to the ionized component as a load. Consequently, the neutral transport properties, such as density and pressure profiles, are related to the ionized component.

We consider a radially expanding plasma and seek self-similar solutions in time where the time evolution is described by the dimensionless evolution function y(t). For this purpose, it is most convenient to think of a Lagrangian fluid description, and consider the position vector of a given laminar flow fluid element $\vec r(t)$. Under self-similarity, the radial profile is time invariant in terms of the radial label $\eta=r(t)/y(t)$, which has the dimension of r. Physically, $\eta$ is the Lagrangian radial position of a fixed fluid element. With a finite plasma, the domain of $\eta$ is bounded by mass conservation

$$0 < \eta_{\rm int} < \eta < \eta_{\rm ext}.$$ (7)

As for the plasma velocity, we consider self-similar structures derived from a spherically symmetric radial velocity which can be written as

$$\vec v = {{\rm d}\vec r(t)\over {\rm d}t}~ =\left\{\eta{{\rm d}y\over {\rm d}t}, 0, 0\right\}~ = \{v, 0, 0\}.$$ (8)

Our self-similar parameter $\eta$, defined through the Lagrangian fluid label, explicitly represents the fluid velocity by the time evolution function y(t). We note that y(0)=1 is the initial condition, so that $\eta=r(0)$ corresponds to the initial radial position. This initial plasma sphere, $(\eta_{\rm int},\eta_{\rm ext})$, corresponds to the distribution in space of the ejected fast wind from the hydrogen layer, together with the entrained magnetic field lines, throughout the entire strong pulsation period of the AGB star. This initial sphere then expands self-similarly in accordance with the evolution function y(t) such that $r(t)=\eta y(t)$. The evolution function will be solved self-consistently with respect to the spatial structures of the plasma. Since we are considering an isotropic radial plasma flow, we might wonder why structures other than unidimensional ones could arise. Indeed, for a purely hydrodynamic system without instabilities, only unidimensional configurations can emerge, with radially dependent mass density $\rho$ and plasma pressure p. Nevertheless, for MHD systems, complex structures with inhomogeneous mass density and plasma pressure will emerge, because of the interactions between the conducting plasma and the inherently two- or three-dimensional magnetic field.

The independent variables are now transformed from $(r,\theta,\phi,t)$ to $(\eta,\theta,\phi,y)$. We then determine the explicit dependence of y on each one of the physical variables with this radial velocity using functional analysis. First, making use of Eqs. (8), (1) becomes

$${\partial\rho\over\partial t} +{1\over r^2}{\partial\over\partial r}\left(r^2v\rho\right) \left({\partial\rho\over\partial t} +v{\partial\rho\over\partial r}\right) +\rho\left({\partial v\over\partial r}+{2v\over r}\right)$$  = 

$${\partial\rho\over\partial y}{{\rm d}y\over {\rm d}t} +{3\rho\ove... ...partial\rho\over\partial y}+{3\rho\over y}\right) {{\rm d}y\over {\rm d}t}~=~0.$$ = (9a)

To reach the second equality, we note that the first bracket in the first equality corresponds to the total time derivative of an Eulerian fluid element which amounts to the time derivative of a Lagrangian fluid element. As for the second bracket, it can be reduced by using $v={\rm d}r/{\rm d}t=\eta {\rm d}y/{\rm d}t$ and $\partial v/\partial r=(1/y)({\rm d}y/{\rm d}t)$. Solving this equation for the y scaling by separating the time part gives

$$\rho(\vec r,t) ={1\over y^3}\bar\rho(\eta,\theta,\phi) .$$ (9b)

Here, we have used the barred variables to denote the spatial part, where the temporal part is explicitly solved in terms of the evolution function. As for Eq. (6), with $\alpha_{0}F=(p/\rho^{\gamma})$ where  $\alpha_{0}$ is a constant that carries the physical dimension so that F is a dimensionless function, it follows

$${\partial F\over\partial t}+v{\partial F\over\partial r} = {\partial F\over\partial y}{{\rm d}y\over {\rm d}t} =0,$$ (10a)

$$\left({p\over\rho^{\gamma}}\right) = \alpha_{0}F(\vec r,t) = {1\over y^0}\alpha_{0}\bar F(\eta,\theta,\phi) .$$ (10b)

As for Eq. (3), with the aid of Eq. (5), the magnetic fields are

$${\partial B_{r}\over\partial t} +v{1\over r^{2}}{\partial\over\pa... ...B_{r}\over\partial t} +v{\partial B_{r}\over\partial r}\right)+{2v\over r}B_{r}$$

$$\qquad \qquad={\partial B_{r}\over\partial y}{{\rm d}y\over {\rm ... ...rtial B_{r}\over\partial y}+{2B_{r}\over y}\right) {{\rm d}y\over {\rm d}t} =0,$$ (11a)

$$B_{r}(\vec r,t) ={1\over y^2}\bar B_{r}(\eta,\theta,\phi) ,$$ (11b)

$${\partial B_{\theta}\over\partial t} +{1\over r}{\partial\over\pa... ...heta}\over\partial r}\right) +{1\over r}B_{\theta}{\partial\over\partial r}(rv)$$

$$\qquad \qquad= {\partial B_{\theta}\over\partial y}{{\rm d}y\over... ...theta}\over\partial y}+{2B_{\theta}\over y}\right) {{\rm d}y\over {\rm d}t} =0,$$ (12a)

$$B_{\theta}(\vec r,t)~ =~{1\over y^2}\bar B_{\theta}(\eta,\theta,\phi) ,$$ (12b)

$${\partial B_{\phi}\over\partial t} +{1\over r}{\partial\over\partial r}( rvB_{\phi}) =0,$$ (13a)

$$B_{\phi}(\vec r,t) ={1\over y^2}\bar B_{\phi}(\eta,\theta,\phi) .$$ (13b)

Making use of Eq. (9b), we derive the plasma pressure from Eq. (10b)

$$p(\vec r,t) {1\over y^{3\gamma}}\alpha_{0}\bar F(\eta,\theta,\phi) \bar\rho^{\gamma}(\eta,\theta,\phi)$$ =

$${1\over y^{3\gamma}}\alpha_{0}\bar F\bar\rho^{\gamma} = {1\over y^{3\gamma}}\bar p(\eta,\theta,\phi) .$$ = (14)

Since Eqs. (12a) and (13a) are of the same form, we conclude that, under self-similarity, $\bar B_{\phi}$ is a linear function of  $\bar B_{\theta}$ with

$$\bar B_{\phi}~=~k\bar B_{\theta} .$$ (15)

Making use of Eq. (4) to eliminate the current density in Eq. (2), we obtain the momentum equation which has three components. The $\phi$, $\theta$, and r components are respectively

$$\bar B_{\theta} \left[{\partial\over\partial\theta}(k\bar B_{\the... ...rtial\phi} -\sin\theta{\partial\over\partial\eta}(\eta k\bar B_{\theta})\right]$$

$$\qquad\qquad -y^{4-3\gamma}{\partial\over\partial\phi}(\mu\bar p) = 0 ,$$ (16)

$$k\bar B_{\theta} \left[{\partial\over\partial\theta}(k\bar B_{\th... ...rtial\theta} -\sin\theta{\partial\over\partial\eta}(\eta\bar B_{\theta})\right]$$

$$\qquad\qquad +y^{4-3\gamma}\sin\theta{\partial\over\partial\theta}(\mu\bar p) =0,$$ (17)

$$k\bar B_{\theta} \left[{1\over\eta\sin\theta}{\partial\bar B_{r}\... ...eta\bar B_{\theta}) -{1\over\eta}{\partial\bar B_{r}\over\partial\theta}\right]$$

$$\qquad\qquad -y^{4-3\gamma}{\partial\over\partial\eta}(\mu\bar p) =\mu\bar\rho y^2{{\rm d}^2y\over {\rm d}t^2} +\mu\bar\rho{GM\over\eta^2} \cdot$$ (18)

We have reduced the general set of time-dependent ideal MHD equations, Eqs. (1)-(6), to a set of self-similar equations with appropriate time scalings, Eqs. (7)-(13). The general ideal MHD set has nonlinear terms of convective type  $(\vec v\cdot\nabla)$. By using the fluid label description, the  $(\vec v\cdot\nabla)$ convective terms are absorbed in the Lagrangian time derivative representation. The structure of the nonlinear terms, absorbed in the Lagrangian fluid label formulation, will appear in the $\eta$ profile of the system.

4 Monopolar jet structures

After this self-similar formulation, we have to solve Eqs. (16)-(18) for the self-similar configurations. For this, we first separate the radial variable from the other two variables by writing

$$\bar B_{r}(\eta,\theta,\phi) =A_{0}R(\eta)\tilde B_{r}(\theta,\phi) ,$$ (19)

where we have used tilded variables for the angular part, and likewise for  $\bar B_{\theta}$ and $\bar B_{\phi}=k\bar B_{\theta}$. Furthermore, we take

$$\bar p(\eta,\theta,\phi) = p_{0}R^2(\eta)\tilde p(\theta,\phi) .$$ (20)

Here, $R(\eta)$, $\tilde B_{r}(\theta,\phi)$, and $\tilde p(\theta,\phi)$ are dimensionless functions, and A₀ and p₀ carry the dimensions of magnetic field and pressure respectively. Since magnetic pressure is a quadratic quantity of the magnetic field, we have used this to write the plasma pressure as $R^2(\eta)$. We take A₀=1 for unit amplitude magnetic fields, such that p₀ is relative to this amplitude. Specifically, we take

$$R(\eta)=(a\eta)^{-n},$$ (21)

to represent a power law field decaying with distance, where a is a normalizing parameter of $\eta$. Considering

$$(4-3\gamma)~=~0,$$

Eqs. (16), (17) are respectively

$$\tilde B_{\theta} \left[{\partial\over\partial\theta}(k\tilde B_{\theta}\sin\theta)... ...theta}\over\partial\phi}\right] -\tilde B_{r}(n-1)\sin\theta k\tilde B_{\theta}$$

$$\qquad -{\partial\over\partial\phi} \left({1\over 2}\tilde B^2_{r}+\mu p_{0}\tilde p\right) =0,$$ (22)

$$k\tilde B_{\theta} \left[{\partial\over\partial\theta}(k\tilde B_... ...\theta}\over\partial\phi}\right] +\tilde B_{r}(n-1)\sin\theta \tilde B_{\theta}$$

$$\qquad\qquad +\sin\theta{\partial\over\partial\theta} \left({1\over 2}\tilde B^2_{r}+\mu p_{0}\tilde p\right) =0,$$ (23)

We now separate the azimuthal dependence by writing

$$\tilde B_{r}(\theta,\phi) =\Theta_{r}(\theta)\Phi(\phi) ,$$ (24a)

$$\tilde B_{\theta}(\theta,\phi) = \Theta(\theta)\Phi(\phi) ,$$ (24b)

$$\tilde p(\theta,\phi) = \Theta^2_{p}(\theta)\Phi^2(\phi) ,$$ (24c)

where $\tilde B_{\phi}=k\tilde B_{\theta}$. These functional dependences give Eq. (22) as

$$\Theta{\partial\over\partial\theta}(k\Theta\sin\theta) -\Theta_{r... ...{0}\Theta^2_{p}\right)\right] {1\over\Phi}{\partial\Phi\over\partial\phi} \cdot$$ (25a)

Considering the separation constant im, such that

$${\partial\Phi\over\partial\phi}~=~{\rm i}m\Phi ,$$ (26b)

$$\Phi(\phi)~=~{\rm e}^{+{\rm i}m\phi} = \cos m\phi + {\rm i}\sin m\phi,$$ (25)

we then have

$$\Theta{\partial\over\partial\theta}(\Theta\sin\theta) -\Theta_{r}... ...ver k}m\left[\Theta^2-\left(\Theta^2_{r}+2\mu p_{0}\Theta^2_{p}\right)\right] ,$$ (26)

which identifies $k=+\rm i$. Following the same procedures, Eq. (23) reads

$$\Theta{\partial\over\partial\theta}(\Theta\sin\theta) -\Theta_{r}... ...\partial\over\partial\theta} \left(\Theta^2_{r}+2\mu p_{0}\Theta^2_{p}\right) .$$ (27)

The left sides of these two equations are the same, which allows the right sides to be equated to give

$$\sin\theta{\partial\Theta^2_{*}\over\partial\theta}~ =~-2m\Theta^2_{*} ,$$ (28a)

$$\Theta^2_{*} = \left(\Theta^2_{r}+2\mu p_{0}\Theta^2_{p}\right).$$ (28b)

To solve for $\Theta^2_{*}$, we integrate Eq. (28a) to get

$$\ln \Theta^2_{*} -2m\int {{\rm d}\theta\over\sin\theta} =-2m\int {\sin\theta {\rm d}\theta\over\sin^2\theta}$$ =

$$+\ln \left({1+x\over 1-x}\right)^m,$$ = (29)

where we have multiplied and divided the right side by $\sin\theta$ to implement the integration, and $x=\cos\theta$. This solution of $\Theta^2_{*}$ is singular at x=+1 for positive m, and x=-1 for negative m. Such a solution results in jet features on the magnetic field lines and plasma density. We note that  $\Theta^2_{*}$ is a composite variable that contains  $\Theta^2_{r}$ of the radial magnetic field and $\Theta^2_{p}$ of the plasma pressure. This composite variable $\Theta^2_{*}$, therefore, reflects the partition of the magnetic energy and the plasma thermal energy through p₀.

To get $\Theta$, instead of solving either Eqs. (26) or (27), we make use of Eq. (5) to get

$$\nabla\cdot\vec B ={1\over y^3}\left({1\over\eta}R(\eta)\right) \... ...B_{\theta}) +{\partial\over\partial\phi}(k\tilde B_{\theta})\right]\Bigg\} = 0,$$ (30a)

$${\partial\over\partial\theta}(\sin\theta\Theta) = m\Theta+(n-2)\sin\theta\Theta_{r}.$$ (30b)

Subsituting Eqs. (30b) into (26) and (27) respectively gives

$$\sin\theta\Theta\Theta_{r}~=~m\Theta^2_{*} ,$$ (31a)

$$2\Theta\Theta_{r} =-{\partial\over\partial\theta}\Theta^2_{*} .$$ (31b)

It can be shown readily that these two equations are equivalent to Eq. (28a). We note that Eq. (30b) couples $\Theta$ with  $\Theta_{r}$. To decouple these two functions, we consider the special case of

n=2. (32)

Coincidentally, although this gives an inverse squared dependence on distance of $R(\eta)=(a\eta)^{-2}$ by Eq. (21), this is not the motivation for choosing n=2. The basic motivation is to decouple Eq. (30b) from $\Theta_{r}$. The case of $n\neq 2$ is discussed in Sect. 7. Multiplying by $\sin\theta$ and defining $P(\theta)=\sin\theta\Theta$, Eq. (30b) can be integrated to give

$$P(x) = \left({1-x\over 1+x}\right)^{m/2} ,$$ (33a)

$$\Theta(x) = {(1-x)^{(m-1)/2}\over (1+x)^{(m+1)/2}},$$ (33b)

$$\Theta_{r}(x)~=~m\left({1+x\over 1-x}\right)^{3m/2},$$ (33c)

$$2\mu p_{0}\Theta^2_{p}(x) = \left[\left({1+x\over 1-x}\right)^{m}-m^2\left({1+x\over 1-x}\right)^{3m}\right] > 0 ,$$ (33d)

where Eq. (33c) is obtained from Eq. (31a), and $\Theta ^{2}_{p}(x)$ can be recovered from $\Theta ^{2}_{*}(x)$ of Eq. (28b). From these solutions, we see that B_r and p are singular at x=+1 through $\Theta _{r}(x)$ and $\Theta ^{2}_{p}(x)$, whereas $B_{\theta}$ and $B_{\phi}$ are singular at x=-1 through $\Theta (x)$. A singular magnetic field by itself, such as a dipole field at the center, is acceptable. To avoid unphysical results, we require the magnetic flux be finite even though the magnetic field is infinite. In other words, we require the singularities be integrable in x, which demands the power of the singularities be less than unity. By inspection of the terms, we conclude that m<1/3. Let us take, as an example,

$$m={1\over 4},$$ (34)

to show $\Theta ^{2}_{*}(x)$, $\Theta _{r}(x)$, and $\Theta (x)$ in Figs. 1-3 respectively. In Fig. 1, the amplitude of the physical quantity $\Theta ^{2}_{*}(x)$ is plotted against the polar angle $\theta$ so that the projections of the plot on the x and y axes give $\Theta^2_{*}(x)\sin\theta$ and $\Theta^2_{*}(x)\cos\theta$ respectively. In this polar plot, the magnitude of the physical variable is given by the radial length from the origin. The y and x axes in the figure stand for the polar z axis and the projected radius on the equatorial plane of the planetary nebula.

	Figure 1: The function $\Theta ^{2}_{*}(x)$ with $x=\cos\theta$ against the polar angle $\theta$ with projections on the vertical and horizontal axes $\Theta^2_{*}(x)\cos\theta$ and $\Theta^2_{*}(x)\sin\theta$ respectively, indicating an integrable singularity of finite flux at x=+1 for a jet structure.
Open with DEXTER

	Figure 2: The function $\Theta _{r}(x)$ with $x=\cos\theta$ in a polar plot indicating an integrable singularity at x=+1 for a jet structure.
Open with DEXTER

Figures 1 and 2 show an integrable singularity at x=+1, or $\theta=0$. Furthermore, Fig. 3 shows a weaker singularity at x=+1 than the one of Fig. 2. However, Fig. 3 shows another singularity at x=-1, while Fig. 2 is regular at that location. As for $\Theta ^{2}_{p}(x)$, the second term of Eq. (33d) exceeds the first term when

$$m\left({1+x\over 1-x}\right)^m>1,$$

or as x gets very close to unity. Writing $x=(1-\delta)$, we get $(\delta/2)^m<m$, or $(\delta/2)<m^{1/m}$. With Eq. (34), we have $\delta=(1/2)^7$. Consequently, $\Theta ^{2}_{p}(x)$ gets smaller as x approaches unity, and it vanishes at $x=(1-\delta)$. Beyond this point, it becomes negative, giving a jet structure as in Fig. 4. This jet structure is along the radial magnetic field given by $\Theta _{r}(x)$. This negative sign can be absorbed in the $\Phi^2(\phi)$ factor of Eq. (24c) as a phase shift by writing

$$\tilde p(\theta,\phi) =\vert\Theta^2_{p}(\theta)\vert({\rm i}... ...)\big\vert\big\vert{\rm e}^{+{\rm i}(m\phi+\pi/2)}\big\vert^2.$$

	Figure 3: The function $\Theta (x)$ with $x=\cos\theta$ indicating integrable singularities at x=+1 and at x=-1 for knot structures.
Open with DEXTER

5 Wind momentum and energy

As for the radial component, with $(4-3\gamma)=0$, Eq. (18) reads

$$\mu\bar\rho\eta y^2{{\rm d}^2y\over {\rm d}t^2} +\mu\bar\rho{GM\o... ...rtial\tilde B_{r}\over\partial\phi} +{\partial B_{r}\over\partial\theta}\right]$$

$$~\qquad\quad+2n\left({1\over\eta}R^2\right)\left(\mu p_{0}\tilde ... ...\left[-m\Theta_{r}+\sin\theta{\partial\Theta_{r}\over\partial\theta}\right]\Phi$$

$$~ \qquad\quad+2n\left({1\over\eta}R^2\right)\mu p_{0}(\Theta_{p}\Phi)^2.$$ (35)

With $\alpha$ as the separation constant, we have

$$\left({1\over\eta}R^2\right) \left \{\tilde B_{\theta} \left[-m\T... ...er\partial\theta}\right] \Phi +2n\mu p_{0}\left(\Theta_{p}\Phi\right)^2\right\}$$

$$\qquad \qquad=\mu\bar\rho\left({GM\over\eta^2}+\alpha\eta\right) ,$$ (36a)

$$\qquad \qquad {{\rm d}^2y\over {\rm d}t^2}~=~{\alpha\over y^2} \cdot$$ (36b)

We write the mass density as

$$\bar\rho(\eta,\theta,\phi)~ =~\rho_{0}\textbf{R}(\eta)\tilde\rho(\theta,\phi),$$ (37)

where $\rho_{0}$ carries the dimension of mass density, and $\textbf{R}(\eta)$ and $\tilde\rho(\theta,\phi)$ are dimensionless functions. We can identify immediately from Eq. (35a) that

$$\textbf{R}(\eta)~ =~{R^2\over (GM/\eta+\alpha\eta^2)} ,$$ (38a)

$$\mu\rho_{0}\tilde\rho \Bigg\{\tilde B_{\theta} \left[-m\Theta_{r}+\sin\theta{\partial\Theta_{r}\over\partial\theta}\right] \Phi +2n\mu p_{0}(\Theta_{p}\Phi)^2\Bigg\}$$ =

$$\left[2n\mu p_{0}\Theta^2_{p}-2m\Theta\Theta_{r}\right]\Phi^2 .$$ = (38b)

The second equation reads

$$\mu\rho_{0}\tilde\rho = \Bigg \{n\left[\left({1+x\over 1-x}\right... ...t)^{3m}\right] -2m^2{(1+x)^{(2m-1)/2}\over (1-x)^{(2m+1)/2}}\Bigg \}\Phi^2 >0 .$$ (39)

This mass density has a singularity at x=+1. With m=1/4, we note that the numerator of the last term in Eq. (38), (1+x)^(2m-1)/2, has a negative power. This gives a singularity at x=-1. Since the power of this singularity is -1/4, it is also integrable. The mass density distribution is shown in Fig. 5 with a jet structure in the x=-1 direction, because of the radial magnetic field. If we consider m<0, Figs. 1-5 would be inverted. Consequently, the jet structures would be on both sides of the polar axis.

	Figure 4: The function $\Theta ^{2}_{p}(x)$ with $x=\cos\theta$ indicating a plasma pressure jet in a very narrow cone about x=+1.
Open with DEXTER

	Figure 5: The function $\tilde\rho(x)$ with $x=\cos\theta$ indicating a mass density jet in a very narrow cone about x=+1.
Open with DEXTER

For the time evolution function of Eq. (36b), we multiply over by d $y/{\rm d}t$ to get

$$\left ({{\rm d}y\over {\rm d}t}\right)^2 =2\left(H-{\alpha\over y}\right) > 0 ,$$ (40a)

where H is an integration constant which stands for the energy density of the system. This equation is subject to the initial condition of y(0)=1 by definition of the Lagrangian fluid label. The solution of the evolution function depends on the nature of the physical system, represented by $\alpha$ and H. Considering that the stellar wind is launched above the escape velocity together with the magnetic fields, H would be positive. This implies that planetary nebulae are open systems that expand indefinitely. Let us consider $\alpha$ positive, such that $\textbf{R}(\eta)$ is a monotonically decreasing function. Since the right side of Eq. (40a) has to be positive, we conclude that $H>\alpha>0$. This gives the plasma velocity a slow start with $({\rm d}y/{\rm d}t)^2=2(H-\alpha)$ that finishes with a fast terminal velocity

$$\left({{\rm d}y\over {\rm d}t}\right)^2=2H,$$ (40b)

which is qualitatively compatible with observations of high momentum and energy (Bujarrabal et al. 2001).

We have modelled the stellar wind by an isotropic radial flow which could be driven by the stellar magnetic field with unspecified topology. Driven by self-organization, the interactions between the conducting plasma and the magnetic fields steer the MHD system to a highly inhomogeneous spatial distribution with plasma jets in the polar directions, although the plasma velocity is radially symmetric. This redistribution in space among the magnetic energy, plasma kinetic energy, and plasma thermal energy is the basic mechanism that drives the jet momentum, expressed through the time evolution function. For this reason, we are dealing with a spatially focusing effect, where energies are rearranged in space, instead of a plasma accelerating mechanism.

6 Knot and filament structures

With the momentum equation solved, the magnetic fields are given by

$$\bar B_{r}=+R(\eta)\Theta_{r}\cos m\phi,$$ (41a)

$$\bar B_{\theta} = +R(\eta)\Theta\cos m\phi ,$$ (41b)

$$\bar B_{\phi} = -R(\eta)\Theta\sin m\phi .$$ (41c)

The magnetic field lines are given by

$${\Theta_{r}\cos m\phi\over a{\rm d}\eta}~ =~{\Theta\cos m\ph... ...\theta} =-{\Theta\sin m\phi\over a\eta\sin\theta \rm d\phi},$$ (42a)

which can be written as

$${a{\rm d}\eta\over a\eta} -m{(1+x)^{2m}\over (1-x)^{2m}}{\rm d}x$$ =

$$-m{(1+x)^{2m+1}\over (1-x)^{2m-1}}{{\rm d}\phi\over\tan m\phi} \cdot$$ = (42b)

With $\eta_{0}$ as an integration constant, the first equality gives

$$\ln{a\eta\over a\eta_{0}} =-m\int {(1+x)^{2m}\over (1-x)^{2m}} {\rm d}x,$$ (43a)

which can be integrated numerically, as is shown in Fig. 6. Since the singularity is integrable, $a\eta$ is finite at x=+1. This result in the radial $\eta$ can then be converted to radial position r(t) through the self-consistently determined evolution function y(t) by $r(t)=\eta y(t)$.

The second equality corresponds to

$${{\rm d}\theta\over\sin\theta} =-{ {\rm d}\phi\over\tan m\phi},$$

and can be integrated to give

$${(1+x)^{m}\over (1-x)^{m}} =K(\sin m\phi)^2,$$ (43b)

where K is an integration constant. With the interval of x between (-1,+1), or $\theta$ between $(\pi,0)$, the left side of Eq. (43b), (1+x)^m/(1-x)^m labelled on the left axis, covers an interval $(0,\infty)$ with positive m, and is plotted in Fig. 7 against x between (-1,+1) labelled on the bottom axis. The right side, $K(\sin m\phi)^2$ labelled on the right axis, is also plotted against $m\phi$ on the top axis, with a scale between 0 and $\pi/2$, in the same figure. In order to view the mapping between the right side and the left side, we assign a large constant K. With m=1/4 and K=3, as x departs from -1, or $\theta$ from $\pi$, $m\phi$ departs from 0, and it maps a root of $\phi$. As x approaches +1, or $\theta$ approaches 0, $m\phi$ reaches $\pi/2$, which takes $\phi$ over $(0,2\pi)$. On the return path of the field lines, x decreases from +1 back to -1, bringing $m\phi$ from $\pi/2$ to $\pi$ along the descending branch of $K(\sin m\phi)^2$, which takes $\phi$ over $2\pi,4\pi$. This descending branch, which is the continuation of Fig. 7, is not shown.

	Figure 6: The magnetic field line $a\eta -\theta$ dependence.
Open with DEXTER

	Figure 7: The magnetic field line $\phi -\theta$ dependence.
Open with DEXTER

To summarize, the field lines starting at x=-1 and $a\eta=+1$ go through the x=(-1,+1,-1) cycle once, with $a\eta=(+1,+2.3,+1)$, while completing the $m\phi=(0,\pi)$ cycle once, covering $\phi=(0,4\pi)$, before closing on themselves again. This generates helical field lines, as shown in Fig. 8, on the surface of revolution of Fig. 6. Since the singularities at $x=\pm 1$ are integrable, and also because of the circulating nature of the fields  $\bar B_{\theta}$ and  $\bar B_{\phi}$, the magnetic field lines converge to $a\eta=1$ at x=-1 axis and to $a\eta=2.3$ at x=+1 axis, as shown in Fig. 6. These locations correspond to point-symmetric magnetic knots, where the field strength is infinite. The filaments correspond to the helical magnetic field lines in space, as shown in Fig. 8. These same field lines have a different shape when they are viewed at different orientations, such as in Figs. 9 and 10. These field lines in terms of $\eta$ expand in space with radial position r(t) according to $r(t)=\eta y(t)$. Further field lines can be generated with, for example, m=1/3.5=2/7=4/14. In this case, the field lines starting at x=-1 and $a\eta=+1$ go through the x=(-1,+1,-1) cycle once, with $a\eta=(+1,+2.3,+1)$, while completing the $m\phi=(0,\pi/2,\pi)$ cycle once, covering $\phi=(0,7\pi/4,7\pi/2)$. Since the lowest $2\pi$ multiple of the $\phi$ cycle is four, the field lines have to complete four cycles of x and of $m\phi$ to make $\phi$ cover $(0,7\pi,14\pi)$, such that the field lines can close on themselves again.

	Figure 8: The three-dimensional magnetic field lines, wound on the surface of revolution of Fig. 6, viewed parallel to the x-y plane at about 45 degrees.
Open with DEXTER

	Figure 9: The three-dimensional magnetic field lines, wound on the surface of revolution of Fig. 6, viewed down the z axis.
Open with DEXTER

A range of $B_{\phi}$ in space corresponds to a range of a. With a range of $a\eta_{0}=1$ of Eq. (43a), we have a distribution of $\eta_{0}$, and we can fill the space with shells of field lines of some m<1/3 up to $a\eta_{\rm ext}$ of Eq. (7), and generate a sequence of magnetic knots on the downward polar axis, $\theta=\pi$. Despite the radial component at x=+1 giving the plasma jet structures, the field lines are closed at a line segment $a\eta=2.3$ on the upward polar axis, $\theta=0$, because of the circulating nature of the meridian and azimuthal components. With m<0, the mirror images of the knots and field lines can be superimposed on those with m>0, generating more lines of knots and concentric shells of magnetic field lines.

7 General n case, symmetric knots

We now solve Eq. (30b) for an arbitrary n. Combining Eqs. (30b) and (31a), we get

$$\left(1-x^2\right){{\rm d}\over {\rm d}x}P^2(x)+2mP^2(x)= -2m(n-2) \left(1-x^2\right)\Theta^2_{*}(x) ,$$ (44a)

where $P(x)=(1-x^2)^{1/2}\Theta(x)$ and $\Theta ^{2}_{*}(x)$ is given by Eq. (29). If the boundary condition is known, this equation can be integrated numerically to get

$$P^2(x) = -2m\int^{x}_{-1}\left[{P^2(x)\over (1-x^2)} +(n-2)\Theta^2_{*}(x)\right]{\rm d}x +P^2(-1).$$ (44b)

To obtain the boundary condition at x=-1, we note that the right side of Eq. (44a) is equal to -2m(n-2)(1+x)^m+1/(1-x)^m-1. Since m<1 but positive, this term vanishes at x=-1 and at x=+1 as well. As a result, in the neighborhood of $x=\pm 1$, P²(x) is described by the homogeneous version of Eq. (44a), where the right side is null. The homogeneous solution can be solved readily as

$$P^2(x)=\left({1-x\over 1+x}\right)^{m},$$ (45)

which provides the needed boundary condition. The self-similar functions are, therefore, given by

$$\Theta(x) = {P(x)\over\sin\theta},$$ (46a)

$$\Theta_{r}(x) = {m\Theta^2_{*}(x)\over P(x)},$$ (46b)

$$2\mu p_{0}\Theta^2_{p}(x) = \left[\Theta^2_{*}(x)-\Theta^2_{r}(x)\right] > 0 .$$ (46c)

The magnetic field lines are described by

$$\ln{a\eta\over a\eta_{0}} =-m\int {\Theta^2_{*}(x)\over P^2(x)} {\rm d}x,$$ (47a)

$${{\rm d}\theta\over\sin\theta} =-{ {\rm d}\phi\over\tan m\phi}\cdot$$ (47b)

The second equation remains the same, while the first equation can be integrated numerically to get the field lines. With n=3, the corresponding self-similar functions of Figs. 2-6 are evaluated anew, and are presented in Figs. 11-15. From Figs. 12 and 15, we can see that the source term in Eq. (44a) makes $\Theta (x)$ and the $(\eta-x)$ mapping more symmetric. The corresponding magnetic field lines of Figs. 8-10 are also shown in Figs. 16-18, with knots, $x=\pm 1$ and $a\eta=+1$, symmetrically placed at equal distances from the center.

	Figure 10: The three-dimensional magnetic field lines, wounded on the surface of revolution of Fig. 6, viewed along the x axis.
Open with DEXTER

	Figure 11: The function $\Theta _{r}(x)$, with n=3, indicating an integrable singularity at x=+1 for a jet structure.
Open with DEXTER

	Figure 12: The function $\Theta (x)$, with n=3, indicating a more symmetric structure.
Open with DEXTER

	Figure 13: The function $\Theta ^{2}_{p}(x)$, with n=3, indicating a plasma pressure jet in a very narrow cone about x=+1.
Open with DEXTER

	Figure 14: The function $\tilde\rho(x)$ with n=3 indicating a mass density jet in a very narrow cone about x=+1.
Open with DEXTER

	Figure 15: The magnetic field line $a\eta -\theta$ dependence, with n=3, showing a symmetric structure.
Open with DEXTER

	Figure 16: The three-dimensional magnetic field lines, with n=3, viewed parallel to the x-y plane at about 45 degrees.
Open with DEXTER

	Figure 17: The three-dimensional magnetic field lines, with n=3, viewed down the z axis.
Open with DEXTER

	Figure 18: The three-dimensional magnetic field lines, with n=3, viewed along the x axis.
Open with DEXTER

8 Discussions and conclusions

We have presented an analysis that describes planetary nebula morphologies as three-dimensional self-organized structures. Analytically, because of self-organization, the morphologies cannot be linked to the central AGB star in a deterministic way, although numerically we can always simulate turbulence from the beginning to reach the end states, given enough computing power and time, as in meteorological models of hurricane formation. Because of our mostly Newtonian scientific tradition with a direct and explicit link between cause and consequence, we are not accustomed to see phenomena from a self-organizational perspective. The situation is somewhat analogous to dealing with uncertainty principle and duality of light in quantum mechanics. In order to get the end configurations without specifying the initial conditions, we have used self-similar MHD solutions with self-consistent evolution functions to represent the self-organized states emerging from turbulence.

In our model, the key role of plasma pressure is represented in the momentum equation, Eq. (2), which is solved together with the magnetic field. Although the plasma velocity is isotropically radial, due to self-organization, complex three-dimensional structures can arise because of the interactions between magnetic fields and the conducting plasma. These structures correspond to the spatial redistribution of plasma and magnetic fields, which amounts to focusing of magnetic energy, plasma thermal energy, and plasma momentum in space. With this self-organization model through self-similar configurations, we have revealed that magnetic fields and plasma are structured in a narrow cone along the polar axis. The magnetic field strength, plasma pressure, and mass density are singular there but integrable. The plasma pressure and mass density along the jet decrease monotonically with distance according to a power law. The resulting plasma jets are under the local pressure of the magnetic field, and, as a result, they can carry large amounts of momentum and energy. Some examples of these energetic jets appear in the M2-9 Twin Jet Nebula, CRL 2688 Egg Nebula, NGC 3242, NGC 6826, NGC 7009. The magnetic fields have their lines of force circulating in space forming filaments that converge at given locations on the ploar axis, giving point-symmetric knots. The radial field, which is also singular on the axis, preserves the divergence-free nature of the magnetic knots. Examples of knots can be found in the M2-9 Twin Jet Nebula, NGC 5307, and filaments in MyCn 18 Hourglass Nebula, NGC 6543 Cat's Eye Nebula, NGC 2392 Eskimo Nebula, M2-9 Twin Jet Nebula, NGC 6543.

In a recent publication (Tsui & Serbeto 2007) of an axisymmetric two-dimensional model, extragalactic polar jets are described from the same self-organizational perspective through an eigenvalue equation with regular eigenfunctions. Jets are driven by plasma pressure to progressively collimate along the polar axis. In particular, most of the ejected mass is distributed incompact plasma tori along the polar jets. The fact that jets, knots, and other sturctures are ubiquitous and reproducible phenomena on stellar and galactic scales has led us to believe that there is an underlying order by self-organization. Without the astrophysical observations, these configurations would be just exercises of self-similar mathematics. Under the matching morphologies, we may see planetary nebulae as the result of self-organization, rather than being direct consequences of the AGB central engine. Self-similarity with self-consistent evolution functions is sufficient to reveal the structures.

Acknowledgements

The author is deeply grateful to Dr. B. C. Low for the inspiring thoughts and physical insights of self-similar solutions, and to Prof. Akira Hasegawa for the very essential concept of self-organization in fluids and plasmas.

References

• Balick, B., Rugers, M., Terzian, Y., & Ghengalur, J. N. 1993, ApJ, 411, 778 [NASA ADS] [CrossRef] (In the text)
• Bian, F. Y., & Lou, Y. Q. 2005, MNRAS, 363, 1315 [NASA ADS] [CrossRef] (In the text)
• Biskamp, D. 1993, Nonlinear Magnetohydrodynamics (Cambridge: Cambridge University Press), Chap. 2 (In the text)
• Blandford, R. D., & Payne, D. G. 1982, MNRAS, 199, 883 [NASA ADS] (In the text)
• Bujarrabal, V., Castro-Carrizo, A., Alcolea, J., & Contreras, C. S. 2001, A&A, 377, 868 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
• Chevalier, R. A., & Luo, D. 1994, ApJ, 421, 225 [NASA ADS] [CrossRef] (In the text)
• Dijkstra, C., & Speck, A. K. 2006, ApJ, 651, 288 [NASA ADS] [CrossRef] (In the text)
• Dyson, J. E., & de Vries, J. 1972, A&A, 20, 223 [NASA ADS] (In the text)
• Gardiner, T. A., & Frank, A. 2001, ApJ, 557, 250 [NASA ADS] [CrossRef] (In the text)
• Garcia-Segura, G. 1997, ApJ, 489, L189 [NASA ADS] [CrossRef] (In the text)
• Garcia-Segura, G., & Lopez, J. L. 2000, ApJ, 544, 336 [NASA ADS] [CrossRef] (In the text)
• Garcia-Segura, G., Lopez, J. L., & Franco, J. 2005, ApJ, 618, 919 [NASA ADS] [CrossRef] (In the text)
• Gibson, S. E., & Low, B. C. 1998, ApJ, 493, 460 [NASA ADS] [CrossRef] (In the text)
• Giuliani, J. L., Jr. 1982, ApJ, 256, 624 [NASA ADS] [CrossRef] (In the text)
• Hasegawa, A. 1985, Adv. Phys., 34, 1 [NASA ADS] [CrossRef] (In the text)
• Jordan, S., Werner, K., & OToole, S. J. 2005, A&A, 432, 273 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
• Kondoh, Y., Serizawa, S., Nakano, A., Takahashi, T., & Van Dam, J. W. 2004, Phys. Rev. E, 70, 066312 [NASA ADS] [CrossRef] (In the text)
• Kahn, F. D., & West, K. A. 1985, MNRAS, 212, 837 [NASA ADS] (In the text)
• Kwok, S., Purton, C. R., & FitzGerald, P. M. 1978, ApJ, 219, L125 [NASA ADS] [CrossRef] (In the text)
• Lopez, J. A., Meaburn, J., & Palmer, J. W. 1993, ApJ, 415, L135 [NASA ADS] [CrossRef] (In the text)
• Lou, Y. Q., & Shen, Y. 2004, MNRAS, 348, 717 [NASA ADS] [CrossRef] (In the text)
• Lou, Y. Q., & Gao, Y. 2006, MNRAS, 373, 1610 [NASA ADS] [CrossRef] (In the text)
• Lou, Y. Q., & Wang, W. G. 2006, MNRAS, 372, 885 [NASA ADS] [CrossRef]
• Low, B. C. 1982a, ApJ, 254, 796 [NASA ADS] [CrossRef] (In the text)
• Low, B. C. 1982b, ApJ, 261, 351 [NASA ADS] [CrossRef]
• Low, B. C. 1984a, ApJ, 281, 381 [NASA ADS] [CrossRef] (In the text)
• Low, B. C. 1984b, ApJ, 281, 392 [NASA ADS] [CrossRef]
• Mastrodemos, N., & Morris, M. 1999, ApJ, 523, 357 [NASA ADS] [CrossRef] (In the text)
• Matt, S., Balick, B., Winglee, R., & Goodson, A. 2000, ApJ, 545, 965 [NASA ADS] [CrossRef] (In the text)
• Mellema, G. 1996, Ap&SS, 245, 239 [NASA ADS] [CrossRef]
• Mellema, G., & Frank, A. 1997, MNRAS, 292, 795 [NASA ADS] (In the text)
• Mellema, G., Eulderink, F., & Icke, V. 1991, A&A, 252, 718 [NASA ADS] (In the text)
• Miranda, L. F., & Solf, J. 1992, A&A, 260, 397 [NASA ADS] (In the text)
• Morris, M. 1981, ApJ, 249, 572 [NASA ADS] [CrossRef] (In the text)
• Osherovich, V. A., Farrugia, C. J., & Burlaga, L. F. 1993, J. Geophys. Res., 98, 13225 [NASA ADS] [CrossRef] (In the text)
• Osherovich, V. A., Farrugia, C. J., & Burlaga, L. F. 1995, J. Geophys. Res., 100, 12307 [NASA ADS] [CrossRef] (In the text)
• Pascoli, G. 1993, J. Astrophys. Astr., 14, 65 [NASA ADS] [CrossRef] (In the text)
• Pelletier, G., & Pudritz, R. E. 1992, ApJ, 394, 117 [NASA ADS] [CrossRef] (In the text)
• Soker, N., & Livio, M. 1994, ApJ, 421, 219 [NASA ADS] [CrossRef] (In the text)
• Tsui, K. H. 2006, Phys. Plasmas, 13, 072102 [NASA ADS] [CrossRef] (In the text)
• Tsui, K. H. 2008, Astron. Astrophys., 482, 793 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
• Tsui, K. H., & Serbeto, A. 2007, ApJ, 658, 794 [NASA ADS] [CrossRef] (In the text)
• Tsui, K. H., & Tavares, M. D. 2005, J. Atmos. Solar-Terr. Phys., 67, 1691 [NASA ADS] [CrossRef] (In the text)
• Tsui, K. H., Navia, C. E., Robba, M. B., Carneiro, L. T., & Emelin, S. E. 2006, Phys. Plasmas, 13, 113503 [NASA ADS] [CrossRef]
• Yoshida, Z., & Mahajan, S. M. 2002, Phys. Rev. Lett., 88, 095001 [NASA ADS] [CrossRef] (In the text)
• Zhang, C. Y., & Kwok, S. 1998, ApJ&SS, 117, 341 [NASA ADS] (In the text)
• Zhu, S. P., Horiuchi, R., Sato, T., & Complexity Simulation Group, 1995, Phys. Rev. E, 51, 6047 [NASA ADS] [CrossRef] (In the text)