Magnetic equilibria of relativistic axisymmetric stars: The impact of flow constants

At the full non-linear level, we work with a general stationary and axisymmetric metric $g_{\mu\nu}$ in spherical-like coordinates $\{t,r,\theta,\varphi\}$ with signature $\{-,+,+,+\}$. Later on it will be convenient to introduce a static and spherically symmetric ‘background’, where the magnetic field is treated perturbatively (see Sec. IV). Here and in Sec. III, MHD dynamics are treated non-linearly, though the fluid is assumed to be perfect, infinitely conducting, and barotropic.

The stress-energy tensor appearing in the Einstein equations contains both fluid and electromagnetic components. The fluid piece is written as

$$T^{\mu\nu}_{\rm fluid}=\left(\epsilon+p\right)u^{\mu}u^{\nu}+pg^{\mu\nu},$$

(1)

where $\epsilon,p,u^{\mu}$ denote the fluid energy density, pressure, and four-velocity, respectively. The former parameter can be expressed in terms of the rest-mass density and specific internal energy $\tilde{\epsilon}$ as

$$\epsilon=\rho(1+\tilde{\epsilon}).$$

(2)

With regard to these densities, there is a certain degree of notational confusion in the literature. In Refs. Ciolfi et al. (2009, 2010), for instance, $\rho$ is used to denote the energy density, which makes it appear at first glance that the resulting GS equation disagrees with others in the literature Ioka and Sasaki (2003, 2004), though this is not the case (cf. Sec. IV.5). The notation used here is the typical one when working with the Tolman-Oppenheimer-Volkoff (TOV) equations, for example. For a barotropic fluid, the energy density and rest-mass density are related through the first law of thermodynamics, $d(\epsilon/\rho)=-pd(1/\rho)$.

Expression (1) corresponds to the stress-energy of a single fluid, where we have $\rho=mn$ for baryon mass, $m$, and number density, $n$. In reality, neutron star matter comprises multiple particle species ($e^{-},p,n,\ldots$) and a multi-fluid description is needed (see, for instance, Refs. Prix (2005); Glampedakis et al. (2012)). The distinction between particle constituents is important when considering the dynamics of magnetised stars, as the four-current is entirely sourced by the flow of the electrons relative to the other charged particles in the limit of charge neutrality. The single fluid description is based on the additional assumptions of a negligible electron inertia and a common velocity for the rest of the particles (charged and uncharged).

The flow constants we introduce in the next section that arise from the induction equation should thus be strictly associated with the electron number density and not the global one. This point is discussed further in Sec. III.3. For a single fluid in the ideal MHD limit, the electromagnetic stress-energy tensor is given by

$$T^{\alpha\beta}_{\rm EM}=g^{\alpha\beta}\frac{B^{2}}{8\pi}+\frac{1}{4\pi}\left(B^{2}u^{\alpha}u^{\beta}-B^{\alpha}B^{\beta}\right),$$

(3)

where $B^{\mu}$ is the magnetic field four-vector. The electric and magnetic components are generally given as $E_{\mu}=F_{\mu\nu}u^{\nu}$ and $B_{\alpha}=\tfrac{1}{2}\sqrt{-g}\epsilon_{\alpha\beta\gamma\delta}u^{\beta}F^{\gamma\delta}$ for Levi-Civita symbol $\epsilon_{\alpha\beta\gamma\delta}$ and Faraday tensor $F_{\mu\nu}=A_{\nu,\mu}-A_{\mu,\nu}$. Here $A_{\mu}$ is the vector potential, whose existence is implied by Maxwell’s equations, $\nabla_{[\alpha}F_{\mu\nu]}=0$, and $E_{\mu}=0$ defines the ideal MHD condition of infinite conductivity.

A result of immediate interest concerns a theorem by Carter Carter (1969, 1973). While a strictly azimuthal flow $(u^{i}=0)$ implies that one can consider a circular spacetime, $g_{ti}=g_{\varphi i}=0$, inclusion of meridional flows or mixed poloidal-toroidal fields generally forces these extra off-diagonal terms to be non-zero and means one cannot work with the Papapetrou metric Papapetrou (1966). Intuitively, if one considers a radial flow ($u^{r}\neq 0$), for instance, it is clear that the physical dynamics are not symmetric under time reversal if the system is also rotating. Similarly, for a mixed magnetic field there will generally be a poloidal current and the same logic applies (Frieben and Rezzolla, 2012). Spacetimes with purely toroidal fields are still circular however, since one can then exploit the azimuthal coordinate/gauge freedom Gourgoulhon et al. (2011). Additionally, one can always set $g_{r\theta}=0$ regardless.

The necessity of non-circular components is one key reason why constructing mixed poloidal-toroidal equilibria at a non-linear level is so difficult in GR Bocquet et al. (1995); it is only recently that such equilibria have been numerically constructed Uryū et al. (2019) and evolved Tsokaros et al. (2022) self-consistently¹¹1It has been argued that even for mixed, virial-strength fields, the degree of non-circularity is of order $\lesssim 10^{-3}$ Oron (2002); Pili et al. (2017). Whether this remains true for all physical setups is unclear owing to the numerical nature of the statement, but for computational purposes a circular approximation may be reasonable. (see also Refs. Lovelace et al. (1986); Fujisawa et al. (2013); Ofengeim and Gusakov (2018) for some Newtonian solutions). Aside from mathematical complexity, these extra pieces enrich the physical system as discussed in Sec. VI.

In this section we list the formulae for the flow constants, as derived in the BO papers Bekenstein and Oron (1978, 1979) and the aforementioned references. A proper derivation is relegated to Appendix A owing to the length of the calculations involved. Each of the constants here are derived from the conservation law $\nabla_{\mu}T^{\mu\nu}=\nabla_{\mu}\left(T^{\mu\nu}_{\rm fluid}+T^{\mu\nu}_{\rm EM}\right)=0$ together with the Einstein equations $R_{\mu\nu}-\tfrac{1}{2}Rg_{\mu\nu}=8\pi T_{\mu\nu}$ of a general stationary and axisymmetric spacetime.

