Non-rigid precession of magnetic stars

Lander, S. K.; Jones, D. I.

doi:10.1093/mnras/stx349

Abstract

Stars are, generically, rotating and magnetized objects with a misalignment between their magnetic and rotation axes. Since a magnetic field induces a permanent distortion to its host, it provides effective rigidity even to a fluid star, leading to bulk stellar motion that resembles free precession. This bulk motion is, however, accompanied by induced interior velocity and magnetic field perturbations, which are oscillatory on the precession time-scale. Extending previous work, we show that these quantities are described by a set of second-order perturbation equations featuring cross-terms scaling with the product of the magnetic and centrifugal distortions to the star. For the case of a background toroidal field, we reduce these to a set of differential equations in radial functions, and find a method for their solution. The resulting magnetic field and velocity perturbations show complex multipolar structure and are strongest towards the centre of the star.

stars: interiors, stars: magnetic field, stars: oscillations, stars: rotation

1 INTRODUCTION

Two of the most important pieces of stellar physics are their magnetic fields and rotation. They have a ubiquity across many classes of star that are otherwise governed by very different interior and exterior physics: main-sequence stars, white dwarfs, neutron stars or combinations of these in binary systems. A fundamental problem is the rich variety of ways in which rotation and magnetism interact, and the effects on stellar properties and evolution.

In this paper, we focus on one of the simplest situations within the broad class of magneto-rotational stellar phenomena: the dynamics of a rigidly rotating star with a frozen-in dynamically stable magnetic field,¹ symmetric about some axis. In the case where the inclination angle χ between the magnetic and rotational axes is zero, i.e. the two axes coincide, the star is stationary. This widely studied situation has the advantage of being simple, but the disadvantage of having limited applicability to astrophysics: observations indicate that stars have a wide distribution of inclination angles (Schmidt & Norsworthy 1991; Tauris & Manchester 1998; Donati & Landstreet 2009; Rookyard, Weltevrede & Johnston 2015). A non-zero χ is essential, for example, to explain the pulsed emission seen from many neutron stars.

