EDP Sciences
Free Access
Issue
A&A
Volume 514, May 2010
Article Number A38
Number of page(s) 11
Section The Sun
DOI https://doi.org/10.1051/0004-6361/200913723
Published online 07 May 2010
A&A 514, A38 (2010)

Three-dimensional solutions of the magnetohydrostatic equations: rigidly rotating magnetized coronae in cylindrical geometry

N. Al-Salti[*] - T. Neukirch - R. Ryan

School of Mathematics and Statistics, University of St. Andrews, St. Andrews KY16 9SS, UK

Received 23 November 2009 / Accepted 1 February 2010

Abstract
Context. Solutions of the magnetohydrostatic (MHS) equations are very important for modelling astrophysical plasmas, such as the coronae of magnetized stars. Realistic models should be three-dimensional, i.e., should not have any spatial symmetries, but finding three-dimensional solutions of the MHS equations is a formidable task.
Aims. We present a general theoretical framework for calculating three-dimensional MHS solutions outside massive rigidly rotating central bodies, together with example solutions. A possible future application is to model the closed field region of the coronae of fast-rotating stars.
Methods. As a first step, we present in this paper the theory and solutions for the case of a massive rigidly rotating magnetized cylinder, but the theory can easily be extended to other geometries, We assume that the solutions are stationary in the co-rotating frame of reference. To simplify the MHS equations, we use a special form for the current density, which leads to a single linear partial differential equation for a pseudo-potential U. The magnetic field can be derived from U by differentiation. The plasma density, pressure, and temperature are also part of the solution.
Results. We derive the fundamental equation for the pseudo-potential both in coordinate independent form and in cylindrical coordinates. We present numerical example solutions for the case of cylindrical coordinates.

Key words: magnetic fields - magnetohydrodynamics (MHD) - stars: magnetic field - stars: coronae - stars: activity

1 Introduction

Finding three-dimensional (3D) solutions to the magnetohydrostatic (MHS) equations, i.e. solutions with no spatial symmetry, is a formidable task. Only very few analytical solutions are known and even using numerical methods for calculating 3D MHS solutions is usually far from straightforward (e.g. Wiegelmann & Neukirch 2006; Wiegelmann et al. 2007).

Analytic solutions have been found for a number of different cases[*]. If external forces like the gravitational force or the centrifugal force can be neglected, force balance between the Lorentz force and the pressure gradient has to be achieved. A small number of analytical solutions in Cartesian and cylindrical coordinates have been found (e.g. Kaiser & Salat 1996; Woolley 1977; Kaiser et al. 1995; Woolley 1976; Kaiser & Salat 1997; Salat & Kaiser 1995; Shivamoggi 1986), and some of them have been generalized to include field-aligned incompressible flows (Petrie & Neukirch 1999).

The case where external forces cannot be neglected is often the most relevant for astrophysical applications. In particular, 3D solutions of the MHS equations in the presence of an external gravitational field have been found for this case, both in Cartesian (e.g. Petrie & Neukirch 2000; Neukirch 1997; Low 1992,1984; Neukirch & Rastätter 1999; Low 1982,1993a,b,1985) and in spherical coordinates (e.g. Bogdan & Low 1986; Osherovich 1985b; Neukirch 1995; Osherovich 1985a).

A systematic method for calculating a special class of 3D MHS equilibria for the case when external forces have to included, has been developed in a series of papers by Low (1992,1993a,b,2005,1991,1985) and Bogdan & Low (1986). The method is applicable to all external forces derived from a potential and assumes a special form for the electric current density to allow analytical progress. In the simplest possible case, the MHS equations reduce to a linear partial differential equation for the magnetic field. It has been shown that in Cartesian and spherical geometry, the fundamental equation is very similar to a Schrödinger equation (Neukirch & Rastätter 1999; Neukirch 1995). Therefore standard methods such as expansion in terms of orthogonal function systems (e.g. Rudenko 2001) or Green's functions (e.g. Petrie & Neukirch 2000) can be used for finding solutions, and this method has been used to model, for example, the solar corona (e.g. Zhao & Hoeksema 1994; Gibson et al. 1996; Zhao & Hoeksema 1993; Ruan et al. 2008; Gibson & Bagenal 1995; Zhao et al. 2000) and stellar coronae (e.g. Lanza 2008).

While the method has been mainly used to find 3D MHS solutions for the case of an external gravitational potential, Low (1991) has also developed the method for rigidly rotating systems subject to centrifugal forces. For those cases, the solutions are only stationary in the frame of reference that is rotating with the same angular velocity as the system itself. Recently, Neukirch (2009) has presented a couple of 3D MHS solutions for rigidly rotating magnetospheres in cylindrical geometry, again using the simplest case leading to a linear differential equation for the magnetic field.

In the present contribution we extend the theory to the case where both gravitational and centrifugal force are taken into account. This case is, for example, relevant to the coronal structure of fast-rotating stars (e.g. Jardine 2004; ud-Doula et al. 2006; Jardine & van Ballegooijen 2005; Townsend & Owocki 2005; Jardine & Unruh 1999; Ryan et al. 2005; Townsend et al. 2005), in particular to the closed field line region. Often, however, potential magnetic fields are used for models derived from stellar surface data (e.g. Morin et al. 2008; Jardine et al. 1999,2002,2001; Donati et al. 2008,2006), but there is observational evidence of some measured surface magnetic fields being non-potential (e.g. Hussain et al. 2002). Recently, Mackay & van Ballegooijen (2006) and Yeates et al. (2008) have developed a numerical technique to produce sequences of quasi-static non-linear force-free equilibria from time series of observed magnetograms. While this technique has so far only been applied to the Sun, it could in principle also be applied to other stars if magnetic field data with a sufficiently high time cadence are obtained. The theory presented in this paper could improve the potential field models and, in its simplest form, is not computationally more demanding than potential field models.