Although we use the term ‘flow’, Eq. (107) from Appendix A shows that for any conserved $Q$ (i.e. $u^{\alpha}Q_{,\alpha}=0$) one necessarily has that $B^{\alpha}Q_{,\alpha}=0$. Introducing the magnetic potential $\psi=A_{\varphi}$, which can be shown to satisfy $B^{\alpha}\psi_{,\alpha}=0$, we express the flow constants in a more familiar form using $\psi$. For circular spacetimes, it is always possible to introduce a zero-angular-momentum observer such that the poloidal field admits a Chandrasekhar-like description with $B^{r}\propto\psi_{,\theta}$ and $B^{\theta}\propto\psi_{,r}$ up to metric factors Bocquet et al. (1995). Although the decomposition is not as simple in general, we still have that $B^{i}\propto A_{\varphi}$ by definition and thus $B^{i}=\mathcal{O}(\psi)$ to leading order since $u^{i}A_{t,i}=0$ from the ideal MHD condition $u^{\nu}F_{t\nu}=0$ (see also Sec. 2.1.2 in C08).

Adopting the notation of IS03 and IS04, we have the following relations amongst the five flow constants ${\cal C}(\psi),\Omega(\psi),E(\psi),L(\psi),$ and $D(\psi)$:

$$E=-\left(\mu+\frac{B^{2}}{4\pi n}\right)u_{t}-\frac{{\cal C}}{4\pi}(u_{t}+\Omega u_{\varphi})B_{t},$$

(4)

$$L=\left(\mu+\frac{B^{2}}{4\pi n}\right)u_{\varphi}+\frac{{\cal C}}{4\pi}(u_{t}+\Omega u_{\varphi})B_{\varphi},$$

(5)

$$D=E-\Omega L,$$

(6)

where $\mu$ is the specific enthalpy (this is defined in Appendix A and should not be confused with the chemical potential unless the temperature is everywhere zero). Some of these constants can be used to express the magnetic field $B^{\mu}$ and meridional flow $u^{i}$ ($i=\{r,\theta\}$) via

$$B^{\mu}=-{\cal C}n\left[\,(u_{t}+\Omega u_{\varphi})u^{\mu}+\delta_{t}^{\mu}+\Omega\delta^{\mu}_{\varphi}\,\right],$$

(7)

and

$$B^{2}=-{\cal C}n(B_{t}+\Omega B_{\varphi}).$$

(8)

Using (8) in Eqs. (4) and (5) one finds the alternative expressions

	$\displaystyle E$	$\displaystyle=-\mu u_{t}+\frac{{\cal C}\Omega}{4\pi}(u_{t}B_{\varphi}-u_{\varphi}B_{t}),$		(9)
	$\displaystyle L$	$\displaystyle=\mu u_{\varphi}+\frac{{\cal C}}{4\pi}(u_{t}B_{\varphi}-u_{\varphi}B_{t}).$		(10)

It is illuminating to decompose (7) into components, viz.

	$\displaystyle B^{i}$	$\displaystyle=-{\cal C}n(u_{t}+\Omega u_{\varphi})u^{i},$		(11)
	$\displaystyle B^{\varphi}$	$\displaystyle=-{\cal C}n\left[\Omega+(u_{t}+\Omega u_{\varphi})u^{\varphi}\right],$		(12)
	$\displaystyle B^{t}$	$\displaystyle=-{\cal C}n\left[1+(u_{t}+\Omega u_{\varphi})u^{t}\right].$		(13)

Following IS04 and BO79, these constants can be interpreted as follows. $E(\psi)$ and $L(\psi)$ represent the total enthalpy of the plasma and specific enthalpic angular momentum of the magnetic field, respectively; $D(\psi)$ symbolises the specific plasma energy according to a corotating observer and $\Omega(\psi)$ relates to the angular velocity²²2Note that what is called $\Omega$ here generally differs from the actual angular velocity $u^{\varphi}/u^{t}$, with agreement only in the limit that the toroidal field vanishes; see Eq. (4.42) of Ref. Gourgoulhon et al. (2011), noting that what we call $\Omega$ is denoted there by $\omega$ and $-A$ by BO. Because of its hybrid nature, $\Omega$ could also be considered ‘induction-like’. of the fluid. These constants can be thought of as ‘Bernoulli-like’, in the sense that there is a simple but non-trivial hydrodynamic limit; the constant $D$ in particular can be directly identified with the Bernoulli integral (see Appendices). The remaining constant ${\cal C}(\psi)$ – representing the magnetic field strength relative to the magnitude of meridional flow – is instead ‘induction-like’. The existence of this constant is a direct consequence of the ideal assumption $E_{\mu}=0$ and becomes trivial in the non-magnetic limit. Moreover, expression (11) shows that in the limit of vanishing meridional flow, $u^{i}\rightarrow 0$, the poloidal field generally vanishes unless ${\cal C}\rightarrow\infty$.

Flow constants analogous to the ones found by BO78 in the context of GRMHD equilibria are known to exist in Newtonian axisymmetric-stationary systems. These are discussed, for instance, in Mestel’s textbook Mestel (1999) but their origins can be traced back to older work Prendergast (1956); Chandrasekhar (1956); Woltjer (1959). Here we summarise the expressions relating these constants with the various hydromagnetic parameters (adopting standard vectorial notation); a detailed derivation of these expressions can be found in Appendix B.

The defining property of a Newtonian flow constant $Q$ is $\mathbf{v}\cdot\bm{\nabla}Q=0$. It is easy to show (see Appendix B) that the same flow constant is a function $Q(\psi)$, where the magnetic scalar potential $\psi$ is defined via the poloidal field component

$$\mathbf{B}_{\rm p}=\bm{\nabla}\psi\times\bm{\nabla}\varphi.$$

(14)

