A combined theory for magnetohydrodynamic equilibria with anisotropic pressure and magnetic shear

We present a new approach to the theory of magnetohydrodynamic equilibria with anisotropic pressure, magnetic shear and translational/rotational invariance. This approach involves combining two existing formalisms in order to eliminate their individual weaknesses. The theoretical aspects of the method are explored in detail along with numerical solutions which make use of the method. Eventually, this method could be applied to model various plasma systems, such as planetary magnetospheres.


Introduction
Magnetohydrodynamic (MHD) equilibria can be used to model a wide variety of physical processes, with applications varying from astrophysical systems such as magnetospheres and the sun, to laboratory plasmas such as those confined within tokamaks and stellarators. Not only are they useful as a starting point for dynamical studies, but the slow evolution of certain plasma systems can be modelled by sequences of equilibria-the so-called quasi-static approximation (see, for example, [1]).
One can often assume a scalar pressure is present in the force balance equation, however occasionally one must include a more general anisotropic pressure tensor. For instance, in particularly dense plasmas where particles have a small mean free path, the presence of collisions provides a mechanism to keep the plasma close to isotropy. However, in a diffuse or collisionless plasma (or indeed a plasma in the presence of a strong magnetic field), the pressure may have different components along and across the magnetic field and thus the inclusion of an anisotropic pressure tensor may be necessary. Alternatively, an anisotropic pressure tensor could arise from a driving process such as the solar wind. As for an astrophysical example, it has been suggested that pressure anisotropies might play a role in models of planetary magnetospheres, such as Jupiter (e.g. [2,3]) and the Earth (e.g. [4][5][6][7][8]). Anisotropic pressure also occurs frequently in laboratory plasma systems, where various heating mechanisms in tokamaks and stellarators can produce large pressure anisotropies [9]: formulations of these systems have been discussed by [10] and [11] amongst others. Thus it is clear that the physical systems where pressure anisotropy plays a major role are both vast in number and varied in application. This underpins the great importance of being able to understand and describe such a system in mathematically consistent way-indeed that is the main aim of this paper.
In the formulations referenced above one often assumes a gyrotropic pressure tensor, which effectively splits the pressure into two components, parallel and perpendicular to the magnetic field. The parallel pressure can then be described as a function of both the magnetic flux function and the magnetic field strength. In symmetric equilibria (i.e. equilibria with translational or rotational invariance), the shear field component of the magnetic field can then be derived from the parallel pressure and some arbitrary function of the magnetic flux function. This leads to an expression for the shear field in terms of the flux function and the magnetic field strength (e.g. [12,13]). We call this the total field formalism.
This has a couple of significant drawbacks. Firstly, one introduces a coupling problem where an equation for the shear field is intrinsically linked to an equation for the total magn etic field, and vice versa. This prevents any progress since one needs to know the shear before finding the total magnetic field, yet we need the total magnetic field in order to find the shear. Secondly, we have a fundamental problem when we prescribe the shear field in terms of the magnetic field strength. One finds that it is possible that the shear field can become larger than the magnetic field strength, which would lead to a contradiction in the definition of the shear field. In fact, there is no method in the total field formalism which guarantees that this problem is avoided a priori.
In a previous paper [14], the authors discussed the benefits of an approach which considers quantities as functions of the magnetic flux function and the poloidal magnetic field strength-we call this the poloidal formalism. In this formalism one specifies a quantity called the effective parallel pressure, which is some combination of the parallel pressure, its derivatives and the shear field. It is then possible to solve the resulting Grad-Shafranov equation to find the magnetic flux function. This approach leads to a consistent mathematical formulation of the equilibrium problem (specifically the shear field can never become larger than the total magnetic field strength), whilst also avoiding the difficulties of the implicit coupling problem and poorly defined shear field mentioned above.
However, there remain some compelling reasons to use the total field formalism instead of the poloidal formalism. Firstly, it makes more sense from a physical perspective to specify the parallel pressure rather than the effective pressure. Moreover, when using this method with rotational invariance, the scale factors involved mean that one must first satisfy a non-linear partial differential equation before specifying the pressures. Therefore, there is a good motivation to consider switching between these formalisms, which we call the combined approach: we can ensure a sensible definition of the pressure in terms of the total magnetic field strength, whilst avoiding the implicit coupling problem, shear field contradictions and the difficulties involving scale factors. For these reasons, we suggest that the combined formalism should be used as the new standard method of computing anisotropic pressure MHD equilibria. Using the combined approach, any physical MHD equilibrium system with anisotropic pressure can be formulated with a sensible definition of the parallel pressure, without ever running into problems with self-consistency.
The combined approach leads to quite a curious situation-a single-valued function in terms of the magnetic field strength can become a multi-valued function in terms of the poloidal magnetic field strength. This property can be traced back to regions where the shear field becomes larger than the total magnetic field strength in the total field formalism. Since these regions cannot exist in the poloidal formalism, they must manifest in some other formnamely the appearance of multiple branches. In this paper, we deal with the mathematical problems that arise from the appearance of multiple branches in quantities such as the parallel pressure and shear field. We will also discuss the idea of branch switching, where quantities can change characteristics by moving between these different branches. This can only happen if certain conditions are satisfied, which we will discuss in some detail. These restrictions also allow us to make deductions about the global magnetic field structure of a domain from knowledge of only a few key points. We would like to emphasise that the branches, branch switching and bifurcations we talk about in this paper are different to those occurring when discussing sequences of equilibria. Our branch switching may occur solely for a fixed set of parameters in a given equilibrium problem.
The structure of this paper is as follows. In section 2 we will outline the basic theory of the problem in Cartesian coordinates. Section 3 will outline the equivalent theory in spherical coordinates with rotational invariance, highlighting the implications of the extra scale factors. Section 4 will explain the drawbacks of the poloidal formalism and the need for a combination of both approaches in order to make progress. This section will also highlight the technical details of the mathematical approach in some generality, along with an example of how to implement the necessary techniques. Then, in section 5 we will discuss branch switching, and the constraints that must be satisfied within a hypothetical domain that allow this to occur. In this section we also discuss the implications of observations at local points to the global field structure, and provide an simplified example of an equilibrium which involves branch switching. Following this, in section 6, we present some first numerical results. We finish with some concluding remarks in section 7.