As in Neukirch (2009), we will present the theory in a general form following Low (1991), but then investigate the somewhat artificial, but illustrative case of a massive rigidly rotating central cylinder. This is done merely for mathematical convenience as it is much easier to impose boundary conditions on a cylindrical boundary. A full solution of the problem would also include a stellar wind on open field line regions and the need to solve a free boundary problem to determine the transition from open to closed field regions. A solution to this problem is beyond the scope of the present paper and we neglect flows altogether. Instead, for the solutions presented in this paper, we impose boundary conditions on an imaginary outer boundary, similar to the source surface used for potential field models. For this case, we determine solutions using standard numerical methods. In future work, one could as a first step towards solving the full problem, assume that the open field line regions are potential and thus try to determine the boundary between open and closed field line regions.

The paper is organized as follows. In Sect. 2 we present a brief derivation of the underlying theory, followed by illustrative example solutions in Sect. 3. We conclude the paper with a summary and discussion in Sect. 4.

2 Theory

2.1 Coordinate-independent theory

Before moving on to the special case of a massive rigidly rotating cylinder, we briefly outline the basic theory in a coordinate-independent form. In this way, the equations derived below are applicable to other cases as well, such as massive rigidly rotating spheres or ellipsoids (stars), or even synchronously rotating double stars, using e.g. the Roche potential.

We basically follow Low (1991) in our outline and refer the reader to his paper for more details (see also Neukirch 2009). The MHS equations in the co-rotating frame of reference are given by (see e.g. Mestel 1999)

                         $\displaystyle {\vec j} \times {\vec B} -\nabla p -\rho \nabla V = 0,$ (1)
    $\displaystyle \nabla \times {\vec B} = \mu_0 {\vec j},$ (2)
    $\displaystyle \nabla \cdot {\vec B} = 0,$ (3)

where ${\vec B}$ is the magnetic field, ${\vec j}$ is the current density, p is the pressure, $\rho$ is the plasma density and V is the combined centrifugal and gravitational potential. Assuming

\begin{displaymath}\mu_0{\vec j} = \nabla F \times \nabla V,
\end{displaymath} (4)

with F a free function one finds from the force balance Eq. (1) that

\begin{displaymath}p(\varpi,\phi,z) = p(F,V),
\end{displaymath} (5)

and

    $\displaystyle \left(\frac{\partial p}{\partial F}\right)_V = -\frac{1}{\mu_0}({\vec B} \cdot \nabla V ),$ (6)
    $\displaystyle \rho = - \left(\frac{\partial p}{\partial V}\right)_F + \frac{1}{\mu_0} ({\vec B} \cdot \nabla F).$ (7)

Further progress can be made by making an appropriate choice for the free function F. Choosing

\begin{displaymath}F(\varpi, \phi,z) = \kappa(V){\vec B}\cdot\nabla V,
\end{displaymath} (8)

with $\kappa(V)$ a free function, as suggested by Low (1991), leads to a linear relation between magnetic field and current density. In this case, we have the following expression

\begin{displaymath}p = p_0(V) - \frac{1}{2 \mu_0} \kappa(V) ({\vec B}\cdot\nabla V)^2
\end{displaymath} (9)

for the pressure. Here, p0(V) is an arbitrary function which represents a hydrostatic background atmosphere. For the density we find

\begin{displaymath}\rho =-\frac{{\rm d} p_0}{{\rm d}V} + \frac{1}{2\mu_0}\; \fra...
...\mu_0}\;\kappa(V) {\vec B}\cdot\nabla ({\vec B}\cdot\nabla V).
\end{displaymath} (10)

An expression for the plasma temperature can be obtained if we assume that the plasma satisfies the equation of state of an ideal gas,

\begin{displaymath}T= \frac{\mu p}{R \rho},
\end{displaymath} (11)

where R is the universal gas constant and $\mu$ is the mean molecular weight.

By integrating Ampère's law (2) one finds the magnetic field to be

\begin{displaymath}{\vec B} = \nabla U +
\frac{\kappa(V) }{1- \kappa(V) (\nabla V)^2} (\nabla U\cdot\nabla V )\nabla V,
\end{displaymath} (12)

where the free function U appears due to the integration. The pseudo-potential U is determined by substituting (12) into $\nabla\cdot {\vec B} =0$:

\begin{displaymath}\nabla \cdot \left(\nabla U + \frac{\kappa(V) }{1- \kappa(V) (\nabla V)^2} (\nabla U\cdot\nabla V )\nabla V \right) =0.
\end{displaymath} (13)

Equation (13) is a single partial differential equation for the pseudo-potential U and is the fundamental equation for the linear case of the theory presented here. An alternative form of this equation is

\begin{displaymath}\nabla \cdot ({{\bf M}}\cdot\nabla U) = 0,
\end{displaymath} (14)

with the 3 $\times$ 3 matrix ${{\bf M}}$ defined as

\begin{displaymath}{{\bf M}} = {{\bf I}} +\frac{\kappa(V) }{1- \kappa(V) (\nabla V)^2} \nabla V~ \nabla V.
\end{displaymath} (15)

Here ${{\bf I}}$ is the 3 $\times$ 3 unit matrix. Equation (14) is particularly useful if $\nabla V$has more than one non-vanishing component. This would, for example, be the case for the combined gravitational and centrifugal potential outside a massive rigidly rotating sphere. In such a case, Eq. (13) is usually not separable.

2.2 Cylindrical geometry

To illustrate how the theory presented above can be used, we treat in the present paper the somewhat artificial, but mathematically simpler case of a cylinder of radius R, infinite length and uniform mass per unit length M, rotating rigidly with angular velocity $\Omega$ about its symmetry axis. We use a co-rotating cylindrical coordinate system $\varpi$, $\phi$, z with the z-axis aligned with the rotation axis. The external gravitational potential (normalized to 0 at $\varpi=R$) of such a cylinder is given by

\begin{displaymath}\Psi = 2 G M \ln (\varpi/R),
\end{displaymath} (16)

and the combined potential V by

\begin{displaymath}V = - \frac{\Omega^2}{2} \varpi^2 + 2 \mbox{G} M \ln (\varpi/R).
\end{displaymath} (17)

Using Eq. (12), the components of ${\vec B}$ in cylindrical coordinates are

                           $\displaystyle B_\varpi = \frac{1}{1-\kappa(V) (V^\prime)^2} \frac{\partial U}{\partial \varpi},$ (18)
    $\displaystyle B_\phi = \frac{1}{\varpi}\frac{\partial U}{\partial \phi},$ (19)
    $\displaystyle B_z = \frac{\partial U}{\partial z},$ (20)