As in GR, the Newtonian constants can be grouped according to their ‘equation of origin’, namely the ideal induction [Eq. (137)] or the Euler equation [Eq. (136)]. In the first group we find the flow constants ${\cal C}_{\rm N}$ and $\Omega_{\rm N}$ which, as their notation suggests, are analogous to the GR ones ${\cal C}$ and $\Omega$. These appear in the following relation between the magnetic field and fluid velocity

$$\mathbf{B}={\cal C}_{\rm N}\rho\left(\mathbf{v}-\varpi\Omega_{\rm N}\bm{\hat{\varphi}}\right),$$

(15)

where $\varpi=r\sin\theta$. In particular, the poloidal/meridional component of this expression is

$$\mathbf{B}_{\rm p}={\cal C}_{\rm N}\rho\mathbf{v}_{\rm p}.$$

(16)

In the second group we find the following three Bernoulli-like flow constants which result by projecting the Euler equation along the $\mathbf{v}$ and $\mathbf{B}$ vectors

	$\displaystyle E_{\rm N}$	$\displaystyle=\frac{1}{2}v^{2}+\Phi+h-\frac{\Omega_{\rm N}{\cal C}_{\rm N}}{4\pi}\varpi B_{\varphi},$		(17)
	$\displaystyle L_{\rm N}$	$\displaystyle=\varpi\left(v_{\varphi}-\frac{{\cal C}_{\rm N}}{4\pi}B_{\varphi}\right),$		(18)
	$\displaystyle D_{\rm N}$	$\displaystyle=\frac{1}{2}v^{2}+\Phi+h-\varpi\Omega_{\rm N}v_{\varphi},$		(19)

where $h$ is the fluid enthalpy and $\Phi$ is the gravitational potential.

One of the most surprising consequences of the flow constant MHD formalism is the theorem derived in BO79 under a circular spacetime assumption (see their Sec. VI for more context): “the absence of meridional circulation implies the vanishing of the toriodal [sic] magnetic field.”

The validity of this theorem hinges on the simultaneous presence of the meridional flow and the flow constant ${\cal C}$ in the ‘induction-originated’ equation (7). Indeed, setting $u^{i}=0$ in (11), and in combination with a non-vanishing poloidal field $B^{i}\neq 0$, leads to the requirement ${\cal C}\to\infty$. This divergence would cause the hydrodynamical constants $E,L$ to diverge too, unless [see Eqs. (4), (5)]

	$\displaystyle B_{t}$	$\displaystyle=g_{tt}B^{t}+g_{t\varphi}B^{\varphi}+g_{ti}B^{i}=0,$		(20)
	$\displaystyle B_{\varphi}$	$\displaystyle=g_{\varphi\varphi}B^{\varphi}+g_{t\varphi}B^{t}+g_{\varphi i}B^{i}=0.$		(21)

The determinant of this linear system is non-vanishing in the absence of horizons, and as a consequence the unique solution is the trivial one, $B^{t}=B^{\varphi}=0$.

As an additional intuitive argument supporting the theorem’s conclusion, consider how the toroidal field appears to a locally inertial observer. The relevant tetrad component reads $B_{(\varphi)}\propto L(\psi)/{\cal C}(\psi)$ up to metric factors (e.g. Ioka and Sasaki (2004); Ciolfi et al. (2009); see Sec. V). Since ${\cal C}\to\infty$ in the $u^{i}\to 0$ limit, keeping the above non-zero requires that $L\to\infty$ at an equally fast rate if the observer is to witness a toroidal field. This is clearly incompatible with the hydrodynamic limit, and thus we are forced to conclude that the toroidal field must vanish to maintain consistency.

Therefore a stationary-axisymmetric MHD equilibrium without meridional flow must be purely poloidal as well as purely spacelike. As far as astrophysically relevant MHD equilibria are concerned this is not an acceptable situation as it is well known that purely poloidal fields are generically unstable (e.g. Lasky et al. (2011); Ciolfi and Rezzolla (2012)). This means that a stable mixed poloidal-toroidal configuration should be accompanied by some amount of meridional flow.

As an aside we note that a vanishing meridional flow is still compatible with a purely toroidal field (which is also generically unstable Frieben and Rezzolla (2012); Kiuchi et al. (2011)). Under these circumstances, setting $u^{i}=B^{i}=0$ in (11) does not imply a divergent flow constant ${\cal C}$. Moreover, the vanishing poloidal field implies $F_{r\varphi}=F_{tr}=0$ and again $\Omega\neq u^{\varphi}/u^{t}$ (see Appendix A). The two non-vanishing field components then read $B^{t}={\cal C}n(u^{\varphi}-\Omega u^{t})u_{\varphi}$ and $B^{\varphi}=-{\cal C}n(u^{\varphi}-\Omega u^{t})u_{t}$; these expressions can be compared with those found in Ref. Kiuchi and Yoshida (2008), which studied nonlinear toroidal equilibria.

Without presenting details, which are carefully laid out by IS03, the Euler equation can be shown to lead to the GS equation³³3The following conversion rule should be applied when adapting formulae from IS to our notation: $\mu_{\rm IS}=\mu/m$, $Q_{\rm IS}=Q/m$ where $Q=\{{\cal C},L,E,D\}$.

	$\displaystyle 0=$	$\displaystyle J^{\varphi}-\Omega J^{t}+\frac{1}{{\cal C}\sqrt{-g}}\left[\,(\mu u_{r})_{,\theta}-(\mu u_{\theta})_{,r}\,\right]$		(22)
		$\displaystyle+nu^{\varphi}\left(L^{\prime}-\Lambda{\cal C}^{\prime}\right)-nu^{t}\left[\,E^{\prime}-\Lambda({\cal C}\Omega)^{\prime}\,\right],$

where

$$\Lambda=\frac{1}{4\pi}(u_{t}B_{\varphi}-u_{\varphi}B_{t})$$

(23)

