# Multipolar electromagnetic fields around neutron stars: general-relativistic vacuum solutions

Jérôme Pétri

► To cite this version:

Jérôme Pétri. Multipolar electromagnetic fields around neutron stars: general-relativistic vacuum solutions. Monthly Notices of the Royal Astronomical Society, 2017, 472 (3), pp.3304-3336. 10.1093/mn-ras/stx2147. hal-01645477

HAL Id: hal-01645477

https://hal.science/hal-01645477v1

Submitted on 29 Aug 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Multipolar electromagnetic fields around neutron stars: general-relativistic vacuum solutions

J. Pétri*

Observatoire astronomique de Strasbourg, Université de Strasbourg, CNRS, UMR 7550, F-67000 Strasbourg, France

Accepted 2017 August 16. Received 2017 August 16; in original form 2016 October 24

## ABSTRACT

Magnetic fields inside and around neutron stars are at the heart of pulsar magnetospheric activity. Strong magnetic fields are responsible for quantum effects, an essential ingredient to produce leptonic pairs and the subsequent broad-band radiation. The variety of electromagnetic field topologies could lead to the observed diversity of neutron star classes. Thus, it is important to include multipolar components to a presumably dominant dipolar magnetic field. Exact analytical solutions for these multipoles in Newtonian gravity have been computed in recent literature. However, flat space-time is not adequate to describe physics in the immediate surroundings of neutron stars. We generalize the multipole expressions to the strong gravity regime by using a slowly rotating metric approximation such as the one expected around neutron stars. Approximate formulae for the electromagnetic field including frame dragging are computed from which we estimate the Poynting flux and the braking index. Corrections to leading order in compactness and spin parameter are presented. As far as spin-down luminosity is concerned, it is shown that frame dragging remains irrelevant. For high-order multipoles starting from the quadrupole, the electric part can radiate more efficiently than the magnetic part. Both analytical and numerical tools are employed.

**Key words:** magnetic fields – methods: analytical – methods: numerical – stars: general – stars: neutron – pulsars: general.

## 1 INTRODUCTION

Our understanding of pulsar physics has benefitted from recent advances in multiwavelength observational campaigns as well as from developments of new numerical tools able to investigate carefully its magnetosphere and the subsequent radiation mechanisms from a theoretical point of view. As the quality and quantity of data increases regularly, theoreticians are forced to improve the physics of their models in conjunction with the precision of their predictions. Detailed analysis of pulsar phase-resolved spectra and polarization properties requires accurate models dealing with all possible perturbations departing from a too simple description of pulsar magnetospheres. General relativity belongs to one of these additional physical ingredients compulsory to investigate properly neutron star electrodynamics. This fact became clear in the past years. Indeed efficient pair creation in force-free magnetospheres seems to require frame-dragging effects close to the magnetic poles (Philippov et al. 2015). Pétri (2016a) investigated in depth general-relativistic force-free magnetospheres. In the same vain Ruiz, Paschalidis & Shapiro (2014) looked at the spin-down luminosity and attempted to match neutron star exterior described in the force-free regime to its interior described by relativistic magnetohydrodynamic, see also Paschalidis & Shapiro (2013). Konno & Kojima (2000) investigated the impact of general-relativistic corrections to curvature radiation. Oscillations of neutron stars in general relativity were also of interest (Kojima & Hosonuma 2000). Morozova, Ahmedov & Zanotti (2010) studied the influence of neutron star oscillations in general relativity on the plasma density in the magnetosphere for a aligned rotator. Curvature of space-time and frame-dragging effects on the surrounding electromagnetic field was already emphasized by Beskin (1990) and Muslimov & Tsygan (1992).

Multipolar components in the vicinity of the neutron star surface are useful to introduce small-scale structures in the magnetic topology and therefore enhance the curvature of magnetic field lines. This finding could boost the pair creation scenario by increasing the multiplicity factor as required from observations of the Crab pulsar and nebula for instance. The long-term evolution of multipolar fields still remains controversial. According to Gourgouliatos & Cumming (2014), Hall effect and Ohmic dissipation tend to decay faster higher multipoles.

<footnote>* E-mail: jerome.petri@astro.unistra.fr</footnote>
Small-scale magnetic structures stability in neutron star crust evolving according to Hall drift has been investigated by Marchant et al. (2014). It can lead to stable and unstable configurations down to the Ohmic dissipation scale. On the other hand, magnetic multipoles could be hidden inside the neutron star interior leaving only a dipole outside until their re-emergence from the stellar crust after a fallback disc episode lasting for several thousands of years (Igoshev, Elfritz & Popov 2016). Multipolar fields confined in the crust certainly exist and must become visible outside neutron stars. They are able to explain several important aspects of their underlying physical engine. It is therefore compulsory to add such topologies to any neutron star magnetosphere.

Another important observation is that any rotating field can be expanded into multipolar components. Thus, Bonazzola, Mottez & Heyvaerts (2015) and Pétri (2015) showed how to compute exact analytical solutions to multipoles in closed analytical form for flat space-time. Extension to general relativity is highly desirable. The simplest case is a rotating dipole for which Kojima, Matsunaga & Okita (2004) gave an approximate solution in general relativity comparing their results with previous analytical work of Rezzolla, Ahmedov & Miller (2001). An analytical estimate for the dipole spin-down has been given by Rezzolla & Ahmedov (2004). Recently, these authors extended their analysis of oscillations by adding damping due to heating and Joule dissipation (Rezzolla & Ahmedov 2016). In this paper, we show how to extend those results to any multipole and to any order of accuracy.

Maxwell equations in presence of strong gravity using the 3+1 formalism are used to solve for arbitrary electromagnetic field configurations in vacuum as presented in Section 2. Exact solutions for static multipolar magnetic fields in Schwarzschild background metric in terms of hypergeometric functions are reminded in Section 3. In the same gravitational field, exact analytical solutions are found in terms of local confluent Heun functions as explained in Section 4. Frame dragging is included in an approximate fashion by numerically solving the system of elliptic equations (Helmholtz equations) for the unknown fields as exposed in Section 5. Our analytical treatment is supported by time-dependent numerical simulations of Maxwell equations in general relativity as presented in Section 6. Conclusions are drawn in Section 7.

## 2 MAXWELL EQUATIONS AND THEIR SOLUTIONS

In this section, the general formalism to solve Maxwell equations in general-relativistic vacuum is reviewed. For the 3+1 split of space-time, we follow the conventions and definitions given by Pétri (2013).

### 2.1 Split of space-time metric

We split space-time into a 3+1 foliation such that the background metric $g_{ik}$ is expressed as

$$ds^2 = g_{ik} dx^i dx^k = \alpha^2 c^2 dt^2 - \gamma_{ab} (dx^a + \beta^a c dt) (dx^b + \beta^b c dt) \quad (1)$$

where $x^i = (ct, x^a)$, $t$ is the time coordinate or universal time and $x^a$ some associated space coordinates. The metric signature is given by $(+, -, -, -)$. $\alpha$ is the lapse function, $\beta^a$ the shift vector and $\gamma_{ab}$ the spatial metric of absolute space. By convention, Latin letters from $a$ to $h$ are used for the components of vectors in absolute space, in the range $\{1, 2, 3\}$, whereas Latin letters starting from $i$ are used for four-dimensional vectors and tensors, in the range $\{0, 1, 2, 3\}$. A fiducial observer (FIDO) is defined by its four-velocity $n^i$ such that

$$n^i = \frac{dx^i}{d\tau} = \frac{c}{\alpha} (1, -\beta) \quad (2a)$$

$$n_i = (\alpha c, 0). \quad (2b)$$

This vector is orthogonal to the hypersurface of constant time coordinate $\Sigma_t$. Its proper time $\tau$ is measured according to $d\tau = \alpha dt$. The relation between the determinants of the space-time metric $g$ and the pure spatial metric $\gamma$ is given by $\sqrt{-g} = \alpha \sqrt{\gamma}$. For a slowly rotating neutron star, the lapse function is

$$\alpha = \sqrt{1 - \frac{R_s}{r}} \quad (3)$$

and the shift vector

$$c \beta = -\omega r \sin \vartheta e_{\varphi} \quad (4a)$$

$$\omega = \frac{R_s a c}{r^3}. \quad (4b)$$

$R_s$ is the Schwarzschild radius of the neutron star.

We use a spherical coordinate system $(r, \vartheta, \varphi)$ and an orthonormal spatial basis $(e_r, e_\vartheta, e_\varphi)$. The metric of a slowly rotating neutron star remains close to the usual flat space, except for the radial direction. Indeed, the components of the spatial metric are given in Boyer-Lindquist coordinates by

$$\gamma_{ab} = \begin{pmatrix} \alpha^{-2} & 0 & 0 \\ 0 & r^2 & 0 \\ 0 & 0 & r^2 \sin^2 \vartheta \end{pmatrix}. \quad (5)$$
For this slow-rotation approximation, the spatial metric does not depend on the spin frequency of the massive body but only on its mass *M* through *α*. This justifies the slow-rotation approximation. The spin parameter *a* is related to the angular momentum *J* by *J* = *Ma* *c*. It follows that *a* has units of a length and should satisfy *a* ≤ *R*<sub>s</sub>/2. Introducing the moment of inertia *I*, we also have *J* = *I*Ω. In the special case of a homogeneous and uniform neutron star interior with spherical symmetry, the moment of inertia reads

$$I = \frac{2}{5} M R^2. \tag{6}$$

Thus, the spin parameter can be expressed as

$$\frac{a}{R_s} = \frac{2}{5} \frac{R}{R_s} \frac{R}{r_L}. \tag{7}$$

$r_L$ is the light-cylinder radius.