with

\begin{displaymath}V^\prime = \frac{{\rm d} V}{{\rm d}\varpi}\cdot
\end{displaymath} (21)

Defining

\begin{displaymath}\xi(\varpi) = \kappa(V) (V^\prime)^2,
\end{displaymath} (22)

one can rewrite Eq. (13) as

\begin{displaymath}\frac{1}{\varpi}\frac{\partial}{\partial \varpi}
\left( \frac...
...2 U}{\partial \phi^2} +
\frac{\partial^2 U}{\partial z^2} = 0,
\end{displaymath} (23)

which is the fundamental equation to be solved.

The pressure and density are given by

\begin{displaymath}p = p_0(V) - \frac{1}{2\mu_0} \kappa(V) {V^\prime}^2 B_\varpi^2
\end{displaymath} (24)

and
                               $\displaystyle \rho$ = $\displaystyle - \frac{{\rm d} p_0}{{\rm d} V} + \frac{1}{2\mu_0} \frac{{\rm d} ...
...^\prime}^2 B_\varpi^2 +
\frac{1}{\mu_0}\kappa(V) \; V^{\prime\prime} B_\varpi^2$  
    $\displaystyle + ~\frac{1}{\mu_0}\kappa(V)\; V^\prime\; {\vec B}\cdot\nabla B_\varpi.$ (25)

\begin{figure}
\par\includegraphics[width=9cm,clip]{13723fig1.eps}
\end{figure} Figure 1:

The combined potential $V(\varpi )$ for a co-rotation radius $\varpi _{\rm co}=4.0$.

Open with DEXTER

While formally, Eqs. (18) to (25) are identical to the case $M\to 0$ (only centrifugal forces) investigated by Neukirch (2009), an important difference is that the combined potential (17) does not have a one-to-one mapping to the radial coordinate $\varpi$ as the gravitational potential or the centrifugal potential on their own have. Instead, $V(\varpi )$ has a maximum (see Fig. 1) at the co-rotation radius given by

\begin{displaymath}\varpi_{co} = \frac{\sqrt{2GM}}{\vert\Omega\vert}\cdot
\end{displaymath} (26)

A test particle in a circular and planar orbit around the cylinder would have an orbital angular velocity which is equal to $\Omega$ so that it would co-rotate with the cylinder. More importantly, though, for a rigidly rotating plasma on the cylindrical surface with radius equalling the co-rotation radius, the outward centrifugal force is exactly balancing the inward gravitational force, as the expression

\begin{eqnarray*}-\rho\nabla V = -\rho V^\prime {\vec e}_\varpi
\end{eqnarray*}


for the combination of the two forces shows. Since $V^\prime$ vanishes at  $\varpi_{\rm co}$, the combined force is zero for $\varpi=\varpi_{\rm co}$. For distances from the cylinder larger than the co-rotation radius the centrifugal force will be bigger than the gravitational force ( $V^\prime <0)$and thus the combined force will be pointing outward. Obviously, overall force balance will have to include the Lorentz force and pressure gradient. The Lorentz force is crucial to be able to obtain force balance beyond the co-rotation radius.

In Neukirch (2009) the expression $\kappa(V) {V^\prime}^2$ was generally replaced by a function  $\xi(\varpi)$. Due to the one-to-one mapping between the centrifugal potential and the radial variable $\varpi$, it was possible to choose the function  $\xi(\varpi)$ instead of $\kappa(V)$. This is not generally possible for the combined potential discussed here. Although defining a function $\xi(\varpi) = \kappa(V) {V^\prime}^2$ is of course possible, choosing  $\xi(\varpi)$ instead of $\kappa(V)$ will generally lead to problems, for example to possible singularities of $\kappa(V)$ at the co-rotation radius, because

\begin{displaymath}\kappa(V(\varpi)) = \frac{\xi(\varpi)}{{V^\prime}^2(\varpi)}\cdot
\end{displaymath} (27)

Obviously the denominator vanishes at $\varpi_{\rm co}$ and $\kappa(V)$ will only be non-singular if  $\xi(\varpi)$ goes to zero with the same or a higher power of  $\varpi-\varpi_{\rm co}$ at the co-rotation radius. This excludes any simple choice such as $\xi(\varpi)=\xi_0$ = const., which was one of the examples used in Neukirch (2009).

Singularities in $\kappa(V)$ in turn lead to singularities of density and temperature as well. If we express the density in terms of  $\xi(\varpi)$ instead of $\kappa(V)$ we first obtain

                                 $\displaystyle \rho$ = $\displaystyle - \frac{ \mbox{d} p_0}{\mbox{d} V} +\frac{1}{2\mu_0} \frac{\mbox{...
...^\prime}^2 B_\varpi^2 + \frac{1}{\mu_0}\kappa(V) \; V^{\prime\prime} B_\varpi^2$  
    $\displaystyle + ~\frac{1}{\mu_0}\kappa(V)\; V^\prime\; {\vec B}\cdot\nabla B_\varpi.$ (28)

Using

\begin{eqnarray*}\frac{{\rm d} \xi}{{\rm d}\varpi} & = & \frac{{\rm d}}{{\rm d}\...
...x{d} V} {V^\prime}^3 + 2 \kappa(V) V^\prime \; V^{\prime\prime},
\end{eqnarray*}


we can rewrite Eq. (28) in the form

\begin{displaymath}\rho = - \frac{{\rm d} p_0}{{\rm d} V} +
\frac{1}{V^\prime}\...
...ime}\frac{1}{\mu_0} \xi(\varpi)({\vec B}\cdot\nabla B_\varpi),
\end{displaymath} (29)

which makes the possible singularity at $\varpi=\varpi_{\rm co}$ ( $V^\prime = 0$) obvious. The pressure is always non-singular since

\begin{displaymath}p = p_0(V) - \frac{1}{2\mu_0} \kappa(V) {V^\prime}^2 B_\varpi^2
= p_0(V) - \frac{1}{2\mu_0} \xi(\varpi) B_\varpi^2.
\end{displaymath} (30)