and $J^{\mu}=\nabla_{\nu}F^{\mu\nu}/4\pi$ is the four-current. Many studies of equilibria concern themselves only with Eq. (22), discarding other information under the assumption that the other ideal MHD equations are trivial. This is not really the case, however, as the number of free functions in (22) make conceptually obvious. As pointed out in Sec. II.2, the flow constant ${\cal C}$ arises directly from the induction equation and thus key physical points are missed if one considers (22) in isolation (see also Sec. III.3). One such point was described in Sec. III.1, namely that setting the meridional flow to zero but keeping a mixed poloidal-toroidal field is inconsistent, even though such a choice leaves (22) well behaved. Another concerns the nature of Alfvén points (i.e. where the magnetic energy density matches the kinetic energy density of the flow) and the stellar surface, which we discuss in Sec. III.4.

Although we work within a relativistic framework in this paper, many of the issues we discuss carry over at a Newtonian level. This is exemplified by the presence of BO-like flow constants in Newtonian systems (Sec. II.3). Owing to the complexity of multifluid dynamics in GR, especially when convective motions and mixed magnetic fields are permitted, we opt to describe some physical features related to induction here in a Newtonian language. In particular, the fact that ‘induction-like’ flow constants are only associated with the electron fluid while the others are associated with the bulk applies to both frameworks. This bares on the physical interpretation of various limits, as discussed throughout, some of which have been addressed in detail by Prix Prix (2005).

Our objective in this section, therefore, is to provide a careful rederivation of the induction-originated ${\cal C}_{\rm N}$ and $\Omega_{\rm N}$ flow constants by accounting for the distinct proton and electron fluid velocities. These are denoted, respectively, as $\mathbf{v}_{\rm c}$ and $\mathbf{v}_{\rm e}$. The induction equation itself can be derived from Faraday’s law, $\bm{\nabla}\times\mathbf{E}=-\partial_{t}\mathbf{B}/c$, and the Euler equation for the electron fluid with the inertial terms omitted, i.e.

$$\mathbf{E}=-\frac{1}{c}\mathbf{v}_{\rm e}\times\mathbf{B}+\bm{\nabla}(\mu_{\rm e}+m_{\rm e}\Phi).$$

(24)

The resulting equation is

$$\partial_{t}\mathbf{B}=\bm{\nabla}\times(\mathbf{v}_{\rm e}\times\mathbf{B}).$$

(25)

The manipulations of the ‘single-fluid’ induction equation described in Appendix B can be exactly repeated if we change $\mathbf{v}\to\mathbf{v}_{\rm e}$ and $\rho\to\rho_{\rm e}$. In terms of the poloidal/meridional and toroidal/azimuthal components of the vectors $\mathbf{B},\mathbf{v}$,

$$\mathbf{B}=\mathbf{B}_{\rm p}+B_{\varphi}\bm{\hat{\varphi}},\qquad\mathbf{v}_{\rm e}=\mathbf{v}_{\rm ep}+v_{e\varphi}\bm{\hat{\varphi}},$$

(26)

we find

$$\mathbf{B}_{\rm p}={\cal C}_{\rm N}\rho_{\rm e}\mathbf{v}_{\rm ep},$$

(27)

where ${\cal C}_{\rm N}$ is constant along poloidal field lines, i.e.

$$\mathbf{B}_{\rm p}\cdot\bm{\nabla}{\cal C}_{\rm N}=0.$$

(28)

In axisymmetry this is equivalent to

$$\mathbf{v}_{\rm e}\cdot\bm{\nabla}{\cal C}_{\rm N}=0,$$

(29)

which shows that ${\cal C}_{\rm N}$ is also conserved with respect to the electron flow.

The poloidal field can be written in a divergence-free form in terms of the scalar potential $\psi(r,\theta)$,

$$\mathbf{B}_{\rm p}=\bm{\nabla}\psi\times\bm{\nabla}\varphi.$$

(30)

In this parametrisation the poloidal field lines lie in constant $\psi$ surfaces and therefore the same must be true for the flow lines of $\mathbf{v}_{\rm ep}$. As a consequence, ${\cal C}_{\rm N}={\cal C}_{\rm N}(\psi)$.

Further manipulation of Eq. (25) leads to the flow constant $\Omega_{\rm N}(\psi)$ as part of an expression for the toroidal field

$$B_{\varphi}={\cal C}_{\rm N}\rho_{\rm e}\left(v_{\varphi}^{\rm e}-\varpi\Omega_{\rm N}\right).$$

(31)

In combination with (27), this allows us to write the Newtonian counterpart of Eq. (7) as

$$\mathbf{B}={\cal C}_{\rm N}\rho_{\rm e}\left(\mathbf{v}_{\rm e}-\varpi\Omega_{\rm N}\bm{\hat{\varphi}}\right).$$

(32)

The final step consists in expressing $\mathbf{v}_{\rm e}$ in terms of the total current $\mathbf{J}$ and the proton velocity (which in non-superfluid matter can be taken to coincide with the neutron fluid velocity). Making use of the assumed charge neutrality of the system, $n_{\rm e}=n_{\rm c}$, we have

$$\mathbf{v}_{\rm e}=\mathbf{v}_{\rm c}-\mathbf{J}/n_{\rm e},$$

(33)

which then implies that

$$\mathbf{B}={\cal C}_{\rm N}\rho_{\rm e}\left(\mathbf{v}_{\rm c}-\varpi\Omega_{\rm N}\bm{\hat{\varphi}}\right)-{\cal C}_{\rm N}m_{\rm e}\mathbf{J}.$$

(34)