The dynamics of a star with non-zero χ – an ‘oblique rotator’ – is still not well understood. The classic work on the topic is by Mestel and collaborators (Mestel & Takhar 1972; Mestel et al. 1981), as discussed in Section 2. The essential idea, put forward by Spitzer (1958), is that the stellar distortion associated with the magnetic field causes it to undergo bulk motion like that of rigid-body free precession. Mestel & Takhar (1972) argue that internally this bulk motion has to be supported by time-varying and non-axisymmetric velocity and magnetic fields. Formulated in generality, this problem is intractable by analytic methods; accordingly, Mestel and collaborators made some major simplifications to produce solutions for their early studies. Various subsequent papers have adopted their ideas and solutions (see e.g. Wasserman 2003; Dall'Osso, Shore & Stella 2009; Lasky & Glampedakis 2016), but none have attempted to extend them.

The goal of the present paper is to build on the pioneering ideas of Mestel and collaborators, but to formulate them in a more rigorous and general way. The practical difference is that we are forced to perform perturbation analysis to a higher order than their work, with an accompanying increase in the complexity of the algebra. This paper is entirely concerned with the theory of this problem, and the ultimate desideratum is to obtain solutions for the perturbed velocity and magnetic fields of a fluid star with a non-zero χ. In order to keep the majority of the calculation analytically tractable, we make two key simplifications: the background magnetic field is assumed to be purely toroidal, and we employ a polytropic equation of state. For numerical reasons, we have also been forced to neglect the outermost layer of the star in our solutions. In future work, we will study the implication of our results for observed phenomena, like the distribution and time-evolution of inclination angles amongst stars, and the damping of precession.

This paper is structured as follows. In Section 2, we summarize existing models of the internal dynamics of oblique rotators, providing a critique of their applicability and suggesting how to extend them. In Section 3, we formulate the oblique-rotator problem using a perturbation scheme in two small quantities, the centrifugal and magnetic distortions, ε_α and ε_B, respectively; we show that the velocity and magnetic-field perturbations are of the same perturbative order as cross-terms involving the interaction of the stellar rotation and background magnetic field, and so finding these entails the solution of the order-ε_αε_B perturbation equations. In Section 4, we systematically present solutions to all the lower order problems: the zeroth-order, order-ε_α and order-ε_B equations. Section 5 then addresses the problem of solving the second-order equations, through a series of stages. In Section 5.2, we reduce the full system of governing vector equations to poloidal and toroidal scalar equations. Following this, we perform a spherical-harmonic decomposition of the perturbed magnetic field (Section 5.3), yielding an infinite system of differential equations (DEs) in radial functions associated with each spherical harmonic. Superficially these appear to be ordinary differential equations (ODEs), but are, in fact, a more complicated system of differential-algebraic equations (DAEs), which cannot be solved with ODE methods. None the less, we find and present a method to solve the equations for the magnetic radial functions, when truncating at the fourth multipole (Section 5.4). In Section 6, we find closed-form expressions for the perturbed velocity field in terms of the magnetic radial functions. After this, we present results for the magnetic and velocity perturbations in Section 7. Finally, we discuss our findings and their applicability to different classes of star in Section 8, and summarize in Section 9.

2 GENERATION OF PRECESSIONAL MOTION IN A MAGNETIzED FLUID STAR

Precession is a rigid-body effect, and superficially one would not expect a fluid star (or the fluid region of a neutron star) to be able to sustain such motion. In fact, we will adopt a stellar model that is isothermal and non-convective, so we would expect only circular motion of fluid elements about the star's rotation axis. However, as argued by Spitzer (1958), a magnetic field |${\boldsymbol B}$| can provide effective rigidity to a star, since it provides a distortion (symmetric about the magnetic-field axis) away from the star's spherically symmetric hydrostatic state (Chandrasekhar & Fermi 1953). If the star is now rotated about a different axis from the magnetic one, it will undergo a secondary rotation² about the magnetic axis, so that the bulk motion of the star resembles precession.

In this section, we begin by recapitulating the work of Mestel & Takhar (1972), who put this conceptual picture on a rigorous footing and derived the form of the secondary rotation. In doing so, we introduce notation that differs in many cases from theirs. We then describe their arguments as to why the motion is not exactly described by rigid-body free precession, and the approaches taken to finding the non-rigid internal dynamics that maintain the bulk precessional motion in both the original paper of Mestel & Takhar (1972) and a rather later follow-up study (Mestel et al. 1981). We conclude the discussion with a critique of these approaches, and describe our strategy for obtaining a more detailed solution that we believe includes the key physics missing from their work. Finally, we discuss the ordering of characteristic time-scales required for non-rigid precession to occur, and compare the various classes of magnetic star to which our analysis may apply.

We model a star as a fluid ball rotating uniformly at angular frequency α about an axis |${\boldsymbol e}_z^{(\alpha )}$|⁠, with a frozen-in magnetic field symmetric about some axis |${\boldsymbol e}_z^{(B)}$|⁠. The axis |${\boldsymbol e}_z^{(B)}$| is misaligned by some angle χ from the primary rotation axis |${\boldsymbol e}_z^{(\alpha )}$|⁠. At different points in our calculation, we will need to refer to both the symmetry axis of primary rotation and that of the magnetic field, and so we form right-handed triads |$({{\boldsymbol e}_x,{\boldsymbol e}_y^{(B)},{\boldsymbol e}_z^{(B)}})$| and |$({\boldsymbol e}_x,{\boldsymbol e}_y^{(\alpha )},{\boldsymbol e}_z^{(\alpha )})$| associated with these magnetic and rotational axes, with the |${\boldsymbol e}_x$| vector being instantaneously common to both (see the left-hand side of Fig. 1). In addition, we denote the spherical polar coordinate system referred to the |${\boldsymbol e}_z^{(B)}$|-triad by (r, θ, ϕ); for the rest of this section, we shall work exclusively in this coordinate system (right-hand side of Fig. 1).

$Left-hand side: the magnetic and rotational triads; we assume ${\boldsymbol e}_y^{(B)},{\boldsymbol e}_z^{(B)},{\boldsymbol e}_y^{(\alpha )}$ and ${\boldsymbol e}_z^{(\alpha )}$ are coplanar. Right-hand side: the ${\boldsymbol e}_z^{(B)}$-triad and its spherical polar coordinate system.$

Figure 1.

Left-hand side: the magnetic and rotational triads; we assume |${\boldsymbol e}_y^{(B)},{\boldsymbol e}_z^{(B)},{\boldsymbol e}_y^{(\alpha )}$| and |${\boldsymbol e}_z^{(\alpha )}$| are coplanar. Right-hand side: the |${\boldsymbol e}_z^{(B)}$|-triad and its spherical polar coordinate system.

Open in new tab Download slide

2.1 The bulk motion of the star

Unless otherwise specified, the rest of this section follows the reasoning of Mestel & Takhar (1972), but sometimes explained in different terms and with different notation. We first recall that a static, unmagnetized star composed of homogeneous fluid would have a spherically symmetric density field ρ₀(r). Including rotation alone adds on a small extra term³ δρ_α(r, θ, ϕ) corresponding to a centrifugal bulge; similarly, the density distribution of a non-rotating magnetized fluid ball could be written as ρ(r, θ) = ρ₀(r) + δρ_B(r, θ) to take account of magnetic distortions δρ_B. Hence, for a rotating, magnetized star, we may write the density of an element at the point (r, θ, ϕ), and at some instant in time, as

\begin{equation} \rho (r,\theta ,\phi ) = \rho _0(r)+\delta \rho _{\rm B}(r,\theta )+\delta \rho _\alpha (r,\theta ,\phi ), \end{equation}

(1)

where – for now – we have neglected cross-terms ∼δρ_αδρ_B for being higher order than the other density components.

The density field of a star rotating with angular velocity |$\alpha {\boldsymbol e}_z^{(\alpha )}$| has the angular momentum vector

\begin{equation} {{\boldsymbol J}}_\alpha = \int \rho {{\boldsymbol r}}\times (\alpha {\boldsymbol e}_z^{(\alpha )}\times {{\boldsymbol r}}) \ \mathrm{d}V. \end{equation}

(2)

However, this alone does not give an invariant angular momentum orientated along the |${\boldsymbol e}_z^{(\alpha )}$| direction, as the |${\boldsymbol e}_y^{(\alpha )}$|-component of (2) is non-zero:

\begin{equation} {\boldsymbol J}_\alpha \cdot {\boldsymbol e}_y^{(\alpha )} = -\alpha \int ({\boldsymbol e}_y^{(\alpha )}\cdot {\boldsymbol r})({\boldsymbol e}_z^{(\alpha )}\cdot {\boldsymbol r})\delta \rho _{\rm B} \ \mathrm{d}V\!, \end{equation}

(3)

where the contributions from ρ₀ and δρ_α vanish by symmetry. To yield an invariant angular momentum, we require an additional rotation ω about the magnetic axis |${\boldsymbol e}_z^{(B)}$| with an associated angular momentum |${\boldsymbol J}_{B}$|⁠, such that |$({\boldsymbol J}_\alpha +{\boldsymbol J}_B)\cdot {\boldsymbol e}_y^{(\alpha )} = 0$|⁠, i.e.

\begin{equation} {\boldsymbol J}_{B}\cdot {\boldsymbol e}_y^{(\alpha )} = \int \rho {\boldsymbol r}\times (\omega {\boldsymbol e}_z^{(B)}\times {\boldsymbol r})\cdot {\boldsymbol e}_y^{(\alpha )} \ \mathrm{d}V = \alpha \int ({\boldsymbol e}_y^{(\alpha )}\cdot {\boldsymbol r})({\boldsymbol e}_z^{(\alpha )}\cdot {\boldsymbol r})\delta \rho _{\rm B} \ \mathrm{d}V. \end{equation}

(4)

Working in spherical polar coordinates with |${\boldsymbol r}= r(\sin \theta \cos \phi \ {\boldsymbol e}_x + \sin \theta \sin \phi \ {\boldsymbol e}_y^{(B)} + \cos \theta \ {\boldsymbol e}_z^{(B)})$|⁠, and writing |${\boldsymbol e}_y^{(\alpha )} = \cos \chi \ {\boldsymbol e}_y^{(B)}+\sin \chi \ {\boldsymbol e}_z^{(B)}$|⁠, |${\boldsymbol e}_z^{(\alpha )} = -\sin \chi \ {\boldsymbol e}_y^{(B)}+\cos \chi \ {\boldsymbol e}_z^{(B)}$| and dV = r²sin θ drdθdϕ, we now evaluate the integral (3) in the |$({\boldsymbol e}_x,{\boldsymbol e}_y^{(B)},{\boldsymbol e}_z^{(B)})$| triad to give

\begin{equation} {\boldsymbol J}_\alpha \cdot {\boldsymbol e}_y^{(\alpha )} = 2\pi \alpha \sin \chi \cos \chi \int \!\!\!\int \delta \rho _{\rm B} P^0_2 (\cos \theta ) r^4\sin \theta \ \mathrm{d}r \mathrm{d}\theta , \end{equation}

(5)

where |$P_l^m(\cos \theta )$| denotes an associated Legendre polynomial of angular index l and azimuthal index m (which, in this equation, is the ordinary Legendre polynomial |$P_2\equiv P_2^0 = \frac{1}{2}(3\cos ^2\theta - 1)$|⁠). We evaluate the |${\boldsymbol e}_y^{(\alpha )}$|-component of |${\boldsymbol J}_{B}$| in a similar fashion to give

\begin{equation*} {\boldsymbol J}_{B}\cdot {\boldsymbol e}_y^{(\alpha )} = \int \rho {\boldsymbol r}\times (\omega {\boldsymbol e}_z^{(B)}\times {\boldsymbol r})\cdot {\boldsymbol e}_y^{(\alpha )}\ \mathrm{d}V\nonumber = I_0 \omega \sin \chi ,\end{equation*}

where |$I_0\equiv \frac{8\pi }{3}\int \rho _0 r^4\ \mathrm{d}r$| is the moment of inertia of the spherically symmetric density field ρ₀; here the two density perturbations are regarded as negligible parts of ρ in comparison with ρ₀. We now use equations (5) and (6), and the requirement |$({\boldsymbol J}_\alpha +{\boldsymbol J}_{B})\cdot {\boldsymbol e}_y^{(\alpha )} = 0$| to find the precession frequency:

\begin{equation} \omega = -\frac{2\pi \alpha \cos \chi }{I_0} \int \!\!\!\int \delta \rho _{\rm B} P_2(\cos \theta )r^4\sin \theta \ \mathrm{d}r\mathrm{d}\theta . \end{equation}

(6)

Rewriting this expression using the xx and zz components of the moment-of-inertia tensor I_jk ≡ ∫ρ(r²δ_jk − x_jx_k) dV referred to the |${\boldsymbol e}_z^{(B)}$|-triad yields

\begin{equation} \omega = \alpha \cos \chi \ \frac{I_{zz}-I_{xx}}{I_0} \equiv \alpha \epsilon _B\cos \chi , \end{equation}

(7)

the usual rigid-body result (Landau & Lifshitz 1976); this is not surprising since we have not yet put any fluid (non-rigid) physics into the calculation.

2.2 Deviation from rigid-body precession due to internal fluid dynamics

The result at the end of the previous section suggests that the macroscopic dynamics of a rotating magnetized fluid body should resemble free precession; however, the fluid is clearly not a rigid body. This presents a question as to what degree the magnetized fluid can be regarded as rigid and hence how similar the motion of a magnetized fluid is to conventional rigid-body precession. Mestel and collaborators sought to answer this by considering the microscopic dynamics: the effect of precession on individual fluid elements and the induced non-rigid velocity field.

Since we will need to distinguish between different frames of reference here, we define for brevity the ‘α-frame’ to be the one comoving with the star's primary rotation (at frequency α) and the ‘ω-frame’ to be the co-precessing frame — i.e. the rigid-body precession frame characterized by the superimposed rotations α and ω.

We wish to investigate the deviation of a rotating magnetized fluid star from free precession. If the fluid moved rigidly, then each fluid element would be stationary as viewed by the co-precessing observer in the ω-frame. Since we do not expect exact rigid-body precession here, let us define the Lagrangian displacement |$\xi$| to be the change in position of a fluid element in the co-precessing frame, with its time derivative |$\dot{\xi }$| giving the velocity of the element as viewed from the ω-frame. Following Mestel & Takhar (1972), we will sometimes refer to this higher order velocity field as ‘the |$\xi$|-motions’. On viewing the star in the inertial frame, we will then see that the motion of a fluid element |${\boldsymbol V}_{\mathrm{inertial}}$| is a vector sum of three characteristic velocities: the primary stellar rotation α about the rotation axis, the slower nutation ω about the magnetic axis and the extra velocity field |$\dot{\xi }$|⁠:

\begin{equation} {\boldsymbol V}_{\mathrm{inertial}} = \alpha {\boldsymbol e}_z^{(\alpha )}\times {\boldsymbol r}+ \omega {\boldsymbol e}_z^{(B)}\times {\boldsymbol r}+ \dot{\xi }. \end{equation}

(8)

An intuitive explanation for the existence of this extra velocity field, presented in Mestel & Takhar (1972), is as follows. Consider the motion of a fluid element in the α-frame; see the left-hand panel of Fig. 2. In the unmagnetized case, the element undergoes only the primary rotation α and so is stationary in the corotating frame. From Section 2.1, we anticipate that the addition of a misaligned magnetic field will cause the star to precess, and a fluid element in the α-frame will therefore undergo a slow secondary rotation (with frequency ω) about the magnetic axis. In doing so, however, the fluid element will be moved through regions of differing density. Since the background density ρ₀ is spherical and the magnetic distortion δρ_B is symmetric about its axis, the density difference will be entirely due to the centrifugal bulge δρ_α. If we move to the ω-frame (Fig. 2, right-hand panel), on a macroscopic level, the entire centrifugal bulge would rotate at a rate ω about the |${\boldsymbol e}_z^{(B)}$|-axis.

$Two views of a precessing fluid star. Left-hand side: the dynamics in the α-frame, i.e. the frame rigidly rotating with rate α; right-hand side: the dynamics in the co-precessing ω-frame. We show only the centrifugal contribution to the distortion here, so that the stellar surface (the solid black line) and the isopycnic surfaces (the dashed black lines) are spheroidal about the ${\boldsymbol e}_z^{(\alpha )}$ axis. Without a magnetic field, a fluid element would be stationary in the α-frame; however, the magnetic field induces a slow precessional motion, superimposed on the normal stellar rotation. This motion will cause a fluid element (the filled red circle) in the α-frame to rotate about the magnetic axis ${\boldsymbol e}_z^{(B)}$ with period 2π/ω. A higher order velocity field $\dot{\xi }$ is needed to prevent the fluid element from travelling through regions of varying density (crossing the density contours as shown). An alternative viewpoint of the process can be obtained from the ω-frame, in which the whole centrifugal bulge would be rotated at rate ω around the ${\boldsymbol e}_z^{(B)}$ axis were it not for the $\dot{\xi }$ field.$

Figure 2.

Two views of a precessing fluid star. Left-hand side: the dynamics in the α-frame, i.e. the frame rigidly rotating with rate α; right-hand side: the dynamics in the co-precessing ω-frame. We show only the centrifugal contribution to the distortion here, so that the stellar surface (the solid black line) and the isopycnic surfaces (the dashed black lines) are spheroidal about the |${\boldsymbol e}_z^{(\alpha )}$| axis. Without a magnetic field, a fluid element would be stationary in the α-frame; however, the magnetic field induces a slow precessional motion, superimposed on the normal stellar rotation. This motion will cause a fluid element (the filled red circle) in the α-frame to rotate about the magnetic axis |${\boldsymbol e}_z^{(B)}$| with period 2π/ω. A higher order velocity field |$\dot{\xi }$| is needed to prevent the fluid element from travelling through regions of varying density (crossing the density contours as shown). An alternative viewpoint of the process can be obtained from the ω-frame, in which the whole centrifugal bulge would be rotated at rate ω around the |${\boldsymbol e}_z^{(B)}$| axis were it not for the |$\dot{\xi }$| field.

Open in new tab Download slide

At this point, we may argue for the existence of a field of non-rigid |$\xi$|-motions. In the α-frame, fluid elements undergoing rigid-body free precession (i.e. with |$\dot{\xi }= {0}$|⁠) would be forced to sustain large-density variations over one precession period. For this reason, there will be a restoring force on each fluid element that acts to return it to its original density; hence, the precessional motion described above cannot be completely rigid. Equivalently, on a macroscopic level and in the ω frame, one would expect the global effect of the internal |$\xi$|-motions to be a restoration of the star to an instantaneous stationary equilibrium.

At this point. Mestel & Takhar (1972) used this intuitive picture to argue for the functional form of δρ_α, the piece of the centrifugal bulge that is time-varying in the ω-frame, by simply performing a rotation of the coordinates with respect to the rotation axis.⁴ They described a double-perturbation formalism in the two small parameters ε_α and ε_B – the centrifugal and magnetic distortions to the star – respectively. Crucially, they argued that one can complete the calculation merely by going to first order in each of the perturbations – i.e. using the zeroth-order, |$\mathcal {O}(\epsilon _\alpha )$| and |$\mathcal {O}(\epsilon _B)$| equations, whilst neglecting higher order cross-terms whose scaling is ε_αε_B or higher. However, at this level, one only has a single equation to fix the three spatial components of |$\xi$|⁠: the continuity equation, which in time-integrated form is |$\delta \rho _\alpha = -\nabla \cdot (\rho _0\xi )$|⁠. To obtain a second condition, they specialized to motions with |$\nabla \cdot \xi = 0$|⁠, effectively appealing to an additional buoyancy force associated with stratification; this clearly already reduces the generality of their results. Despite this, one extra relation is still required to close the system. In the first of the two papers from Mestel and collaborators (Mestel & Takhar 1972), they chose ξ_ϕ = 0, whilst in the second one, they demanded minimization of the kinetic energy of |$\xi$| (Mestel et al. 1981).

2.3 The need for second-order perturbations

Although Mestel and collaborators claim that their results are likely to be qualitatively correct, they clearly used ad hoc assumptions to close their system of equations. This should not have been necessary, since the original problem of a rotating fluid star was perfectly well defined. In particular, it is odd that the magnetic field – by dint of which the precession was possible to start with – does not enter directly at all. Instead, in our approach, we follow the original steps of Mestel & Takhar (1972) and set up a perturbation scheme in the small parameters ε_α and ε_B, but then rigorously write out all the resulting systems of equations. By doing so, we show that the non-rigid response of the fluid to the precessional motion – encoded in the velocity field |$\dot{\xi }$| and a perturbed magnetic field |$\delta {\boldsymbol B}$| – only enters at order ε_αε_B. The resulting hierarchy of equations – at the zeroth, first and second perturbative order – is unsurprisingly more complex than that considered in previous work, but we believe that these equations contain the minimal information required to get a reliable solution.

2.4 Key parameters for precessing magnetic stars

Before we begin our detailed modelling of non-rigid precession, let us pause to define the ordering of time-scales necessary for this to occur. First, by inverting equation (7), we immediately see that the precession period T_ω = 2π/ω will always be significantly longer than the primary-rotation period T_α = 2π/α:

\begin{equation} T_\omega = \frac{T_\alpha }{\epsilon _B\cos \chi } \gg T_\alpha . \end{equation}

(9)

Note that the time-scale T_ω must also be the oscillation period of the star's non-rigid response to precession, encoded in the perturbed fluid velocity |$\dot{\xi }$| and magnetic field |$\delta {\boldsymbol B}$|⁠. In order for the star to precess in a manner analogous to a rigid body, it is natural to assume that the characteristic magnetic-mode crossing time-scale should be short compared with the precession period. There is some subtlety in defining a suitable magnetic-mode time-scale, however: in non-rotating magnetic stars, the relevant time-scale is the Alfvén-mode crossing time-scale τ_A, but rotation splits each Alfvén mode into a pair of co- and counter-rotating magneto-inertial modes (Lander, Jones & Passamonti 2010). If rotation is relatively unimportant, in the sense that τ_A/T_α is small, then these modes are virtually Alfvén modes; if τ_A/T_α is large, then one branch effectively becomes a pure inertial mode, whereas the other branch is virtually stationary as viewed in the rotating frame, and so its oscillation period is far longer than τ_A. This agrees with the result of Levin & D'Angelo (2004) that rotation increases magnetic-mode time-scales.

For simplicity, let us assume that the magneto-inertial crossing time-scale is approximately that of a pure Alfvén wave τ_A, in which case the criterion for precession is

\begin{equation} \tau _A\sim \frac{R_*}{B/\sqrt{\rho }}\sim \frac{\mathcal {M}^{1/2}}{R_*^{1/2}B} \le T_\omega , \end{equation}

(10)

where R_* is the radius and |$\mathcal {M}$| the mass of the star. We will see in the table below that this criterion is very easily satisfied, so the inequality will also hold even for more rapidly rotating stars where τ_A is increased by an order of magnitude or more.

A number of different classes of star are thought to harbour long-lived magnetic fields and to have significant rotation rates. These are all candidates for undergoing the non-rigid precession mechanism upon which this paper is focused. To orientate the reader, Table 1 gives some typical parameters for these stars, including estimates for τ_A and T_ω. These are made by assuming that the volume-averaged magnetic field strength is equal to the surface value B_surf as inferred from observations. In reality, the interior field could be considerably stronger (especially if it is dominated by a toroidal component), so the absolute values we report for τ_A and T_ω should be taken with caution – however, the ordering τ_A ≪ T_ω will hold regardless. In addition, our estimate of T_ω for a typical Ap/Bp star is in broad agreement with the results of the calculations of Mestel et al. (1981) and Nittmann & Wood (1981).

Table 1.

Open in new tab

Key parameters and time-scales for different classes of potentially precessing magnetic stars.

Class of star	Mass \|$\mathcal {M}$\|/g	Radius R_*/cm	B_surf/G	T_α	Estimated τ_A	Estimated T_ω/yr
O star¹	7 × 10³⁴	7 × 10¹¹	10³	15 d	10 yr	5 × 10⁷
Ap/Bp star²	4 × 10³³	10¹¹	10⁴	1 d	0.6 yr	3 × 10⁵
Magnetic white dwarf³	2 × 10³³	10⁹	10⁸	1 h	0.4 h	3000
Pulsar⁴	3 × 10³³	10⁶	10¹²	0.1 s	50 s	2000
Magnetar⁴	3 × 10³³	10⁶	10¹⁵	10 s	0.05 s	0.2

Class of star	Mass \|$\mathcal {M}$\|/g	Radius R_*/cm	B_surf/G	T_α	Estimated τ_A	Estimated T_ω/yr
O star¹	7 × 10³⁴	7 × 10¹¹	10³	15 d	10 yr	5 × 10⁷
Ap/Bp star²	4 × 10³³	10¹¹	10⁴	1 d	0.6 yr	3 × 10⁵
Magnetic white dwarf³	2 × 10³³	10⁹	10⁸	1 h	0.4 h	3000
Pulsar⁴	3 × 10³³	10⁶	10¹²	0.1 s	50 s	2000
Magnetar⁴	3 × 10³³	10⁶	10¹⁵	10 s	0.05 s	0.2

References: ¹Donati et al. (2002)

²Landstreet (1992)

³Ferrario, de Martino & Gänsicke (2015)

⁴Turolla, Zane & Watts (2015)

Table 1.

Open in new tab

Key parameters and time-scales for different classes of potentially precessing magnetic stars.

Class of star	Mass \|$\mathcal {M}$\|/g	Radius R_*/cm	B_surf/G	T_α	Estimated τ_A	Estimated T_ω/yr
O star¹	7 × 10³⁴	7 × 10¹¹	10³	15 d	10 yr	5 × 10⁷
Ap/Bp star²	4 × 10³³	10¹¹	10⁴	1 d	0.6 yr	3 × 10⁵
Magnetic white dwarf³	2 × 10³³	10⁹	10⁸	1 h	0.4 h	3000
Pulsar⁴	3 × 10³³	10⁶	10¹²	0.1 s	50 s	2000
Magnetar⁴	3 × 10³³	10⁶	10¹⁵	10 s	0.05 s	0.2

Class of star	Mass \|$\mathcal {M}$\|/g	Radius R_*/cm	B_surf/G	T_α	Estimated τ_A	Estimated T_ω/yr
O star¹	7 × 10³⁴	7 × 10¹¹	10³	15 d	10 yr	5 × 10⁷
Ap/Bp star²	4 × 10³³	10¹¹	10⁴	1 d	0.6 yr	3 × 10⁵
Magnetic white dwarf³	2 × 10³³	10⁹	10⁸	1 h	0.4 h	3000
Pulsar⁴	3 × 10³³	10⁶	10¹²	0.1 s	50 s	2000
Magnetar⁴	3 × 10³³	10⁶	10¹⁵	10 s	0.05 s	0.2

References: ¹Donati et al. (2002)

²Landstreet (1992)

³Ferrario, de Martino & Gänsicke (2015)

⁴Turolla, Zane & Watts (2015)

3 FLUID PRECESSION AS A SECOND-ORDER PERTURBATION PROBLEM

3.1 Full equations

For our non-axisymmetric stellar model, we need to consider a very general form of the standard equations of motion. First, we have the Euler equation, referred to a frame that rotates non-uniformly with angular velocity |${\boldsymbol \Omega }(t)$|⁠:

\begin{equation} \frac{\mathrm{\partial} {\boldsymbol v}}{\mathrm{\partial} t} + ({\boldsymbol v} \cdot \nabla ) {\boldsymbol v} + 2 {\boldsymbol \Omega } \times {\boldsymbol v} + \frac{\mathrm{d}{\boldsymbol {\Omega }}}{\mathrm{d}t} \times {\boldsymbol r} + {\boldsymbol \Omega } \times ({\boldsymbol \Omega } \times {\boldsymbol r}) =- \nabla H - \nabla \Phi + \frac{1}{4\pi \rho } (\nabla \times {\boldsymbol B}) \times {\boldsymbol B}, \end{equation}

(11)

where Φ is the gravitational potential and |${\boldsymbol v}$| the fluid velocity (at this stage, we have not chosen a specific rotating frame, so this is not yet the same thing as |$\dot{\xi }$|⁠). Note the presence of the Euler force|$\frac{\mathrm{d}{\boldsymbol \Omega }}{\mathrm{d}{t}}\times {\boldsymbol r}$| (present because we are allowing for non-uniform rotation), in addition to the more familiar Coriolis and centrifugal force terms |$2\boldsymbol \Omega \times {\boldsymbol v}$| and |$\boldsymbol \Omega \times (\boldsymbol \Omega \times {\boldsymbol r})$|⁠. As we will deal only with barotropic fluids, with equations of state of the form P = P(ρ), we have chosen to use enthalpy H in preference to the pressure P; the two quantities are related by

\begin{equation} H({\boldsymbol r}) \equiv \int _{\hat{P} = 0}^{\hat{P} = P({\boldsymbol r})} \frac{\mathrm{d}\hat{P}}{\rho (\hat{P})}. \end{equation}

(12)

This choice of variable will make the perturbation equations simpler. Together with the Euler equations, we also have the continuity, Poisson and induction equations, the equation of state and the solenoidal constraint on the magnetic field:

\begin{eqnarray} \frac{\mathrm{\partial} \rho }{\mathrm{\partial} t} = -\nabla \cdot (\rho {\boldsymbol v}), \end{eqnarray}

(13)

\begin{eqnarray} \nabla ^2\Phi = 4\pi G \rho , \end{eqnarray}

(14)

\begin{eqnarray} \frac{\mathrm{\partial} {\boldsymbol B}}{\mathrm{\partial} t} = \nabla \times ({\boldsymbol v} \times {\boldsymbol B}), \end{eqnarray}

(15)

\begin{eqnarray} H = H(\rho ), \end{eqnarray}

(16)

\begin{eqnarray} \nabla \cdot {\boldsymbol B}= 0. \end{eqnarray}

(17)

Using the logic of the previous section, we argue that the motion of our star is close to that of a freely precessing rigid body. To some approximation, then, the motion must then consist of two superimposed rotations, one at rate α about |${\boldsymbol e}_z^{(\alpha )}$|⁠, and the other at a much slower rate ω = αε_Bcos χ about the magnetic axis |${\boldsymbol e}_z^{(B)}$|⁠, with |${\boldsymbol e}_z^{(B)}$| tracing out a cone of half-angle χ about |${\boldsymbol e}_z^{(\alpha )}$| at a rate α. However, from the point of view of an observer rotating with the body, the centrifugal bulge of size ε_α rotates about |${\boldsymbol e}_z^{(B)}$|⁠, again in a cone of half-angle χ, at the slow precession frequency ω; recall Fig. 2. It is this slow density wave that produces an Eulerian density perturbation δρ_α(t) (see below for the precise meaning), which in turn induces the |$\xi$|-motions that encode the non-rigid response.

3.2 Perturbative scheme

In order to derive sets of perturbation equations, we first need to choose which frame to work in. There are three natural options:

The inertial frame, |$\boldsymbol \Omega = 0$|⁠. This is conceptually simple, but in this frame, both the mass and magnetic fields are time-varying.
The α-frame, |${\boldsymbol \Omega } = \alpha {\boldsymbol e}_z^{(\alpha )}$|⁠. In this frame, both the mass and magnetic fields are stationary to the lowest order, but individual fluid elements move in large circles about |${\boldsymbol e}_z^{(B)}$|⁠.
The ω-frame, |${\boldsymbol \Omega } = \alpha {\boldsymbol e}_z^{(\alpha )} + \omega {\boldsymbol e}_z^{(B)}$|⁠. In this frame, the magnetic field is stationary to the lowest order, but the |${\boldsymbol e}_z^{(\alpha )}$| axis moves about |${\boldsymbol e}_z^{(B)}$| with a period 2π/ω. Consequently, the mass field is also time-varying on the ω time-scale, inducing the |$\xi$|-motions.

We will work in the ω-frame. This is the frame that is closest to the ‘body frame’ to which the equations of rigid bodies are conventionally referred. One key advantage of this is that the only velocity component in the Euler equation is |$\dot{\xi }$|⁠. In the ω-frame, |$\boldsymbol \Omega$| is the actual angular velocity of the body itself, which we know from our free precession ansatz. With respect to the magnetic-field |${\boldsymbol e}_z^{(B)}$| triad, it is given by

\begin{equation} {\boldsymbol \Omega }(t) = \omega {\boldsymbol e}_z^{(B)} + \alpha [\cos \chi {\boldsymbol e}_z^{(B)} + \sin \chi ({\boldsymbol e}_x^{(B)}\cos \omega t - {\boldsymbol e}_y^{(B)}\sin \omega t) ]. \end{equation}

(18)

From this point on, we will refer all quantities to the |${\boldsymbol e}_z^{(B)}$| triad rather than that associated with the primary rotation axis |${\boldsymbol e}_z^{(\alpha )}$|⁠; and since there is no longer any ambiguity in the axes, we will also drop the (B) superscripts on |$({\boldsymbol e}_x,{\boldsymbol e}_y,{\boldsymbol e}_z)$|⁠.

We wish to set up a perturbative scheme in which we can describe the |$\xi$|-motions. We will follow Mestel and collaborators by expanding about a non-rotating and unmagnetized spherical background, assuming that terms related to the stellar rotation and magnetic field are small. The size of these terms can be defined by the centrifugal and magnetic ellipticities ε_α and ε_B, dimensionless quantities which scale with the ratio of their respective energies to the gravitational binding energy of the spherical background star:

\begin{eqnarray} \epsilon _\alpha \sim \frac{\Omega ^2 R_*^3}{G\mathcal {M}}, \end{eqnarray}

(19)

\begin{eqnarray} \epsilon _B \sim \frac{B^2 R_*^4}{G\mathcal {M}^2}. \end{eqnarray}

(20)

Although we will assume that both ε_α and ε_B are separately small, we do not make any assumption about their relative size, so the order |$\mathcal {O}$| of each is formally different (even though they could be numerically comparable). In particular, we will treat second-order perturbations in the most general case, where

\begin{equation} \mathcal {O}(\epsilon _\alpha ^2)\ne \mathcal {O}(\epsilon _\alpha \epsilon _B) \ne \mathcal {O}(\epsilon _B^2). \end{equation}

(21)

Finally, one could envisage a different perturbative scheme with χ as the small parameter, and perturbations being performed about a background aligned-rotator model, but we prefer to be able to allow for an arbitrary degree of misalignment between the rotation and magnetic axes.

In obvious notation, we can now write any given quantity (e.g. the density) as a perturbative expansion of the form

\begin{equation} \rho = \rho _0 + \delta \rho _{\alpha } + \delta \rho _{B} + \delta \rho _{\alpha \alpha } + \delta \rho _{\alpha B} + \delta \rho _{BB} + \ldots . \end{equation}

(22)

This enables us to expand all the terms on the right-hand side of the Euler equation. Note that for the perturbative expansion for the magnetic field itself, which is of the form

\begin{equation} {\boldsymbol B} = {\boldsymbol B}_0 + \delta {\boldsymbol B}_\alpha + \delta {\boldsymbol B}_B + \ldots , \end{equation}

(23)

half-integer powers of ε_B will occur, for example,

\begin{equation} {\boldsymbol B}_0 \sim \epsilon _B^{1/2}, \hspace{28.45274pt} \delta {\boldsymbol B}_\alpha \sim \epsilon _B^{1/2} \epsilon _\alpha , \hspace{28.45274pt} \delta {\boldsymbol B}_B \sim \epsilon _B^{3/2}, \ldots . \end{equation}

(24)

For the left-hand side of the Euler equation (11), we need to think about the assumptions already made for our precessional-like motion; these fix the leading order scalings of various quantities. First, the secondary rotation ω scales as

\begin{equation} \omega \sim \epsilon _B \alpha \sim \epsilon _\alpha ^{1/2} \epsilon _B. \end{equation}

(25)

The angular velocity itself can be written as

\begin{equation} {\boldsymbol \Omega }(t) = \omega {\boldsymbol e}_z + {\boldsymbol \Omega }_\alpha , \end{equation}

(26)

where

\begin{equation} {\boldsymbol \Omega }_\alpha \equiv \alpha [\cos \chi {\boldsymbol e}_z + \sin \chi ({{\boldsymbol e}_x} \cos \omega t - {{\boldsymbol e}_y} \sin \omega t) ], \end{equation}

(27)

so that the two pieces of |$\boldsymbol \Omega$| have scalings:

\begin{eqnarray} \omega {\boldsymbol e}_z \sim \omega \sim \epsilon _\alpha ^{1/2} \epsilon _B, \end{eqnarray}

(28)

\begin{eqnarray} {\boldsymbol \Omega }_\alpha \sim \alpha \sim \epsilon _\alpha ^{1/2}. \end{eqnarray}

(29)

We assume that to leading order, the displacements |$\xi$| are sourced by the motion of the centrifugal bulge, so that

\begin{equation} \xi \sim \epsilon _\alpha R_* \sim \epsilon _\alpha . \end{equation}

(30)

To leading order, all time derivatives are due to quantities varying on the time-scale of ω⁻¹, so that, symbolically,

\begin{equation} \frac{\mathrm{\partial} }{\mathrm{\partial} t} \rightarrow \omega \sim \epsilon _\alpha ^{1/2} \epsilon _B. \end{equation}

(31)

In particular, the leading-order piece of the fluid velocity must then scale as

\begin{equation} \dot{\xi }\sim \omega \xi \sim \epsilon _\alpha ^{3/2} \epsilon _B. \end{equation}

(32)

We can then use these results to write down the scalings of the leading-order parts of the first three terms on the left-hand side of equation (11):

\begin{eqnarray} \frac{\mathrm{\partial} {\dot{\xi }}}{\mathrm{\partial} t} \sim \omega ^2 \xi \sim \epsilon _\alpha ^2 \epsilon _B^2, \end{eqnarray}

(33)

\begin{eqnarray} ({\dot{\xi }} \cdot \nabla ) {\dot{\xi }} \sim (\dot{\xi })^2 \sim \epsilon _\alpha ^3 \epsilon _B^2, \end{eqnarray}

(34)

\begin{eqnarray} {\boldsymbol \Omega } \times {\dot{\xi }} \sim \alpha (\epsilon _\alpha \epsilon _B \alpha ) \sim \epsilon _\alpha ^{2} \epsilon _B. \end{eqnarray}

(35)

We will find below that all these leading-order pieces are of sufficiently high order that they will not be needed in our analysis.

Given that we have prescribed the exact form of |${\boldsymbol \Omega }(t),$| we can compute the fourth and fifth terms on the left-hand side of equation (11) exactly, but for now, let us just note their scalings with ε_α and ε_B. For the fourth term on the left-hand side of equation (11), we have a leading-order piece that scales as

\begin{equation} \frac{\mathrm{d}{\boldsymbol {\Omega }}}{\mathrm{d}t} \times {\boldsymbol r} \sim \alpha \omega R_* \sim \epsilon _\alpha \epsilon _B, \end{equation}

(36)

which will be relevant. For the fifth term on the left-hand side of the Euler equation, there will be a leading-order piece of the form

\begin{equation} {\boldsymbol \Omega }_\alpha \times ({\boldsymbol \Omega }_\alpha \times {\boldsymbol r}) \sim \alpha ^2 \sim \epsilon _\alpha , \end{equation}

(37)

which will be needed, and the pieces given by

\begin{equation} (\omega {\boldsymbol e_{z}}) \times ({\boldsymbol \Omega _\alpha } \times {\boldsymbol r}) + {\boldsymbol \Omega _\alpha } \times [(\omega {\boldsymbol e_{z}}) \times {\boldsymbol r}] \sim \omega \alpha \sim \epsilon _\alpha \epsilon _B, \end{equation}

(38)

which will also be needed, and a piece

\begin{equation} (\omega {\boldsymbol e_{z}}) \times [(\omega {\boldsymbol e_{z}}) \times {\boldsymbol r}] \sim \omega ^2 \sim \epsilon _\alpha \epsilon _B^2, \end{equation}

(39)

which will turn out to be of too high an order to be important in this paper.

Finally, to find the equation-of-state relations needed to close each system of perturbation equations, we perform a Taylor expansion of H = H(ρ) about the point ρ = ρ₀:

\begin{eqnarray} H(\rho ) & =& H(\rho _0)+(\rho -\rho _0)\frac{\mathrm{d}{H}}{\mathrm{d}{\rho }}+\frac{(\rho -\rho _0)^2}{2!}\frac{\mathrm{d}{^2H}}{\mathrm{d}{\rho ^2}}+\ldots \nonumber\\ & =& H(\rho _0)+\delta \rho _\alpha \frac{\mathrm{d}{H}}{\mathrm{d}{\rho }}+\delta \rho _{\rm B}\frac{\mathrm{d}{H}}{\mathrm{d}{\rho }} +\delta \rho _{\alpha B}\frac{\mathrm{d}{H}}{\mathrm{d}{\rho }} + \delta \rho _\alpha \delta \rho _{\rm B}\frac{\mathrm{d}{^2H}}{\mathrm{d}{\rho ^2}}+\ldots \end{eqnarray}

(40)

Now comparing this with H = H₀ + δH_α + δH_B + δH_αB + …, we can read off the relations for different perturbative orders.

We can now insert all of these perturbative expansions into the full set of equations, to give sets of perturbation equations to be solved simultaneously. It is convenient to label these sets according to the order in ε_α and ε_B to which terms are retained in the Euler equation.

3.3 Zeroth-order equations

To zeroth order in ε_α and ε_B, the Euler equation simply gives the hydrostatic force balance of a non-rotating unmagnetized star:

\begin{equation} 0 = -\nabla H_0 - \nabla \Phi _0. \end{equation}

(41)

To the same order, we have Poisson's equation

\begin{equation} \nabla ^2 \Phi _0 = 4\pi G \rho _0, \end{equation}

(42)

and the equation of state

\begin{equation} H_0 = H_0(\rho _0). \end{equation}

(43)

The solution is static and spherical, so that only the radial component of equation (41) is non-trivial, leaving us with three equations in the three unknowns H₀, Φ₀, ρ₀, each of which will depend only upon the radial coordinate r.

3.4 Order-ε_B equations

To order ε_B, the Euler equation gives the perturbation to the spherical star caused by a magnetic field |${\boldsymbol B}_0$|⁠:

\begin{equation} 0 = -\nabla \delta H_B - \nabla \delta \Phi _B + \frac{1}{4 \pi \rho _0} (\nabla \times {\boldsymbol B}_0) \times {\boldsymbol B}_0. \end{equation}

(44)

To the same order, we have Poisson's equation

\begin{equation} \nabla ^2 \delta \Phi _B = 4\pi G \delta \rho _{\rm B}, \end{equation}

(45)

the equation of state

\begin{equation} \delta H_B = \frac{\mathrm{d}{H}}{\mathrm{d}{\rho }} \delta \rho _{\rm B}, \end{equation}

(46)

and the solenoidal constraint

\begin{equation} \nabla \cdot {\boldsymbol B}_0 = 0. \end{equation}

(47)

We have six equations in the six unknowns (B₀)_r, (B₀)_θ, (B₀)_ϕ, δH_B, δΦ_B, δρ_B. By assumption, the leading-order piece of the magnetic field, |${\boldsymbol B}_0$|⁠, is static and axisymmetric, so all these unknown quantities will be functions only of r and θ.

3.5 Order-ε_α equations

To order ε_α, the Euler equation gives the perturbation to the spherical star caused by the rotation |$\boldsymbol \Omega _\alpha$|⁠:

\begin{equation} {\boldsymbol \Omega }_\alpha \times ({\boldsymbol \Omega }_\alpha \times {\boldsymbol r}) = -\nabla \delta H_\alpha - \nabla \delta \Phi _\alpha . \end{equation}

(48)

To the same order, we have Poisson's equation

\begin{equation} \nabla ^2 \delta \Phi _\alpha = 4\pi G \delta \rho _\alpha , \end{equation}

(49)

and the equation of state

\begin{equation} \delta H_\alpha = \frac{\mathrm{d}{H}}{\mathrm{d}{\rho }} \delta \rho _\alpha . \end{equation}

(50)

As can be expected given the form of equation (27) above, the centrifugal force term on the left-hand side of the above perturbed Euler equation is time varying, and has no particular symmetry, so we can expect the unknown quantities (δH_α, δρ_α, δΦ_α) to depend upon all four coordinates (t, r, θ, ϕ). The centrifugal term on the left-hand side of equation (48) can easily be written as the gradient of a scalar, so that the equation reduces to a single scalar equation, leaving three equations in the three unknowns.

3.6 Order-ε_αε_B equations

To order ε_αε_B, the Euler equation gives the perturbation to the spherical star caused by the interaction between the rotation and the magnetic field. Explicitly, we find

\begin{eqnarray} \frac{\mathrm{d}{\boldsymbol \Omega _\alpha }}{\mathrm{d}t} \times {\boldsymbol r} + \omega {\boldsymbol e}_z \times ({\boldsymbol \Omega _\alpha } \times {\boldsymbol r}) + {\boldsymbol \Omega _\alpha } \times (\omega {\boldsymbol e}_z \times {\boldsymbol r}) &=& -\nabla \delta H_{\alpha B} - \nabla \delta \Phi _{\alpha B} - \frac{\delta \rho _\alpha }{4\pi \rho _0^2} (\nabla \times {\boldsymbol B}_0) \times {\boldsymbol B}_0\nonumber\\ && + \, \frac{1}{4\pi \rho _0} [ (\nabla \times \delta {\boldsymbol B}_\alpha ) \times {\boldsymbol B}_0 + (\nabla \times {\boldsymbol B}_0) \times \delta {\boldsymbol B}_\alpha ]. \end{eqnarray}

(51)

To the same order, we have Poisson's equation

\begin{equation} \nabla ^2 \delta \Phi _{\alpha B} = 4\pi G \delta \rho _{\alpha B}, \end{equation}

(52)

and the equation of state

\begin{equation} \delta H_{\alpha B} = \frac{\mathrm{d}{H}}{\mathrm{d}{\rho }} \delta \rho _{\alpha B} + \frac{\mathrm{d}{^2 H}}{\mathrm{d}{\rho ^2}} \delta \rho _{\alpha }\delta \rho _{\rm B}. \end{equation}

(53)

To close these equations, we can make use of the continuity and the induction equations. The leading-order non-zero part of the continuity equation contains terms that scale as |$\epsilon _\alpha ^{3/2} \epsilon _B$|⁠:

\begin{equation} \frac{\mathrm{\partial} \delta \rho _\alpha }{\mathrm{\partial} t} + \nabla \cdot (\rho _0 \dot{\xi }) = 0. \end{equation}

(54)

The leading-order non-zero part of the induction equation contains terms that scale as |$\epsilon _\alpha ^{3/2} \epsilon _B^{3/2}$|⁠:

\begin{equation} \frac{\mathrm{\partial} \delta {\boldsymbol B}_\alpha }{\mathrm{\partial} t} - \nabla \times (\dot{\xi }\times {\boldsymbol B}_0) = 0. \end{equation}

(55)

Together, these form a set of nine equations in the nine unknowns |$(\delta \rho _{\alpha B}, \delta H_{\alpha B}, \delta \Phi _{\alpha B},\dot{\xi }, \delta {\boldsymbol B}_\alpha$|⁠). One could add the further equation

\begin{equation} \nabla \cdot \delta {\boldsymbol B}_\alpha = 0, \end{equation}

(56)

but this is redundant, as the induction equation guarantees that this constraint is preserved by the evolution, so provided the initial data are divergence-free, the solution at later times will be too. Since this is the only perturbation of |${\boldsymbol B}$| we will need to consider, we will drop its α subscript in future.

3.7 Equations of order |$\epsilon _\alpha ^2$|⁠, |$\epsilon _B^2$| and higher

As discussed earlier, there are three different sets of second-order perturbation equations: one set with |$\mathcal {O}(\epsilon _\alpha ^2)$| terms, another with |$\mathcal {O}(\epsilon _B^2)$| terms and a third with |$\mathcal {O}(\epsilon _\alpha \epsilon _B)$| terms. We saw from the previous section that the latter set encodes the non-rigid motions in which we are interested. The reason for neglecting the other two is not that they are of a higher order – in fact, |$\mathcal {O}(\epsilon _\alpha \epsilon _B)$| quantities will never numerically be the largest of the three – but simply that they contain more mundane information.

The order-|$\epsilon _\alpha ^2$| equations merely describe a correction δρ_αα to the result of calculating the centrifugal bulge from the order-ε_α equations. This correction scales with α⁴ and therefore does not become important except for very rapidly rotating stellar models. Analogously, the order-|$\epsilon _B^2$| equations encode a correction δρ_BB, scaling with B⁴, to the magnetic distortion as given by the order-ε_B calculation. Here, a stronger statement may be made: there is no known star for which the magnetic energy is anywhere near large enough to warrant including this higher order correction (Reisenegger 2009).

Although we will not consider any higher order quantities than those discussed, we note that the density perturbation δρ_αB should be time-dependent, and therefore through a higher order analogue of the continuity equation (54) will induce a velocity perturbation of a higher order than |$\dot{\xi }$|⁠. Similarly, δρ_αα induces a velocity field, but merely a rapid-rotation correction to |$\dot{\xi }$|⁠. Finally, δρ_BB is stationary, so the associated velocity correction is identically zero.

4 SOLUTION OF THE ZEROTH- AND FIRST-ORDER EQUATIONS

Before tackling the second-order |$\mathcal {O}(\epsilon _\alpha \epsilon _B)$| system of equations we wish to solve, we first need results from the three lower order systems of equations: the zeroth-order, |$\mathcal {O}(\epsilon _B)$| and |$\mathcal {O}(\epsilon _\alpha )$| sets. The first two are readily found in the literature, but the latter calculation is non-standard and so we report it in full.

We will specialize to the case of a γ = 2 polytropic equation of state. This has the advantage that virtually the entire calculation may be carried out analytically, and is the simplest case for which the fluid motion is compressible. Although this polytrope better mimics the pressure–density distribution of neutron stars than for other classes of star, we will argue later that many features of our solutions are likely to be applicable to a generic magnetic star. A detailed discussion of the applicability, and limitations, of our model may be found in Section 8.

4.1 The zeroth-order equations

The density profile of a γ = 2 polytrope in hydrostatic equilibrium is a classic result (Chandrasekhar 1939), and is given by

\begin{equation} \rho _0 = \rho _0(r) = \rho _{\rm c}\frac{R_*\sin \left({{ {\frac{\pi r}{R_*}}}}\right)}{\pi r}, \end{equation}

(57)

where ρ_c is the central density.

4.2 The order-ε_B equations

The study of magnetic-field distributions in axisymmetric stars also has a considerable pedigree. For purely poloidal and poloidal-toroidal fields, one needs to solve for both the magnitude and direction of the field, encapsulated in a magnetic streamfunction (Chandrasekhar & Prendergast 1956; Monaghan 1976). Purely toroidal fields are the simplest of all, however, as the direction is known (azimuthal, i.e. along |${\boldsymbol e}_\phi$|⁠) and one need only determine the functional form of the magnitude. One such solution for the toroidal field is

\begin{equation} {\boldsymbol B}_0 = B_\phi {\boldsymbol e}_\phi = \Lambda \rho _0 r\sin \theta {\boldsymbol e}_\phi , \end{equation}

(58)

where Λ is a constant governing the strength of the field (Roxburgh 1966). This corresponds to a Lorentz force per unit mass |$(\nabla \times {\boldsymbol B}_0)\times {\boldsymbol B}_0/4\pi \rho _0$|⁠, which is the gradient of the following scalar (Lander & Jones 2009):

\begin{equation} M = -\frac{\Lambda ^2}{4\pi } r^2\sin ^2\theta \rho _0 = -\frac{\Lambda ^2 R_*\rho _c}{4\pi ^2} r\sin \left({ {\frac{\pi r}{R_*}}}\right)\sin ^2\theta . \end{equation}

(59)

Later we will also need to relate the ellipticity ε_B to the coefficient Λ. It seems, however, that a closed-form expression for ε_B is only possible for an incompressible star. For our γ = 2 polytrope, we will need to be content with using the scaling of ε_B from equation (20), together with the fact that the mass and average magnetic field 〈B〉 scale as follows:

\begin{equation} \mathcal {M}\sim \rho _{\rm c} R_*^3,\ \ \langle B\rangle \sim \Lambda \rho _{\rm c} R_*. \end{equation}

(60)

The scaling of the ellipticity with stellar parameters is then

\begin{equation} \epsilon _B \sim \frac{\langle B\rangle {}^2 R_*^4}{G\mathcal {M}^2} \sim \frac{\Lambda ^2}{G}. \end{equation}

(61)

ε_B and Λ²/G are related by a dimensionless constant of proportionality k_B, which must be determined numerically. Using the code from Lander & Jones (2009), we find that for a γ = 2 polytrope,

\begin{equation} \epsilon _B = k_B \frac{\Lambda ^2}{G} \approx -0.019 \frac{\Lambda ^2}{G}, \end{equation}

(62)

although the constant depends only weakly on the equation of state; it is again −0.019 for a polytrope with γ = 5/3, and is −0.015 for γ = 4/3.

This completes the description of order-ε_B quantities needed as input for the second-order perturbation equations later. In particular, although the quantities δρ_B, δH_B or δΦ_B appear directly in the order-ε_B equations, we will not need explicit forms for these.

4.3 The order-ε_α equations

In the ω-frame, the centrifugal distortion δρ_α is neither stationary nor axisymmetric. It is sourced by the centrifugal force |$\boldsymbol \Omega _\alpha \times (\boldsymbol \Omega _\alpha \times {\boldsymbol r})$|⁠, whose form we know from our free-precession ansatz on |$\boldsymbol \Omega _\alpha$|⁠. Our first task is to write this force as the gradient of a scalar function.

4.3.1 The centrifugal force

We will work in spherical polar coordinates, and begin by converting the Cartesian expression for |$\boldsymbol \Omega _\alpha$| from equation (27) into spherical-polar form:

\begin{equation} \boldsymbol \Omega _\alpha = \alpha \lbrace [\cos \chi \cos \theta +\sin \chi \sin \theta \cos (\omega t+\phi )]{\boldsymbol e}_r + [-\cos \chi \sin \theta +\sin \chi \cos \theta \cos (\omega t+\phi )]{\boldsymbol e}_\theta - \sin \chi \sin (\omega t+\phi ){\boldsymbol e}_\phi \rbrace . \end{equation}

(63)

We use this to calculate the r, θ and ϕ components of the centrifugal force. To express this force as the gradient of some scalar F_α, we compare the three components of the force with the general expression for ∇F_α. After suitable integrations, we can find F_α explicitly and thus write the centrifugal force as

\begin{equation} \boldsymbol \Omega _\alpha \times (\boldsymbol \Omega _\alpha \times {\boldsymbol r}) = \nabla \left\lbrace \frac{\alpha ^2 r^2}{4} \left[ -2\sin ^2\theta + \sin ^2\chi (3\sin ^2\theta -2) + \sin (2\chi )\sin (2\theta )\cos \lambda + \sin ^2\chi \sin ^2\theta \cos (2\lambda ) \right] \right \rbrace , \end{equation}

(64)

which reduces to the standard result for a stationary rotating star in the limit χ → 0, as expected. In equation (64), we introduced a time-shifted azimuthal coordinate

\begin{equation} \lambda \equiv \phi +\omega t \end{equation}

(65)

to simplify the expression. Note that we may freely replace the original azimuthal coordinate ϕ with the new one λ, since it only acts to redefine the origin of time in our equations, and since derivatives are unaffected:

\begin{equation} \frac{\mathrm{\partial} {f(\theta ,\lambda )}}{\mathrm{\partial} {\phi }} = \frac{\mathrm{\partial} {f}}{\mathrm{\partial} {\lambda }}\frac{\mathrm{\partial} {\lambda }}{\mathrm{\partial} {\phi }} = \frac{\mathrm{\partial} {f}}{\mathrm{\partial} {\lambda }}\frac{\mathrm{\partial} { }}{\mathrm{\partial} {\phi }}(\omega t+\phi ) = \frac{\mathrm{\partial} {f}}{\mathrm{\partial} {\lambda }} \Rightarrow \frac{\mathrm{\partial} { }}{\mathrm{\partial} {\phi }} = \frac{\mathrm{\partial} { }}{\mathrm{\partial} {\lambda }}. \end{equation}

(66)

The stationary background field |${\boldsymbol B}_0$| discussed in the previous section may, of course, be written equivalently in terms of λ, i.e. |${\boldsymbol B}_0 = B_\phi {\boldsymbol e}_\phi = B_\lambda {\boldsymbol e}_\lambda$|⁠.

4.3.2 The |$\mathcal {O}(\epsilon _\alpha )$| system as a single Helmholtz equation

Having found an expression for the perturbed centrifugal force, we now turn to the set of |$\mathcal {O}(\epsilon _\alpha )$| equations. We wish to reduce the three equations in three variables to a single equation in one variable. Although it is equivalent, in principle, to work with a final equation in either δρ_α, δH_α or δΦ_α, we choose the latter – as only in this case is it straightforward to impose the necessary boundary conditions at the stellar surface. For a polytropic equation of state, the enthalpy may be written as

\begin{equation} H = \frac{\gamma P}{(\gamma -1)\rho } = \frac{\gamma k\rho ^{\gamma -1}}{\gamma -1}. \end{equation}

(67)

Substituting this into the |$\mathcal {O}(\epsilon _\alpha )$| Euler equation, we have

\begin{equation} \boldsymbol \Omega _\alpha \times (\boldsymbol \Omega _\alpha \times {\boldsymbol r}) = -k\gamma \delta \rho _\alpha \nabla (\rho ^{\gamma -2}) -k\gamma \rho ^{\gamma -2}\nabla \delta \rho _\alpha - \nabla \delta \Phi _\alpha . \end{equation}

(68)

This is clearly particularly simple for our chosen polytropic index, γ = 2. Together with the result of (64), the equation becomes

\begin{equation} \nabla \Big \lbrace \frac{\alpha ^2 r^2}{4} \Big [ -2\sin ^2\theta + \sin ^2\chi (3\sin ^2\theta -2) + \sin (2\chi )\sin (2\theta )\cos \lambda + \sin ^2\chi \sin ^2\theta \cos (2\lambda ) \Big ] \Big \rbrace = -2k\nabla \delta \rho _\alpha -\nabla \delta \Phi _\alpha , \end{equation}

(69)

whose first integral is

\begin{equation} C+\frac{\alpha ^2 r^2}{4} \Big [ -2\sin ^2\theta + \sin ^2\chi (3\sin ^2\theta -2) + \sin (2\chi )\sin (2\theta )\cos \lambda + \sin ^2\chi \sin ^2\theta \cos (2\lambda ) \Big ] = -2k\delta \rho _\alpha -\delta \Phi _\alpha , \end{equation}

(70)

where C is an integration constant. Finally, we use the Poisson equation in the above to replace δρ_α, to give

\begin{equation} {\left({\nabla ^2+\frac{2\pi G}{k}}\right)\delta \Phi _\alpha = -\frac{2\pi G}{k}\Big \lbrace C + \frac{\alpha ^2 r^2}{4} \Big [ -2\sin ^2\theta + \sin ^2\chi (3\sin ^2\theta -2) + \sin (2\chi )\sin (2\theta )\cos \lambda + \sin ^2\chi \sin ^2\theta \cos (2\lambda ) \Big ] \Big \rbrace .} \end{equation}

(71)

This is an inhomogeneous Helmholtz equation, which we solve next, beginning with the homogeneous solution.

4.3.3 Homogeneous solution

The general solution of the homogeneous Helmholtz equation

\begin{equation} {\left({\nabla ^2+\frac{2\pi G}{k}}\right)\delta \Phi _\alpha ^{\rm homog} = 0} \end{equation}

(72)

is a standard result (Arfken & Weber 2005), which, in spherical polar coordinates, may be written as the infinite sum

\begin{equation} \delta \Phi _\alpha ^{\rm homog} = \sum \limits _{l = 0}^{\infty } \sum \limits _{m = -l}^{l} b_l^m j_l\left({ {\sqrt{\frac{2\pi G}{k}}r}}\right) Y_l^m(\theta ,\lambda ), \end{equation}

(73)

where |$b_l^m$| are constants, j_l are spherical Bessel functions and |$Y_l^m$| are spherical harmonics.

4.3.4 Particular solution

We begin by decomposing the various angular pieces of the centrifugal potential (64) into spherical harmonics:

\begin{eqnarray} 1 = 2\sqrt{\pi } Y_0^0, \end{eqnarray}

(74)

\begin{eqnarray} &&{\sin ^2\theta = \frac{2}{3}\left({1 - P_2^0}\right) = -\frac{4\sqrt{\pi }}{3}\left({Y_0^0-\frac{1}{\sqrt{5}}Y_2^0}\right),} \end{eqnarray}

(75)

\begin{eqnarray} &&{\sin (2\theta )\cos \lambda = \sin (2\theta )\left({\frac{\mathrm{e}^{i\lambda }}{2}+\frac{\mathrm{e}^{-i\lambda }}{2}}\right) = -\frac{1}{3}P_2^1 \mathrm{e}^{i\lambda } + 2P_2^{-1}\mathrm{e}^{-i\lambda } = -2\sqrt{\frac{2\pi }{15}}\left({Y_2^1 - Y_2^{-1}}\right),} \end{eqnarray}

(76)

\begin{eqnarray} &&{\sin ^2\theta \cos (2\lambda ) = \sin ^2\theta \left({\frac{\mathrm{e}^{2i\lambda }}{2}+\frac{\mathrm{e}^{-2i\lambda }}{2}}\right) = \frac{1}{6}P_2^2 \mathrm{e}^{2i\lambda } + 4P_2^{-2}\mathrm{e}^{-2i\lambda } = 4\sqrt{\frac{\pi }{30}}\left({Y_2^2 + Y_2^{-2}}\right).} \end{eqnarray}

(77)

Using the above relations and equation (64), the right-hand side of the inhomogeneous Helmholtz equation (71) may, therefore, be written as a sum over l = 0, 2 and m = 0, ±2:

\begin{equation} {-\frac{2\pi ^{3/2} G\alpha ^2 r^2}{k} \Bigg [ 2\left({\frac{C}{\alpha ^2 r^2}-\frac{1}{3}}\right)Y_0^0 +\frac{1}{\sqrt{5}}\left({\frac{2}{3}-\sin ^2\chi }\right)Y_2^0 -\frac{\sin (2\chi )}{\sqrt{30}}(Y_2^1-Y_2^{-1}) +\frac{\sin ^2\chi }{\sqrt{30}}(Y_2^2-Y_2^{-2}) \Bigg ]} \end{equation}

(78)

and so we expect the same angular structure for the particular solution |$\delta \Phi _\alpha ^{\rm PS}$|⁠:

\begin{equation} \delta \Phi _\alpha ^{\rm PS} = \delta \Phi _0^0(r)Y_0^0 + \delta \Phi _2^0(r)Y_2^0 + \delta \Phi _2^1(r)Y_2^1 + \delta \Phi _2^{-1}(r)Y_2^{-1} + \delta \Phi _2^2(r)Y_2^2 + \delta \Phi _2^{-2}(r)Y_2^{-2}. \end{equation}

(79)

Now, using

\begin{equation} \nabla ^2(f(r)Y_l^m) = f(r)\nabla ^2 Y_l^m + Y_l^m\left({\frac{\mathrm{d}{^2}}{\mathrm{d}{r^2}}+\frac{2}{r}\frac{\mathrm{d}{ }}{\mathrm{d}{r}}}\right)f(r) \end{equation}

(80)

and

\begin{equation} \nabla ^2 Y_l^m = -\frac{l(l+1)}{r^2}Y_l^m, \end{equation}

(81)

we rewrite the left-hand side of (71):

\begin{equation} {\left({\nabla ^2+\frac{2\pi G}{k}}\right)\delta \Phi _\alpha ^{\rm PS} = \sum \limits _{l = 0}^2\sum _{m = -l}^l \left({\frac{\mathrm{d}{^2}}{\mathrm{d}{r^2}}+\frac{2}{r}\frac{\mathrm{d}{ }}{\mathrm{d}{r}}+\frac{2\pi G}{k}-\frac{l(l+1)}{r^2}}\right)\delta \Phi _l^m(r)Y_l^m.} \end{equation}

(82)

Next, we solve for the radial functions |$\delta \Phi _l^m$| by equating the different |$Y_l^m$| terms from the left-hand and right-hand sides of the inhomogeneous Helmholtz equation, i.e. the terms from (78) with their counterparts from equation (82). Before continuing, we recall that our final aim is to calculate time-dependent quantities in a precessing, magnetized fluid star: the perturbations to the velocity and magnetic field, which oscillate over a precession time-scale. However, two of the six terms in δΦ_α are stationary:

\begin{equation} \frac{\mathrm{\partial} {\delta \Phi _0^0}}{\mathrm{\partial} {t}} = \frac{\mathrm{\partial} {\delta \Phi _2^0}}{\mathrm{\partial} {t}} = 0, \end{equation}

(83)

so we need not solve for these terms, and will drop them in the rest of the calculation. Next, we turn to the time-dependent terms. Equating components of |$Y_2^1$|⁠, we have

\begin{equation} {\left({\frac{\mathrm{d}{^2}}{\mathrm{d}{r^2}}+\frac{2}{r}\frac{\mathrm{d}{ }}{\mathrm{d}{r}}+\frac{2\pi G}{k}-\frac{6}{r^2}}\right)\delta \Phi _2^1 = \frac{2\sin (2\chi )\pi ^{3/2} G\alpha ^2 r^2}{\sqrt{30} k}.} \end{equation}

(84)

To solve this, we make the ansatz that |$\delta \Phi _2^1$| is a quadratic function of r, finding that its constant and linear pieces are zero, so that

\begin{equation} \delta \Phi _2^1 = \sqrt{\frac{\pi }{30}} \sin (2\chi )\alpha ^2 r^2. \end{equation}

(85)

A similar procedure applied to the other components gives

\begin{equation} \delta \Phi _2^{-1} = -\delta \Phi _2^1 = -\sqrt{\frac{\pi }{30}} \sin (2\chi )\alpha ^2 r^2,\ \ \delta \Phi _2^2 = \delta \Phi _2^{-2} = -\sqrt{\frac{\pi }{30}} \sin ^2\chi \alpha ^2 r^2. \end{equation}

(86)

4.3.5 Full interior solution

The full solution for the perturbed gravitational potential inside the star |$\delta \Phi _\alpha ^{\rm int}$| – with yet-to-be-determined coefficients |$b_l^m$| – is given by the sum of the particular solution and the homogeneous solution. For the latter, we need only consider the terms in the infinite sum with the same multipolar order as the time-dependent pieces in the particular solution (again, since δΦ_α is time-dependent by ansatz):

\begin{eqnarray} \delta \Phi _\alpha ^{\rm int} &=& \left({b_2^1j_2\left({{\sqrt{\frac{2\pi G}{k}}}r}\right) + \sqrt{\frac{\pi }{30}}\sin (2\chi )\alpha ^2 r^2 }\right)Y_2^1 + \left({b_2^{-1}j_2\left({{\sqrt{\frac{2\pi G}{k}}}r}\right) - \sqrt{\frac{\pi }{30}}\sin (2\chi )\alpha ^2 r^2 }\right)Y_2^{-1}\nonumber \\ &&+\, \left({b_2^2j_2\left({{\sqrt{\frac{2\pi G}{k}}}r}\right) - \sqrt{\frac{\pi }{30}}\sin ^2\chi \alpha ^2 r^2 }\right)Y_2^2 + \left({b_2^{-2}j_2\left({{\sqrt{\frac{2\pi G}{k}}}r}\right) - \sqrt{\frac{\pi }{30}}\sin ^2\chi \alpha ^2 r^2 }\right)Y_2^{-2}. \end{eqnarray}

(87)

4.3.6 Exterior solution and boundary conditions

In order to fix the constants |$b_l^m$| for the interior gravitational potential, we need to match the interior and exterior gravitational potentials appropriately at the stellar surface. The exterior gravitational potential |$\delta \Phi _\alpha ^{\rm ext}$| is governed by Laplace's equation:

\begin{equation} \nabla ^2\delta \Phi _\alpha ^{\rm ext} = 0. \end{equation}

(88)

Discarding unphysical terms that diverge at infinity, and keeping only those multipoles that feature in the interior solution, we are left with

\begin{equation} \delta \Phi _\alpha ^{\rm ext} = \frac{1}{r^3}\left({ c_2^1 Y_2^1 + c_2^{-1} Y_2^{-1} + c_2^2 Y_2^2 + c_2^{-2} Y_2^{-2} }\right). \end{equation}

(89)

We need to match the four multipoles above to the corresponding ones in the interior solution at the stellar surface R_*, such that

\begin{eqnarray} \delta \Phi _\alpha ^{\rm int}|_{r = R*} =\delta \Phi _\alpha ^{\rm ext}|_{r = R*}, \end{eqnarray}

(90)

\begin{eqnarray} \nabla \delta \Phi _\alpha ^{\rm int}|_{r = R*} =\nabla \delta \Phi _\alpha ^{\rm ext}|_{r = R*}. \end{eqnarray}

(91)

Since the angular pieces of both interior and exterior solutions are the same |$Y_l^m$| functions, imposing the boundary conditions reduces to matching the radial functions of interior and exterior solutions, and the radial derivatives of these. To perform this matching, we now use the fact that ε_α ≪ 1 and |ε_B| ≪ 1, so that it is legitimate to assume the stellar surface for the time-independent background is spherical. For a γ = 2 polytrope in hydrostatic equilibrium, the stellar surface is located at a radius (Chandrasekhar 1939) of

\begin{equation} R_* = \sqrt{\frac{\pi k}{2G}}. \end{equation}

(92)

Matching the interior and exterior solutions at this radius, we find that the four constants in the interior gravitational potential solution are

\begin{eqnarray} b_2^1 = -b_2^{-1} = -\sqrt{\frac{5}{6}}\frac{\pi ^{3/2}k\alpha ^2}{2G}\sin (2\chi ), \end{eqnarray}

(93)

\begin{eqnarray} b_2^2 = b_2^{-2} = \sqrt{\frac{5}{6}}\frac{\pi ^{3/2}k\alpha ^2}{2G}\sin ^2\chi . \end{eqnarray}

(94)

4.3.7 The first-order perturbation in the density

Finally, we may use Poisson's equation and our expression for |$\delta \Phi _\alpha ^{\rm int}$| to find the time-dependent piece of the centrifugal bulge, and therefore of the mass distribution of our precessing star:

\begin{equation} \delta \rho _\alpha = -\frac{5\pi \alpha ^2}{16 G} j_2\left({{ {\frac{\pi r}{R_*}}}}\right) \left({ \sin (2\chi )\sin (2\theta )\cos \lambda + \sin ^2\chi \sin ^2\theta \cos (2\lambda ) }\right)\!, \end{equation}

(95)

where we have used the relation |$R_* = \sqrt{\pi k/2G}$| to simplify the result and replace all instances of the polytropic constant k. The angular dependence of this result is in agreement with that of Mestel & Takhar (1972), but unlike their expression, we have calculated the explicit radial dependence (for the specific case of a γ = 2 polytrope).

5 SOLUTION OF THE SECOND-ORDER EQUATIONS

We are interested in the time-dependent velocity and magnetic fields that are generated in a precessing magnetized fluid star. Looking at the system of second-order equations (51)–(55), however, we see that higher order fluid terms (e.g. δH_αB) are also present. These quantities represent the tiny distortion to the star induced by the perturbed Lorentz force, which we are not interested in solving for. We can remove them and simplify the system of equations by taking the curl of the Euler equation (51); the equation of state then becomes redundant, as does the Poisson equation. The centrifugal-force terms from the Euler equation also vanish, since these may be written as the gradient of a scalar:

\begin{equation} (\omega {\boldsymbol e}_z)\times (\boldsymbol \Omega _\alpha \times {\boldsymbol r}) + \boldsymbol \Omega _\alpha \times ((\omega {\boldsymbol e}_z)\times {\boldsymbol r}) = \nabla \Big [\frac{\alpha \omega r^2}{2} (\sin \chi \sin (2\theta )\cos \lambda -2\cos \chi \sin ^2\theta )\Big ]. \end{equation}

(96)

The remaining set of perturbation equations reads:

\begin{eqnarray} \nabla \times \Big \lbrace \frac{\mathrm{d}{\boldsymbol \Omega _\alpha }}{\mathrm{d}{t}} \times {\boldsymbol r} \Big \rbrace = \nabla \times \Big \lbrace - \frac{\delta \rho _\alpha }{4\pi {\rho _0}^2} (\nabla \times {\boldsymbol B}_0) \times {\boldsymbol B}_0 + \frac{1}{4\pi \rho _0} [ (\nabla \times \delta {\boldsymbol B}) \times {\boldsymbol B}_0 + (\nabla \times {\boldsymbol B}_0) \times \delta {\boldsymbol B}] \Big \rbrace , \end{eqnarray}

(97)

\begin{eqnarray} \frac{\mathrm{\partial} {\delta \rho _\alpha }}{\mathrm{\partial} {t}} = -\nabla \cdot (\rho _0\dot{\xi }), \end{eqnarray}

(98)

\begin{eqnarray} \frac{\mathrm{\partial} {\delta {\boldsymbol B}}}{\mathrm{\partial} {t}} = \nabla \times ({\dot{\xi }} \times {\boldsymbol B}_0), \end{eqnarray}

(99)

\begin{eqnarray} \nabla \cdot \delta {\boldsymbol B}= 0. \end{eqnarray}

(100)

Now, the curled Euler equation has only the three components of |$\delta {\boldsymbol B}$| as unknown variables. Because each term of the equation is divergence-free, however (since the divergence of a curl is zero), it may be expressed as the sum of a poloidal and a toroidal piece. Thus, this curled Euler equation only involves two degrees of freedom, not three; it must be reducible to one scalar equality governing the various poloidal components, and another for the toroidal components. Together with the solenoidal nature of |$\delta {\boldsymbol B}$|⁠, we get a well-defined system of three equations in three unknowns (the components of |$\delta {\boldsymbol B}$|⁠):

\begin{eqnarray} \Big [\nabla \times \Big \lbrace \frac{\mathrm{d}{\boldsymbol \Omega _\alpha }}{\mathrm{d}{t}} \times {\boldsymbol r} +\frac{\delta \rho _\alpha }{4\pi \rho _0^2} (\nabla \times {\boldsymbol B}_0) \times {\boldsymbol B}_0\Big \rbrace \Big ]_{\rm {poloidal}} = \Big [\nabla \times \Big \lbrace \frac{1}{4\pi \rho _0} [ (\nabla \times \delta {\boldsymbol B}) \times {\boldsymbol B}_0 + (\nabla \times {\boldsymbol B}_0) \times \delta {\boldsymbol B}] \Big \rbrace \Big ]_{\rm {poloidal}}, \end{eqnarray}

(101)

\begin{eqnarray} \Big [\nabla \times \Big \lbrace \frac{\mathrm{d}{\boldsymbol \Omega _\alpha }}{\mathrm{d}{t}} \times {\boldsymbol r} +\frac{\delta \rho _\alpha }{4\pi \rho _0^2} (\nabla \times {\boldsymbol B}_0) \times {\boldsymbol B}_0\Big \rbrace \Big ]_{\rm {toroidal}} = \Big [\nabla \times \Big \lbrace \frac{1}{4\pi \rho _0} [ (\nabla \times \delta {\boldsymbol B}) \times {\boldsymbol B}_0 + (\nabla \times {\boldsymbol B}_0) \times \delta {\boldsymbol B}] \Big \rbrace \Big ]_{\rm {toroidal}}, \end{eqnarray}

(102)

\begin{eqnarray} \nabla \cdot \delta {\boldsymbol B}= 0. \end{eqnarray}

(103)

In the full set of Euler equations, the divergence-free nature of |$\delta {\boldsymbol B}$| did not count as an equation, merely a constraint, giving an initial condition |$\nabla \cdot \delta {\boldsymbol B}(t = 0) = 0$|⁠. The induction equation would then ensure |$\nabla \cdot \delta {\boldsymbol B}= 0$| for all times. In the above system, we have no induction equation, so |$\nabla \cdot \delta {\boldsymbol B}= 0$| is elevated to the role of an equation in its own right. This is also the case in, for example, hydromagnetic equilibrium equations; cf. Section 3.4.

Now, assume that we have solved the above for |$\delta {\boldsymbol B}$|⁠. We have two equations left over – the continuity and induction equations:

\begin{eqnarray} \frac{\mathrm{\partial} {\delta \rho _\alpha }}{\mathrm{\partial} {t}} = -\nabla \cdot (\rho _0\dot{\xi }), \end{eqnarray}

(104)

\begin{eqnarray} \frac{\mathrm{\partial} {\delta {\boldsymbol B}}}{\mathrm{\partial} {t}} = \nabla \times ({\dot{\xi }} \times {\boldsymbol B}_0). \end{eqnarray}

(105)

These represent three equations (not four, since the induction equation again only represents two equations – a poloidal and a toroidal degree of freedom) in three unknowns: the components of the velocity |$\dot{\xi }$|⁠. At this stage, having already found the perturbed magnetic field |$\delta {\boldsymbol B}$|⁠, one then has enough information to use the continuity and induction equations to obtain a solution for |$\dot{\xi }$|⁠; see Section 6.

The above gives us a strategy to solve the second-order perturbation equations, but to proceed, we need to find explicit forms of the schematically written equations (101) and (102). To this end, we first describe a general decomposition of a solenoidal field.

5.1 Preliminaries on solenoidal fields

If a vector field |${\boldsymbol F}$| is divergence-free, it has one fewer degree of freedom than its number of dimensions – in the case of a three-dimensional star, it has two degrees of freedom. These are the poloidal and toroidal components of the vector field. These two components may be written in terms of the gradients of scalar functions⁵ Ψ and ϒ, in what is sometimes known as the Mie representation of the vector field (see e.g. Backus 1986):

\begin{equation} {\bf F} = \underbrace{\nabla \times (\Psi {\boldsymbol r})}_{{\rm toroidal}} + \underbrace{\nabla \times (\nabla \Upsilon \times {\boldsymbol r})}_{{\rm poloidal}}. \end{equation}

(106)

Note that unlike the special case of axisymmetric fields, we cannot generally orientate our axes so that the toroidal unit vector is aligned with the azimuthal direction. Now, if we have an equation involving a number of different poloidal and toroidal vector fields, we can separate the equation into one equality governing the toroidal fields and another for the poloidal fields. Furthermore, we can reduce each of these equations to a relation involving just the poloidal/toroidal scalar functions, e.g. if |${\boldsymbol P}_1, {\boldsymbol P}_2, {\boldsymbol P}_3$| are three poloidal fields, then

\begin{eqnarray} {\boldsymbol P}_1 + {\boldsymbol P}_2 = {\boldsymbol P}_3 &\Rightarrow & \nabla \times (\nabla \Upsilon _1\times {\boldsymbol r}) + \nabla \times (\nabla \Upsilon _2\times {\boldsymbol r}) = \nabla \times (\nabla \Upsilon _3\times {\boldsymbol r})\nonumber \\ &\Rightarrow &\nabla \times (\nabla [\Upsilon _1+\Upsilon _2-\Upsilon _3]\times {\boldsymbol r}) = 0\nonumber \\ &\Rightarrow &\Upsilon _1+\Upsilon _2-\Upsilon _3 = \mathcal {A}, \end{eqnarray}

(107)

using the distributive properties of the cross product and curl. |$\mathcal {A}$| represents the gauge freedom that one may add on any terms to the solution, which satisfy

\begin{equation} \nabla \times (\nabla \mathcal {A}\times {\boldsymbol r}) = 0. \end{equation}

(108)

One may equate the scalar functions of toroidal fields in an analogous manner; in this case one may add to the solution any term |$\mathcal {B}$| that satisfies

\begin{equation} \nabla \mathcal {B}\times {\boldsymbol r}= 0. \end{equation}

(109)

5.2 An explicit toroidal–poloidal split of the curled Euler equation

Our first task is to find explicit expressions for the poloidal and toroidal terms in the curled Euler equations (101) and (102), by comparing the terms with the Mie representation (106).

First, note that the time derivative of the perturbed angular velocity may be written as the gradient of a scalar:

\begin{equation} \frac{\mathrm{d}{\boldsymbol \Omega _\alpha }}{\mathrm{d}{t}} = \nabla \left({-\alpha \omega \sin \chi r\sin \theta \sin \lambda }\right) = \nabla \left({-\alpha ^2\sin \chi \cos \chi \epsilon _B r\sin \theta \sin \lambda }\right)\!. \end{equation}

(110)

The curl of the Euler force |$\nabla \times ({\rm d}\boldsymbol \Omega _\alpha /{\rm d}t \times {\boldsymbol r})$| is therefore a purely poloidal vector.

The next term from (101) and (102) involves δρ_α and the background Lorentz force. We may simplify this using the result that the Lorentz force divided by the density is the gradient of a scalar M:

\begin{equation} \nabla \times \Big \lbrace \frac{\delta \rho _\alpha }{4\pi {\rho _0^2}} (\nabla \times {\boldsymbol B}_0) \times {\boldsymbol B}_0 \Big \rbrace = \nabla \times \left({\frac{\delta \rho _\alpha }{\rho _0}\nabla M}\right) = \nabla \left({\frac{\delta \rho _\alpha }{\rho _0}}\right) \times \nabla M. \end{equation}

(111)

Next, expand the two gradients in components:

\begin{eqnarray} \nabla \left({\frac{\delta \rho _\alpha }{\rho _0}}\right) \times \nabla M &=& \left({ \frac{\mathrm{\partial} {(\delta \rho _\alpha /\rho _0)}}{\mathrm{\partial} {r}}{\boldsymbol e}_r + \frac{1}{r}\frac{\mathrm{\partial} {(\delta \rho _\alpha /\rho _0)}}{\mathrm{\partial} {\theta }}{\boldsymbol e}_\theta + \frac{1}{r\sin \theta }\frac{\mathrm{\partial} {(\delta \rho _\alpha /\rho _0)}}{\mathrm{\partial} {\lambda }}{\boldsymbol e}_\lambda }\right) \times \left({ \frac{\mathrm{\partial} {M}}{\mathrm{\partial} {r}}{\boldsymbol e}_r + \frac{1}{r}\frac{\mathrm{\partial} {M}}{\mathrm{\partial} {\theta }}{\boldsymbol e}_\theta }\right)\nonumber \\ &=& -\frac{1}{r^2\sin \theta }\frac{\mathrm{\partial} {M}}{\mathrm{\partial} {\theta }}\frac{\mathrm{\partial} {(\delta \rho _\alpha /\rho _0)}}{\mathrm{\partial} {\lambda }}{\boldsymbol e}_r + \frac{1}{r\sin \theta }\frac{\mathrm{\partial} {M}}{\mathrm{\partial} {r}}\frac{\mathrm{\partial} {(\delta \rho _\alpha /\rho _0)}}{\mathrm{\partial} {\lambda }}{\boldsymbol e}_\theta + \frac{1}{r} \left({ \frac{\mathrm{\partial} {M}}{\mathrm{\partial} {\theta }}\frac{\mathrm{\partial} {(\delta \rho _\alpha /\rho _0)}}{\mathrm{\partial} {r}} -\frac{\mathrm{\partial} {M}}{\mathrm{\partial} {r}}\frac{\mathrm{\partial} {(\delta \rho _\alpha /\rho _0)}}{\mathrm{\partial} {\theta }} }\right){\boldsymbol e}_\lambda . \end{eqnarray}

(112)

We now have to extract the poloidal and toroidal contributions to this vector. To do so, first note that a toroidal field cannot have any component in the |${\boldsymbol e}_r$|-direction, and therefore the |${\boldsymbol e}_r$| component of any solenoidal vector |${\boldsymbol S}$| must come from the vector's poloidal piece. Next, we use vector identities to identify the relation between the poloidal-field scalar ϒ and the r-component of |${\boldsymbol S}$|⁠:

\begin{eqnarray} \nabla \times (\nabla \Upsilon \times {\boldsymbol r}) &=& (\nabla \cdot {\boldsymbol r})\nabla \Upsilon -{\boldsymbol r}\nabla ^2\Upsilon +({\boldsymbol r}\cdot \nabla )\nabla \Upsilon -(\nabla \Upsilon \cdot \nabla ){\boldsymbol r}\nonumber \\ &=& 2\nabla \Upsilon -{\boldsymbol r}\nabla ^2\Upsilon +r\frac{\mathrm{\partial} { }}{\mathrm{\partial} {r}}(\nabla \Upsilon ), \end{eqnarray}

(113)

and therefore

\begin{equation} S_r = \big [ 2\nabla \Upsilon -{\boldsymbol r}\nabla ^2\Upsilon +r\frac{\mathrm{\partial} { }}{\mathrm{\partial} {r}}(\nabla \Upsilon )\big ]_r = -\frac{1}{r\sin \theta }\frac{\mathrm{\partial} { }}{\mathrm{\partial} {\theta }}\left({\sin \theta \frac{\mathrm{\partial} {\Upsilon }}{\mathrm{\partial} {\theta }}}\right) - \frac{1}{r\sin ^2\theta }\frac{\mathrm{\partial} {^2\Upsilon }}{\mathrm{\partial} {\lambda ^2}}. \end{equation}

(114)

Decomposing ϒ and S_r into sums of spherical harmonics, |$\Upsilon = \sum _{l,m} \Upsilon _l^m Y_l^m\,,\ S_r = \sum _{l,m} G_l^m Y_l^m$|⁠, we can then identify the spherical-harmonic components of ϒ from those of S_r by using

\begin{equation} S_r = -\sum \limits _{l = 0}^\infty \sum \limits _{m = -l}^l \Upsilon _l^m(r) r\nabla ^2 Y_l^m = \sum \limits _{l = 0}^\infty \sum \limits _{m = -l}^l \frac{l(l+1) \Upsilon _l^m(r)}{r} Y_l^m, \end{equation}

(115)

i.e.

\begin{equation} \Upsilon _l^m = \frac{r}{l(l+1)} G_l^m\ \ (\forall \ l,m). \end{equation}

(116)

To apply these results to the vector ∇(δρ_α/ρ₀) × ∇M, we need an explicit expression for its r-component. Using the expressions from (57), (59), (95), we have

\begin{equation} \Big [\nabla \left({\frac{\delta \rho _\alpha }{\rho _0}}\right)\times \nabla M \Big ]_r = -\frac{1}{r^2\sin \theta }\frac{\mathrm{\partial} {M}}{\mathrm{\partial} {\theta }}\frac{\mathrm{\partial} {(\delta \rho _\alpha /\rho _0)}}{\mathrm{\partial} {\lambda }} = \frac{5\Lambda ^2\alpha ^2}{32G} j_2\left({{ {\frac{\pi r}{R_*}}}}\right) \sin (2\theta ) \big ( \sin (2\chi )\cos \theta \sin \lambda + \sin ^2\chi \sin \theta \sin (2\lambda ) \big ). \end{equation}

(117)

In order to use (116) to find the radial functions |$\Upsilon _l^m$| of the poloidal-field scalar, we must convert the angular terms in the above expression to spherical harmonics:

\begin{eqnarray} \sin (2\theta )\cos \theta \sin \lambda = 2\left(\sin \theta -\sin ^3\theta \right)\sin \lambda &=& -\frac{1}{5i}\left( {\frac{2}{3}}P_3^1+P_1^1\right) \left(\mathrm{e}^{i\lambda }-\mathrm{e}^{-i\lambda }\right)\nonumber \\ &=& -\frac{1}{5i} \left( P_1^1\mathrm{e}^{i\lambda }+2P_1^{-1}\mathrm{e}^{-i\lambda } + {\frac{2}{3}}P_3^1\mathrm{e}^{i\lambda }+8P_3^{-1}\mathrm{e}^{-i\lambda }\right)\nonumber \\ &=& -\frac{2\sqrt{\pi }}{5i}\left({ {\sqrt{\frac{2}{3}}}(Y_1^1+Y_1^{-1}) + {\frac{4}{\sqrt{21}}}(Y_3^1+Y_3^{-1}) }\right). \end{eqnarray}

(118)

\begin{eqnarray} \sin (2\theta )\sin \theta \sin (2\lambda ) &=& -\frac{1}{15i}P_3^2(\mathrm{e}^{2i\lambda }-\mathrm{e}^{-2i\lambda })\nonumber \\ &=& -\frac{1}{15i} \left(P_3^2\mathrm{e}^{2i\lambda }-120P_3^{-2}\mathrm{e}^{-2i\lambda }\right)\nonumber \\ &=& \frac{4\sqrt{\pi }}{15i}\sqrt{\frac{30}{7}}\left( Y_3^2-Y_3^{-2}\right). \end{eqnarray}

(119)

Substituting these relations into (117), we can use (116) to show that

\begin{eqnarray} \Upsilon _1^1 = \Upsilon _1^{-1} = \frac{i\pi ^{1/2}\Lambda ^2\alpha ^2\sin (2\chi )}{16\sqrt{6}G} r j_2\left({\frac{\pi r}{R_*}}\right), \end{eqnarray}

(120)

\begin{eqnarray} \Upsilon _3^1 = \Upsilon _3^{-1 } = \frac{i\pi ^{1/2}\Lambda ^2\alpha ^2\sin (2\chi )}{48\sqrt{21}G} r j_2\left({\frac{\pi r}{R_*}}\right), \end{eqnarray}

(121)

\begin{eqnarray} \Upsilon _3^2 = -\Upsilon _3^{-2} = -\frac{i\pi ^{1/2}\sqrt{10}\Lambda ^2\alpha ^2\sin ^2\chi }{96\sqrt{21}G} r j_2\left({\frac{\pi r}{R_*}}\right). \end{eqnarray}

(122)

Finally, then, the scalar function ϒ is given by

\begin{eqnarray} \Upsilon &=& \Upsilon _1^1 Y_1^1 + \Upsilon _1^{-1} Y_1^{-1} + \Upsilon _3^1 Y_3^1 + \Upsilon _3^{-1} Y_3^{-1} + \Upsilon _3^2 Y_3^2 + \Upsilon _3^{-2} Y_3^{-2} \nonumber \\ &=& \frac{5\alpha ^2\Lambda ^2}{192G} r j_2\left({ {\frac{\pi r}{R_*}}}\right) \left({ \sin (2\chi )\sin \theta (1+\cos ^2\theta )\sin \lambda + \sin ^2\chi \sin ^2\theta \cos \theta \sin (2\lambda ) }\right) + \mathcal {A}, \end{eqnarray}

(123)

where |$\mathcal {A}$| represents the gauge freedom in the solution; see Section 5.1. From this, we may now calculate the complete poloidal vector |$\nabla \times (\nabla \Upsilon \times {\boldsymbol r})$|⁠. The toroidal vector is then given by the difference |${\boldsymbol T}$| between the full expression (112) and the poloidal vector:

\begin{equation} \nabla \Psi \times {\boldsymbol r}= {\boldsymbol T} \equiv \nabla \left({\frac{\delta \rho _\alpha }{\rho _0}}\right)\times \nabla M - \nabla \times (\nabla \Upsilon \times {\boldsymbol r}), \end{equation}

(124)

where here – and for the rest of the derivation – we suppress the lengthy explicit expressions for ∇(δρ_α/ρ₀) × ∇M and its poloidal component. The θ and λ components of this equation give us expressions for Ψ:

\begin{eqnarray} \frac{1}{\sin \theta }\frac{\mathrm{\partial} {\Psi }}{\mathrm{\partial} {\lambda }} = T_\theta \Rightarrow \Psi = \int \sin \theta T_\theta \ \mathrm{d}\lambda + f(r,\theta ), \end{eqnarray}

(125)

\begin{eqnarray} -\frac{\mathrm{\partial} {\Psi }}{\mathrm{\partial} {\theta }} = T_\lambda \Rightarrow \Psi = -\int T_\lambda \ \mathrm{d}\theta + g(r,\lambda ). \end{eqnarray}

(126)

From this, we identify a function h(r, θ, λ) common to both expressions for Ψ, and deduce that the full expression for Ψ is

\begin{equation} \Psi = h(r,\theta ,\lambda )+f(r,\theta )+g(r,\lambda ). \end{equation}

(127)

After considerable rearrangement and simplification, we find that Ψ is given by

\begin{eqnarray} \Psi &=& \frac{5\alpha ^2\Lambda ^2}{384G} \Big [ \left({ {\frac{\pi r}{R_*}}}j_1-j_2\right) (2-3\sin ^2\theta ) + 6\left({1+{ {\frac{\pi r}{R_*}}}\cot \left({ {\frac{\pi r}{R_*}}}\right)}\right)j_2\sin ^2\theta \Big ] \left(\sin (2\chi )\sin (2\theta )\cos \lambda +\sin ^2\chi \sin ^2\theta \cos (2\lambda )\right)\nonumber \\ &&+\, \frac{5\alpha ^2\Lambda ^2\sin ^2\chi }{3072G} \left(j_2-{ {\frac{\pi r}{R_*}}}j_1\right) + \mathcal {B}. \end{eqnarray}

(128)

We now exploit the ‘gauge freedom’ of Ψ by choosing the function |$\mathcal {B}$| as follows:

\begin{equation} \mathcal {B} = -\frac{5\alpha ^2\Lambda ^2\sin ^2\chi }{3072G}\left(j_2-{ {\frac{\pi r}{R_*}}}j_1\right). \end{equation}

(129)

This clearly satisfies the requirement that |$\nabla \mathcal {B}\times {\boldsymbol r}= 0$| (from Section 5.1), since |$\nabla \mathcal {B}$| must be parallel to |${\boldsymbol r}$|⁠. Therefore, we may add this on to the original expression for Ψ, cancelling the final term:

\begin{equation} \Psi +\mathcal {B} = \frac{5\alpha ^2\Lambda ^2}{384G} \Big [ \left({ {\frac{\pi r}{R_*}}}j_1{-}j_2\right)(2{-}3\sin ^2\theta ) + 6\left({1+{ {\frac{\pi r}{R_*}}}\cot \left({ {\frac{\pi r}{R_*}}}\right)}\right)j_2\sin ^2\theta \Big ] (\sin (2\chi )\sin (2\theta )\cos \lambda +\sin ^2\chi \sin ^2\theta \cos (2\lambda ))+\mathcal {C}. \end{equation}

(130)

Note that one is still free to add additional terms |$\mathcal {C}$| satisfying the condition |$\nabla \mathcal {C}\times {\boldsymbol r}= 0$|⁠. To summarize, we have found scalar functions – ϒ from (123) and |$\Psi +\mathcal {B}$| from (130) – representing the two degrees of freedom of |$\nabla \times [\delta \rho _\alpha (\nabla \times {\boldsymbol B}_0)\times {\boldsymbol B}_0/4\pi \rho _0^2]$|⁠, one of the three pieces of the perturbed Lorentz force. We can now replace the schematic poloidal and toroidal pieces of this quantity in equations (101) and (102) with explicit results in terms of scalar functions.

The final terms in the curled Euler equations are the two other pieces of the perturbed Lorentz force. Again, we know that the curl of these force terms must be solenoidal and therefore expressible in terms of another pair of poloidal and toroidal scalars:

\begin{equation} \nabla \times \Big \lbrace \frac{(\nabla \times \delta {\boldsymbol B})\times {\boldsymbol B}_0+(\nabla \times {\boldsymbol B}_0)\times \delta {\boldsymbol B}}{4\pi \rho _0}\Big \rbrace = \nabla \times (\Psi _{\delta \mathcal {L}}{\boldsymbol r}) + \nabla \times (\nabla \Upsilon _{\delta \mathcal {L}}\times {\boldsymbol r}). \end{equation}

(131)

These scalars involve the unknown |$\delta {\boldsymbol B}$|⁠, but they can none the less now be calculated explicitly, since we know all the other terms in equations (101) and (102). Equating the terms under the curl, using the results (110), (123) and (130),

\begin{eqnarray} \Upsilon _{\delta \mathcal {L}} &=& -\alpha \omega \sin \chi r\sin \theta \sin \lambda +\Upsilon \nonumber \\ &=& \frac{\alpha ^2\Lambda ^2\sin \chi }{G} r\sin \theta \sin \lambda \Big \lbrace {\frac{5}{192}}j_2\left({ {\frac{\pi r}{R_*}}}\right) \big [2\cos \chi (1+\cos ^2\theta )+\sin \chi \sin (2\theta )\cos \lambda \big ] -k_B \cos \chi \Big \rbrace + \mathcal {A}, \end{eqnarray}

(132)

\begin{eqnarray} \Psi _{\delta \mathcal {L}} &=& \Psi +\mathcal {B}\nonumber \\ &=& \frac{5\alpha ^2\Lambda ^2}{384G} \Big [ \left({ {\frac{\pi r}{R_*}}}j_1-j_2\right) (2-3\sin ^2\theta ) + 6\left({1+{ {\frac{\pi r}{R_*}}}\cot \left({ {\frac{\pi r}{R_*}}}\right)}\right)j_2\sin ^2\theta \Big ] \left(\sin (2\chi )\sin (2\theta )\cos \lambda +\sin ^2\chi \sin ^2\theta \cos (2\lambda )\right) + \mathcal {C},\nonumber\\ \end{eqnarray}

(133)

where we have used the scaling of ε_B taken from (62) to simplify the first equation.

5.3 Equations for |$\delta {\boldsymbol B}$|

With our expressions for the two scalar functions |$\Upsilon _{\delta \mathcal {L}}$| and |$\Psi _{\delta \mathcal {L}}$|⁠, we now have an explicit expression for the curl of the perturbed Lorentz force – see equation (131). Inverting the curl operator from this equation gives us an expression for the two terms involving the unknown |$\delta {\boldsymbol B}$|⁠, in terms of the known quantities |$\Upsilon _{\delta \mathcal {L}}$| and |$\Psi _{\delta \mathcal {L}}$|⁠. In doing so, however, we have to allow for the possibility of this force containing an additional, and unknown, irrotational term ∇Ξ:

\begin{equation} \frac{(\nabla \times \delta {\boldsymbol B})\times {\boldsymbol B}_0+(\nabla \times {\boldsymbol B}_0)\times \delta {\boldsymbol B}}{4\pi \rho _0} = \Psi _{\delta \mathcal {L}}{\boldsymbol r}+ \nabla \Upsilon _{\delta \mathcal {L}}\times {\boldsymbol r}+\nabla \Xi . \end{equation}

(134)

Note that the gauge-freedom functions |$\mathcal {A}$| and |$\mathcal {C}$| may be absorbed into ∇Ξ, since the addition of these terms to |$\Psi _{\delta \mathcal {L}}$| and |$\Upsilon _{\delta \mathcal {L}}$| produces the following contribution to the right-hand side of equation (134):

\begin{eqnarray} \mathcal {C}{\boldsymbol r}+ \nabla \mathcal {A}\times {\boldsymbol r}, \end{eqnarray}

(135)

and these terms are curl-free by the definitions of |$\mathcal {A}$| and |$\mathcal {C}$|⁠. If we compare (134) with the original second-order Euler equation (51), we see that Ξ must include all the information about the second-order fluid perturbations (e.g. δH_αB), which was thrown away by taking the curl of the Euler equation.

Equation (134) gives three scalar equations in four unknowns – the three components of |$\delta {\boldsymbol B}$| and the new scalar Ξ. The set of equations is then closed with the divergence-free condition on |$\delta {\boldsymbol B}$|⁠. We are, at last, in a position to obtain a set of DEs featuring |$\delta {\boldsymbol B}$| explicitly.

Our equations for the perturbed field involve vector operations on |$\delta {\boldsymbol B}$| and |${\boldsymbol B}_0$|⁠. Operations of this form are generally simplest, performed in a vector-spherical-harmonic basis, and so we write |$\delta {\boldsymbol B}$| as an infinite sum over these basis vectors, and |${\boldsymbol B}_0$| in its (known) vector-spherical-harmonic form:

\begin{eqnarray} \delta {\boldsymbol B} = \sum \limits _{l = 1}^\infty \sum \limits _{m = -l,m\ne 0}^{l} \left({ U_l^mY_l^m{\boldsymbol e}_r + V_l^m\nabla Y_l^m + W_l^m{\boldsymbol e}_r\times \nabla Y_l^m }\right), \end{eqnarray}

(136)

\begin{eqnarray} {\boldsymbol B}_0 = -2\sqrt{\frac{\pi }{3}}\Lambda r^2\rho _0 {\boldsymbol e}_r\times \nabla Y_1^0. \end{eqnarray}

(137)

Recall that we are looking for time-dependent solutions, ∼e^imλ = e^{i(mϕ + ωt)}, so we explicitly exclude m = 0 terms from the sum for |$\delta {\boldsymbol B}$| (though our reasoning is applicable to this case too). In what follows, we will use ∑_{l, m} as shorthand for the summation used for |$\delta {\boldsymbol B}$| above. First, we look at the divergence-free condition for |$\delta {\boldsymbol B}$| with the above decomposition:

\begin{eqnarray} 0 = \nabla \cdot \delta {\boldsymbol B}&=&\sum \limits _{l,m} \left\lbrace \nabla (U_l^mY_l^m)\cdot {\boldsymbol e}_r + U_l^m Y_l^m\nabla \cdot {\boldsymbol e}_r + \nabla V_l^m\cdot \nabla Y_l^m+V_l^m\nabla ^2 Y_l^m +\nabla \cdot [W_l^m\nabla \times (Y_l^m{\boldsymbol e}_r)]\right\rbrace \nonumber \\ &=&\sum \limits _{l,m} \left[ \nabla (U_l^mY_l^m)\cdot {\boldsymbol e}_r + U_l^m Y_l^m\nabla \cdot {\boldsymbol e}_r + V_l^m\nabla ^2 Y_l^m \right] \nonumber \\ &=&\sum \limits _{l,m} \left({U_l^m{}^{\prime } + \frac{2}{r}U_l^m - \frac{l(l+1)}{r^2}V_l^m}\right)Y_l^m, \end{eqnarray}

(138)

where |$U_l^m{}^{\prime } = (U_l^m)^{\prime } = \mathrm{d}U_l^m/\mathrm{d}r$|⁠. Now, equation (138) implies that

\begin{equation} V_l^m = \frac{r^2 U_l^m {}^{\prime } + 2rU_l^m}{l(l+1)}. \end{equation}

(139)

Note that this result shows how a solenoidal vector field loses one degree of freedom. Next, we rewrite the curl of |$\delta {\boldsymbol B}$|⁠, using standard vector identities:

\begin{eqnarray} \nabla \times \delta {\boldsymbol B}&=&\sum \limits _{l,m} \left[ \nabla (U_l^m Y_l^m)\times {\boldsymbol e}_r + \nabla V_l^m\times \nabla Y_l^m + \nabla W_l^m\times ({\boldsymbol e}_r\times \nabla Y_l^m) + W_l^m\nabla \times ({\boldsymbol e}_r\times \nabla Y_l^m)\right]\nonumber \\ &=& \sum \limits _{l,m} \left\lbrace U_l^m\nabla Y_l^m\times {\boldsymbol e}_r + V_l^m{}^{\prime }{\boldsymbol e}_r\times \nabla Y_l^m - \nabla Y_l^m(\nabla W_l^m\cdot {\boldsymbol e}_r)\right.\nonumber \\ \left.+\, W_l^m\left[{\boldsymbol e}_r\nabla ^2 Y_l^m-\nabla Y_l^m(\nabla \cdot {\boldsymbol e}_r) + (\nabla Y_l^m\cdot \nabla ){\boldsymbol e}_r-({\boldsymbol e}_r\cdot \nabla )\nabla Y_l^m\right]\right\rbrace \nonumber \\ &=& \sum \limits _{l,m} \left[ -\frac{l(l+1)W_l^m}{r^2} Y_l^m{\boldsymbol e}_r - W_l^m{}^{\prime }\nabla Y_l^m + (V_l^m{}^{\prime }-U_l^m){\boldsymbol e}_r\times \nabla Y_l^m \right]. \end{eqnarray}

(140)

Comparing this result from earlier ones in this section, we see that this is an embodiment of the known result that the curl of a poloidal field is toroidal, and vice versa. We use the above expression to calculate one of the unknown pieces of the perturbed Lorentz force:

\begin{equation} (\nabla \times \delta {\boldsymbol B})\times {\boldsymbol B}_0 = -\Lambda \rho _0 \sum \limits _{l,m} \bigg \lbrace \big [ W_l^m{}^{\prime }\sin \theta Y_l^m{}_{,\theta } +(V_l^m{}^{\prime }-U_l^m)Y_l^m{}_{,\lambda } \big ]{\boldsymbol e}_r + l(l+1)W_l^mY_l^m\nabla (\cos \theta ) \bigg \rbrace , \end{equation}

(141)

where we have used the following vector identities

\begin{eqnarray} {\boldsymbol e}_r\times ({\boldsymbol e}_r\times \nabla Y_1^0) = -\nabla Y_1^0= -\frac{\sqrt{3/\pi }}{2}\nabla (\cos \theta ), \end{eqnarray}

(142)

\begin{eqnarray} \nabla Y_l^m\times ({\boldsymbol e}_r\times \nabla Y_1^0) = {\boldsymbol e}_r \nabla Y_l^m\cdot \nabla Y_1^0 = -\frac{\sqrt{3/\pi }}{2r^2}\sin \theta Y_l^m{}_{,\theta }{\boldsymbol e}_r, \end{eqnarray}

(143)

\begin{eqnarray} ({\boldsymbol e}_r\times \nabla Y_l^m)\times ({\boldsymbol e}_r\times \nabla Y_1^0) = \nabla Y_l^m\times \nabla Y_1^0 = \frac{\sqrt{3/\pi }}{2r^2} Y_l^m{}_{,\lambda }{\boldsymbol e}_r. \end{eqnarray}

(144)

Very similar relations may be used to calculate the other unknown piece of the perturbed Lorentz force:

\begin{eqnarray} (\nabla \times {\boldsymbol B}_0)\times \delta {\boldsymbol B}&=& -\Lambda \rho _0 \sum \limits _{l,m} \left[ (V_l^m Y_l^m{}_{,\lambda } + W_l^m\sin \theta Y_l^m{}_{,\theta })\left({\frac{2}{r}+\frac{\rho ^{\prime }_0}{\rho _0}}\right){\boldsymbol e}_r +2W_l^m\cos \theta \nabla Y_l^m\right.\nonumber \\ \, \left.+\, \left({2r+\frac{r^2\rho ^{\prime }_0}{\rho _0}}\right)U_l^mY_l^m{\boldsymbol e}_r\times \nabla (\cos \theta ) - 2V_l^m\cos \theta {\boldsymbol e}_r\times \nabla Y_l^m \right]. \end{eqnarray}

(145)

We are now in a position to reduce our original system of equations – (134) plus |$\nabla \cdot \delta {\boldsymbol B}= 0$| – to three equations (per l, m) in the unknown radial functions |$U_l^m$|⁠, |$W_l^m$| and |$\Xi _l^m$|⁠, where the latter is defined through |$\Xi = \sum _{l,m} \Xi _l^m Y_l^m$|⁠. To find this simpler set of equations, we take the expressions for the perturbed Lorentz force from equations (141) and (145), and use the divergence-free condition (138) to eliminate all |$V_l^m$| in favour of |$U_l^m$|⁠. We also rewrite the right-hand side of (134) slightly so that it is manifestly in vector-spherical-harmonic form. The resulting set of equations is

\begin{eqnarray} 0&=& \sum \limits _{l,m}\bigg \lbrace \frac{4\pi }{\Lambda }\big [ (r\tilde{\Psi }_l^m+\Xi _l^m{}^{\prime })Y_l^m{\boldsymbol e}_r + \Xi _l^m\nabla Y_l^m - r\tilde{\Upsilon }_l^m{\boldsymbol e}_r\times \nabla Y_l^m \big ] + \left[ W_l^m{}^{\prime }+\left({\frac{2}{r}+\frac{\rho ^{\prime }_0}{\rho _0}}\right)W_l^m \right]\sin \theta Y_l^m{}_{,\theta }{\boldsymbol e}_r \nonumber \\ \, +\, \left[ r^2 U_l^m{}^{\prime \prime } + \left({6+\frac{r\rho ^{\prime }_0}{\rho _0}}\right)rU_l^m{}^{\prime } +\left({6-l(l+1)+\frac{2r\rho ^{\prime }_0}{\rho _0}}\right)U_l^m \right]\frac{Y_l^m{}_{,\lambda }}{l(l+1)}{\boldsymbol e}_r +\left({2r+\frac{r^2\rho ^{\prime }_0}{\rho _0}}\right)U_l^mY_l^m{\boldsymbol e}_r\times \nabla (\cos \theta ) \nonumber \\ \,\, -\, \frac{2}{l(l+1)}(r^2 U_l^m{}^{\prime } + 2rU_l^m)\cos \theta {\boldsymbol e}_r\times \nabla Y_l^m +2W_l^m\cos \theta \nabla Y_l^m + l(l+1)W_l^m Y_l^m\nabla (\cos \theta ) \bigg \rbrace . \end{eqnarray}

(146)

In the above, |$\tilde{\Psi }_l^m$| and |$\tilde{\Upsilon }_l^m$| are the radial functions in spherical-harmonic decompositions of |$\Psi _{\delta \mathcal {L}}$| and |$\Upsilon _{\delta \mathcal {L}}$|⁠:

\begin{equation} \Psi _{\delta \mathcal {L}} = \sum _{l,m} \tilde{\Psi }_l^m(r)Y_l^m\,,\ \Upsilon _{\delta \mathcal {L}} = \sum _{l,m} \tilde{\Upsilon }_l^m(r)Y_l^m, \end{equation}

(147)

defined with tildes to distinguish them from the |$\Psi _l^m,\Upsilon _l^m$| defined in the previous section. The explicit form of the functions |$\tilde{\Psi }_l^m,\tilde{\Upsilon }_l^m$| is given in Appendix A.

Note that equation (146) is not expressed in a vector-spherical-harmonic basis, as it contains product terms like |$Y_l^m\nabla (\cos \theta )$|⁠. Writing a term like this in a vector-spherical-harmonic basis appears to require an infinite sum in itself, which would have to be performed for each l, m in the original (infinite) sum. Instead of facing this unappealing double infinite sum, we look at the spherical polar coordinate components of equation (146). Our aim will be to convert these into three ODEs, per l, m, in the three unknown radial functions. As we will see later, the angular dependence can be removed using the spherical harmonic orthogonality relation, but to do so, we need to cast all angular dependence of the components of (146) in the form of spherical harmonics. To this end, we will need three key identities:

\begin{eqnarray} \sin \theta Y_l^m{}_{,\theta } = lQ_{l+1}Y_{l+1}^m - (l+1)Q_lY_{l-1}^m , \end{eqnarray}

(148)

\begin{eqnarray} \cos \theta Y_l^m = Q_{l+1}Y_{l+1}^m + Q_l Y_{l-1}^m , \end{eqnarray}

(149)

\begin{eqnarray} Y_l^m{}_{,\lambda } = im Y_l^m, \end{eqnarray}

(150)

where we have defined

\begin{equation} Q_l\equiv \sqrt{\frac{(l-m)(l+m)}{(2l-1)(2l+1)}}. \end{equation}

(151)

Using these relations, it is clear that all angular terms in the r-component of equation (146) may be rewritten accordingly. The θ- and λ-components of the equation, however, involve combinations like |$\cos \theta Y_l^m{}_{,\theta }$|⁠, which cannot be readily simplified (we encounter a similar problem as for the vector-spherical-harmonic basis, where one needs to invoke an infinite sum). To avoid this problem, we leave the r-component of (146) unchanged, but multiply its θ and λ components by sin θ, giving us the set:

\begin{eqnarray} 0 &=& \sum \limits _{l,m}\bigg \lbrace \frac{4\pi }{\Lambda } (r\tilde{\Psi }_l^m+\Xi _l^m{}^{\prime }) Y_l^m + \frac{im}{l(l+1)}\left[ r^2 U_l^m{}^{\prime \prime } + \left({6+\frac{r\rho ^{\prime }_0}{\rho _0}}\right)rU_l^m{}^{\prime } + \left({6-l(l+1)+\frac{2r\rho ^{\prime }_0}{\rho _0}}\right)U_l^m \right] Y_l^m \nonumber \\ \, +\,\left[ W_l^m{}^{\prime } + \left({2+\frac{r\rho ^{\prime }_0}{\rho _0}}\right)\frac{W_l^m}{r} \right]\sin \theta Y_l^m{}_{,\theta } \bigg \rbrace , \end{eqnarray}

(152)

\begin{eqnarray} 0 = \sum \limits _{l,m}\bigg \lbrace\frac{4\pi }{\Lambda } \left({\frac{\Xi _l^m}{r}\sin \theta Y_l^m{}_{,\theta } + im\tilde{\Upsilon }_l^m Y_l^m}\right) + \frac{2im}{l(l+1)} (r U_l^m{}^{\prime } +2U_l^m) \cos \theta Y_l^m + [2\cos \theta \sin \theta Y_l^m{}_{,\theta }-l(l+1)\sin ^2\theta Y_l^m]\frac{W_l^m}{r} \bigg \rbrace , \end{eqnarray}

(153)

\begin{eqnarray} 0 &=& \sum \limits _{l,m}\bigg \lbrace \frac{4\pi }{\Lambda } \left({ \frac{im\Xi _l^m}{r} Y_l^m - \tilde{\Upsilon }_l^m\sin \theta Y_l^m{}_{,\theta }}\right) - \left({2+\frac{r\rho ^{\prime }_0}{\rho _0}}\right)U_l^m\sin ^2\theta Y_l^m -\frac{2}{l(l+1)}(rU_l^m{}^{\prime }+2U_l^m)\cos \theta \sin \theta Y_l^m{}_{,\theta }+ \frac{2im W_l^m}{r}\cos \theta Y_l^m \bigg \rbrace .\nonumber\\ \end{eqnarray}

(154)

The latter two equations in the above set may now also be rewritten, by the use of the following auxiliary relations (which may be derived from 148 and 149):

\begin{eqnarray} \cos \theta \sin \theta Y_l^m{}_{,\theta } = l Q_{l+1}Q_{l+2}Y_{l+2}^m + [lQ_{l+1}^2-(l+1)Q_l^2]Y_l^m - (l+1)Q_{l-1}Q_lY_{l-2}^m, \end{eqnarray}

(155)

\begin{eqnarray} \sin ^2\theta Y_l^m = Y_l^m-\cos ^2\theta Y_l^m = -Q_{l+1}Q_{l+2}Y_{l+2}^m+ (1-Q_l^2-Q_{l+1}^2)Y_l^m - Q_{l-1}Q_lY_{l-2}^m. \end{eqnarray}

(156)

With these relations, we are now able to eliminate all instances of trigonometric functions and derivatives of |$Y_l^m$| in our set of equations:

\begin{eqnarray} 0 &=& \sum \limits _{l,m} \bigg \lbrace \frac{4\pi }{\Lambda } (r\tilde{\Psi }_l^m+\Xi _l^m{}^{\prime }) Y_l^m + \frac{im}{l(l+1)}\left[ r^2 U_l^m{}^{\prime \prime } + \left({6+\frac{r\rho ^{\prime }_0}{\rho _0}}\right)rU_l^m{}^{\prime } + \left({6-l(l+1)+\frac{2r\rho ^{\prime }_0}{\rho _0}}\right)U_l^m \right] Y_l^m \nonumber \\ \, + \, \left[ W_l^m{}^{\prime } + \left({2+\frac{r\rho ^{\prime }_0}{\rho _0}}\right)\frac{W_l^m}{r} \right] \Big [ lQ_{l+1}Y_{l+1}^m - (l+1)Q_lY_{l-1}^m \Big ] \bigg \rbrace , \end{eqnarray}

(157)

\begin{eqnarray} 0&=& \sum \limits _{l,m} \bigg \lbrace\frac{4\pi }{\Lambda } \left\lbrace \frac{\Xi _l^m}{r}[lQ_{l+1}Y_{l+1}^m-(l+1)Q_l Y_{l-1}^m] + im\tilde{\Upsilon }_l^m Y_l^m \right\rbrace + \frac{2im}{l(l+1)} (r U_l^m{}^{\prime }+2U_l^m)(Q_{l+1}Y_{l+1}^m+Q_l Y_{l-1}^m) \nonumber \\ \, +\, \frac{W_l^m}{r}\Big \lbrace l(l+3)Q_{l+1}Q_{l+2}Y_{l+2}^m + \big [ l(l+3)Q_{l+1}^2 + (l-2)(l+1)Q_l^2 - l(l+1) \big ]Y_l^m + (l-2)(l+1)Q_{l-1}Q_l Y_{l-2}^m\Big \rbrace \bigg \rbrace , \end{eqnarray}

(158)

\begin{eqnarray} 0&=& \sum \limits _{l,m} \bigg \lbrace\frac{4\pi }{\Lambda } \left\lbrace \frac{im\Xi _l^m}{r} Y_l^m - \tilde{\Upsilon }_l^m[lQ_{l+1}Y_{l+1}^m-(l+1)Q_l Y_{l-1}^m] \right\rbrace + U_l^m\bigg \lbrace \left({\frac{2(l-1)}{(l+1)}+\frac{r\rho ^{\prime }_0}{\rho _0} }\right)Q_{l+1}Q_{l+2} Y_{l+2}^m\nonumber \\ \, -\, \left[ 2\left({1-\frac{(l+2)Q_l^2}{l}-\frac{(l-1)Q_{l+1}^2}{(l+1)}}\right) + (1-Q_l^2-Q_{l+1}^2)\frac{r\rho ^{\prime }_0}{\rho _0} \right]Y_l^m + \left({\frac{2(l+2)}{l}+\frac{r\rho ^{\prime }_0}{\rho _0}}\right)Q_{l-1}Q_l Y_{l-2}^m\bigg \rbrace \nonumber \\ \, -\, \frac{2}{l(l+1)} r U_l^m{}^{\prime }\Big \lbrace lQ_{l+1}Q_{l+2}Y_{l+2}^m + [lQ_{l+1}^2-(l+1)Q_l^2]Y_l^m -(l+1)Q_{l-1}Q_l Y_{l-2}^m \Big \rbrace + \frac{2im W_l^m}{r} (Q_{l+1}Y_{l+1}^m+Q_lY_{l-1}^m) \bigg \rbrace . \end{eqnarray}

(159)

Note that in the above, spherical-harmonic terms with l < |m| are automatically excluded from the sum by the form of their Q_l pre-factors (which are zero for l = ±m). Next, we will eliminate all the angular terms in the above, by using the spherical-harmonic orthogonality relation:

\begin{equation} \int Y_{\hat{l}}^{\hat{m}} Y_l^m{}^*\ \sin \theta \mathrm{d}\theta \mathrm{d}\lambda = \delta _{\hat{l}l} \delta _{\hat{m}m}. \end{equation}

(160)

Let us relabel the indices l, m in equations (157), (158) and (159) to |$\hat{l},\hat{m}$|⁠, to match the notation of the above relation. Multiplying these three equations by |$Y_l^m{}^*$|⁠, and integrating and using the above relation, then picks out individual elements in the sums. This reduces the three infinite sums to three equalities (without summation) for each value of m ≠ 0 and l ≥ 1. These equalities are now just DEs in the radial coordinate, as the angular dependence has dropped out:

\begin{eqnarray} 0 &=& \frac{4\pi }{\Lambda } (r\tilde{\Psi }_l+\Xi ^{\prime }_l) + \frac{im}{l(l+1)}\left[ r^2 U^{\prime \prime }_l + \left({6+\frac{r\rho ^{\prime }_0}{\rho _0}}\right)rU^{\prime }_l + \left({6-l(l+1)+\frac{2r\rho ^{\prime }_0}{\rho _0}}\right)U_l \right]\nonumber \\ + \,(l-1)Q_l \left[ W^{\prime }_{l-1}+\left({2+\frac{r\rho ^{\prime }_0}{\rho _0}}\right)\frac{W_{l-1}}{r} \right] - (l+2)Q_{l+1} \left[ W^{\prime }_{l+1}+\left({2+\frac{r\rho ^{\prime }_0}{\rho _0}}\right)\frac{W_{l+1}}{r} \right], \end{eqnarray}

(161)

\begin{eqnarray} 0 &=& \frac{4\pi }{\Lambda } \left[ (l-1)Q_l\frac{\Xi _{l-1}}{r}-(l+2)Q_{l+1}\frac{\Xi _{l+1}}{r} + im\tilde{\Upsilon }_l\right] + \frac{2imQ_l}{(l-1)l} (rU^{\prime }_{l-1}+2U_{l-1}) + \frac{2imQ_{l+1}}{(l+1)(l+2)}(rU^{\prime }_{l+1}+2U_{l+1})\nonumber \\ +\, (l-2)(l+1)Q_{l-1}Q_l \frac{W_{l-2}}{r} + \big [l(l+3)Q_{l+1}^2+(l-2)(l+1)Q_l^2-l(l+1)\big ] \frac{W_l}{r} + l(l+3)Q_{l+1}Q_{l+2}\frac{W_{l+2}}{r}, \end{eqnarray}

(162)

\begin{eqnarray} 0&=& \frac{4\pi }{\Lambda }\left[ \frac{im\Xi _l}{r} - (l-1)Q_l \tilde{\Upsilon }_{l-1}+(l+2)Q_{l+1}\tilde{\Upsilon }_{l+1} \right] - \frac{2}{(l-1)}Q_{l-1}Q_l r U^{\prime }_{l-2} + Q_{l-1}Q_l\left({ \frac{2(l-3)}{(l-1)}+\frac{r\rho ^{\prime }_0}{\rho _0} }\right) U_{l-2}\nonumber \\ -\, \frac{2}{l(l+1)} [lQ_{l+1}^2-(l+1)Q_l^2] r U^{\prime }_l - \bigg [ 2\left({1-\frac{(l+2)Q_l^2}{l}-\frac{(l-1)Q_{l+1}^2}{(l+1)}}\right) + (1-Q_l^2-Q_{l+1}^2)\frac{r\rho ^{\prime }_0}{\rho _0} \bigg ] U_l\nonumber \\ +\, \frac{2Q_{l+1}Q_{l+2}}{(l+2)}r U^{\prime }_{l+2} + Q_{l+1}Q_{l+2}\left({ \frac{2(l+4)}{(l+2)}+\frac{r\rho ^{\prime }_0}{\rho _0} }\right)U_{l+2} + 2im Q_l \frac{W_{l-1}}{r} + 2imQ_{l+1}\frac{W_{l+1}}{r}. \end{eqnarray}

(163)

The above equations form an infinite set for different values of l and m. Individual equations couple together U and W functions with different angular indices l, but terms with different azimuthal index m decouple, allowing us to suppress m. Since the angular structure of the source terms ϒ, Ψ in equations (161)–(163) only includes terms with azimuthal index |m| = 1 and 2 (see Appendix A), we expect the solution to reflect this, i.e. |$U_l^m = W_l^m = 0$| for |m| ≥ 3 (and as before, we exclude the stationary m = 0 terms). In addition, in the case of an orthogonal rotator (χ = π/2), only the |m| = 2 terms survive.

For brevity, we have shown equations for general l, but by doing so, some of the above equations as written feature U and W functions with l < |m|, and in particular with negative l (e.g. the U_l − 2 terms when l = m = 1). The origin of these terms, however, is in l < |m| spherical harmonics from (157) to (159) – and these are identically zero (see the notes after these equations). Accordingly, equations (161)–(163) and those which follow them in this section are presented on the understanding that any instances of U or W functions with l < |m| are undefined and should be excluded. This is then also consistent with the original sum (136) in which the U and W functions first appeared.

To find the perturbed magnetic field |$\delta {\boldsymbol B}$|⁠, we need only the functions U_l and W_l, and are not interested in solving for the other unknown set of functions Ξ_l. We next use (163) to eliminate all instances of Ξ_l terms in equations (161) and (162). After simplification, the result is the following two coupled DEs for each value of l:

\begin{eqnarray} 0 &=& \frac{4\pi }{\Lambda } \left[ mr\tilde{\Psi }_l - i(l-1)Q_l(r\tilde{\Upsilon }^{\prime }_{l-1}+\tilde{\Upsilon }_{l-1}) + i(l+2)Q_{l+1}(r\tilde{\Upsilon }^{\prime }_{l+1}+\tilde{\Upsilon }_{l+1}) \right]\nonumber \\ +\, Q_{l-1}Q_l\bigg \lbrace -\frac{2}{(l-1)}r^2\bar{U}^{\prime \prime }_{l-2} + \left[ \frac{2(l-5)}{(l-1)}+\frac{r\rho ^{\prime }_0}{\rho _0}\right]r\bar{U}^{\prime }_{l-2} + \left[ \frac{2(l-3)}{(l-1)}+\frac{r\rho ^{\prime }_0}{\rho _0}+r\left({\frac{r\rho ^{\prime }_0}{\rho _0}}\right)^{\prime }\right]\bar{U}_{l-2}\bigg \rbrace \nonumber \\ +\, \frac{(m^2+2(l+1)Q_l^2-2lQ_{l+1}^2)}{l(l+1)}r^2\bar{U}^{\prime \prime }_l + \left[\frac{6m^2}{l(l+1)}-2+\frac{2(l+4)Q_l^2}{l}+\frac{2(l-3)Q_{l+1}^2}{(l+1)} +\left({\!\frac{m^2}{l(l+1)}\!-\!1\!+\!Q_l^2\!+\!Q_{l+1}^2\!}\right)\!\frac{r\rho ^{\prime }_0}{\rho _0} \right]r\bar{U}^{\prime }_l\nonumber \\ +\, \left[ \frac{6m^2}{l(l+1)}-m^2 -2+\frac{2(l+2)Q_l^2}{l}+\frac{2(l-1)Q_{l+1}^2}{(l+1)} + \left({\frac{2m^2}{l(l+1)}-1+Q_l^2+Q_{l+1}^2}\right)\frac{r\rho ^{\prime }_0}{\rho _0} - (1-Q_l^2-Q_{l+1}^2)r\left({\frac{r\rho ^{\prime }_0}{\rho _0}}\right)^{\prime } \right]\bar{U}_l \nonumber \\ +\, Q_{l+1}Q_{l+2}\bigg \lbrace \frac{2}{(l+2)}r^2\bar{U}^{\prime \prime }_{l+2} + \left[ \frac{2(l+6)}{(l+2)}+\frac{r\rho ^{\prime }_0}{\rho _0} \right]r\bar{U}^{\prime }_{l+2} + \left[ \frac{2(l+4)}{(l+2)}+\frac{r\rho ^{\prime }_0}{\rho _0}+r\left({\frac{r\rho ^{\prime }_0}{\rho _0}}\right)^{\prime } \right]\bar{U}_{l+2} \bigg \rbrace \nonumber \\ +\, (l-3)mQ_l rX^{\prime }_{l-1} + mQ_l\left[ (3l-5)+(l-1)\frac{r\rho ^{\prime }_0}{\rho _0} \right]X_{l-1} - (l+4)mQ_{l+1} rX^{\prime }_{l+1} - mQ_{l+1}\left[ (3l+8)+(l+2)\frac{r\rho ^{\prime }_0}{\rho _0} \right]X_{l+1}, \end{eqnarray}

(164)

\begin{eqnarray} 0 &=& \frac{4\pi i}{\Lambda }\left\lbrace -(l-1)(l-2)Q_{l-1}Q_l\tilde{\Upsilon }_{l-2} + [ m^2+(l-1)(l+1)Q_l^2+l(l+2)Q_{l+1}^2 ]\tilde{\Upsilon }_l - (l+2)(l+3)Q_{l+1}Q_{l+2}\tilde{\Upsilon }_{l+2} \right\rbrace \nonumber \\ -\, (l-1)Q_{l-2}Q_{l-1}Q_l\left[ \frac{2}{(l-2)}r\bar{U}^{\prime }_{l-3}-\left({\!\frac{2(l-4)}{(l-2)}+\frac{r\rho ^{\prime }_0}{\rho _0}\!}\right)\!\bar{U}_{l-3} \right] + \frac{Q_l}{l}\bigg \lbrace 2\left[ \frac{m^2}{(l\!-\!1)}+lQ_{l-1}^2-(l\!-\!1)Q_l^2+(l\!+\!2)Q_{l+1}^2 \right]r\bar{U}^{\prime }_{l-1} \nonumber \\ +\, \bigg [ \frac{4m^2}{(l\!-\!1)} -2l\big [l-1-(l\!+\!1)Q_{l-1}^2\big ] + 2(l\!-\!2)\big [(l\!-\!1)Q_l^2 \!-\! (l\!+\!2)Q_{l+1}^2\big ] - l\big [ (l\!-\!1)(1\!-\!Q_{l-1}^2\!-\!Q_l^2) + (l\!+\!2)Q_{l+1}^2 \big ]\frac{r\rho ^{\prime }_0}{\rho _0} \bigg ]\bar{U}_{l-1} \bigg \rbrace \nonumber \\ +\, \frac{Q_{l+1}}{(l+1)}\bigg \lbrace 2\left[ \frac{m^2}{(l+2)}+(l-1)Q_l^2-(l+2)Q_{l+1}^2+(l+1)Q_{l+2}^2 \right]r\bar{U}^{\prime }_{l+1} + \bigg [ \frac{4m^2}{(l+2)} +2(l+1)\big [ l+2-lQ_{l+2}^2 \big ]\nonumber \\ +\,2(l+3)[(l-1)Q_l^2 - (l+2)Q_{l+1}^2] + (l+1)\big [ (l+2)(1-Q_{l+1}^2-Q_{l+2}^2) + (l-1)Q_l^2 \big ]\frac{r\rho ^{\prime }_0}{\rho _0} \bigg ]\bar{U}_{l+1} \bigg \rbrace \nonumber \\ -\, (l+2)Q_{l+1}Q_{l+2}Q_{l+3}\left[ \frac{2}{(l+3)}r\bar{U}^{\prime }_{l+3}+\left({\frac{2(l+5)}{(l+3)}+\frac{r\rho ^{\prime }_0}{\rho _0}}\right)\bar{U}_{l+3} \right] + mQ_{l-1}Q_l[ (l-2)(l+1)-2(l-1) ]X_{l-2}\nonumber \\ +\, m\left\lbrace -l(l+1) + [ (l-2)(l+1)-2(l-1) ]Q_l^2 + [ l(l+3)+2(l+2) ]Q_{l+1}^2 \right\rbrace X_l + mQ_{l+1}Q_{l+2}[ l(l+3)+2(l+2) ]X_{l+2}, \end{eqnarray}

(165)

where we have defined

\begin{equation} X_l\equiv \frac{W_l}{r}, \bar{U}_l\equiv iU_l. \end{equation}

(166)

The former definition gives a toroidal function with the same dimension as the poloidal functions U_l; the latter produces a set of DEs where all coefficients are real (note that the |$\tilde{\Upsilon }_l$| functions are all imaginary, so |$i\tilde{\Upsilon }_l$| is real).

5.3.1 Analysis of equations

In principle, both equations (164) and (165) must be solved for all l > |m|, but upon closer inspection, we see that half of the equations are trivial. To see this, note that only the l = 2 and 4 terms of |$\tilde{\Psi }_l^m$| and the l = 1 and 3 terms of |$\tilde{\Upsilon }_l^m$| are non-zero (⁠|$\tilde{\Upsilon }^{\pm 2}_1$| is also zero, as expected since l < |m|). In addition, the equations only couple U_l to U_l + 2, U_l + 4, …, with the same being true for X_l. The result is that (164) for l = 2, 4, 6, … couples to (165) for l = 1, 3, 5, …, and to all the source terms. These are clearly the relevant equations to solve for our problem. On the other hand, equation (164) for odd l couples to (165) for even l, but since there are no terms sourcing the variables involved in these equations, we may take

\begin{equation} X_{{\rm even }\,\,l} = U_{{\rm odd }\,\,l} = 0. \end{equation}

(167)

Next, we argue that the terms with negative m are equal to plus-or-minus the corresponding terms with positive m (as usual for problems involving spherical-harmonic decompositions).

Although equations (164) and (165) are rather lengthy, it is enough to note that they take the schematic forms:

\begin{eqnarray} 0 = {\rm terms\, of\, the\, form }\,\lbrace m\tilde{\Psi }_l^m,\tilde{\Upsilon }_l^m,U_l^m,mX_l^m\rbrace , \end{eqnarray}

(168)

\begin{eqnarray} 0 = {\rm terms\, of\, the\, form}\,\lbrace \tilde{\Upsilon }_l^m,m^2\tilde{\Upsilon }_l^m,U_l^m,mX_l^m\rbrace , \end{eqnarray}

(169)

respectively (where we have restored the suppressed m superscripts on the variables). From the forms of the source terms, (A4) and (A5), we see that

\begin{equation} \tilde{\Upsilon }_l^1 = \tilde{\Upsilon }_l^{-1},\tilde{\Upsilon }_l^2 = -\tilde{\Upsilon }_l^{-2}, \tilde{\Psi }_l^1 = -\tilde{\Psi }_l^{-1},\tilde{\Psi }_l^2 = \tilde{\Psi }_l^{-2}, \end{equation}

(170)

and therefore, we expect

\begin{equation} U_l^1 = U_l^{-1},U_l^2 = -U_l^{-2}, X_l^1 = -X_l^{-1},X_l^2 = X_l^{-2}. \end{equation}

(171)

We have now accumulated a number of results – specifically those of equations (139), (166), (167) and (171) – which allow us to simplify the original expression (136) for |$\delta {\boldsymbol B}$| and rewrite it in terms of the new variables used in our system of DEs (164) and (165):

\begin{eqnarray} \delta {\boldsymbol B}&=& \sum \limits _{l=1}^\infty \bigg \lbrace-i\bar{U}_{2l}^1(Y_{2l}^1+Y_{2l}^{-1}){\boldsymbol e}_r - i\frac{(r^2\bar{U}_{2l}^1{}^{\prime }+2r\bar{U}_{2l}^1)}{l(l+1)}\nabla (Y_{2l}^1+Y_{2l}^{-1}) + rX_{2l-1}^1{\boldsymbol e}_r\times \nabla (Y_{2l-1}^1-Y_{2l-1}^{-1}) \nonumber \\ \, -\,i\bar{U}_{2l}^2(Y_{2l}^2-Y_{2l}^{-2}){\boldsymbol e}_r - i\frac{(r^2\bar{U}_{2l}^2{}^{\prime }+2r\bar{U}_{2l}^2)}{l(l+1)}\nabla (Y_{2l}^2-Y_{2l}^{-2}) + rX_{2l+1}^2{\boldsymbol e}_r\times \nabla (Y_{2l+1}^2+Y_{2l+1}^{-2}) \bigg \rbrace . \end{eqnarray}

(172)

Since the quantities |$\bar{U}_l^m,X_l^m,(Y_l^1-Y_l^{-1}),(Y_l^2+Y_l^{-2})$| are all real, whilst |$(Y_l^1+Y_l^{-1})$| and |$(Y_l^2-Y_l^{-2})$| are imaginary, it is clear that |$\delta {\boldsymbol B}$| is also real, as expected.

5.4 Differential algebraic equations; conditions at the centre

The first of the above pair of equations, (164), is second order in |$\bar{U}_l$|⁠, but the second (165) is only first order. If we were to define a new variable

\begin{equation} Z_l\equiv \bar{U}^{\prime }_l \end{equation}

(173)

and substitute this back into the pair of equations above, the first equation would reduce to first order, and the second would be algebraic, with no derivatives. Thus, instead of a conventional system of ODEs, we have a system of DAEs. This is an unwelcome result, as DAEs are, in many senses, more difficult to solve than ODEs. Extra care is needed in choosing suitable boundary conditions for the differentiated variables, to ensure that they are consistent – i.e. that they satisfy any algebraic relations stemming from the DAE system. Furthermore, although initial-value problems for DAEs have been relatively well explored, we have a boundary-value problem, since we will have conditions to impose both at the centre and outer boundary. See e.g. Petzold (1982) for more discussion of the difficulties of solving DAEs.

Despite these problems, we have been able to find a solution method, albeit one with limited applicability. For this method, we chose to work with equation (164) in its original form, but coupled to the r-derivative of equation (165):

\begin{eqnarray} 0 &=& \frac{4\pi i}{\Lambda }\left\lbrace -(l-1)(l-2)Q_{l-1}Q_l\tilde{\Upsilon }^{\prime }_{l-2} + [ m^2+(l-1)(l+1)Q_l^2+l(l+2)Q_{l+1}^2 ]\tilde{\Upsilon }^{\prime }_l - (l+2)(l+3)Q_{l+1}Q_{l+2}\tilde{\Upsilon }^{\prime }_{l+2} \right\rbrace \nonumber \\ -\, (l-1)Q_{l-2}Q_{l-1}Q_l\left[ \frac{2}{(l-2)}r\bar{U}^{\prime \prime }_{l-3} -\left({\frac{2(l-5)}{(l-2)}+\frac{r\rho ^{\prime }_0}{\rho _0}}\right)\bar{U}^{\prime }_{l-3}-\left({\frac{r\rho ^{\prime }_0}{\rho _0}}\right)^{\prime } \bar{U}_{l-3} \right] + \frac{Q_l}{l}\bigg \lbrace 2\bigg [ \frac{m^2}{(l-1)}+lQ_{l-1}^2-(l-1)Q_l^2\nonumber \\ +\,(l+2)Q_{l+1}^2 \bigg ]r\bar{U}^{\prime \prime }_{l-1} + \bigg [ \frac{6m^2}{(l-1)} -2l(l-1) +2l(l+2)Q_{l-1}^2 + 2(l-3)(l-1)Q_l^2 - 2(l-3)(l+2)Q_{l+1}^2 \nonumber \\ -\, l\big [ (l-1)(1-Q_{l-1}^2-Q_l^2) + (l+2)Q_{l+1}^2 \big ]\frac{r\rho ^{\prime }_0}{\rho _0} \bigg ]\bar{U}^{\prime }_{l-1} - l\big [ (l-1)(1-Q_{l-1}^2-Q_l^2) + (l+2)Q_{l+1}^2 \big ]\left({\frac{r\rho ^{\prime }_0}{\rho _0}}\right)^{\prime }\bar{U}_{l-1} \bigg \rbrace \nonumber \\ +\, \frac{Q_{l+1}}{(l+1)}\bigg \lbrace 2\left[ \frac{m^2}{(l+2)}+(l-1)Q_l^2-(l+2)Q_{l+1}^2+(l+1)Q_{l+2}^2 \right]r\bar{U}^{\prime \prime }_{l+1} + \bigg [ \frac{6m^2}{(l+2)} +2(l+1)(l+2)+2(l-1)(l+4)Q_l^2\nonumber \\ -\,2(l+2)(l+4)Q_{l+1}^2-(l-1)(l+1)Q_{l+2}^2 + (l+1)\big [ (l+2)(1-Q_{l+1}^2-Q_{l+2}^2) + (l-1)Q_l^2 \big ]\frac{r\rho ^{\prime }_0}{\rho _0} \bigg ]\bar{U}^{\prime }_{l+1} \nonumber \\ +\, (l+1)\big [ (l+2)(1-Q_{l+1}^2-Q_{l+2}^2) + (l-1)Q_l^2 \big ]\left({\frac{r\rho ^{\prime }_0}{\rho _0}}\right)^{\prime }\bar{U}_{l+1} \bigg \rbrace - (l+2)Q_{l+1}Q_{l+2}Q_{l+3}\bigg [ \frac{2}{(l+3)}r\bar{U}^{\prime \prime }_{l+3}\nonumber \\ +\,\left({\frac{2(l+6)}{(l+3)}+\frac{r\rho ^{\prime }_0}{\rho _0}}\right)\bar{U}^{\prime }_{l+3}+\left({\frac{r\rho ^{\prime }_0}{\rho _0}}\right)^{\prime }\bar{U}_{l+3} \bigg ] + mQ_{l-1}Q_l[ (l-2)(l+1)-2(l-1) ]X^{\prime }_{l-2}\nonumber \\ +\, m\left\lbrace -l(l+1) + [ (l-2)(l+1)-2(l-1) ]Q_l^2 + [ l(l+3)+2(l+2) ]Q_{l+1}^2 \right\rbrace X^{\prime }_l + mQ_{l+1}Q_{l+2}[ l(l+3)+2(l+2) ]X^{\prime }_{l+2}. \end{eqnarray}

(174)

Now we really do have a conventional system of coupled ODEs, where both equations (per l, m) are second order in U_l and first order in X_l – and may be solved with conventional numerical methods. However, the set of solutions to this new system of equations is larger than those that solve our original problem. In particular, if we integrate equation (174), we recover the original equation plus some arbitrary integration constant. Plugging a solution to the differentiated equation back into the original system thus generically results in an substantial r-independent error; see Appendix B.

In order to fix the integration constant, we need information from the original equation (165). We note that since all derivatives are pre-multiplied by r, if we evaluate (165) at the centre of the star, these are zero, assuming sufficient regularity of the |$\bar{U}$| and X functions. In addition, the source terms |$\tilde{\Psi }_l^m,\tilde{\Upsilon }_l^m$| are zero. We are thus left with a set of algebraic equations in U_l(0) and X_l(0) (one per value of l):

\begin{eqnarray} 0 &=& \frac{(l-4)(l-1)}{(l-2)}Q_{l-2}Q_{l-1}Q_l\bar{U}_{l-3}(0) + \frac{Q_l}{l}\left[ \frac{2m^2}{(l-1)}-l(l-1)+l(l+1)Q_{l-1}^2+(l-2)(l-1)Q_l^2 -(l-2)(l+2)Q_{l+1}^2\right]\bar{U}_{l-1}(0) \nonumber \\ +\, \frac{Q_{l+1}}{(l+1)}\left[ \frac{2m^2}{(l+2)}+(l+1)(l+2)+(l-1)(l+3)Q_l^2 -(l+2)(l+3)Q_{l+1}^2 -l(l+1)Q_{l+2}^2\right]\bar{U}_{l+1}(0)\nonumber \\ -\, \frac{(l+2)(l+5)}{(l+3)}Q_{l+1}Q_{l+2}Q_{l+3}\bar{U}_{l+3}(0) + mQ_{l-1}Q_l\left[\frac{(l-2)(l+1)}{2}-l+1\right]X_{l-2}(0)\nonumber \\ +\, m\left\lbrace -\frac{l(l+1)}{2}+\left[ \frac{(l-2)(l+1)}{2}-l+1 \right]Q_l^2 +\left[ \frac{l(l+3)}{2}+l+2 \right]Q_{l+1}^2 \right\rbrace X_l(0) + mQ_{l+1}Q_{l+2}\left[\frac{l(l+3)}{2}+l+2\right]X_{l+2}(0). \end{eqnarray}

(175)

In numerical solutions of systems of DEs like (164) and (174), one truncates the infinite set of equations at some finite value of the angular index l_max. For a particular l_max, equation (175) may be ‘solved’ – it does not give numerical values for all the variables at the centre, of course, since this requires information from the full equations including the source terms, but it does allow one to express these variables in terms of each other. For example, for l_max = 4 and m = 1, one can show that equation (175) implies that

\begin{eqnarray} \bar{U}_4^1(0) = -\frac{55}{6\sqrt{6}}\bar{U}_2^1(0), \end{eqnarray}

(176)

\begin{eqnarray} X_3^1(0) = -\frac{7}{12}\sqrt{\frac{35}{2}}\bar{U}_2^1(0). \end{eqnarray}

(177)

Note that for this example, and in fact for all l_max, X₁(0) drops out of equation (175) completely, as its pre-factor is always zero.

Unfortunately, applying the same method to higher values of l_max does not give solutions to the original system of DAEs. We believe this is because the relations obtained at the centre from (175) become more complex – relating, for example, all quantities to a combination of |$\bar{U}_2(0)$|and|$\bar{U}_4(0)$|⁠, instead of |$\bar{U}_2(0)$| alone. As a result, it seems that the centre conditions do not give sufficient information to fix the integration constant for higher l_max.

It would naturally have been more satisfactory to go to higher values of l_max to check the convergence of the solution in the limit l → ∞. However, the low-l, m nature of the source terms gives us reason to believe that the full l → ∞ sum will be dominated by lower multipoles, and so we anticipate that our results will be representative of the full, untruncated, solution. We will later find some hints of how the higher l components behave, in Section 7.2.

5.5 Exterior solution

To understand what boundary conditions are appropriate, let us first recall the physical meaning of the quantities in our equations. From (136) and (139), we see that the U_l functions are associated with the poloidal component of |$\delta {\boldsymbol B}$|⁠, and the W_l functions with the toroidal component. We may, therefore, regard the two unknown quantities in the above coupled DEs as being |$\delta {\boldsymbol B}_{{\rm pol}}$| and |$\delta {\boldsymbol B}_{{\rm tor}}$|⁠. Note that by eliminating the Ξ_l functions, we have produced a set of equations with no dependence on second-order fluid perturbations like δH_αB; see the discussion after equation (134). This means that we only need to impose boundary conditions on the perturbed magnetic field, rather than on the fluid too.

We assume that the exterior of the star is vacuum; it may be more realistic to model a charge-carrying magnetosphere, but the aim of this study is to isolate the physics of a quasi-rigidly precessing magnetized fluid star. Assuming a vacuum exterior means there are no particles to carry an electric current |$\delta {\boldsymbol J}$|⁠. Using Ampère's law, this gives us

\begin{equation} \nabla \times \delta {\boldsymbol B}= 4\pi \delta {\boldsymbol J} = 0\ \ {\rm (}r>R_*{\rm )}. \end{equation}

(178)

From (140), we require that each of the three components of the vector-spherical-harmonic decomposition of |$\nabla \times \delta {\boldsymbol B}$| be zero outside the star, i.e.

\begin{equation} -\frac{l(l+1)W_l^m}{r^2} = -W_l^m{}^{\prime } = V_l^m{}^{\prime } - U_l^m = 0\ \ {\rm (}r>R_*{\rm )}. \end{equation}

(179)

As for the equations for the interior, the exterior equations decouple for different azimuthal index, and so we may suppress the m indices. Using the above result, and eliminating |$V_l^m$| in favour of |$U_l^m$| using (139), the equations for the exterior are

\begin{eqnarray} 0 = r^2U^{\prime \prime }_l + 4rU^{\prime }_l+[2-l(l+1)]U_l \ \ {\rm (}r>R_*{\rm )}, \end{eqnarray}

(180)

\begin{eqnarray} 0 = W_l\ \ {\rm (}r>R_*{\rm )}. \end{eqnarray}

(181)

The latter equation immediately tells us that there is no exterior toroidal field. The former relation is a Cauchy–Euler equation, which we may solve by first making the ansatz that U_l = r^p; plugging this into equation (180) gives us an indicial equation, which can be factorized to give the values of p:

\begin{equation} 0 = p^2 + 3p + [2-l(l+1)]p = (p+l+2)(p-l+1) \Rightarrow p = -l-2,p = l-1. \end{equation}

(182)

Thus the general exterior solution |$U^{{\rm ext}}_l$| for U_l is

\begin{equation} U^{{\rm ext}}_l(r) = u^{\rm ext}_l r^{-l-2}+\hat{u}^{\rm ext}_l r^{l-1}. \end{equation}

(183)

The latter term in this expression diverges at infinite radius, however, and so is unphysical. Therefore we set |$\hat{u}^{\rm ext}_l = 0$|⁠, leaving us with

\begin{equation} U^{{\rm ext}}_l(r) = u^{\rm ext}_l r^{-l-2}. \end{equation}

(184)

Note that this gives us the expected physical result: the dipole component (l = 1) of a poloidal field should decay as 1/r³, the quadrupole as 1/r⁴, and so on. Differentiating the above, we obtain a condition on U_l that does not involve the unknown constants |$u_l^{\rm ext}$|⁠:

\begin{equation} (U^{\rm ext}_l)^{\prime } = -(l+2)u_l^{\rm ext}r^{-(l+3)} = -\frac{(l+2)}{r}U^{\rm ext}_l \Rightarrow rU^{\prime }_l + (l+2)U_l = 0\ \ {\rm (}r>R_*{\rm )}. \end{equation}

(185)

5.6 Internal solution behaviour as outer boundary → R_*

Within our model, the most natural place to impose boundary conditions matching the interior and exterior solutions would be at the stellar surface R_*. A closer look at equations (164) and (165), however, indicates problems with doing so. In particular, the equations feature terms of the form |$\rho ^{\prime }_0/\rho _0$| and the radial derivative thereof, both of which diverge as r → R_* for a γ = 2 polytrope.⁶

The divergent terms originate from the perturbed Lorentz force; see equation (51). Of these three terms, two have a denominator of ρ₀ and one has a denominator of |$\rho _0^2$|⁠; we therefore need corresponding factors in the numerators to avoid the Lorentz force diverging when ρ₀ → 0. The second of the three Lorentz-force terms from (51) involves a cross-product with |${\boldsymbol B}_0$|⁠, which, in turn, scales with ρ₀, and so this term is indeed well behaved. The first and third terms will diverge, however, unless |$\nabla \times {\boldsymbol B}_0\rightarrow 0$| at least as fast as ρ₀ – which it does not; it contains one term that instead scales with |$\rho ^{\prime }_0$| (this can be seen by direct evaluation of the curl of equation 58, or by applying relation (140) to |${\boldsymbol B}_0$| in the form given in equation 137), and |$\rho ^{\prime }_0$| is non-zero at R_*. As originally noted by Mestel et al. (1981), this divergence is a problem before even reaching the surface, since the perturbed Lorentz force – an |$\mathcal {O}(\epsilon _\alpha \epsilon _B)$| quantity – will eventually become numerically larger than the |$\mathcal {O}(\epsilon _\alpha )$| centrifugal force from equation (27) at some radius r < R_*. At this point, our perturbative ordering breaks down, and our scheme and equations become invalid.

To see how this problem would manifest itself in practice, we solved the system of equations (164) and (174) up to an outer radius R_out < R_*, and investigated the behaviour of the solutions in the limit R_out → R_*. As R_out was increased, the amplitude of the radial functions |$\bar{U},X$| was seen to diverge. At the inner and outer boundaries, however, the solutions were well behaved.

The simplest resolution to this problem is to place an outer boundary R_out somewhere beneath the stellar surface. For the stellar model studied in Mestel et al. (1981), this was chosen to be 0.99R_*. Our numerical results from varying the cut-off radius, on the other hand, showed that the solutions were relatively insensitive to the exact cut-off value provided R_out ≲ 0.93R_*; accordingly, we will place the outer boundary for the calculation at a value R_out = 0.9R_*. For main-sequence stars and white dwarfs, this is an arbitrary choice, and a discussion of its validity is given at the end of this paper, in Section 8. For neutron stars, however, we can give a more quantitative argument for placing the outer boundary at 0.9R_*: it is (approximately) the radius at which the star's fluid core gives way to a solid, elastic crust, and so provides a natural cut-off for our fluid calculation. Next, we discuss how our solutions might be affected by elastic forces in a neutron-star crust, if we were to extend our modelling to include this region.

5.7 Conditions in the crust of a neutron star

To extend our analysis to include the effects of a neutron-star crust, we would need to include elastic-force terms in our perturbation equations. This is not straightforward, however, as the resulting terms have a separate perturbative scaling from those in our problem. In addition, one has to define what the relaxed state of the crust is – i.e. the state in which the stresses are zero (and therefore the crust is described by fluid terms alone). In a real neutron star, the crust's elastic stress pattern is likely to be complicated and to depend on the star's rotational and seismic history.

For the purposes of this paper, we just wish to make a qualitative assessment of how elasticity might modify our solution; for this, there are two limiting cases. In one limit, the crust is so strong that it can maintain its spheroidal centrifugal distortion even as the star's instantaneous rotation vector is substantially misaligned from the primary rotation vector (recall Fig. 2); the crust would therefore undergo rigid-body free precession with no additional |$\xi$|-motions. In the other limit, the crust is so weak that elastic forces have negligible contribution at all perturbative orders we consider. We will argue next that a neutron-star crust falls between these two limits. A separate issue to consider is whether the stresses could be large enough to exceed the elastic yield limit of the crust. Finally, we discuss the coupling between the crust and core in a magnetized precessing neutron star.

5.7.1 Importance of elastic force terms

We begin by assessing the size of elastic force terms. For the sake of simplicity, let us define the crust as being relaxed in the presence of its magnetic distortion; this is not crucial for the argument, but makes it clearer. In this case, the combined background and |$\mathcal {O}(\epsilon _B)$| equations will look the same as in the fluid case, since no elastic terms will appear. If we then subtract force terms of these two orders from the full, unperturbed Euler equation (11), we are left with an equation containing |$\mathcal {O}(\epsilon _\alpha )$| terms, |$\mathcal {O}(\epsilon _\alpha \epsilon _B)$| terms and the elastic-force terms per unit mass |$\delta {\boldsymbol F}_{\rm el}$|⁠, whose leading-order pieces scale as

\begin{equation} \delta {\boldsymbol F}_{\rm el}\sim \frac{\mu }{\rho }\nabla ^2\xi \sim \frac{\mu }{\rho }\frac{\xi }{R_*^2}\sim \frac{\mu \epsilon _\alpha }{\rho R_*}, \end{equation}

(186)

where μ is the shear modulus. To find the perturbative order of this quantity, we non-dimensionalize it by dividing by |$P R^2/\mathcal {M}$|⁠:

\begin{equation} \frac{\delta {\boldsymbol F}_{\rm el}}{P R_*^2/\mathcal {M}} \sim \frac{\mu }{P\rho }\frac{\mathcal {M}}{R_*^3} \sim \frac{\mu }{P}\epsilon _\alpha . \end{equation}

(187)

Since this quantity scales with a new dimensionless parameter, μ/P, it prevents us from separating out the |$\mathcal {O}(\epsilon _\alpha )$| and |$\mathcal {O}(\epsilon _\alpha \epsilon _B)$| equations as before. The ratio μ/P ≲ 0.01 for neutron-star crusts – so one can think of the elastic forces as representing a numerically small correction to the |$\mathcal {O}(\epsilon _\alpha )$| equations, and therefore a small correction to our existing δρ_α solution, i.e. equation (95). On the other hand, they represent a potentially dominant piece of the |$\mathcal {O}(\epsilon _\alpha \epsilon _B)$| equations, and therefore could drastically change the solution for |$\delta {\boldsymbol B}$|⁠.

The conclusion we draw from this is that the crust is weak enough that the centrifugal bulge will always remain approximately symmetrical about the primary-rotation axis, but strong enough that our solutions for |$\dot{\xi }$| and |$\delta {\boldsymbol B}$| could be radically altered in the crustal region. Since we are obliged to place an outer boundary some distance within the star anyway, to ensure the validity of our perturbative scheme, we find it natural simply to exclude the crust from our modelling altogether, and place an outer boundary at the crust–core interface, r = 0.9R_*.

5.7.2 Exceeding the crustal breaking strain

Like any elastic medium, a neutron-star crust has a yield strain σ_max, beyond which value it ceases to respond elastically to additional strain and instead breaks.⁷ Recent molecular-dynamics simulations of neutron-star crustal matter (Horowitz & Kadau 2009) indicate that this happens at a very large value, σ_max ∼ 0.1. We need to compare this with the strain tensor induced by the |$\xi$|-motions, |$\nabla \xi$|⁠:

\begin{equation} \nabla \xi \sim \frac{\epsilon _\alpha R_*}{R_*} \sim \epsilon _\alpha \sim 0.2\left({\frac{\alpha }{{\rm 1 kHz}}}\right)^2, \end{equation}

(188)

where the right-hand side is an estimate of ε_α from Jones & Andersson (2002). Comparing this with σ_max, we conclude that the crust will not break unless the rotation frequency exceeds ∼700 Hz, about half the Keplerian frequency for a neutron star. Rotation rates more rapid than this value start to violate ε_α ≪ 1, one of our original ansätze, anyway – so this case is outside the confines of our model. We may therefore assume that the crust never breaks, without an additional loss of generality.

5.7.3 Crust–core coupling

This paper focusses on fluid stars whose only ‘rigidity’ comes from their magnetic fields. By contrast, the solid crust of a neutron star can, of course, undergo conventional rigid-body free precession without the aid of any magnetic field – provided it has some permanent distortion misaligned from the rotation axis. If we temporarily ignore the magnetic rigidity of the core, this fluid region will be unable to precess, and the resulting dynamics of the neutron star will be a precessing crust atop a rigidly rotating core.

In fact, the star's magnetic field would generally thread the crust and core, and couple together the two regions on an Alfvén time-scale τ_A. In the case where |$\tau _{\rm A} \lesssim \sqrt{T_\alpha T_\omega }$|⁠, Levin & D'Angelo (2004) showed that the core's motion is strongly coupled to that of the crust, so that the whole star would precess as one.⁸ Thus, in general, a neutron star's precession should combine two effects: the non-rigid precession on which this paper is focused (related to the magnetic distortion of the fluid) and a coupling to the rigid-body precession of the crust (related only to the crust's permanent distortion). We close by noting that for our particular stellar model, the latter effect may not actually occur: purely toroidal magnetic field lines are always non-radial and so no field line will cross from the crust to the core.

5.8 Boundary conditions for our model

Our system of equations, (164) and (174), is second-order in |$\bar{U}_l$| and first-order in X_l. This means that for each |$\bar{U}_l$|⁠, we need to impose two boundary conditions, and for each X_l, we need to impose one boundary condition.

We begin with boundary conditions at the centre, where we are obliged to use certain relations to fix the ‘integration constant’ in our problem, and thus to ensure our results are solutions of the original problem as encoded in equations (164) and (165) (see Section 5.4). For l_max = 4 and m = 1, the centre conditions (176) and (177) fix |$\bar{U}^1_4$| and |$X^1_3$| in terms of |$\bar{U}_2^1$|⁠, but the value of |$\bar{U}_2^1$| itself at the centre is unknown. Let us therefore simply demand that |$\bar{U}_2^1$| tends smoothly to some constant at the centre, by imposing the condition

\begin{equation} \frac{\mathrm{d}{\bar{U}_2^1}}{\mathrm{d}{r}}(r = 0) = 0. \end{equation}

(189)

So far |$X_1^1$| is unfixed; it cannot be written in terms of the other quantities at the centre, and its actual value there is unknown. We also cannot evade the problem as done above for |$\bar{U}_2^1$|⁠, by imposing |$X_1^1{}^{\prime }(0) = 0$|⁠, since our system of equations is only first order in X. Instead, we will impose an outer boundary condition on |$X^1_1$|⁠.

For l_max = 4 and m = 2, the only X-function is |$X_3^2$|⁠. The method of Section 5.4 gives |$X_3^2(0)$| in terms of the other two radial functions at the centre, |$U_2^2(0)$| and |$U_4^2(0)$|⁠:

\begin{equation} X_3^2(0) = -\frac{2}{15\sqrt{7}}U^2_2(0)+\frac{207}{550}\sqrt{\frac{3}{7}}U^2_4(0). \end{equation}

(190)

Again, we ensure that the |$\bar{U}$| functions of this expression approach some constant value at the centre by imposing

\begin{equation} \frac{\mathrm{d}{U_2^2}}{\mathrm{d}{r}}(r = 0) = \frac{\mathrm{d}{U_4^2}}{\mathrm{d}{r}}(r = 0) = 0. \end{equation}

(191)

We now move to the conditions at the outer boundary, corresponding to the crust-core boundary for a neutron star; here, we would ideally like to ensure that the interior perturbed magnetic field matches smoothly to its exterior counterpart. Mathematically speaking, this matching is not necessary, as steps in the non-radial field components across the boundary may be matched consistently using a current sheet (steps in the radial component would, however, violate |$\nabla \cdot {\boldsymbol B}= 0$|⁠). Physically, however, such a current sheet is undesirable, as it must be balanced by additional forces beyond those included in our problem, or otherwise would result in a net acceleration acting tangential to the boundary and oscillating on the precession time-scale.

For the continuity of the poloidal component, we need to match the interior U_l and V_l functions smoothly to the exterior, but by equation (139) |$V_l\sim U^{\prime }_l$|⁠, meaning that we need the continuity of U_land its radial derivative across the boundary R_out. This means we can safely differentiate the expression (184) at the boundary, as we were able to do outside the star, and extend the result (185) to the outer boundary:

\begin{equation} R_{\rm out} U^{\prime }_l(R_{\rm out})+(l+2)U_l(R_{\rm out}) = 0. \end{equation}

(192)

We have two boundary conditions to specify for each U_l, and since we have only specified one at the centre for each of them, we may indeed impose equation (192), ensuring the continuity of the poloidal field.

For the continuity of the toroidal component, we see from (136) that it suffices to match the interior and exterior W_l functions (or equivalently, the X_l functions) at r = R_out = 0.9R_*. But the exterior condition is X_l = 0, so these functions must also be zero at the outer boundary:

\begin{equation} X_l(R_{\rm out}) = 0. \end{equation}

(193)

For the X_l functions, we only have the freedom to impose a single boundary condition, however, and for |$X_3^1$| and |$X_3^2$|⁠, this must be the one at the centre, so that we pick out the correct DAE solution. Only for |$X_1^1$|⁠, which does not enter the centre conditions, are we able to enforce (193) at the outer boundary. In our results, we will see that the l = 3 X functions do indeed not vanish at the outer boundary, although their values there are numerically small in comparison with the typical value of the X functions in the interior. If one insisted on smooth interior–exterior matching for the toroidal field, it should be possible to solve the equations with the additional boundary condition (193) for allX, but at the expense of making the system into an eigenvalue problem.

The boundary conditions used in our solution of the ODE system, (164) and (174), are summarized in Table 2.

Table 2.

Open in new tab

Summary of boundary conditions for our problem, truncated at l = l_max = 4. There is one fewer boundary condition for m = 2, since X₁ does not exist for this case.

m	Centre	Crust–core boundary
	\|$U^{\prime }_2(0) = 0$\|	\|$R_{\rm out}U^{\prime }_2(R_{\rm out})+4U_2(R_{\rm out}) = 0$\|
1	\|$U_4(0) = -\frac{55}{6\sqrt{6}}U_2(0)$\|	\|$R_{\rm out}U^{\prime }_4(R_{\rm out})+6U_4(R_{\rm out}) = 0$\|
	\|$X_3(0) = -\frac{7}{12}\sqrt{\frac{35}{2}}U_2(0)$\|	X₁(R_out) = 0
	\|$U^{\prime }_2(0) = 0$\|	\|$R_{\rm out}U^{\prime }_2(R_{\rm out})+4U_2(R_{\rm out}) = 0$\|
2	\|$U^{\prime }_4(0) = 0$\|	\|$R_{\rm out}U^{\prime }_4(R_{\rm out})+6U_4(R_{\rm out}) = 0$\|
	\|$X_3(0) = -\frac{2}{15\sqrt{7}}U_2(0)+\frac{207}{550}\sqrt{\frac{3}{7}}U_4(0)$\|	–

m	Centre	Crust–core boundary
	\|$U^{\prime }_2(0) = 0$\|	\|$R_{\rm out}U^{\prime }_2(R_{\rm out})+4U_2(R_{\rm out}) = 0$\|
1	\|$U_4(0) = -\frac{55}{6\sqrt{6}}U_2(0)$\|	\|$R_{\rm out}U^{\prime }_4(R_{\rm out})+6U_4(R_{\rm out}) = 0$\|
	\|$X_3(0) = -\frac{7}{12}\sqrt{\frac{35}{2}}U_2(0)$\|	X₁(R_out) = 0
	\|$U^{\prime }_2(0) = 0$\|	\|$R_{\rm out}U^{\prime }_2(R_{\rm out})+4U_2(R_{\rm out}) = 0$\|
2	\|$U^{\prime }_4(0) = 0$\|	\|$R_{\rm out}U^{\prime }_4(R_{\rm out})+6U_4(R_{\rm out}) = 0$\|
	\|$X_3(0) = -\frac{2}{15\sqrt{7}}U_2(0)+\frac{207}{550}\sqrt{\frac{3}{7}}U_4(0)$\|	–

Table 2.

Open in new tab

Summary of boundary conditions for our problem, truncated at l = l_max = 4. There is one fewer boundary condition for m = 2, since X₁ does not exist for this case.

m	Centre	Crust–core boundary
	\|$U^{\prime }_2(0) = 0$\|	\|$R_{\rm out}U^{\prime }_2(R_{\rm out})+4U_2(R_{\rm out}) = 0$\|
1	\|$U_4(0) = -\frac{55}{6\sqrt{6}}U_2(0)$\|	\|$R_{\rm out}U^{\prime }_4(R_{\rm out})+6U_4(R_{\rm out}) = 0$\|
	\|$X_3(0) = -\frac{7}{12}\sqrt{\frac{35}{2}}U_2(0)$\|	X₁(R_out) = 0
	\|$U^{\prime }_2(0) = 0$\|	\|$R_{\rm out}U^{\prime }_2(R_{\rm out})+4U_2(R_{\rm out}) = 0$\|
2	\|$U^{\prime }_4(0) = 0$\|	\|$R_{\rm out}U^{\prime }_4(R_{\rm out})+6U_4(R_{\rm out}) = 0$\|
	\|$X_3(0) = -\frac{2}{15\sqrt{7}}U_2(0)+\frac{207}{550}\sqrt{\frac{3}{7}}U_4(0)$\|	–

m	Centre	Crust–core boundary
	\|$U^{\prime }_2(0) = 0$\|	\|$R_{\rm out}U^{\prime }_2(R_{\rm out})+4U_2(R_{\rm out}) = 0$\|
1	\|$U_4(0) = -\frac{55}{6\sqrt{6}}U_2(0)$\|	\|$R_{\rm out}U^{\prime }_4(R_{\rm out})+6U_4(R_{\rm out}) = 0$\|
	\|$X_3(0) = -\frac{7}{12}\sqrt{\frac{35}{2}}U_2(0)$\|	X₁(R_out) = 0
	\|$U^{\prime }_2(0) = 0$\|	\|$R_{\rm out}U^{\prime }_2(R_{\rm out})+4U_2(R_{\rm out}) = 0$\|
2	\|$U^{\prime }_4(0) = 0$\|	\|$R_{\rm out}U^{\prime }_4(R_{\rm out})+6U_4(R_{\rm out}) = 0$\|
	\|$X_3(0) = -\frac{2}{15\sqrt{7}}U_2(0)+\frac{207}{550}\sqrt{\frac{3}{7}}U_4(0)$\|	–

6 THE |$\xi$|-MOTIONS

We have derived a set of DEs whose solution gives us |$\delta {\boldsymbol B}$|⁠. Together with the simple solution for δρ_α (see equation 95), this gives us enough information to solve for the perturbed velocity field |$\dot{\xi }$|⁠, using the perturbed continuity and induction equations (104) and (105). In this section, we show that |$\dot{\xi }$| may be found explicitly in terms of the magnetic functions |$U_l^m$| and |$W_l^m$|⁠, through relatively simple algebraic manipulation of equations (104) and (105).

The first task is to unwrap the curl operator in the induction equation, for which we use the fact that the magnetic field must be expressible as a Mie representation (106):

\begin{equation} \nabla \times (\dot{\xi }\times {\boldsymbol B}_0) = \delta \dot{{\boldsymbol B}}= \nabla \times (\nabla \dot{\mathcal {P}}\times {\boldsymbol r}+\dot{\mathcal {T}}{\boldsymbol r}) \Rightarrow \dot{\xi }\times {\boldsymbol B}_0 = \nabla \dot{\mathcal {P}}\times {\boldsymbol r}+\dot{\mathcal {T}}{\boldsymbol r}+\nabla \dot{\mathcal {G}}, \end{equation}

(194)

where |$\mathcal {P},\mathcal {T}$| are poloidal and toroidal scalar functions, respectively, and |$\mathcal {G}$| is a third scalar function representing the freedom that – upon equating terms under a curl operator – the resulting equation is only defined up to the addition of curl-free terms |$\nabla \mathcal {G}$|⁠.

We first need, therefore, to rewrite the general expression for |$\delta {\boldsymbol B}$| from (136) in its Mie representation. Note that we could, in fact, commence with |$\delta {\boldsymbol B}$| in the simpler form given in (172), which is specific to our problem. We prefer, however, to consider the general case, so that parts of the analysis here may be directly applied to other problems (for example, precession of a magnetic star where |${\boldsymbol B}_0$| is a poloidal or mixed poloidal-toroidal field).

To find the toroidal scalar |$\mathcal {T}$|⁠, decomposed in spherical harmonics as |$\mathcal {T} = \sum _{l,m}\mathcal {T}_l^m Y_l^m$|⁠, we need only rearrange the toroidal component |$\delta {\boldsymbol B}_{\rm tor}$| from equation (136):

\begin{equation} \delta {\boldsymbol B}_{\rm tor} = \sum \limits _{l,m} W_l^m{\boldsymbol e}_r\times \nabla Y_l^m = -\sum \limits _{l,m} \nabla \times \left({\frac{W_l^m Y_l^m}{r}{\boldsymbol r}}\right) = \sum \limits _{l,m} \nabla \times \left({\mathcal {T}_l^m Y_l^m {\boldsymbol r}}\right), \end{equation}

(195)

using the product rule. In the above we have, as before, defined

\begin{equation} \sum \limits _{l,m}\equiv \sum \limits _{l = 1}^\infty \sum \limits _{m = -l,m\ne 0}^{l} \end{equation}

(196)

for brevity.

For the poloidal component |$\delta {\boldsymbol B}_{\rm pol}$| we begin by writing |$\mathcal {P} = \sum _{l,m} \mathcal {P}_l^m Y_l^m$| and look at one arbitrary (l, m)-component in the sum for |$\nabla \times (\nabla \mathcal {P}\times {\boldsymbol r})$|⁠:

\begin{eqnarray} \nabla \times [\nabla (\mathcal {P}_l^m Y_l^m)\times {\boldsymbol r}] &=& \nabla \mathcal {P}_l^m\times (\nabla Y_l^m\times {\boldsymbol r}) + \mathcal {P}_l^m\nabla \times (r\nabla Y_l^m\times {\boldsymbol e}_r)\nonumber \\ &=& r\mathcal {P}_l^m{}^{\prime }\nabla Y_l^m + \mathcal {P}_l^m\nabla Y_l^m + r\mathcal {P}_l^m\nabla \times (\nabla Y_l^m\times {\boldsymbol e}_r)\nonumber \\ &=& (r\mathcal {P}_l^m{}^{\prime }+\mathcal {P}_l^m)\nabla Y_l^m + \frac{l(l+1)\mathcal {P}_l^m}{r}Y_l^m{\boldsymbol e}_r, \end{eqnarray}

(197)

using similar algebraic steps to those in equation (140). Now comparing the above result with the expression obtained by substituting (139) into (136):

\begin{equation} \delta {\boldsymbol B}_{\rm pol} = \sum \limits _{l,m} \left[ \frac{(r^2 U_l^m{}^{\prime }+2rU_l^m)}{l(l+1)}\nabla Y_l^m + U_l^m Y_l^m{\boldsymbol e}_r \right], \end{equation}

(198)

we see that the poloidal radial functions in the Mie representation are given by

\begin{equation} \mathcal {P}_l^m = \frac{rU_l^m}{l(l+1)}. \end{equation}

(199)

We are now able to write |$\dot{\xi }\times {\boldsymbol B}_0$| in vector-spherical-harmonic form:

\begin{eqnarray} \dot{\xi }\times {\boldsymbol B}_0 &=& \frac{\mathrm{\partial} { }}{\mathrm{\partial} {t}}\left({\nabla \mathcal {P}\times {\boldsymbol r}+\mathcal {T}{\boldsymbol r}+\nabla \mathcal {G}}\right)\nonumber \\ & =& \frac{\mathrm{\partial} { }}{\mathrm{\partial} {t}}\sum \limits _{l,m}\left[r\mathcal {P}_l^m\nabla Y_l^m\times {\boldsymbol e}_r +r\mathcal {T}_l^m Y_l^m{\boldsymbol e}_r +\nabla (\mathcal {G}_l^m Y_l^m) \right]\nonumber \\ &=& \sum \limits _{l,m}\left[ (\mathcal {G}_l^m{}^{\prime }-W_l^m) \frac{\mathrm{\partial} {Y_l^m}}{\mathrm{\partial} {t}}{\boldsymbol e}_r +\mathcal {G}_l^m\nabla \left({\frac{\mathrm{\partial} {Y_l^m}}{\mathrm{\partial} {t}}}\right) + \frac{r^2 U_l^m}{l(l+1)}\nabla \left({\frac{\mathrm{\partial} {Y_l^m}}{\mathrm{\partial} {t}}}\right)\times {\boldsymbol e}_r \right] \nonumber \\ &=& \omega \sum \limits _{l,m}\left[im\left({ (\mathcal {G}_l^m{}^{\prime }-W_l^m) Y_l^m{\boldsymbol e}_r +\mathcal {G}_l^m\nabla Y_l^m + \frac{r^2 U_l^m}{l(l+1)}\nabla Y_l^m\times {\boldsymbol e}_r }\right) \right]\!, \end{eqnarray}

(200)

where we have used the result

\begin{equation} \frac{\mathrm{\partial} {Y_l^m(\theta ,\lambda )}}{\mathrm{\partial} {t}} = \frac{\mathrm{\partial} {Y_l^m}}{\mathrm{\partial} {\lambda }}\frac{\mathrm{\partial} {\lambda }}{\mathrm{\partial} {t}} = \frac{\mathrm{\partial} {Y_l^m}}{\mathrm{\partial} {\lambda }}\frac{\mathrm{\partial} {}}{\mathrm{\partial} {t}}(\omega t+\phi ) = im\omega Y_l^m \end{equation}

(201)

and have, predictably, defined |$\mathcal {G}_l^m$| as the radial terms in a spherical-harmonic expansion of |$\mathcal {G}$|⁠, i.e. |$\mathcal {G} = \sum _{l,m}\mathcal {G}_l^m Y_l^m$|⁠. Now, since we have a toroidal background field |${\boldsymbol B}_0 = B_0{\boldsymbol e}_\lambda$|⁠,

\begin{equation} \dot{\xi }\times {\boldsymbol B}_0 = B_0(\dot{\xi }_\theta {\boldsymbol e}_r-\dot{\xi }_r{\boldsymbol e}_\theta ). \end{equation}

(202)

We now use equations (200) and (202) to equate the r, θ and λ components of |$\dot{\xi }\times {\boldsymbol B}_0$|⁠:

\begin{eqnarray} {}[\dot{\xi }\times {\boldsymbol B}_0]_r =\ B_0\dot{\xi }_\theta \ = \omega \sum \limits _{l,m} [imY_l^m(\mathcal {G}_l^m{}^{\prime }-W_l^m)], \end{eqnarray}

(203)

\begin{eqnarray} {}[\dot{\xi }\times {\boldsymbol B}_0]_\theta = -B_0\dot{\xi }_r = \omega \sum \limits _{l,m} \left[ \frac{1}{r\sin \theta } \left({ -\frac{m^2 r^2 U_l^m Y_l^m}{l(l+1)} + im\mathcal {G}_l^m\sin \theta Y_{l,\theta }^m}\right) \right], \end{eqnarray}

(204)

\begin{eqnarray} [\dot{\xi }\times {\boldsymbol B}_0]_\lambda =\ \ \ 0\ \ \ \ = \omega \sum \limits _{l,m} \left[ -\frac{1}{r\sin \theta } \left({ \frac{imr^2U_l^m\sin \theta Y_{l,\theta }^m}{l(l+1)} + m^2\mathcal {G}_l^m Y_l^m}\right) \right]. \end{eqnarray}

(205)

First, by the orthogonality of the e^imλ terms within spherical harmonics, we know that equation (205) must hold for each m separately, so that this equation corresponds to a set of equalities involving sums over l alone. In each of these equalities it is permissible to take m outside the sum, and in particular to divide by m², so that

\begin{equation} \sum \limits _{l\ge |m|}\mathcal {G}_l^m Y_l^m = -\sum \limits _{l\ge |m|}\frac{ir^2 U_l^m\sin \theta Y_l^m{}_{,\theta }}{l(l+1)m} = -\sum \limits _{l\ge |m|}\frac{ir^2 U_l^m}{l(l+1)m} [lQ_{l+1}Y_{l+1}^m-(l+1)Q_l Y_{l-1}^m], \end{equation}

(206)

using equation (148). Differentiating this relation with respect to r gives |$(\mathcal {G}_l^m Y_l^m)_{,r} = \mathcal {G}_l^m{}^{\prime } Y_l^m$|⁠; we substitute this result back into equation (203) to eliminate |$\mathcal {G}_l^m$| and get an explicit expression for the θ-component of |$\dot{\xi }$|⁠:

\begin{equation} \dot{\xi }_\theta = -\frac{\omega }{B_0}\sum \limits _{l,m}\bigg \lbrace imW_l^m Y_l^m - \frac{(r^2 U_l^m{}^{\prime }+2rU_l^m)}{l(l+1)}[lQ_{l+1}Y_{l+1}^m-(l+1)Q_l Y_{l-1}^m] \bigg \rbrace . \end{equation}

(207)

Next, with a view to eliminate |$\mathcal {G}$| from equation (204), we use (206) to show that

\begin{equation} \sum \limits _{l\ge |m|}\mathcal {G}_l^m Y_{l,\theta }^m = \sum \limits _{l\ge |m|} (\mathcal {G}_l^m Y_l^m)_{,\theta } = -\sum \limits _{l\ge |m|}\frac{ir^2 U_l^m}{l(l+1)m}(\sin \theta Y_{l,\theta }^m)_{,\theta } = -\sum \limits _{l\ge |m|}\frac{ir^2 U_l^m}{l(l+1)m} [ lQ_{l+1}Y_{l+1}^m{}_{,\theta }-(l+1)Q_l Y_{l-1}^m{}_{,\theta } ]. \end{equation}

(208)

Now applying (148) to the spherical harmonic derivatives on the right-hand side of the above equation, we can show that

\begin{equation} \sum \limits _{l\ge |m|}\mathcal {G}_l^m\sin \theta Y_l^m{}_{,\theta } = -\sum \limits _{l\ge |m|}\frac{ir^2 U_l^m}{l(l+1)m}\Big \lbrace l(l+1)Q_{l+1}Q_{l+2}Y_{l+2}^m - [ l(l+2)Q_{l+1}^2+(l-1)(l+1)Q_l^2 ]Y_l^m + l(l+1)Q_{l-1}Q_l Y_{l-2}^m \Big \rbrace . \end{equation}

(209)

Substituting this result into (204) yields an explicit expression for the r-component of |$\dot{\xi }$|⁠:

\begin{equation} \dot{\xi }_r = -\frac{\omega }{B_0}\sum \limits _{l,m} \bigg \lbrace \frac{rU_l^m}{l(l+1)\sin \theta } \Big \lbrace l(l+1)Q_{l+1}Q_{l+2}Y_{l+2}^m - [m^2+l(l+2)Q_{l+1}^2+(l-1)(l+1)Q_l^2 ]Y_l^m + l(l+1)Q_{l-1}Q_l Y_{l-2}^m \Big \rbrace \bigg \rbrace . \end{equation}

(210)

The induction equation has given us expressions for |$\dot{\xi }_r$| and |$\dot{\xi }_\theta$|⁠. To find the azimuthal component |$\dot{\xi }_\lambda$|⁠, we now use the continuity equation:

\begin{equation} -\delta \dot{\rho }_\alpha = \nabla \cdot (\rho _0\dot{\xi }) = \rho ^{\prime }_0\dot{\xi }_r + \rho _0 \left[ \frac{(r^2\dot{\xi }_r)_{,r}}{r^2} + \frac{(\sin \theta \dot{\xi }_\theta )_{,\theta }}{r\sin \theta } + \frac{\dot{\xi }_{\lambda ,\lambda }}{r\sin \theta } \right], \end{equation}

(211)

which we rearrange to give

\begin{equation} \dot{\xi }_{\lambda ,\lambda } = -\frac{r\sin \theta \delta \dot{\rho }_\alpha }{\rho _0} - \frac{r\sin \theta \rho ^{\prime }_0\dot{\xi }_r}{\rho _0} -\frac{(r^2\dot{\xi }_r)_{,r}\sin \theta }{r} - (\sin \theta \dot{\xi }_\theta )_{,\theta }. \end{equation}

(212)

The middle two terms on the right-hand side of this equation may be written down immediately, whilst the other two need minor algebraic manipulation. First, differentiating (95) with respect to time yields:

\begin{equation} \delta \dot{\rho }_\alpha = \frac{5\pi \alpha ^2\omega }{8G} j_2\left({ {\frac{\pi r}{R_*}}}\right) \sin \theta [ \sin (2\chi )\cos \theta \sin \lambda + \sin ^2\chi \sin \theta \sin (2\lambda ) ]. \end{equation}

(213)

The final term on the right-hand side of (212) is

\begin{eqnarray} -(\sin \theta \dot{\xi }_\theta )_{,\theta } &=& \omega \frac{\mathrm{\partial} { }}{\mathrm{\partial} {\theta }}\left({\frac{\sin \theta }{B_0}}\right)\sum \limits _{l,m}\bigg \lbrace imW_l^m Y_l^m - \frac{(r^2 U_l^m{}^{\prime }+2rU_l^m)}{l(l+1)}[lQ_{l+1}Y_{l+1}^m-(l+1)Q_l Y_{l-1}^m] \bigg \rbrace \nonumber \\ &&+\, \frac{\omega }{B_0}\sum \limits _{l,m}\bigg \lbrace imW_l^m\sin \theta Y_l^m{}_{,\theta } - \frac{(r^2 U_l^m{}^{\prime }+2rU_l^m)}{l(l+1)}[lQ_{l+1}\sin \theta Y_{l+1}^m{}_{,\theta }-(l+1)Q_l\sin \theta Y_{l-1}^m{}_{,\theta }] \bigg \rbrace . \end{eqnarray}

(214)

However, the first right-hand-side term of this vanishes, since

\begin{equation} \frac{\mathrm{\partial} { }}{\mathrm{\partial} {\theta }}\left({\frac{\sin \theta }{B_0}}\right) = \frac{\mathrm{\partial} { }}{\mathrm{\partial} {\theta }}\left({\frac{1}{\Lambda r\rho _0}}\right) = 0 \end{equation}

(215)

for our background field, given by equation (58). The result, upon a further application of (148), is

\begin{eqnarray} -(\sin \theta \dot{\xi }_\theta )_{,\theta } &=& \frac{\omega }{B_0}\sum \limits _{l,m}\bigg \lbrace imW_l^m [ lQ_{l+1}Y_{l+1}^m-(l+1)Q_l Y_{l-1}^m ] - \frac{(r^2 U_l^m{}^{\prime }+2rU_l^m)}{l(l+1)}\Big \lbrace l(l+1)Q_{l+1}Q_{l+2}Y_{l+2}^m \nonumber \\ \, -\, [ l(l+2)Q_{l+1}^2+(l-1)(l+1)Q_l^2 ]Y_l^m + l(l+1)Q_{l-1}Q_l Y_{l-2}^m \Big \rbrace \bigg \rbrace . \end{eqnarray}

(216)

Now evaluating the right-hand side of (212), we find that a large number of the terms involving |$U_l^m$| functions cancel, and we are left with

\begin{eqnarray} \dot{\xi }_{\lambda ,\lambda } &=& -\frac{5\pi \alpha ^2\omega }{8G}\frac{r j_2({ {\frac{\pi r}{R_*}}})}{\rho _0} \sin ^2\theta [ \sin (2\chi )\cos \theta \sin \lambda + \sin ^2\chi \sin \theta \sin (2\lambda ) ]\nonumber \\ +\, \frac{\omega }{B_0}\sum \limits _{l,m} \Big \lbrace imW_l^m [lQ_{l+1}Y_{l+1}^m-(l+1)Q_l Y_{l-1}^m] -\frac{m^2}{l(l+1)}(r^2 U_l^m{}^{\prime }+2rU_l^m)Y_l^m \Big \rbrace . \end{eqnarray}

(217)

We now integrate this expression with respect to λ to give

\begin{eqnarray} \dot{\xi }_\lambda &=& \frac{5\pi \alpha ^2\omega }{16G}\frac{r j_2({ {\frac{\pi r}{R_*}}})}{\rho _0} \sin ^2\theta [ 2\sin (2\chi )\cos \theta \cos \lambda + \sin ^2\chi \sin \theta \cos (2\lambda ) ]\nonumber \\ &&+\, \frac{\omega }{B_0}\sum \limits _{l,m} \Big \lbrace W_l^m [lQ_{l+1}Y_{l+1}^m-(l+1)Q_l Y_{l-1}^m] +\frac{im}{l(l+1)}(r^2 U_l^m{}^{\prime }+2rU_l^m)Y_l^m \Big \rbrace . \end{eqnarray}

(218)

Finally, by following the reasoning of Section 5.3.1, we rewrite the three components of |$\dot{\xi }$| (equations 207, 210 and 218) in a manifestly real-valued form, with sums over l alone:

\begin{eqnarray} \dot{\xi }_r &=& -\frac{\omega r}{B_0\sin \theta }\sum \limits _{l=1}^\infty \bigg \lbrace \frac{U_{2l}^1}{2l(2l+1)} \Big \lbrace 2l(2l+1)Q_{2l+1}Q_{2l+2}(Y_{2l+2}^1+Y_{2l+2}^{-1}) - [1+2l(2l+2)Q_{2l+1}^2+(2l-1)(2l+1)Q_{2l}^2 ](Y_{2l}^1+Y_{2l}^{-1}) \nonumber \\ \,+\, 2l(2l+1)Q_{2l-1}Q_{2l} (Y_{2l-2}^1+Y_{2l-2}^{-1}) \Big \rbrace + \frac{U_{2l}^2}{2l(2l+1)} \Big \lbrace 2l(2l+1)Q_{2l+1}Q_{2l+2}(Y_{2l+2}^2-Y_{2l+2}^{-2}) \nonumber \\ -\, [4+2l(2l+2)Q_{2l+1}^2+(2l-1)(2l+1)Q_{2l}^2 ](Y_{2l}^2-Y_{2l}^{-2}) + 2l(2l+1)Q_{2l-1}Q_{2l} (Y_{2l-2}^2-Y_{2l-2}^{-2}) \Big \rbrace \bigg \rbrace , \end{eqnarray}

(219)

\begin{eqnarray} \dot{\xi }_\theta &=& -\frac{\omega }{B_0}\sum \limits _{l=1}^\infty \bigg \lbrace i W_{2l-1}^1 (Y_{2l-1}^1+Y_{2l-1}^{-1}) - \frac{(r^2 U_{2l}^1{}^{\prime }+2rU_{2l}^1)}{2l(2l+1)}[2lQ_{2l+1}(Y_{2l+1}^1+Y_{2l+1}^{-1})-(2l+1)Q_{2l}(Y_{2l-1}^1+Y_{2l-1}^{-1})]\nonumber \\ \,+\,2i W_{2l-1}^2 (Y_{2l-1}^2-Y_{2l-1}^{-2}) - \frac{(r^2 U_{2l}^2{}^{\prime }+2rU_{2l}^2)}{2l(2l+1)}[2lQ_{2l+1}(Y_{2l+1}^2-Y_{2l+1}^{-2})-(2l+1)Q_{2l}(Y_{2l-1}^2-Y_{2l-1}^{-2})]\bigg \rbrace , \end{eqnarray}

(220)

\begin{eqnarray} \dot{\xi }_\lambda &=& \ \frac{5\pi \alpha ^2\omega }{16G}\frac{r j_2({ {\frac{\pi r}{R_*}}})}{\rho _0} \sin ^2\theta [ 2\sin (2\chi )\cos \theta \cos \lambda + \sin ^2\chi \sin \theta \cos (2\lambda ) ]\nonumber \\ +\, \frac{\omega }{B_0}\sum \limits _{l = 1}^\infty \bigg \lbrace W_{2l-1}^1 [ (2l-1)Q_{2l}(Y_{2l}^1-Y_{2l}^{-1})-2lQ_{2l-1}(Y_{2l-2}^1-Y_{2l-2}^{-1})] + \frac{i(r^2 U_{2l}^1{}^{\prime }+2rU_{2l}^1)}{2l(2l+1)}(Y_{2l}^1-Y_{2l}^{-1}) \nonumber \\ &&+\, W_{2l-1}^2 [ (2l-1)Q_{2l}(Y_{2l}^2+Y_{2l}^{-2}) - 2lQ_{2l-1}(Y_{2l-2}^2+Y_{2l-2}^{-2})] + \frac{i(r^2 U_{2l}^2{}^{\prime }+2rU_{2l}^2)}{l(2l+1)}(Y_{2l}^2+Y_{2l}^{-2}) \bigg \rbrace . \end{eqnarray}

(221)

7 RESULTS

Finally, we come to the results of solving equations (164) and (174), with the boundary conditions given in Table 2. Having recast the original system of DAEs into a simpler system of ODEs, the numerical solution then becomes fairly straightforward, subject to a few numerical subtleties that we detail in Appendix B. In what follows, the radial functions are of course time-independent, but the results for |$\delta {\boldsymbol B}$| and |$\dot{\xi }$| are presented for t = 0. Over a period of 2π/ω, these three-dimensional figures would rotate around the z-axis.

7.1 The poloidal and toroidal radial functions

We recall that the fundamental quantities to solve for were the toroidal radial functions |$X_l^m\equiv W_l^m/r$| and the poloidal radial functions |$\bar{U}_l^m\equiv iU_l^m$|⁠. From these, one can find the perturbed magnetic field |$\delta {\boldsymbol B}$|⁠, by using equation (172), and the perturbed velocity field – the |$\xi$|-motions – by using equations (219)–(221). We begin by plotting these poloidal and toroidal radial functions, in Fig. 3. All of these functions appear to tend to zero at the origin, suggesting that the central boundary conditions described in Section 5.4 were actually unnecessary and that simply imposing the vanishing of all variables would have been sufficient. Having tried this, however, we found that the resulting functions were not solutions of the original system of DAEs, (164) and (165). Numerically, the U and X functions have tiny – but non-zero – values at the centre, and fixing the correct ratios between them (those of Table 2) is, in fact, crucial for getting a solution to the DAEs.

$Plots of the radial functions of $\delta {\boldsymbol B}$ against dimensionless radius $\hat{r}$, for lmax = 4, inclination angle χ = π/4, dimensionless primary rotation frequency $\hat{\alpha }= 0.1$ and dimensionless background magnetic field strength $\hat{\Lambda }= 0.1$ (see equations 229–231 for formulae to redimensionalize these to physical stellar parameters). Varying χ only changes the amplitude of the results, and the relative size of the m = 1 and 2 components; in the limit of an orthogonal rotator, only the m = 2 functions are non-zero.$

Figure 3.

Plots of the radial functions of |$\delta {\boldsymbol B}$| against dimensionless radius |$\hat{r}$|⁠, for l_max = 4, inclination angle χ = π/4, dimensionless primary rotation frequency |$\hat{\alpha }= 0.1$| and dimensionless background magnetic field strength |$\hat{\Lambda }= 0.1$| (see equations 229–231 for formulae to redimensionalize these to physical stellar parameters). Varying χ only changes the amplitude of the results, and the relative size of the m = 1 and 2 components; in the limit of an orthogonal rotator, only the m = 2 functions are non-zero.

Open in new tab Download slide

The m = 1 functions in Fig. 3 are numerically larger than their m = 2 counterparts for the chosen inclination angle of χ = π/4; only very close to the orthogonal-rotator case (χ = π/2) do the m = 2 perturbations become larger than those for m = 1. Perhaps the other most significant feature to note is that the |$\bar{U}_2^1$| function is a factor of 8 smaller than the |$\bar{U}_4^1$|⁠, i.e. the quadrupolar perturbed magnetic field is significantly smaller than the l = 4 component. This contradicts one's experience from simpler calculations – in axisymmetric magnetic equilibria, for example, the solutions feature odd l and are dominantly dipolar, with the power in the multipoles falling off monotonically (Ciolfi et al. 2009; Lander, Andersson & Glampedakis 2012).

The solutions plotted in Fig. 3 are numerical, with no closed-form expressions, and so are of limited applicability. For this reason, we next explore the scaling of the solutions with rotation rate, field strength and inclination angle, and present polynomial approximations to the radial functions |$\bar{U}_l^m$| and |$W_l^m$|⁠.

From equations (164) and (165), the source terms for the radial functions |$U^m_l,X^m_l$| are of the form |$\tilde{\Psi }_l^m/\Lambda$| and |$\tilde{\Upsilon }_l^m/\Lambda$|⁠; comparing with equations (A4) and (A5), we see that the m = 1 source terms are proportional to α²Λsin (2χ) and that the m = 2 source terms are proportional to α²Λsin ²χ. We therefore expect the same scaling for the radial functions |$U^m_l,X^m_l$|⁠, which we have confirmed by solving the equations for different values of α, Λ and χ (plots omitted for brevity).

To avoid very lengthy approximations, we will use four significant figures in the pre-factor for each power of r. For each radial function, we will employ the minimum number of powers of r required so that the deviation from the actual solution – defined as the difference between the approximate and actual solution, divided by the maximum value attained by the actual solution – never exceeds 2 per cent (it is typically below 1 per cent, in fact). With these demands, and using the ‘Fit’ function within mathematica, we find the following polynomial approximations to |$\bar{U}_l^m$| and |$W_l^m$|⁠, where we have terminated the angular sum at l_max = 4 as usual:

\begin{eqnarray} \bar{U}_2^1(\hat{r}) = \hat{\alpha }^2\hat{\Lambda }\sin (2\chi ) \hat{r}^2 (-7.616 + 20.24 \hat{r} - 13.84 \hat{r}^2 + 0.9159 \hat{r}^3 - 5.827 \hat{r}^4 + 6.186 \hat{r}^5), \end{eqnarray}

(222)

\begin{eqnarray} \bar{U}_4^1(\hat{r})= \hat{\alpha }^2\hat{\Lambda }\sin (2\chi ) \hat{r}^2 (121.4 - 639.5 \hat{r} + 1548 \hat{r}^2 - 2062 \hat{r}^3 + 1430 \hat{r}^4 - 397.7 \hat{r}^5), \end{eqnarray}

(223)

\begin{eqnarray} W_1^1(\hat{r}) = \hat{\alpha }^2\hat{\Lambda }\sin (2\chi ) \hat{r}^2 (188.6 - 1410 \hat{r} + 4340 \hat{r}^2 - 7531 \hat{r}^3 + 8342 \hat{r}^4 - 5643 \hat{r}^5 + 1720 \hat{r}^6), \end{eqnarray}

(224)

\begin{eqnarray} W_3^1(\hat{r}) = \hat{\alpha }^2\hat{\Lambda }\sin (2\chi ) \hat{r}^3 (81.58 - 474.4 \hat{r} + 1245 \hat{r}^2 - 1819 \hat{r}^3 + 1409 \hat{r}^4 - 442.6 \hat{r}^5), \end{eqnarray}

(225)

\begin{eqnarray} \bar{U}_2^2(\hat{r}) = \hat{\alpha }^2\hat{\Lambda }\sin ^2\chi \hat{r} (0.2287 + 0.3832 \hat{r} - 1.611 \hat{r}^2 + 2.844 \hat{r}^3 - 1.779 \hat{r}^4), \end{eqnarray}

(226)

\begin{eqnarray} \bar{U}_4^2(\hat{r}) = \hat{\alpha }^2\hat{\Lambda }\sin ^2\chi \hat{r}^2 (-0.082\,30 + 1.079 \hat{r} - 2.798 \hat{r}^2 + 4.073 \hat{r}^3 - 2.286 \hat{r}^4), \end{eqnarray}

(227)

\begin{eqnarray} W_3^2(\hat{r}) = \hat{\alpha }^2\hat{\Lambda }\sin ^2\chi \hat{r}^3 (0.0595 - 0.2124 \hat{r} + 0.9521 \hat{r}^2 - 0.8182 \hat{r}^3). \end{eqnarray}

(228)

Note that the number of terms required for an accurate approximation of each function differ, and that the minimum and maximum powers featured in the above set of relations are r and r⁸, respectively. Constant (i.e. r⁰) terms are unnecessary, since the U and X functions are very small at the origin, and the W functions are exactly zero (since |$W_l^m = rX_l^m$|⁠, and the |$X_l^m$| functions are finite at r = 0); see Fig. 3.

These functions may be redimensionalized by multiplying by the appropriate combinations of G, ρ_c and R_*, to give values for a particular star. For a neutron star with mass 1.4 M_⊙ and radius R_* = 10 km, one can redimensionalize the above functions to chosen values of rotation rate and field strength using the replacements

\begin{eqnarray} \hat{\alpha } = \frac{\alpha }{\sqrt{G\rho _{\rm c}}} = \frac{\alpha }{12100\,\,{\rm rad\,\, s}^{-1}} = \frac{\alpha /2\pi }{1920\,{\rm Hz}}, \end{eqnarray}

(229)

\begin{eqnarray} \hat{\Lambda } = \frac{\Lambda }{\sqrt{G}} = \frac{\langle B\rangle }{7.94\times 10^{16}\,{\rm G}}, \end{eqnarray}

(230)

\begin{eqnarray} \hat{r} = \frac{r}{R_*} = \frac{r}{10^6\,{\rm cm}}, \end{eqnarray}

(231)

where we have used the central density ρ_c = 2.19 × 10¹⁵ g cm⁻³ of a γ = 2 polytrope with the above physical mass and radius (Lander & Jones 2009), and have rewritten Λ in terms of the average field strength 〈B〉 by combining equation (62) with the ellipticity formula (38) from Lander & Jones (2009).

7.2 Behaviour of higher multipoles at the centre

As discussed in Section 5.4, we are limited to truncating our solutions at a maximum angular index l = l_max = 4, and we anticipate that higher multipoles will gradually become negligible. However, we have already seen from Fig. 3 that the magnitude of the multipolar functions does not decrease monotonically, and so we would like some idea of the relative strength of higher multipoles.

Returning to our system of DAEs, (164) and (165), we note a few features. They couple together X functions over an l-range of 4, from X_l − 2 to X_l + 2, and |$\bar{U}$| functions over an l-range of 6, so it is natural to expect solutions to have significant contributions from many multipoles. The highest multipoles in the source terms, on the other hand, are an l = 4 piece of |$\tilde{\Psi }$| and an l = 3 piece of |$\tilde{\Upsilon }$| (see Section A). From the way the source terms enter the DAEs, this means that the l = 5 system of equations is the last to contain any explicit source, through the |$\tilde{\Upsilon }_{l-2}$| term of equation (165). For all higher values of l, the DAEs merely relate different |$\bar{U}$| and X functions to one another, and the higher l functions are only non-zero by virtue of coupling to their lower l counterparts.

In lieu of solving the full system of DAEs for high l_max, we may evaluate the equations at the centre of the star – as done for the second DAE, (165), in Section 5.4. For r = 0, all source terms vanish (as in the full DAEs at high l), and we simply have a system of algebraic simultaneous equations relating the values of different functions at the centre, |$\bar{U}_l(0)$| and X_l(0). When truncating for odd values of l_max, the equations may be ‘solved’ to give each quantity as a fraction of X₁(0), thus giving some idea of the way different multipoles couple together, and their behaviour at high l. For lower values of l_max, the results show some variation with the exact truncation value, so we terminate at l_max = 101, well after the point at which any variation can be seen.

We plot the results of this experiment in Fig. 4, where we see that the quantities X_l(0)/X₁(0) and |$\bar{U}_l(0)/X_1(0)$| decay in an oscillatory manner with increasing l, alternating between negative and positive values, rather than in a more predictable monotonic manner. We also see that the U_l(0) quantities decay far more slowly, with even the l = 24 multipole being more than 1 per cent of the value of X₁(0); by contrast, every X_l(0) quantity with l > 3 is below 1 per cent of X₁(0).

$Setting the radius to zero in the original system of DAEs, (164) and (165), results in a system of simultaneous algebraic equations. We truncate these at lmax = 101 to solve them, plotting the results for l ≤ 25 here. As solutions, we find the relative magnitudes of different multipoles $\bar{U}_l(0),X_l(0)$ at the centre, as a fraction of the value of the dipole piece X1(0) there (blue crosses). This suggests that the U functions fall off more slowly with increasing l than the X functions. For reference, the maximum values of U2, U4 and X3 from our full solution are plotted as ratios of the maximum value of X1 (filled red circles).$

Figure 4.

Setting the radius to zero in the original system of DAEs, (164) and (165), results in a system of simultaneous algebraic equations. We truncate these at l_max = 101 to solve them, plotting the results for l ≤ 25 here. As solutions, we find the relative magnitudes of different multipoles |$\bar{U}_l(0),X_l(0)$| at the centre, as a fraction of the value of the dipole piece X₁(0) there (blue crosses). This suggests that the U functions fall off more slowly with increasing l than the X functions. For reference, the maximum values of U₂, U₄ and X₃ from our full solution are plotted as ratios of the maximum value of X₁ (filled red circles).

Open in new tab Download slide

The only rigorous statement that one can make about Fig. 4 is, of course, that it shows the solutions to the DAEs with r set to zero. Nonetheless, we believe that it also captures the essence of how higher multipoles couple together in our problem. Evidence in favour of this is the fact that it predicts the magnitudes of |$\bar{U}_2^{\rm max}/X_1^{\rm max}$|⁠, |$\bar{U}_4^{\rm max}/X_1^{\rm max}$| and |$X_3^{\rm max}/X_1^{\rm max}$|⁠, from our l_max = 4 solutions, in the correct order, and in particular the fact that |$\bar{U}^1_2$| is far smaller than |$\bar{U}_4^1$| and has the opposite sign.

7.3 Magnetic-field perturbations

Having found solutions for U_l and W_l, we may now insert these into equation (172) to find |$\delta {\boldsymbol B}$|⁠. The resulting field configuration is complicated, and we have found it clearest to separate it into its poloidal and toroidal components. First, in Fig. 5, we plot the magnitudes of these components. Broadly speaking, the perturbed magnetic field is strongest towards the centre, and the toroidal component is about double the strength of the poloidal one. For most values of the inclination angle, the m = 1 component dominates over m = 2, so the near-aligned (χ = π/16) and near-orthogonal (χ = 7π/16) results are visibly similar. For χ = π/2, when the magnetic and rotation axes are orthogonal, there is only an m = 2 component. However, this situation is actually stationary (and therefore does not correspond to precession); from equation (7), we see that ω = 0 for χ = π/2, and consequently (by equation 201), |$\mathrm{\partial} (Y_l^m(\theta ,\lambda ))/\mathrm{\partial} t = 0$|⁠.

Figure 5.

Three-dimensional contour plots showing the magnitude of the perturbed magnetic field, separated into its poloidal (top row) and toroidal (bottom row) components. Note that viewing angles differ by π/2 for the poloidal and toroidal plots. The results are symmetric about the z = 0 plane, so we show only the Northern hemisphere. Results are shown for three different values of inclination angle χ. Contours are colour coded so that the blue contour represents a field strength double that of the light brown contour, green is double the value of blue, and red is double the value of green. For a fixed inclination angle, the contour colouring is fixed, so that poloidal and toroidal components may be compared directly, but different inclination angles use different keys. The m = 1 field vanishes for the case of an orthogonal rotator (χ = π/2), leaving only the m = 2 component, and resulting in a very different-looking field (which is actually time-independent). However, the m = 2 component is generally weak, so that even the near-orthogonal configuration (χ = 7π/16) is similar to the near-aligned model (χ = π/16).

Open in new tab Download slide

Next, we show the direction of the magnetic field, in Fig. 6. The clearest visualization we could find was to plot the poloidal and toroidal unit vectors over a nested sequence of constant-radius shells with r/R_* = 0.2, 0.4, 0.6 and 0.8. Comparing the field across consecutive shells should, hopefully, allow the reader to visualize the full three-dimensional field-line structure. The toroidal component is clearer, since (by definition) its field lines run along constant-r shells. For the poloidal component, note that the blue and white patches are in roughly the same position on each contour, meaning that a given field line moves towards the centre of the star through the blue regions, then bends back and comes out through a neighbouring white region to its ‘starting’ location, thus forming a closed loop (as dictated by the requirement that |$\nabla \cdot \delta {\boldsymbol B}= 0$|⁠).

Figure 6.

The direction of the perturbed magnetic field, plotted along a series of hemispherical shells beginning near the centre and ending near the surface. From left to right, the shells represent r/R_* = 0.2, 0.4, 0.6 and 0.8, respectively. The upper row of plots shows the unit poloidal vector and the lower row the unit toroidal vector. Since arrows into and out of contours are hard to distinguish, we use colour coding for the strength of the radial field component: blue represents a region of field lines pointing towards the centre of the star, whilst white indicates outwards-pointing field lines. This only applies to the poloidal-field plots, as toroidal fields are, by definition, non-radial, and so their arrows are all tangential to the hemispherical shells. Each poloidal/toroidal plot is viewed from the same (x, y) position as in Fig. 5, but from above – i.e. for a z > 0 position, unlike the z < 0 vantage points of Fig. 5. To save space, the axes and bounding box are only drawn for the left-hand plots. Finally, we take χ = π/16 for all plots, though there is little variation of the field's direction with inclination angle.

Open in new tab Download slide

7.4 |$\xi$|-motions

At last we come to the |$\xi$|-motions: the non-rigid velocity field present in a precessing fluid star, and the original goal of this paper. We recall that these may now be calculated from our solutions |$\bar{U}_l,X_l$| using equations (219)–(221). As for |$\delta {\boldsymbol B}$|⁠, we separate the information about |$\dot{\xi }$| into two figures, one showing the magnitude of |$\dot{\xi }$| and the other its direction. The contours of |$\dot{\xi }$|⁠, shown in Fig. 7, represent smaller differences than those used above for |$\delta {\boldsymbol B}$| – the speed |$|\dot{\xi }|$| only differs by a factor of |$\sqrt{2}$| between neighbouring contours, and therefore only varies by a factor of 2 across most of the stellar interior. Like the perturbed toroidal field from the lower panels in Fig. 5, the speed is greatest in a torus centred on the origin and aligned along the x = 0 plane. For the near-aligned case (χ = π/16), |$\dot{\xi }$| is approximately reflection symmetric about the x = 0 plane, but as the inclination angle is increased, the speed lowers more quickly in the x < 0 region (χ = 7π/16 plot). Finally, as mentioned in the previous section, |$\dot{\xi }= {0}$| for the limiting case of χ = π/2.

$Three-dimensional contour plots showing the magnitude of $\dot{\xi }$, for χ = π/16 (left-hand plot) and χ = 7π/16 (right-hand plot). The details are the same as for Fig. 5, except that each contour represents a value $\sqrt{2}$ times that of the last, meaning that the innermost contour has twice the value of the outermost.$

Figure 7.

Three-dimensional contour plots showing the magnitude of |$\dot{\xi }$|⁠, for χ = π/16 (left-hand plot) and χ = 7π/16 (right-hand plot). The details are the same as for Fig. 5, except that each contour represents a value |$\sqrt{2}$| times that of the last, meaning that the innermost contour has twice the value of the outermost.

Open in new tab Download slide

The direction of |$\dot{\xi }$| is shown in Fig. 8. Close to the centre and to the outer boundary, the velocity is roughly tangential to constant-r shells, but in between, it has a large radial component. As for the poloidal-field direction plots, the blue and white patches are in the same location for each of the four shells in Fig. 8, meaning that there is a meridional circulation of fluid, with streamlines moving towards the centre of the star through the blue regions and back out through the white regions. Unlike the magnetic-field case, however, these streamlines do not need to close, since |$\nabla \cdot \dot{\xi }\ne 0$| (i.e. the flow is compressible). Not restricting to incompressible flow was, in fact, the place at which our approach first diverged from that of Mestel & Takhar (1972); recall the discussion of Section 2. For this reason, it is logical to close this section by looking at the behaviour of |$\nabla \cdot \dot{\xi }$| itself.

$The direction of the $\dot{\xi }$ vector field on constant-radius shells of r/R* = 0.2, 0.4, 0.6 and 0.8 (left to right), for χ = π/16. See the caption of Fig. 6 for details.$

Figure 8.

The direction of the |$\dot{\xi }$| vector field on constant-radius shells of r/R_* = 0.2, 0.4, 0.6 and 0.8 (left to right), for χ = π/16. See the caption of Fig. 6 for details.

Open in new tab Download slide

In Fig. 9, we plot the divergence of |$\nabla \cdot \dot{\xi }$|⁠, normalized by |$|\dot{\xi }|/r$| to give a dimensionless quantity. In a barotropic star, one cannot assume incompressible motions, and indeed we see that |$\nabla \cdot \dot{\xi }$| deviates greatly from zero, especially in the outer regions. Thermal- or composition-gradient stratification in real stars would provide some additional buoyancy force to suppress compressible motions, but for a given model, one still needs to account for such a buoyancy force consistently, and to verify it is sufficiently strong to impose |$\nabla \cdot \dot{\xi }= 0$|⁠. Fig. 9 gives an idea of how and where such forces would have to act to ensure |$\nabla \cdot \dot{\xi }= 0$|⁠.

$Divergence of $\dot{\xi }$, normalized by the speed $|\dot{\xi }|$ divided by the radius, to make a dimensionless quantity. We show the behaviour of this quantity for λ = π/2 and three values of θ, as labelled. This plot demonstrates that for barotropic stellar models, the higher order motions in a precessing fluid star are, in fact, highly compressible. It also gives an idea of how strong buoyancy forces would need to be in order to produce the incompressible motions discussed in Mestel & Takhar (1972).$

Figure 9.

Divergence of |$\dot{\xi }$|⁠, normalized by the speed |$|\dot{\xi }|$| divided by the radius, to make a dimensionless quantity. We show the behaviour of this quantity for λ = π/2 and three values of θ, as labelled. This plot demonstrates that for barotropic stellar models, the higher order motions in a precessing fluid star are, in fact, highly compressible. It also gives an idea of how strong buoyancy forces would need to be in order to produce the incompressible motions discussed in Mestel & Takhar (1972).

Open in new tab Download slide

8 DISCUSSION

This paper has been concerned with a detailed description of the dynamics of a rotating magnetized star, in the general ‘oblique rotator’ case when the magnetic and rotational axes are misaligned by some angle 0 < χ < π/2. The star's bulk motion is that of free precession, but sustaining this in a fluid star requires additional dynamics in the interior: time-dependent velocity and magnetic fields. The concept of non-rigid precession was first described by Mestel & Takhar (1972), who produced simplified models of the perturbed field |$\dot{\xi }$| using first-order perturbation equations. A new calculation of |$\dot{\xi }$| was presented in Mestel et al. (1981), who also used this, in turn, to calculate magnetic field perturbations, given a poloidal background field. Nittmann & Wood (1981) performed the same analysis for a toroidal background field.

In this paper, we argue that a consistent solution for both |$\dot{\xi }$| and |$\delta {\boldsymbol B}$| can only come from solving a higher order set of equations than those previously considered. A large part of the work presented here involves a detailed formalism to describe these perturbations (Section 3), algebraic manipulation to reduce them to a set of DEs in radial functions (Section 5) and then solving these (Section 7 and Appendix B).

In order to make the problem tractable, we have had to make a number of assumptions, and now need to consider how these affect the validity of our model when applied to different classes of star. None the less, we emphasize that our model contains the minimum essential physics for non-rigid precession, and so everything that follows concerns additional effects. Even when our assumptions are unjustifiable, then, the work presented here can be regarded as a ‘background’ model upon which to build.

The most major simplification we have made is assuming a toroidal background magnetic field. Purely toroidal fields, like purely poloidal fields, suffer rapidly growing dynamical instabilities – and so cannot be realistic models of stellar magnetic fields (Tayler 1973), which must instead feature both poloidal and toroidal components. However, in the absence of an additional stabilizing mechanism, there are indications that even mixed poloidal–toroidal fields are generically unstable (Lander & Jones 2012; Mitchell et al. 2015). On the other hand, toroidal-dominated stable equilibria have been found for stellar models with an additional thermal-buoyancy force (Braithwaite 2009). The issue of stable magnetic equilibria is not fully settled, but if the toroidal component is generically the dominant one in stars, our results may indeed be a reasonable first approximation to their precessional dynamics.

For a stellar model with a purely poloidal field, we believe that an analogous analysis to ours could be carried out. However, whilst a toroidal field has a component only in the direction of one spherical–polar coordinate (the azimuthal one, λ), a poloidal field's depends on the other two (r and θ), suggesting that spherical–polar coordinates will not be appropriate for studying this problem. Instead, one can work in a coordinate system fitted to poloidal-field lines (Bernstein et al. 1958), in which case all the steps of our toroidal-field calculation should be tractable for a poloidal field too. There is, however, no coordinate system in which a mixed poloidal–toroidal field becomes similarly simple (see e.g. Tayler 1980); for a stellar model with such a background field, it seems inevitable that a fully numerical approach would be required.

Neutron stars: Most of the assumptions made in this paper are typical ones in the context of neutron-star modelling. This is certainly not equivalent to saying that they are all justifiable on physical grounds, however, so each needs further attention. The γ = 2 polytrope we adopt is quite suitable for a single-fluid model of a neutron star, but the cores of mature neutron stars are better described by a two-fluid equation of state in which the neutrons are superfluid; simple magnetic equilibrium models for mature neutron-star cores may be found in which each fluid component obeys its own polytropic relation (Glampedakis, Andersson & Lander 2012; Lander et al. 2012). The neutron superfluid will decouple from the precessional dynamics we consider in this paper, and so only the charged component is relevant. This component is, however, better described by a polytropic index γ ≈ 3/2.

It is reasonable to assume that no additional buoyancy forces enter the precessional dynamics of a mature neutron star. A neutron star in equilibrium has a density-dependent fraction of neutrons to protons. When the neutrons are not superfluid, displacing a fluid element results in chemical reactions that convert between neutrons and protons to restore the equilibrium fraction of protons and neutrons, and this represents a buoyancy force associated with the composition-gradient stratification (Reisenegger & Goldreich 1992). However, these reactions are strongly suppressed (along with the associated buoyancy) when the neutrons are superfluid. In addition, neutron stars become isothermal very early on in their life, so there is negligible thermal-gradient stratification. Finally, we have already argued (in Section 5.7) that crustal elasticity and the matching between crust and core could seriously modify our solutions in the outermost region, 0.9R_* ≤ r ≤ R_*, and so we considered it wisest to exclude this region altogether, and place the outer boundary for our solutions at a radius of 0.9R_*.

The most significant neutron-star physics we have failed to account for is the superconducting nature of the stellar core. Since most of this region is believed to form a type-II superconductor, however, the magnetic field will still penetrate the core and produce a macroscopically uniform magnetic force (Baym, Pethick & Pines 1969). The conceptual picture should thus remain intact: the magnetic field lends rigidity to the fluid, inducing bulk precession but at the expense of requiring additional internal motions and magnetic-field perturbations. Quantitatively, however, the nature of these perturbations will undoubtedly be different – and qualitatively we would expect them to scale with the superconducting magnetic force⁹ ∝H_c1B ∼ 10¹⁵B instead of the normal Lorentz force ∝B² (Easson & Pethick 1977; Wasserman 2003).

The velocity field |$\dot{\xi }$| represents a ‘macroscopic’ coupling between the star's rotation and magnetic field. In a neutron-star core with type-II superconducting protons and superfluid neutrons, however, there is a competing ‘microscopic’ coupling mechanism: the interaction of rotational neutron vortices with magnetic proton fluxtubes (Ruderman, Zhu & Chen 1998; Link 2003). Ascertaining the relative importance of these two mechanisms would be an important step towards a quantitative picture of non-rigid neutron-star precession.

White dwarfs: Unlike neutron stars, white dwarfs form a broad class of stars, with differences in chemical composition and the physics governing their different regions. The major sub-division of white dwarfs of relevance to this paper is between the magnetic and non-magnetic white dwarfs. Like neutron stars, magnetic white dwarfs have field strengths spanning several orders of magnitude, but with an upper limit of ≲10⁹ G. Many of these strongly magnetized objects have, like magnetars, long rotation periods, although there are some rapidly rotating and highly magnetized white dwarfs (Ferrario et al. 2015) – the ideal combination for relatively short precession periods.

The equation of state for a white dwarf is a degenerate Fermi gas. In the high-density limit, this reduces to a γ = 4/3 polytrope, and in the low-density limit to a γ = 5/3 polytrope (Nauenberg 1972). For lower densities, then, the γ = 2 relation adopted here is probably not an egregious approximation to white dwarf matter. Imposing an outer boundary below the surface is also justifiable, since the degenerate core of white dwarfs gives way abruptly to a low-mass non-degenerate envelope region – although the envelope is only ∼30 km (Althaus et al. 2010), meaning that a more appropriate radius to impose the outer boundary at would be ∼0.997R_* (taking care that the perturbative scheme is still applicable at this radius; see Section 5.6).

In this work, we have assumed a fluid region that extends to the stellar centre, which is not always the case for white dwarfs: As they become older and colder, white dwarfs begin to form solid cores (Van Horn 1968), with the crystallization proceeding from the centre outwards and solidifying the whole star over a time-scale of ∼10⁹ yr (Iben & Tutukov 1984). Elasticity will then affect the induced velocity and magnetic-field perturbations in the solid core; for a sufficiently rigid core lattice, these perturbations would become negligible, and the core's motion would then be well described by standard rigid-body free precession. Finally, more detailed modelling of white dwarf precession should account for the fact that these stars are non-barotropic, with stratification due to both composition and thermal gradients.

Non-degenerate (main-sequence) stars: With some caveats, our analysis applies to objects along the upper part of the main sequence of non-degenerate stars. Upper main-sequence stars have higher masses and more simple field topologies than other non-degenerate stars, and since their envelopes appear to be in radiative equilibrium (rather than convective), dynamo action seems unlikely. The conclusion is that the magnetic fields in these envelopes are likely to be in dynamically stable equilibrium states. Of these stars, the most highly magnetized are the chemically peculiar A- and B-type stars, with field strengths up to ∼10⁴ G. Upper main-sequence stars with strong magnetic fields tend to rotate far more slowly than their weakly magnetized counterparts, such as the normal A and B stars (Donati & Landstreet 2009).

In the radiative envelopes of main-sequence stars, our analysis is broadly applicable, but as for white dwarfs, it should be modified by the inclusion of a thermal buoyancy force. This will constrain the form of |$\dot{\xi }$| and |$\delta {\boldsymbol B}$|⁠, but in a way that can only be confirmed by including thermal effects in a self-consistent matter in the whole analysis. As for partially crystallized white dwarfs, our analysis will not be valid all the way to the centre of main-sequence stars. In this case, the reason is the presence of an inner core region where vigorous convection disrupts and regenerates the magnetic field – and so our assumption of a large-scale stable background magnetic field becomes invalid here. Moving to the outer boundary, there is no distinct outer region of the radiative zone where one would expect the |$\xi$|-motions to stop. Our approach requires us to cut off the solution before reaching the stellar surface R_* (see Section 5.6), but as long as we choose an outer boundary close to R_*, we will only be neglecting a small, low-density region of the star.

For all classes of star, the internal velocity and magnetic fields solved for in this paper may have significant effects on their evolution and dynamics. In particular, since they couple the star's rotation and magnetic fields, the dissipation of these motions could cause the evolution of the inclination angle (Spitzer 1958). For toroidal fields - the case considered here - the dissipation will tend to cause the star to become an orthogonal rotator, χ → π/2; for poloidal fields, χ → 0 (Mestel & Takhar 1972). For the more realistic case of mixed poloidal–toroidal fields the evolution of χ will be dictated by the dominant field component,¹⁰ and so observations of this evolution could allow one to infer details of the star's internal magnetic field. Although we do not intend to tackle a fully self-consistent evolution of dissipative processes, a qualitative discussion of this and applications to astrophysics will be presented in a follow-up paper.

9 SUMMARY

We have presented a formalism and solution method for models of rotating magnetic stars with misaligned rotation and magnetic axes – sometimes known as oblique rotators. We began with arguably the simplest necessary pieces of physics to describe this problem: a rigid primary rotation of the star, and a large-scale m = 0 dipolar toroidal magnetic field. Despite this simplicity, we have found that the dynamics of such stars involve internal velocity and magnetic-field perturbations with highly complex geometries. The perturbed velocity and magnetic field are both stronger towards the centre of the star, and the toroidal perturbed field is a factor of 2 more intense than the poloidal. There are indications that whilst the perturbed toroidal field is dominated by its dipole (l = 1) piece, the poloidal component of the magnetic-field perturbation may have significant contributions from many higher multipoles – and, in fact, its l = 4 component is stronger than its lowest order multipole, a quadrupole (l = 2) piece. The velocity field has significant azimuthal and meridional components, and deviates strongly from being solenoidal, i.e. the internal motions are highly compressible. The secular dissipation of these perturbations will cause the star to evolve over time into either an aligned or orthogonal rotator, and so the distribution of inclination angles in stars is likely to be connected with their internal magnetic fields.

Acknowledgments

We acknowledge support from STFC via grant number ST/M000931/1, and from NewCompStar (a COST-funded Research Networking Programme). We thank the anonymous referee for a number of helpful suggestions for improving this paper.

1

That is, a field that does not require continuous dynamo action to maintain it, in the way that the Sun's field does.

2

referred to as ‘nutation’ in some studies.

3

This would be axisymmetric in the |${\boldsymbol e}_z^{(\alpha )}$|-coordinate system, but is non-axisymmetric in coordinates referred to the axis |${\boldsymbol e}_z^{(B)}$|⁠.

4

By a more systematic analysis later in this paper, we calculate δρ_α, and its angular dependence does indeed agree with Mestel & Takhar (1972).

5

The poloidal-field scalar function is conventionally denoted by Φ; we use ϒ instead, to avoid confusion with the gravitational potential.

6

These quantities also diverge for γ = 4/3 – often used to model the pressure–density relation for main-sequence stars – and for γ = 5/3, a typical white dwarf model. This already hints at a more general problem with the model.

7

‘Break’ in this context is likely to be a plastic deformation rather than a brittle fracture (Jones 2003).

8

More accurately, only the charged component of the core (protons and electrons) would couple to the crust, not the neutron superfluid component.

9

H_c1 is the lower critical field for a type-II superconductor and is associated with the microscopic field along fluxtubes.

10

The possibility of no evolution of χ can probably be dismissed, since it would require very fine tuning of the field components.

REFERENCES

Althaus

L. G.

Córsico

A. H.

Isern

J.

García-Berro

E.

,

2010

,

A&AR

,

18

,

471

Crossref

Search ADS

Arfken

G. B.

Weber

H. J.

,

2005

,

Mathematical Methods for Physicists

.

Elsevier Academic Press

,

Oxford

Google Scholar

Google Preview

OpenURL Placeholder Text

WorldCat

Backus

G.

,

1986

,

Rev. Geophys.

,

24

,

75

Crossref

Search ADS

Baym

G.

Pethick

C.

Pines

D.

,

1969

,

Nature

,

224

,

673

Crossref

Search ADS

Bernstein

I. B.

Frieman

E. A.

Kruskal

M. D.

Kulsrud

R. M.

,

1958

,

Phil. Trans. R. Soc. A

,

244

,

17

Braithwaite

J.

,

2009

,

MNRAS

,

397

,

763

Crossref

Search ADS

Chandrasekhar

S.

,

1939

,

An Introduction to the Study of Stellar Structure

.

Univ. Chicago Press

,

Chicago

Google Scholar

Google Preview

OpenURL Placeholder Text

WorldCat

Chandrasekhar

S.

Fermi

E.

,

1953

,

ApJ

,

118

,

116

Crossref

Search ADS

Chandrasekhar

S.

Prendergast

K. H.

,

1956

,

Proc. Natl. Acad. Sci. USA

,

42

,

5

Crossref

Search ADS

Ciolfi

R.

Ferrari

V.

Gualtieri

L.

Pons

J. A.

,

2009

,

MNRAS

,

397

,

913

Crossref

Search ADS

Dall'Osso

S.

Shore

S. N.

Stella

L.

,

2009

,

MNRAS

,

398

,

1869

Crossref

Search ADS

Donati

J.-F.

Landstreet

J. D.

,

2009

,

ARA&A

,

47

,

333

Crossref

Search ADS

Donati

J.-F.

Babel

J.

Harries

T. J.

Howarth

I. D.

Petit

P.

Semel

M.

,

2002

,

MNRAS

,

333

,

55

Crossref

Search ADS

Easson

I.

Pethick

C. J.

,

1977

,

Phys. Rev. D

,

16

,

275

Crossref

Search ADS

Ferrario

L.

de Martino

D.

Gänsicke

B. T.

,

2015

,

Space Sci. Rev.

,

191

,

111

Crossref

Search ADS

Glampedakis

K.

Andersson

N.

Lander

S. K.

,

2012

,

MNRAS

,

420

,

1263

Crossref

Search ADS

Horowitz

C. J.

Kadau

K.

,

2009

,

Phys. Rev. Lett.

,

102

,

191102

Crossref

Search ADS

PubMed

Iben

I.

Tutukov

A. V.

,

1984

,

ApJ

,

282

,

615

Crossref

Search ADS

Jones

P. B.

,

2003

,

ApJ

,

595

,

342

Crossref

Search ADS

Jones

D. I.

Andersson

N.

,

2002

,

MNRAS

,

331

,

203

Crossref

Search ADS

Landau

L. D.

Lifshitz

E. M.

,

1976

,

Mechanics

.

Elsevier Butterworth-Heinemann

,

Oxford

Google Scholar

Google Preview

OpenURL Placeholder Text

WorldCat

Lander

S. K.

Jones

D. I.

,

2009

,

MNRAS

,

395

,

2162

Crossref

Search ADS

Lander

S. K.

Jones

D. I.

,

2012

,

MNRAS

,

424

,

482

Crossref

Search ADS

Lander

S. K.

Jones

D. I.

Passamonti

A.

,

2010

,

MNRAS

,

405

,

318

Lander

S. K.

Andersson

N.

Glampedakis

K.

,

2012

,

MNRAS

,

419

,

732

Crossref

Search ADS

Landstreet

J. D.

,

1992

,

A&AR

,

4

,

35

Crossref

Search ADS

Lasky

P. D.

Glampedakis

K.

,

2016

,

MNRAS

,

458

,

1660

Crossref

Search ADS

Levin

Y.

D'Angelo

C.

,

2004

,

ApJ

,

613

,

1157

Crossref

Search ADS

Link

B.

,

2003

,

Phys. Rev. Lett.

,

91

,

10

Crossref

Search ADS

Mestel

L.

Takhar

H. S.

,

1972

,

MNRAS

,

156

,

419

Crossref

Search ADS

Mestel

L.

Nittmann

J.

Wood

W. P.

Wright

G. A. E.

,

1981

,

MNRAS

,

195

,

979

Crossref

Search ADS

Mitchell

J. P.

Braithwaite

J.

Reisenegger

A.

Spruit

H.

Valdivia

J. A.

Langer

N.

,

2015

,

MNRAS

,

447

,

1213

Crossref

Search ADS

Monaghan

J. J.

,

1976

,

Ap&SS

,

40

,

385

Crossref

Search ADS

Nauenberg

M.

,

1972

,

ApJ

,

175

,

417

Crossref

Search ADS

Nittmann

J.

Wood

W. P.

,

1981

,

MNRAS

,

196

,

491

Crossref

Search ADS

Petzold

L.

,

1982

,

SIAM J. Sci. Stat. Comput.

,

3

,

367

Crossref

Search ADS

Reisenegger

A.

,

2009

,

A&A

,

499

,

557

Crossref

Search ADS

Reisenegger

A.

Goldreich

P.

,

1992

,

ApJ

,

395

,

240

Crossref

Search ADS

Rookyard

S. C.

Weltevrede

P.

Johnston

S.

,

2015

,

MNRAS

,

446

,

3356

Crossref

Search ADS

Roxburgh

I. M.

,

1966

,

MNRAS

,

132

,

347

Crossref

Search ADS

Ruderman

M.

Zhu

T.

Chen

K.

,

1998

,

ApJ

,

492

,

267

Crossref

Search ADS

Schmidt

G. D.

Norsworthy

J. E.

,

1991

,

ApJ

,

366

,

270

Crossref

Search ADS

Spitzer

L.

,

1958

,

Lehnert

B.

,

Proc. IAU Symp. 6, Electromagnetic Phenomena in Cosmical Physics

.

Cambridge Univ. Press

,

Cambridge

, .

169

Google Scholar

Google Preview

OpenURL Placeholder Text

WorldCat

Tauris

T. M.

Manchester

R. N.

,

1998

,

MNRAS

,

298

,

625

Crossref

Search ADS

Tayler

R. J.

,

1973

,

MNRAS

,

161

,

365

Crossref

Search ADS

Tayler

R. J.

,

1980

,

MNRAS

,

191

,

151

Turolla

R.

Zane

S.

Watts

A. L.

,

2015

,

Rep. Prog. Phys.

,

78

,

116901

Crossref

Search ADS

PubMed

Van Horn

H. M.

,

1968

,

ApJ

,

151

,

227

Crossref

Search ADS

Wasserman

I.

,

2003

,

MNRAS

,

341

,

1020

Crossref

Search ADS

APPENDIX A: SPHERICAL-HARMONIC DECOMPOSITIONS OF SOURCE TERMS

Since we solve equations for a given m, we need the source terms decomposed into |$Y_l^m$| pieces too. We use the following relations between trigonometric functions and spherical harmonics |$Y_l^m(\theta ,\lambda )$|⁠:

\begin{eqnarray} &&{\sin (\theta )\sin \lambda = i\sqrt{\frac{2\pi }{3}}(Y_1^1+Y_1^{-1}),} \end{eqnarray}

(A1)

\begin{eqnarray} &&{\sin ^2\theta \sin (2\theta )\cos \lambda = \frac{8}{21}\sqrt{\frac{\pi }{5}}(Y_4^1-Y_4^{-1}) - \frac{8}{7}\sqrt{\frac{2\pi }{15}}(Y_2^1-Y_2^{-1}),} \end{eqnarray}

(A2)

\begin{eqnarray} &&{\sin ^4\theta \cos (2\lambda ) = -\frac{4}{21}\sqrt{\frac{2\pi }{5}}(Y_4^2+Y_4^{-2}) + \frac{4}{7}\sqrt{\frac{6\pi }{5}}(Y_2^2+Y_2^{-2}),} \end{eqnarray}

(A3)

together with those given in equations (76), (77), (118) and (119). These allow us to decompose the scalar functions |$\Upsilon _{\delta \mathcal {L}}$| and |$\Psi _{\delta \mathcal {L}}$| from equations (132) and (133), which enter as source terms in our final system of equations:

\begin{equation} \Upsilon _{\delta \mathcal {L}} = \frac{i\alpha ^2\Lambda ^2 r}{16G}\sqrt{\frac{\pi }{6}}\left\lbrace \sin (2\chi )\left[j_2\left({ {\frac{\pi r}{R_*}}}\right)-16k_{B}\right](Y_1^1+Y_1^{-1}) + \frac{1}{3\sqrt{7}} j_2\left({ {\frac{\pi r}{R_*}}}\right) \left[ \sqrt{2}\sin (2\chi )(Y_3^1+Y_3^{-1}) - \sqrt{5}\sin ^2\chi (Y_3^2-Y_3^{-2}) \right] \right\rbrace \end{equation}

(A4)

and

\begin{eqnarray} \Psi _{\delta \mathcal {L}} &=& \frac{\alpha ^2\Lambda ^2}{336G}\sqrt{\frac{5\pi }{6}} \bigg \lbrace -\sin (2\chi )\left[ { {\frac{\pi r}{R_*}}}j_1\left({ {\frac{\pi r}{R_*}}}\right) +11 j_2\left({ {\frac{\pi r}{R_*}}}\right) +12{ {\frac{\pi r}{R_*}}}\cot \left({ {\frac{\pi r}{R_*}}}\right)j_2\left({ {\frac{\pi r}{R_*}}}\right) \right]\left(Y_2^1-Y_2^{-1}\right)\nonumber \\ \,+\, \sin ^2\chi \left[ -2{ {\frac{\pi r}{R_*}}}j_1\left({ {\frac{\pi r}{R_*}}}\right) +20 j_2\left({ {\frac{\pi r}{R_*}}}\right) +18{ {\frac{\pi r}{R_*}}}\cot \left({ {\frac{\pi r}{R_*}}}\right)j_2\left({ {\frac{\pi r}{R_*}}}\right)\right]\left(Y_2^2+Y_2^{-2}\right)\nonumber \\ \,+\, \sqrt{6}\sin \left(2\chi \right)\left[-{ {\frac{\pi r}{R_*}}}j_1\left({ {\frac{\pi r}{R_*}}}\right)+3j_2\left({ {\frac{\pi r}{R_*}}}\right)+2{ {\frac{\pi r}{R_*}}}\cot \left({ {\frac{\pi r}{R_*}}}\right)j_2\right]\left(Y_4^1-Y_4^{-1}\right)\nonumber \\ \,-\, \sqrt{3}\sin ^2\chi \left[-{ {\frac{\pi r}{R_*}}}j_1\left({ {\frac{\pi r}{R_*}}}\right)+3j_2\left({ {\frac{\pi r}{R_*}}}\right)+2{ {\frac{\pi r}{R_*}}}\cot \left({ {\frac{\pi r}{R_*}}}\right)j_2\right]\left(Y_4^2+Y_4^{-2}\right)\bigg \rbrace . \end{eqnarray}

(A5)

APPENDIX B: NUMERICAL SOLUTION AND ERROR ESTIMATION

We use the NDSolve command in the software package mathematica in order to solve our system of equations. NDSolve is not capable of solving DAEs as boundary-value problems, as we would need, so we instead convert our original equations into a system of ODEs by working with equations (164) and (174). We implement the boundary conditions as summarized in Table 2. The full solution to our precession problem should be an infinite sum of harmonics, but we terminate at a value of l_max = 4, as we have failed to find a successful method of solving the equations for arbitrary truncation values; implementation of the techniques of Section 5.4 to l_max > 4 typically results in a solution that does not satisfy the highest l equation to acceptable accuracy. Finally, we non-dimensionalize the variables in our equations by suitable combinations of G, ρ_c and R_* (these three quantities contain instances of mass, length and time, so together they may be used to non-dimensionalize any other quantity). This allows us to work with quantities whose magnitude does not deviate excessively from unity, unlike the unwieldy physical parameters for stellar models. It also allows one to scale results easily to parameters appropriate for models of main-sequence stars, white dwarfs and neutron stars.

The problems discussed up to this point have been mathematical in origin, but the numerical solution has also required some care. We have found that the numerical algorithm typically requires a 25-digit precision or higher for success. Finally, mathematica returns zero-divided-by-zero errors with the equations in their original form. To avoid these, certain quantities need to be evaluated at a point slightly displaced from the centre. To minimize the introduced error, and to prevent it from affecting conditions at the outer boundary (the crust–core boundary), we define the function

\begin{equation} \mathcal {D}(r)\equiv {\mathfrak d} (r-R_{\rm out}), \end{equation}

(B1)

where |$\mathfrak {d}$| is some small constant, and in the equations in the mathematica notebook, we make the replacement |$r\mapsto r+\mathcal {D}(r)$| in the quantities |$rU^{\prime }_l(r),\Upsilon ^{\prime }_l(r),(r\rho ^{\prime }_0/\rho _0)^{\prime }$|⁠. At the end of this section, we verify that the introduced error is small and converges to zero as |$\mathfrak {d}\rightarrow 0$|⁠.

B1 Verification of DAE solution method

In Section 5.4, we presented a method to solve the system of DAEs (164) and (165). We differentiated the second of these DAEs in order to get a system of ODEs (164) and (174); however, the set of solutions to these ODEs is larger than the set of solutions to the original DAEs, and so a typical solution to the ODEs would differ from a solution to the DAEs by some ‘integration constant’ factor. The method of Section 5.4 was designed to pick out the correct integration constant, so let us now check that our solutions do indeed satisfy the original DAEs.

We substitute our numerical solution |$\mathcal {S}\equiv \lbrace X_l,\bar{U}_l\rbrace$| back into equations (164) and (165). To produce a normalized error estimate, we recall that the X_l functions are only non-zero for odd l, and the |$\bar{U}_l$| functions are non-zero for even l. Since we solve (165) (in differentiated form) for odd l and (164) for even l, this means that for a particular value of l, we can associate one equation with one particular function, e.g. equation (165) for l = 1 may be associated with X₁. We produce normalized errors |$\mathfrak {e}_l(r)$| using this observation:

\begin{equation} \mathfrak {e}_l(r)\equiv \left\lbrace \begin{array}{@{}l@{\quad }l@{}}\displaystyle {\frac{\mathcal {S}\,{\rm substituted\, into\, equation\, (164)\,}}{\max {(\bar{U}_l)}}\,{\rm (for\, even \,}l{\rm )}},\\ \displaystyle {\frac{\mathcal {S}\,{\rm substituted\, into\, equation\, (165)\,}}{\max {(X_l)}}\,{\rm (for\, odd\, }l{\rm )}}. \end{array}\right. \end{equation}

(B2)

In Fig. B1, we show the resulting errors, whose small size indicate that we have successfully obtained a solution to the DAE system. For comparison, we also perform a consistency check by substituting the solution back into the ODEs (164) and (174); since these were the equations actually solved by mathematica, the resulting error is very small.

$Comparison of errors $\mathfrak {e}_l$ (see equation (B2)) for different boundary conditions. We convert our original DAEs (164) and (165) into ODEs which mathematica is capable of solving by differentiating the latter equation, giving equation (174). The left-hand plots show the result of substituting the obtained solution back into (164) and (174); this is merely a consistency check of mathematica’s algorithm, and the resulting error is predictably small (regardless of which set of boundary conditions we use). We then solve (164) and (174) together with a naïve set of boundary conditions (Table B1) that do not contain information about the original DAEs. We substitute the solution obtained with these boundary conditions back into the original DAEs (164) and (165), and show the resulting error in the central plots. For the lower ones of these plots, corresponding to the equation that was differentiated, a large, constant error can be seen. This reflects the fact that a solution to (174) will typically differ by some integration constant from a solution of the original equation (165). To fix this constant at the correct value (of zero), we need to use the more careful boundary conditions of Table 2. With these, as shown in the right-hand plots, the error becomes small.$

Figure B1.

Comparison of errors |$\mathfrak {e}_l$| (see equation (B2)) for different boundary conditions. We convert our original DAEs (164) and (165) into ODEs which mathematica is capable of solving by differentiating the latter equation, giving equation (174). The left-hand plots show the result of substituting the obtained solution back into (164) and (174); this is merely a consistency check of mathematica’s algorithm, and the resulting error is predictably small (regardless of which set of boundary conditions we use). We then solve (164) and (174) together with a naïve set of boundary conditions (Table B1) that do not contain information about the original DAEs. We substitute the solution obtained with these boundary conditions back into the original DAEs (164) and (165), and show the resulting error in the central plots. For the lower ones of these plots, corresponding to the equation that was differentiated, a large, constant error can be seen. This reflects the fact that a solution to (174) will typically differ by some integration constant from a solution of the original equation (165). To fix this constant at the correct value (of zero), we need to use the more careful boundary conditions of Table 2. With these, as shown in the right-hand plots, the error becomes small.

Open in new tab Download slide

If we had not worried about fixing the integration constant correctly, we might instead have adopted a set of boundary conditions like those of Table B1, which are natural-looking ones to use to ensure that solutions are well behaved at the centre and do not give rise to a current sheet at the outer boundary. Fig. B1 clearly shows that the resulting solutions fail to satisfy the original DAE system.

Table B1.

Open in new tab

An example set of superficially sensible boundary conditions on the ODE system (164) and (174), which are none the less not acceptable, as they do not ensure that the resulting solution is also a solution of the original DAEs (164) and (165).

m	Centre	Crust–core boundary
	\|$\bar{U}^{\prime }_2(0) = 0$\|	\|$R_{\rm out}\bar{U}^{\prime }_2(R_{\rm out})+4\bar{U}_2(R_{\rm out}) = 0$\|
1	\|$\bar{U}_4^{\prime }(0) = 0$\|	\|$R_{\rm out}\bar{U}^{\prime }_4(R_{\rm out})+6\bar{U}_4(R_{\rm out}) = 0$\|
		X₁(R_out) = 0
		X₃(R_out) = 0
	\|$\bar{U}^{\prime }_2(0) = 0$\|	\|$R_{\rm out}\bar{U}^{\prime }_2(R_{\rm out})+4\bar{U}_2(R_{\rm out}) = 0$\|
2	\|$\bar{U}_4^{\prime }(0) = 0$\|	\|$R_{\rm out}\bar{U}^{\prime }_4(R_{\rm out})+6\bar{U}_4(R_{\rm out}) = 0$\|
		–
		X₃(R_out) = 0

m	Centre	Crust–core boundary
	\|$\bar{U}^{\prime }_2(0) = 0$\|	\|$R_{\rm out}\bar{U}^{\prime }_2(R_{\rm out})+4\bar{U}_2(R_{\rm out}) = 0$\|
1	\|$\bar{U}_4^{\prime }(0) = 0$\|	\|$R_{\rm out}\bar{U}^{\prime }_4(R_{\rm out})+6\bar{U}_4(R_{\rm out}) = 0$\|
		X₁(R_out) = 0
		X₃(R_out) = 0
	\|$\bar{U}^{\prime }_2(0) = 0$\|	\|$R_{\rm out}\bar{U}^{\prime }_2(R_{\rm out})+4\bar{U}_2(R_{\rm out}) = 0$\|
2	\|$\bar{U}_4^{\prime }(0) = 0$\|	\|$R_{\rm out}\bar{U}^{\prime }_4(R_{\rm out})+6\bar{U}_4(R_{\rm out}) = 0$\|
		–
		X₃(R_out) = 0

Table B1.

Open in new tab

An example set of superficially sensible boundary conditions on the ODE system (164) and (174), which are none the less not acceptable, as they do not ensure that the resulting solution is also a solution of the original DAEs (164) and (165).

m	Centre	Crust–core boundary
	\|$\bar{U}^{\prime }_2(0) = 0$\|	\|$R_{\rm out}\bar{U}^{\prime }_2(R_{\rm out})+4\bar{U}_2(R_{\rm out}) = 0$\|
1	\|$\bar{U}_4^{\prime }(0) = 0$\|	\|$R_{\rm out}\bar{U}^{\prime }_4(R_{\rm out})+6\bar{U}_4(R_{\rm out}) = 0$\|
		X₁(R_out) = 0
		X₃(R_out) = 0
	\|$\bar{U}^{\prime }_2(0) = 0$\|	\|$R_{\rm out}\bar{U}^{\prime }_2(R_{\rm out})+4\bar{U}_2(R_{\rm out}) = 0$\|
2	\|$\bar{U}_4^{\prime }(0) = 0$\|	\|$R_{\rm out}\bar{U}^{\prime }_4(R_{\rm out})+6\bar{U}_4(R_{\rm out}) = 0$\|
		–
		X₃(R_out) = 0

m	Centre	Crust–core boundary
	\|$\bar{U}^{\prime }_2(0) = 0$\|	\|$R_{\rm out}\bar{U}^{\prime }_2(R_{\rm out})+4\bar{U}_2(R_{\rm out}) = 0$\|
1	\|$\bar{U}_4^{\prime }(0) = 0$\|	\|$R_{\rm out}\bar{U}^{\prime }_4(R_{\rm out})+6\bar{U}_4(R_{\rm out}) = 0$\|
		X₁(R_out) = 0
		X₃(R_out) = 0
	\|$\bar{U}^{\prime }_2(0) = 0$\|	\|$R_{\rm out}\bar{U}^{\prime }_2(R_{\rm out})+4\bar{U}_2(R_{\rm out}) = 0$\|
2	\|$\bar{U}_4^{\prime }(0) = 0$\|	\|$R_{\rm out}\bar{U}^{\prime }_4(R_{\rm out})+6\bar{U}_4(R_{\rm out}) = 0$\|
		–
		X₃(R_out) = 0

B2 Convergence

To prevent numerical errors at the centre, we evaluated some quantities at a small radial distance – given by equation (B1) – from their correct location. Our final check is to make sure that the introduced error is small, and converges to zero in the limit |$\mathfrak {d}\rightarrow 0$|⁠; we do this in Fig. B2. Here we have used the correct boundary conditions, from Table 2; note that the large ‘integration constant’ error resulting from using the boundary conditions of Table B1 does not converge away as |$\mathfrak {d}\rightarrow 0$|⁠.

$Convergence of the error in solutions as $\mathfrak {d}\rightarrow 0$. For $\mathfrak {d} = 10^{-5}$, the noisy inherent error present from mathematica’s solution becomes visible; i.e. the error introduced from evaluating quantities slightly away from the origin reduces to the level of intrinsic error from mathematica’s solution method. Note that the actual physical results for $\bar{U}_l,X_l$ for these three values of $\mathfrak {d}$ are indistinguishable.$

Figure B2.

Convergence of the error in solutions as |$\mathfrak {d}\rightarrow 0$|⁠. For |$\mathfrak {d} = 10^{-5}$|⁠, the noisy inherent error present from mathematica’s solution becomes visible; i.e. the error introduced from evaluating quantities slightly away from the origin reduces to the level of intrinsic error from mathematica’s solution method. Note that the actual physical results for |$\bar{U}_l,X_l$| for these three values of |$\mathfrak {d}$| are indistinguishable.

Open in new tab Download slide

Download all slides

m	Centre	Crust–core boundary
	\|$U^{\prime }_2(0) = 0$\|	\|$R_{\rm out}U^{\prime }_2(R_{\rm out})+4U_2(R_{\rm out}) = 0$\|
1	\|$U_4(0) = -\frac{55}{6\sqrt{6}}U_2(0)$\|	\|$R_{\rm out}U^{\prime }_4(R_{\rm out})+6U_4(R_{\rm out}) = 0$\|
	\|$X_3(0) = -\frac{7}{12}\sqrt{\frac{35}{2}}U_2(0)$\|	X₁(R_out) = 0
	\|$U^{\prime }_2(0) = 0$\|	\|$R_{\rm out}U^{\prime }_2(R_{\rm out})+4U_2(R_{\rm out}) = 0$\|
2	\|$U^{\prime }_4(0) = 0$\|	\|$R_{\rm out}U^{\prime }_4(R_{\rm out})+6U_4(R_{\rm out}) = 0$\|
	\|$X_3(0) = -\frac{2}{15\sqrt{7}}U_2(0)+\frac{207}{550}\sqrt{\frac{3}{7}}U_4(0)$\|	–

m	Centre	Crust–core boundary
	\|$U^{\prime }_2(0) = 0$\|	\|$R_{\rm out}U^{\prime }_2(R_{\rm out})+4U_2(R_{\rm out}) = 0$\|
1	\|$U_4(0) = -\frac{55}{6\sqrt{6}}U_2(0)$\|	\|$R_{\rm out}U^{\prime }_4(R_{\rm out})+6U_4(R_{\rm out}) = 0$\|
	\|$X_3(0) = -\frac{7}{12}\sqrt{\frac{35}{2}}U_2(0)$\|	X₁(R_out) = 0
	\|$U^{\prime }_2(0) = 0$\|	\|$R_{\rm out}U^{\prime }_2(R_{\rm out})+4U_2(R_{\rm out}) = 0$\|
2	\|$U^{\prime }_4(0) = 0$\|	\|$R_{\rm out}U^{\prime }_4(R_{\rm out})+6U_4(R_{\rm out}) = 0$\|
	\|$X_3(0) = -\frac{2}{15\sqrt{7}}U_2(0)+\frac{207}{550}\sqrt{\frac{3}{7}}U_4(0)$\|	–

m	Centre	Crust–core boundary
	\|$U^{\prime }_2(0) = 0$\|	\|$R_{\rm out}U^{\prime }_2(R_{\rm out})+4U_2(R_{\rm out}) = 0$\|
1	\|$U_4(0) = -\frac{55}{6\sqrt{6}}U_2(0)$\|	\|$R_{\rm out}U^{\prime }_4(R_{\rm out})+6U_4(R_{\rm out}) = 0$\|
	\|$X_3(0) = -\frac{7}{12}\sqrt{\frac{35}{2}}U_2(0)$\|	X₁(R_out) = 0
	\|$U^{\prime }_2(0) = 0$\|	\|$R_{\rm out}U^{\prime }_2(R_{\rm out})+4U_2(R_{\rm out}) = 0$\|
2	\|$U^{\prime }_4(0) = 0$\|	\|$R_{\rm out}U^{\prime }_4(R_{\rm out})+6U_4(R_{\rm out}) = 0$\|
	\|$X_3(0) = -\frac{2}{15\sqrt{7}}U_2(0)+\frac{207}{550}\sqrt{\frac{3}{7}}U_4(0)$\|	–

m	Centre	Crust–core boundary
	\|$U^{\prime }_2(0) = 0$\|	\|$R_{\rm out}U^{\prime }_2(R_{\rm out})+4U_2(R_{\rm out}) = 0$\|
1	\|$U_4(0) = -\frac{55}{6\sqrt{6}}U_2(0)$\|	\|$R_{\rm out}U^{\prime }_4(R_{\rm out})+6U_4(R_{\rm out}) = 0$\|
	\|$X_3(0) = -\frac{7}{12}\sqrt{\frac{35}{2}}U_2(0)$\|	X₁(R_out) = 0
	\|$U^{\prime }_2(0) = 0$\|	\|$R_{\rm out}U^{\prime }_2(R_{\rm out})+4U_2(R_{\rm out}) = 0$\|
2	\|$U^{\prime }_4(0) = 0$\|	\|$R_{\rm out}U^{\prime }_4(R_{\rm out})+6U_4(R_{\rm out}) = 0$\|
	\|$X_3(0) = -\frac{2}{15\sqrt{7}}U_2(0)+\frac{207}{550}\sqrt{\frac{3}{7}}U_4(0)$\|	–

m	Centre	Crust–core boundary
	\|$\bar{U}^{\prime }_2(0) = 0$\|	\|$R_{\rm out}\bar{U}^{\prime }_2(R_{\rm out})+4\bar{U}_2(R_{\rm out}) = 0$\|
1	\|$\bar{U}_4^{\prime }(0) = 0$\|	\|$R_{\rm out}\bar{U}^{\prime }_4(R_{\rm out})+6\bar{U}_4(R_{\rm out}) = 0$\|
		X₁(R_out) = 0
		X₃(R_out) = 0
	\|$\bar{U}^{\prime }_2(0) = 0$\|	\|$R_{\rm out}\bar{U}^{\prime }_2(R_{\rm out})+4\bar{U}_2(R_{\rm out}) = 0$\|
2	\|$\bar{U}_4^{\prime }(0) = 0$\|	\|$R_{\rm out}\bar{U}^{\prime }_4(R_{\rm out})+6\bar{U}_4(R_{\rm out}) = 0$\|
		–
		X₃(R_out) = 0

Month:	Total Views:
February 2017	4
March 2017	4
April 2017	14
May 2017	2
June 2017	18
July 2017	12
August 2017	3
September 2017	2
October 2017	4
November 2017	3
December 2017	1
January 2018	6
February 2018	1
March 2018	4
May 2018	1
June 2018	1
July 2018	2
August 2018	2
September 2018	3
December 2018	3
February 2019	5
March 2019	10
May 2019	1
June 2019	3
August 2019	8
September 2019	1
October 2019	8
November 2019	8
December 2019	5
January 2020	1
February 2020	5
April 2020	3
May 2020	4
June 2020	11
July 2020	5
August 2020	6
September 2020	4
October 2020	8
November 2020	3
December 2020	5

Article Contents

Non-rigid precession of magnetic stars Free

Abstract

1 INTRODUCTION

2 GENERATION OF PRECESSIONAL MOTION IN A MAGNETIzED FLUID STAR

2.1 The bulk motion of the star

2.2 Deviation from rigid-body precession due to internal fluid dynamics

2.3 The need for second-order perturbations

2.4 Key parameters for precessing magnetic stars

3 FLUID PRECESSION AS A SECOND-ORDER PERTURBATION PROBLEM

3.1 Full equations

3.2 Perturbative scheme

3.3 Zeroth-order equations

3.4 Order-εB equations

3.5 Order-εα equations

3.6 Order-εαεB equations

3.7 Equations of order |$\epsilon _\alpha ^2$|⁠, |$\epsilon _B^2$| and higher

4 SOLUTION OF THE ZEROTH- AND FIRST-ORDER EQUATIONS

4.1 The zeroth-order equations

4.2 The order-εB equations

4.3 The order-εα equations

4.3.1 The centrifugal force

4.3.2 The |$\mathcal {O}(\epsilon _\alpha )$| system as a single Helmholtz equation

4.3.3 Homogeneous solution

4.3.4 Particular solution

4.3.5 Full interior solution

4.3.6 Exterior solution and boundary conditions

4.3.7 The first-order perturbation in the density

5 SOLUTION OF THE SECOND-ORDER EQUATIONS

5.1 Preliminaries on solenoidal fields

5.2 An explicit toroidal–poloidal split of the curled Euler equation

5.3 Equations for |$\delta {\boldsymbol B}$|

5.3.1 Analysis of equations

5.4 Differential algebraic equations; conditions at the centre

5.5 Exterior solution

5.6 Internal solution behaviour as outer boundary → R*

5.7 Conditions in the crust of a neutron star

5.7.1 Importance of elastic force terms

5.7.2 Exceeding the crustal breaking strain

5.7.3 Crust–core coupling

5.8 Boundary conditions for our model

6 THE |$\xi$|-MOTIONS

7 RESULTS

7.1 The poloidal and toroidal radial functions

7.2 Behaviour of higher multipoles at the centre

7.3 Magnetic-field perturbations

7.4 |$\xi$|-motions

8 DISCUSSION

9 SUMMARY

Acknowledgments

REFERENCES

APPENDIX A: SPHERICAL-HARMONIC DECOMPOSITIONS OF SOURCE TERMS

APPENDIX B: NUMERICAL SOLUTION AND ERROR ESTIMATION

B1 Verification of DAE solution method

B2 Convergence

Citations

Views

Altmetric

Email alerts

Astrophysics Data System

Citing articles via

Latest

Most Read

Most Cited

This Feature Is Available To Subscribers Only

Non-rigid precession of magnetic stars

3.4 Order-ε_B equations

3.5 Order-ε_α equations

3.6 Order-ε_αε_B equations

4.2 The order-ε_B equations

4.3 The order-ε_α equations

5.6 Internal solution behaviour as outer boundary → R_*