But even if a singularity of $\kappa(V)$ and $\rho$ could be avoided for a suitable choice of  $\xi(\varpi)$ (going through 0 quadratically at  $\varpi_{\rm co}$), the inverse mapping from $\varpi$ to V would not be well-defined across  $\varpi_{\rm co}$, and therefore a given function  $\xi(\varpi)/{V^\prime}^2(\varpi)$ cannot generally be expressed as a function of V. This has to borne in mind when we considering Eq. (23) in which  $\xi(\varpi)$ should merely be regarded as an abbreviation for  $\kappa(V) {V^\prime}^2$, but not as an independent free function as, for example, in Neukirch (2009).

As already stated above, Eq. (23) is the fundamental equation that has to be solved. It is straightforward to see that for

\begin{displaymath}1-\kappa(V) {V^\prime}^2 = 1 -\xi(\varpi) \gtrless 0,
\end{displaymath} (31)

Eq. (23) is either elliptic or hyperbolic, respectively, while having singularities at any radius  $\varpi_{\rm s}$ where $\kappa(V) {V^\prime}^2=1$. Obviously, none of the singularities coincides with the co-rotation radius (at  $\varpi=\varpi_{\rm co}$ we have $1-\kappa(V) {V^\prime}^2 = 1$). The singularities, i.e. the transition from an elliptic to a hyperbolic equation or vice versa can only occur for $\kappa(V)$ positive. For example, assuming for simplicity that $\kappa>0$ is constant, the critical points occur at

\begin{displaymath}\varpi _{s_{\pm }}^{2}=\frac{\varpi _{\rm co}^{2}}{2}\left[ \...
...\varpi _{\rm co}^{2}}{(2GM)^2\kappa
}+2\right) ^{2}-4}\right].
\end{displaymath} (32)

The critical point defined by $\varpi _{\rm s_{+}}^{2}$ is beyond the co-rotation radius, whereas the one given by  $\varpi _{\rm s_{-}}^{2}$ is within the co-rotation radius (see Fig. 2). The inner singularity lies inside the central cylinder, if

\begin{displaymath}(2GM)^2 \kappa < \frac{\varpi_{\rm co}^4R^2}{(\varpi_{\rm co}^2-R^2)^2},
\end{displaymath}

where of course $R^2<\varpi_{\rm co}^2$. One should note, however, that Eq. (32) only applies if $\kappa$ = const. If $\kappa$ depends on V (and thus on $\varpi$), even the possibility of more than two critical points exists in principle.

\begin{figure}
\par\includegraphics[width=9cm,clip]{13723fig2.eps}
\vspace*{5mm}
\end{figure} Figure 2:

The critical points $\varpi _{\rm s_{+}}^{2}$ (blue) and $\varpi _{\rm s_{-}}^{2}$ (red) for a co-rotation radius $\varpi _{\rm co}=2.0$.

Open with DEXTER

The case $\kappa(V)$ positive generally corresponds to a stretching of magnetic field compared to a potential field ($\kappa=0$), as can be seen from Eq. (18). A thought experiment where one starts with $\kappa=0$ (potential field) and then slowly increases $\kappa$ shows that the radial component of the magnetic field will increase due to the decrease of $1-\kappa {V^\prime}^2$, if one assumes that to lowest order U does not change too rapidly with changing $\kappa$. Furthermore it is relatively straightforward to see that the radial component of the Lorentz force will be directed inwards for magnetic fields which have the same general behaviour as a dipole field close to the equator (z=0). This is exactly what is expected of stretched magnetic fields acting to confine plasma pulled away from the cylinder by the centrifugal force.

In this paper we will only consider Eq. (23) for cases where it is elliptic. We shall follow a common approach used in solar and stellar applications and in addition to the inner cylindrical boundary define an artificial outer boundary. This is similar to the source surface used in many global potential magnetic field models of the solar corona (the so-called potential field source surface or PFSS models). It should be noted, however, that because the magnetic fields calculated in the present paper are non-potential, we impose slightly different boundary conditions from those usually imposed on a source surface when potential fields are used.

3 Solution methods and example solutions

In this section we discuss possible solution methods and a few illustrative example solutions. We first nondimensionalize all quantities and equations using $\widetilde{\mathbf{\nabla }}=R\mathbf{\nabla}$, $\widetilde{\varpi}=\frac{\varpi}{R}$, $\widetilde{{\vec B}}=\frac{{\vec B}}{B_{0}}$, $\widetilde{V}=\frac{V}{2 G M}$, $\widetilde{{\vec j}}=\frac{{\vec j}}{(B_{0}/\mu_{0} R)}$, $\widetilde{p}=\frac{p}{B_{0}^{2}/\mu_{0}}$, and $ \widetilde{\rho}=\frac{\rho}{B_{0}^{2}/(2 \mu_{0} G M)}$, where B0 is a typical magnetic field value. The dimensionless combined centrifugal and gravitational potential is given by

\begin{displaymath}V = - \frac{1}{2 \varpi_{\rm co}^2} (\varpi^2 -\varpi_{\rm co}^2 \ln (\varpi^2)),
\end{displaymath} (33)

where, $\varpi_{\rm co}= \frac{\sqrt{2GM}/\vert\Omega\vert}{R}$ is the dimensionless co-rotation radius and the cylinder radius in these dimensionless coordinates is 1.

\begin{figure}
\par\includegraphics[width=18.35cm,clip]{13723fig3.eps}
\vspace*{6mm}
\end{figure} Figure 3:

Field line plots for the example solution. The left panel shows a side view, the right panel a view along the z-axis. The colours on the boundary represent the radial magnetic field component, $B_\varpi $ on that boundary.

Open with DEXTER

3.1 Separation of variables

In the case when it is elliptic, Eq. (23) is very similar to Laplace's equation and admits separable solutions (see also Neukirch 2009) of the form

\begin{displaymath}U(\varpi,\phi,z) = F_{mk}(\varpi) \exp(\mbox{i} m \phi) \exp(\mbox{i} k z).
\end{displaymath} (34)

If we substitute (34) into (23) we find that the radial function  $F_{mk}(\varpi)$ satisfies the equation