Assuming $n_{\rm e}=0$ at the surface, this expression can easily accommodate a non-vanishing magnetic field as long as $\mathbf{J}\neq 0$, thus alleviating any pathological behaviour at the surface (see related discussion in the following section). In practice, this may not be an issue since astrophysical neutron stars are endowed with a non-vacuum magnetosphere in which $n_{\rm e}$ is non-vanishing Goldreich and Julian (1969). Eq. (34) may also invalidate the BO result applying in the limit $\mathbf{v}_{\rm cp}\to 0$ (see the final paragraph of Appendix B), as ${\cal C}_{{\rm N}}\to\infty$ is no longer necessary to maintain $\mathbf{B}_{p}\neq 0$ if the poloidal current $\mathbf{J}_{\rm p}$ does not vanish.

As a final point, although entropy gradients have been ignored here, one generally has $\mathbf{v}\cdot\nabla s=0$ in the single fluid case, meaning that $s=s(\psi)$. Yoshida et al. (2012) argue that this forces a vanishing meridional flow else convectively unstable regions with $\nabla p\cdot\nabla s<0$ exist in a realistic, stratified star. This is clearly in conflict with the BO conclusion, though the above considerations, which naturally carry over to GR, may remedy the situation.

In the context of twisted-torus models where one considers the GS equation in isolation, the vanishing density at the stellar surface generally poses no problem to regularity. In fact, in the exterior of a static star where $n=0$ and there is a purely poloidal ‘test’ field, Eq. (22) decouples with respect to the angular coordinates and admits an analytical solution in terms of hypergeometric functions Konno et al. (1999). Although rotating stars will not be surrounded by vacuum in reality Goldreich and Julian (1969), the vacuum GS equation necessarily describes a force-free magnetosphere (see Ref. Pétri (2016) for a discussion of ‘force-free’ in the GR context).

When considering non-zero velocity fields, the limit $n\to 0$ needs special attention as a result of equations (7) and (8). The two subequations which are most relevant here are written again for convenience, viz.

	$\displaystyle B^{i}$	$\displaystyle=-{\cal C}n\left(u_{t}+\Omega u_{\varphi}\right)u^{i},$		(35)
	$\displaystyle B^{\varphi}$	$\displaystyle=-{\cal C}n\left[\Omega+(u_{t}+\Omega u_{\varphi})u^{\varphi}\right].$		(36)

The toroidal field can be assumed to be confined in the stellar interior, rendering Eq. (36) trivial at the surface and beyond, without placing any restriction on ${\cal C}$. The poloidal equations are clearly problematic though at $n=0$, unless ${\cal C}\to\infty$ for every $\psi$ curve that crosses the surface or one carefully treats the multi-fluid dynamics, as in Sec. III.3. This and a related problem can be described in terms of the poloidal Alfvén Mach number⁴⁴4It is somewhat unintuitive that $M_{\rm A}$ does not depend explicitly on the magnetic field, though this is because the poloidal flow scales with the poloidal magnetic field, as evident from expression (35).,

$$M^{2}_{\rm A}=-\frac{4\pi\mu}{{\cal C}^{2}n\left(g_{tt}+2\Omega g_{t\phi}+\Omega^{2}g_{\phi\phi}\right)}.$$

(37)

From Eqs. (5) and (8), it can be shown that $u_{\varphi}$ tends to diverge as $M_{\rm A}\to 1$ (i.e. at the Alfvén surface) unless $L$ and $\Omega$ are proportional to each other there (Ioka and Sasaki, 2004). This problem is discussed in the Newtonian context by Mestel Mestel (1968) and Ogilvie Ogilvie (2016) (see also section VI C in Ref. Gourgoulhon et al. (2011)). Unfortunately, this does not really help reduce the pool of flow constants since proportionality is only required in a particular limit. The more obvious problem is that for a finite value of $M_{\rm A}$ to be maintained one requires ${\cal C}^{2}n>0$ everywhere, which would seem to demand ${\cal C}\sim\psi^{-1}$ with a field confined to the interior such that $\psi\to 0$ towards the surface (see also Sec. IV). Although this choice was the one adopted in IS04, it is clearly unrealistic.

There are several possibilities worth considering in order to avoid this conundrum:

• (i)

  The ideal MHD approximation (upon which the BO relations are based) breaks down in the limit $n\to 0$. This is somewhat expected as, in reality, it is meaningless to discuss a current or velocity gradient without matter (e.g., manipulations of the ideal MHD assumption $E_{\mu}=F_{\mu\nu}u^{\nu}=0$ become dubious; see Appendix A).

• (ii)

  A non-relaxed crust alters (or invalidates) some of the BO relations. That is, the term $T^{\mu\nu}_{\rm shear}$, defined through the Carter-Quintana equations Carter and Quintana (1972), may also enter into the total stress-energy tensor. This term will naturally adjust the boundary conditions and flow constant structure (see Ref. Kojima et al. (2022) for some results in the Newtonian context). The presence of an ocean further complicates the dynamics.

• (iii)

  The function ${\cal C}(\psi)$ behaves in such a way that it tends to infinity everywhere except in an equatorial torus. Suppose that $\psi$ attains a maximum $\psi_{c}$ at the stellar surface, and we write ${\cal C}={\tilde{\mathcal{C}}}/\left[\psi(\psi-\psi_{c})\Theta(\psi-\psi_{c})\right]$ for constant ${\tilde{\mathcal{C}}}$ with the meridional flow confined to the region $\psi>\psi_{c}$. Outside of this torus (e.g. at the surface) we have ${\cal C}\to\infty$ but $u^{i}\to 0$ and so $B^{i}$ from Eq. (35) can be non-zero even if $n\to 0$ there, while inside both $u^{i}$ and ${\cal C}$ are finite with $n>0$ and so $B^{i}$ is well behaved there also. Although there are mathematical issues with this approach related to the undefined limit of the Heaviside function as $\psi\to\psi_{c}$, one has to bear in mind that the use of discontinuous functions is a convenient way of modelling sharp yet continuous gradients of the real physical system.

• (iv)

  The neutron star is unlikely to be surrounded by vacuum in reality as mentioned previously, arguably rendering this issue moot. Care must be taken in this case to choose appropriate boundary conditions, for instance using the Goldreich-Julian number density $n_{\rm GJ}\sim B_{\mu}u^{\mu}$ Goldreich and Julian (1969). If $n\to 0$ anywhere in the exterior universe however, this problem reappears at that interface. Note this is related to point (i), as the ideal MHD approximation itself may not be valid in a tenuous magnetosphere.