Basic theory with translational symmetry
The basic equations for an MHD equilibrium with anisotropic pressure are as follows. Firstly, the magnetic field B must satisfy the solenoidal constraint Secondly, we have an expression for the current density j from a form of Ampère's law: where µ 0 is the magnetic permeability of free space. Lastly we have a force balance equation, which we take to be where B is the magnetic field, j is the current density and P is the anisotropic pressure tensor. We assume P is gyrotropic and thus has the form 2 where I is the unit dyadic and =| | B B is the magnitude of the magnetic field. In Cartesian geometry with translational invariance in the y-direction we can express B in terms of a magnetic flux function A and shear field B y . Specifically, we have that In the isotropic case, one finds that the pressure and shear field must both be functions of the magnetic flux function. This is not the case for anisotropic pressure equilibria, however if one assumes that these quantities are functions of the magnetic flux function and the magnetic field strength, (1) reduces to the following equations which we will call the total field formalism: The first of these is a Grad-Shafranov type equation for the magnetic flux function. The second is a constraint equation which ties the shear field profile to the choice of parallel pressure. The function F on the right hand side of (3b) is some arbitrary function of A (see [13]). One can interpret F as the shear field profile one would have specified in the isotropic case. This can be seen intuitively by substituting an isotropic pressure ( / ∥ ∂ ∂ = P B 0) into (3b), where one sees that in the isotropic limit B y = F. Lastly, (3c) gives a relation which ties together the two components of pressure.
The problem with the formalism above is that (3a) and (3b) are implicitly coupled. In order to solve (3a), one must know B y before being able to determine B. However, in order to determine B y from (3b) we must know B. We call this the implicit coupling problem, and it poses some difficulty when trying to solve (3a) and (3b) in their current form. There are also situations where the shear field becomes larger than the total magnetic field strength, which is an obvious contradiction. This stems from (3b), which implicitly defines the shear field as a function of the total magnetic field strength via the parallel pressure. As an example, take a parallel pressure of the form / ∥ µ = P BB 0 0 and let F = B 0 , for some typical background field strength B 0 . Then (3b) gives the shear field as /( ) = − B BB B B y 0 0 . Thus, in the region 0 < B < 2B 0 , the shear field is larger than the total magnetic field strength. This problem is also present in the parallel pressure function ( ) /( ) ∥ = + P AB B exp 1, shown in figure 1 for F(A) = A 2 , which is frequently used since it corresponds to a bi-Maxwellian distribution function. Any region coloured red exhibits the unphysical property of | |> B B y . Whilst using the total field formalism, it is impossible to ensure beforehand that, after solving (3a), one does not end up in the red region of figure 1 at some point in the solution domain.
An alternative formalism is therefore required. In [14], the equilibrium problem is tackled under the assumption that the pressures are not functions of A and B (as in the total field formalism), but are functions of A and B p , where =|∇ | B A p is the poloidal magnetic field strength. One can also introduce a quantity ∥ P , the effective parallel pressure, which simplifies the resulting equations significantly. We list these equations for reference below, which we shall refer to as the poloidal formalism: This poloidal formulation is reminiscent of the 2-D problem: if the shear field is equal to zero in (3a), one arrives at (4a) with ∥ P and B in place of ∥ P and B p . The implicit coupling problem has also been avoided. If one specifies the effective parallel pressure, the shear field plays no role in (4a), and thus can be solved entirely independently of the other equations. The shear field can then be immediately found in terms of B p from (4b). We should also note that since B y is not defined in terms of B, there is no possibility of the shear becoming larger than the total magnetic field strength. Thus the poloidal formalism avoids the contradictions that occur in the total field formalism.