\begin{displaymath}\frac{1}{\varpi}\frac{{\rm d}}{{\rm d} \varpi}
\left( \frac{\...
...\right) -
\left( \frac{m^2}{\varpi^2} + k^2\right) F_{mk} = 0.
\end{displaymath} (35)

This ordinary second order differential equation will have two linearly independent solutions,  $F^{(1)}_{mk}(\varpi)$ and  $F^{(2)}_{mk}(\varpi)$, say. Since the partial differential equation for U is linear, the solutions for different m and k may be superposed to generate other solutions, and the most general form of a solution of (23) is

                                            $\displaystyle U(\varpi,\phi,z) = \sum_{m=-\infty}^\infty \exp(\mbox{i} m \phi)$  
    $\displaystyle \qquad\quad \times ~\int_{-\infty}^\infty {\rm d}k [A_m(k) F_{mk}^{(1)}(\varpi) +
B_m(k) F_{mk}^{(2)}(\varpi)] \exp(\mbox{i} k z).$ (36)

Here the Am(k) and Bm(k) are complex coefficients, which are determined by the boundary conditions, e.g. Dirichlet or von Neumann conditions in the elliptic case.

As already discussed above, we are not allowed to choose the function  $\xi(\varpi)$in the present case, if the domain includes the co-rotation radius. However, to illustrate the method and for use as a test case for the numerical method used later, we show a few plots for an analytical solution which can be obtained for $\xi(\varpi)=\xi_0$ = const. (see e.g. Neukirch 2009). In this case the general solutions to Eq. (35) are given by

\begin{displaymath}F_{mk}(\varpi) = A_{m}(k) \; I_\nu(k\sqrt{1-\xi_0} \;\varpi) +
B_{m}(k) \; K_\nu(k\sqrt{1-\xi_0} \;\varpi),
\end{displaymath} (37)

where $I_\nu(x)$ and $K_\nu(x)$ are modified Bessel functions (Abramowitz & Stegun 1965), $\nu = m \sqrt{1-\xi_0}$. Am(k) and Bm(k) are constants which would usually be determined by the boundary conditions.

For this illustrative example we have chosen the parameter values $\xi_0=3/4$, m=2 and $k=\pi/5$. This choice of parameters leads to $\sqrt{1-\xi_0} = 1/2$ and $\nu = m \sqrt{1-\xi_0} = 1$, and we set Am(k)=Bm(k)=0, except for  $B_{\pm 2}(k)$ which we choose so that the pseudo-potential is given by

\begin{displaymath}U=B_0 {K}_1(\pi \;\varpi/10)\sin(2\phi) \sin(\pi \; z/5).
\end{displaymath} (38)

The magnetic field components are then given by

                                      $\displaystyle B_\varpi = -\frac{2\pi}{5} B_0 [{K}_0(\pi \;\varpi/10) + \frac{10}{\pi \; \varpi} {K}_1(\pi \;\varpi/10) ]$  
    $\displaystyle \qquad\;~ \times ~\sin(2\phi) \sin(\pi \; z/5),$  
    $\displaystyle B_\phi = \frac{2\;B_0}{\varpi} {K}_1(\pi \;\varpi/10)\cos(2\phi)\sin(\pi \; z/5),$  
    $\displaystyle B_z = \frac{\pi\; B_0}{5} {K}_1(\pi \;\varpi/10)\sin(2\phi)\cos(\pi \; z/5).$ (39)

The reason we choose $K_\nu$ instead of $I_\nu$ is that $K_\nu$ decreases with increasing argument, which means that the magnetic field strength decreases with increasing distance from the z-axis. In Fig. 3 we show a 3D plot of magnetic field lines from two different viewing angles. The boundary colours represent the radial magnetic field component, $B_\varpi $. The non-symmetric nature of the magnetic field is obvious from the plot.

\begin{figure}
\par\includegraphics[width=18cm,clip]{13723fig4a.eps}\par\includegraphics[width=18cm,clip]{13723fig4b.eps}
\end{figure} Figure 4:

Cross section plots of the pressure ( top panels) and density ( bottom panels) variation in the xz-plane at y=0.5, y= 1, and y= 2.5, respectively.

Open with DEXTER

\begin{figure}
\par\includegraphics[width=18cm,clip]{13723fig5a.eps}\par\includegraphics[width=18cm,clip]{13723fig5b.eps}
\end{figure} Figure 5:

Examples of pressure ( top panels) and density isosurfaces ( bottom panels).

Open with DEXTER

\begin{figure}
\par\includegraphics[width=18cm,clip]{13723fig6.eps}
\vspace*{5mm}
\end{figure} Figure 6:

3D file line plot ( left panel) and view along the z-axis ( right) for the numerical solution of the $\xi (\varpi )=\xi _0=3/4$ case. The similarity of the plots with Fig. 3 is obvious.

Open with DEXTER

The pressure is given by

\begin{displaymath}p=p_0(V) - \frac{\xi_0}{2} B_\varpi^2 = p_0(V) - \frac{3}{4} B_\varpi^2,
\end{displaymath} (40)

which as discussed above is non-singular at the co-rotation radius. The density, however, is given by

\begin{displaymath}\rho = -\frac{{\rm d} p_0}{{\rm d} V} + \frac{\xi_0}{\; V^\pr...
...{4(\varpi_{\rm co}^2-\varpi^2)} {\vec B}\cdot \nabla B_\varpi,
\end{displaymath} (41)

and here the singularity at the co-rotation radius $\varpi_{\rm co}$ is obvious. We therefore consider this solution only for $\varpi < \varpi_{\rm co}$. It should be noted that when  $\xi(\varpi)$ is chosen directly, the value of the co-rotation radius affects the solution only through the presence of $V^\prime$ in the density. As the co-rotation radius is the only parameter in which the angular velocity $\Omega$ appears, choosing  $\xi(\varpi)$ instead of $\kappa(V)$ basically eliminates the rotation rate from the problem. Again this is a feature of the solutions which is not necessarily wanted if one wants to study the effect of increasing $\Omega$ on the solutions. Plots of the pressure and density contours and isosurfaces are shown in Figs. 4 and 5. Note that in the plots we only show the deviations from the cylindrically symmetrical background pressure and density.

3.1.1 Numerical solutions of equation (23)

In general, finding analytical solutions of Eqs. (23) or (35) will be impossible even for simple choices of the function $\kappa(V)$. Thus, numerical methods will have to be used to find solutions. Since Eq. (23) is a simple linear partial differential equation, standard numerical methods can be used to solve it. In the present paper we used an adaptive mesh finite element method from the COMSOL Multiphysics 3.4 package with MATLAB to solve Eq. (23).

\begin{figure}
\par\includegraphics[width=9cm,clip]{13723fig7a.eps}\par\includeg...
...}\par\includegraphics[width=9cm,clip]{13723fig7c.eps}
\vspace*{3mm}
\end{figure} Figure 7:

Magnetic field lines plots for the three cases of a) aligned rotator, b) oblique rotator and  c) displaced dipole.