The following sections consider the leading-order behaviour of the flow constants, essentially through a power-counting analysis. A similar procedure was undertaken by IS04, however they assumed that rotation is always sub-leading and that $\psi\to 0$ at the stellar surface (cf. Sec. III.4). One goal of this section is to reexamine their results when these conditions are relaxed.

Firstly, it is instructive to define perturbative when considering an expansion in terms of magnetic parameters. If, as is typical in numerical GR investigations, we further consider ‘units’ of length such that the star has unit mass, $M_{\star}=1$, the basic measure of magnetic flux becomes the dimensionless number $[\text{G cm}^{2}]=1.39\times 10^{-30}(1.4M_{\odot}/\tilde{M})$ for physical-units mass $\tilde{M}$. Given that the magnetic flux is related to the local field strength through $\psi\approx BR_{\star}^{2}$, where $R_{\star}$ denotes the stellar radius, we see that expansions in powers of $\psi$ are well behaved ($\psi\ll 1$) provided that $(B/\text{G})(R_{\star}/\text{cm})^{2}\times 1.39\times 10^{-30}(1.4M_{\odot}/\tilde{M})\ll 1$, which requires $B\ll 7.2\times 10^{17}\,$G for $R_{\star}=10^{6}\,$cm and $\tilde{M}=1.4M_{\odot}$. This same limit can be independently derived by demanding that the magnetic pressure at the center of the star is much less than the minimum central hydrostatic pressure. This limit is physically equivalent to a magnetic-to-gravitational-binding energy ratio much below unity.

Although we consider expansions in $\psi$, it is not obvious that this is the most natural approach when considering a rotating star. In reality, most if not all neutron stars are both ‘slow’ (rotating much slower than the break-up limit) and ‘weakly magnetised’ (as defined above), and it may be more appropriate to consider a double expansion in both $u^{\varphi}$ and $\psi$ simultaneously in the form of a magnetic-Hartle-Thorne scheme (see Ref. Steil (2018) for such an approach). We ignore such complications here (though see Footnote 5 below).

Here we treat the magnetic field as a perturbation on a static star, with the magnetic stream function $\psi$ used as a perturbative bookkeeping parameter. We assume the following scaling for the magnetic field components

$$B^{i},B^{\varphi}={\cal O}(\psi),\quad B^{t}={\cal O}(\psi^{\gamma}),$$

(39)

for some $\gamma\geq 1$. To be more precise, the scaling $B^{i}={\cal O}(\psi)$ is a direct consequence of the definition of $\psi=A_{\varphi}$ (see Sec. II.2) and we impose a comparable toroidal field to maintain linearity (cf. Sec. IV.5). This ‘3+1’ scaling accommodates the presence of a mixed poloidal-toroidal field to leading order, while also taking into account the post-Newtonian character of $B^{t}$.

Metric perturbations appear first at ${\cal O}(\psi^{2})$ since $T^{\alpha\beta}_{\rm EM}\sim B^{2}$ Ioka and Sasaki (2003), i.e.

$$g_{\mu\nu}=g_{\mu\nu}^{\rm TOV}+{\cal O}(\psi^{2}).$$

(40)

As discussed earlier, a magnetic field may induce fluid motion in an otherwise static configuration. We suppose that azimuthal flow appears at some sub-leading order,

$$u^{\varphi}={\cal O}(\psi^{\eta}),\qquad\Omega={\cal O}(\psi^{\eta}),$$

(41)

for $\eta\geq 1$, together with $u^{t}={\cal O}(1)$ set by the wind equation $u^{t}u_{t}=-1$ at leading order. The meridional flow scaling is taken to be

$$u^{i}={\cal O}(\psi^{\beta}),$$

(42)

for some $\beta\geq\eta$. The flow constants $E,L$, and $D$ are hydrodynamical. That is, they are present even as $\psi\to 0$. We Taylor-expand these through

	$\displaystyle E$	$\displaystyle=E_{0}+E_{1}\psi+E_{2}\psi^{2}+{\cal O}(\psi^{3}),$		(43)
	$\displaystyle L$	$\displaystyle=L_{0}+L_{1}\psi+L_{2}\psi^{2}+{\cal O}(\psi^{3}),$		(44)

where $E_{0}$ and $L_{0}$ are non-magnetic parameters, and the expansion for $D$ follows from Eq. (6). Last but not least, the flow constant ${\cal C}$ is assumed to scale as

$${\cal C}={\cal O}(\psi^{\alpha}),$$

(45)

where $\alpha\neq 0$, since induction-like constants should not enter at background order. The above scalings are not independent; inserting them in Eqs. (11), (12), (13) and (8), we obtain

$$\psi\sim nu_{t}\psi^{\alpha+\beta},$$

(46)

$$\psi\sim nu_{t}u^{\varphi}\psi^{\alpha}\sim n\,\psi^{\alpha+\eta},$$

(47)

$$\psi^{\gamma}\sim n(\Omega u^{t}-u^{\varphi})\psi^{\alpha+\eta},$$

(48)

and

$$\psi^{2}\sim n(B_{t}+\Omega B_{\varphi})\psi^{\alpha}\sim n\psi^{\alpha+k},\,k=\mbox{min}[\gamma,\eta+1],$$

(49)

respectively, where the symbol $\sim$ is used to relate parameters with a non-trivial $\psi$-scaling. The balance of $\psi$-powers in the first two expressions leads to

$$\alpha=1-\eta,\qquad\beta=\eta.$$

(50)

Subsequently, Eqs. (48) and (49) lead to

$$\gamma\geq 1+\eta,\qquad\Omega u^{t}-u^{\varphi}={\cal O}(\psi^{\gamma-1})\lesssim{\cal O}(\psi^{\eta}).$$