For the remainder of this paper, we will use this expression to relate the spin parameter intervening in the metric to the spin frequency of the neutron star. From the above expression, note that the parameter *a*/*R*<sub>s</sub> remains usually small.

## 2.2 Maxwell equations in general relativity

Maxwell equations in 3+1 notation take a traditional form close to the one known in flat space-time except that, as the reader should keep in mind, the three-dimensional space is curved and differential operators defined according to the spatial metric $\gamma_{ab}$. In vacuum, the system reads

$$\nabla \cdot \mathbf{B} = 0 \tag{8a}$$

$$\nabla \times \mathbf{E} = -\frac{1}{\sqrt{\gamma}} \partial_t(\sqrt{\gamma} \mathbf{B}) \tag{8b}$$

$$\nabla \cdot \mathbf{D} = 0 \tag{8c}$$

$$\nabla \times \mathbf{H} = \frac{1}{\sqrt{\gamma}} \partial_t(\sqrt{\gamma} \mathbf{D}). \tag{8d}$$

The three-dimensional vector fields are not independent, they are related by two important constitutive relations, namely

$$\varepsilon_0 \mathbf{E} = \alpha \mathbf{D} + \varepsilon_0 c \beta \times \mathbf{B} \tag{9a}$$

$$\mu_0 \mathbf{H} = \alpha \mathbf{B} - \frac{\beta \times \mathbf{D}}{\varepsilon_0 c} \tag{9b}$$

$\varepsilon_0$ is the vacuum permittivity and $\mu_0$ the vacuum permeability. The curvature of absolute space is taken into account by the lapse function factor $\alpha$ in the first term on the right-hand side and the frame-dragging effect is included in the second term, the cross-product between the shift vector $\beta$ and the fields. Space curvature and frame dragging have therefore an imprint on the constitutive relations equation (9). From the auxiliary fields ($\mathbf{E}$, $\mathbf{H}$), we get the Poynting flux through a sphere of radius *r* by computing the two-dimensional integral on this sphere by

$$L = \int_{\omega} E \wedge H r^2 d\omega \tag{10}$$

where $d\omega$ is the infinitesimal solid angle and $\omega$ the full sky angle of 4 $\pi$ sr. This integral can be computed analytically in the asymptotic regime of large distances given in the wave zone. The Poynting flux should then be interpreted as the power radiated as seen by a distant observer and not the intrinsic spin-down as measured on the neutron star surface. Indeed, due to gravitational redshift, energy is degraded when photons propagate from the surface to the observer and this affects also the measured power in the wave zone. Our spin-down luminosity is computed according to the normalization done by this distant observer. Actually our distant observer is located at infinity, thus very faraway from the light cylinder and from any perceptible gravitational effect. In this paper, by distant observer we mean someone making measurements at spatial infinity.

## 2.3 General solution to Maxwell equations

It is formally possible to give arbitrary solutions to Maxwell equations for divergencelessness electric and magnetic fields in vacuum in a stationary regime. Indeed their expansion into vector spherical harmonics reads, neglecting a possible monopolar $\ell = 0$ contribution,

$$\mathbf{D}(r, \vartheta, \varphi, t) = \sum_{\ell=1}^{\infty} \sum_{m=-\ell}^{\ell} (\nabla \times [f_{\ell,m}^{\text{D}}(r) \mathbf{\Phi}_{\ell,m}] + g_{\ell,m}^{\text{D}}(r) \mathbf{\Phi}_{\ell,m}) e^{-i m \Omega t} \tag{11a}$$

$$\mathbf{B}(r, \vartheta, \varphi, t) = \sum_{\ell=1}^{\infty} \sum_{m=-\ell}^{\ell} (\nabla \times [f_{\ell,m}^{\text{B}}(r) \mathbf{\Phi}_{\ell,m}] + g_{\ell,m}^{\text{B}}(r) \mathbf{\Phi}_{\ell,m}) e^{-i m \Omega t} \tag{11b}$$
which correspond to stationary solutions expressed in the frame of a distant observer. The $\Phi_{\ell,m}$ are vector spherical harmonics defined and introduced in recent works by Pétri (2013). The functions $g_{\ell,m}^D$ and $g_{\ell,m}^B$ are related to the function $f_{\ell,m}^B$ and $f_{\ell,m}^D$ according to Maxwell equations by a linear scaling. Indeed, there exists a simple algebraic relation between $g_{\ell,m}^D$ and $f_{\ell,m}^B$ on one side, and between $g_{\ell,m}^B$ and $f_{\ell,m}^D$ on the other side such that

$$ \alpha g_{\ell,m}^D = +i \epsilon_0 m \tilde{\omega} f_{\ell,m}^B \quad (11c) $$

$$ \alpha g_{\ell,m}^B = -i \mu_0 m \tilde{\omega} f_{\ell,m}^D \quad (11d) $$

where $\tilde{\omega} = \Omega - \omega$ is the rotation rate as measured by a local observer and corrected for frame dragging. After substitution in Maxwell equations, we get inhomogeneous Helmholtz equations for the potentials $f_{\ell,m}^D$ and $f_{\ell,m}^B$. Indeed as shown by Pétri (2013), introducing the Helmholtz operator in curved space-time by

$$ W_{\ell,m}[f] = \alpha^2 \left[ \frac{1}{r} \frac{\partial}{\partial r} \left( \alpha^2 \frac{\partial}{\partial r} (r f) \right) - \frac{\ell (\ell+1)}{r^2} f \right] + m^2 \frac{\tilde{\omega}^2}{c^2} f \quad (12) $$

the potentials $f_{\ell,m}^D$ must satisfy

$$ W_{\ell}[f_{\ell,m}^{D}] = 3 \epsilon_0 \alpha^2 \frac{\omega}{r} \left[ f_{\ell-1,m}^{B} \sqrt{(\ell-1)(\ell+1)} J_{\ell,m} - f_{\ell+1,m}^{B} \sqrt{\ell(\ell+2)} J_{\ell+1,m} \right] \quad (13a) $$

where $J_{\ell,m} = \sqrt{\frac{\ell^2 - m^2}{4 \ell^2 - 1}}$ and similarly for the magnetic field $f_{\ell,m}^B$

$$ W_{\ell}[f_{\ell,m}^{B}] = -3 \mu_{0} \alpha^{2} \frac{\omega}{r} \left[ f_{\ell-1,m}^{D} \sqrt{(\ell-1)(\ell+1)} J_{\ell,m} - f_{\ell+1,m}^{D} \sqrt{\ell(\ell+2)} J_{\ell+1,m} \right]. \quad (13b) $$

The boundary conditions on the neutron star surface are imposed similarly to their Newtonian counterparts. We impose continuity of the radial component of the magnetic field **B** and continuity of the tangential component of the electric field **D**. Let us remind that inside a perfectly conducting star, the rotation of the plasma induces an electric field **E** satisfying

$$ E + r \Omega \sin\vartheta e_{\varphi} \times B = 0. \quad (14) $$

Thus, the electric field as measured by a FIDO is given by

$$ D = -\epsilon_0 \frac{\tilde{\omega}}{\alpha} r \sin\vartheta e_{\varphi} \times B. \quad (15) $$

The electromagnetic field symbolized by (D, B) has to verify the jump conditions across the plasma/vacuum interface as in flat space-time. Explicitly, the normal component $B^\mu$ and the tangential components ($D^\mu, D^\nu$) have to be continuous across the stellar surface. By construction, the component $D^\mu$ remains zero in the exterior vacuum, as it is inside the star. Using the divergencelessness expansion of (D, B) in equation (11), the coefficients for the electric field $f_{\ell,m}^D$ are related to those of the magnetic field $f_{\ell,m}^B$ by

$$ \alpha^2 \partial_r (r f_{\ell,m}^D) = \epsilon_0 r \tilde{\omega} \left[ \sqrt{(\ell+1)(\ell-1)} J_{\ell,m} f_{\ell-1,m}^B - \sqrt{\ell(\ell+2)} J_{\ell+1,m} f_{\ell+1,m}^B \right]. \quad (16) $$

The detailed derivation of these boundary conditions using recurrence properties of scalar spherical harmonics can be found in Pétri (2013). This last expression shows that the electric field strength is proportional to $\tilde{\omega}$ which contains the frame-dragging effect. This leads to a lowering of the actual rotation rate of the star as seen by a local observer on the surface. Thus, frame dragging decreases the electric field intensity and therefore also the spin-down luminosity corrections due to rotation of space-time. However, as shown later in this work, for realistic neutron star parameters, these corrections remain negligible. A second correction is induced by the space curvature and implies an additional factor $\alpha^{-2}$ compared to flat space-time. The constants of integration are magnified but this effect is sometimes completely cancelled by the general-relativistic spherical Hankel functions when $\ell = m$ as shown in the numerical results in Section 5. For other multipoles with $\ell > m$, compensation is only partial. These conclusions are discussed in detail in the numerical approximate solution Section 5.

The Helmholtz operator is conveniently rewritten by introducing the tortoise coordinate $r_*$ such that

$$ r_* = r + R_s \ln \left( \frac{r}{R_s} - 1 \right). \quad (17) $$

We are looking for solutions describing outgoing waves that reduce to $e^{i k_m r}$ in flat space-time, thus we introduce another unknown field $u_{\ell,m}^{B/D}$ such that

$$ f_{\ell,m}^{B/D}(r) = u_{\ell,m}^{B/D}(r) \frac{e^{i k_m r_*}}{r} \quad (18) $$

where $k_m = m/r_L$. The curved space-time Helmholtz operator in terms of these new dependent variables $u_{\ell,m}$ becomes

$$ W_{\ell,m}[f_{\ell,m}] = \left\{ \frac{1}{r^3} \left(1 - \frac{R_s}{r}\right) \left[ r(r - R_s) u_{\ell,m}'' + (2i k_m r^2 + R_s) u_{\ell,m}' - \ell(\ell+1)u_{\ell,m} \right] + k_m^2 \left[ \left(1 - \frac{\omega}{\Omega}\right)^2 - 1 \right] \frac{u_{\ell,m}}{r} \right\} e^{i k_m r_*} \quad (19) $$

$'$ and $''$ are the first and second derivatives respectively.
In terms of $u_{\ell,m}^{B/D}$, the inhomogeneous Helmholtz equations become for the electric field

$$ \frac{1}{r} \left[ r (r - R_s) u_{\ell,m}^{D''} + (2i k_m r^2 + R_s) u_{\ell,m}^{D'} - \ell (\ell+1) u_{\ell,m}^{D} \right] + k_m^2 \left[ \left(1 - \frac{\omega}{\Omega}\right)^2 - 1 \right] r \frac{u_{\ell,m}^{D}}{\alpha^2} = 3 \epsilon_0 \omega \left[ u_{\ell-1,m}^{B} \sqrt{(\ell-1)(\ell+1)} J_{\ell,m} - u_{\ell+1,m}^{B} \sqrt{\ell(\ell+2)} J_{\ell+1,m} \right] \quad (20a) $$

and similarly for the magnetic field

$$ \frac{1}{r} \left[ r (r - R_s) u_{\ell,m}^{B''} + (2i k_m r^2 + R_s) u_{\ell,m}^{B'} - \ell (\ell+1) u_{\ell,m}^{B} \right] + k_m^2 \left[ \left(1 - \frac{\omega}{\Omega}\right)^2 - 1 \right] r \frac{u_{\ell,m}^{B}}{\alpha^2} = -3 \mu_0 \omega \left[ u_{\ell-1,m}^{D} \sqrt{(\ell-1)(\ell+1)} J_{\ell,m} - u_{\ell+1,m}^{D} \sqrt{\ell(\ell+2)} J_{\ell+1,m} \right]. \quad (20b) $$

For one single multipole with potential specified at the surface by $f_{\ell,m}^{B}(R)$, the boundary conditions for the electric field reduces to

$$ \alpha^2 \partial_r (u_{\ell-1,m}^D e^{i k_m r_s}) = -\epsilon_0 r \tilde{\omega} \sqrt{(\ell-1)(\ell+1)} J_{\ell,m} f_{\ell,m}^B \quad (21a) $$

$$ \alpha^2 \partial_r (u_{\ell+1,m}^D e^{i k_m r_s}) = +\epsilon_0 r \tilde{\omega} \sqrt{\ell(\ell+2)} J_{\ell+1,m} f_{\ell,m}^B \quad (21b) $$

to be evaluated at the surface for $r=R$.

## 2.4 Wave zone and Poynting flux

The Poynting flux of a rotating multipole can be most easily computed in the asymptotic flat space-time at very large distance. Because of energy conservation law, the flux leaving the star must reach infinity, there is no absorption layer in between. In the wave zone, the expressions (11) can be drastically reduced by the fact that the potential functions behave asymptotically as $f_{\ell,m}(r) \approx u_{\ell,m}^{\infty} e^{i k_m r}/r$ where $\lim_{r \to +\infty} u_{\ell,m}(r) = u_{\ell,m}^{\infty}$. Neglecting the axisymmetric mode decreasing much faster, like $r^{-(\ell+1)}$, the electromagnetic field becomes in the limit of large distances $r \gg r_L$

$$ B_w = \sum_{\ell \ge 1, m \ne 0} -i \frac{e^{i m (k r - \Omega t)}}{r} k_m \left( u_{\ell,m}^{B,\infty} \Psi_{\ell,m} + \mu_0 c u_{\ell,m}^{D,\infty} \Phi_{\ell,m} \right) \quad (22a) $$

$$ D_w = \sum_{\ell \ge 1, m \ne 0} -i \frac{e^{i m (k r - \Omega t)}}{r} k_m \left( u_{\ell,m}^{D,\infty} \Psi_{\ell,m} - \epsilon_0 c u_{\ell,m}^{B,\infty} \Phi_{\ell,m} \right) \quad (22b) $$

$$ = \epsilon_0 c \mathbf{B}_w \wedge \mathbf{n}. \quad (22c) $$

Equation (22c) shows that the solution behaves as a monochromatic plane wave propagating in the radial direction $n = e_r$ at frequency $\Omega$. The time-averaged Poynting flux is therefore

$$ S = \frac{D_w \wedge B_w^*}{2 \mu_0 \epsilon_0} \quad (23) $$

where $B_w^*$ is the complex conjugate of $B_w$. Integrating the radial component of the Poynting vector along the solid angle $\omega$, we get the power radiated, using the orthonormality of the vector spherical harmonics, such that

$$ L = \int_{\omega} S^2 r^2 d\omega = \frac{c}{2 \mu_0} \sum_{\ell \ge 1, m \ne 0} k_m^2 \left( |u_{\ell,m}^{B,\infty}|^2 + \mu_0^2 c^2 |u_{\ell,m}^{D,\infty}|^2 \right). \quad (24) $$

The spin-down luminosity $L$ is independent of the radius as it should from the energy conservation law. Equation (24) represents the most general expression for the magneto-multipole losses from an arbitrary multipole magnetic field in general relativity. As soon as the constants of integration $u_{\ell,m}^{D/B,\infty}$ are known, the full stationary electromagnetic field is determined and its subsequent properties such as spin-down and magnetic topology. Our main goal is to fix these constants either with some analytical argument or more accurately via numerical integration of the elliptic problems related to the Helmholtz equations in curved space-time. We next recall the exact analytical solutions of the static multipoles in general relativity and then look for the stationary rotating multipoles solutions found numerically.

## 3 EXACT STATIC MULTIPOLE SOLUTIONS

Finding explicit exact solutions to the rotating multipole problem is difficult because there is no known analytical solution to Helmholtz equations in general relativity in the Schwarzschild metric. Nevertheless, in the static limit of non-rotating neutron stars, it is possible to find exact close formula for the multipoles to any order.

From the expansion into vector spherical harmonics equations (11), each function $f_{\ell,m}^{D/B}$ has to satisfy a homogeneous second-order linear differential equation in Schwarzschild space-time such that

$$ \frac{\alpha^2}{r} \partial_r (\alpha^2 \partial_r (r f_{\ell,m})) - \alpha^2 \frac{\ell (\ell+1)}{r^2} f_{\ell,m} = 0. \quad (25) $$
Introducing the new unknown function $\phi_{\ell, m} = r f_{\ell, m}$, we get the simple differential equation

$$ \partial_r(\alpha^2 \partial_r(\phi_{\ell,m})) - \frac{\ell(\ell+1)}{r^2} \phi_{\ell,m} = 0 \qquad (26) $$

to be solved with appropriate boundary conditions, namely vanishing potentials at spatial infinity. Moreover, introducing the normalized inverse radial coordinate by $x = R_s/r$, the functions $\phi_{\ell, m}$ will be solution of

$$ x^2 (1-x) \phi''_{\ell,m} + x (2-3x) \phi'_{\ell,m} - \ell(\ell+1) \phi_{\ell,m} = 0 \qquad (27) $$

which reduces to the hypergeometric differential equation by the change of unknown function $\phi_{\ell, m} = x^\ell v_{\ell, m}$ to

$$ x (1-x) v''_{\ell,m} + (2(\ell+1) - (2\ell+3)x) v'_{\ell,m} - \ell(\ell+2) v_{\ell,m} = 0. \qquad (28) $$

Setting the parameters $a = \ell$, $b = \ell + 2$ and $c = 2(\ell + 1)$, we indeed find the standard form of the hypergeometric differential equation (Olver 2010). The only solution vanishing at infinity is

$$ \phi_{\ell,m}(x) = C x^{\ell/2} F_1(\ell, \ell+2, 2(\ell+1), x) \qquad (29) $$

which is the expression found by Muslimov & Tsygan (1986). C is a constant of integration imposed by the boundary condition on the neutron star surface. The solution does not depend on the quantum number $m$ in the static limit. The constant C is chosen such that the asymptotic value of the potentials converge to their flat space-time counterpart at large distance. This represents our normalization of the magnetic field strength throughout the paper. We now give the exact analytical expression for the first four multipoles $\ell \in \{1, 2, 3, 4\}$. Useful expressions for the magnetic field components can be found in Appendix A. For brevity, here we only give exact and approximate expressions for the potentials $f_{\ell,m}^B$.

## 3.1 The magnetic dipole $\ell = 1$

We start our discussion with the general-relativistic magnetic dipole. Introducing the vector spherical harmonics expansion, the static aligned dipole frozen into the neutron star is conveniently written as

$$ f_{1,0}^{B} = 2 \sqrt{6\pi} \frac{B R^3}{R_s^2} \left[ \frac{\ln(1-x)}{x} + 1 + \frac{x}{2} \right] \approx -2 \sqrt{\frac{2\pi}{3}} \frac{B R^3}{r^2} \left[ 1 + \frac{3}{4} \frac{R_s}{r} \right] \qquad (30) $$

(recall that $x = R_s/r$). This solution already found by Ginzburg & Oznoy (1964) asymptotes to the flat space-time field at very large distances $r \gg R_s$. B is the characteristic strength of the magnetic field.

## 3.2 The magnetic quadrupole $\ell = 2$

Let us perform the same expansion to the general-relativistic magnetic quadrupole such that it can be expressed inside the star by

$$ f_{2,0}^{B} = -\frac{8}{3} \sqrt{\frac{10\pi}{3} \frac{B R^4}{R_s^3}} \frac{x(x(x+6)-24)+6(3x-4)\log(1-x)}{x^2} \qquad (31a) $$

$$ \approx -4 \sqrt{\frac{2\pi}{15}} \frac{B R^4}{r^3} \left[ 1 + \frac{4}{3} \frac{R_s}{r} \right]. \qquad (31b) $$

The static (2,1)-quadrupole frozen into the neutron star is conveniently written with the normalization

$$ f_{2,1}^{B} = -\frac{1}{2} \sqrt{\frac{3}{2}} f_{2,0}^{B}. \qquad (32) $$

The static (2,2)-quadrupole frozen into the neutron star is conveniently written with the normalization

$$ f_{2,2}^{B} = \sqrt{\frac{3}{2}} f_{2,0}^{B}. \qquad (33) $$

## 3.3 The magnetic hexapole $\ell = 3$

Next the magnetic hexapole existing inside the neutron star is written as

$$ f_{3,0}^{B} = \frac{5}{2} \sqrt{\frac{7\pi}{3} \frac{B R^5}{R_s^4}} \frac{12(6x^2 - 20x + 15)\log(1-x) + x(x(x+12)-150)+180}{x^3} \qquad (34a) $$

$$ \approx -2 \sqrt{\frac{\pi}{21}} \frac{B R^5}{r^4} \left[ 1 + \frac{15}{8} \frac{R_s}{r} \right]. \qquad (34b) $$
The static (3,1)-quadrupole frozen is normalized according to

$$f_{3,1}^{B} = -\frac{32}{\sqrt{3}} f_{3,0}^{B}. \qquad (35)$$

The static (3,2)-quadrupole frozen is normalized according to

$$f_{3,2}^{B} = \sqrt{\frac{2}{15}} f_{3,0}^{B}. \qquad (36)$$

The static (3,3)-quadrupole is normalized to

$$f_{3,3}^{B} = \frac{8}{\sqrt{5}} f_{3,0}^{B}. \qquad (37)$$

### 3.4 The magnetic octopole $ℓ = 4$

To finish this discussion about the lowest order general-relativistic multipoles, we give the exact solution to the octopole field for any azimuthal mode $m$. Inside the neutron star, the magnetic field is described a priori by a general octopolar expansion such that the magnetic potential useful for any $m$ is given with appropriate normalization by

$$f_{4,0}^{B} = \frac{336}{5} \sqrt{\frac{\pi}{5}} \frac{B R^6}{R_s^5 x^4} \left[ x(x(x(3x(x+20)-1570)+4620)-3360) + 60(5x(2(x-6)x+21)-56)\log(1-x) \right] \qquad (38a)$$

$$\approx 8 \sqrt{\frac{\pi}{5}} \frac{B R^6}{r^5} \left[ 1 + \frac{12}{5} \frac{R_s}{r} \right]. \qquad (38b)$$

The static (4,1)-quadrupole is conveniently written with the normalization

$$f_{4,1}^{B} = \frac{8}{\sqrt{5}} f_{4,0}^{B}. \qquad (39)$$

The static (4,2)-quadrupole is conveniently written with the normalization

$$f_{4,2}^{B} = \frac{1}{\sqrt{10}} f_{4,0}^{B}. \qquad (40)$$

The static (4,3)-quadrupole is conveniently written with the normalization

$$f_{4,3}^{B} = \sqrt{\frac{5}{7}} f_{4,0}^{B}. \qquad (41)$$

The static (4,4)-quadrupole is conveniently written with the normalization

$$f_{4,4}^{B} = \sqrt{\frac{10}{7}} f_{4,0}^{B}. \qquad (42)$$

All these expressions for the first multipoles will be very useful for forthcoming general-relativistic time-dependent simulations of Maxwell equations to look deeper into stationary solutions of the electromagnetic field in vacuum as well as in force-free magnetospheres.

## 4 EXACT ROTATING MULTIPOLE FIELDS WITHOUT FRAME DRAGGING

Frame dragging induces a coupling between different components rendering an analytical solution difficult to find. This is seen in the right-hand side term in equations (20) where terms in $\omega$ appear. However, in order to catch the essentials of general-relativistic effects on a rotating multipole, we start with a Schwarzschild background metric, thus neglecting frame dragging. We justify this approximation a posteriori when solving numerically the system of equations (20), at least for realistic neutron star parameters. This allows a treatment very similar to flat space-time as shown below.

### 4.1 General treatment

From the expansion into vector spherical harmonics, each function $f_{ℓ,m}^{D/B}$ has to satisfy a scalar wave equation in Schwarzschild space-time such that

$$\frac{\alpha^2}{r} \partial_r (\alpha^2 \partial_r (r f)) + \left( k_m^2 - \alpha^2 \frac{\ell(\ell+1)}{r^2} \right) f = 0. \qquad (43)$$

Using the tortoise coordinate $r_*$ and the new unknown function $\phi = rf$ we arrive at

$$\partial_{r_*}^2 \phi + \left( k_m^2 - \alpha^2 \frac{\ell(\ell+1)}{r^2} \right) \phi = 0. \qquad (44)$$
We are looking for solutions describing outgoing waves that reduce to $e^{i k_m r}$ in flat space-time, thus we introduce another unknown field $u$ as
$$ \phi = u e^{i k_m r_s}. \quad (45) $$

Therefore, $u$ is solution of the ordinary differential equation

$$ \alpha^2 u'' + \left( 2 i k_m + \frac{R_s}{r^2} \right) u' - \frac{\ell (\ell + 1)}{r^2} u = 0 \quad (46) $$

or written more explicitly

$$ r (r - R_s) u'' + (2 i k_m r^2 + R_s) u' - \ell (\ell + 1) u = 0. \quad (47) $$

Outgoing waves at infinity imposes a finite value for $u$ such that $\lim_{r \to +\infty} |u(r)| = u^\infty < \infty$. Changing to normalized radius defined by $z = r/R_s$ and wavenumber by $\epsilon_m = k_m R_s = R_s/r_L \ll 1$, this equation then becomes

$$ z (z - 1) u''(z) + (1 + 2 i \epsilon_m z^2) u'(z) - \ell (\ell + 1) u(z) = 0. \quad (48) $$

Equation (48) is known as the confluent Heun equation. It has two regular singular points located, respectively, at the origin of the coordinate system $z = 0$ and at the Schwarzschild radius $z = 1$ and one irregular singular point at spatial infinity $z = +\infty$ (Olver 2010). In the standard non-symmetrical canonical form given by Ronveaux & Arscott (1995), the confluent Heun equation is usually summarized by

$$ \frac{d^2 w}{dz^2} + \left( 4 p + \frac{\gamma}{z} + \frac{\delta}{z-1} \right) \frac{dw}{dz} + \frac{4 p \alpha z - \sigma}{z(z-1)} w = 0. \quad (49) $$

In the case of equation (48), the parameters are defined by

$$ 4 p = 2 i \epsilon_m \quad (50a) $$

$$ \gamma = -1 \quad (50b) $$

$$ \delta = 1 + 2 i \epsilon_m \quad (50c) $$

$$ \alpha = 0 \quad (50d) $$

$$ \sigma = \ell (\ell + 1). \quad (50e) $$

Local solutions to the confluent Heun equation around the regular singular points $z = 0, 1$ can be expressed as a Frobenius series expansion. Unfortunately, such expansions do not extend to $z = +\infty$ because their radius of converge is at most equal to the distance to the next singular point, thus a convergence around the point $z_0$ only within $|z - z_0| < 1$. The solution called radial function and denoted by $\text{Hc}^{(r)}(p, \alpha, \gamma, \delta, \sigma; z)$ remaining finite at infinite radius (Cook & Zalutskiy 2014) is defined by

$$ \lim_{z \to +\infty} z^{\alpha} \text{Hc}^{(r)}(p, \alpha, \gamma, \delta, \sigma; z) = 1. \quad (51) $$

The two independent local solutions at infinity are therefore

$$ \text{Hc}^{(r)}(p, \alpha, \gamma, \delta; z) \quad (52a) $$

$$ e^{-4 p z} \text{Hc}^{(r)}(-p, -\alpha + \gamma + \delta, \gamma, \delta, \sigma - 4 p \gamma; z). \quad (52b) $$

In our case, $\alpha = 0$ therefore $\lim_{z \to +\infty} \text{Hc}^{(r)}(p, \alpha, \gamma, \delta, \sigma; z) = 1$ which is indeed compatible with the boundary condition we impose on $u_{\ell, m}$. These constants should remain finite and different from zero. The other solution behaves like $e^{-4 p z} z^{\alpha - \gamma - \delta} = e^{-2 i \epsilon_m (z + \ln z)} = e^{-2 i k_m (r + R_s \ln(r/R_s))}$ which is the solution for ingoing wave from infinity that we discard in our analysis. The asymptotic series at infinity given as Thomé solution for the function $\text{Hc}^{(r)}(p, \alpha, \gamma, \delta, \sigma; z)$ is (Philipp & Perlick 2015)

$$ \text{Hc}^{(r)}(p, \alpha, \gamma, \delta, \sigma; z) = \sum_{n=0}^{+\infty} a_n z^{-n} \quad (53) $$

where the coefficients $a_n$ satisfy a three-term recurrence

$$ f_n a_{n+1} + g_n a_n + h_n a_{n-1} = 0 \quad (54a) $$

$$ f_n = -2 i k R_s (n + 1) \quad (54b) $$

$$ g_n = n (n + 1) - \ell (\ell + 1) \quad (54c) $$
$$h_n = -(n-1)(n+1) \qquad (54d)$$

leading unfortunately to the divergence of the series proposed in equation (53) (Ronveaux & Arscott 1995). Divergent series are not useless, they can approximate to good accuracy some functions having a divergent series as formal solution (Bender & Orszag 1999). Computations is performed by an optimal truncation of the series. Nevertheless, for the coefficients given by the three-term recurrence in equation (54), it is not easy to find the explicit form of the $a_n$. Moreover, the convergence is efficient only for $z \gg 1$. However, to impose boundary condition on the star, accurate values are also required for $z \ll 1$. Eventually, we have to resort to numerical methods to find accurate solutions in the whole integration domain.

In Newtonian gravity, $R_s = 0$ and the three-term recurrence simplifies into a two-term recurrence because $f_n = 0$ and it is finite, stopping whenever $n = \ell$. The solutions are polynomials in $z^{-1}$ and corresponds to the standard spherical Hankel functions $h_{\ell}^{(1)}(k_m r)$ as used by Pétri (2015). In the general-relativistic static limit, solutions are found by setting $\epsilon_m = 0$ in equation (48) which means $k_m = 0$ and therefore no rotation. In that case, the confluent Heun equation reduces to the hypergeometric differential equation of Section 3.

Series solutions of the confluent Heun equation (48) are found with Frobenius method and a judicious change of the independent variable (Leaver 1986). Because of the two regular singular points located, respectively, at $z = 0$ and $1$ and the irregular singular point at $z = +\infty$ we use a Jaffé transform to expand the solution around the regular singular point $z = 1$. According to Jaffé, the transform reads

$$x = \frac{z-1}{z} \qquad (55)$$

In that way the regular singular point $z = 0$ is rejected at $x = -\infty$, the regular singular point $z = 1$ at a value $x = 0$ and the irregular singular point $z = +\infty$ at $x = 1$. The differential equation then reads

$$x(1-x)^2 u''(x) + (2i \epsilon_m + (1-x)(1-3x)) u'(x) - \ell(\ell+1)u(x) = 0. \qquad (56)$$

This technique was used by Kearney, Kegeles & Cohen (1978) to solve electrodynamics problems in Schwarzschild metric applied to black holes and neutron stars. We look for series solution around $x = 0$ such that $u(x) = \sum_{n=0}^{+\infty} a_n x^n$. We know that the series is convergent up to the nearest singular point thus up to $x = 1$ which corresponds to $z = +\infty$. The solution is therefore valid in whole space outside the Schwarzschild radius. Equating the monomes of same order $x^n$, we arrive at the following recurrence relations

$$\ell(\ell+1)a_0 = (1+2i\epsilon_m)a_1 \qquad (57a)$$

$$(n+1)(n+1+2i\epsilon_m)a_{n+1} = [2n(n+1)+\ell(\ell+1)]a_n - (n+1)(n-1)a_{n-1}. \qquad (57b)$$

For very large $n$ with $n \gg 1$, the relation reduces to $a_{n+1} = 2a_n - a_{n-1}$ therefore the ratio $a_{n+1}/a_n = 2 - a_{n-1}/a_n$ leading to $\lim_{n \to +\infty} a_{n+1}/a_n = 1$. Thus, the radius of convergence of the series is indeed $\rho = 1$. To find a solution to our problem, we need to check that the series is also convergent at the point $x = 1$ that is on the circle of convergence. There is no general theorem on convergence on this circle and the problem is a delicate mathematical question with no definite answer. Actually the series does not converge at $x = 1$. It is thus impossible to satisfy the boundary condition at infinity with this series expansion. We leave the solution with the notation $\text{Hc}^{(r)}(p, \alpha, \gamma, \delta, \sigma; r/R_s)$ and do not give explicit expressions for it.

Close to the neutron star, we have $\epsilon z^2 \ll 1$ so that we can neglect this term in front of $u'$. The solution therefore reduces to the general-relativistic static multipole as described in Section 3. Let us write $\mathcal{H}_{\ell}^{(1)}(k_m r)$, the solution to the confluent Heun equation that remains finite at infinity. The solution corresponding to a spherical outgoing wave, that with $m > 0$, is

$$u_{\ell,m}^{\text{D/B}} = a_{\ell,m}^{\text{D/B}} \text{Hc}^{(r)}(p, \alpha, \gamma, \delta, \sigma; z). \qquad (58)$$

The constant of integration $a_{\ell,m}^{\text{B}}$ is determined from the boundary conditions on the neutron star surface noting that

$$u_{\ell m}^{\text{B}}(R) = R f_{\ell m}^{\text{B}}(R) e^{-i k_m r_{\ell m}(R)}. \qquad (59)$$

Moreover, we introduced the new function

$$\mathcal{H}_{\ell}^{(1)}(z) = \frac{\text{Hc}^{(r)}(p, \alpha, \gamma, \delta; z)}{z} e^{i k_m r_{\ell}} \qquad (60)$$

as a generalization of the spherical Hankel functions to the curved metric of Schwarzschild type. The constant of integration for the magnetic part is therefore

$$a_{\ell,m}^{\text{B}} = \frac{f_{\ell,m}^{\text{B}}(R)}{\mathcal{H}_{\ell}^{(1)}(k_m R)}. \qquad (61)$$

The constant of integration for the electric part obtained from the continuity of the tangential component of the electric field gives

$$\alpha^2 a_{\ell,m}^{\text{D}} \partial_r (r \mathcal{H}_{\ell}^{(1)}(k_m r)) \Big|_{r=R} = \epsilon_0 R \Omega \left[ \sqrt{(\ell+1)(\ell-1)} J_{\ell,m} f_{\ell-1,m}^{\text{B}}(R) - \sqrt{\ell(\ell+2)} J_{\ell+1,m} f_{\ell+1,m}^{\text{B}}(R) \right]. \qquad (62)$$

From this definition, the electric and magnetic potentials are deduced according to

$$f_{\ell,m}^{\text{D/B}}(r) = a_{\ell,m}^{\text{D/B}} \mathcal{H}_{\ell}^{(1)}(k_m r). \qquad (63)$$
<table><caption>**Table 1.** Spin-down contribution from the magnetic and electric multipoles in the limit of point multipoles ($R/r_L \to 0$). The last column gives the percentage of spin-down luminosity attributed to the dominant electric multipole.</caption><thead><tr><th>Multipole (ℓ, m)</th><th>C<sub>ℓ,m</sub></th><th>C<sub>ℓ-1,m</sub></th><th>C<sub>ℓ+1,m</sub></th><th>C<sub>ℓ-1,m</sub>/S<sub>ℓ,m</sub></th></tr></thead><tbody><tr><td>(1, 1)</td><td>1</td><td>-</td><td>0</td><td>0</td></tr><tr><td>(2, 1)</td><td><sup>1</sup>&frasl;<sub>9</sub></td><td><sup>3</sup>&frasl;<sub>5</sub></td><td>0</td><td><sup>27</sup>&frasl;<sub>32</sub> ≈ 84 per cent</td></tr><tr><td>(2, 2)</td><td><sup>64</sup>&frasl;<sub>9</sub></td><td>-</td><td>0</td><td>0</td></tr><tr><td>(3, 1)</td><td><sup>1</sup>&frasl;<sub>225</sub></td><td><sup>16</sup>&frasl;<sub>315</sub></td><td>0</td><td><sup>80</sup>&frasl;<sub>87</sub> ≈ 92 per cent</td></tr><tr><td>(3, 2)</td><td><sup>256</sup>&frasl;<sub>225</sub></td><td><sup>138</sup>&frasl;<sub>63</sub></td><td>0</td><td><sup>25</sup>&frasl;<sub>39</sub> ≈ 64 per cent</td></tr><tr><td>(3, 3)</td><td><sup>729</sup>&frasl;<sub>25</sub></td><td>-</td><td>0</td><td>0</td></tr><tr><td>(4, 1)</td><td><sup>1</sup>&frasl;<sub>11025</sub></td><td><sup>1</sup>&frasl;<sub>567</sub></td><td>0</td><td><sup>175</sup>&frasl;<sub>184</sub> ≈ 95 per cent</td></tr><tr><td>(4, 2)</td><td><sup>1024</sup>&frasl;<sub>11025</sub></td><td><sup>1024</sup>&frasl;<sub>2835</sub></td><td>0</td><td><sup>35</sup>&frasl;<sub>44</sub> ≈ 80 per cent</td></tr><tr><td>(4, 3)</td><td><sup>6561</sup>&frasl;<sub>1225</sub></td><td><sup>27</sup>&frasl;<sub>5</sub></td><td>0</td><td><sup>245</sup>&frasl;<sub>488</sub> ≈ 50 per cent</td></tr><tr><td>(4, 4)</td><td><sup>1048576</sup>&frasl;<sub>11025</sub></td><td>-</td><td>0</td><td>0</td></tr></tbody></table>

The functions $\mathcal{H}_{\ell}^{(1)}(k_m r)$ are well defined as solutions to a linear boundary value problem but are still ill studied by Mathematicians. Confluent Heun functions have many applications in many areas of physics due to their general expression but as a drawback, they are hard to study. There is no efficient implementation in any mathematical software so far due to their intrinsic complexity.

## 4.2 Solution for one multipole

Useful exact solutions are given for particular multipole fields with fixed numbers $(\ell, m)$. Let us assume that inside the star, the magnetic field is solely represented by the function $f_{\ell,m}^B(r)$. Then, the only non-vanishing magnetic field coefficient is given by equation (61) (we discard axisymmetric cases $m=0$ which represent static solutions). Moreover, the two non-vanishing electric field coefficients are (if $\ell=1$ only one solution exists, see the dipole case below) for $m > 0$

$$ \alpha^2 a_{\ell+1,m}^D \left. \partial_r (r \mathcal{H}_{\ell+1}^{(1)}(k_m r)) \right|_{r=R} = \varepsilon_0 R \Omega \sqrt{\ell(\ell+2)} J_{\ell+1,m} f_{\ell,m}^B(R) \quad (64a) $$

$$ \alpha^2 a_{\ell-1,m}^D \left. \partial_r (r \mathcal{H}_{\ell-1}^{(1)}(k_m r)) \right|_{r=R} = -\varepsilon_0 R \Omega \sqrt{(\ell-1)(\ell+1)} J_{\ell,m} f_{\ell,m}^B(R). \quad (64b) $$

We conclude that the solution is fully specified by the three constants of integration $(a_{\ell,m}^B, a_{\ell+1,m}^D, a_{\ell-1,m}^D)$ as already demonstrated in flat space-time. The Poynting flux associated with this particular solution is for the relevant cases $m > 0$ (it vanishes for axisymmetric cases $m=0$)

$$ P_{\ell,m} = \frac{c}{2\mu_0} \left[ |a_{\ell,m}^B|^2 + \mu_0^2 c^2 (|a_{\ell-1,m}^D|^2 + |a_{\ell+1,m}^D|^2) \right] \quad (65a) $$

$$ = \frac{c|f_{\ell,m}^B(R)|^2}{2\mu_0} S_{\ell,m} \quad (65b) $$

$$ S_{\ell,m} = \frac{1}{|\mathcal{H}_{\ell}^{(1)}(k_m R)|^2} + \frac{R^2}{r_L^2} \left( \frac{\ell(\ell+2) J_{\ell+1,m}^2}{|\partial_r(r \mathcal{H}_{\ell+1}^{(1)}(k_m r))|_R^2} + \frac{(\ell-1)(\ell+1) J_{\ell,m}^2}{|\partial_r(r \mathcal{H}_{\ell-1}^{(1)}(k_m r))|_R^2} \right). \quad (65c) $$

For the special case $\ell=1$, the constant $a_{\ell-1,m}^D$ does not exist.

Expression (65) generalizes the Newtonian formula found in Pétri (2015). At this point, it is opportune to stress that the Poynting flux arises mainly not from the single magnetic multipole $(\ell, m)$ depicted by the constant $a_{\ell,m}^B$, but from the lowest order multipole $\ell'$, let it be magnetic or electric. Because of the boundary conditions imposed on the neutron star surface where electric multipoles of order $\ell \pm 1$ are induced from the magnetic multipole of order $(\ell, m)$, there are cases where the electric multipole is dominant compared to the magnetic multipole $(\ell, m)$. For the magnetic dipole with quantum numbers $(\ell, m) = (1, 1)$, the only electric multipole has quantum numbers $(\ell', m) = (2, 1)$, thus represents a higher order multipole and therefore remains negligible to first approximation in spin rate. The Poynting flux $P_{\ell,m}$ is mainly attributed to the constant of integration $a_{\ell,1}^B$ that is the magnetic dipole.

Any electric multipole must satisfy the requirement $m \le \ell' = \ell \pm 1$. Moreover electric multipoles are dominant if their quantum number satisfies $\ell' < \ell$, thus the only choice is $\ell' = \ell - 1 < \ell$ and this requires $\ell \ge m + 1$. We conclude that for single magnetic multipoles given by $(\ell \ge m + 1, m)$, the electric multipole of order $(\ell - 1, m)$ dominates the spin-down rate. Let us prove this fact in the flat space-time geometry where exact analytical formulae for any ratio $R/r_L$ have been found. Table 1 summarizes the weights associated with each multipole and the percentage of luminosity attributed to the electric multipole (last column). Note that the normalized total luminosity is split into contributions from different multipoles according to $S_{\ell,m} = C_{\ell,m} + C_{\ell-1,m} + C_{\ell+1,m}$. $C_{\ell,m}$ represents the contribution from the magnetic multipole through the constant $a_{\ell,m}^B$, $C_{\ell-1,m}$ represents the contribution from the electric multipole through the constant $a_{\ell-1,m}^D$ and $C_{\ell+1,m}$ represents the contribution from the magnetic multipole through the constant $a_{\ell+1,m}^B$.
represents the contribution from the electric multipole through the constant $a_{\ell+1,m}^D$. This table emphasizes three points. First, the electric multipole ($\ell + 1, m$) never contributes to the spin-down in the limit of a point multipole. Second, the magnetic multipole ($\ell = m, m$) is the only relevant contributor to the spin-down. Third, for any multipole ($\ell > m, m$), the contribution from the electric multipole ($\ell - 1, m$) can be substantial and is always larger than the luminosity emanating from the magnetic multipole ($\ell, m$). Consequently, apart for the magnetic dipole, in order to estimate the spin-down luminosity of a single magnetic multipole of order $\ell \ge 2$, we need to take care of the associated electric multipole to properly deduce the energy loss rate. The discussion focused on the Newtonian case, but the same obviously applies to a general-relativistic magnetic multipole as will be shown from the numerical solutions.

This apparently misleading results arises because although a factor $(R/r_L)^2$ appears in front of the electric multipoles, it is compensated by the lower order of one of these electric multipoles containing a spherical Hankel function behaving like $|\partial_r(r \mathcal{H}_{\ell-1}^{(1)}(k_m r))|_R^{-2}$. As a result, $(R/r_L)^2 |\partial_r(r \mathcal{H}_{\ell-1}^{(1)}(k_m r))|_R^{-2}$ is of the same order in $R/r_L$ as the magnetic multipole $|\mathcal{H}_{\ell}^{(1)}(k_m R)|^{-2}$ for $R/r_L \ll 1$. Actually, we showed that contrary to being negligible, this electric multipole is dominant for the energy loss rate as soon as $\ell > m$. Indeed, in flat space-time, the spherical Hankel functions can be expanded to lowest order when $k_m R \ll 1$. In that case, the contributions to spin-down luminosities are

$$S_{\ell,m} = m^{2\ell} (k R)^{2\ell+2} \left[ \frac{m^2}{((2\ell-1)!!)^2} + \frac{(\ell+1) J_{\ell,m}^2}{(\ell-1)((2\ell-3)!!)^2} \right] \quad (66)$$

where the first term accounts for the magnetic multipole of order ($\ell, m$) and the second term accounts for the electric multipole of order ($\ell - 1, m$) (for $\ell > 1$). We neglect terms of higher order in the product $kR$, especially those arising from the ($\ell + 1, m$) electric multipole. The term in bracket in front of $R_L^2$ in equation (65) is not necessarily a higher order correction of the spin-down in the point multipole limit compared to the first term in $S_{\ell,m}$. It is actually exactly of the same order as the associated magnetic multipole.

There exist no simple general-relativistic expression for the spherical Hankel functions $\mathcal{H}_{\ell}^{(1)}(k_m R)$. Nevertheless, following Rezzolla & Ahmedov (2004), we could adopt the flat space-time counterpart as a good approximation taking into account the gravitational redshift of the spin frequency at the surface by noting $\Omega_R = \Omega/\alpha_R$, thus replacing $\mathcal{H}_{\ell}^{(1)}(k_m R) = h_{\ell}^{(1)}(m \Omega_R R/c)$. Using the asymptotic expression for $h_{\ell}^{(1)}(x)$ for $x \ll 1$, the major contribution to the spin-down luminosity should be

$$P_{\ell,m} = \frac{c|f_{\ell,m}^{\text{B}}(R)|^2}{2 \mu_0} \frac{1}{(2(\ell-1)!!)^2} \left( \frac{m \Omega R}{\alpha_R c} \right)^{2\ell+2} \quad (67)$$

Compared to flat space-time, we note two correcting factors. The first one arising from the amplification of the magnetic field strength at the stellar surface through the factor $|f_{\ell,m}^{\text{B}}(R, R_s)/f_{\ell,m}^{\text{B}}(R, R_s = 0)|^2$ and the second from the rotation frequency shift through a factor $\alpha_R^{-(2\ell+2)}$. These are the results already discussed by Rezzolla & Ahmedov (2004). The general-relativistic luminosity increase is therefore completely determined by the Schwarzschild radius $R_s$ of the star. There are no corrections including $r_L$. Unfortunately as we show in this paper and as was already found in Pétri (2013) by time-dependent numerical simulations, this estimate is not correct. For extremely slow-rotation rate, when $r_L \to +\infty$, the general-relativistic spin-down reduces to the Newtonian expression, their ratio tending to one.

Good estimates of the Poynting flux are given by the strength of the electromagnetic field at the light cylinder. Note however that there is no boundary condition imposed on this light cylinder. This fictitious surface represents the transition between the near quasi-static regime and the asymptotic wave zone. Any multipolar electromagnetic field is completely and uniquely determined by the outgoing wave solution at large distances and by fixing the field value at the stellar surface. We therefore enforce to boundary conditions as required by any elliptic partial differential equation such as the stationary wave equation leading to Helmholtz equation. If the light cylinder is rejected to large distances, with our normalization, the magnetic dipole geometry looks very similar to flat space-time with no gravitational perturbations and the relevant radiating field is close to its Minkowski expression. This agrees with our previous semi-analytical results and with the spin-down rate extracted from time-dependent numerical simulations.

Strictly speaking, in general relativity, the light cylinder is located at a surface where the local speed reach the speed of light. Taking into account our 3+1 foliation of the metric, we found (Ruiz et al. 2014)

$$\gamma_{ab} (v^a + \beta^a c) (v^b + \beta^b c) = \alpha^2 c^2. \quad (68)$$

Neglecting the shift vector at the location of the light cylinder, to first order the light cylinder radius is shifted inwards down to

$$r_{\text{L}}^{\text{GR}} = r_{\text{L}} \left( 1 - \frac{R_s}{2 r_{\text{L}}} \right). \quad (69)$$

In practice, these corrections can be neglected and we always assumed that $r_{\text{L}}^{\text{GR}} \approx r_{\text{L}} = c/\Omega$. An estimate of the spin-down contribution from a multipole of order $\ell$ can be guessed from the fact that the static magnetic field decreases like

$$B_{\ell} = B \frac{R^{\ell+2}}{r^{\ell+2}} \quad (70)$$

At the light cylinder, electric and magnetic fields have equal normalized intensities such that $E_{\ell} = c B_{\ell}$ and the Poynting flux is given by

$$S_{\ell} = \frac{E_{\ell} B_{\ell}}{\mu_{0}} = \frac{c}{\mu_{0}} B^{2} \frac{R^{2 \ell+4}}{c^{2 \ell+4}} \Omega^{2 \ell+4} \quad (71)$$
<figure><img src="image_3.png" /><figcaption>Figure 1. Comparison of the spherical Hankel function for $ℓ = 1$ in flat and Schwarzschild space-time for $R/r_L = 0.1$ and several ratio $R_s/R$ as shown in the legend. Real and imaginary parts are shown separately.</figcaption></figure>

The total radial flux through a surface of radius $r_L$ is

$$L_{\ell} = 4 \pi r_{L}^{2} S_{\ell} = \frac{4 \pi}{\mu_{0} c^{2} \ell + 1} B^{2} R^{2 \ell + 4} \Omega^{2 \ell + 2} \quad (72)$$

Apart from a constant numerical factor different from $4\pi$, these estimates are very close to the exact expressions. They show the typical behaviour as $\Omega^{2\ell+2}$ valid for a multipole of order $\ell$. These straightforward calculations point out the relation between the field intensities at the stellar surface, their intensities at the light cylinder end the related Poynting flux. No boundary conditions are enforced at the light cylinder. The only assumption is that the wave zone starts at this location with the actual value of the field in this region.

The above reasoning applied to general-relativistic rotating dipoles and multipoles lead to the conclusion that when the light cylinder is rejected to infinity, that is very slowly rotating stars, the spin-down tends to its Minkowski rate, except for subtle case when $\ell > m$ as will be demonstrated later.

To support our conclusion on firm basis, we show very accurate numerical solutions of spherical Hankel functions in general relativity up to $\ell = 4$. They are thoroughly investigated in the next subsection.

## 4.3 Generalized spherical Hankel functions

The wave equations in the Schwarzschild metric and in the slow-rotation metric are solved numerically to high accuracy with an expansion on to rational Chebyshev polynomials. Here, we give some useful numerical approximations for these general-relativistic spherical Hankel functions for different compactness $R/R_s$ and spin rate $R/r_L$. The results of our integration of the boundary value problems are given by the functions shown in Fig. 1 for the mode $\ell = 1$, in Fig. 2 for the mode $\ell = 2$, in Fig. 3 for the mode $\ell = 3$ and in Fig. 4 for the mode $\ell = 4$. To ease comparison, we use the same boundary conditions on both end of the integration interval. The solution must vanish at infinity and matches the flat spherical Hankel functions on the left boundary, $H_{\ell}^{(1)}(R/r_L) = h_{\ell}^{(1)}(R/r_L)$.

With these same boundary conditions, the curved spherical Hankel functions suffer from an amplitude decrease for large distances compared to their flat space-time equivalents. The decrease gets stronger for larger compactness. This means that taking the flat function as an exact solution to the curved space-time problem but adapting the inner boundary conditions according to general-relativistic corrections overestimates the Poynting flux as measured by a distant observer. A correct answer requires a careful integration of the curved space-time wave equation as presented in this paper.

The evaluation of the spin-down luminosity requires the value of the functions $u^{D/B}$ at infinity and denoted by $u_{\infty}^{D/B}$. As a comparison between flat and curved space-time, we have shown in Table 2 these values for the same parameters $R/R_s$ and $R/r_L$. These values are useful to compute the energy loss from a magnetic multipole of order ($\ell = m, m$), but not for the other multipoles. The reason has been exposed previously.
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Comparison of the spherical Hankel function for $\ell = 2$ in flat and Schwarzschild space-time for $R/r_L = 0.1$ and several ratio $R_s/R$ as shown in the legend. Real and imaginary parts are shown separately.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 3.** Comparison of the spherical Hankel function for $\ell = 3$ in flat and Schwarzschild space-time for $R/r_L = 0.1$ and several ratio $R_s/R$ as shown in the legend. Real and imaginary parts are shown separately.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Comparison of the spherical Hankel function for $\ell = 4$ in flat and Schwarzschild space-time for $R/r_L = 0.1$ and several ratio $R_s/R$ as shown in the legend. Real and imaginary parts are shown separately.</figcaption></figure>

<table><caption>**Table 2.** Limit of the spherical Hankel functions for several multipoles $\ell$. The table gives the values of $u_{\ell,1}^{B,\infty}$. $R_s/R = 0$ represents flat space-time, $R_s/R = \{0.2, 0.4\}$ (Schw) Schwarzschild metric and $R_s/R = \{0.2, 0.4\}$ (SRNS) the slowly rotating neutron star metric.</caption><thead><tr><th rowspan="2">ℓ</th><th rowspan="2">R<sub>s</sub>/R = 0</th><th colspan="2">R<sub>s</sub>/R = 0.2</th><th colspan="2">R<sub>s</sub>/R = 0.4</th></tr><tr><th>(Schw)</th><th>(SRNS)</th><th>(Schw)</th><th>(SRNS)</th></tr></thead><tbody><tr><td>1</td><td>-1</td><td>-0.872+0.0675*i</td><td>-0.8719+0.0676*i</td><td>-0.732+0.0926*i</td><td>-0.731+0.0926*i</td></tr><tr><td>2</td><td>+i</td><td>0.0700+0.768*i</td><td>0.0700+0.768*i</td><td>0.0853+0.551*i</td><td>0.0853+0.551*i</td></tr><tr><td>3</td><td>+1</td><td>0.682-0.0674*i</td><td>0.682-0.0674*i</td><td>0.424-0.0724*i</td><td>0.424-0.0723*i</td></tr><tr><td>4</td><td>-i</td><td>-0.0635-0.609*i</td><td>-0.0634-0.609*i</td><td>-0.0599-0.330*i</td><td>-0.0599-0.329*i</td></tr></tbody></table>

We also integrated the wave equation with frame dragging and did not find significant deviation from Schwarzschild solutions as long as $R/r_L \ll 1$. The solution in Schwarzschild space-time remain valid to very good accuracy for a slowly rotating neutron star metric. They are not represented here as they overlap to previous plots.

For concreteness, we switch now to explicit application of low-order multipole solutions in order to get more physical insight into their properties. For the remainder of this paper, we focus on some illuminating cases such as the low-order multipoles: dipole, quadrupole, hexapole and octopole.

## 5 APPROXIMATE ROTATING MULTIPOLE SOLUTIONS

Including frame-dragging effects into the picture of rotating multipoles renders the problem analytically infeasible. The system of equations satisfied by the transformed potentials $u_{\ell,m}^{D/B}$ has to be solved numerically. Confluent Heun equations and its generalization are solved by spectral methods, expanding the unknown functions on to Chebyshev polynomials. More precisely, elliptic equations satisfied by the potentials are efficiently solved by a radial expansion on to rational Chebyshev functions as defined and explained in Boyd (2001). Boundary conditions taking into account frame dragging and space curvature have been given in equations (21) and (59). Moreover, these transformed potentials should converge to a finite non-vanishing value at infinite radius.

In the following subsections, we given accurate numerical solutions to the electromagnetic field of a single rotating multipole up to octupole order $\ell = 4$. To emphasize the effect of space curvature and frame dragging, we compute two kind of solutions. The first one assumes a background Schwarzschild metric thus without dragging and the second one includes this dragging adopting the slowly rotating approximation. Actually, we will show that except for unrealistically short periods with $P < 1$ ms, this frame dragging can safely be ignored as far as the Poynting flux is concerned.
<table><caption>**Table 3.** Normalized spin-down luminosity for the perpendicular rotating dipole in Schwarzschild background. The analytical approximation equation (75) is shown in the last line (RA).</caption><thead><tr><th>R/r<sub>L</sub> / R<sub>s</sub>/R</th><th>0</th><th>0.1</th><th>0.2</th><th>0.3</th><th>0.4</th><th>0.5</th></tr></thead><tbody><tr><td>0.01</td><td>0.9999</td><td>1.003</td><td>1.006</td><td>1.009</td><td>1.012</td><td>1.015</td></tr><tr><td>0.02</td><td>0.9996</td><td>1.005</td><td>1.012</td><td>1.018</td><td>1.025</td><td>1.032</td></tr><tr><td>0.05</td><td>0.9975</td><td>1.013</td><td>1.029</td><td>1.045</td><td>1.063</td><td>1.081</td></tr><tr><td>0.1</td><td>0.9901</td><td>1.020</td><td>1.053</td><td>1.088</td><td>1.125</td><td>1.166</td></tr><tr><td>0.2</td><td>0.9615</td><td>1.018</td><td>1.083</td><td>1.156</td><td>1.239</td><td>1.336</td></tr><tr><td>RA</td><td>1.000</td><td>1.444</td><td>2.171</td><td>3.434</td><td>5.800</td><td>10.70</td></tr></tbody></table>

Because of frame-dragging effects that couple $f_{\ell,m}^B$ modes to $f_{\ell,m}^D$ modes according to the right-hand side terms in equation (13), there is an infinity of modes $(\ell, m)$ that are induced by a single rotating multipole. However, because these modes remain much weaker than the fundamental single multipole imposed by the star, we neglect them in the subsequent treatment. Time-dependent numerical simulations at the end of this paper taking into account many multipoles $(\ell' > \ell, m')$ will confirm our expectations.

## 5.1 Dipole solution

The rotating dipole has been extensively investigated in flat space-time. If space-time curvature is including, the transformed potentials have to satisfy the following generalization of the Helmholtz equation for a stationary wave. The two relevant functions are $u_{2,1}^D$ and $u_{1,1}^B$ verifying the coupled system

$$ \frac{1}{r} \left[ r (r - R_s) u_{2,1}^{D''} + (2 i k r^2 + R_s) u_{2,1}^{D'} - 6 u_{2,1}^{D} \right] + k^2 \left[ \left(1 - \frac{\omega}{\Omega}\right)^2 - 1 \right] r \frac{u_{2,1}^{D}}{\alpha^2} = 3 \sqrt{\frac{3}{5}} \epsilon_0 \omega u_{1,1}^B \quad (73a) $$

$$ \frac{1}{r} \left[ r (r - R_s) u_{1,1}^{B''} + (2 i k r^2 + R_s) u_{1,1}^{B'} - 2 u_{1,1}^{B} \right] + k^2 \left[ \left(1 - \frac{\omega}{\Omega}\right)^2 - 1 \right] r \frac{u_{1,1}^{B}}{\alpha^2} = 3 \sqrt{\frac{3}{5}} \mu_0 \omega u_{2,1}^{D}. \quad (73b) $$

Without frame dragging that is when $\omega = 0$, these equations decouple in two scalar ordinary differential equations as is the case for Minkowski space-time. In any case, the electric part has to satisfy the following boundary condition

$$ \alpha^2 \partial_t (u_{2,1}^D e^{i k r_s}) = +\epsilon_0 r \tilde{\omega} \sqrt{\frac{3}{5}} f_{1,1}^B. \quad (74) $$

The solution is completely and uniquely determined by the aforementioned equations and boundary conditions. In all subsequent results, as a check of the numerical implementation and high accuracy of our algorithm, we compare the output of the numerical approximation with the exact analytical solution expressed by spherical Hankel functions that are given by Pétri (2015).

For a dipole rotating in a Schwarzschild background metric, thus neglecting frame dragging, the spin-down luminosity for different stellar compactness parameters $\Xi = R_s/R$ and different periods expressed by the ratio $a = R/r_L$ is shown in Table 3. The column with $R_s/R = 0$ corresponds to the Minkowski solution and serves as a reference to prove the high precision of the method. Indeed, we found at least 10 digits of precision compared to the analytical expression. This clearly points out the very good set of basis functions used to solve the system. As a general trend, the spin-down monotonically increases with compactness for fixed period reaching up to 1.336 times the reference value. As a comparison with the analytical expression found by Rezzolla & Ahmedov (2004), the corresponding general-relativistic enhancement is shown in the last line of the table. It depends on $R_s/R$ but remains independent of the rotation rate so no need to specify $R/r_L$. The formula is

$$ \frac{L_{\text{GR}}}{L_{\text{Newt}}} = \left( \frac{f_{1,0}^{B}(R_s, R)}{f_{1,0}^{B}(R_s = 0, R)} \right)^2 \left( 1 - \frac{R_s}{R} \right)^{-2} \quad (75) $$

The first correction is induced by the magnetic field amplification and the second by gravitational redshift of the spin frequency. It should be a good guess for a point dipole, thus for $R/r_L \ll 1$. A quick comparison with the first line of the table with $R/r_L = 0.01$ clearly emphasizes the mismatch between both estimates. Rezzolla & Ahmedov (2004) are severely overestimating the Poynting flux for very compact stars.

When frame dragging is switch on, the Poynting flux decreases slightly compared to the previous case as seen in Table 4. This trend is easily explained by the boundary conditions for the electric field which contain the term $\tilde{\omega}$ that is a decrease in the effective rotation rate of the star as measured by a local observer explaining a weaker radiating electric field contribution. However, for realistic neutron star parameters, we would expect at most $R/r_L \le 0.1$ and $R_s/R \le 0.4$. In that parameter range, the discrepancy between Tables 3 and 4 remains less than 1 per cent. Thus, for practical purposes when discussing the spin-down luminosities, frame-dragging effects are irrelevant. The plots in Fig. 5 give a synthetic overview of the full set of results showing the increase in spin-down luminosity compared to flat space-time given by the ratio

$$ S(a, \Xi) = \frac{L_{\text{GR}}}{L_{\text{flat}}}. \quad (76) $$
<table><caption>**Table 4.** Normalized spin-down luminosity for the perpendicular rotating dipole in the slowly rotating metric background.</caption><thead><tr><th>R/r<sub>L</sub> / R<sub>s</sub>/R</th><th>0</th><th>0.1</th><th>0.2</th><th>0.3</th><th>0.4</th><th>0.5</th></tr></thead><tbody><tr><th>0.01</th><td>0.9999</td><td>1.003</td><td>1.006</td><td>1.009</td><td>1.012</td><td>1.015</td></tr><tr><th>0.02</th><td>0.9996</td><td>1.005</td><td>1.012</td><td>1.018</td><td>1.025</td><td>1.031</td></tr><tr><th>0.05</th><td>0.9975</td><td>1.012</td><td>1.028</td><td>1.045</td><td>1.062</td><td>1.080</td></tr><tr><th>0.1</th><td>0.9901</td><td>1.020</td><td>1.052</td><td>1.086</td><td>1.123</td><td>1.163</td></tr><tr><th>0.2</th><td>0.9615</td><td>1.017</td><td>1.079</td><td>1.149</td><td>1.229</td><td>1.321</td></tr></tbody></table>

<figure><img src="image_7.png" /><figcaption>**Figure 5.** Spin-down luminosity of the orthogonal dipole in Schwarzschild and slowly rotating neutron star metric approximation.</figcaption></figure>

Each curve corresponds to the ratio between a compactness $\Xi > 0$ and the reference case $\Xi = 0$. There is always an increase in the Poynting flux whatever the period $a$ and compactness $\Xi$. Solid lines with the same colour correspond to a same compactness. The Schwarzschild luminosity always lies above the slowly rotating neutron star metric. This is also valid for the higher multipoles shown in the following paragraphs.

To emphasize the difference between the numerical solution and some analytical approximate expressions, we plot the real and imaginary parts of the radial dependence of the magnetic field for a dipole in Schwarzschild space-time with the parameters $R/r_L = 0.1$ and $R/R_s = 2$. The discrepancies are highlighted in Fig. 6. The numerical solution is shown in red solid line (num in the legend). The expression deduced according to Rezzolla & Ahmedov (2004) does not reproduce the solution to good accuracy, green solid line (RA in the legend). Actually, we found that the Newtonian expression, in solid blue line ($h_1^{(1)}$ in the legend), assuming variations according to spherical Hankel functions $h_p^{(1)}(k_m r)$, regardless of any general-relativistic effect such as gravitation redshift or stellar boundary condition, gives the best approximation to the numerical solution. Corrections are only significant for unrealistically high rotation rates $R/r_L \lesssim 1$. Our conclusions remain valid for a slowly rotating space-time metric as long as $R/r_L \ll 1$. This explains why the general-relativistic luminosity equals approximately the Minkowski results. For a distant observer, there is no distinction between general-relativistic or Minkowski metric. Obviously, for an observer sitting on the star, the measurements would be very different.

## 5.2 Quadrupole solution

The rotating quadrupole has been less extensively investigated, even in flat space-time. The three relevant transformed potentials for the $m = 1$ mode are $u_{1,1}^D$, $u_{3,1}^D$ and $u_{2,1}^B$. They verify the coupled system for the three unknown functions

$$ \frac{1}{r} \left[ r(r-R_s)u_{1,1}^{D''} + (2ikr^2+R_s)u_{1,1}^{D'} - 2u_{1,1}^D \right] + k^2 \left[ \left(1-\frac{\omega}{\Omega}\right)^2 - 1 \right] r \frac{u_{1,1}^D}{\alpha^2} = -3 \sqrt{\frac{3}{5}} \epsilon_0 \omega u_{2,1}^B \quad (77a) $$

$$ \frac{1}{r} \left[ r(r-R_s)u_{2,1}^{B''} + (2ikr^2+R_s)u_{2,1}^{B'} - 6u_{2,1}^B \right] + k^2 \left[ \left(1-\frac{\omega}{\Omega}\right)^2 - 1 \right] r \frac{u_{2,1}^B}{\alpha^2} = -3 \mu_0 \omega \left( \sqrt{\frac{3}{5}} u_{1,1}^D - \frac{8}{\sqrt{35}} u_{3,1}^D \right) \quad (77b) $$

$$ \frac{1}{r} \left[ r(r-R_s)u_{3,1}^{D''} + (2ikr^2+R_s)u_{3,1}^{D'} - 12u_{3,1}^D \right] + k^2 \left[ \left(1-\frac{\omega}{\Omega}\right)^2 - 1 \right] r \frac{u_{3,1}^D}{\alpha^2} = 3 \sqrt{\frac{8}{35}} \epsilon_0 \omega u_{2,1}^B. \quad (77c) $$
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Comparison between the numerical, red solid line and the analytical expressions, green and blue solid lines, for the radial dependence of the magnetic field profile. The parameters are $R/r_L = 0.1$ and $R/R_s = 2$ for Schwarzschild space-time.</figcaption></figure>

<table><caption>**Table 5.** Normalized spin-down luminosity for the $m = 1$ rotating quadrupole in Schwarzschild background.</caption><thead><tr><th>R/r<sub>L</sub> / R<sub>s</sub>/R</th><th>0</th><th>0.1</th><th>0.2</th><th>0.3</th><th>0.4</th><th>0.5</th></tr></thead><tbody><tr><th>0.01</th><td>1.000</td><td>1.164</td><td>1.382</td><td>1.680</td><td>2.105</td><td>2.742</td></tr><tr><th>0.02</th><td>1.000</td><td>1.168</td><td>1.391</td><td>1.697</td><td>2.132</td><td>2.788</td></tr><tr><th>0.05</th><td>1.001</td><td>1.181</td><td>1.420</td><td>1.749</td><td>2.221</td><td>2.934</td></tr><tr><th>0.1</th><td>1.007</td><td>1.206</td><td>1.474</td><td>1.846</td><td>2.386</td><td>3.214</td></tr><tr><th>0.2</th><td>1.031</td><td>1.272</td><td>1.604</td><td>2.081</td><td>2.795</td><td>3.929</td></tr></tbody></table>

The appropriate boundary conditions for the electric part are given by

$$ \alpha^2 \partial_t (u_{1,1}^D e^{i k r_s}) = -\epsilon_0 r \tilde{\omega} \sqrt{\frac{3}{5}} f_{2,1}^B \quad (78a) $$

$$ \alpha^2 \partial_t (u_{3,1}^D e^{i k r_s}) = +\epsilon_0 r \tilde{\omega} \frac{8}{\sqrt{35}} f_{2,1}^B. \quad (78b) $$

The spin-down luminosity obtained from the numerical integration of this coupled system is summarized in Table 5 for the Schwarzschild metric. Here again, as a check the flat space-time results are compared to analytical expressions and the precision remains better than 10 digits. The luminosity increase with respect to compactness is more sensitive than for the dipole $m = 1$ mode up to 4 times larger than the flat space-time quadrupole. If frame dragging is switch on, the decrease in luminosity can be quantified by comparing Table 5 with Table 6. We always observe a lowering of the spin-down rate because of the boundary conditions imposed on the electric part for the same reason as for the dipole. Frame dragging decreases the effective stellar rotation rate for a local observer and thus the expected electric field strength contribution to the radiating field. The plots in Fig. 7 give a synthetic compilation of the results showing the increase in spin-down luminosity compared to flat space-time reference values.

For the $m = 2$, the two relevant transformed potentials are $u_{3,2}^D$ and $u_{2,2}^B$. They verify the coupled system of two unknown functions

$$ \frac{1}{r} \left[ r (r - R_s) u_{2,2}^{D''} + (2 i k_2 r^2 + R_s) u_{2,2}^{B'} - 6 u_{2,2}^B \right] + k_2^2 \left[ \left(1 - \frac{\omega}{\Omega}\right)^2 - 1 \right] r \frac{u_{2,2}^B}{\alpha^2} = 6 \sqrt{\frac{2}{7}} \mu_0 \omega u_{3,2}^D \quad (79a) $$

$$ \frac{1}{r} \left[ r (r - R_s) u_{3,2}^{D''} + (2 i k_2 r^2 + R_s) u_{3,2}^{D'} - 12 u_{3,2}^D \right] + k_2^2 \left[ \left(1 - \frac{\omega}{\Omega}\right)^2 - 1 \right] r \frac{u_{3,2}^D}{\alpha^2} = 6 \frac{\sqrt{14}}{7} \epsilon_0 \omega u_{2,2}^B. \quad (79b) $$
<table><caption>**Table 6.** Normalized spin-down luminosity for the *m* = 1 rotating quadrupole in the slowly rotating background metric.</caption><thead><tr><th>R/r<sub>L</sub></th><th>0</th><th>0.1</th><th>0.2</th><th>0.3</th><th>0.4</th><th>0.5</th></tr></thead><tbody><tr><th>0.01</th><td>1.000</td><td>1.136</td><td>1.312</td><td>1.545</td><td>1.864</td><td>2.323</td></tr><tr><th>0.02</th><td>1.000</td><td>1.140</td><td>1.321</td><td>1.560</td><td>1.888</td><td>2.362</td></tr><tr><th>0.05</th><td>1.001</td><td>1.152</td><td>1.348</td><td>1.607</td><td>1.965</td><td>2.483</td></tr><tr><th>0.1</th><td>1.007</td><td>1.176</td><td>1.397</td><td>1.693</td><td>2.105</td><td>2.708</td></tr><tr><th>0.2</th><td>1.031</td><td>1.238</td><td>1.513</td><td>1.892</td><td>2.434</td><td>3.250</td></tr></tbody></table>

<figure><img src="image_9.png" /><figcaption>**Figure 7.** Spin-down luminosity of the orthogonal quadrupole for the mode *m* = 1 in Schwarzschild and slowly rotating neutron star metric approximation.</figcaption></figure>

<table><caption>**Table 7.** Normalized spin-down luminosity for the *m* = 2 rotating quadrupole in Schwarzschild background.</caption><thead><tr><th>R/r<sub>L</sub></th><th>0</th><th>0.1</th><th>0.2</th><th>0.3</th><th>0.4</th><th>0.5</th></tr></thead><tbody><tr><th>0.01</th><td>0.9998</td><td>1.006</td><td>1.012</td><td>1.0182</td><td>1.025</td><td>1.031</td></tr><tr><th>0.02</th><td>0.9994</td><td>1.012</td><td>1.024</td><td>1.0387</td><td>1.051</td><td>1.065</td></tr><tr><th>0.05</th><td>0.9966</td><td>1.028</td><td>1.061</td><td>1.0967</td><td>1.133</td><td>1.172</td></tr><tr><th>0.1</th><td>0.9866</td><td>1.048</td><td>1.117</td><td>1.1924</td><td>1.277</td><td>1.372</td></tr><tr><th>0.2</th><td>0.9468</td><td>1.064</td><td>1.205</td><td>1.3755</td><td>1.585</td><td>1.849</td></tr></tbody></table>

<table><caption>**Table 8.** Normalized spin-down luminosity for the *m* = 2 rotating quadrupole in the slowly rotating background metric.</caption><thead><tr><th>R/r<sub>L</sub></th><th>0</th><th>0.1</th><th>0.2</th><th>0.3</th><th>0.4</th><th>0.5</th></tr></thead><tbody><tr><th>0.01</th><td>0.9998</td><td>1.006</td><td>1.012</td><td>1.018</td><td>1.025</td><td>1.031</td></tr><tr><th>0.02</th><td>0.9994</td><td>1.012</td><td>1.024</td><td>1.037</td><td>1.051</td><td>1.065</td></tr><tr><th>0.05</th><td>0.9966</td><td>1.027</td><td>1.060</td><td>1.094</td><td>1.131</td><td>1.170</td></tr><tr><th>0.1</th><td>0.9866</td><td>1.047</td><td>1.114</td><td>1.188</td><td>1.269</td><td>1.361</td></tr><tr><th>0.2</th><td>0.9468</td><td>1.059</td><td>1.193</td><td>1.353</td><td>1.548</td><td>1.788</td></tr></tbody></table>

The boundary conditions apply on the electric field such that

$$ \alpha^2 \partial_r (u_{3,2}^D e^{i k_2 r_s}) = 2 \sqrt{\frac{2}{7}} \varepsilon_0 r \tilde{\omega} f_{2,1}^B \quad (80a) $$

Table 7 summarizes results for the Schwarzschild metric and Table 8 allows comparison with the slowly rotating neutron star metric. Compactness increases again the spin-down rate although it is less pronounced than for the *m* = 1 mode. The plots in Fig. 8 give a synthetic compilation of the results showing the increase in spin-down luminosity compared to flat space-time.
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Spin-down luminosity of the orthogonal rotator quadrupole for the mode $m = 2$ in Schwarzschild and slowly rotating neutron star metric approximation.</figcaption></figure>

<table><caption>**Table 9.** Normalized spin-down luminosity for the $m = 1$ rotating hexapole in Schwarzschild background.</caption><thead><tr><th>R<sub>s</sub>/R<br>R/r<sub>L</sub></th><th>0</th><th>0.1</th><th>0.2</th><th>0.3</th><th>0.4</th><th>0.5</th></tr></thead><tbody><tr><td>0.01</td><td>0.9999</td><td>1.188</td><td>1.441</td><td>1.791</td><td>2.296</td><td>3.066</td></tr><tr><td>0.02</td><td>0.9999</td><td>1.192</td><td>1.450</td><td>1.807</td><td>2.324</td><td>3.113</td></tr><tr><td>0.05</td><td>0.9999</td><td>1.203</td><td>1.477</td><td>1.860</td><td>2.415</td><td>3.268</td></tr><tr><td>0.1</td><td>0.9998</td><td>1.222</td><td>1.525</td><td>1.951</td><td>2.578</td><td>3.552</td></tr><tr><td>0.2</td><td>0.9994</td><td>1.259</td><td>1.622</td><td>2.148</td><td>2.943</td><td>4.216</td></tr></tbody></table>

## 5.3 Hexapole solution

Higher order multipoles $ℓ ≥ 2$ are treated in the same way as a dipole or quadrupole. There is no increase in the complexity of the algorithm or in the numerical computation of the Poynting flux. However, the possible parameter space for the geometric configuration of the magnetic field augments with $ℓ$ because there is an increasing number of modes $m$ such that $m ≤ ℓ$. High-order multipoles are useful to investigate small-scale structure on the neutron star surface and also to look for off-centred dipole solutions (Pétri 2016b). In order to quantitatively fix the modification brought by general relativity, we decide two give results up to the octupole.

The required equations for the $m = 1$ hexapole mode involving $u_{2,1}^D, u_{4,1}^D, u_{3,1}^B$ are

$$ \frac{1}{r} \left[ r (r - R_s) u_{2,1}^{D''} + (2 i k r^2 + R_s) u_{2,1}^{D'} - 6 u_{2,1}^D \right] + k^2 \left[ \left(1 - \frac{\omega}{\Omega}\right)^2 - 1 \right] r \frac{u_{2,1}^D}{\alpha^2} = -3 \frac{8}{\sqrt{35}} \epsilon_0 \omega u_{3,1}^B \quad (81a) $$

$$ \frac{1}{r} \left[ r (r - R_s) u_{3,1}^{B''} + (2 i k r^2 + R_s) u_{3,1}^{B'} - 12 u_{3,1}^B \right] + k^2 \left[ \left(1 - \frac{\omega}{\Omega}\right)^2 - 1 \right] r \frac{u_{3,1}^B}{\alpha^2} = -3 \mu_0 \omega \left( \frac{8}{\sqrt{35}} u_{2,1}^D - \frac{5}{\sqrt{7}} u_{4,1}^D \right) \quad (81b) $$

$$ \frac{1}{r} \left[ r (r - R_s) u_{4,1}^{D''} + (2 i k r^2 + R_s) u_{4,1}^{D'} - 20 u_{4,1}^D \right] + k^2 \left[ \left(1 - \frac{\omega}{\Omega}\right)^2 - 1 \right] r \frac{u_{4,1}^D}{\alpha^2} = 3 \frac{5}{\sqrt{7}} \epsilon_0 \omega u_{3,1}^B. \quad (81c) $$

The boundary conditions for the electric field are

$$ \alpha^2 \partial_t (u_{2,1}^D e^{i k r_s}) = -\frac{8}{\sqrt{35}} \epsilon_0 r \tilde{\omega} f_{3,1}^B \quad (82a) $$

$$ \alpha^2 \partial_t (u_{4,1}^D e^{i k r_s}) = +\frac{5}{\sqrt{7}} \epsilon_0 r \tilde{\omega} f_{3,1}^B. \quad (82b) $$

Applying the same procedure as before, Table 9 shows the luminosity we obtained for the Schwarzschild metric. The increase in luminosity reaches values up to 4. Comparison with the slowly rotating neutron star metric is always useful and given in Table 10. The plots in Fig. 9 give a synthetic compilation of the results showing the increase in spin-down luminosity compared to flat space-time.

If we consider the $m = 2$ hexapole mode, the appropriate functions are $u_{2,2}^D, u_{4,2}^D, u_{3,2}^B$ that satisfy

$$ \frac{1}{r} \left[ r (r - R_s) u_{2,2}^{D''} + (2 i k_2 r^2 + R_s) u_{2,2}^{D'} - 6 u_{2,2}^D \right] + k_2^2 \left[ \left(1 - \frac{\omega}{\Omega}\right)^2 - 1 \right] r \frac{u_{2,2}^D}{\alpha^2} = -6 \sqrt{\frac{2}{7}} \epsilon_0 \omega u_{3,2}^B \quad (83a) $$
<table><caption>**Table 10.** Normalized spin-down luminosity for the *m* = 1 rotating hexapole in a slowly rotating metric background.</caption><thead><tr><th>R/r<sub>L</sub></th><th>0</th><th>0.1</th><th>0.2</th><th>0.3</th><th>0.4</th><th>0.5</th></tr></thead><tbody><tr><th>0.01</th><td>0.9999</td><td>1.154</td><td>1.354</td><td>1.622</td><td>1.994</td><td>2.536</td></tr><tr><th>0.02</th><td>0.9999</td><td>1.158</td><td>1.363</td><td>1.637</td><td>2.018</td><td>2.574</td></tr><tr><th>0.05</th><td>0.9999</td><td>1.168</td><td>1.388</td><td>1.684</td><td>2.096</td><td>2.701</td></tr><tr><th>0.1</th><td>0.9998</td><td>1.186</td><td>1.432</td><td>1.765</td><td>2.234</td><td>2.929</td></tr><tr><th>0.2</th><td>0.9994</td><td>1.221</td><td>1.520</td><td>1.935</td><td>2.534</td><td>3.445</td></tr></tbody></table>

<figure><img src="image_11.png" /><figcaption>**Figure 9.** Spin-down luminosity of the orthogonal rotator hexapole for the mode *m* = 1 in Schwarzschild and slowly rotating neutron star metric approximation.</figcaption></figure>

<table><caption>**Table 11.** Normalized spin-down luminosity for the *m* = 2 rotating hexapole in Schwarzschild background.</caption><thead><tr><th>R/r<sub>L</sub></th><th>0</th><th>0.1</th><th>0.2</th><th>0.3</th><th>0.4</th><th>0.5</th></tr></thead><tbody><tr><th>0.01</th><td>0.9999</td><td>1.136</td><td>1.317</td><td>1.568</td><td>1.931</td><td>2.483</td></tr><tr><th>0.02</th><td>0.9998</td><td>1.143</td><td>1.334</td><td>1.598</td><td>1.980</td><td>2.564</td></tr><tr><th>0.05</th><td>0.9992</td><td>1.163</td><td>1.385</td><td>1.692</td><td>2.140</td><td>2.831</td></tr><tr><th>0.1</th><td>0.9972</td><td>1.197</td><td>1.471</td><td>1.861</td><td>2.440</td><td>3.356</td></tr><tr><th>0.2</th><td>0.9899</td><td>1.261</td><td>1.653</td><td>2.242</td><td>3.178</td><td>4.771</td></tr></tbody></table>

$$ \frac{1}{r} \left[ r(r - R_s) u_{3,2}^{B''} + (2i k_2 r^2 + R_s) u_{3,2}^{B'} - 12 u_{3,2}^{B} \right] + k_2^2 \left[ \left(1 - \frac{\omega}{\Omega}\right)^2 - 1 \right] r \frac{u_{3,2}^{B}}{\alpha^2} = -6 \mu_0 \omega \left( \frac{\sqrt{14}}{7} u_{2,2}^{D} - \frac{\sqrt{35}}{7} u_{4,2}^{D} \right) \quad (83b) $$

$$ \frac{1}{r} \left[ r(r - R_s) u_{4,2}^{D''} + (2i k_2 r^2 + R_s) u_{4,2}^{D'} - 20 u_{4,2}^{D} \right] + k_2^2 \left[ \left(1 - \frac{\omega}{\Omega}\right)^2 - 1 \right] r \frac{u_{4,2}^{D}}{\alpha^2} = 6 \sqrt{\frac{5}{7}} \epsilon_0 \omega u_{3,2}^{B}. \quad (83c) $$

The boundary conditions for the electric field are

$$ \alpha^2 \partial_r(u_{2,2}^D e^{i k_2 r_*}) = -2 \sqrt{\frac{2}{7}} \epsilon_0 r \tilde{\omega} f_{3,2}^B \quad (84a) $$

$$ \alpha^2 \partial_r(u_{4,2}^D e^{i k_2 r_*}) = 2 \sqrt{\frac{5}{7}} \epsilon_0 r \tilde{\omega} f_{3,2}^B. \quad (84b) $$

The accurate values of the spin-down for Schwarzschild metric are given in Table 11 and for the slowly rotating neutron star in Table 12. The plots in Fig. 10 give a synthetic compilation of the results showing the increase in spin-down luminosity compared to flat space-time.

Finally, the *m* = 3 hexapole mode requires the solution of $u_{4,3}^D, u_{3,3}^B$ satisfying

$$ \frac{1}{r} \left[ r(r - R_s) u_{3,3}^{B''} + (2i k_3 r^2 + R_s) u_{3,3}^{B'} - 12 u_{3,3}^{B} \right] + k_3^2 \left[ \left(1 - \frac{\omega}{\Omega}\right)^2 - 1 \right] r \frac{u_{3,3}^{B}}{\alpha^2} = \sqrt{15} \mu_0 \omega u_{4,3}^{D} \quad (85a) $$
<table><caption>**Table 12.** Normalized spin-down luminosity for the *m* = 2 rotating hexapole in a slowly rotating metric background.</caption><thead><tr><th>R/r<sub>L</sub> / R<sub>s</sub>/R</th><th>0</th><th>0.1</th><th>0.2</th><th>0.3</th><th>0.4</th><th>0.5</th></tr></thead><tbody><tr><th>0.01</th><td>0.9999</td><td>1.112</td><td>1.257</td><td>1.450</td><td>1.717</td><td>2.107</td></tr><tr><th>0.02</th><td>0.9998</td><td>1.118</td><td>1.272</td><td>1.477</td><td>1.761</td><td>2.176</td></tr><tr><th>0.05</th><td>0.9992</td><td>1.139</td><td>1.320</td><td>1.562</td><td>1.900</td><td>2.396</td></tr><tr><th>0.1</th><td>0.9972</td><td>1.170</td><td>1.399</td><td>1.711</td><td>2.154</td><td>2.816</td></tr><tr><th>0.2</th><td>0.9899</td><td>1.228</td><td>1.557</td><td>2.029</td><td>2.737</td><td>3.863</td></tr></tbody></table>

<figure><img src="image_12.png" /><figcaption>**Figure 10.** Spin-down luminosity of the orthogonal rotator hexapole for the mode *m* = 2 in Schwarzschild and slowly rotating neutron star metric approximation.</figcaption></figure>

<table><caption>**Table 13.** Normalized spin-down luminosity for the *m* = 3 rotating hexapole in Schwarzschild background.</caption><thead><tr><th>R/r<sub>L</sub> / R<sub>s</sub>/R</th><th>0</th><th>0.1</th><th>0.2</th><th>0.3</th><th>0.4</th><th>0.5</th></tr></thead><tbody><tr><th>0.01</th><td>0.9998</td><td>1.009</td><td>1.018</td><td>1.028</td><td>1.038</td><td>1.048</td></tr><tr><th>0.02</th><td>0.9992</td><td>1.018</td><td>1.037</td><td>1.057</td><td>1.078</td><td>1.100</td></tr><tr><th>0.05</th><td>0.9955</td><td>1.042</td><td>1.093</td><td>1.148</td><td>1.207</td><td>1.270</td></tr><tr><th>0.1</th><td>0.9821</td><td>1.076</td><td>1.183</td><td>1.305</td><td>1.447</td><td>1.614</td></tr><tr><th>0.2</th><td>0.9296</td><td>1.108</td><td>1.335</td><td>1.630</td><td>2.020</td><td>2.551</td></tr></tbody></table>

<table><caption>**Table 14.** Normalized spin-down luminosity for the *m* = 3 rotating quadrupole in a slowly rotating metric background.</caption><thead><tr><th>R/r<sub>L</sub> / R<sub>s</sub>/R</th><th>0</th><th>0.1</th><th>0.2</th><th>0.3</th><th>0.4</th><th>0.5</th></tr></thead><tbody><tr><th>0.01</th><td>0.9998</td><td>1.009</td><td>1.018</td><td>1.028</td><td>1.038</td><td>1.048</td></tr><tr><th>0.02</th><td>0.9992</td><td>1.018</td><td>1.037</td><td>1.057</td><td>1.078</td><td>1.099</td></tr><tr><th>0.05</th><td>0.9955</td><td>1.042</td><td>1.092</td><td>1.146</td><td>1.204</td><td>1.266</td></tr><tr><th>0.1</th><td>0.9821</td><td>1.074</td><td>1.178</td><td>1.297</td><td>1.433</td><td>1.591</td></tr><tr><th>0.2</th><td>0.9296</td><td>1.100</td><td>1.313</td><td>1.586</td><td>1.940</td><td>2.409</td></tr></tbody></table>

$$ \frac{1}{r} \left[ r (r - R_s) u_{4,3}^{D''} + (2 i k_3 r^2 + R_s) u_{4,3}^{D'} - 20 u_{4,3}^{D} \right] + k_3^2 \left[ \left(1 - \frac{\omega}{\Omega}\right)^2 - 1 \right] r \frac{u_{4,3}^{D}}{\alpha^2} = 3 \sqrt{\frac{5}{3}} \epsilon_0 \omega u_{3,3}^{B} \quad (85b) $$

supplemented with the boundary conditions for the electric field as

$$ \alpha^2 \partial_t (u_{4,3}^D e^{i k_3 r_s}) = \sqrt{\frac{5}{3}} \epsilon_0 r \tilde{\omega} f_{3,3}^B. \quad (86) $$

Comparison of Table 13 with Table 14 demonstrates again the decrease of spin-down induced by the frame-dragging effect. The plots in Fig. 11 give a synthetic compilation of the results showing the increase in spin-down luminosity compared to flat space-time.
<figure><img src="image_13.png" /><figcaption>**Figure 11.** Spin-down luminosity of the orthogonal rotator hexapole for the mode *m* = 3 in Schwarzschild and slowly rotating neutron star metric approximation.</figcaption></figure>

<table><caption>**Table 15.** Normalized spin-down luminosity for the *m* = 1 rotating octupole in Schwarzschild background.</caption><thead><tr><th>R<sub>s</sub>/R<br>R/r<sub>L</sub></th><th>0</th><th>0.1</th><th>0.2</th><th>0.3</th><th>0.4</th><th>0.5</th></tr></thead><tbody><tr><th>0.01</th><td>0.9999</td><td>1.199</td><td>1.468</td><td>1.844</td><td>2.390</td><td>3.231</td></tr><tr><th>0.02</th><td>0.9999</td><td>1.204</td><td>1.481</td><td>1.867</td><td>2.429</td><td>3.295</td></tr><tr><th>0.05</th><td>0.9998</td><td>1.215</td><td>1.509</td><td>1.920</td><td>2.522</td><td>3.455</td></tr><tr><th>0.1</th><td>0.9992</td><td>1.234</td><td>1.556</td><td>2.013</td><td>2.689</td><td>3.748</td></tr><tr><th>0.2</th><td>0.9971</td><td>1.270</td><td>1.653</td><td>2.212</td><td>3.060</td><td>4.427</td></tr></tbody></table>

<table><caption>**Table 16.** Normalized spin-down luminosity for the *m* = 1 rotating octupole in a slowly rotating metric background.</caption><thead><tr><th>R<sub>s</sub>/R<br>R/r<sub>L</sub></th><th>0</th><th>0.1</th><th>0.2</th><th>0.3</th><th>0.4</th><th>0.5</th></tr></thead><tbody><tr><th>0.01</th><td>0.9999</td><td>1.161</td><td>1.373</td><td>1.656</td><td>2.053</td><td>2.636</td></tr><tr><th>0.02</th><td>0.9999</td><td>1.166</td><td>1.384</td><td>1.677</td><td>2.086</td><td>2.689</td></tr><tr><th>0.05</th><td>0.9998</td><td>1.177</td><td>1.410</td><td>1.724</td><td>2.166</td><td>2.818</td></tr><tr><th>0.1</th><td>0.9992</td><td>1.195</td><td>1.454</td><td>1.807</td><td>2.307</td><td>3.053</td></tr><tr><th>0.2</th><td>0.9971</td><td>1.229</td><td>1.542</td><td>1.980</td><td>2.614</td><td>3.585</td></tr></tbody></table>

## 5.4 Octupole solution

Eventually, the octupole can be solved following the same lines. For the *m* = 1 octupole mode, the coupled system reads for $u_{3,1}^D$, $u_{5,1}^D$, $u_{4,1}^B$

$$ \frac{1}{r} \left[ r (r - R_s) u_{3,1}^{D''} + (2 i k r^2 + R_s) u_{3,1}^{D'} - 12 u_{3,1}^{D} \right] + k^2 \left[ \left(1 - \frac{\omega}{\Omega}\right)^2 - 1 \right] r \frac{u_{3,1}^{D}}{\alpha^2} = -3 \frac{5}{\sqrt{7}} \epsilon_0 \omega u_{4,1}^{B} \quad (87a) $$

$$ \frac{1}{r} \left[ r (r - R_s) u_{4,1}^{B''} + (2 i k r^2 + R_s) u_{4,1}^{B'} - 20 u_{4,1}^{B} \right] + k^2 \left[ \left(1 - \frac{\omega}{\Omega}\right)^2 - 1 \right] r \frac{u_{4,1}^{B}}{\alpha^2} \\ = -3 \mu_0 \omega \left( \frac{5}{\sqrt{7}} f_{3,1}^{D} - \frac{8}{\sqrt{11}} f_{5,1}^{D} \right) \quad (87b) $$

$$ \frac{1}{r} \left[ r (r - R_s) u_{5,1}^{D''} + (2 i k r^2 + R_s) u_{5,1}^{D'} - 30 u_{5,1}^{D} \right] + k^2 \left[ \left(1 - \frac{\omega}{\Omega}\right)^2 - 1 \right] r \frac{u_{5,1}^{D}}{\alpha^2} = 3 \frac{8}{\sqrt{11}} \epsilon_0 \omega u_{4,1}^{B} \quad (87c) $$

with the boundary conditions for the electric field

$$ \alpha^2 \partial_r(u_{3,1}^D e^{i k r_*}) = -\frac{5}{\sqrt{7}} \epsilon_0 r \tilde{\omega} f_{4,1}^B \quad (88a) $$

$$ \alpha^2 \partial_r(u_{5,1}^D e^{i k r_*}) = +\frac{8}{\sqrt{11}} \epsilon_0 r \tilde{\omega} f_{4,1}^B \quad (88b) $$

Performing the numerical integration results are shown in Table 15 for Schwarzschild metric and in Table 16 for the slowly rotating neutron star metric. The plots in Fig. 12 give a synthetic compilation of the results showing the increase in spin-down luminosity compared to flat space-time.
<figure><img src="image_14.png" /><figcaption>**Figure 12.** Spin-down luminosity of the orthogonal rotator octopole for the mode $m = 1$ in Schwarzschild and slowly rotating neutron star metric approximation.</figcaption></figure>

<table><caption>**Table 17.** Normalized spin-down luminosity for the $m = 2$ rotating octopole in Schwarzschild background.</caption><thead><tr><th>R/r<sub>L</sub></th><th>0</th><th>0.1</th><th>0.2</th><th>0.3</th><th>0.4</th><th>0.5</th></tr></thead><tbody><tr><td>0.01</td><td>0.9999</td><td>1.171</td><td>1.403</td><td>1.727</td><td>2.197</td><td>2.923</td></tr><tr><td>0.02</td><td>0.9998</td><td>1.179</td><td>1.422</td><td>1.760</td><td>2.254</td><td>3.017</td></tr><tr><td>0.05</td><td>0.9991</td><td>1.200</td><td>1.475</td><td>1.863</td><td>2.434</td><td>3.327</td></tr><tr><td>0.1</td><td>0.9967</td><td>1.235</td><td>1.567</td><td>2.046</td><td>2.770</td><td>3.928</td></tr><tr><td>0.2</td><td>0.9870</td><td>1.299</td><td>1.757</td><td>2.458</td><td>3.583</td><td>5.516</td></tr></tbody></table>

<table><caption>**Table 18.** Normalized spin-down luminosity for the $m = 2$ rotating octopole in a slowly rotating metric background.</caption><thead><tr><th>R/r<sub>L</sub></th><th>0</th><th>0.1</th><th>0.2</th><th>0.3</th><th>0.4</th><th>0.5</th></tr></thead><tbody><tr><td>0.01</td><td>0.9999</td><td>1.140</td><td>1.322</td><td>1.568</td><td>1.911</td><td>2.416</td></tr><tr><td>0.02</td><td>0.9998</td><td>1.147</td><td>1.340</td><td>1.598</td><td>1.960</td><td>2.494</td></tr><tr><td>0.05</td><td>0.9991</td><td>1.168</td><td>1.390</td><td>1.690</td><td>2.115</td><td>2.745</td></tr><tr><td>0.1</td><td>0.9967</td><td>1.200</td><td>1.474</td><td>1.852</td><td>2.396</td><td>3.222</td></tr><tr><td>0.2</td><td>0.9870</td><td>1.259</td><td>1.642</td><td>2.199</td><td>3.048</td><td>4.416</td></tr></tbody></table>

For the $m = 2$ octopole mode, the coupled system reads for $u_{3,2}^D, u_{5,2}^D, u_{4,2}^B$

$$ \frac{1}{r} \left[ r (r - R_s) u_{3,2}^{D''} + (2 i k_2 r^2 + R_s) u_{3,2}^{D'} - 12 u_{3,2}^D \right] + k_m^2 \left[ \left(1 - \frac{\omega}{\Omega}\right)^2 - 1 \right] r \frac{u_{3,2}^D}{\alpha^2} = -6 \frac{\sqrt{35}}{7} \epsilon_0 \omega u_{4,2}^B \quad (89a) $$

$$ \frac{1}{r} \left[ r (r - R_s) u_{4,2}^{B''} + (2 i k r^2 + R_s) u_{4,2}^{B'} - 20 u_{4,2}^B \right] + k_2^2 \left[ \left(1 - \frac{\omega}{\Omega}\right)^2 - 1 \right] r \frac{u_{4,2}^B}{\alpha^2} = -3 \mu_0 \omega \left( 2 \sqrt{\frac{5}{7}} u_{3,2}^D - 2 \sqrt{\frac{14}{11}} u_{5,2}^D \right) \quad (89b) $$

$$ \frac{1}{r} \left[ r (r - R_s) u_{5,2}^{D''} + (2 i k_2 r^2 + R_s) u_{5,2}^{D'} - 30 u_{5,2}^D \right] + k_2^2 \left[ \left(1 - \frac{\omega}{\Omega}\right)^2 - 1 \right] r \frac{u_{5,2}^D}{\alpha^2} = 6 \sqrt{\frac{14}{11}} \epsilon_0 \omega u_{4,2}^B \quad (89c) $$

with the boundary conditions

$$ \alpha^2 \partial_t (u_{3,2}^D e^{i k_2 r_*}) = -2 \sqrt{\frac{5}{7}} \epsilon_0 r \tilde{\omega} f_{4,2}^B \quad (90a) $$

$$ \alpha^2 \partial_t (u_{5,2}^D e^{i k_2 r_*}) = 2 \sqrt{\frac{14}{11}} \epsilon_0 r \tilde{\omega} f_{4,2}^B. \quad (90b) $$

Integration results are shown in Table 17 for Schwarzschild metric and in Table 18 for the slowly rotating neutron star metric. The plots in Fig. 13 give a synthetic compilation of the results showing the increase in spin-down luminosity compared to flat space-time.
<figure><img src="image_15.png" /><figcaption>**Figure 13.** Spin-down luminosity of the orthogonal rotator octopole for the mode *m* = 2 in Schwarzschild and slowly rotating neutron star metric approximation.</figcaption></figure>

<table><caption>**Table 19.** Normalized spin-down luminosity for the *m* = 3 rotating octopole in Schwarzschild background.</caption><thead><tr><th>R/r<sub>L</sub></th><th>0</th><th>0.1</th><th>0.2</th><th>0.3</th><th>0.4</th><th>0.5</th></tr></thead><tbody><tr><th>0.01</th><td>0.9999</td><td>1.1145</td><td>1.267</td><td>1.479</td><td>1.787</td><td>2.260</td></tr><tr><th>0.02</th><td>0.9996</td><td>1.124</td><td>1.291</td><td>1.522</td><td>1.857</td><td>2.372</td></tr><tr><th>0.05</th><td>0.9976</td><td>1.1547</td><td>1.364</td><td>1.656</td><td>2.084</td><td>2.749</td></tr><tr><th>0.1</th><td>0.9906</td><td>1.1996</td><td>1.488</td><td>1.901</td><td>2.525</td><td>3.532</td></tr><tr><th>0.2</th><td>0.9634</td><td>1.2745</td><td>1.738</td><td>2.465</td><td>3.677</td><td>5.865</td></tr></tbody></table>

<table><caption>**Table 20.** Normalized spin-down luminosity for the *m* = 3 rotating octopole in a slowly rotating metric background.</caption><thead><tr><th>R/r<sub>L</sub></th><th>0</th><th>0.1</th><th>0.2</th><th>0.3</th><th>0.4</th><th>0.5</th></tr></thead><tbody><tr><th>0.01</th><td>0.9999</td><td>1.094</td><td>1.216</td><td>1.378</td><td>1.604</td><td>1.935</td></tr><tr><th>0.02</th><td>0.9996</td><td>1.104</td><td>1.239</td><td>1.418</td><td>1.667</td><td>2.030</td></tr><tr><th>0.05</th><td>0.9976</td><td>1.132</td><td>1.307</td><td>1.540</td><td>1.866</td><td>2.344</td></tr><tr><th>0.1</th><td>0.9906</td><td>1.175</td><td>1.421</td><td>1.758</td><td>2.241</td><td>2.973</td></tr><tr><th>0.2</th><td>0.9634</td><td>1.241</td><td>1.638</td><td>2.227</td><td>3.148</td><td>4.682</td></tr></tbody></table>

For the *m* = 3 octopole mode, the coupled system reads for $u_{3,3}^D$, $u_{5,3}^D$, $u_{4,3}^B$

$$ \frac{1}{r} \left[ r(r - R_s) u_{3,3}^{D''} + (2i k_3 r^2 + R_s) u_{3,3}^{D'} - 12 u_{3,3}^D \right] + k_m^2 \left[ \left(1 - \frac{\omega}{\Omega}\right)^2 - 1 \right] r \frac{u_{3,3}^D}{\alpha^2} = -\sqrt{15} \epsilon_0 \omega u_{4,3}^B \quad (91a) $$

$$ \frac{1}{r} \left[ r(r - R_s) u_{4,3}^{B''} + (2i k_3 r^2 + R_s) u_{4,3}^{B'} - 20 u_{4,3}^B \right] + k_3^2 \left[ \left(1 - \frac{\omega}{\Omega}\right)^2 - 1 \right] r \frac{u_{4,3}^B}{\alpha^2} = -3 \mu_0 \omega \left( \sqrt{\frac{5}{3}} u_{3,3}^D - 8 \sqrt{\frac{2}{33}} u_{5,3}^D \right) \quad (91b) $$

$$ \frac{1}{r} \left[ r(r - R_s) u_{5,3}^{D''} + (2i k_3 r^2 + R_s) u_{5,3}^{D'} - 30 u_{5,3}^D \right] + k_3^2 \left[ \left(1 - \frac{\omega}{\Omega}\right)^2 - 1 \right] r \frac{u_{5,3}^D}{\alpha^2} = 24 \sqrt{\frac{2}{33}} \epsilon_0 \omega u_{4,3}^B \quad (91c) $$

with the boundary conditions

$$ \alpha^2 \partial_r(u_{3,3}^D e^{i k_3 r_*}) = -\sqrt{\frac{5}{3}} \epsilon_0 r \tilde{\omega} f_{4,3}^B \quad (92a) $$

$$ \alpha^2 \partial_r(u_{5,3}^D e^{i k_3 r_*}) = 8 \sqrt{\frac{2}{33}} \epsilon_0 r \tilde{\omega} f_{4,3}^B. \quad (92b) $$

Integration results are shown in Table 19 for Schwarzschild metric and in Table 20 for the slowly rotating neutron star metric. The plots in Fig. 14 give a synthetic compilation of the results showing the increase in spin-down luminosity compared to flat space-time.
<figure><img src="image_16.png" /><figcaption>**Figure 14.** Spin-down luminosity of the orthogonal rotator octopole for the mode *m* = 3 in Schwarzschild and slowly rotating neutron star metric approximation.</figcaption></figure>

<table><caption>**Table 21.** Normalized spin-down luminosity for the *m* = 4 rotating octopole in Schwarzschild background.</caption><thead><tr><th>R/r<sub>L</sub> / R<sub>s</sub>/R</th><th>0</th><th>0.1</th><th>0.2</th><th>0.3</th><th>0.4</th><th>0.5</th></tr></thead><tbody><tr><th>0.01</th><td>0.9997</td><td>1.012</td><td>1.025</td><td>1.038</td><td>1.051</td><td>1.065</td></tr><tr><th>0.02</th><td>0.9990</td><td>1.024</td><td>1.050</td><td>1.077</td><td>1.106</td><td>1.135</td></tr><tr><th>0.05</th><td>0.9942</td><td>1.057</td><td>1.127</td><td>1.202</td><td>1.285</td><td>1.377</td></tr><tr><th>0.1</th><td>0.9773</td><td>1.104</td><td>1.253</td><td>1.429</td><td>1.641</td><td>1.898</td></tr><tr><th>0.2</th><td>0.9118</td><td>1.152</td><td>1.478</td><td>1.929</td><td>2.572</td><td>3.520</td></tr><tr><th>0.5</th><td>0.5417</td><td>0.891</td><td>1.548</td><td>2.887</td><td>5.906</td><td>13.75</td></tr></tbody></table>

<table><caption>**Table 22.** Normalized spin-down luminosity for the *m* = 4 rotating octopole in a slowly rotating metric background.</caption><thead><tr><th>R/r<sub>L</sub> / R<sub>s</sub>/R</th><th>0</th><th>0.1</th><th>0.2</th><th>0.3</th><th>0.4</th><th>0.5</th></tr></thead><tbody><tr><th>0.01</th><td>0.9997</td><td>1.012</td><td>1.025</td><td>1.038</td><td>1.051</td><td>1.065</td></tr><tr><th>0.02</th><td>0.9990</td><td>1.024</td><td>1.050</td><td>1.077</td><td>1.105</td><td>1.134</td></tr><tr><th>0.05</th><td>0.9942</td><td>1.057</td><td>1.125</td><td>1.199</td><td>1.281</td><td>1.370</td></tr><tr><th>0.1</th><td>0.9773</td><td>1.101</td><td>1.245</td><td>1.415</td><td>1.617</td><td>1.860</td></tr><tr><th>0.2</th><td>0.9118</td><td>1.140</td><td>1.444</td><td>1.856</td><td>2.427</td><td>3.241</td></tr></tbody></table>

For the *m* = 4 octopole mode, the coupled system reads for $u_{5,4}^D, u_{4,4}^B$

$$ \frac{1}{r} \left[ r (r - R_s) u_{4,4}^{B''} + (2 i k_4 r^2 + R_s) u_{4,4}^{B'} - 20 u_{4,4}^B \right] + k_4^2 \left[ \left(1 - \frac{\omega}{\Omega}\right)^2 - 1 \right] r \frac{u_{4,4}^B}{\alpha^2} = 6 \sqrt{\frac{6}{11}} \mu_0 \omega u_{5,4}^D \quad (93a) $$

$$ \frac{1}{r} \left[ r (r - R_s) u_{5,4}^{D''} + (2 i k_4 r^2 + R_s) u_{5,4}^{D'} - 30 u_{5,4}^D \right] + k_4^2 \left[ \left(1 - \frac{\omega}{\Omega}\right)^2 - 1 \right] r \frac{u_{5,4}^D}{\alpha^2} = 6 \sqrt{\frac{6}{11}} \varepsilon_0 \omega u_{4,4}^B \quad (93b) $$

with the boundary conditions

$$ \alpha^2 \partial_t (u_{5,4}^D e^{i k_4 r_s}) = 2 \sqrt{\frac{6}{11}} \varepsilon_0 r \tilde{\omega} f_{4,4}^B \quad (94a) $$

Integration results are shown in Table 21 for Schwarzschild metric and in Table 22 for the slowly rotating neutron star metric. The plots in Fig. 15 give a synthetic compilation of the results showing the increase in spin-down luminosity compared to flat space-time.

## 5.5 Contributions to spin-down luminosity

As we saw in the Newtonian case, the highest contribution to the luminosity does not necessarily correspond to the magnetic multipole of order ($\ell, m$) supposed to be anchored in the neutron star. We showed thanks to exact analytical solutions expressed in terms of spherical Hankel functions that for $\ell > m$, it is always the electric multipole of order ($\ell - 1, m$) that contributes most to the energy loss. This trend should also be observed in the general-relativistic case. Thus, for some special configurations of magnetic field, rotation and gravity, we computed separately the contribution from the magnetic and the electric multipoles to check the relevance of both part to the total Poynting flux. The
<figure><img src="image_17.png" /><figcaption>**Figure 15.** Spin-down luminosity of the orthogonal rotator octopole for the mode $m = 4$ in Schwarzschild and slowly rotating neutron star metric approximation.</figcaption></figure>

<table><caption>**Table 23.** Normalized spin-down contribution from the magnetic multipole $B-(\ell, m)$ and electric multipole $E-(\ell \pm 1, m)$ for a rotating multipole of order $(\ell, m)$ in the flat, Schwarzschild (Schw) and slowly rotating (SRNS) metric. The sign of $E-(\ell \pm 1, m)$ is chosen such that $\ell \pm 1 \ge m$ for the lowest possible $\ell$. The plus sign applies only when $\ell = m$. The parameter used are $R/r_L = 0.1$ and $R_s/R = 0.5$ (if relevant).</caption><thead><tr><th rowspan="2">Multipole</th><th colspan="2">Flat</th><th colspan="2">Schw</th><th colspan="2">SRNS</th></tr><tr><th>B-(ℓ, m)</th><th>E-(ℓ ± 1, m)</th><th>B-(ℓ, m)</th><th>E-(ℓ ± 1, m)</th><th>B-(ℓ, m)</th><th>E-(ℓ ± 1, m)</th></tr></thead><tbody><tr><td>(1,1)</td><td>0.9901</td><td>1.667e-6</td><td>1.166</td><td>1.505e-6</td><td>1.163</td><td>1.874e-6</td></tr><tr><td>(2,1)</td><td>0.1557</td><td>0.8521</td><td>0.1832</td><td>3.031</td><td>0.1832</td><td>2.524</td></tr><tr><td>(2,2)</td><td>0.9866</td><td>2.026e-6</td><td>1.372</td><td>2.414e-6</td><td>1.361</td><td>2.948e-6</td></tr><tr><td>(3,1)</td><td>0.08029</td><td>0.9195</td><td>0.09436</td><td>3.457</td><td>0.09434</td><td>2.834</td></tr><tr><td>(3,2)</td><td>0.3561</td><td>0.6411</td><td>0.4942</td><td>2.862</td><td>0.492</td><td>2.324</td></tr><tr><td>(3,3)</td><td>0.9821</td><td>1.901e-6</td><td>1.614</td><td>2.827e-6</td><td>1.591</td><td>3.406e-6</td></tr><tr><td>(4,1)</td><td>0.04884</td><td>0.9504</td><td>0.05736</td><td>3.691</td><td>0.05745</td><td>2.996</td></tr><tr><td>(4,2)</td><td>0.2033</td><td>0.7933</td><td>0.2817</td><td>3.647</td><td>0.2810</td><td>2.941</td></tr><tr><td>(4,3)</td><td>0.4915</td><td>0.4990</td><td>0.8056</td><td>2.726</td><td>0.7976</td><td>2.176</td></tr><tr><td>(4,4)</td><td>0.9773</td><td>1.706e-6</td><td>1.898</td><td>3.097e-6</td><td>1.860</td><td>3.692e-6</td></tr></tbody></table>

results are summarized in the Table 23. For the magnetic dipole $(\ell, m) = (1, 1)$, the dipolar component radiates very efficiently, it represents almost 100 per cent of the energy losses, the electric quadrupolar part is irrelevant, contributing to one part to a million. Adding space curvature and frame-dragging effects does not change our conclusion. A very good estimate of the spin-down requires only the knowledge of the magnetic radiating part. For the magnetic quadrupole $(\ell, m) = (2, 1)$, the situation is quite opposite. Indeed, the quadrupolar magnetic component radiates only about 15 per cent of the total losses, it represents a negligible part of the dynamics. However, the electric dipolar part is highest, contributing to about 85 per cent. Adding space curvature and frame-dragging effects does not drastically change the magnetic contribution, but the electric part is enhanced by a factor around 3.5. In this configuration, a very good estimate of the spin-down requires only the knowledge of the electric radiating part. Frame dragging decreases slightly the flux. For the magnetic quadrupole $(\ell, m) = (2, 2)$, we retrieve the behaviour reminiscent of the magnetic dipole. Indeed the quadrupolar magnetic component radiates almost exclusively, whereas the electric hexapolar part is completely negligible. Space curvature and frame-dragging effects enhance the magnetic radiating part slightly. Going briefly through the hexapolar and octopolar magnetic fields, the previous findings apply in the same way. For the mode $\ell = m$, the magnetic part is highest and the electric part remains negligible. General relativity slightly increases the luminosity of this magnetic multipole radiation. For all other cases with $\ell > m$, the electric part radiates at the highest level and the magnetic part just adds some corrections. However, general relativity enhances the losses from the electric multipole, multipliying the rate by a factor of several units. A corollary of this effect is a drastic difference between flat space-time and Schwarzschild space-time Poynting flux for point multipoles. General-relativistic expectations do not converge to Newtonian energy losses when $R/r_L \to 0$. To my knowledge, this subtle non-cancellation has never been noticed in the literature.

## 5.6 Fitting expressions

As a summary of all the above results, we give fitting expressions for the spin-down luminosity for each multipole. Only for the modes $\ell = m$, the spin-down luminosity tends to the flat space-time value when $R/r_L \to 0$. For the other modes, it seems that the Poynting flux stay larger than its flat space-time counterpart even for a point multipole.
<table><caption>**Table 24.** Best fit of the spin-down luminosity for all multipoles, frame dragging included.</caption><thead><tr><th>Multipole (ℓ, m)</th><th>l<sub>1</sub></th><th>l<sub>2</sub></th><th>l<sub>3</sub></th><th>l<sub>4</sub></th><th>l<sub>5</sub></th><th>l<sub>6</sub></th><th>Max error</th></tr></thead><tbody><tr><td>(1,1)</td><td>-0.00114</td><td>3.206</td><td>-2.989</td><td>-0.0009536</td><td>0.1251</td><td>10.69</td><td>0.001606</td></tr><tr><td>(2,1)</td><td>0.7208</td><td>1.706</td><td>-5.975</td><td>3.627</td><td>11.08</td><td>23.95</td><td>0.05797</td></tr><tr><td>(2,2)</td><td>-0.002069</td><td>6.542</td><td>-8.615</td><td>0.005343</td><td>-0.4638</td><td>42.52</td><td>0.007425</td></tr><tr><td>(3,1)</td><td>0.7691</td><td>1.517</td><td>-5.751</td><td>4.359</td><td>12.14</td><td>28.33</td><td>0.06858</td></tr><tr><td>(3,2)</td><td>0.5325</td><td>4.368</td><td>-20.76</td><td>3.041</td><td>16.33</td><td>97.29</td><td>0.09678</td></tr><tr><td>(3,3)</td><td>-0.007845</td><td>10.3</td><td>-20.15</td><td>0.02878</td><td>-2.698</td><td>108.1</td><td>0.02187</td></tr><tr><td>(4,1)</td><td>0.7479</td><td>2.176</td><td>-8.628</td><td>4.641</td><td>15.48</td><td>18.94</td><td>0.07213</td></tr><tr><td>(4,2)</td><td>0.6308</td><td>4.201</td><td>-24.79</td><td>3.98</td><td>20.82</td><td>111.6</td><td>0.1204</td></tr><tr><td>(4,3)</td><td>0.382</td><td>8.505</td><td>-47.93</td><td>2.564</td><td>15.63</td><td>225.2</td><td>0.1508</td></tr><tr><td>(4,4)</td><td>-0.03599</td><td>15.14</td><td>-42.8</td><td>0.08634</td><td>-8.143</td><td>226.3</td><td>0.05161</td></tr></tbody></table>

The approximate formula used to numerically fit the data for the luminosity increase includes terms up to quadratic order such that the model function is

$$ \frac{L_{GR}}{L_{flat}} = S(a, \Xi) = 1 + \Xi (l_1 + l_2 a + l_3 a^2) + \Xi^2 (l_4 + l_5 a + l_6 a^2). \quad (95) $$

The choice of this formal expression is dictated by the fact that for zero compactness $\Xi = 0$ the spin-down luminosity converges to the flat space-time limit thus $\frac{L_{GR}}{L_{flat}} = 1$. We then added corrections to first and second order in compactness $\Xi$ therefore also keeping second-order corrections in the spin parameter $a$. The accuracy of these fits dependent strongly on the multipole we try to adjust. Best results are obtained for the dipole. Nevertheless, all fits give accuracy better than 20 per cent. The list of coefficients and the maximal error in the fit are summarized in Table 24. We notice that for the cases $\ell = m$, the constant $l_1$ and $l_4$ almost vanish. This means that in the limit $a \to 0$ that is for very slow rotators, the general-relativistic luminosity tends to the flat space-time limit. This conclusion does not hold for the other multipoles satisfying $\ell \neq m$. There is always a larger spin-down in general relativity compared to Newtonian gravity.

## 5.7 Braking index

As a diagnostic of multipole fields, we computed the Poynting flux of single multipoles labelled by the mode $(\ell, m)$ taking into account the finite size of the star, space curvature and frame dragging. The braking index $n$ is another interesting related quantity which describes the efficiency of electromagnetic braking by the law $\dot{\Omega} = -K \Omega^n$ where $K$ is a constant depending on boundary conditions on the neutron star. For magnetic multipolar point sources of order $\ell$, we know that $n_\ell = 2\ell + 1$ (Krolik 1991; Pétri 2015). It is well known for instance that the braking index for magnetodipole losses is equal to $n = 3$ within small corrections in power of $R/r_L$. However, in general, the braking index can differ from this fiducial value if the size of the star is taken into account. What about corrections due to general relativity? To answer this question, we compute the braking index derived from the spin-down luminosity expressed in terms of the dimensionless parameter $a = R/r_L$ according to

$$ n_{GR} = \frac{a}{L_{GR}} \frac{dL_{GR}}{da} - 1 = \frac{d \ln L_{GR}}{d \ln a} - 1. \quad (96) $$

This is the general formula to compute the braking index in any case, knowing the luminosity of the star with respect to the spin normalized by the parameter $a$. Taking into account the general-relativistic correction factor compared to flat space-time, we get

$$ n_{GR} = \frac{d \ln S}{d \ln a} + \frac{d \ln L_{flat}}{d \ln a} - 1 = n_{flat} + \frac{d \ln S}{d \ln a} \quad (97) $$

where $n_{flat}$ is the braking index obtained from flat space-time multipoles, the one computed in Pétri (2015). General-relativity clearly introduces some corrections to the standard value expected from Newtonian gravity. From the values of the fitting coefficients, numerical applications show that the variation in the braking index remains bounded to small increase not higher than 0.3–0.4. General-relativistic effects cannot account for large discrepancies between measured braking indexes of pulsars and the fiducial value of the point dipole equal to $n = 3$.

## 6 TIME-DEPENDENT SIMULATIONS

As a final consistency check of our results, we performed some time-dependent numerical simulations of the vacuum Maxwell equations in general relativity for a magnetic dipole with $(\ell, m) = (1, 1)$, a magnetic quadrupole with $(\ell, m) = (2, 1)$ and $(\ell, m) = (2, 2)$ and a magnetic hexapolar with $(\ell, m) = (3, 1)$, $(\ell, m) = (3, 2)$ and $(\ell, m) = (3, 2)$. The relevant parameters of the simulations are a spin rate of $R/r_L = 0.1$ and use of the three metrics: Minkowski ($R_s/R = 0$), Schwarzschild ($R_s/R = 0.5$) and slowly rotating neutron star ($R_s/R = 0.5$). The pseudo-spectral discontinuous Galerkin technique using the 3+1 formalism of general relativity has been exposed in depth in Pétri (2013) and is not reproduced here.
<table><caption>**Table 25.** Normalized spin-down extracted from the time-dependent numerical simulations for $R/r_L = 0.1$ and $R_s/R = 0.5$ (if applicable).</caption><thead><tr><th>Multipole</th><th>Flat</th><th>Schw</th><th>SRNS</th></tr></thead><tbody><tr><td>(1,1)</td><td>0.987</td><td>1.163</td><td>1.160</td></tr><tr><td>(2,1)</td><td>1.006</td><td>3.209</td><td>2.704</td></tr><tr><td>(2,2)</td><td>0.987</td><td>1.373</td><td>1.362</td></tr><tr><td>(3,1)</td><td>0.995</td><td>3.583</td><td>3.022</td></tr><tr><td>(3,2)</td><td>0.997</td><td>3.446</td><td>2.906</td></tr><tr><td>(3,3)</td><td>0.982</td><td>1.580</td><td>1.558</td></tr></tbody></table>

A summary of our results are given in Table 25. These values are in perfect agreement with our semi-analytical treatment resolving the stationary problem with Helmholtz equations and expansion into rational Chebyshev functions of the generalized spherical Hankel functions. Although the numerical simulations take into account many magnetic multipoles due to coupling by frame dragging, in practice only the lowest magnetic multipole is required for a good approximate estimate.

The good agreement gives us confidence about the consistency of our analytical, semi-analytical and numerical approaches to investigate some properties of general-relativistic rotating magnetic multipoles.

## 7 CONCLUSION

Computing general-relativistic extensions of the Deutsch solution including space curvature and frame-dragging effects is easily achieved numerically to very high precision thanks to appropriate change of the coordinate system and a judicious expansion of the unknown functions on to rational Chebyshev functions. These spectral methods even allow us to get reasonable approximate analytical solutions to the Helmholtz equations in the slow-rotation metric approximation. These expressions are very handy to deduce important properties of rotating multipoles in general relativity. In this spirit, we estimated the corrections to the electromagnetic spin-down luminosity as well as to the braking index of pulsars in the parameter range of compactness and periods relevant to neutron star electrodynamics. General relativity increases the spin-down luminosity with respect to Newtonian gravity although frame dragging does not significantly impact on this estimates. Also for realistic pulsar parameters, the braking index remains essentially unaffected by gravity perturbations.

For slowly rotating neutron stars, general-relativistic effects remain small. They have only little impact on the Poynting flux and the related braking index. Nevertheless, close to the neutron star, the magnetic topology and strength can significantly differ from its Newtonian counterpart, as measured locally by an observer sitting in the gravitational field. Braking indexes larger than 3 could be explained by multipoles $ℓ > 1$ generated by the rotating dipole. Gravitational wave emission is another relativistic alternative but not directly connected to the magnetic field. In the Newtonian case, rotation already deforms the shape of the star, making it oblate at the surface and modifying the boundary conditions and therefore producing similar effects. In a real neutron star, general relativity and rotation add up to change the braking index and Poynting flux. Such corrections are also left for future work.

Future time-dependent numerical simulations of rotating multipoles in general relativity should confirm our estimates given in this paper. As an example of multipolar electromagnetic field, we plan to extend our method to a rotating off-centred dipole in general relativity. Investigation of pulse profiles and related polarization characteristics will offer a valuable insight into the magnetic topology close to the neutron star surface in the strong gravity regime. Last but not least, the effect of plasma charges and currents circulating inside the magnetosphere must be carefully included into the picture of magnetic multipoles to point out the effect of frame dragging which will become sensitive to first order in spin frequency. Only numerical simulations can cope with this tremendous task of finding stationary solutions to the filled magnetosphere.

## ACKNOWLEDGEMENTS

I am grateful to the referees for their comments and suggestions to L. Rezzolla for stimulating discussions and to A. Igoshev for further references on multipolar fields. This work has been supported by the French National Research Agency (ANR) through the grant no. ANR-13-JS05-0003-01 (project EMPERE). It also benefitted from the computational facilities available at Equip@Meso from Université de Strasbourg.

## REFERENCES

Bender C., Orszag S., 1999, Advanced Mathematical Methods for Scientists and Engineers I: Asymptotic Methods and Perturbation Theory. Springer, New York

Beskin V. S., 1990, Sov. Astron. Lett., 16, 286

Bonazzola S., Mottez F., Heyvaerts J., 2015, A&A, 573, A51

Boyd J. P., 2001, Chebyshev and Fourier Spectral Methods. Springer-Verlag, Berlin

Cook G. B., Zalutskiy M., 2014, Phys. Rev. D, 90, 124021

Ginzburg V. L., Ozernoy L. M., 1964, Zh. Eksp. Teor. Fiz., 47, 1030

Gourgouliatos K. N., Cumming A., 2014, Phys. Rev. Lett., 112, 171101
Igoshev A. P., Elfritz J. G., Popov S. B., 2016, MNRAS, 462, 3689

Kearney M. W., Kegeles L. S., Cohen J. M., 1978, Ap&SS, 56, 129

Kojima Y., Hosonuma M., 2000, Phys. Rev. D, 62, 044006

Kojima Y., Matsunaga N., Okita T., 2004, MNRAS, 348, 1388

Konno K., Kojima Y., 2000, Prog. Theor. Phys., 104, 1117

Krolik J. H., 1991, ApJ, 373, L69

Leaver E. W., 1986, J. Math. Phys., 27, 1238

Marchant P., Reisenegger A., Valdivia J. A., Hoyos J. H., 2014, ApJ, 796, 94

Morozova V. S., Ahmedov B. J., Zanotti O., 2010, MNRAS, 408, 490

Muslimov A. G., Tsygan A. I., 1986, SvA, 30, 567

Muslimov A. G., Tsygan A. I., 1992, MNRAS, 255, 61

Olver F. W. J., 2010, NIST Handbook of Mathematical Functions. National Institute of Standards and Technology (U.S.), Cambridge Univ. Press, Cambridge

Paschalidis V., Shapiro S. L., 2013, Phys. Rev. D, 88

Pétri J., 2013, MNRAS, 433, 986

Pétri J., 2015, Multipolar Electromagnetic Fields around Neutron Stars. EWASS-2015, La Laguna, Ténérife

Pétri J., 2016a, MNRAS, 455, 3779

Pétri J., 2016b, MNRAS, 463, 1240

Philipp D., Perlick V., 2015, preprint (arXiv:1503.08101)

Philippov A. A., Cerutti B., Tchekhovskoy A., Spitkovsky A., 2015, ApJ, 815, L19

Rezzolla L., Ahmedov B. J., 2004, MNRAS, 352, 1161

Rezzolla L., Ahmedov B. J., 2016, MNRAS, 459, 4144

Rezzolla L., Ahmedov B. J., Miller J. C., 2001, MNRAS, 322, 723

Ronveaux A., Arscott F. M. eds, 1995, Heun's Differential Equations. Oxford Univ. Press, Oxford

Ruiz M., Paschalidis V., Shapiro S. L., 2014, Phys. Rev. D, 89, 084045

# APPENDIX A: THE FIRST GENERAL-RELATIVISTIC MULTIPOLES

## A1 The magnetic dipole $ℓ = 1$

We start our discussion with the general-relativistic magnetic dipole. Introducing the vector spherical harmonics expansion, the static aligned dipole frozen into the neutron star is conveniently written as

$$f_{1,0}^{B} = 2 \sqrt{6\pi} \frac{B R^3}{R_s^2} \left[ \frac{\ln(1-x)}{x} + 1 + \frac{x}{2} \right] \approx -2 \sqrt{\frac{2\pi}{3}} \frac{B R^3}{r^2} \left[ 1 + \frac{3}{4} \frac{R_s}{r} \right] \quad (A1)$$

(recall that $x = R_s/r$). This solution already found by Ginzburg & Ozernoy (1964) asymptotes to the flat space-time field at very large distances $r \gg R_s$. The aligned dipolar magnetic field components are given in an orthonormal basis by

$$B^i = -3 \frac{B R^3}{R_s^3} \cos \vartheta \mathcal{L}_1(x) \approx 2 \frac{B R^3}{r^3} \cos \vartheta \left[ 1 + \frac{3}{4} \frac{R_s}{r} \right] \quad (A2a)$$

$$B^{\hat{\vartheta}} = 3 \frac{B R^3}{R_s^3} \sin \vartheta \mathcal{T}_1(x) \approx \frac{B R^3}{r^3} \sin \vartheta \left[ 1 + \frac{R_s}{r} \right] \quad (A2b)$$

$$B^{\varphi} = 0 \quad (A2c)$$

where we introduced the longitudinal and transversal parts by

$$\mathcal{L}_1(x) = x(x+2) + 2 \ln(1-x) \approx -\frac{2}{3} x^3 \left[ 1 + \frac{3}{4} x + o(x^2) \right] \quad (A3a)$$

$$\mathcal{T}_1(x) = \frac{(2-x)x + 2(1-x)\ln(1-x)}{\sqrt{1-x}} \approx \frac{1}{3} x^3 \left[ 1 + x + o(x^2) \right]. \quad (A3b)$$

The static perpendicular dipole frozen into the neutron star is conveniently written with the normalization

$$f_{1,1}^{B} = -\sqrt{2} f_{1,0}^{B} \quad (A4)$$

meaning inclining the previous aligned dipole to 90° with respect to the rotation axis leading to the magnetic field components

$$B^i = -3 \frac{B R^3}{R_s^3} e^{i \varphi} \sin \vartheta \mathcal{L}_1(x) \approx 2 \frac{B R^3}{r^3} e^{i \varphi} \sin \vartheta \left[ 1 + \frac{3}{4} \frac{R_s}{r} \right] \quad (A5a)$$

$$B^{\vartheta} = -3 \frac{B R^3}{R_s^3} e^{i \varphi} \cos \vartheta \mathcal{T}_1(x) \approx -\frac{B R^3}{r^3} e^{i \varphi} \cos \vartheta \left[ 1 + \frac{R_s}{r} \right] \quad (A5b)$$

$$B^{\varphi} = -3 \frac{B R^3}{R_s^3} i e^{i \varphi} \mathcal{T}_1(x) \approx -\frac{B R^3}{r^3} i e^{i \varphi} \left[ 1 + \frac{R_s}{r} \right]. \quad (A5c)$$
Gravitational corrections to first order expressed by the coefficient $R_s/r$ are shown to better grab the increase in magnetic field components. The amplification is different depending on the component under consideration thus the field line topology is also modified with respect to a flat space-time dipole.

## A2 The magnetic quadrupole $ℓ = 2$

Let us perform the same expansion to the general-relativistic magnetic quadrupole such that it can be expressed inside the star by

$$f_{2,0}^{B} = -\frac{8}{3} \sqrt{\frac{10\pi}{3} \frac{BR_s^4}{R_s^3}} \frac{x(x(x+6)-24)+6(3x-4)\log(1-x)}{x^2} \qquad (A6a)$$

$$\approx -4 \sqrt{\frac{2\pi}{15}} \frac{BR_s^4}{r^3} \left[ 1 + \frac{4}{3} \frac{R_s}{r} \right]. \qquad (A6b)$$

The quadrupolar magnetic field components for the axisymmetric mode $m = 0$ are given in an orthonormal basis by

$$B^{\hat{\rho}} = -\frac{10}{3} \frac{BR_s^4}{R_s^4} (3 \cos 2\vartheta + 1) \mathcal{L}_2(x) \approx \frac{BR_s^4}{r^4} (3 \cos 2\vartheta + 1) \left[ 1 + \frac{4}{3} \frac{R_s}{r} \right] \qquad (A7a)$$

$$B^{\hat{\vartheta}} = 20 \frac{BR_s^4}{R_s^4} \sin 2\vartheta \mathcal{T}_2(x) \approx 2 \frac{BR_s^4}{r^4} \sin 2\vartheta \left[ 1 + \frac{3}{2} \frac{R_s}{r} \right] \qquad (A7b)$$

$$B^{\hat{\varphi}} = 0 \qquad (A7c)$$

where we introduced the longitudinal and transversal parts by

$$\mathcal{L}_2(x) = -\frac{x(x(x+6)-24)+6(3x-4)\log(1-x)}{x} = -\frac{3}{10} x^4 \left[ 1 + \frac{4}{3} x + o(x^2) \right] \qquad (A8a)$$

$$\mathcal{T}_2(x) = \frac{x((x-12)x+12)+6(x-2)(x-1)\log(1-x)}{x\sqrt{1-x}} = \frac{1}{10} x^4 \left[ 1 + \frac{3}{2} x + o(x^2) \right]. \qquad (A8b)$$

The static (2,1)-quadrupole frozen into the neutron star is conveniently written with the normalization

$$f_{2,1}^{B} = -\frac{1}{2} \sqrt{\frac{3}{2}} f_{2,0}^{B} \qquad (A9)$$

giving the $m = 1$ mode by

$$B^{\hat{\rho}} = -5 \frac{BR_s^4}{R_s^4} e^{i\varphi} \sin 2\vartheta \mathcal{L}_2(x) \approx \frac{3}{2} \frac{BR_s^4}{r^4} e^{i\varphi} \sin 2\vartheta \left[ 1 + \frac{4}{3} \frac{R_s}{r} \right] \qquad (A10a)$$

$$B^{\hat{\vartheta}} = -10 \frac{BR_s^4}{R_s^4} e^{i\varphi} \cos 2\vartheta \mathcal{T}_2(x) \approx -\frac{BR_s^4}{r^4} e^{i\varphi} \cos 2\vartheta \left[ 1 + \frac{3}{2} \frac{R_s}{r} \right] \qquad (A10b)$$

$$B^{\hat{\varphi}} = -10 \frac{BR_s^4}{R_s^4} i e^{i\varphi} \cos \vartheta \mathcal{T}_2(x) \approx -\frac{BR_s^4}{r^4} i e^{i\varphi} \cos \vartheta \left[ 1 + \frac{3}{2} \frac{R_s}{r} \right]. \qquad (A10c)$$

The static (2,2)-quadrupole frozen into the neutron star is conveniently written with the normalization

$$f_{2,2}^{B} = \sqrt{\frac{3}{2}} f_{2,0}^{B} \qquad (A11)$$

leading to the $m = 2$ mode

$$B^{\hat{\rho}} = -10 \frac{BR_s^4}{R_s^4} e^{2i\varphi} \sin^2 \vartheta \mathcal{L}_2(x) \approx 3 \frac{BR_s^4}{r^4} e^{2i\varphi} \sin^2 \vartheta \left[ 1 + \frac{4}{3} \frac{R_s}{r} \right] \qquad (A12a)$$

$$B^{\hat{\vartheta}} = -10 \frac{BR_s^4}{R_s^4} e^{2i\varphi} \sin 2\vartheta \mathcal{T}_2(x) \approx -\frac{BR_s^4}{r^4} e^{2i\varphi} \sin 2\vartheta \left[ 1 + \frac{3}{2} \frac{R_s}{r} \right] \qquad (A12b)$$

$$B^{\hat{\varphi}} = -20 \frac{BR_s^4}{R_s^4} i e^{2i\varphi} \sin \vartheta \mathcal{T}_2(x) \approx -2 \frac{BR_s^4}{r^4} i e^{2i\varphi} \sin \vartheta \left[ 1 + \frac{3}{2} \frac{R_s}{r} \right]. \qquad (A12c)$$
A3 The magnetic hexapole ℓ = 3

Next the magnetic hexapole existing inside the neutron star is written as

$$
f_{3,0}^{B} = \frac{5}{2} \sqrt{\frac{7\pi}{3} \frac{BR^5}{R_s^4}} \frac{12(6x^2 - 20x + 15) \log(1-x) + x(x(x+12)-150) + 180}{x^3} \quad (A13a)
$$

$$
\approx -2 \sqrt{\frac{\pi}{21}} \frac{BR^5}{r^4} \left[ 1 + \frac{15}{8} \frac{R_s}{r} \right]. \qquad (A13b)
$$

The hexapolar magnetic field components for the axisymmetric mode m = 0 are given in an orthonormal basis by

$$
B^r = -\frac{35B R^5}{4R_s^5} \cos \vartheta (5 \cos^2 \vartheta - 3) \mathcal{L}_3(x) \approx \frac{B R^5}{r^5} \cos \vartheta (5 \cos^2 \vartheta - 3) \left[ 1 + \frac{15}{8} \frac{R_s}{r} \right] \quad (\text{A14a})
$$

$$
B^{\theta} = \frac{105 B R^5}{16 R_s^5} (\sin \vartheta + 5 \sin 3 \vartheta) \mathcal{T}_3(x) \approx \frac{3}{16} \frac{B R^5}{r^5} \cos \vartheta (5 \cos^2 \vartheta - 3) \left[ 1 + 2 \frac{R_s}{r} \right] \quad (\text{A14b})
$$

$$
B^{\phi} = 0 \qquad (\text{A14c})
$$

where we introduced the longitudinal and transversal parts by

$$
\mathcal{L}_3(x) = \frac{(12(6x^2 - 20x + 15) \log(1-x) + x(x(x+12) - 150) + 180)}{x^2} = -\frac{4}{35} x^5 \left[ 1 + \frac{15}{8} x + o(x^2) \right] \quad (\text{A15a})
$$

$$
\mathcal{T}_3(x) = \frac{(2-x)x((x-30)x+30) + 12(1-x)((x-5)x+5)\log(1-x)}{x^2\sqrt{1-x}} = \frac{1}{35}x^5[1+2x+o(x^2)]. \quad (\text{A15b})
$$

The static (3,1)-quadrupole frozen is normalized according to

$$
f_{3,1}^{B} = - \frac{32}{\sqrt{3}} f_{3,0}^{B} \qquad (A16)
$$

such that the m = 1 mode becomes

$$
B^r = -35 \frac{BR^5}{R_s^5} e^{i\varphi} (\sin\vartheta + 5 \sin 3\vartheta) \mathcal{L}_3(x) \approx 4 \frac{BR^5}{r^5} e^{i\varphi} (\sin\vartheta + 5 \sin 3\vartheta) \left[1 + \frac{15}{8} \frac{R_s}{r}\right] \quad (\text{A17a})
$$

$$
B^{\theta} = -\frac{35B R^5}{R_s^5} e^{i\psi}(\cos\vartheta + 15\cos 3\vartheta) \mathcal{T}_3(x) \approx -\frac{B R^5}{r^5} e^{i\psi}(\cos\vartheta + 15\cos 3\vartheta) \left[1 + 2 \frac{R_s}{r}\right] \quad (\text{A17b})
$$

$$
B^{\phi} = -\frac{70B R^5}{R_s^5} i e^{i\psi}(5\cos 2\vartheta + 3) \mathcal{T}_3(x) \approx -2\frac{B R^5}{r^5} i e^{i\psi}(5\cos 2\vartheta + 3) \left[1 + 2\frac{R_s}{r}\right]. \quad (\text{A17c})
$$

The static (3,2)-quadrupole frozen is normalized according to

$$
f_{3,2}^{B} = \sqrt{\frac{2}{15}} f_{3,0}^{B} \qquad (A18)
$$

such that the m = 2 mode gives

$$
B^r = -\frac{35 B R^5}{4 R_s^5} e^{2i\psi} \sin^2 \vartheta \cos \vartheta \mathcal{L}_3(x) \approx \frac{B R^5}{r^5} e^{2i\psi} \sin^2 \vartheta \cos \vartheta \left[1 + \frac{15}{8} \frac{R_s}{r}\right] \quad (\text{A19a})
$$

$$
B^{\theta} = \frac{35 B R^5}{16 R_s^5} e^{2i\psi} (\sin\vartheta - 3\sin 3\vartheta) \mathcal{T}_3(x) \approx \frac{B R^5}{16 r^5} e^{2i\psi} (\sin\vartheta - 3\sin 3\vartheta) \left[1 + 2 \frac{R_s}{r}\right] \quad (\text{A19b})
$$

$$
B^{\phi} = -\frac{35 B R^5}{4 R_s^5} i e^{2i\psi} \sin 2\vartheta \mathcal{T}_3(x) \approx -\frac{B R^5}{4 r^5} i e^{2i\psi} \sin 2\vartheta \left[1 + 2 \frac{R_s}{r}\right]. \quad (\text{A19c})
$$

The static (3,3)-quadrupole is normalized to

$$
f_{3,3}^{B} = \frac{8}{\sqrt{5}} f_{3,0}^{B} \qquad (A20)
$$

thus we get for the m = 3 mode

$$
B^r = \frac{35 B R^5}{R_s^5} e^{3i\psi} \sin^3 \vartheta \mathcal{L}_3(x) \approx -\frac{4B R^5}{r^5} e^{3i\psi} \sin^3 \vartheta \left[1 + \frac{15}{8} \frac{R_s}{r}\right] \quad (A21a)
$$

$$
B^{\theta} = \frac{105 B R^5}{R_s^5} e^{3i\psi} \sin^2 \vartheta \cos \vartheta \mathcal{T}_3(x) \approx \frac{3B R^5}{r^5} e^{3i\psi} \sin^2 \vartheta \cos \vartheta \left[1 + 2 \frac{R_s}{r}\right] \quad (A21b)
$$

$$
B^{\phi} = \frac{105 B R^5}{R_s^5} i e^{3i\psi} \sin^2 \vartheta \mathcal{T}_3(x) \approx \frac{3B R^5}{r^5} i e^{3i\psi} \sin^2 \vartheta \left[1 + 2 \frac{R_s}{r}\right]. \quad (A21c)
$$
A4 The magnetic octopole $ℓ = 4$

To finish this discussion about the lowest order general-relativistic multipoles, we give the exact solution to the octopole field for any azimuthal mode $m$. Inside the neutron star, the magnetic field is described a priori by a general octopolar expansion such that the magnetic potential useful for any $m$ is given with appropriate normalization by

$$f_{4,0}^{B} = \frac{336}{5} \sqrt{\frac{\pi}{5}} \frac{B R^6}{R_s^5 x^4} [x(x(3x(x+20)-1570)+4620)-3360) + 60(5x(2(x-6)x+21)-56)\log(1-x)] \quad (A22a)$$

$$\approx 8 \sqrt{\frac{\pi}{5}} \frac{B R^6}{r^5} \left[ 1 + \frac{12}{5} \frac{R_s}{r} \right] \quad (A22b)$$

The octopolar magnetic field components for the axisymmetric mode $m = 0$ are given in an orthonormal basis by

$$B^{\hat{\varphi}} = -\frac{63 B R^6}{20 R_s^6} (9 + 20 \cos 2\vartheta + 35 \cos 4\vartheta) \mathcal{L}_4(x) \quad (A23a)$$

$$\approx -\frac{3}{8} \frac{B R^6}{r^6} (9 + 20 \cos 2\vartheta + 35 \cos 4\vartheta) \left[ 1 + \frac{12}{5} \frac{R_s}{r} \right] \quad (A23b)$$

$$B^{\hat{\vartheta}} = -\frac{63 B R^6}{R_s^6} (2 \sin 2\vartheta + 7 \sin 4\vartheta) \mathcal{T}_4(x) \approx -\frac{3}{2} \frac{B R^6}{r^6} (2 \sin 2\vartheta + 7 \sin 4\vartheta) \left[ 1 + \frac{5}{2} \frac{R_s}{r} \right] \quad (A23c)$$

$$B^{\psi} = 0 \quad (A23d)$$

where we introduced the longitudinal and transversal parts by

$$\begin{aligned} \mathcal{L}_4(x) &= \frac{x(3360 - x(x(3x(x+20)-1570)+4620)) - 60(5x(2(x-6)x+21)-56)\log(1-x)}{x^3} \\ &= \frac{5}{42} x^6 \left[ 1 + \frac{12}{5} x + o(x^2) \right] \end{aligned} \quad (A24a)$$

$$\begin{aligned} \mathcal{T}_4(x) &= \frac{x(x(x(3x-190)+1030)-1680)+840}{x^3 \sqrt{1-x}} + 60(x-2)(x-1)(x-7)x+7)\log(1-x) \\ &= -\frac{1}{42} x^6 \left[ 1 + \frac{5}{2} x + o(x^2) \right]. \end{aligned} \quad (A24b)$$

The static (4,1)-quadrupole is conveniently written with the normalization

$$f_{4,1}^{B} = \frac{8}{\sqrt{5}} f_{4,0}^{B} \quad (A25)$$

For the $m = 1$ mode, we get

$$B^{\hat{\varphi}} = \frac{63 B R^6}{5 R_s^6} e^{i\varphi} \sin 2\vartheta (1 + 7 \cos 2\vartheta) \mathcal{L}_4(x) \approx \frac{3 B R^6}{2 r^6} e^{i\varphi} \sin 2\vartheta (1 + 7 \cos 2\vartheta) \left[ 1 + \frac{12}{5} \frac{R_s}{r} \right] \quad (A26a)$$

$$B^{\hat{\vartheta}} = \frac{126 B R^6}{5 R_s^6} e^{i\varphi} (\cos 2\vartheta + 7 \cos 4\vartheta) \mathcal{T}_4(x) \approx -3 \frac{B R^6}{5 r^6} e^{i\varphi} (\cos 2\vartheta + 7 \cos 4\vartheta) \left[ 1 + \frac{5}{2} \frac{R_s}{r} \right] \quad (A26b)$$

$$B^{\psi} = \frac{63 B R^6}{5 R_s^6} i e^{i\psi} (9 \cos \vartheta + 7 \cos 3\vartheta) \mathcal{T}_4(x) \approx -3 \frac{B R^6}{10 r^6} i e^{i\psi} (9 \cos \vartheta + 7 \cos 3\vartheta) \left[ 1 + \frac{5}{2} \frac{R_s}{r} \right]. \quad (A26c)$$

The static (4,2)-quadrupole is conveniently written with the normalization

$$f_{4,2}^{B} = \frac{1}{\sqrt{10}} f_{4,0}^{B} \quad (A27)$$
For the $m = 2$ mode, we get

$$B^r = -\frac{63B R^6}{20R_s^6} e^{2i\varphi}(3 + 4 \cos 2\vartheta - 7 \cos 4\vartheta) \mathcal{L}_4(x) \quad (A28a)$$

$$\approx -3 \frac{B R^6}{8r^6} e^{2i\varphi}(3 + 4 \cos 2\vartheta - 7 \cos 4\vartheta) \left[ 1 + \frac{12}{5} \frac{R_s}{r} \right] \quad (A28b)$$

$$B^\vartheta = \frac{63B R^6}{5R_s^6} e^{2i\varphi}(2 \sin 2\vartheta - 7 \sin 4\vartheta) \mathcal{T}_4(x) \approx -3 \frac{B R^6}{10r^6} e^{2i\varphi}(2 \sin 2\vartheta - 7 \sin 4\vartheta) \left[ 1 + \frac{5}{2} \frac{R_s}{r} \right] \quad (A28c)$$

$$B^\phi = -\frac{63B R^6}{5R_s^6} i e^{2i\varphi}(3 \sin \vartheta + 7 \sin 3\vartheta) \mathcal{T}_4(x) \approx \frac{3B R^6}{10r^6} i e^{2i\varphi}(3 \sin \vartheta + 7 \sin 3\vartheta) \left[ 1 + \frac{5}{2} \frac{R_s}{r} \right]. \quad (A28d)$$

The static (4,3)-quadrupole is conveniently written with the normalization

$$f_{4,3}^{B} = \sqrt{\frac{5}{7}} f_{4,0}^{B} \quad (A29)$$

For the $m = 3$ mode, we get

$$B^r = \frac{63}{2} \frac{B R^6}{R_s^6} e^{3i\varphi} (2 \sin 2\vartheta - 4 \sin 4\vartheta) \mathcal{L}_4(x) \approx \frac{15}{4} \frac{B R^6}{r^6} e^{3i\varphi} (2 \sin 2\vartheta - 4 \sin 4\vartheta) \left[ 1 + \frac{12}{5} \frac{R_s}{r} \right] \quad (A30a)$$

$$B^\vartheta = 126 \frac{B R^6}{R_s^6} e^{3i\varphi} (\cos 2\vartheta - \cos 4\vartheta) \mathcal{T}_4(x) \approx -3 \frac{B R^6}{r^6} e^{3i\varphi} (\cos 2\vartheta - \cos 4\vartheta) \left[ 1 + \frac{5}{2} \frac{R_s}{r} \right] \quad (A30b)$$

$$B^\phi = 189 \frac{B R^6}{R_s^6} i e^{3i\varphi} (\cos \vartheta - \cos 3\vartheta) \mathcal{T}_4(x) \approx -\frac{9}{2} \frac{B R^6}{r^6} i e^{3i\varphi} (\cos \vartheta - \cos 3\vartheta) \left[ 1 + \frac{5}{2} \frac{R_s}{r} \right]. \quad (A30c)$$

The static (4,4)-quadrupole is conveniently written with the normalization

$$f_{4,4}^{B} = \sqrt{\frac{10}{7}} f_{4,0}^{B} \quad (A31)$$

For the $m = 4$ mode, we get

$$B^r = -126 \frac{B R^6}{R_s^6} e^{4i\varphi} \sin^4 \vartheta \mathcal{L}_4(x) \approx -15 \frac{B R^6}{r^6} e^{4i\varphi} \sin^4 \vartheta \left[ 1 + \frac{12}{5} \frac{R_s}{r} \right] \quad (A32a)$$

$$B^\vartheta = -504 \frac{B R^6}{R_s^6} e^{4i\varphi} \sin 3\vartheta \cos \vartheta \mathcal{T}_4(x) \approx 12 \frac{B R^6}{r^6} e^{4i\varphi} \sin 3\vartheta \cos \vartheta \left[ 1 + \frac{5}{2} \frac{R_s}{r} \right] \quad (A32b)$$

$$B^\phi = -504 \frac{B R^6}{R_s^6} i e^{4i\varphi} \sin^3 \vartheta \mathcal{T}_4(x) \approx 12 \frac{B R^6}{r^6} i e^{4i\varphi} \sin^3 \vartheta \left[ 1 + \frac{5}{2} \frac{R_s}{r} \right]. \quad (A32c)$$

All these expressions for the first multipoles will be very useful for forthcoming general-relativistic time-dependent simulations of Maxwell equations to look deeper into stationary solutions of the electromagnetic field in vacuum as well as in force-free magnetospheres.

This paper has been typeset from a TeX/LATEX file prepared by the author.