Open with DEXTER

To check the accuracy of our numerical method, we have first solved Eq. (23) for the constant $\xi$ case presented in Sect. 3.1, using the same parameter values, as well as boundary conditions that are consistent with the analytical solution. We solve Eq. (23) for U on a numerical domain, which is bounded by an inner cylinder of radius 1, an outer cylinder of radius $\varpi_{\rm o}=6$, and which extends from -5 to 5 in the z-direction. The outer boundary  $\varpi_{\rm o}$ is assumed to be smaller than the co-rotation radius in this case to avoid singularities in the density.

The exact boundary conditions used for this case are

  • $B_\varpi (1,\phi,z) = B_{\varpi, {\rm analytical}}(1,\phi,z)$ on the inner boundary, where $B_{\varpi, {\rm analytical}}(\varpi,\phi,z)$ is the expression given in (39);

  • $U(6,\phi,z)= U_{\rm analytical}(6,\phi,z)$ on the outer boundary, where $U_{\rm analytical}(\varpi,\phi,z)$ is given by Eq. (38) and;

  • ${{\vec n}\cdot {\vec B}}= 0$ at $z=\pm 5$.
The mesh size used for this calculation consists of 261 996 elements.

A magnetic field line plot for the numerical solution obtained is shown in Fig. 6, which shows good agreement with the analytical solution. The only noticeable difference is the structure of the field lines towards $z=\pm 5$ which is due to the effect of the boundary condition at $z=\pm 5$ for the numerical solution.

Having thus convinced ourselves that the numerical tool gives satisfactory results, we have considered the simplest possible choice of $\kappa(V)$, which is

\begin{displaymath}\kappa(V) = \kappa_0 = {\rm const.}
\end{displaymath} (42)

as an example for a case where $\kappa$ is chosen directly. We have calculated numerical solutions to Eq. (23) using as boundary conditions on the surface of the central cylinder ($\varpi=1$) a magnetic dipole field ( ${\vec B}_{\rm dip}$) for the three cases of the
(a)
magnetic dipole moment at the origin and aligned with the rotation axis (aligned rotator);

(b)
magnetic dipole moment at the origin, but inclined with respect to the rotation axis (oblique rotator) and;

(c)
magnetic dipole moment not located at the origin and inclined with respect to the rotation axis (displaced dipole).
These three cases are similar to the cases discussed by Neukirch (2009).

We use the same domain and mesh size as in the previous example, with the outer boundary conditions given by $U(6,\phi,z)=U_{\rm o}(6,\phi,z)$, where $U(\varpi,\phi,z)$ satisfies $\mathbf{\nabla}U_{\rm o}={\vec B}_{\rm dip}$ and ${{\vec n}\cdot {\vec B}}= 0$ at $z=\pm 5$.

By choosing $\kappa(V)$, the density remains non-singular at the co-rotation radius. Hence, we can now calculate solutions in a domain including and extending beyond the co-rotation radius.

Numerical solutions for the case of $\kappa=\frac{3}{4}$ and $\varpi_{\rm co}=3.5$ are illustrated in Figs. 7-10, where the letters (a)-(c) above each plot indicate the three different boundary conditions mentioned above. For the oblique rotator case the magnetic dipole moment is in the xz-plane at an angle of  $\frac{\pi}{4}$ with the x-axis. For the displaced dipole case the magnetic moment is again in the xz-plane, but now at x=0.3, making an angle of  $\frac{\pi}{4}$ with the x-axis. In the plots showing pressure and density, we only show the 3D deviation from background pressure and density. It turns out that both the 3D pressure deviation and the 3D density deviation are negative, which means that these terms will reduce any background pressure and density to lower values.

Figure 7 shows magnetic field line plots, with the colour contours on the central cylinder indicating the strength of the radial magnetic field component, $B_\varpi $, for the three different boundary conditions. As is to be expected, the change of boundary conditions have a clear effect on the structure of the magnetic field, which is clearly symmetric for the aligned rotator case, but becomes non-symmetric for the other two cases.

This is also visible in Fig. 8, where we show isosurfaces of the 3D deviation of the pressure from the background pressure. One can see that one has smaller isosurfaces close to the inner boundary where the magnetic field (and thus the pressure deviation) is strong, whereas the isosurfaces become more extended as one moves away from the cylinder and the magnetic field becomes weaker.

It can be clearly seen that for the case of aligned rotator (top panels in Fig. 8), the pressure isosurfaces are symmetric, whereas this symmetry is broken for the other two cases. In particular, the symmetry with respect to the x-axis and the z-axis is broken, but for the oblique rotator case (middle panels in Fig. 8) a notion of symmetry about the dipole axis remains. The least symmetric case, at least in terms of pressure isosurfaces, is the displaced dipole case (lower panels in Fig. 8).

\begin{figure}
\par\includegraphics[width=18cm,clip]{13723fig8a.eps}\par\include...
...\par\includegraphics[width=18cm,clip]{13723fig8c.eps}
\vspace*{2mm}
\end{figure} Figure 8:

Pressure isosurface plots for the aligned rotator case ( top panels), the oblique rotator case ( middle panels) and the displaced dipole case ( bottom panels). The transition from rotational symmetric isosurfaces in the aligned rotator case to asymmetric isosurfaces in the other two case can be seen very clearly.