(51)

The resulting scalings ${\cal C}\sim\psi^{-1}$ and $u^{i},u^{\varphi}\sim\psi^{2}$ are consistent with the analysis of IS04 if we set $\eta=2$. More generally, however, any $\eta>1$ with $\alpha=1-\eta<0$ is viable, though demanding integer solutions to permit (Laurent/Taylor) expansions forces $\eta=2$ and $\alpha=-1$. The remaining ‘hydrodynamical’ Eqs. (9), (10) are compatible with them and offer no new information about the scalings. It is also easy to verify that

$$g_{tt}^{\rm TOV}(u^{t})^{2}=-1+{\cal O}(\psi^{2}).$$

(52)

One last remark concerns the non-magnetic ($\psi\to 0$) limit of the MHD equations. The balance of $\psi$-powers implemented here means that all of them should automatically behave smoothly in that limit (for example, the above purely magnetic equations reduce to a trivial $0=0$).

For a static star, the above analysis thus determines the leading-order behaviour of the flow constants irrespective of the true system’s initial conditions.

Much of the previous section’s results depend on the crucial assumption of a non-rotating star in the limit $\psi\to 0$. More astrophysically relevant is the scenario where the star is rigidly rotating in the non-magnetic limit and $g_{t\varphi}\neq 0$ (though cf. Footnote 5 below). In such a case,

$$u^{\varphi}={\cal O}(1),\qquad\Omega={\cal O}(1).$$

(53)

As rigid rotation may preclude comparable energy partitions, we instead allow for a more general scaling of the toroidal field through

$$B^{\varphi}={\cal O}(\psi^{\lambda}),$$

(54)

together with⁵⁵5Note that if $A_{t}\sim u^{i}$ and $u^{\varphi}=\mathcal{O}(1)$, the arguments outlined at the beginning of Sec. II.2 would appear to imply $B^{i}={\cal O}(\psi^{\chi})$ with $\chi=\text{min}[\beta,1]$. This, however, forces $\lambda=0$ if $\beta\leq 1$ which leads to a contradiction. Physically speaking though, this simply echoes the sentiment that it does not really make sense to talk about a non-magnetic star with circulating charged particles (i.e. the ‘background’ does not truly exist). On the other hand, if the fluid is composed only of uncharged matter, then the flow constant ${\cal C}$ does not exist; cf. Sec. III.3. $B^{i}={\cal O}(\psi)$. As before, we assume

$$u^{i}={\cal O}(\psi^{\beta}),\qquad{\cal C}={\cal O}(\psi^{\alpha}).$$

(55)

Upon inserting these in Eqs. (11), (12), and (13) we find

	$\displaystyle\psi$	$\displaystyle\sim nu_{t}\psi^{\alpha+\beta},$		(56)
	$\displaystyle\psi^{\lambda}$	$\displaystyle\sim nu_{t}u^{\varphi}\psi^{\alpha}\sim n\psi^{\alpha},$		(57)
	$\displaystyle\psi^{\gamma}$	$\displaystyle\sim n(\Omega u^{t}-u^{\varphi})\psi^{\alpha},$		(58)

from which we easily deduce that

$$\alpha=\lambda,\qquad\alpha+\beta=1,\qquad\gamma=\alpha.$$

(59)

If we have comparable poloidal and toroidal field strengths, then we find that $\alpha=1,~{}\beta=0$. That is,

$$u^{i}={\cal O}(1),\qquad{\cal C}={\cal O}(\psi),$$

(60)

for $\lambda=1$. These scalings could be preempted from Eq. (50) with $\eta=0$, though are markedly different to those obtained for a non-rotating star in Sec. IV.2. In particular, a leading-order meridional flow is an uncomfortable result as it implies a strong deviation from rigid rotation in the non-magnetic limit. This impasse prompts us to explore non-standard scalings for the toroidal field component, in combination with a sub-leading and non-integer-scaling meridional flow, $\beta>0$. This requirement gives $\lambda=1-\beta<1$, which entails a toroidally-dominated MHD equilibrium unless $\psi$ itself is non-perturbative (i.e. $\psi^{\lambda}>\psi$ for $0<\psi,\lambda<1$). Some discussion on this point is given in Sec. VI.

Here we consider a perturbative and integer expansion of the GS equation (22) without $\mathcal{O}(1)$ rotation. Based on the previously obtained scalings and writing

$${\cal C}={{\tilde{\mathcal{C}}}}/{\psi},$$

(61)

where ${\tilde{\mathcal{C}}}$ is a constant parameter, we have

$$J^{\varphi}\sim{\cal O}(\psi),\,\,\Omega J^{t}\lesssim{\cal O}(\psi^{2}),\,\,\Lambda=\frac{1}{4\pi}u_{t}B_{\varphi}+{\cal O}(\psi^{5}),$$

(62)

$$\frac{1}{{\cal C}\sqrt{-g}}\left[\,(\mu u_{r})_{,\theta}-(\mu u_{\theta})_{,r}\,\right]\sim\frac{\psi^{\beta-\alpha}}{\tilde{{\cal C}}}\sim{\cal O}(\psi^{3}),$$

(63)

	$\displaystyle nu^{\varphi}\left(L^{\prime}-\Lambda{\cal C}^{\prime}\right)$	$\displaystyle=nu^{\varphi}\left(L_{1}+2L_{2}\psi-\alpha\Lambda\tilde{{\cal C}}\psi^{\alpha-1}\right)$		(64)
		$\displaystyle=n\tilde{{\cal C}}u^{\varphi}\frac{\Lambda}{\psi^{2}}+{\cal O}(\psi^{2}),$

and

	$\displaystyle-nu^{t}\left[E^{\prime}-\Lambda({\cal C}\Omega)^{\prime}\right]$	$\displaystyle=-nu^{t}\Big{[}E_{1}+2E_{2}\psi$		(65)
		$\displaystyle-\frac{\Lambda\tilde{{\cal C}}}{\psi}\left(\Omega^{\prime}-\Omega\psi^{-1}\right)\Big{]}+{\cal O}(\psi^{2}).$