Basic theory with rotational symmetry
We now discuss the total field and poloidal formalisms with respect to equilibria which have rotational symmetry. Using spherical geometry and assuming invariance in the φ-direction we can express B in terms of a magnetic flux function A and shear field φ B , i.e.
Again, making the assumption that the pressures are functions of the magnetic flux function and the magnetic field strength, (1) reduces to the following equations: where, These equations are similar to the Cartesian case, however we now notice the presence of scale factors. This formalism also suffers from the same problems intrinsic to the Cartesian case: the implicit coupling problem is still present and (5b) still allows the possibility that the shear becomes larger than the total magnetic field strength.
We proceed in the same way as in the Cartesian case, but we must account for the scale factors inherent to spherical geometry. In [14] it is noted that if the pressure is written as a function of only two variables (A and B p ), then a simplification with an effective pressure is not possible. However, we can make a similar simplification if we assume that the pressures are a function of three variables: the magnetic flux function (A), the poloidal magnetic field strength ( , and the scale factor R. This leads to the following set of equations: Since we have assumed the parallel pressure is a function of three variables instead of two, we must have some uniqueness constraint. The appropriate constraint, which is derived in appendix, takes the form of a non-linear partial differential equation. In order for (6a)-(6e) to be valid, the effective parallel pressure must satisfy,

A combined approach
For each coordinate system we now have two formalisms, one based on the total magnetic field strength and one based on the poloidal magnetic field strength. For clarity, we will refer to these as the total field and poloidal formalisms respectively. The poloidal formalism eliminates both major flaws of the total field formalism and allows progress both analytically and numerically. However, on closer inspection, there are some compelling reasons why one may prefer to use the total field formalism. Firstly, it makes sense from a physical perspective to specify the parallel pressure rather than the effective pressure when solving the equilibrium problem. Since the effective pressure is a combination of the actual parallel pressure, the magnetic shear and the anisotropy, it is rather difficult to decide whether a given effective pressure is descriptive of a real physical system. One can, of course, work backwards from an effective pressure to recover the parallel pressure and magnetic shear, but one does not know from the outset what these quantities will be. This can prove to be a fairly substantial problem from a mathematical modelling point of view. For these reasons, it would be better to specify the parallel pressure, as is done in the total field formalism, in order to ensure that the equilibrium problem has physical relevance.
There is also the benefit that one can ensure the positivity of the pressures a priori in the total field formalism. When considering pressures as a function of A and B, one normally specifies the parallel pressure and then derives the perpendicular pressure from (3c). One can easily show that in order for both pressures to remain positive for all values of the total magnetic field strength then the parallel pressure must be given by In the poloidal formalism, however, it is not so easy to guarantee the positivity of both pressures. One can specify a positive parallel pressure without too much difficulty (at least in Cartesian coordinates), but when it comes to deriving the perpendicular pressure from (4d) or (6d), it is not clear what sort of contribution will come from the M y or φ M factors. That is to say, one can check whether both pressures are positive after specifying ∥ P , but there is not a formula a priori which will guarantee that both pressures are positive. This can lead to an awkward process of guessing forms of ∥ P , finding that the pressures could become negative, and then repeating until one finds a valid pressure by luck rather than judgement.
We must also not forget the uniqueness constraint in the problem with rotational invariance using the poloidal formalism. In the Cartesian case, we could simply specify an effective parallel pressure as a function of A and B p , and then solve (4a)-(4e) without fear of contradiction. In the rotationally invariant case, one cannot simply specify ∥ P as a function of the three variables A, B p and R, since there is no guarantee that this will satisfy (7). There are some analytic solutions to (7) which the authors have found by assuming certain forms for ∥ P , but these are generally unphysical and become negative or imaginary in certain parts of the R-z plane. A numerical approach is therefore required to solve the uniqueness constraint. Whilst it is possible to proceed in this way, this method lacks simplicity. Moreover, often one requires accurate derivatives of ∥ P in order to solve (6a) numerically, thus not having an analytic expression for ∥ P is quite a drawback. Thus we find ourselves in a situation where both formalisms have subtle difficulties associated with them. The key point, however, is that all of the problems of the total field formalism are solved by the poloidal formalism, and vice versa. We therefore propose a method which combines both formalisms. One should start by specifying a parallel pressure in terms of A and B, as in the total field formalism. In this way, the pressures can be guaranteed to be positive and (if we have rotational invariance) the uniqueness constraint (7) is automatically satisfied. This also means that we do not have to concern ourselves with specifying an effective parallel pressure that may not have physical relevance. Once the parallel pressure has been specified, we should then convert it into the poloidal formalism. This is a reasonably straightforward process which will make use of either (3b) or (5b), depending on the symmetry. Now that the parallel pressure is given in terms of the poloidal magnetic field strength, we can continue to solve the problem via (4a)-(4e) or (6a)-(6e).
We now arrive at the main point of this paper. When we convert between the formalisms as outlined above, one finds that a typical parallel pressure in the total field formalism can lead to a parallel pressure which is multi-valued in the poloidal formalism. Before we consider the implications of the appearance of multi-valued quantities, we should first discuss how and why they arise in the first place. Essentially, this is due to the invalid regions that exist in the total field formalism (the red regions in figure 1, for example). Since these regions cannot exist in the poloidal formalism, something must happen when we convert between the two systems. Specifically, the invalid regions 'roll up' and are forced into the complex plane In the process of pushing the invalid regions into the complex plane, distortions arise, hence the appearance of the multi-valued functions. Thus, the multi-valued functions are not simply a mathematical curiosity, but arise from the contradictory nature of the total field formalism itself.
We now describe the appearance of the multi-valued functions more rigorously. Here, it is easier to consider what happens to the shear field when we switch formalisms, and so we will focus on that quantity for now. Similarly, for simplicity, we will refer to the Cartesian casebut we should note that the following analysis also holds in the rotationally invariant case if we fix the radial coordinate R.
Recall that in the total field formalism, since the parallel pressure is given in terms of A and B, the shear field can also be written in terms of those variables, from an inversion of (3b). Let us fix A and consider lines of constant B y and B in the B y -B plane (see figure 2(a)). Note how constant B y (horizontal) lines start on the diagonal B y = B, otherwise we would have | |> B B y . Now consider how these lines transform when we switch into the poloidal formalism. If B y = k for some constant k in the B y -B plane, then we still have B y = k in the B y -B p plane, but now these horizontal lines start at It is the transformation of the vertical lines to circles that is responsible for the appearance of multi-valued quantities. A single valued curve can be bent back on itself by this transformation. This is shown in figure 2, where the black line is a typical single valued curve (in this ), which transforms into a semi-circle offset from the origin when we move into the B y -B p plane.
We now understand why these multi-valued quantities arise-the next logical question is to ask when they arise. To answer this we must consider the reverse grid transformation, shown in figure 3. This shows how a square grid in the B y -B p plane is transformed when we move back to the B y -B plane. Just as before, the horizontal lines remain unchanged. The vertical lines however, transform into hyperbolic curves, with Therefore, if one has a shear field given as a function of A and B, it is possible to check the quantity S for any sign changes. These indicate that the shear field must bend back on itself in B y -B p space. This covers branches that are connected to each other. It is also possible that entirely disconnected branches can arise. These will occur when the shear becomes greater than the total magnetic field strength, and then returns below that threshold. In the B y -B p plane this will correspond to a branch disappearing and reappearing at = B 0 p with different values of B y . This is shown graphically in figure 4. We can deduce precisely when this will happen by solving (3b) for B y and setting up the inequality | |> B B y . It is now useful to introduce the concept of the critical points of a parallel pressure. These critical points are values of the magnetic flux function where the general shape of the shear field (as a function of B p ) changes structure. This could be where a branch disappears, for instance, or a branch straightens out to become single valued, or even two distinct branches joining together. The latter two possibilities are shown in figures 5 and 6.

Example
As an example we take the parallel pressure function given in [5] which corresponds to an underlying bi-Maxwellian distribution function, If we choose F = A 2 for the free function in (3b) (an arbitrary choice made for illustrative purposes), then the regions in which the shear field becomes larger than the total magnetic field strength are given by the inequality We should also consider the quantity S, defined by (8), which we find to be From (10) and (11), we can deduce the critical points of (9) (recall that the critical points describe significant changes in the structure of the shear field). We find that there exist two critical points. The first is the real solution of the equation , which we denote c 1 = −0.7034..., and the second occurs at c 2 = 0.  For A < c 1 the region described by (10) spans from B = 0 to some positive value B = B c , after which all values of B do not satisfy (10). Thus, we should expect only a single branch of the shear field for A < c 1 . We then check S for A < c 1 and find that it does not change sign, hence the branch is only single-valued here. Indeed this is what we see when we graph the shear field in figure 7(a).
For A > c 1 , but ≠ A c 2 , the region described by (10) spans between two positive values of B. Thus, for these values of A we should expect precisely two branches. Again, we also check for sign changes in S. One finds that for these values of A prior to the region described by (10), there is a sign change in S. After the region described described by (10), there are no sign changes in S. We deduce that, for A > c 1 , but ≠ A c 2 , one of the branches bends back on itself, and the other branch does not. We see this behaviour in figures 7(b) and (c).
For the singular case where A = c 2 , the inequality (10) is never satisfied. In fact, when A = c 2 , we have that F = 0. This implies, from (3b), that either B y = 0 or For the latter we find that = = B B 0.465 57... 0 . Thus, we find that for A = c 2 , we either have B y = 0 or we lie on the circle radius B 0 centred on the origin in the B y -B p plane. This is shown in figure 7(d).
One should compare the plots in figure 7(a) with the total field plot of the parallel pressure in figure 1. The top single-valued branch in figure 7 corresponds to a parallel pressure in the green region of figure 1, and the lower branch corresponds to a parallel pressure in the blue region. That is why the lower branch only appears above the first critical point-the blue region in figure 1 vanishes for A < c 1 . We also note that this correspondence between the two formalisms means that the lower branches in figure 7 are always firehose unstable, in contrast to the upper branches which are stable to the firehose instability.

Branch switching
Now that we are aware of the possibility of multiple branches existing when B y is expressed as a function of A and B p , it makes sense to ask when can we switch between these branches. In general, the equilibrium must satisfy some rather specific conditions in order for branch switching to occur. The types of branch switching allowed will vary greatly depending on the precise choice of parallel pressure function and F. In this section we will restrict ourselves to continuous branch switching, however we note that it is possible for branch switching to occur across a discontinuity, although we leave that topic for further study in a future paper.
Suppose we have a function ( ) B A B , . It is also possible that for a certain value of A the bifurcation point will disappear, such as when A passes through a critical point. This could be due to either the branch portion straightening out, or the branch disappearing altogether. Therefore the bifurcation point will typically only be defined between two values of the magnetic flux function. These values can be deduced by examining the critical points of the parallel pressure as described above. We now introduce a parameter σ, where Again, this σ is only defined for the same A values as * B p . The contours of σ are necessary for determining when continuous branch switching can occur, and what specific types of switching are allowed. Crucially we should consider the contour σ = 1. Along this contour, the poloidal magnetic field equals * B p and potential branch switching can occur. The key point is that σ = 1 must be either a local minimum or a local maximum of σ in order to switch branches continuously. If σ was found to be higher on one side of the contour than the other, that would imply that we have passed a point where the shear field no longer exists (at least not on the branch which bends back upon itself). There is one exception to this, and that is when we cross σ = 1 at the same time as a critical value of A. When this happens it is possible that several branches join together and we can potentially switch onto one of these.
Indeed, we must also consider the contours where A is equal to one of the critical values for other reasons. These critical values determine a change in the general structure of the shear field, and crossing them implies certain constraints which must be taken into account. For instance, it is relatively common for certain branches to not exist past some critical values. In that case, for a smooth transition one requires that we stay on a branch that exists in both regions.
Considering the types of branch switching that are possible in a specific configuration allows us to make some deductions about the global field structure based on observations of a few specific points. For instance, the appearance of a null point will often signify that we are on a specific branch and thus the global structure of the field must reflect this. Alternatively, we might know the value of A in the entire domain, and thus, by analysing critical points, we can often deduce the types of shear field allowed without much in the way of further calculations.

A qualitative example
Consider the plot in figure 8, which shows a typical cut of the shear field given in the previous example in section 4.1. We can split B y into three branches, which we will call branches I-III corresponding to the solid, dashed and dot-dashed lines respectively. Note that branch III is entirely disconnected from branches I and II, which meet at the bifurcation point mentioned above (the black cross). Branch III exists for all values of A and ⩾ B 0 p , while branches I and II only exist below the bifurcation point and for A > c 1 .
Suppose we are on branch I (so A > A c ). Clearly, for a fixed value of A, the poloidal magnetic field must always be less than or equal to the value at the bifurcation point, otherwise the shear field would have to jump to branch III, creating a discontinuity. Also, if the poloidal magnetic field is always strictly less than the value at the bifurcation point, we are confined to the first branch since any move to a different branch would cause a jump in the shear field. The interesting case is, therefore, when the maximum value of the poloidal magnetic field is precisely equal to the value at the bifurcation point or, equivalently, when σ = 1 . It is only when this occurs that we can smoothly switch from branch I to branch II. We must also ensure that we never pass below A = A c , for at this point the bottom branches disappear. We show this diagrammatically in figure 9 for a hypothetical domain D in the x-z plane, where one must pass through the contour σ = 1 when it is a local maximum.
Another possible transition occurs when the contour of σ = 1 precisely overlaps with the contour F = 0. When this happens, across that contour the shear field is at the triple point of the pitchfork bifurcation in figure 7(d). Only at this point can we transition from branch I or II to branch III (or vice versa). Again, we show how this might occur in a hypothetical domain, D, in figure 10. The shear field is on branch I on the left hand side of the domain, and as one travels through to the right hand side of the domain we cross a contour of σ = 1 at the same time as a contour of A = A 0 = 0. After this point the shear field has switched onto branch III.
We can also discuss the implications of observations of specific structures inside a hypothetical domain. Suppose there exists a null point in the domain. One possibility is that the shear field is zero everywhere in the domain and we have a purely 2D solution. Otherwise, in order to preserve continuity, there must exist a region around the null point where we are on branch I, and thus the magnetic flux must be greater than the first critical point, A c . In this region we also know that the shear field must be negative (a property of branch I), and thus (since F is positive) the equilibrium must be firehose unstable. Suppose, in the same domain, we also measured a point where there was a positive value of the shear field. That point must be on branch III. Therefore, there must be a contour of σ = 1 and A = 0 which precisely overlap in order for us to have switched between branches (as in the example in figure 10). In another domain, we might observe that the magnetic flux function is less than the first critical value everywhere. In this case, we can automatically deduce that the shear field can only ever be positive, since we must be on branch III. Furthermore, since the third branch corresponds to values of B > B c , we know that there cannot exist any null points in this domain.

An explicit 1.5D example
It is relatively difficult to actually construct a full 2.5D example of the qualitative branch switching process that is described in the previous section. Indeed, if computing solutions numerically, one must keep track of all contours at all times, and be able to identify the appearance of local extrema in σ in order to switch onto different branches. However, it is simpler to consider a 1.5D example. In these cases we treat the magnetic flux as a function of x alone. This seems like a step back, since the entire point of the past few chapters has been to construct 2.5D solutions, however, considering 1.5D solutions is the first step in constructing examples which can actually exhibit branch switching.
In this example, consider a normalized effective parallel pressure as a dome in A-B p space, i.e.
where ∥ P has been normalized by /µ B 0 2 0 and B has been normalized by B 0 . This form of the effective parallel pressure was chosen in order to create a simplified model of the type of multi-valued effective parallel pressure that one may expect from the analysis in the previous sections. In particular, the dome structure is representative of the lower branches one recovers when the parallel pressure is given by (9), via a bi-Maxwellian distribution function (see also figures 7(c) and 8). We show a plot of ∥ P in figure 11, where one can see the individual branches as the red and blue surfaces. There are two critical points associated with ∥ P , which are the values = ± A 1. In the region between the critical points,the effective parallel pressure is as shown in figure 11, but beyond these critical points, both branches disappear. Thus the bifurcation point is only defined for −1 < A < 1, and is given by This allows us to find the quantity σ, which is defined for the same values of A as the bifurcation point, and is given by When σ = 1, we lie on the dome equator, and in order for branch switching to occur, we require that σ = 1 is a local maximum at some point in the domain.
It remains to solve the Grad-Shafranov equation for translationally invariant equilibria with a shear field component (4a), which upon substitution of the pressure (and after normalization) becomes To proceed, we assume the magnetic flux is a function of x alone. We then solve (15) using the computer software package Maple, which uses a Runge-Kutta-Fehlberg method to find A(x) after specifying the values of A(0) and ( ) ′ A 0 . This involves solving (15) by taking the positive sign of the pressure (i.e. on the upper portion of the dome). Eventually the solution will approach the dome equator (i.e. σ = 1), where we can either bounce back (remaining on the upper portion of the dome), or continue into the lower portion of the dome. In the second case, our solution will have exhibited branch switching.
In figure 12(a) we show plots of the magnetic flux function and magnetic field strength for a solution which remains on the upper portion of the dome, where no branch switching has occurred. The dotted line represents a point where we are on the dome equator, i.e. σ = 1 and B p is equal to the value at the bifurcation point * B p . However, figure 12(b) shows the solution which has exhibited branch switching. As the solution passes through the bifurcation point B p , we move onto the lower portion of the dome defined by ∥ P . We see that for a short distance after the switch from upper to lower branch, there is not much difference between the two cases. However, the behaviour of both the magnetic flux function and magnetic field strength changes drastically thereafter. The consequences on the value of the pressure are shown in figure 13. In the case without branch switching ( figure 13(a)), the pressure is bounded below by 1 since we are confined to the upper portion of the dome. In the case with branch switching ( figure 13(b)), we pass onto the lower portion of the dome, giving substantially lower pressures than the other case.
The effects are perhaps most dramatic on the shear field, which is shown for both cases in figure 14. When we have no branch switching ( figure 14(a)) the shear field is always positive remaining steady around 0.5, apart from dropping to zero briefly when we approach the possible switching location. In the case with branch switching ( figure 14(b)) we see a completely different shear field profile. The shear field changes sign after the switch point, but also increases by just over an order of magnitude.

Some numerical results
The analysis presented in section 4 is not only of theoretical use. In this section we show that one can actually apply the theory and obtain numerical solutions to the force balance equation with rotational symmetry and a shear field component of the magnetic field. We present a solution to the equilibrium problem (1) with a similar parallel pressure to that described in the example sections 4.1 and 5.1. In this paper we will focus only on one solution, rather than carrying out an extensive parameter study. This section is purely to show that it is now possible to find numerical solutions to the equilibrium problem by using the combined approach given in section 4.
Firstly we must normalize all of the quantities in (1). Since we are not focussing on modelling a particular physical system, we normalize all units by typical quantities: the magnetic field is normalized by a typical B 0 ; lengths by a typical lengthscale L 0 ; pressures by /µ B 0 2 0 ; and the magnetic flux by B L 0 2 . We then consider a normalized pressure of the form, The parameter k can be thought of as an anisotropy parameter (when k = 0 we have an isotropic pressure). As in previous examples, the anisotropic /( ) + B B k factor in the parallel pressure corresponds to an underlying bi-Maxwellian distribution function, and the isotropic A 2 dependence of the parallel pressure is a commonly used function when it comes to modelling isotropic planetary magnetospheres. The choice of a piecewise function introduces a cut-off field line, that is a field line after which the pressure is set to zero (in this case the cut off field line is A = 0.5). The cut-off field line is necessary to prevent the build up of current on the boundaries of the domain. We specify a similar cut-off for the free function F in (5b), which we choose to be This choice of F (which is arbitrary) leads to a shear field profile which is qualitatively the same as in figure 8: three branches, two of which join together and one of which is distinct. There exists one critical point, which is the value of the cut-off field line, below which branches I and II do not exist. We will make the assumption that we are always on branch III (the distinct branch), which ensures that there exists a single value of the shear field for all values of B p . Thus we avoid the difficulties associated with branch switching in 2.5D which were discussed in section 5. We then use the continuation code described above to increase the anisotropy parameter from k = 0 to k = 0.1. The solution is fixed at the boundaries (1 < r < 50 and θ π < < 0 ), where it takes on the value of a dipole field. Figure 15(a) shows values of the magnetic flux function in a region near the origin. It consists of a dipole field that has been stretched out along the equatorial plane. We also see that, in this example, we do not have to concern ourselves with branch switching, since it is impossible to switch from branch III to any other branch without passing through a local minimum at A = 0.5, a feature that is not present in figure 15(a). The poloidal magnetic field strength is shown in figure 15(b) with a logarithmic scale. In general the poloidal magnetic field strength is quite low, increasing towards the origin, but in the equatorial plane there is a region where the poloidal magnetic field strength drops to near zero. Now we consider the quantity φ b , shown in figure 16(a). We see that the shear field is strongest in the equatorial plane, peaking at around = φ b 0.25, and quickly drops to zero elsewhere in the domain. Now, recall that the function F is the value that the shear field would have taken in the isotropic case. Thus, if we plot the quantity / φ b F, (i.e. the shear field normalized by what the shear field would have been in the isotropic case), we can obtain a good picture of how much the shear field deviates from what we would expect in an isotropic regime. This is shown in figure 16(b). We see that in most parts of the domain, the shear field is equal to the free function F, however there is a region in the equatorial plane (1 < R < 4 and / θ π | | < 8) where we can see some large differences. At its peak in this region, the shear field is about three and a half times as large as it would have been in the isotropic case. This corresponds to the region of low poloidal magnetic field in figure 15(b). A common approximation in these sorts of models is to assume that the shear field is precisely equal to the function F (e.g. [13,31,32]), which would imply that / φ b F is unity everywhere in the domain. Thus, apart from demonstrating how the method can be used in general, this example also demonstrates that setting the shear field equal to F(A) may not be an approximation which is generally valid.

Conclusion
In a previous paper [14], it was demonstrated that defining quantities, such as the pressure and shear field, in terms of the magnetic field strength can lead to problems. Firstly, it is possible that the shear field can take on values larger than the total magnetic field strength, a clear contradiction. Secondly, the equations which arise are coupled in such a way that any progress, analytic or numeric, is not possible.
One can avoid these problems by using the poloidal formalism where such problems cannot arise. This formalism treats quantities as functions of the poloidal magnetic field strength, rather than the total magnetic field strength. The implicit coupling problem vanishes, and it is impossible for the shear field to become larger than the total magnetic field strength.
However, it soon becomes clear that there are compelling reasons not to define quantities in this manner. The effective pressure is a rather nebulous quantity to define, and it is difficult to ensure the positivity of pressures that are derived from it. Moreover, with rotational invariance, even specifying the parallel pressure as a function of the poloidal magnetic field strength can prove impossible, considering a non-linear partial differential equation must be satisfied. Therefore, in this paper, we propose a combined approach. One should specify quantities in terms of the total magnetic field strength, as in the total field formalism, and then convert into the poloidal formalism. In this way, we avoid the problems inherent to both formalisms. We can specify pressures in a meaningful way, assure their positivity and even avoid the uniqueness constraint with rotational invariance, whilst not needing to worry about the implicit coupling problem, or shear field components becoming too large.
At this point we encounter a curious situation, namely that single-valued quantities in the total field can become multi-valued when we convert them into the poloidal formalism. We have shown the conditions under which these multi-valued functions arise, and how to predict their appearance through analysis of various quantities related to the parallel pressure.
The appearance of multi-valued solutions then begs the question, under what conditions can we switch between branches? To answer this one must consider the quantity σ, introduced in section 5, and the crucial contour σ = 1. We must also consider the critical points of the parallel pressure-values of A where there is a restructuring of the shear field profile. We have also shown that by considering observations taken at a few points, one can deduce global properties of an equilibrium. In the example case shown, the mere appearance of a null point provides a large amount of information concerning the surrounding region. Indeed, we showed that when we observe a null point and a point of positive shear, there must be an overlap of contours A = c 2 and σ = 1-a property of our example that would have been difficult to predict otherwise. In section 5.2 we were able to construct a simplified 1.5D example equilibrium which exhibited branch switching. This example serves as an indication of the types of equilibria which may be found in future studies with branch switching in 2.5D.
We conclude with an indication of what work could be done in the future. The theory outlined in this paper can, for example, be used to conduct numerical simulations in spherical coordinates with rotational symmetry and magnetic shear. We show a brief example of this in section 6, where we note that a common simplification in the literature (namely that = φ b F) may not be valid in certain situations. A possible application for this could be fast rotating magnetospheres, where it has been suggested that pressure anisotropies may play a role [2]. Yet there are still some problems that need to be overcome. In our numerical calculations we were restricted to a single valued branch that existed for all values of the magnetic flux function. It would be interesting indeed to see some simulations of the theory discussed in section 5, where we switch between branches.

Appendix. A derivation of the poloidal formalism with rotational invariance
Starting from the rotationally invariant equations where quantities depend on A and B (i.e. (5a) and (5b)), we now try to resolve the problem by considering functions in terms of . A problem arises, however, where we cannot consider functions only in terms of A and B p (as we could in the translationally symmetric case) due to the presence of geometric scale factors. This is seen by considering B, the magnitude of the magnetic field. Recall, in spherical coordinates with axial symmetry, The presence of R in the definition of B means that when we change from considering quantities as functions of A and B to A and B p , we must not neglect the R dependence. We therefore consider quantities as functions of three variables: A, B p and R. The parallel pressure derivatives are related by the following equations: The R dependence of ( ) ∥ P A B R , , p is not entirely arbitrary because we have changed the parallel pressure from being a function of two variables to a function of three variables. In fact, the dependence is governed by the following equation which is found by combining (A.2) and (A.3), , , . We now use the poloidal formalism, where ∥ P is understood as a function of the three variables A, B p and R, to rewrite (5a) and (5b) as Note that (A.7) and (A.8) are equivalent formulations of (5b). We can combine (A.7) with the right hand side of (A.6) to arrive at the simplification In the translationally invariant case, it was possible to construct an effective pressure ( ∥ P ) which reduced the Grad-Shafranov equation with shear field to an equivalent Grad-Shafranov equation without shear field. It is possible to do this in the rotationally invariant case as well, however there are some difficulties when it comes to specifying ∥ P as we will see later. We define the effective parallel pressure as Note that the form of the effective pressure is given differently in (A.10) and (6c). These two equations are equivalent, and one can convert between the two by using (6b). Now consider the R derivatives of ∥ P , which will give a condition on the allowed form of the effective parallel pressure, namely