Open with DEXTER

\begin{figure}
\par\includegraphics[width=18cm,clip]{13723fig9a.eps}\par\include...
...\par\includegraphics[width=18cm,clip]{13723fig9c.eps}
\vspace*{5mm}
\end{figure} Figure 9:

Variation of the pressure deviation from the background pressure in the xz-plane at y=0 (through the central cylinders), y=2, and y=3.5 (touching the co-rotation cylinder), respectively. Shown is the logarithm of the pressure. The increasing asymmetry from top to bottom is obvious.

Open with DEXTER

\begin{figure}
\par\includegraphics[width=17.4cm,clip]{13723fig10a.eps}\par\incl...
...fig10b.eps}\par\includegraphics[width=17.4cm,clip]{13723fig10c.eps}
\end{figure} Figure 10:

Variation of the density deviation from the background pressure in the xz-plane at y=0 (through the central), y= 2, and y= 3.5 (touching the co-rotation cylinder), respectively. Shown is the logarithm of the density. The density deviation is largest close to the cylinder.

Open with DEXTER

Figures 9 and 10 show cross section plots of the variation of the 3D pressure and density deviations from the background pressure and density in planes parallel to the xz-plane for different y-values. These plots show that there is some symmetry of the pressure and density deviations about the y-axis for all three cases, but clearly show symmetry about the x and z-axes only for the aligned rotator case. The intersection between the planes shown in Fig. 9 and the cylindrical surface with radius equal to the co-rotation radius coincides with the dark vertical features in the plots. In the rightmost panels, the plane basically touches the co-rotation cylinder and thus one only sees a single broad vertical feature. It can be easily seen (e.g. from Eq. (24)) that the total pressure is equal to the background pressure at the co-rotation radius, i.e. p=p0(V) at $\varpi=\varpi_{\rm co}$. The dark vertical features in Fig. 9 thus correspond to a vanishing 3D pressure deviation, whereas no corresponding feature exists for the 3D density deviation (Fig. 10). The pressure and density cross-section plots confirm the increasing degree of asymmetry when going from the aligned rotator case over the oblique rotator case to the displace dipole moment case.

4 Summary and discussion

We have presented a relatively simple (semi-)analytical approach which allows the modelling of 3D rigidly rotating magnetized coronae or magnetospheres around massive central objects. In the present paper we have restricted our analysis for illustrative purposes to the simpler, but less realistic case of cylindrical geometry. The possibility of extending the theory to other geometries will be discussed below.

The theory contains free functions $\kappa(V)$ and p0(V), where, in the case presented in this paper, V is the combined gravitational and centrifugal potential in the co-rotating frame of reference. Whereas the function $\kappa(V)$ implicitly determines the current density in the corona, p0(V) is an independent background pressure. Alternatively the derivative  ${\rm d} p_0/{\rm d} V= -\rho_0(V)$ can be chosen, where $\rho_0(V)$is a background density. The background pressure can then be determined by integration, if an equation of state and/or a temperature profile is assumed.

The function $\kappa(V)$ appears in the theory in the combination  $\kappa(V) {V^\prime}^2$ and, in the cylindrical geometry used in the present paper, a new function $\xi$ of the radial coordinate, $\varpi$, can be defined as $\xi(\varpi) = \kappa(V) {V^\prime}^2$. As has been shown before (Neukirch 2009) for the case of V being just the centrifugal potential (no gravitational force), analytical solutions of the theory can in principle be found if  $\xi(\varpi)$is chosen to have a convenient form. However, for the case of a combined centrifugal and gravitational potential as presented in this paper, the direct choice of a function  $\xi(\varpi)$ instead of deriving it from a chosen function $\kappa(V)$ generally leads to singularities, in particular of the density, at the co-rotation radius.

One can avoid these problems by choosing $\kappa(V)$ instead of  $\xi(\varpi)$. In this case, however, the fundamental equation is usually too complicated to allow for analytical solutions to be found, but the equation is still sufficiently simple that standard numerical methods can be used to solve it. We have presented an example of an analytical solution to be able to test our numerical method, and the numerical solution shows good agreement with the analytical solution on its domain of validity inside the co-rotation radius. We have then presented numerical solutions for the case $\kappa(V) = \kappa_0= {\rm const.}$ for three different types of boundary conditions on the surface of the central cylinder: a magnetic dipole field generated by a dipole moment located at the origin, aligned with the rotation axis (aligned rotator), a magnetic dipole field generated by a dipole moment located at the origin, but at an angle with the rotation axis (oblique rotator) and a magnetic dipole moment displaced from the origin, with the dipole moment not aligned with the rotation axis. These three cases were used to illustrate the transition from a rotationally symmetric corona to an asymmetric corona for the simple geometry of a magnetic dipole field.

A similar theory can also be developed for rotating spherical massive bodies. The combined gravitational and centrifugal potential for a body of mass M0 whose rotation axis is aligned with the z-axis has the form (using spherical coordinates r$\theta$ and $\phi$)

\begin{displaymath}V(r,\theta) = -\frac{1}{2} \Omega^2 r^2 \sin^2 \theta - \frac{G M_0}{r}\cdot
\end{displaymath} (43)

Due to the dependence of V on two of the coordinates in this case, Eq. (13) has a much more complicated form since

    $\displaystyle B_r = \left[1+f\left(\frac{\partial V }{\partial r} \right)^2\rig...
...l V}{\partial \theta}\left(\frac{1}{r}\frac{\partial U}{\partial \theta}\right)$ (44)
    $\displaystyle B_\theta = \frac{f}{r}\frac{\partial V}{\partial r} \frac{\partia...
...\right)^2\right]\left(\frac{1}{r}\frac{\partial U}{\partial \theta}\right)\cdot$ (45)

where

\begin{displaymath}f = \frac{\kappa(V)}{1-\kappa(V) (\nabla V)^2}\cdot
\end{displaymath} (46)