The above imply that (22) takes the form

$$\mathcal{O}(\psi^{2})=J^{\varphi}+n\tilde{{\cal C}}\frac{u^{\varphi}\Lambda}{\psi^{2}}-nu^{t}\left[E_{1}+2E_{2}\psi-\Omega_{2}\Lambda\tilde{{\cal C}}\right],$$

(66)

with $\Omega(\psi)=\Omega_{2}\psi^{2}$. Upon inserting the approximate $\Lambda$ from Eq. (62), one finally obtains the $\mathcal{O}(\psi)$ GS equation

	$\displaystyle 0=$	$\displaystyle J^{\varphi}+\frac{n}{4\pi}\Big{\{}\tilde{{\cal C}}\frac{u^{\varphi}u_{t}B_{\varphi}}{\psi^{2}}$		(67)
		$\displaystyle-u^{t}\left[4\pi\left(E_{1}+2E_{2}\psi\right)-\Omega_{2}u_{t}B_{\varphi}\tilde{{\cal C}}\right]\Big{\}}.$

By using expression (10), one can also rewrite this equation using the flow constant $L$ (modulo some factors) in exchange for $u_{t}B_{\varphi}$; see Sec. V.

The first interesting aspect of Eq. (67) we point out is that the term $\sim nu^{t}E_{1}$ is ${\cal O}(1)$, and therefore we should expect $E_{1}=0$. That is to say, a strict power-counting argument appears to preclude the possibility of $E_{1}\neq 0$, else in the non-magnetic limit (67) unphysically implies that $0=nu^{t}$. The choice $E_{1}\neq 0$ is however made in the literature (where $E_{i+1}$ are usually rebranded as $c_{i}$), as discussed in the next section.

A second point concerns the explicit appearance of $\Omega$, which is absent in the linearised equation presented by IS04 [see their Eq. (64)]. We believe this is because these authors implicitly assume $\eta>2$ rather than $\eta=2$ as written in their Eq. (58), meaning that all $\Omega$ terms are ignored as being higher order. [Note also the typo in their Eq. (60) in the expansion of $D(\psi)$.]

In this section, we discuss the perturbative GS equation and compare various approaches towards it found in the literature. Ignoring $\Omega$ corrections at all orders, writing $E^{\prime}(\psi)\propto F(\psi)$ and $L_{0}/{\tilde{\mathcal{C}}}\propto\zeta$, one can eventually write the linear GS equation in coordinate form as

$$0=J_{\varphi}-{\zeta^{2}}e^{-2\nu}\psi-r^{2}\sin^{2}\theta(\epsilon+p)F(\psi),$$

(68)

with $F^{\prime\prime}(\psi)=0$. We write Eq. (68) in this way so as to match the $\mathcal{O}(\psi)$ equation presented by C08 [i.e. their Eq. (57)]. Sometimes the parameter $\zeta$ is upgraded to a general function of $\psi$, making (68) non-linear (see below).

Strict linearisation of the GS equation limits the possible set of equilibria matching to an exterior poloidal electrovacuum, with the only options being (i) the entire field is confined within the star, (ii) $B^{\varphi}$ is discontinuous and there is a non-zero surface current, or (iii) $\zeta=0$ and $B^{\varphi}=0$ identically. To see this, note that if one matches to an exterior where $\psi\neq 0$ then $\zeta^{2}\psi\neq 0$ at the surface unless $\zeta=0$ everywhere or $\zeta(R_{\star})=0$ imposed by hand discontinuously. Option (i) is that selected by IS04 and others (e.g. Haskell et al. (2008)), while (ii) is chosen by C08.

While neither of these options characterise a realistic system, they are the only choices available for a mixed poloidal-toroidal field at strictly linear order. An alternative is to consider non-linear choices for $\zeta^{2}\psi$ and/or $F(\psi)$, allowing a toroidal field that decays regularly towards the boundary (à la twisted torus), such as in Refs. Ciolfi et al. (2009, 2010); Ciolfi and Rezzolla (2013). However, if one were to take $\zeta\propto\psi$ (for instance), quantities explicitly related to $u^{i}$ should appear in the GS equation since meridional flow enters at $\mathcal{O}(\psi^{3})$. Omitting these terms thus implies that power counting is not handled self-consistently, presenting a conceptual problem. At quadratic order one could avoid the appearance of meridional terms (i.e. taking $\zeta\propto\sqrt{\psi}$), though this again likely rules out a regular toroidal component. Moreover, metric perturbations should be accounted for as these appear at $\mathcal{O}(\psi^{2})$, though this is often circumvented with the Cowling approximation.

A related issue surrounds $F(\psi)$; C08 and others (e.g. Ciolfi et al. (2009, 2010); Glampedakis et al. (2012); Ciolfi and Rezzolla (2013)) Taylor-expand this function as

$$F(\psi)=c_{0}+c_{1}\psi+{\cal O}(\psi^{2}),$$

(69)

though declare $c_{0}={\cal O}(B)$. In C08, the authors study the simpler system (68) with $c_{1}=0$, which admits decoupled multipoles. In C09, by contrast, both $c_{0}$ and $c_{1}$ are kept and the multipole expansion becomes non-trivial (cf. Section V). In either case, however, it seems inconsistent to keep $c_{0}\neq 0$ since this would imply one does not recover $0=0$ in the non-magnetic limit, $\psi\rightarrow 0$.

The above two points show that there are subtle issues involved in the construction of magnetic equilibria, aside from those related to induction explored previously. In the next section we resolve the second of these points, meaning that we numerically build equilibria where $c_{1}\neq 0$ but $c_{0}=0$, doing so in the flow-constant language. Solving the non-linear system self-consistently with a twisted-torus configuration, meridional flow, and metric backreactions is left to future work.