Threedimensional solutions of the magnetohydrostatic equations: rigidly rotating magnetized coronae in cylindrical geometry
N. AlSalti^{}  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
threedimensional, i.e., should not have any spatial symmetries,
but finding threedimensional solutions of the MHS equations is a
formidable task.
Aims. We present a general theoretical framework for calculating
threedimensional 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
fastrotating 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
corotating 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
pseudopotential 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
pseudopotential 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 threedimensional (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 fieldaligned 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 fastrotating stars (e.g. Jardine 2004; udDoula 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 nonpotential (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 quasistatic nonlinear forcefree 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 Coordinateindependent theory
Before moving on to the special case of a massive rigidly rotating cylinder, we briefly outline the basic theory in a coordinateindependent 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 corotating frame of reference are given by (see e.g. Mestel 1999)
where is the magnetic field, is the current density, p is the pressure, is the plasma density and V is the combined centrifugal and gravitational potential. Assuming
(4) 
with F a free function one finds from the force balance Eq. (1) that
(5) 
and
Further progress can be made by making an appropriate choice for the free function F. Choosing
with 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
for the pressure. Here, p_{0}(V) is an arbitrary function which represents a hydrostatic background atmosphere. For the density we find
An expression for the plasma temperature can be obtained if we assume that the plasma satisfies the equation of state of an ideal gas,
where R is the universal gas constant and is the mean molecular weight.
By integrating Ampère's law (2) one finds the magnetic field to be
where the free function U appears due to the integration. The pseudopotential U is determined by substituting (12) into :
Equation (13) is a single partial differential equation for the pseudopotential U and is the fundamental equation for the linear case of the theory presented here. An alternative form of this equation is
with the 3 3 matrix defined as
(15) 
Here is the 3 3 unit matrix. Equation (14) is particularly useful if has more than one nonvanishing 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
about its symmetry axis. We use a corotating cylindrical coordinate system ,
,
z with the zaxis aligned with the rotation axis. The external gravitational potential (normalized to 0 at )
of such a cylinder is given by
(16) 
and the combined potential V by
Using Eq. (12), the components of in cylindrical coordinates are
with
(21) 
Defining
one can rewrite Eq. (13) as
which is the fundamental equation to be solved.
The pressure and density are given by
and
Figure 1: The combined potential for a corotation radius . 

Open with DEXTER 
While formally, Eqs. (18) to (25) are identical to the case
(only centrifugal forces) investigated by Neukirch (2009), an important difference is that the combined potential (17) does not have a onetoone mapping to the radial coordinate
as the gravitational potential or the centrifugal potential on their own have. Instead, has a maximum (see Fig. 1) at the corotation radius given by
(26) 
A test particle in a circular and planar orbit around the cylinder would have an orbital angular velocity which is equal to so that it would corotate with the cylinder. More importantly, though, for a rigidly rotating plasma on the cylindrical surface with radius equalling the corotation radius, the outward centrifugal force is exactly balancing the inward gravitational force, as the expression
for the combination of the two forces shows. Since vanishes at , the combined force is zero for . For distances from the cylinder larger than the corotation radius the centrifugal force will be bigger than the gravitational force ( 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 corotation radius.
In Neukirch (2009) the expression
was generally replaced by a
function
.
Due to the onetoone mapping between the centrifugal potential and the radial variable ,
it was possible to choose the function
instead of .
This is not generally possible for the combined potential discussed here. Although defining a function
is of course possible, choosing
instead of
will generally lead to problems, for example to possible singularities of
at the corotation radius, because
(27) 
Obviously the denominator vanishes at and will only be nonsingular if goes to zero with the same or a higher power of at the corotation radius. This excludes any simple choice such as = const., which was one of the examples used in Neukirch (2009).
Singularities in
in turn lead to singularities of density and temperature as well. If we express the density in terms of
instead of
we first obtain
Using
we can rewrite Eq. (28) in the form
which makes the possible singularity at ( ) obvious. The pressure is always nonsingular since
But even if a singularity of and could be avoided for a suitable choice of (going through 0 quadratically at ), the inverse mapping from to V would not be welldefined across , and therefore a given function cannot generally be expressed as a function of V. This has to borne in mind when we considering Eq. (23) in which should merely be regarded as an abbreviation for , 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
(31) 
Eq. (23) is either elliptic or hyperbolic, respectively, while having singularities at any radius where . Obviously, none of the singularities coincides with the corotation radius (at we have ). The singularities, i.e. the transition from an elliptic to a hyperbolic equation or vice versa can only occur for positive. For example, assuming for simplicity that is constant, the critical points occur at
The critical point defined by is beyond the corotation radius, whereas the one given by is within the corotation radius (see Fig. 2). The inner singularity lies inside the central cylinder, if
where of course . One should note, however, that Eq. (32) only applies if = const. If depends on V (and thus on ), even the possibility of more than two critical points exists in principle.
Figure 2: The critical points (blue) and (red) for a corotation radius . 

Open with DEXTER 
The case positive generally corresponds to a stretching of magnetic field compared to a potential field (), as can be seen from Eq. (18). A thought experiment where one starts with (potential field) and then slowly increases shows that the radial component of the magnetic field will increase due to the decrease of , if one assumes that to lowest order U does not change too rapidly with changing . 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 socalled potential field source surface or PFSS models). It should be noted, however, that because the magnetic fields calculated in the present paper are nonpotential, 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
,
,
,
,
,
,
and
,
where B_{0} is a typical magnetic field value. The dimensionless combined centrifugal and gravitational potential is given by
where, is the dimensionless corotation radius and the cylinder radius in these dimensionless coordinates is 1.
Figure 3: Field line plots for the example solution. The left panel shows a side view, the right panel a view along the zaxis. The colours on the boundary represent the radial magnetic field component, 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
If we substitute (34) into (23) we find that the radial function satisfies the equation
This ordinary second order differential equation will have two linearly independent solutions, and , 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
(36) 
Here the A_{m}(k) and B_{m}(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
in
the present
case, if the domain includes the corotation 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
= const. (see e.g. Neukirch 2009). In this case the general solutions to Eq. (35) are given by
(37) 
where and are modified Bessel functions (Abramowitz & Stegun 1965), . A_{m}(k) and B_{m}(k) are constants which would usually be determined by the boundary conditions.
For this illustrative example we have chosen the parameter values ,
m=2 and .
This choice of parameters leads to
and
,
and we set
A_{m}(k)=B_{m}(k)=0, except for
which we choose so that the pseudopotential is given by
The magnetic field components are then given by
The reason we choose instead of is that decreases with increasing argument, which means that the magnetic field strength decreases with increasing distance from the zaxis. 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, . The nonsymmetric nature of the magnetic field is obvious from the plot.
Figure 4: Cross section plots of the pressure ( top panels) and density ( bottom panels) variation in the xzplane at y=0.5, y= 1, and y= 2.5, respectively. 

Open with DEXTER 
Figure 5: Examples of pressure ( top panels) and density isosurfaces ( bottom panels). 

Open with DEXTER 
Figure 6: 3D file line plot ( left panel) and view along the zaxis ( right) for the numerical solution of the case. The similarity of the plots with Fig. 3 is obvious. 

Open with DEXTER 
The pressure is given by
which as discussed above is nonsingular at the corotation radius. The density, however, is given by
and here the singularity at the corotation radius is obvious. We therefore consider this solution only for . It should be noted that when is chosen directly, the value of the corotation radius affects the solution only through the presence of in the density. As the corotation radius is the only parameter in which the angular velocity appears, choosing instead of 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 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 . 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).
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 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 , and which extends from 5 to 5 in the zdirection. The outer boundary is assumed to be smaller than the corotation radius in this case to avoid singularities in the density.
The exact boundary conditions used for this case are
 on the inner boundary, where is the expression given in (39);
 on the outer boundary, where is given by Eq. (38) and;
 at .
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 which is due to the effect of the boundary condition at for the numerical solution.
Having thus convinced ourselves that the numerical tool gives
satisfactory results, we have considered the simplest possible choice
of ,
which is
(42) 
as an example for a case where is chosen directly. We have calculated numerical solutions to Eq. (23) using as boundary conditions on the surface of the central cylinder () a magnetic dipole field ( ) 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).
We use the same domain and mesh size as in the previous example, with the outer boundary conditions given by , where satisfies and at .
By choosing , the density remains nonsingular at the corotation radius. Hence, we can now calculate solutions in a domain including and extending beyond the corotation radius.
Numerical solutions for the case of and are illustrated in Figs. 710, 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 xzplane at an angle of with the xaxis. For the displaced dipole case the magnetic moment is again in the xzplane, but now at x=0.3, making an angle of with the xaxis. 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, , 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 nonsymmetric 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 xaxis and the zaxis 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).
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 
Figure 9: Variation of the pressure deviation from the background pressure in the xzplane at y=0 (through the central cylinders), y=2, and y=3.5 (touching the corotation cylinder), respectively. Shown is the logarithm of the pressure. The increasing asymmetry from top to bottom is obvious. 

Open with DEXTER 
Figure 10: Variation of the density deviation from the background pressure in the xzplane at y=0 (through the central), y= 2, and y= 3.5 (touching the corotation 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 xzplane for different yvalues. These plots show that there is some symmetry of the pressure and density deviations about the yaxis for all three cases, but clearly show symmetry about the x and zaxes only for the aligned rotator case. The intersection between the planes shown in Fig. 9 and the cylindrical surface with radius equal to the corotation radius coincides with the dark vertical features in the plots. In the rightmost panels, the plane basically touches the corotation 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 corotation radius, i.e. p=p_{0}(V) at . 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 crosssection 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 and p_{0}(V), where, in the case presented in this paper, V is the combined gravitational and centrifugal potential in the corotating frame of reference. Whereas the function implicitly determines the current density in the corona, p_{0}(V) is an independent background pressure. Alternatively the derivative can be chosen, where 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 appears in the theory in the combination and, in the cylindrical geometry used in the present paper, a new function of the radial coordinate, , can be defined as . 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 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 instead of deriving it from a chosen function generally leads to singularities, in particular of the density, at the corotation radius.
One can avoid these problems by choosing instead of . 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 corotation radius. We have then presented numerical solutions for the case 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 M_{0} whose rotation axis is aligned with the zaxis has the form (using spherical coordinates r,
and )
Due to the dependence of V on two of the coordinates in this case, Eq. (13) has a much more complicated form since
(44)  
(45) 
where
(46) 
Both B_{r} and depend on and 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 . 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
The nature of Eq. (13) is determined by the signs of the eigenvalues of the real and symmetric matrix . If all eigenvalues have the same sign, Eq. (13) is elliptic, otherwise it is hyperbolic. A straightforward calculation shows that , as given in (47), has a double eigenvalue 1 and that the third eigenvalue is given by , 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.
AcknowledgementsN.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 (MTRNCT2006035484).
References
 Abramowitz, M., & Stegun, I. A. 1965, Handbook of Mathematical Functions (New York: Dover Publications) [Google Scholar]
 Bogdan, T. J., & Low, B. C. 1986, ApJ, 306, 271 [NASA ADS] [CrossRef] [Google Scholar]
 Donati, J.F., Howarth, I. D., Jardine, M. M., et al. 2006, MNRAS, 370, 629 [NASA ADS] [CrossRef] [Google Scholar]
 Donati, J.F., Morin, J., Petit, P., et al. 2008, MNRAS, 390, 545 [NASA ADS] [CrossRef] [Google Scholar]
 Gibson, S. E., & Bagenal, F. 1995, J. Geophys. Res., 100, 19865 [NASA ADS] [CrossRef] [Google Scholar]
 Gibson, S. E., Bagenal, F., & Low, B. C. 1996, J. Geophys. Res., 101, 4813 [NASA ADS] [CrossRef] [Google Scholar]
 Hussain, G. A. J., van Ballegooijen, A. A., Jardine, M., & Collier Cameron, A. 2002, ApJ, 575, 1078 [NASA ADS] [CrossRef] [Google Scholar]
 Jardine, M. 2004, A&A, 414, L5 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Jardine, M., & Unruh, Y. C. 1999, A&A, 346, 883 [NASA ADS] [Google Scholar]
 Jardine, M., & van Ballegooijen, A. A. 2005, MNRAS, 361, 1173 [NASA ADS] [CrossRef] [Google Scholar]
 Jardine, M., Barnes, J. R., Donati, J.F., & Collier Cameron, A. 1999, MNRAS, 305, L35 [NASA ADS] [CrossRef] [Google Scholar]
 Jardine, M., Collier Cameron, A., Donati, J.F., & Pointer, G. R. 2001, MNRAS, 324, 201 [NASA ADS] [CrossRef] [Google Scholar]
 Jardine, M., Collier Cameron, A., & Donati, J.F. 2002, MNRAS, 333, 339 [NASA ADS] [CrossRef] [Google Scholar]
 Kaiser, R., & Salat, A. 1996, , 77, 3133 [Google Scholar]
 Kaiser, R., & Salat, A. 1997, J. Plasma Phys., 57, 425 [NASA ADS] [CrossRef] [Google Scholar]
 Kaiser, R., Salat, A., & Tataronis, J. A. 1995, Phys. Plasmas, 2, 1599 [NASA ADS] [CrossRef] [Google Scholar]
 Lanza, A. F. 2008, A&A, 487, 1163 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Low, B. C. 1982, ApJ, 263, 952 [NASA ADS] [CrossRef] [Google Scholar]
 Low, B. C. 1984, ApJ, 277, 415 [NASA ADS] [CrossRef] [Google Scholar]
 Low, B. C. 1985, ApJ, 293, 31 [NASA ADS] [CrossRef] [Google Scholar]
 Low, B. C. 1991, ApJ, 370, 427 [NASA ADS] [CrossRef] [Google Scholar]
 Low, B. C. 1992, ApJ, 399, 300 [NASA ADS] [CrossRef] [Google Scholar]
 Low, B. C. 1993a, ApJ, 408, 689 [NASA ADS] [CrossRef] [Google Scholar]
 Low, B. C. 1993b, ApJ, 408, 693 [NASA ADS] [CrossRef] [Google Scholar]
 Low, B. C. 2005, ApJ, 625, 451 [NASA ADS] [CrossRef] [Google Scholar]
 Mackay, D. H., & van Ballegooijen, A. A. 2006, ApJ, 641, 577 [NASA ADS] [CrossRef] [Google Scholar]
 Mestel, L. 1999, Stellar Magnetism (Oxford: Clarendon Press) [Google Scholar]
 Morin, J., Donati, J.F., Petit, P., et al. 2008, MNRAS, 390, 567 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Neukirch, T. 1995, A&A, 301, 628 [NASA ADS] [Google Scholar]
 Neukirch, T. 1997, A&A, 325, 847 [NASA ADS] [Google Scholar]
 Neukirch, T. 2009, Geophys. Astrophys. Fluid Dyn., 135, 535 [NASA ADS] [CrossRef] [Google Scholar]
 Neukirch, T., & Rastätter, L. 1999, A&A, 348, 1000 [NASA ADS] [Google Scholar]
 Osherovich, V. A. 1985a, Aust. J. Phys., 38, 975 [NASA ADS] [Google Scholar]
 Osherovich, V. A. 1985b, ApJ, 298, 235 [NASA ADS] [CrossRef] [Google Scholar]
 Petrie, G. J. D., & Neukirch, T. 1999, Geophys. Astrophys. Fluid Dyn., 91, 269 [NASA ADS] [CrossRef] [Google Scholar]
 Petrie, G. J. D., & Neukirch, T. 2000, A&A, 356, 735 [NASA ADS] [Google Scholar]
 Ruan, P., Wiegelmann, T., Inhester, B., et al. 2008, A&A, 481, 827 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Rudenko, G. V. 2001, Sol. Phys., 198, 279 [NASA ADS] [CrossRef] [Google Scholar]
 Ryan, R. D., Neukirch, T., & Jardine, M. 2005, A&A, 433, 323 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Salat, A., & Kaiser, R. 1995, Phys. Plasmas, 2, 3777 [NASA ADS] [CrossRef] [Google Scholar]
 Shivamoggi, B. K. 1986, Q. Appl. Math., 44, 487 [Google Scholar]
 Townsend, R. H. D., & Owocki, S. P. 2005, MNRAS, 357, 251 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Townsend, R. H. D., Owocki, S. P., & Groote, D. 2005, ApJ, 630, L81 [NASA ADS] [CrossRef] [Google Scholar]
 udDoula, A., Townsend, R. H. D., & Owocki, S. P. 2006, ApJ, 640, L191 [NASA ADS] [CrossRef] [Google Scholar]
 Wiegelmann, T., & Neukirch, T. 2006, A&A, 457, 1053 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Wiegelmann, T., Neukirch, T., Ruan, P., & Inhester, B. 2007, A&A, 475, 701 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Woolley, M. L. 1976, J. Plasma Phys., 15, 423 [NASA ADS] [CrossRef] [Google Scholar]
 Woolley, M. L. 1977, J. Plasma Phys., 17, 259 [NASA ADS] [CrossRef] [Google Scholar]
 Yeates, A. R., Mackay, D. H., & van Ballegooijen, A. A. 2008, Sol. Phys., 247, 103 [NASA ADS] [CrossRef] [Google Scholar]
 Zhao, X., & Hoeksema, J. T. 1993, Sol. Phys., 143, 41 [NASA ADS] [CrossRef] [Google Scholar]
 Zhao, X., & Hoeksema, J. T. 1994, Sol. Phys., 151, 91 [NASA ADS] [CrossRef] [Google Scholar]
 Zhao, X. P., Hoeksema, J. T., & Scherrer, P. H. 2000, ApJ, 538, 932 [NASA ADS] [CrossRef] [Google Scholar]
Footnotes
 ... AlSalti^{}
 Permanent address: DOMAS, Sultan Qaboos University, PO Box 36 P.C. 123 Muscat, OMAN.
 ... cases^{}
 We explicitly exclude forcefree fields from this discussion.
All Figures
Figure 1: The combined potential for a corotation radius . 

Open with DEXTER  
In the text 
Figure 2: The critical points (blue) and (red) for a corotation radius . 

Open with DEXTER  
In the text 
Figure 3: Field line plots for the example solution. The left panel shows a side view, the right panel a view along the zaxis. The colours on the boundary represent the radial magnetic field component, on that boundary. 

Open with DEXTER  
In the text 
Figure 4: Cross section plots of the pressure ( top panels) and density ( bottom panels) variation in the xzplane at y=0.5, y= 1, and y= 2.5, respectively. 

Open with DEXTER  
In the text 
Figure 5: Examples of pressure ( top panels) and density isosurfaces ( bottom panels). 

Open with DEXTER  
In the text 
Figure 6: 3D file line plot ( left panel) and view along the zaxis ( right) for the numerical solution of the case. The similarity of the plots with Fig. 3 is obvious. 

Open with DEXTER  
In the text 
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 
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 
Figure 9: Variation of the pressure deviation from the background pressure in the xzplane at y=0 (through the central cylinders), y=2, and y=3.5 (touching the corotation cylinder), respectively. Shown is the logarithm of the pressure. The increasing asymmetry from top to bottom is obvious. 

Open with DEXTER  
In the text 
Figure 10: Variation of the density deviation from the background pressure in the xzplane at y=0 (through the central), y= 2, and y= 3.5 (touching the corotation 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