Both Br and $B_\theta$ depend on $\partial U/ \partial r$ and $\partial U/ \partial \theta$ for this case and this leads to mixed second derivatives in Eq. (13). It is highly unlikely that the resulting Eq. (13) for this case has any analytical solutions, although this still has to be investigated in detail. However, despite its more complicated form, solving Eq. (13) for the spherical case using standard numerical methods as, for example, the ones we have used in this paper, is not generically more difficult than solving the cylindrical case presented above. The main reason for this is that, despite its seemingly more complicated form, the type of Eq. (13) is again determined completely by the sign of the expression  $1-\kappa(V) (\nabla V)^2$. If this term is positive Eq. (13) is elliptic, otherwise it is hyperbolic. This can be seen relatively easily by writing Eq. (13) in the form (14) with

\begin{displaymath}{{\bf M}} = \left(
\begin{array}{ccc}
1 + f \left(\frac{\part...
...partial \theta}\right)^2 & 0 \\
0 & 0 & 1
\end{array}\right).
\end{displaymath} (47)

The nature of Eq. (13) is determined by the signs of the eigenvalues of the real and symmetric matrix ${{\bf M}}$. If all eigenvalues have the same sign, Eq. (13) is elliptic, otherwise it is hyperbolic. A straightforward calculation shows that ${{\bf M}}$, as given in (47), has a double eigenvalue 1 and that the third eigenvalue is given by $1/[1-\kappa(V) (\nabla V)^2]$, which corroborates our statement from above. We can thus conclude that it should be possible to use standard numerical methods for linear elliptic second order partial differential equations to solve Eq. (13) for the spherical case. Preliminary results obtained for the spherical case with the same numerical methods used for the cylindrical case so far confirm this conclusion and it is planned that full account of the spherical case will be given in a future publication.

Acknowledgements
N.A. acknowledges financial support by Sultan Qaboos University, Oman, T.N. acknowledges financial support by the UK's Science and Technology Facilities Council and by the European Commission through the SOLAIRE Network (MTRN-CT-2006-035484).

References

Footnotes

... Al-Salti[*]
Permanent address: DOMAS, Sultan Qaboos University, PO Box 36 P.C. 123 Muscat, OMAN.
... cases[*]
We explicitly exclude force-free fields from this discussion.

All Figures

  \begin{figure}
\par\includegraphics[width=9cm,clip]{13723fig1.eps}
\end{figure} Figure 1:

The combined potential $V(\varpi )$ for a co-rotation radius $\varpi _{\rm co}=4.0$.

Open with DEXTER
In the text

  \begin{figure}
\par\includegraphics[width=9cm,clip]{13723fig2.eps}
\vspace*{5mm}
\end{figure} Figure 2:

The critical points $\varpi _{\rm s_{+}}^{2}$ (blue) and $\varpi _{\rm s_{-}}^{2}$ (red) for a co-rotation radius $\varpi _{\rm co}=2.0$.

Open with DEXTER
In the text

  \begin{figure}
\par\includegraphics[width=18.35cm,clip]{13723fig3.eps}
\vspace*{6mm}
\end{figure} Figure 3:

Field line plots for the example solution. The left panel shows a side view, the right panel a view along the z-axis. The colours on the boundary represent the radial magnetic field component, $B_\varpi $ on that boundary.

Open with DEXTER
In the text

  \begin{figure}
\par\includegraphics[width=18cm,clip]{13723fig4a.eps}\par\includegraphics[width=18cm,clip]{13723fig4b.eps}
\end{figure} Figure 4:

Cross section plots of the pressure ( top panels) and density ( bottom panels) variation in the xz-plane at y=0.5, y= 1, and y= 2.5, respectively.

Open with DEXTER
In the text

  \begin{figure}
\par\includegraphics[width=18cm,clip]{13723fig5a.eps}\par\includegraphics[width=18cm,clip]{13723fig5b.eps}
\end{figure} Figure 5:

Examples of pressure ( top panels) and density isosurfaces ( bottom panels).

Open with DEXTER
In the text

  \begin{figure}
\par\includegraphics[width=18cm,clip]{13723fig6.eps}
\vspace*{5mm}
\end{figure} Figure 6:

3D file line plot ( left panel) and view along the z-axis ( right) for the numerical solution of the $\xi (\varpi )=\xi _0=3/4$ case. The similarity of the plots with Fig. 3 is obvious.

Open with DEXTER
In the text

  \begin{figure}
\par\includegraphics[width=9cm,clip]{13723fig7a.eps}\par\includeg...
...}\par\includegraphics[width=9cm,clip]{13723fig7c.eps}
\vspace*{3mm}
\end{figure} Figure 7:

Magnetic field lines plots for the three cases of a) aligned rotator, b) oblique rotator and  c) displaced dipole.

Open with DEXTER
In the text

  \begin{figure}
\par\includegraphics[width=18cm,clip]{13723fig8a.eps}\par\include...
...\par\includegraphics[width=18cm,clip]{13723fig8c.eps}
\vspace*{2mm}
\end{figure} Figure 8:

Pressure isosurface plots for the aligned rotator case ( top panels), the oblique rotator case ( middle panels) and the displaced dipole case ( bottom panels). The transition from rotational symmetric isosurfaces in the aligned rotator case to asymmetric isosurfaces in the other two case can be seen very clearly.

Open with DEXTER
In the text

  \begin{figure}
\par\includegraphics[width=18cm,clip]{13723fig9a.eps}\par\include...
...\par\includegraphics[width=18cm,clip]{13723fig9c.eps}
\vspace*{5mm}
\end{figure} Figure 9:

Variation of the pressure deviation from the background pressure in the xz-plane at y=0 (through the central cylinders), y=2, and y=3.5 (touching the co-rotation cylinder), respectively. Shown is the logarithm of the pressure. The increasing asymmetry from top to bottom is obvious.

Open with DEXTER
In the text

  \begin{figure}
\par\includegraphics[width=17.4cm,clip]{13723fig10a.eps}\par\incl...
...fig10b.eps}\par\includegraphics[width=17.4cm,clip]{13723fig10c.eps}
\end{figure} Figure 10:

Variation of the density deviation from the background pressure in the xz-plane at y=0 (through the central), y= 2, and y= 3.5 (touching the co-rotation cylinder), respectively. Shown is the logarithm of the density. The density deviation is largest close to the cylinder.

Open with DEXTER
In the text


Copyright ESO 2010

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.