1. Introduction
Twisted magnetospheres of neutron stars are being investigated for the role that the twist plays in the dynamics and radiation of these stars, particularly in the case of magnetars (Thompson, Lyutikov & Kulkarni Reference Thompson, Lyutikov and Kulkarni2002). By twist, it is meant that the magnetic field has a toroidal component that is supported by magnetospheric currents. Among other things, twists can affect the opacity of the magnetosphere (e.g. Viganò et al. Reference Viganò, Pons and Miralles2011) because of the strong currents and particle densities they require, or affect the spin-down rate of the star by modifying the polar cap structure (Glampedakis, Lander & Andersson Reference Glampedakis, Lander and Andersson2014). Untwisting is proposed to be at the origin of at least some of the stellar radiation (e.g. Beloborodov Reference Beloborodov2009; Viganò et al. Reference Viganò, Pons and Miralles2011). Twists represent an energy reservoir which can provide the energy dissipated into giant magnetar flares (Mahlmann et al. Reference Mahlmann, Akgün, Pons, Aloy and Cerdá-Durán2019). Recent work has been dedicated to studying the dynamics of untwisting in three-dimensional configurations, introducing increasingly complex configurations (e.g. Mahlmann et al. Reference Mahlmann, Akgün, Pons, Aloy and Cerdá-Durán2019, Reference Mahlmann, Philippov, Mewes, Ripperda, Most and Sironi2023; Carrasco et al. Reference Carrasco, Viganò, Palenzuela and Pons2019). The literature mostly focuses on dipolar boundary conditions at the stellar surface, and the case of multipolar configurations has had relatively little attention.
The equation describing an aligned non-rotating force-free axisymmetric magnetosphere is the Grad–Shafranov equation, and the so-called pulsar equation (Michel Reference Michel1973) can be seen as its extension to rotating, relativistic, magnetospheres. A few semi-analytical solutions exist for non-rotating twisted magnetospheres (Wolfson Reference Wolfson1995; Thompson et al. Reference Thompson, Lyutikov and Kulkarni2002); see Viganò et al. (Reference Viganò, Pons and Miralles2011) for a review. In particular, self-similar solutions have the property that both poloidal and toroidal components of the magnetic-field decay with the same power law of the radius. These solutions therefore do not asymptotically tend to vacuum. Self-similar solutions allow for multipolar-like fields (Pavan et al. Reference Pavan, Turolla, Zane and Nobili2009), but a single multipole is possible at a given time because the nonlinearity of the equations precludes linear combinations. Much work has recently been done on numerical solutions where the twist is confined, by construction, to a subset of field lines. This presents the advantage of allowing for vacuum fields outside of this region (e.g. Akgün et al. Reference Akgün, Miralles, Pons and Cerdá-Durán2016; Ntotsikas et al. Reference Ntotsikas, Gourgouliatos, Contopoulos and Lander2024, and references therein). Some of this work has been focused on the matching of surface boundary conditions with internal solutions in order to study coupled equilibrium and evolution (Viganò et al. Reference Viganò, Pons and Miralles2011; Fujisawa & Kisaka Reference Fujisawa and Kisaka2014; Glampedakis et al. Reference Glampedakis, Lander and Andersson2014; Pili, Bucciantini & Del Zanna Reference Pili, Bucciantini and Del Zanna2015; Akgün et al. Reference Akgün, Cerdá-Durán, Miralles and Pons2017; Uryu et al. Reference Uryu, Yoshida, Gourgoulhon, Markakis, Fujisawa, Tsokaros, Taniguchi and Zamani2023). In Parfrey, Beloborodov & Hui (Reference Parfrey, Beloborodov and Hui2013), it is shown that beyond a certain critical twist the force-free magnetosphere becomes tearing-mode unstable, forming a current sheet that dissipates the energy through magnetic reconnection. The existence of a critical twist beyond which no stationary solution is possible has been recognised by many authors (e.g. Akgün et al. Reference Akgün, Miralles, Pons and Cerdá-Durán2016, and references therein).
In this work we present a class of semi-analytical solutions to the non-rotating Grad–Shafranov equation akin to a multipolar expansion of the magnetic field. It allows for the choice of an arbitrary set of multipoles and can asymptotically tend to a vacuum dipole. In § 2 we briefly review the mathematical framework, in § 3 we present the class of solutions, in § 4 we present a couple of examples of solutions and we conclude in § 5.
2. Mathematical framework
For a non-rotating star, the force-free condition
$\boldsymbol{j}\times \boldsymbol{B} = 0$
(e.g. Gruzinov Reference Gruzinov2006) leads to
\begin{equation} \boldsymbol{\nabla } \times \boldsymbol{B} = \alpha \boldsymbol{B}, \end{equation}
where
$\alpha$
is a function such that
$\boldsymbol{j} = \alpha \boldsymbol{B}$
is the current density with
$\boldsymbol{B}$
the magnetic field vector.
Taking the divergence of (2.1)
\begin{equation} \boldsymbol{B}\boldsymbol{\cdot }\boldsymbol{\nabla } \alpha = 0, \end{equation}
such that
$\alpha$
is constant along a field line.
In the following we seek cylindrically symmetric solutions, such that
$\partial _\varphi = 0$
where
$\varphi$
is the angle around the symmetry axis. The system is therefore effectively two-dimensional, and in a system of coordinates where one of them labels magnetic-field lines, (2.2) implies that
$\alpha$
depends on that single coordinate. Because of axial symmetry the poloidal part of the magnetic field is itself divergence free. This allows one to express it through a single Euler potential (e.g. Stern Reference Stern1970)
$\mathcal{P}$
\begin{equation} \boldsymbol{B}_{\textrm {p }} = \boldsymbol{\nabla } \mathcal{P} \times \boldsymbol{\nabla } \varphi = \frac{\boldsymbol{\nabla } \mathcal{P} \times \boldsymbol{e}_\varphi }{r\sin \theta } = \frac{r^{-1} \partial _\theta \mathcal{P} \boldsymbol{e}_r - \partial _r \mathcal{P} \boldsymbol{e}_\theta }{ r\sin \theta }, \end{equation}
where
$(r,\theta,\varphi)$
denote spherical coordinates,
$(\boldsymbol{e}_r, \boldsymbol{e}_\theta, \boldsymbol{e}_\varphi)$
is the associated orthonormal base, and
$\mathcal{P}$
labels magnetic-field lines.Footnote
1
It follows that one can write
$\alpha = \alpha (\mathcal{P})$
. The magnetic field reads
$ \boldsymbol{B}= \boldsymbol{B}_{\textrm {p }} + B_{\varphi } \boldsymbol{e}_\varphi$
.
Expending (2.1) in spherical coordinates, we get
\begin{equation} \left (\begin{matrix} \frac{1}{r\sin \theta } \partial _\theta \left (\sin \theta B_{\varphi }\right ) \\[5pt] -\frac{1}{r}\partial _r\left (r B_{\varphi }\right ) \\[5pt] \frac{1}{r} \left (\partial _r \left (r B_{\theta }\right ) - \partial _\theta B_{r}\right ) \end{matrix}\right ) = \alpha \left (\begin{matrix} B_{r} \\[5pt] B_{\theta } \\[5pt] B_{\varphi } \end{matrix}\right ). \end{equation}
From (2.3) we see that the poloidal part of the left-hand side of (2.4) can be written
\begin{eqnarray} \frac{1}{r\sin \theta } \partial _\theta \left (\sin \theta B_{\varphi }\right )\ & = & \ \alpha \frac{\partial _\theta \mathcal{P} }{r^2\sin \theta }, \nonumber \\[5pt] \frac{1}{r}\partial _r\left (r B_{\varphi }\right )\ & = &\ \alpha \frac{\partial _r \mathcal{P}}{ r\sin \theta }. \end{eqnarray}
Since we have seen that
$\alpha$
is a function of
$\mathcal{P}$
, we can introduce the primitive
$A = \int \alpha \mathrm{d} \mathcal{P}$
, and see from (2.5) that
\begin{equation} A = r\sin \theta B_{\varphi }. \end{equation}
Inserting (2.3) into the third component of the force-free equation, (2.4), we obtain the Grad–Shafranov equation
\begin{equation} -\partial _r^2\mathcal{P} - \frac{1-\mu ^2}{r^2} \partial _\mu ^2\mathcal{P} = \alpha (\mathcal{P}) A(\mathcal{P}), \end{equation}
where
$\mu \equiv \cos \theta$
. Boundary conditions must be such that
$\mathcal{P}(\mu =\pm 1) = 0$
, which ensures that the field line going out of either side of the symmetry axis is the same and avoids field-line crossings (Wolfson Reference Wolfson1995).
More generally, the boundary conditions of this equation are set by specifying
$\mathcal{P}$
on the stellar surface, and at infinity. Here and in § 3 we present solutions that are decomposed on a functional basis such that asymptotically
$\mathcal{P} \propto 1/r$
and therefore the external boundary condition is automatically satisfied. The treatment of the inner boundary condition then depends on the class of solutions considered.
It is instructive to look at the example of self-similar solutions to (2.7), see Low & Lou (Reference Low and Lou1990), Lynden-Bell & Boily (Reference Lynden-Bell and Boily1994), Wolfson (Reference Wolfson1995)and Thompson et al. (Reference Thompson, Lyutikov and Kulkarni2002). It consists in using the ansatz
\begin{equation} \mathcal{P} = F(\mu )/r^p \;\; ; \;\; \alpha (\mathcal{P}) = c \mathcal{P}^{1/p}, \end{equation}
where
$c$
is a constant,
$p \geqslant 1$
is a constant index and
$F$
is a function to be determined. An ordinary differential equation for
$F$
is obtained by inserting (2.8) into (2.7). This equation is nonlinear and can be integrated numerically. In order to proceed, one must specify the condition
$F(\mu =1) =\mathcal{P}(\mu = 1) = 0$
and its derivative with respect to
$\mu$
, that is
$F'(\mu =1)$
. The latter corresponds to the radial magnetic field at the pole, as can be seen from (2.3). However, if setting
$(F(1),F'(1))$
determines the initial value problem for numerical integration, there is no guarantee that the boundary value problem
$F(\mu = -1) =\mathcal{P}(\mu = -1) = 0$
will be met as well. In fact, for fixed
$c$
and
$p$
, boundary conditions are verified for a discrete set of
$F'(\mu =1)$
which one must find numerically (for example by a shooting method). Since it is physically more interesting to specify
$F'(\mu =1)$
, one rather varies
$c$
, for which a discrete spectrum of solutions exist. This can be seen as a nonlinear eigenvalue problem of
$c(p)$
. Since (2.7) with the ansatz (2.8) is even,
$\mathcal{P}(-\mu )$
is a solution if
$\mathcal{P}(\mu )$
is, and one can integrate only until the equator at
$\mu = 0$
and replace the boundary condition at
$-1$
by
$\partial _\theta \mathcal{P}(\mu = 0) =0$
(Thompson et al. Reference Thompson, Lyutikov and Kulkarni2002).
3. Solutions
3.1. Ansatz
Multipolar expansions are general solutions of the vacuum electromagnetic field in spherical symmetry (Bonazzola, Mottez & Heyvaerts Reference Bonazzola, Mottez and Heyvaerts2015; Pétri Reference Pétri2015). This motivates us to posit the ansatz
\begin{equation} \mathcal{P} = B_0 R^2 \sum _{i=1} F_i(\mu )\frac{R^i}{r^i}, \end{equation}
where
$R$
is the stellar radius,
$B_0$
gives the scale of the surface magnetic-field intensity and
$F_i$
are functions to be determined.
A heuristic for the induction of the ansatz in (2.8) consists of analysing the dependence on
$r$
on both sides of (2.7). Indeed, once one has posited the ansatz for
$\mathcal{P}$
, it is clear that the left-hand side is
$\propto 1/r^{p+2}$
. From this observation the form proposed for
$\alpha$
in (2.8) appears as the simplest way to balance out the powers of
$r$
. With the current ansatz (3.1), applying this heuristics brings us to postulate the following form for
$\alpha$
and its integral
$A$
:
\begin{eqnarray} \alpha (\mathcal{P}) \ & = &\ R^{-1}\sum _{i=1} (i+1) c_i \left (\frac{\mathcal{P}}{B_0R^2}\right )^i, \end{eqnarray}
\begin{eqnarray} A(\mathcal{P}) \ & = & \ B_0 R \sum _{i=1} c_i \left (\frac{\mathcal{P}}{B_0R^2}\right )^{i+1}, \end{eqnarray}
where
$c_i$
are a priori free coupling constants.
Keeping in mind that
$[\mathcal{P}] = \text{magnetic strength} \times \text{length}^2$
, it follows that the
$F_i$
are dimensionless. Similarly, from (2.1) one sees that
$[\alpha ] = \text{length}^{-1}$
, such that the coupling constants
$c_i$
are also dimensionless. Thereafter we use units such that
$B_0=R=1$
.
3.2. Source term
$\alpha A$
The source term of (2.7) reads
\begin{equation} \alpha A = \sum _{i,j=1} (i+1) c_i c_j \mathcal{P}^{i+j+1} = \sum _{k = 3} \sum _{i=1}^{k-2} c_i c_{\underline {j}} (i+1) \mathcal{P}^{k}, \end{equation}
where
$k = i+j+1$
and
$\underline {j} = k -i-1$
.
We formulate the powers of
$\mathcal{P}$
as a Laurent series in
$r$
\begin{eqnarray} \mathcal{P}^k & = & \sum _{i_1, \ldots ,i_k}\frac{F_{i_1}\ldots F_{i_k}}{r^l} \text{ with } l = \sum _m^k i_m \end{eqnarray}
\begin{eqnarray} & = & \sum _{l=k}\frac{G_{l}^{(k)}}{r^l}, \end{eqnarray}
with
\begin{eqnarray} G_{l}^{(k)} & = & \sum _{i_1+\ldots +i_k = l; i_j \geqslant 1} F_{i_1}\ldots F_{i_k} \end{eqnarray}
\begin{eqnarray} & = & \sum _{i_1 = 1}^{\hat {i}_1}\ldots \sum _{i_{k-1} = 1}^{\hat {i}_{k-1}} F_{i_1}\ldots F_{\underline {i}_k}, \end{eqnarray}
where we have
$\underline {i}_k = l - \sum _{j = 1}^{ k-1} i_j$
and
${\hat {i}}_j = l - \sum _{m=1}^{j-1} i_m - (k-j)$
. In particular, it results that
$G_{k}^{(k)}$
=
$ F_1^k$
and, more generally, one can show the following functional dependence:
\begin{eqnarray} G_{k+m}^{(k)} & = & f\left (F_1, \ldots , F_{m+1}\right ). \end{eqnarray}
This results from
$\forall k, \max (\underline {i}_k) = m+1$
.
We can now decompose the source term into a Laurent series with explicit coefficients. First, inserting (3.5) into (3.4), we obtain
\begin{equation} \alpha A = \sum _{k = 3} \sum _{i=1}^{k-2} c_i c_{\underline {j}} (i+1) \sum _{l=k} \frac{G_{l}^{(k)}}{r^l}. \end{equation}
Second, we swap sums such that
\begin{equation} \alpha A = \sum _{l = 3} \frac{1}{r^l}\sum _{k=3}^{l}\sum _{i=1}^{k-2} c_i c_{\underline {j}} (i+1) G_{l}^{(k)} = \sum _{l = 3} \frac{1}{r^l} [\alpha A]_{l}, \end{equation}
where
$ [\alpha A]_{l} = \sum _{k=3}^{l}\sum _{i=1}^{k-2} c_i c_{\underline {j}} (i+1) G_{l}^{(k)}$
denotes the angular part of the
$\bigcirc (r^{-l})$
term of the function
$\alpha A$
. In order to generate the subsequent orders of the source we see from (3.11) that
\begin{equation} \left [\alpha A \right ]_{l+1} = \left [\alpha A \right ]_{l} (l \rightarrow l+1) + G_{l+1}^{(l+1)} \sum _{i = 1}^{l-1} c_i c_{\underline {j}} (i+1), \end{equation}
where the arrow indicates replacement of the lower indices of the
$G$
functions, and
$\underline {j} = l-i$
.
Here, we explicit the first 5 orders of the source term
\begin{align}\left [\alpha A \right ]_3& = 2c_1^2 G_{3}^{(3)},\\[-8pt]\nonumber\end{align}
\begin{align}\left [\alpha A \right ]_4 &= 2c_1^2 G_{4}^{(3)} + 5c_1 c_2 G_{4}^{(4)},\\[-8pt]\nonumber \end{align}
\begin{align}\left [\alpha A \right ]_5 & = 2c_1^2 G_{5}^{(3)} + 5c_1 c_2 G_{5}^{(4)} + (6c_1 c_3 + 3c_2^2) G_{5}^{(5)},\\[-8pt]\nonumber \end{align}
\begin{align}\left [\alpha A \right ]_6 & = 2c_1^2 G_{6}^{(3)} + 5c_1 c_2 G_{6}^{(4)} + (6c_1 c_3 + 3c_2^2) G_{6}^{(5)} + 7(c_1c_4 + c_2c_3) G_6^{(6)},\\[-8pt]\nonumber \end{align}
\begin{align}\left [\alpha A \right ]_7 & = 2c_1^2 G_{7}^{(3)} + 5c_1 c_2 G_{7}^{(4)} + (6c_1 c_3 + 3c_2^2) G_{7}^{(5)} + 7(c_1c_4 + c_2c_3) G_7^{(6)}\nonumber\\[5pt]&\quad + \left [8(c_1 c_5 + c_2 c_4) + 4c_3^2\right ] G_{7}^{(7)}.\end{align}
3.3. Hierarchy
Inserting (3.11) into the Grad–Shafranov equation (2.7), and solving for each coefficient of the Laurent series in
$r$
, it transforms into a set of ordinary differential equations
\begin{equation} \forall i \geqslant 1, -i(i+1) F_i - (1-\mu ^2) F_i'' = [\alpha A]_{l}\left (\{F_j\}_{j\leq i}\right ), \end{equation}
where
$l=i+2$
and, by virtue of (3.9), the right-hand side only depends on
$F_{j \leq i}$
. As a consequence of the hierarchy of this set of equations, it can be solved iteratively up to a certain truncation order.
The boundary conditions
$\mathcal{P}(\mu =\pm 1) = 0$
translate into
$\forall i\geqslant 1, F_i(\pm 1) = 0$
. Indeed, we need
$\mathcal{P}(\mu =\pm 1) = 0$
for all
$r$
along the symmetry axis which, given the definition (3.1), implies that all
$F_i$
must individually cancel at the boundaries. In practice, we will build solutions by specifying
$\{F_i'(\mu =1)\}_{i\geqslant 1}$
(see below).
When all coupling constants are null,
$c_i=0$
, then each element of (3.18) gives the corresponding vacuum multipole (see Appendix A):
$F_1$
is a dipole,
$F_2$
a quadrupole,
$F_3$
an octupole, and so on.
If
$c_1 = 0$
then all equations are linear (albeit with a potentially cumbersome source term). If the dipole order is present, that is
$F_1 \neq 0$
, then it is a vacuum dipole as its source term is empty. Provided that some other coupling constant is non-zero, then this dipole sources an infinity of higher-order equations. For example, if
$c_2$
is active, then the source term of the equation for
$F_3$
is
$[\alpha A]_{5} = 3c_2^2 G_{5}^{(5)} = 3c_2^2 F_1^5$
, according to (3.15). Since all higher-order terms decay with radius faster that the
$F_1$
term, the solution is asymptotically a vacuum dipole. More generally, for
$c_1=0$
solutions asymptotically tend to the first non-zero vacuum multipole at infinity. Since equations become linear, they can in principle be solved analytically using power series. However, we here mostly use numerical integration as it appears simpler for our purposes.
On the other hand equations are highly nonlinear when
$c_1 \neq 0$
. For example, the source term of the leading-order equation for
$F_1$
is
$2c_1^2 F_1^3$
according to (3.13). As a result, these solutions do not asymptotically connect to vacuum.
3.4. Boundary conditions, current sheet and regular or anti-twist
3.4.1. Boundary conditions
We need to satisfy the boundary conditions given by
$\forall i\geqslant 1, F_i(\pm 1) = 0$
(see § 3.3). In principle, one possibility would be to integrate from the north pole at
$\mu =1$
with
$\{F_i(\mu =1)=0\}$
until the south pole at
$\mu = -1$
and adjust the
$\{F_i'(\mu =1)\}$
, or alternatively the constants
$\{c_i\}$
, such that the second boundary condition
$\{F_i(\mu =-1)=0\}$
be met. This, in fact, would be the generalisation of the method used to obtain self-similar solutions, § 2. However, we have empirically found that this approach leads to the divergence of the
$\{F_i\}$
series even in the simplest cases (e.g. § 4.1).
One can see that for any solution
$\{F_i(\mu )\}_{i\geqslant 1}$
of (3.18) on the interval
$[-1,0[$
,
$\{F_i(-\mu )\}_{i\geqslant 1}$
is a solution on the interval
$]0,1]$
. This even symmetry with respect to the equatorial plane allows us to build a solution that fulfils boundary conditions at the poles and that is defined everywhere except in that plane. To this end, one needs to integrate (3.18) with boundary conditions
$\{F_i(\mu =1)=0, F_i'(\mu =1)\}_{i\geqslant 1}$
in, say, the northern hemisphere (
$\mu \gt 0$
), and then mirror it onto the other hemisphere. Here,
$\{F_i'(\mu =1)\}$
are free parameters. This procedure comes at the price of an equatorial current sheet, as explained below, but makes the satisfaction of the boundary conditions
$\forall i\geqslant 1, F_i(\pm 1) = 0$
straightforward (see § 3.3).
The case where
$c_1$
is the only non-zero constant is particular. Indeed, if all
$F_i'(\mu =1) =0$
except
$F_1'$
, that is the dipolar order, then the hierarchy stops at
$F_1$
(i.e.
$F_{i\gt 1}= 0$
) and we are in the nonlinear case already studied in Low & Lou (Reference Low and Lou1990). In this case, there is a discrete spectrum of values of
$c_1$
(or
$F_1'(\mu =1)$
) for which the solution satisfies the boundary conditions and that is continuous in the equatorial plane. It is possible to add higher orders to such solution, although at the price of a discontinuity at these orders. We expand on this in § 4.2.
3.4.2. Regularly and anti-twisted solutions
Two such solutions are in fact possible: for a given set of coupling constants
$\{c_i\}_{i \geqslant 1}$
in the northern hemisphere one can choose
$\{\pm c_i\}_{i \geqslant 1}$
in the south. Indeed, the right-hand side of (3.18) depends only on products
$c_i c_j$
such that a change of sign of all coupling constants has no effect. However, the sign of
$B_\varphi$
is then reversed across the equator through (2.6). We shall refer to solutions with reversed
$B_{\varphi }$
as anti-twisted solutions and others as regularly twisted.
3.4.3. Solutions with a current sheet
The solution is constructed as follows:
\begin{equation} \mathcal{P} = \left \{\begin{array}{cc} \mathcal{P}_{\textrm {n}}(\mu ) \text{ for } \mu \gt 0 ,\\ \mathcal{P}_{\textrm {s}} = \mathcal{P}_{\textrm {n}}(-\mu ) \text{ for } \mu \lt 0 ,\end{array}\right . \end{equation}
where
$\mathcal{P}_{\textrm {n}}$
is the potential of the northern hemisphere obtained by integration of (3.18) from
$\mu = 1$
to 0, while the southern hemisphere
$\mathcal{P}_{\textrm {s}}$
is obtained by symmetry. As a result of this even symmetry,
$\mathcal{P}(0^{+}) = \mathcal{P}(0^{-})$
and therefore the magnetic component normal to the equator,
$B_\theta \propto \partial _{r} \mathcal{P}$
, is continuous across it (see (2.3)). On the other hand, the radial component
$B_{ r} \propto \partial _{\theta } \mathcal{P}$
, tangent to the equatorial plane, is generally discontinuous since
$\partial _{\theta } \mathcal{P}(0^{+}) = - \partial _{\theta } \mathcal{P}(0^{-})$
. The exception to this rule is when
$B_{ r} = 0$
, such as in the case of a vacuum dipole for example. As a result, the even symmetry produces an azimuthal current sheet with surface current density
\begin{equation} \sigma _{\varphi } = -2B_r(\mu =0^{+}). \end{equation}
For anti-twisted solutions, everything is identical except for the sign of the toroidal component
$B_\varphi$
which reverses across the equator. This is supported by a radial component of the current sheet with density given by
\begin{equation} \sigma _r = \left \{\begin{array}{cc} -2B_r(\mu =0^{+}) \text{ (anti-twisted)},\\ 0 \text{ (regularly twisted)}. \end{array} \right . \end{equation}
3.4.4. Practical consequences
For odd-order
$F_i$
functions, vacuum solutions are unaffected since
$B_r = 0$
in the equatorial plane (vacuum dipole, octupole, etc.). This is because these functions are odd (see Appendix D1). For example, if
$c_1=0$
,
$F_1$
is always a vacuum dipole which does not generate a current sheet. However, particular solutions of higher orders do not share the property
$B_r = 0$
which results in a toroidal current sheet
$\sigma _\varphi$
as seen above. For example, a vacuum dipole sources higher orders with equatorial discontinuities supported by the current sheet. In the case with coupling constants
$c_{i\neq 2}=0$
, the dipole sources octupolar and higher orders,
$F_{i\geqslant 3}$
. This case is detailed in Appendices. As a result, the current sheet is localised in the sense that it is born by components
$1/r^{i\geqslant 5}$
.
Even-order vacuum solutions
$F_i$
are odd functions. As a consequence, even symmetry makes them become what we may call split multipoles. For example, a vacuum quadrupole, which is possible provided that
$c_1=0$
, becomes a split quadrupole. As before, this is supported by the toroidal component of the current sheet. Discontinuities in the particular solutions behave the same.
In order to keep the natural symmetry of the quadrupole and simultaneously satisfy boundary conditions, one would need to impose odd parity to all orders. However, the quadrupole sources higher-order equations in the hierarchy. Particular solutions to these equations, contrary to vacuum ones, will not in general satisfy
$B_\theta (\mu =0)=0$
. As a result, discontinuities arise in the magnetic component normal to the equatorial plane, which cannot be supported by a current sheet and is unphysical.
3.5. Convergence
A necessary condition for convergence of the series defining the potential
$\mathcal{P}$
, (3.1), is that the source term
$\alpha A$
defined by (3.11) itself converges. Indeed one expects that
$F_i/r^l \sim [\alpha A]_{l}/r^l$
with
$l=i+2$
from (3.18). We conjecture that this is a sufficient condition in general as well, without being able to demonstrate it at this point. Nonetheless, all the cases studied in § 4 are numerically shown to be converging. The special case for
$c_{i\neq 2} = 0$
is discussed below and in Appendices.
Since with
$c_{i\neq 2} = 0$
all equations are linear beyond first order, solutions are the sum of a particular and a vacuum solution. Vacuum solutions are proportional to the initial condition
$F_i'(\mu = 1)$
(see Appendices). Therefore, a necessary condition for convergence is that
$F_i'(\mu = 1) \rightarrow _{i\rightarrow \infty } 0$
. For the rest, the source term tending to zero implies the particular solutions also tending to zero.
Concerning the
$F_i$
functions, we have not had any convergence or singularity issue in all the cases that have been studied numerically. We have analytically studied the particular case of the function
$F_3$
with
$c_{i\neq 2}=0$
, § 5, shows that indeed the solution is analytical.
4. Applications
In this section, we first study the simplest case that produces (i) asymptotically a vacuum dipole and (ii) the leading-order toroidal component. Such configurations are obtained for
$c_2 \neq 0$
and
$c_{i\neq 2} = 0$
. In all the cases presented, we present plots corresponding to the anti-twisted solutions unless otherwise mentioned. This choice is only made to limit the number of plots to a minimum. Indeed, the regularly twisted configurations can be straightforwardly deduced from the anti-twisted ones by symmetrising with respect to the equator, and setting the radial current-sheet component to zero, that is
$\sigma _r =0$
.
We also discuss the case
$c_{i\neq 1} = 0$
, where the dipole order is not in vacuum, but rather interacts nonlinearly. This case was already studied in previous literature but without the possibility of adding multipoles. Here we exhibit such an example.
For practical purposes it is convenient to carry out the integration of (3.18) as a function of
$\bar \mu = 1 -\mu$
. In the following we denote by a dot the derivative with respect to
$\bar \mu$
, and unless otherwise stated the boundary conditions are denoted
$\dot F_i = -F_i'(\mu = 1)$
.
In all cases we give the maximum twist, that is the largest azimuthal shift between the two footpoints of a field line,
$\Delta \varphi = \oint _{\text{field line}} \mathrm{d}\varphi$
, and the magnetic helicity
$H = \int \mathcal{A} \cdot \boldsymbol{B} \;\mathrm{d} V$
. We also numerically compute the magnetic energy as
$E = \int B^2 \mathrm{d}V$
, which we compare with the energy of the vacuum multipoles alone (i.e.
$\forall i, c_i = 0$
). Values for energy are given in units of
$2\pi B_0^2R^3$
and values of helicity are in units of
$2\pi B_0^2 R^4$
. The
$2\pi$
factor accounts for the azimuthal integration. For example, a vacuum dipole has exactly
$E_{\mathrm{dip}}=1/3$
for
$B_1 = 1$
.
The python script used to compute the solutions and produce the figures presented in this paper is made available on the dedicated Zenodo repository.Footnote
2
Unless otherwise stated, solutions to (3.18) for the
$F_i$
functions have been obtained numerically using the Runge–Kutta algorithm of order 4 with order 5 for error estimation implemented in the SciPy libraryFootnote
3
4.1. General case with
$c_{i\neq 2} = 0$
The most general solution can be expressed up to
$\bigcirc (1/r^5)$
as
\begin{equation} \boldsymbol{B} = \boldsymbol{B}_{1}^{\textrm {v}} + \boldsymbol{B}_{2}^{\textrm {v, split}} \pm c_2 \frac{B_1^3}{8} \frac{(1-\mu ^2)^{5/2}}{r^4} \boldsymbol{e}_\varphi +\bigcirc \left (\frac{1}{r^5}\right ), \end{equation}
where
$\boldsymbol{B}_{1}^{\textrm {v}}, \boldsymbol{B}_{2}^{\textrm {v, split}}$
are the vacuum dipole and split quadrupole, respectively, which we make explicit in Appendices, (A6) and (A9). Here,
$B_1$
is the dipole strength at the pole. The dipole sources a toroidal field at order
$1/r^4$
through (2.6), while both vacuum fields are purely poloidal. The
$\pm$
sign in front of the toroidal term signifies the potential parity operation between hemispheres: the sign can be reversed across the equator for anti-twisted solutions. If so the solution generates a current sheet in the radial direction with surface density
$\sigma _r = -2 c_2 B_1^3 (1-\mu ^2)^{5/2}/r^4$
, as obtained in § 3.4.3.
The twist of a field line is the difference of azimuth
$\Delta \varphi = \int \mathrm{d}\varphi$
between its two footpoints integrated along the line. In absence of split quadrupole, the twist of a magnetic-field line emerging at colatitude
$\theta$
is straightforward to calculate in the approximation of (4.1)
\begin{eqnarray} \Delta \varphi (\theta ) & = & 2\oint _{\theta }^{\pi /2} \frac{B_{\varphi }}{B_{\theta }} \frac{\mathrm{d} \theta '}{\sin \theta '} \end{eqnarray}
\begin{eqnarray} &&\qquad\qquad = \frac{1}{2} \frac{c_2 B_1^2}{R_*} \cos \theta \sin ^2\theta , \end{eqnarray}
where in those units
$R_*=1$
. The factor 2 in (4.2) implies that this is the full twist from one hemisphere to the other. In the case of the anti-twisted configuration the hemisphere-to-hemisphere twist is by definition 0, but the twist between a pole and the current sheet is half the above value.
Helicity is defined as
$H = \int \mathcal{A} \boldsymbol{\cdot }\boldsymbol{B} \;\mathrm{d} V$
where
$\mathcal{A}$
is the magnetic potential vector. In the present case, one can show that (Appendix D1)
\begin{equation} H = 2\int \frac{\mathcal{P} B_\varphi }{r \sin \theta } \;\mathrm{d} V = \frac{4}{105} c_2 B_1^4, \end{equation}
where the second equality is valid for the approximation of (4.1). As for the twist, this expression is valid only in the regularly twisted configuration, and
$H=0$
in the anti-twisted case.
At order
$1/r^5$
, the equation for
$F_3$
is sourced by
$3c_2 G_5^{(5)} = F_1^5$
, (3.15). As shown in Appendices, this results in a solution with discontinuous
$F_3'$
at the equator due to imposed even parity. This discontinuity is sustained by the toroidal component of the current sheet, § 3.4.3, which therefore also decreases as
$1/r^5$
. Appendices also give an analytical solution of the particular solution for
$F_3$
. Beyond this order, we consider in this paper that it is more economical to use numerical integration. Provided a quadrupole is also present, that is
$F_2\neq 0$
, a toroidal field will also be present at order
$r^{-5}$
expressed by
$\pm c_2F_1^2F_2/r^5\sin \theta$
. More generally, none of the functions
$F_{i\geqslant 3}$
are vacuum solutions, but instead are solutions of (3.18) sourced by combinations of functions of lower order. These higher-order functions,
$F_{i\geqslant 3}$
, are the ones that generate the equatorial discontinuity supported by the toroidal component of the current sheet. In the particular case where
$F_2 = 0$
and
$\forall i \gt 1, \dot F_{2i} = 0$
one can show that all even-order functions are null, that is
$\forall i \gt 1, F_{2i} = 0$
.
4.1.1. Dipole
In figure 1 we show the anti-twisted solution sourced by a vacuum dipole. All multipoles have been computed up to order 30 and are shown in figure 2. Since no quadrupole is present, only odd orders
$F_1, F_3, F_5\ldots$
contribute. We can see in figure 1 that at the surface the dipole component dominates up to a colatitude of
${\sim}60^{\circ }$
and is dominated by higher-order multipoles around the equator. Indeed, one can see in the top-right panel of figure 1 that the radial magnetic component is far from following a dipolar behaviour as it flips sign within the northern hemisphere. This is primarily due to the influence of the octupolar component,
$F_3$
, at the surface, figure 2. However, since the latter component decays as
$1/r^5$
it gives way to the dipolar component generated by
$F_1$
after a few stellar radii. Thus, the solution is asymptotically dipolar.
The largest twist angle is found to be
$\Delta \varphi \simeq 1.20$
rad, close to the approximation of (4.3),
$\Delta \varphi \simeq 1.15$
rad obtained at
$\theta _{\max }\simeq 55$
deg. Similarly, magnetic helicity has
$H \simeq = 0.16$
, compared with the analytical value of 0.23. The magnetic energy is remarkably close to the vacuum-dipole energy, indeed,
$E/E_{\textrm {dip}} \simeq 0.9995$
(numerically significant). Approximately 12 % of the magnetic energy is stored in the toroidal component, mostly from its leading-order contribution (4.1), which implies that the sourced higher multipoles partly cancel the dipolar component such that the total energy actually decreases slightly.
The value of
$c_2 = 6$
is about the largest value that we found to converge, everything else being equal. The series of
$F_i$
displays geometric convergence as can be seen in figure 2, as expected from the arguments given in Appendices. We also show in Appendices that in this case the problem only depends on the parameter
$\rho = 3(\dot F_1/2)^{4}c_2^2$
up to a scaling factor
$\dot F/2$
. This means that for a constant value of
$\rho$
,
$F_{i}/(\dot F_1/2)$
are identical.
Solution for
$\dot F_1= B_1 =1, c_2=6$
between 1 and 5 stellar radii, where
$B_1$
is the dipole field strength at the pole. All quantities are multiplied by
$r^3$
for better visualisation. Left, left hemisphere: poloidal cross-section of the current amplitude
$\alpha |\boldsymbol{B}| r^3$
, normalised by its maximum. Left, right hemisphere: poloidal cross-section of the toroidal field
$B_\varphi r^3$
, normalised by its maximum. Poloidal magnetic-field lines are shown as solid black lines, and vacuum-dipole field lines are shown for comparison as grey dashed lines. Top right: magnetic-field components at the stellar surface,
$B_r$
dashed line,
$B_\theta$
dot-dashed line and
$B_\varphi$
dotted line in units of
$B_1$
. For comparison, vacuum-dipole components are plotted as grey lines with corresponding styles. Bottom right: current-sheet current components
$\sigma _r$
, dashed line, and
$\sigma _\varphi$
, dotted line, rescaled by a factor
$r^3$
.
Components (left) and convergence (right) of the solution to (3.18) for the conditions of figure 1. Left: the components
$F_i$
up to the truncation at order 30. For clarity only the lowest 9 orders are shown in colour and labelled, and the remaining ones are shown in grey. The thick black line represents their sum up to order 30 such that it represents the potential
$\mathcal{P}$
, (3.1), on the stellar surface at
$r=1$
. Right: convergence of the series defining
$\mathcal{P}$
, (3.1), for the conditions of figure 1. The maximum of each
$|F_i(\mu )|$
for
$\mu$
between 0 and 1 is plotted against its order
$i$
, where
$F_i$
is solution of (3.18). Only odd orders are shown, as even orders are all equal to 0.
Same as figure 1 for
$ \dot F_1 =1, \dot F_3 =4, c_2=3.2$
in the anti-twisted configuration.
The components
$F_{i \gt 1}$
in figure 2 are all particular solutions, as a result of the choice of boundary conditions given by
$\dot F_{i \gt 1}=0$
. Indeed, a given solution is a linear combination of a vacuum and a particular solution,
$F_i = F_i^{\textrm {v}} + F_i^{\textrm {p}}$
. The power-series expansion of particular solutions are proportional to
$\bar \mu ^{l} = (1-\mu )^{l}$
where
$l=i+2$
(see Appendices). It follows that particular solutions have derivatives equal to zero at
$\mu = 1$
, that is
$\dot F_i^{\textrm {p}}{} = 0$
. On the other hand, vacuum solutions all have a term
$\bigcirc (\bar \mu )$
which means that the boundary condition
$\dot F_i$
is entirely determined by the vacuum term.
The higher the order of
$F_i$
the more their effect is localised around the equator. This is a direct consequence of the fact that
$F_i \propto \bar \mu ^{i+2}$
, which is visible in figure 2.
4.1.2. Dipole + octupole
In figure 3 we show the anti-twisted configuration of a vacuum dipole with the contribution of a vacuum octupole. We used the values
$F_1' =-1, F_3' =-4 c_2=3.2$
, where
$c_2$
has about the largest value we found to allow for convergence in this case. Figure 4 shows the components separately, and one can clearly distinguish the octupole contribution,
$F_3$
, which slope at the origin is not null, meaning a vacuum contribution is present on top of the particular solution sourced by the dipole. The convergence plot in figure 4 shows a monotonic convergence of
$F_{i\gt 5}$
while the lower orders are dominated by the effect of the values imposed on
$F_1'$
and
$F_3'$
.
The result is a peak of the toroidal field at the surface at higher latitudes, around
$\pm 45^{\circ }$
. These contributions are, however, much more localised than in the dipole case above, as they result mainly from terms of the form
$\propto c_2^2 F_3 F_1^2/r^6$
. More generally, the slow convergence means that many multipoles contribute, especially near the equator. This is particularly visible in the very steep decay of the toroidal component of the current sheet. Beyond a few stellar radii the structure of the field tends to that of figure 1. We note that, in this case as in the others, it is possible to produce a vanishing current sheet on the surface by fitting
$\{\dot F_i\}$
so as to produce a continuous boundary condition, while the current sheet appears at higher altitude. This is particularly what is being done in § 4.3.
The largest twist is here
$\Delta \varphi \simeq 0.78$
rad, and helicity
$H \simeq 0.18$
. Similarly to the dipole case in § 4.1 the magnetic energy is very close to the vacuum energy (dipole + octupole) with
$E/E_{\textrm {dip + oct}} \simeq 0.991$
(numerically significant). Approximately 4 % of the magnetic energy is stored in the toroidal component. However, the absolute amount of energy here represents
$1.28 E_{\textrm {dip}}$
due to the octupolar component. As a result, the absolute toroidal energy is larger than in § 4.1.1 at approximately
$0.14 E_{\textrm {dip}}$
. Although less twisted, this multipolar configuration represents a somewhat larger energy reservoir.
4.2. Non-vacuum dipole + octupole
The case
$c_{i\neq 1} = 0$
has already been partially studied in Low & Lou (Reference Low and Lou1990). Indeed, for
$\dot F_{i\neq 1} = 0$
, it is equivalent to the self-similar model discussed in § 2 for index
$p = 1$
, and
$F_1$
is the only function that is not null. There is a discrete spectrum of values of
$c_1$
for which the function is fully continuous at the equator and does not require a current sheet. For
$c_1 =0$
one has a vacuum dipole, and for higher values the angular dependency shows an increasing number of poles, but with a radial dependency that remains
$r^{-3}$
for all components, including the toroidal one. The difference in this work, is that it is possible to add multipoles, for example by setting
$\dot F_3 \neq 0$
. This then triggers a cascade of higher-order multipoles, similarly to the above cases.
This behaviour is illustrated in figures 5 and 6 where we have used
$\dot F_1=1 ; \dot F_3=2$
at
$\mu = 1$
and
$c_1 = 15.9$
, with all other initial conditions and constants equal to zero. The value of
$c_1$
is determined by numerically solving
$\dot {F}_1(\mu =0, c_1) =0$
, which is the sufficient condition of
$F_1$
to be even with respect to the equator. With this particular value of the constant
$F_1$
is continuous at the equator and does not generate a current sheet. This is why here we show the twisted case instead of the anti-twisted one. A current sheet with only a toroidal component is, however, still supported by the higher orders, mainly
$F_3$
.
The maximum twist in this configuration reaches
$\Delta \varphi \simeq 5.8$
rad. This extreme value is related to the homogeneous radial dependency of all components much more than to the octupolar component. Indeed, it is very close to the value obtained without the additional octupole. Helicity is
$H \simeq 0.088$
.
The magnetic energy compared with the vacuum case (dipole + octupole) is
$E/E_{\textrm {dip+oct}} \simeq 1.05$
. Consistently with the large twist, the toroidal energy is proportionally larger and represents approximately 27 % of the total magnetic energy. Contrary to the two other cases, the energy is larger than in the vacuum case, and the difference between the two is also more significant. We conjecture that this is due to the nonlinear nature of this configuration where, contrary to the two previous cases, the top-level dipole is itself modified compared with vacuum.
Same as figure 1 for
$\dot F_1 =1, \dot F_3 =2, c_1=15.9$
and regular twist.
4.3. Surface boundary conditions, energy and twist
We have so far considered solutions where the value of
$\mathcal{P}$
on the surface of the star was a by-product of the solution rather than a boundary condition. This is because the boundary conditions of the problem are conveyed by the constants
$\{\dot F_i\}$
which are not straightforwardly connected to the surface potential. On the other hand, given that the inner boundary condition is determined by an infinite number of degrees of freedom,
$\{\dot F_i\}$
, we can conjecture that this class of solutions allows for arbitrary conditions at the stellar surface. This is not, for example, the case of self-similar solution in § 2.
Here, we take the opposite route: the surface potential is fixed to the vacuum-dipole value, i.e.
$\mathcal{P} (r=1, \mu ) = \mathcal{D}(\mu )$
where
$\mathcal{D}=F_1^{\textrm {v}}$
is analytically defined by (A3). As before we set
$B_1=\dot F_1^{\textrm {v}}=1$
such that
$\mathcal{D} \sim 1$
. To do so we use a Nelder–Mead minimiser (see footnote 4) in order to fit the
$\{\dot F_i\}$
such that the boundary condition is fulfilled. As a metric, we use
$\chi ^2 = \sum _j (\mathcal{P}(r=1, \mu _j) - \mathcal{D}(\mu _j))^2$
where we take 100 evenly spaced
$\mu _j$
to sample the function. The task is potentially difficult due to the large number of parameters. Fortunately, we found that a limited number of
$\dot F_i$
is sufficient for a reasonable accuracy. In practice, for
$c_2 \lt 4$
we limited ourselves to twenty multipolar orders, but fitted only
$\dot F_{i\leq 10}$
while
$\dot F_{i\gt 10}=0$
. For
$c_2 \geqslant 4$
, we used 28 orders and fitted the first 14. This is because as the coupling constant increases the source terms remain significant at higher orders. The absolute difference with the target boundary condition is
$\max (\|\mathcal{P} - \mathcal{D}\|) \leq 0.0001$
.
Left axis: relative total magnetic energy
$(E - E_{\textrm {dip}})/E_{\textrm {dip}}$
(orange ‘x’ markers), relative toroidal magnetic energy
$E_{\textrm {toro}}/E_{\textrm {dip}}$
(green stars), dipolar component strength with respect to a vacuum dipole
$\dot F_1/\dot F_1^{\textrm {v}}$
(red circles) and helicity
$H$
in units of
$B_0^2R^4$
(purple diamonds). Right axis: maximum twist angle (blue ‘+’ markers). All quantities are computed for a sequence of solutions characterised by the value of
$c_2$
and with
$c_{i\neq 2}=0$
. These solutions share the same boundary condition
$\mathcal{P}(r=1)=\mathcal{D}$
at the stellar surface. Dashed lines represent the twist (blue), helicity (purple) and relative toroidal energy (green) computed using the analytical approximation in (4.1) with the correct boundary conditions
$\dot F_1 = \dot F_1^{\textrm {v}}$
. Dotted lines correspond to a slightly stronger dipole
$\dot F_1 = 1.048\dot F_1^{\textrm { v}}$
.
The results are shown in figure 7, and compared with the analytic approximation of (4.1). The coupling constant
$c_2$
ranges in
$[0, 4.5]$
. For larger values we found that the solutions do not converge with this boundary condition. We found that the maximum twist angle and helicity are nearly linear functions of
$c_2$
and the toroidal energy a quadratic function. This is in agreement with the analytic approximation, which is particularly good for smaller values of
$c_2$
. As seen in other solutions (Akgün et al. Reference Akgün, Miralles, Pons and Cerdá-Durán2016), the dipolar component increases with increasing currents and as such is largely responsible for the increase in total magnetic energy. Indeed, we can see that the total energy increases faster than the toroidal component, while the increase is exactly the same in the analytical model, since the dipole remains constant and there are no additional multipole. In the exact solution, multipolar orders also increasingly contribute to the energy budget. Finally, we stress that, although the dipolar component increases by
${\sim}20\,\%$
compared with vacuum, the surface field is nonetheless that prescribed by the boundary condition, where other components effectively screen the enhanced dipole.
We performed a fit of the toroidal energy with the analytical model in (4.1) in order to find the effective dipolar amplitude
$B_1$
that allows us to best reproduce the data. We found
$B_1 \simeq 1.048$
to work best, and it allows us to reproduce quite accurately not only the toroidal energy but also helicity as well as the maximum twist. Thus, we may say that as far as these quantities are concerned, the exact solution effectively behaves as the analytical approximation for a surface boundary condition given by a vacuum dipole approximately 5 % stronger.
5. Discussion and conclusion
We present in this work a class of solutions to the Grad–Shafranov equation based on a multipolar expansion the angular part of which is determined by a hierarchy of ordinary differential equations as exhibited in (3.18). Solutions are characterised by a set of arbitrarily chosen coupling constants
$\{c_i\}_{i\geqslant 1}$
and initial conditions for the derivatives of the angular functions
$\{\dot F_i\}_{i\geqslant 1}$
. The main limitation in the choice of couplings is the necessary convergence of the series. The multipolar vacuum expansion is recovered when all coupling constants are zero.
In § 4 we present a subset of solutions symmetric about the equatorial plane and supported by an equatorial current sheet (in particular figures 1, 3 and 5). The twist, or equivalently the toroidal component of the magnetic field, can have opposite signs in opposite hemispheres. However, the physical interpretation of such anti-twisted configurations is unclear.
Current sheets have been previously obtained (see in particular the numerical work of Parfrey et al. (Reference Parfrey, Beloborodov and Hui2013) and references therein), and are expected to form past a certain level of twisting. Here, however, the current sheet is part of the stationary solution, albeit presumably dynamically unstable to reconnection. The case of reversed toroidal fields across the equator is reminiscent of the formation of similar internal magnetar configurations due to Hall drift (Viganò et al. Reference Viganò, Pons and Miralles2012).
A case of interest is when the coupling constants are chosen such that the magnetosphere asymptotically tends to vacuum. In such a case, the hierarchy of equations is linear, and can possibly be solved through power-series expansions. However, the complexity of the source terms makes this approach rather unwieldy and one may prefer a numerical solution. To this effect, we share the python script used for this work (see footnote 2).
We show in (4.1) the simplest solution truncated at order
$r^{-5}$
that asymptotically tends to a vacuum dipole. This results in a simple analytical expression where the leading-order toroidal components intervenes at order
$r^{-4}$
and is stronger near the equator. This approximation is valid at distance sufficiently large to neglect multipolar-order contributions, but is also a good approximation for low to moderate values of the coupling constants of solutions with a dipolar inner boundary condition, reproducing well the twist and toroidal magnetic energy, in particular (§ 4.3).
The class of solutions described in this paper has infinity degrees of freedom, both in terms of coupling constants
$c_i$
and initial condition
$\dot F_i$
for the angular functions. Both can in principle be adjusted in order to produce a particular surface field, for example to match an internal solution. The question of whether this freedom allows one to reproduce any desired surface condition, that is whether the multipolar expansion forms a complete set of functions, is open. In addition, this method can become impractical if too many components are needed. So far, our attempts at fitting dipolar boundary conditions in § 4.3 was quite straightforward, involving only a limited number of parameters.
Thus, we could study the evolution of the magnetosphere as a function of current coupling for fixed boundary conditions. As in other works (e.g. Akgün et al. Reference Akgün, Miralles, Pons and Cerdá-Durán2016) we find a
${\sim} 25\,\%$
increase of the dipolar component with increasing currents, as well as enhanced multipolar orders, leading to a similar increase in total energy. Toroidal energy increases less,
${\sim} 15\,\%{-}20\,\%$
reflecting the amplification of poloidal fields. The maximum twist increases quasi-linearly and all our solutions failed to converge for twists little above 1.1 rad. Overall, we could not obtain twists larger than
${\sim} 1.2$
rad for asymptotically dipolar configurations. This is similar to critical values found in previous numerical work on stationary solutions (Akgün et al. Reference Akgün, Miralles, Pons and Cerdá-Durán2016, and references therein) which typically to not exceed
${\sim} 1.5$
rad, but lower than dynamical studies which reach values of
${\sim} 3$
rad (Parfrey et al. Reference Parfrey, Beloborodov and Hui2013; Mahlmann et al. Reference Mahlmann, Philippov, Mewes, Ripperda, Most and Sironi2023). As an exception, the strongly nonlinear solution presented in § 4.2 reaches a twist of 5.3 rad, but this solution is not asymptotically dipolar.
Supplementary data
The scripts implementing the solution described in this paper and for generating all the figures are available at https://doi.org/10.5281/zenodo.18155842.
Acknowledgements
The author thanks Dr F. Mottez (Observatoire de Paris, CNRS, France) for his advice and support, and Dr K. Uryū (University of the Ryukyus, Okinawa, Japan) for helpful discussions about his 2023 paper. The author would like to thank the anonymous referees for their constructive comments that helped improving this work.
Editor Roger Blandford thanks the referees for their advice in evaluating this article.
Funding
This research received no specific grant from any funding agency, commercial or not-for-profit sectors.
Declaration of interests
The author report no conflict of interest.
Appendix A. Vacuum solutions
Vacuum solutions to the Grad–Shafranov equation are usually expressed using Legendre polynomials (see e.g. Viganò et al. Reference Viganò, Pons and Miralles2011). For the purpose of this work we find it insightful to outline a demonstration of these solutions and explicit the first multipoles.
In vacuum (3.18) becomes
\begin{equation} \forall i \geqslant 1, -l(l+1) F_l^{\textrm {v}} - (1-\mu ^2) {F_l^{\textrm {v}}}'' = 0. \end{equation}
We remark that this could directly have been obtained from the Grad–Shafranov equation, (2.7), with an ansatz of the form
$\mathcal{P} = F_l^{\textrm {v}}(\mu )/r^l$
.
Since we specify boundary conditions at
$\mu = 1$
, we will seek a power-series solution of the form
$F_l^{\textrm { v}} = \sum _{i=0} a_i \bar \mu ^i$
, where
$\bar \mu = 1- \mu$
. Indeed, expansions around
$\mu =0$
do not always converge at
$\mu =1$
. Injecting this series into (A1), one finds after a few manipulations that solutions to the equations follow the sequence defined by
\begin{equation} \forall i\gt 0, a_{i+1} = a_i \frac{i(i-1) - l(l+1)}{2i(i+1)}, \end{equation}
and
$a_0 = 0$
. This sequence is initialised by choosing
$a_1$
. One sees that
$a_{i \gt l+1} = 0$
, such that these solutions are finite polynomials. Since
$a_0 =0$
, these solutions satisfy
$F_l^{\textrm {v}}(\bar \mu = 0) = 0$
, as required. It is less obvious what their behaviour is at
$\bar \mu = 2$
(
$\mu = -1$
). A change of variable
$\bar \mu \rightarrow 1-\mu$
shows that these solutions are even for odd values of
$l$
and odd for even values, meaning that
$F_l^{\textrm {v}}(\mu = \pm 1) = 0$
, and all required boundary conditions are matched. It is immediate to see that
$a_1 = -{F^{\textrm {v}}}'(\mu =1) = \dot F^{\textrm {v}}(\bar \mu =0)$
where the dot denotes differentiation with respect to
$\bar \mu$
.
These solutions form a multipole basis. We make explicit the first orders
\begin{align} F_1^{\textrm {v}} & = a_1 \left(\bar \mu - \frac{1}{2}\bar \mu ^2\right) = \frac{a_1}{2}(1-\mu ^2),\end{align}
\begin{align} F_2^{\textrm {v}} & = a_1\bar \mu \left (1 - \frac{3}{2}\bar \mu +\frac{1}{2}\bar \mu ^2\right ) = \frac{a_1}{2}\mu (1-\mu ^2) , \end{align}
\begin{align} F_3^{\textrm {v}} & = a_1\bar \mu \left (1 - 3\bar \mu +\frac{5}{2}\bar \mu ^2 - \frac{5}{8}\bar \mu ^3\right ) \nonumber \\ & = -\frac{a_1}{8}(1-6\mu ^2 + 5\mu ^4). \end{align}
Using (2.3) with
$\mathcal{P} = F_l/r^l$
, we obtain expressions for the vacuum magnetic fields
\begin{align} \boldsymbol{B}_{1}^{\textrm {v}} & = \frac{B_1}{r^3}\left (\cos \theta \boldsymbol{e}_r + \frac{\sin \theta }{2}\boldsymbol{e}_\theta \right ), \end{align}
\begin{align} \boldsymbol{B}_{2}^{\textrm {v}} & = \frac{B_2}{r^4} \left ((2\cos \theta - 1)\boldsymbol{e}_r - \cos \theta \sin \theta \boldsymbol{e}_\theta \right ),\end{align}
\begin{align}\boldsymbol{B}_{3}^{\textrm {v}} & = \frac{B_3}{r^5} \left (\left (\frac{5}{2}\cos ^4\theta - \frac{3}{2}\cos \theta \right)\boldsymbol{e}_r + \left (\frac{\sin ^3}{8}\theta - \frac{\sin \theta }{2}\right)\boldsymbol{e}_\theta \right ),\end{align}
where
$B_1, B_2, B_3$
are the magnetic strengths at the pole corresponding respectively to
$a_1, a_1/2, a_1$
in the notations of (A3)–(A5).
In the case of imposed even parity across the equator, the quadrupole field becomes a split quadrupole which can be written
\begin{equation} \boldsymbol{B}_{2}^{\textrm {v, split}} = \left (H(\mu ) - H(-\mu )\right ) \boldsymbol{B}_{2}^{\textrm {v}} , \end{equation}
where
$H(x)$
is the Heaviside step function such that
$H(x\geqslant 0) = 1$
and
$H(x\lt 0)=0$
.
Appendix B Example of analytical solution for
$F_3$
We consider the case where
$F_1$
is a vacuum dipole, (A3),
$ F_2 = 0$
and
$c_2$
is the only non-zero coupling constant. In this case, the source term of the
$F_3$
equation in (2.7) is
$ [\alpha A]_5 = 3c_2^2 G_5^{(5)} = 3c_2^2 F_1^5$
as can be seen from (3.15).
We solve the equation for
$\phi _3 = F_3/(3c_2^2 \dot F_1(0)^5)$
where we denote
$\dot F$
the derivative of
$F$
with respect to
$\bar \mu = 1-\mu$
\begin{equation} -12 \phi _3 - \bar \mu (2-\bar \mu ) \ddot \phi _3 = \left (\bar \mu - \frac{1}{2}\bar \mu ^2\right )^5. \end{equation}
In order to get a particular solution
$\phi _3^{\rm (p)}$
we proceed exactly as for the vacuum solutions of § 5, by expanding
$ \phi _3^{\textrm {(p)}} = \sum _{i=0} \alpha _i \bar \mu ^i$
, obtaining
\begin{equation} \sum _{i=0}\bar \mu ^i \left [-2(i-1)i \alpha _{i+1} - \alpha _i \left (l(l+1) - i (i-1)\right )\right ] = \left (\bar \mu - \frac{1}{2}\bar \mu ^2\right )^5, \end{equation}
where here
$l = 3$
.
We find the particular solution
\begin{eqnarray} \phi _3^{\textrm {(p)}} &\ =\ & -\frac{8}{15} \bar \mu ^6 + \frac{88}{105}\bar \mu ^7 - \frac{24}{49}\bar \mu ^8 + \frac{113}{882} \bar \mu ^9 - \frac{17}{1323}\bar \mu ^{10} \nonumber \\ & & - \frac{1}{97020}\bar \mu ^{11} + \sum _{i=12}^{\infty }\alpha _i \bar \mu ^i, \end{eqnarray}
where the polynomial part results from matching the first 11 terms with the source term on the right-hand side of (B1), and the series is generated by using the relation in (A2) (replacing
$a$
with
$\alpha$
) for
$i \geqslant 11$
seeded with
$\alpha _{11} = 1/97020$
. The terms up to
$\bigcirc (\bar \mu ^5)$
do not appear here because they correspond to the vacuum solution in (B1). The full solution is thus
$F_3 = F_3^{\textrm {v}} + 3c_2^2 \dot F_1(0)^5 \phi _3^{\textrm { (p)}}$
where the first term is given by (B1) and the second one by (B3). We have checked that (B3) matches the numerical integration used in the main text within numerical error.
We see that the particular solution of (B3) fulfils the boundary condition
$\phi _3^{\textrm { (p)}}(\bar \mu =0) = 0$
and that its first derivative is equally null
$\dot \phi _3^{\textrm {(p)}}(\bar \mu =0) = 0$
, meaning that the slope at the origin is entirely given by the vacuum solution. More generally, it can be shown by induction that the source terms of the equation of order
$i$
are proportional to
$ \bar \mu ^{l}(1 + \bigcirc (\bar \mu ))$
, with
$l= i+2$
, implying that
$\phi _i^{\textrm {(p)}} = \bigcirc (\bar \mu ^l)$
.
On the other hand, this solution is neither symmetric with respect to
$\mu = 0$
nor does it fulfil the second boundary condition
$\phi _3^{\textrm {(p)}}(\bar \mu =2) = 0$
. Indeed, we find numerically
$\phi _3^{\textrm { (p)}}(\bar \mu =2) \simeq 0.057$
. As a consequence, in order to respect the boundary conditions we construct a solution by forcing a symmetry with respect to the equator. In particular, this implies a discontinuity of
$\dot F_3$
at the equator, which results in a discontinuity of the radial magnetic field supported by a toroidal current sheet (see main text).
The series term in (B3) converges for
$\bar \mu \leq 2$
, however, its derivative is not defined at
$\bar \mu = 2$
as its series diverges. Indeed one can show that at
$\bar \mu = 2$
the series is asymptotically given by
$\alpha _i \bar \mu ^i \sim 1/i^2$
, the derivative of which is
${\sim} 1/i$
and thus diverges. From a practical point of view, we note that the polynomial part of (B3) is accurate to
${\sim} 10^{-5}$
, which we expect to be sufficient for most applications.
Appendix C. Convergence
C.1 Discussion of the
$c_{i\neq 2}$
case
A case of interest that we exhibit in § 4 is when
$c_{i\neq 2} = 0$
,
$F_1 \neq 0$
and
$F_{i\gt 1} =0$
. That is, we have a vacuum dipole which sources higher multipoles through the coupling
$c_2$
. Only functions
$\{F_{2i+1}\}_{i\geqslant 1}$
are non-zero, and the fact that
$F_{2i+1}'=0$
means that all solutions are particular solutions without vacuum terms (see main text and Appendices). We show below that the right-hand side of (3.18) is, up to a factor of order 1,
$[\alpha A ]_{l=i+2} \propto \beta _i = \rho ^i F_1'/2$
for
$i \geqslant 3$
, where
$\rho = 3(F_1'/2)^{4}c_2^2$
. In addition, we show that dividing both sides of each equation in (3.18) by
$F_1'/2$
results in
$F_i/F_1'$
depending only on
$\rho$
.
We call
$\chi _i = F_i/(F_i'/2)$
and get for the first orders
\begin{align} \mathcal{L}_1 \chi _1 & = 0, \end{align}
\begin{align} \mathcal{L}_3 \chi _3 & = \rho \chi _1^5, \end{align}
\begin{align} \mathcal{L}_5 \chi _5 & = \rho n_1 \chi _3 \chi _1^4 = \bigcirc \left ( \rho ^2 \right ),\end{align}
\begin{align} \mathcal{L}_7 \chi _7 & = \rho (n_2\chi _5 \chi _1^4 + n_3\chi _3^2\chi _1^3) = \bigcirc \left ( \rho ^3 \right ),\end{align}
where
$\mathcal{L}_i = -i(i+1) - (1-\mu ^2) \partial _{\mu ^2}^2$
is the linear operator on the left-hand side of (3.18). The parameters
$n_i$
are the number of terms involved that can be calculated from (3.7), but they are irrelevant for our current discussion.
By definition
$\chi _1 \sim 1$
. Indeed it is a vacuum solution, (A3), such that
$\max (|\chi _1|) =1$
. It follows that
$\chi _3 \sim \rho$
,
$\chi _5 \sim \rho ^2$
up to a factor of order 1. The fact that particular solutions of
$\chi _i$
are proportional to the amplitude of the right-hand side is immediately apparent by looking for series solutions as we do in Appendices.
We proceed by induction to show that
$\chi _{2i +1} \sim \rho ^i$
. Initialisation has just been shown, therefore, we assume that this is true for all functions up to order
$2i +1$
. At order
$2i+3$
, we see from (3.7) that the right-hand side is
\begin{equation} \rho G_{l+2}^{(5)} = \sum _{i_1+\ldots +i_5 = l+2; i_j \geqslant 1} \chi _{i_1}\ldots \chi _{i_5}, \end{equation}
with
$l =2i+1 +2$
. We see that the only combinations of
$i_1, i_2, i_3, i_4, i_5$
possible can be obtained by taking those present at order
$l$
in
$G_{l}^{(5)}$
and by replacing a single
$i_j$
by
$i_j+2$
. In these terms the largest
$i_j$
is
$2i-1$
for the term
$\chi _{2i-1} \chi _1^4$
. Since
$\chi _{i_j +2 } \sim \rho \chi _{i_j}$
, these terms are indeed of
$\bigcirc (\rho ^{i+1})$
. Thus the induction is proven.
Based on this argument, we expect that the
$\chi _{2i+1}$
can be bounded by a geometric sequence with common ratio
${\sim} \rho$
. Convergence then corresponds to
$\rho \lesssim 1$
, which means
$|c_2| \lesssim 2.3 F_1'(1)^2$
. Numerically, we observe approximate geometric convergence up to
$|c_2| \lesssim 6 F_1'(1)^2$
. It is to be noted, that we have here neglected the increasing number of terms on the right-hand side. However, one can see by looking at (3.7) that, without detailed calculation, this number can be bounded by
$l^4$
. Thus the number of terms increases polynomially, which does not affect geometric convergence.
The series that actually must converge is
$\chi _i/r^l$
. It follows that if
$\chi _{2i+1} = \bigcirc (\rho ^i)$
then
$\chi _i/r^l= \bigcirc ((\rho /r)^i)$
. As a result, if the series converges at the surface, for
$r=1$
, then it converges for any
$r \geqslant 1$
. If
$\rho = \rho _{\max }$
, the largest value for which convergence is allowed at
$r=1$
, then the series is expected to diverge for
$r \lt 1$
. This is the behaviour we numerically observe.
Appendix D. Parity
We want to show that the source terms of even-order functions
$F_{2i}$
are odd and even for odd-order functions
$F_{2i-1}$
in the case
$c_{i\neq 2} =0$
. The source term of the equation for
$F_i$
is
$3c_2^2 G_{i+2}^{(5)}$
for
$i \geqslant 3$
as can be seen from (3.12). For
$i = 3$
,
$G_{5}^{(5)} = F_1^5$
is even, as desired. For
$i=4$
,
$G_{6}^{(5)} \propto F_2 F_1^4$
, which is odd, equally as desired. Equation (3.7) shows that
$G_{l}^{(5)} = \sum _{i_1, .., i_5} F_{i_1}..F_{i_5}$
with
$i_1+..+i_5=l$
and
$l=i+2$
. From (3.9), we know that
$G_{l}^{(5)}$
depends only on function
$F_{j\lt i}$
, such that we can prove by induction the result. Given
$i\gt 2$
, we assume that parity properties are fulfilled for all
$j\lt i$
. Then, for any
$l \gt 4$
, if
$l$
is even then an even number of odd indices is required which means that the number of even indices is odd. It follows that
$G_{l}^{(5)}$
is odd for even
$l$
, as desired. Conversely, for odd
$l$
we need an odd number of odd indices and an even number of even indices, such that
$G_{l}^{(5)}$
is even, as desired. Thus, we have proven the result.
Appendix E. Magnetic helicity
One can check that
\begin{equation} \mathcal{A} = \left (-r\int _{0}^{\theta }B_\varphi \mathrm{d}\theta ',0,\frac{\mathcal{P}}{r\sin \theta }\right )_{(r,\theta , \varphi )} \end{equation}
forms a suitable vector potential for the solution presented in § 3.
Helicity is defined by
$H = \int \mathcal{A} \boldsymbol{\cdot}\boldsymbol{B} \;\mathrm{d} V$
. The radial term of the scalar product gives
\begin{equation} -\int \frac{\partial _\theta \mathcal{P} }{r^2 \sin \theta } r\int _{0}^{\theta }B_\varphi \mathrm{d}\theta ' \;\mathrm{d} V = \int \frac{\mathcal{P} B_\varphi }{r \sin \theta } \;\mathrm{d} V, \end{equation}
where the right-hand side is obtained by integration by parts in
$\theta$
. Adding up the other term of the scalar product we obtain
\begin{equation} H = 2\int \frac{\mathcal{P} B_\varphi }{r \sin \theta }\; \mathrm{d} V. \end{equation}
Open access
