¹¹institutetext: Université de Strasbourg, CNRS, Observatoire astronomique de Strasbourg, UMR 7550, F-67000 Strasbourg, France.
¹¹email: [email protected]

Spheroidal magnetic stars rotating in vacuum

J. Pétri

(Received ; accepted )

Abstract

Context. Gravity shapes stars to become almost spherical because of the isotropic nature of gravitational attraction in Newton’s theory. However, several mechanisms break this isotropy like for instance their rotation generating a centrifugal force, magnetic pressure or anisotropic equations of state. The stellar surface therefore deviates slightly or significantly from a sphere depending on the strength of these anisotropic perturbations.

Aims. In this paper, we compute analytical and numerical solutions of the electromagnetic field produced by a rotating spheroidal star of oblate or prolate nature. This study is particularly relevant for millisecond pulsars for which strong deformations are produced by the rotation or a strong magnetic field, leading to indirect observational signatures of the polar cap thermal X-ray emission.

Methods. First we solve the time harmonic Maxwell equations in vacuum by using oblate and prolate spheroidal coordinates adapted to the stellar boundary conditions. The solutions are expanded in series of radial and angular spheroidal wave functions. Particular emphasize is put on the magnetic dipole radiation. Second, we compute approximate solutions by integrating numerically the time-dependent Maxwell equations in spheroidal coordinates.

Results. We show that the spin down luminosity corrections compared to a perfect sphere are to leading order given by terms involving $(a/r_{\rm L})^{2}$ and $(a/R)^{2}$ where $a$ is the stellar oblateness or prolateness, $R$ the smallest star radius and $r_{\rm L}$ the light-cylinder radius. The corresponding perturbations in the electromagnetic field are only perceptible close to the surface, deforming the polar cap rims. At large distances $r\gg a$ , the solution tends asymptotically to the perfect spherical case of a rotating dipole.

Key Words.:

Magnetic fields – Methods: analytical – Methods: numerical – Stars: general – Stars: rotation – Stars: pulsars: general

1 Introduction

Large celestial bodies are approximately spherical because of the preponderance of gravity against other internal forces and stresses. Gravitation does not favour any direction in the sky and therefore a perfect spherical shape is expected for isotropic materials. However, celestial bodies like stars and molecular clouds are often subject to rotation, producing a centrifugal force $F_{\rm cen}$ breaking the isotropy imposed by gravity $F_{\rm grav}$ . A reference direction appears along the rotation axis. The surface of the body is deformed and can be approximated for instance by an ellipsoidal shape of oblate nature. The strength of this force must be compared to gravity at the surface. A good guess for this ratio is given by the comparison between centrifugal and gravitational forces

\frac{F_{\rm cen}}{F_{\rm grav}}\approx\frac{\Omega^{2}}{\Omega_{k}^{2}}

(1)

$\Omega$ being the rotation rate of the star (assumed to be in solid body rotation) and $\Omega_{k}$ the Keplerian angular frequency at the equator. For a more quantitative description, see Chandrasekhar (1970) discussion about spheroids and ellipsoids of revolution and Horedt (2004) who presents a comprehensive analysis of rotational effects on polytropic stars. Centrifugal forces are particularly important during the violent birth of a neutron star, where they deform the proto-neutron star to an oblate shape. The study proposed in this paper will be relevant to such early infancy of a newly born neutron star.

For fast rotating stars, with rotation rate being a fraction of the equatorial Keplerian frequency, we expect its shape to deviate from a spherical body by a fraction given by eq.(1). Following Zanazzi & Lai (2015), the rotation-induced oblateness of a celestial corps of uniform density is estimated analytically through the parameter

\epsilon=\frac{15}{16\,\pi}\,\frac{\Omega^{2}}{G\,\rho}=\frac{5}{4}\,\frac{\Omega^{2}\,R^{3}}{G\,M}=\frac{5}{4}\,\frac{\Omega^{2}}{\Omega_{k}^{2}}

(2)

defined by the moment of inertia difference between equatorial and polar direction. This is in agreement with the estimate derived in eq.(1). It gives an estimate (likely an overestimate) for the oblateness ratio $a/R$ . Exact analytical solutions exist for a constant density and uniformly rotating axisymmetric ellipsoidal fluid. If the ellipticity is defined as

e=\sqrt{1-\left(\frac{R_{\rm pol}}{R_{\rm eq}}\right)^{2}}

(3)

with $R_{\rm pol}$ and $R_{\rm eq}$ being the polar and equatorial radii of the fluid, frequency $\Omega$ given by Chandrasekhar (1970) reads

\frac{\Omega^{2}}{2\,\pi\,G\,\rho}=\frac{\sqrt{1-e^{2}}}{e^{3}}\,(3-2\,e^{2})\,\arcsin e-3\,\frac{1-e^{2}}{e^{2}}\approx\frac{4}{15}\,e^{2}.

(4)

The approximation is valid for small ellipticities $e\ll 1$ . A more realistic estimate is obtained for instance for the SLy4 equation of state as found from Table 1 of Silva et al. (2021).

The body becomes oblate and impacts its immediate surrounding gravitationally but also electromagnetically if it possesses a magnetic field. We are interested in the latter possibility namely an oblate magnetic star rotating in vacuum. It launches an electromagnetic wave responsible for its magnetic braking. Strongly magnetized neutron stars are of particular interest with this respect because they spin down according to the magnetodipole losses. Millisecond pulsars are the fastest rotating neutron stars known, they become elongated along their equator because of strong centrifugal forces. It is therefore important to understand how their oblate shape perturbs the electromagnetic field and Poynting flux in conjunction with their ageing and spin evolution. So far exact analytical solutions only exist for a perfect sphere as expressed by Deutsch (1955). Extension to multipolar fields have been thoroughly investigated by Pétri (2015), see also Bonazzola et al. (2015). It is our goal to extend these important results to oblate stars by introducing a coordinate system adapted to the stellar surface shape. Consequently oblate spheroidal coordinates are best suited to achieve this goal. These coordinates are among the eleven separable systems (Morse & Feshbach, 1953), leading to well defined solutions for the Laplace and Helmholtz equations. For completeness, we will also consider prolate shapes although these cannot be produced by rotation, but for instance by magnetic pressure. Indeed, observational signatures of such prolateness is witnessed by the torque-free precession of magnetars subject to strong toroidal magnetic fields (Makishima et al., 2016, 2019).

For strongly magnetized systems, the magnetic pressure becomes comparable to the gaseous pressure and deforms its surface to an aspherical shape with no particular symmetry axis. Spheroidal geometries could therefore be used as a first step towards more general configurations. Spheroidal coordinate systems possess the interesting and important properties of fully separation of variable for the Laplace and Helmholtz operators. It is therefore possible to solve analytically time harmonic wave emission and propagation in spheroidal coordinates. Moreover for compact objects, anisotropic equation of states for nuclear matter at very high density like in neutron stars also produces non spherical astronomical objects.

From a mathematical perspective, spheroidal wave functions applied to electromagnetic theory have been extensively discussed in Li et al. (2001). Some early application to light scattering was studied by Asano & Yamamoto (1975). A different approach employing only spheroidal coordinates has been proposed by Zeppenfeld (2009).

In this paper, we compute formally exact analytical solutions to the electromagnetic wave radiation of a stationary rotating star of spheroidal shape, oblate or prolate. In section 2, we recall the useful and important properties of the curvilinear and orthogonal coordinate systems formed by oblate and prolate coordinates, with their metric and natural basis. The separation of variables is presented and the related spheroidal wave functions are introduced. Next in section 3, we compute static solutions for the multipolar magnetic field sustained by a spheroidal object. We will discuss the important question about the normalization of a spheroidal multipole with respect to a spherical multipole taking as the reference solution. Eventually, in section 4 we compute the solution for a rotating spheroid by expansion of the electromagnetic field into spheroidal wave functions. Some useful approximate analytical solutions are presented to the lowest order in oblateness or prolateness. We close our work with accurate numerical results of spheroidal stars rotating in vacuum performed by pseudo-spectral time-dependent simulations in section 5. Conclusions are drawn in section 6.

2 Spheroidal coordinates

We start with a brief overview of both spheroidal coordinate systems, namely oblate and prolate. Let us assume a star with minimal radius $R$ , corresponding to the polar radius for oblate coordinates and to the equatorial radius for prolate coordinates. An oblate shape describes well a self-gravitating gas deformed by its own rotation. For completeness, we also consider prolate shapes although these deformations are not produced by rotation. It is advantageous to adapt the curvilinear coordinate system to the boundary conditions imposed by the gas or the star.

The Cartesian coordinate system is defined by the unit orthonormal basis $(\mathbf{e}_{\rm x},\mathbf{e}_{\rm y},\mathbf{e}_{\rm z})$ with the associated coordinates $(x,y,z)$ . We define the spheroidal coordinates relying on the Cartesian correspondence.

2.1 Oblate spheroidal coordinates

We introduce the oblate spheroidal coordinate system $(\rho,\psi,\varphi)$ such that the Cartesian coordinates $(x,y,z)$ are given by


$\displaystyle x$	$\displaystyle=\eta\,\sin\psi\,\cos\varphi$	(5a)
$\displaystyle y$	$\displaystyle=\eta\,\sin\psi\,\sin\varphi$	(5b)
$\displaystyle z$	$\displaystyle=\rho\,\cos\psi$	(5c)

where we define $\eta=\sqrt{\rho^{2}+a^{2}}$ and $a$ is a real and positive parameter related to the oblateness. The equatorial radius of the surface is then $R_{\rm eq}=\sqrt{R^{2}+a^{2}}$ . The ellipticity is defined by (Shapiro & Teukolsky, 1983)

\epsilon=\frac{R_{\rm eq}-R}{(R_{\rm eq}+R)/2}.

(6)

The natural basis vectors derived from eq.(5) are expressed as


$\displaystyle\eta\,\@vec{e}_{\rho}$	$\displaystyle=\rho\,\sin\psi\,\cos\varphi\,\mathbf{e}_{\rm x}+\rho\,\sin\psi\,\sin\varphi\,\mathbf{e}_{\rm y}+\eta\,\cos\psi\,\mathbf{e}_{\rm z}$	(7a)
$\displaystyle\@vec{e}_{\psi}$	$\displaystyle=\eta\,\cos\psi\,\cos\varphi\,\mathbf{e}_{\rm x}+\eta\,\cos\psi\,\sin\varphi\,\mathbf{e}_{\rm y}-\rho\,\sin\psi\,\mathbf{e}_{\rm z}$	(7b)
$\displaystyle\@vec{e}_{\varphi}$	$\displaystyle=-\eta\,\sin\psi\,\sin\varphi\,\mathbf{e}_{\rm x}+\eta\,\sin\psi\,\cos\varphi\,\mathbf{e}_{\rm y}.$	(7c)

The position vector $\@vec{r}$ is then expanded into

\@vec{r}=\frac{\rho\,\eta^{2}}{\Delta_{o}}\,\@vec{e}_{\rho}+\frac{a^{2}\,\cos\psi\,\sin\psi}{\Delta_{o}}\,\@vec{e}_{\psi}

(8)

with $\Delta_{o}=\rho^{2}+a^{2}\,\cos^{2}\psi$ . The oblate coordinate system is orthogonal, thus leading to a diagonal metric $g$ with coefficients


$\displaystyle g_{\rho\rho}$	$\displaystyle=\frac{\rho^{2}+a^{2}\,\cos^{2}\psi}{\eta^{2}}=\frac{\Delta_{o}}{\eta^{2}}$	(9a)
$\displaystyle g_{\psi\psi}$	$\displaystyle=\rho^{2}+a^{2}\,\cos^{2}\psi=\Delta_{o}$	(9b)
$\displaystyle g_{\varphi\varphi}$	$\displaystyle=\eta^{2}\,\sin^{2}\psi.$	(9c)

Its determinant is

\gamma=\det(g)=(\rho^{2}+a^{2}\,\cos^{2}\psi)^{2}\,\sin^{2}\psi=\Delta_{o}^{2}\,\sin^{2}\psi.

(10)

For the computation of fluxes along the coordinate $\rho$ it is useful to have the surface element vector $d\@vec{\Sigma}$ defined for an orthogonal coordinate system such as the oblate one expressed in terms of the variation of the position vector $\@vec{r}$ along two directions $\@vec{e}_{1}$ and $\@vec{e}_{2}$ such that (with the indices $1$ and $2$ being $\rho$ , $\psi$ or $\varphi$ )

d\@vec{\Sigma}=d\@vec{r}_{1}\wedge d\@vec{r}_{2}=\frac{\partial\@vec{r}}{\partial x_{1}}\wedge\frac{\partial\@vec{r}}{\partial x_{2}}\,dx_{1}\,dx_{2}.

(11)

Defined by its covariant components we get

d\Sigma_{i}=\sqrt{\gamma}\,\epsilon_{i12}\,dx_{1}\,dx_{2}

(12)

$\epsilon_{ijk}$ being the spatial Levi-Civita tensor (Arfken & Weber, 2005). For instance, the radial Poynting flux across a closed surface $\Sigma$ defined by the spheroid $\rho=\rho_{0}$ follows as

L=\oiint_{\Sigma}\@vec{S}\cdot d\@vec{\Sigma}=\oiint_{\Sigma}S^{\rho}\,d\Sigma_{\rho}

(13)

with $d\Sigma_{\rho}=g_{\rho\rho}\,d\Sigma^{\rho}=\Delta_{o}\,\sin\psi\,d\psi\,d\varphi$ and $\@vec{S}$ the Poynting vector. Exactly similar reckoning is performed for prolate coordinates as shown in th next paragraph.

2.2 Prolate spheroidal coordinates

The prolate spheroidal coordinate system $(\rho,\psi,\varphi)$ is given by


$\displaystyle x$	$\displaystyle=\rho\,\sin\psi\,\cos\varphi$	(14a)
$\displaystyle y$	$\displaystyle=\rho\,\sin\psi\,\sin\varphi$	(14b)
$\displaystyle z$	$\displaystyle=\eta\,\cos\psi.$	(14c)

It has the same functional form as eq.(5) except that it inverts $\rho$ and $\eta$ . Now the polar radius of the surface delimiting the gas or the star is then $R_{\rm pol}=\sqrt{R^{2}+a^{2}}$ . The definition of the ellipticity must be changed accordingly such that

\epsilon=\frac{R_{\rm pol}-R}{(R_{\rm pol}+R)/2}.

(15)

The natural basis vectors derived from eq.(14) are expressed as


$\displaystyle\eta\,\@vec{e}_{\rho}$	$\displaystyle=\eta\,\sin\psi\,\cos\varphi\,\mathbf{e}_{\rm x}+\eta\,\sin\psi\,\sin\varphi\,\mathbf{e}_{\rm y}+\rho\,\cos\psi\,\mathbf{e}_{\rm z}$	(16a)
$\displaystyle\@vec{e}_{\psi}$	$\displaystyle=\rho\,\cos\psi\,\cos\varphi\,\mathbf{e}_{\rm x}+\rho\,\cos\psi\,\sin\varphi\,\mathbf{e}_{\rm y}-\eta\,\sin\psi\,\mathbf{e}_{\rm z}$	(16b)
$\displaystyle\mathbf{e}_{\varphi}$	$\displaystyle=-\rho\,\sin\psi\,\sin\varphi\,\mathbf{e}_{\rm x}+\rho\,\sin\psi\,\cos\varphi\,\mathbf{e}_{\rm y}.$	(16c)

The position vector $\@vec{r}$ is given by

\@vec{r}=\frac{\rho\,\eta^{2}}{\Delta_{p}}\,\@vec{e}_{\rho}-\frac{a^{2}\,\cos\psi\,\sin\psi}{\Delta_{p}}\,\@vec{e}_{\psi}

(17)

with $\Delta_{p}=\rho^{2}+a^{2}\,\sin^{2}\psi$ . The prolate coordinate system is also orthogonal and the metric given by


$\displaystyle g_{\rho\rho}$	$\displaystyle=\frac{\rho^{2}+a^{2}\,\sin^{2}\psi}{\eta^{2}}=\frac{\Delta_{p}}{\eta^{2}}$	(18a)
$\displaystyle g_{\psi\psi}$	$\displaystyle=\rho^{2}+a^{2}\,\sin^{2}\psi=\Delta_{p}$	(18b)
$\displaystyle g_{\varphi\varphi}$	$\displaystyle=\rho^{2}\,\sin^{2}\psi.$	(18c)

Its determinant is

\gamma=det(g)=\frac{\rho^{2}}{\eta^{2}}\,(\rho^{2}+a^{2}\,\sin^{2}\psi)^{2}\,\sin^{2}\psi=\frac{\rho^{2}\,\Delta_{p}^{2}}{\eta^{2}}\,\sin^{2}\psi.

(19)

The surface element on a surface $\rho=\rho_{0}$ is in contravariant components

d\Sigma^{\rho}=\eta\,\rho\,\sin\psi\,d\psi\,d\phi.

(20)

In the whole paper, we work in these coordinate systems, using the natural basis either from eq.(7) or from eq.(16) for the components of vector fields. We now discuss the central property of spheroidal coordinates leading to fully separable variables for the scalar Helmholtz equation.

2.3 Separation of oblate variables

The scalar Helmholtz equation, that reads

\Delta W+k^{2}\,W=0

(21)

where $W$ is the unknown scalar field in three dimensions, is well known to be separable in 11 coordinate systems (Morse & Feshbach, 1953). The spheroidal coordinates, being prolate or oblate, belong to these sets. Therefore eq. (21) can be separated in three functions of one independent variable each, a radial part $P(\rho)$ , an angular part $\Psi(\psi)$ and an azimuthal part $\Phi(\varphi)$ . The second order linear differential equations derived from this separation generate the radial and angular spheroidal wave functions (Olver, 2010; Abramowitz & Stegun, 1965).

Writing explicitly the solution with the ansatz $W(\rho,\psi,\varphi)=P(\rho)\,\Psi(\psi)\,\Phi(\varphi)$ , the separation of variable leads to an azimuthal dependence $\Phi(\varphi)\propto e^{i\,m\,\varphi}$ where $m$ is an integer because $\Phi(\varphi)$ must be single-valued in $\varphi\in[0,2\,\pi]$ and to two second order linear differential equations for $\rho\geq 0$ and $\psi\in[0,\pi]$ given respectively by


$\displaystyle\frac{d}{d\rho}\left[(\rho^{2}+a^{2})\,\frac{dP}{d\rho}\right]+\left[k^{2}\,(\rho^{2}+a^{2})+\frac{m^{2}\,a^{2}}{\rho^{2}+a^{2}}\right]\,P$	$\displaystyle=+\lambda\,P$	(22a)
$\displaystyle\frac{1}{\sin\psi}\,\frac{d}{d\psi}\left[\sin\psi\,\frac{d\Psi}{d\psi}\right]-\left[k^{2}\,a^{2}\,\sin^{2}\psi+\frac{m^{2}}{\sin^{2}\psi}\right]\,\Psi$	$\displaystyle=-\lambda\,\Psi$	(22b)

$\lambda$ is a separation constant, being an eigenvalue determined by the boundary conditions. By a change to a new independent variable $z$ , letting $\rho=\pm i\,a\,z$ for eq.(22a) and $z=\cos\psi$ for eq.(22b) both equations reduce to the same Sturm-Liouville problem

\frac{d}{dz}\left[(1-z^{2})\,\frac{df}{dz}\right]+\left[\lambda+\gamma^{2}\,(1-z^{2})-\frac{m^{2}}{1-z^{2}}\right]\,f=0

(23)

with $\gamma^{2}=-k^{2}\,a^{2}<0$ meaning that $\gamma$ is purely imaginary. Note that the angular and radial wave functions are defined in different intervals, $|z|<1$ and $z>0$ respectively.

The most general solutions in each direction are given by


$\displaystyle P(\xi)$	$\displaystyle=p_{1}\,{S_{\ell}^{m}}^{(1)}(i\,\xi,\gamma)+p_{2}\,{S_{\ell}^{m}}^{(2)}(i\,\xi,\gamma)$	(24a)
$\displaystyle\Psi(\psi)$	$\displaystyle=s_{1}\,\textrm{Ps}_{\ell}^{m}(\psi,\gamma^{2})+s_{2}\,\textrm{Qs}_{\ell}^{m}(\psi,\gamma^{2})$	(24b)
$\displaystyle\Phi(\varphi)$	$\displaystyle=h_{1}\cos(m\,\varphi)+h_{2}\sin(m\,\varphi)$	(24c)

with $\xi=\rho/a$ . The radial ${S_{\ell}^{m}}^{(j)}$ ( $j=1,2$ ) and angular $\textrm{Ps}_{\ell}^{m},\textrm{Qs}_{\ell}^{m}$ spheroidal functions are defined in Olver (2010) and normalized according to Meixner & Schäfke (1954). Outgoing wave solutions require to fix $P$ proportional to ${S_{\ell}^{m}}^{(3)}(i\,\xi,\gamma)$ with

{S_{\ell}^{m}}^{(3)}(i\,\xi,\gamma)={S_{\ell}^{m}}^{(1)}(i\,\xi,\gamma)+i\,{S_{\ell}^{m}}^{(2)}(i\,\xi,\gamma)

(25)

which represents a generalization of the spherical Hankel functions $h_{\ell}^{(1)}$ and $h_{\ell}^{(2)}$ . The angular regularity condition imposes $s_{2}=0$ . In complex form, a particular solution for the scalar Helmholtz equation eq.(21) with outgoing wave boundary conditions is

W={S_{\ell}^{m}}^{(3)}(i\,\xi,\gamma)\,\textrm{Ps}_{\ell}^{m}(\psi,\gamma^{2})\,e^{i\,m\,\varphi}.

(26)

The same technique is applied to prolate coordinates as shown below.

2.4 Separation of prolate variables

Following the above lines, a same procedure for prolate coordinates leads to the angular and radial equations identical to eq. (23) but with the important difference that $\gamma^{2}=k^{2}\,a^{2}>0$ , meaning that $\gamma$ is real. The solution for outgoing waves now reads

W={S_{\ell}^{m}}^{(3)}(\zeta,\gamma)\,\textrm{Ps}_{\ell}^{m}(\psi,\gamma^{2})\,e^{i\,m\,\varphi}

(27)

with $\zeta=\eta/a=\sqrt{1+\xi^{2}}$ . The arguments of ${S_{\ell}^{m}}^{(3)}$ are all real whereas for oblate coordinates they are all purely imaginary.

Rotating objects in stationary state leading to vector Helmholtz equations reducible to a set of scalar Helmholtz equations will be studied in section 4. However, first we need to set the background magnetic field of a static magnetized object. This is exposed in the next section.

3 Static spheroidal star

For non-rotating stars, Helmholtz equation reduces to Poisson equation. By introducing a magnetic scalar potential $\phi_{M}$ , the solution to the magnetic field structure in vacuum can be solved like an electrostatic problem (Jackson, 2001). We describe the procedure for any spheroidal multipole and give explicit examples of magnetic monopoles and dipoles in oblate and prolate geometries.

3.1 Oblate magnetic star

The magnetic field $\@vec{B}$ in vacuum outside a static star is curl free and divergenceless. It can be written as the gradient of a magnetic scalar potential $\phi_{M}$ such that

\@vec{B}=-\@vec{\nabla}\phi_{M}.

(28)

The condition $\@vec{\nabla}\cdot\@vec{B}=0$ implies that the scalar potential satisfies Laplace equation

\Delta\phi_{M}=0.

(29)

This equation is fully separable in oblate spheroidal coordinates. Assuming that the inner boundary is located on the surface $\rho=\rho_{0}$ , the general solution is expanded with $\xi=\rho/a$ into

\phi_{M}(\rho,\psi,\varphi)=\sum_{\ell,|m|\leq\ell}\left(a_{\ell}^{m}\,\frac{P_{\ell}^{m}(i\,\xi)}{P_{\ell}^{m}(i\,\xi_{0})}+b_{\ell}^{m}\,\frac{Q_{\ell}^{m}(i\,\xi)}{Q_{\ell}^{m}(i\,\xi_{0})}\right)\,Y_{\ell}^{m}(\psi,\varphi)

(30)

where $P_{\ell}^{m}$ and $Q_{\ell}^{m}$ are the Legendre functions of first and second kind respectively and $Y_{\ell}^{m}$ are the spherical harmonics, see appendix A. The potential is imposed at $\rho=\rho_{0}$ with $\phi_{M}(\rho_{0},\psi,\varphi)=V(\psi,\varphi)$ and must vanish at infinity at $\rho=+\infty$ . Therefore the coefficients $a_{\ell}^{m}$ vanish. Moreover, the coefficients $b_{\ell}^{m}$ are determined by the decomposition of the surface potential into spherical harmonics

V(\psi,\varphi)=\sum_{\ell,|m|\leq\ell}V_{\ell}^{m}\,Y_{\ell}^{m}(\psi,\varphi)

(31)

and then identify $b_{\ell}^{m}=V_{\ell}^{m}$ . The solution for any magnetic potential is therefore

\phi_{M}(\rho,\psi,\varphi)=\sum_{\ell,|m|\leq\ell}V_{\ell}^{m}\,\frac{Q_{\ell}^{m}(i\,\xi)}{Q_{\ell}^{m}(i\,\xi_{0})}\,Y_{\ell}^{m}(\psi,\varphi).

(32)

The magnetic field follows from Eq. (28). The physical components are

B^{\hat{i}}=-\frac{1}{\sqrt{g_{ii}}}\,\partial_{i}\phi_{M}.

(33)

When the stellar shape tends to a perfect sphere, $a$ vanishes and

\lim_{a\to 0}\frac{Q_{\ell}^{m}(i\,\xi)}{Q_{\ell}^{m}(i\,\xi_{0})}=\left(\frac{R}{r}\right)^{\ell+1}

(34)

and we retrieve the expressions for standard magnetic multipoles (Pétri, 2015). Note that in this limit $\rho_{0}=R=R_{\rm eq}$ .

3.1.1 Monopole solution

For completeness and to better understand the impact of the stellar ellipticity on the electromagnetic field, let us start by computing the monopole solution given by the numbers $(\ell,m)=(0,0)$ . Assuming a constant potential $V$ at its surface $\rho=R$ ( $R$ should not be confused with the stellar radius that depends on colatitude $\psi$ ), the magnetic potential reads

\phi_{M}(\rho,\psi,\varphi)=V_{0}^{0}\,\frac{Q_{0}^{0}(i\,\xi)}{Q_{0}^{0}(i\,\xi_{0})}\,Y_{0}^{0}(\psi,\varphi)=V\,\frac{\textrm{arccot }\xi}{\textrm{arccot }\xi_{0}}

(35)

getting rid of spherical harmonic normalization factor by setting $V_{0}^{0}=V\sqrt{4\,\pi}$ . This potential actually depends only on the coordinate $\rho$ . The magnetic field has therefore only a $\rho$ component

B_{\rho}=-\partial_{\rho}\phi_{M}=V\,\frac{a}{(\rho^{2}+a^{2})\,\textrm{arccot }\xi_{0}}.

(36)

Introducing the constant magnetic field strength $B$ at the surface $\rho=R$ by the definition $B=B_{\rho}(R)$ the magnetic surface potential reads

V=\frac{R^{2}+a^{2}}{a}\,B\,\textrm{arccot }\xi_{0}

(37)

and the magnetic field becomes

B_{\rho}=B\,\frac{R^{2}+a^{2}}{\rho^{2}+a^{2}}.

(38)

The physical component where $B$ is the magnetic field strength at the pole ( $\rho=R,\psi=0$ ) is

B_{\hat{\rho}}=B\,\frac{R^{2}+a^{2}}{\sqrt{(\rho^{2}+a^{2})\,(\rho^{2}+a^{2}\,\cos^{2}\psi)}}.

(39)

The physical component at the equator $B_{\hat{\rho}}^{\rm eq}$ is stronger because

B_{\hat{\rho}}^{\rm eq}=B\,\sqrt{1+\frac{a^{2}}{R^{2}}}.

(40)

For $a=0$ , the star becomes perfectly spherical and the field purely radial simplifies into

B_{r}=B\,\frac{R^{2}}{r^{2}}.

(41)

In case of a small deformation with $a\ll 1$ the field expands into

B_{\hat{\rho}}\approx B\,\frac{R^{2}}{\rho^{2}}\,\left[1+\frac{a^{2}}{R^{2}}\,\left(1-\frac{R^{2}}{2\,\rho^{2}}\,(1+\cos^{2}\psi)\right)\right].

(42)

We stress that the $\rho$ component does not correspond to the spherical radial component and as such the monopolar magnetic field has actually two component both contained in the meridional plane. In the spherical coordinate system $(r,\theta,\varphi)$ , $B_{\rho}$ decomposes into a $B_{r}$ and a $B_{\theta}$ component because of the natural basis expressions in eq. (7).

At large distances the field simplifies into

B_{\hat{\rho}}\approx B\,\frac{R^{2}}{\rho^{2}}\,\left[1+\frac{a^{2}}{R^{2}}\right]\approx B\,\frac{R^{2}}{r^{2}}\,\left[1+\frac{a^{2}}{R^{2}}\right]

(43)

showing that the oblateness has still an imprint far away from the surface depicted by the correcting factor $(1+a^{2}/R^{2})$ . This factor corresponds to an increase in the magnetic field strength compared to a spherical dipole with surface field $B$ as measured by a distant observer.

3.1.2 Dipole solution

The procedure to follow for the dipole or for any multipole is very similar to the monopole case. Because of the linearity of the problem, we solve separately for an aligned and an orthogonal dipole. Even in the static case, both solutions are different because the ellipsoid defines new preferred axes breaking the spherical symmetry.

For an oblique rotator, the magnetic potential with parallel $V_{1}^{0}$ and perpendicular $V_{1}^{1}$ contribution is

\phi_{M}(\rho,\psi,\varphi)=V_{1}^{0}\,\frac{Q_{1}^{0}(i\,\xi)}{Q_{1}^{0}(i\,\xi_{0})}\,Y_{1}^{0}(\psi,\varphi)+V_{1}^{1}\,\frac{Q_{1}^{1}(i\,\xi)}{Q_{1}^{1}(i\,\xi_{0})}\,Y_{1}^{1}(\psi,\varphi)

(44)

or again getting rid of spherical harmonic normalization factors with parallel $V_{\parallel}$ and perpendicular $V_{\perp}$ contributions and a possible negative sign for the magnetic field components we write

-\phi_{M}(\rho,\psi,\varphi)=V_{\parallel}\,\frac{\xi\,\textrm{arccot }\xi-1}{\xi_{0}\,\textrm{arccot }\xi_{0}-1}\,\cos\psi+\\ V_{\perp}\,\frac{(\xi^{2}+1)\,\textrm{arccot }\xi-\xi}{(\xi_{0}^{2}+1)\,\textrm{arccot }\xi_{0}-\xi_{0}}\,\sqrt{\frac{\xi_{0}^{2}+1}{\xi^{2}+1}}\sin\psi\,e^{i\,\varphi}.

(45)

The magnetic field is decomposed into an aligned rotator as


$\displaystyle B_{\rho}$	$\displaystyle=\frac{V_{\parallel}}{a}\,\frac{(\xi^{2}+1)\,\operatorname{arccot}\xi-\xi}{(\xi^{2}+1)\,(\xi_{0}\,\operatorname{arccot}\xi_{0}-1)}\,\cos\psi$	(46a)
$\displaystyle B_{\psi}$	$\displaystyle=-V_{\parallel}\,\frac{\xi\,\operatorname{arccot}\xi-1}{\xi_{0}\,\operatorname{arccot}\xi_{0}-1}\,\sin\psi$	(46b)
$\displaystyle B_{\varphi}$	$\displaystyle=0$	(46c)

or by introducing a typical surface magnetic field strength $B_{\parallel}=B_{\rho}(R,\psi=0)$


$\displaystyle B_{\rho}$	$\displaystyle=2\,B_{\parallel}\,\frac{\xi_{0}^{2}+1}{\xi^{2}+1}\,\frac{(\xi^{2}+1)\,\operatorname{arccot}\xi-\xi}{(\xi_{0}^{2}+1)\,\operatorname{arccot}\xi_{0}-\xi_{0}}\,\cos\psi$	(47a)
$\displaystyle B_{\psi}$	$\displaystyle=2\,a\,B_{\parallel}\,\frac{(\xi_{0}^{2}+1)\,(1-\xi\,\operatorname{arccot}\xi)}{(\xi_{0}^{2}+1)\,\operatorname{arccot}\xi_{0}-\xi_{0}}\,\sin\psi$	(47b)
$\displaystyle B_{\varphi}$	$\displaystyle=0$	(47c)

and an orthogonal rotator as


$\displaystyle B_{\rho}$	$\displaystyle=\frac{V_{\perp}}{a}\,\frac{(\xi^{2}+1)\,\xi\,\operatorname{arccot}\xi-\xi^{2}-2}{(\xi_{0}^{2}+1)\,\operatorname{arccot}\xi_{0}-\xi_{0}}\,\frac{\sqrt{\xi_{0}^{2}+1}}{(\xi^{2}+1)^{3/2}}\,\sin\psi\,e^{i\,\varphi}$	(48a)
$\displaystyle B_{\psi}$	$\displaystyle=V_{\perp}\,\frac{(\xi^{2}+1)\,\textrm{arccot }\xi-\xi}{(\xi_{0}^{2}+1)\,\textrm{arccot }\xi_{0}-\xi_{0}}\,\sqrt{\frac{\xi_{0}^{2}+1}{\xi^{2}+1}}\,\cos\psi\,e^{i\,\varphi}$	(48b)
$\displaystyle B_{\varphi}$	$\displaystyle=i\,V_{\perp}\,\frac{(\xi^{2}+1)\,\textrm{arccot }\xi-\xi}{(\xi_{0}^{2}+1)\,\textrm{arccot }\xi_{0}-\xi_{0}}\,\sqrt{\frac{\xi_{0}^{2}+1}{\xi^{2}+1}}\,\sin\psi\,e^{i\,\varphi}.$	(48c)

or by introducing a typical surface magnetic field strength $B_{\perp}=B_{\rho}(R,\psi=\pi/2,\varphi=0)$


$\displaystyle B_{\rho}$	$\displaystyle=2\,B_{\perp}\,\frac{\xi^{2}+2-(\xi^{2}+1)\,\xi\,\operatorname{arccot}\xi}{\xi_{0}^{2}+2-(\xi_{0}^{2}+1)\,\xi_{0}\,\operatorname{arccot}\xi_{0}}\,\left(\frac{\xi_{0}^{2}+1}{\xi^{2}+1}\right)^{3/2}\,\sin\psi\,e^{i\,\varphi}$	(49a)
$\displaystyle B_{\psi}$	$\displaystyle=2\,a\,B_{\perp}\,\frac{(\xi^{2}+1)\,\operatorname{arccot}\xi-\xi}{(\xi_{0}^{2}+1)\,\xi_{0}\,\operatorname{arccot}\xi_{0}-\xi_{0}^{2}-2}\,\frac{(\xi_{0}^{2}+1)^{3/2}}{\sqrt{\xi^{2}+1}}\,\cos\psi\,e^{i\,\varphi}$	(49b)
$\displaystyle B_{\varphi}$	$\displaystyle=2\,i\,a\,B_{\perp}\,\frac{(\xi^{2}+1)\,\operatorname{arccot}\xi-\xi}{(\xi_{0}^{2}+1)\,\xi_{0}\,\operatorname{arccot}\xi_{0}-\xi_{0}^{2}-2}\,\frac{(\xi_{0}^{2}+1)^{3/2}}{\sqrt{\xi^{2}+1}}\,\sin\psi\,e^{i\,\varphi}.$	(49c)

In order to better connect these expressions to the spherical dipole, we introduced the magnetic field strength at the north pole by respectively $2\,B_{\parallel}$ and $2\,B_{\perp}$ for the aligned and orthogonal rotator. Because the metric coefficient $g_{\rho\rho}=1$ at the pole for an oblate coordinate system, the natural component is equal to the physical component, therefore $B_{\hat{\rho}}=B_{\rho}=2\,B_{\parallel}$ for aligned and $B_{\hat{\rho}}=B_{\rho}=2\,B_{\perp}$ for perpendicular rotators.

At large distances, aligned and perpendicular components of an oblate dipole tend respectively to


$\displaystyle B_{\rho}$	$\displaystyle=\frac{4}{3}\,B_{\parallel}\,\frac{a^{3}}{\rho^{3}}\,\frac{\xi_{0}^{2}+1}{(\xi_{0}^{2}+1)\,\operatorname{arccot}\xi_{0}-\xi_{0}}\,\cos\psi$	(50a)
$\displaystyle B_{\psi}$	$\displaystyle=\frac{2}{3}\,B_{\parallel}\,\frac{a^{3}}{\rho^{2}}\,\frac{\xi_{0}^{2}+1}{(\xi_{0}^{2}+1)\,\operatorname{arccot}\xi_{0}-\xi_{0}}\,\sin\psi$	(50b)
$\displaystyle B_{\varphi}$	$\displaystyle=0$	(50c)

for the aligned part and


$\displaystyle B_{\rho}$	$\displaystyle=\frac{8}{3}\,B_{\perp}\,\frac{a^{3}}{\rho^{3}}\,\frac{(\xi_{0}^{2}+1)^{3/2}}{\xi_{0}^{2}+2-(\xi_{0}^{2}+1)\,\xi_{0}\,\operatorname{arccot}\xi_{0}}\,\sin\psi\,e^{i\,\varphi}$	(51a)
$\displaystyle B_{\psi}$	$\displaystyle=\frac{4}{3}\,B_{\perp}\,\frac{a^{3}}{\rho^{2}}\,\frac{(\xi_{0}^{2}+1)^{3/2}}{(\xi_{0}^{2}+1)\,\xi_{0}\,\operatorname{arccot}\xi_{0}-\xi_{0}^{2}-2}\,\cos\psi\,e^{i\,\varphi}$	(51b)
$\displaystyle B_{\varphi}$	$\displaystyle=\frac{4}{3}\,i\,B_{\perp}\,\frac{a^{3}}{\rho^{2}}\,\frac{(\xi_{0}^{2}+1)^{3/2}}{(\xi_{0}^{2}+1)\,\xi_{0}\,\operatorname{arccot}\xi_{0}-\xi_{0}^{2}-2}\,\sin\psi\,e^{i\,\varphi}$	(51c)

for the orthogonal part.

3.2 Prolate magnetic star

Following the same procedure as for the oblate magnetic star we find the magnetic potential $\phi_{M}$ for prolate star as

\phi_{M}(\rho,\psi,\varphi)=\sum_{\ell,|m|\leq\ell}V_{\ell}^{m}\,\frac{Q_{\ell}^{m}(\zeta)}{Q_{\ell}^{m}(\zeta_{0})}\,Y_{\ell}^{m}(\psi,\varphi)

(52)

with $\zeta=\sqrt{1+(\rho/a)^{2}}$ and $\zeta_{0}=\sqrt{1+(R/a)^{2}}$ .

For small prolateness tending to a spherical shape, $a$ tends to zero and the potential simplifies to a multipole like

\lim_{a\to 0}\frac{Q_{\ell}^{m}(\zeta)}{Q_{\ell}^{m}(\zeta_{0})}=\left(\frac{R}{r}\right)^{\ell+1}

(53)

the same expression as in the previous subsection because oblate and prolate tend both to a spherical shape when $a\to 0$ .

3.2.1 Monopole solution

For the monopole, the potential reads explicitly

\phi_{M}(\rho,\psi,\varphi)=V_{0}^{0}\,\frac{Q_{0}^{0}(\zeta)}{Q_{0}^{0}(\zeta_{0})}\,Y_{0}^{0}(\psi,\varphi)=V\,\frac{\textrm{arccoth }\zeta}{\textrm{arccoth }\zeta_{0}}.

(54)

It actually depends only on the $\rho$ coordinate. The magnetic field has therefore only a $\rho$ component

B_{\rho}=-\partial_{\rho}\phi_{M}=V\,\frac{a}{\rho\,\sqrt{\rho^{2}+a^{2}}\,\textrm{arccoth }\zeta_{0}}.

(55)

At the stellar surface $\rho=R$ and $B_{\rho}=B$ thus

V=\frac{R\,\sqrt{R^{2}+a^{2}}}{a}\,B\,\textrm{arccoth }\zeta_{0}

(56)

such that

B_{\rho}=B\,\frac{R}{\rho}\sqrt{\frac{R^{2}+a^{2}}{\rho^{2}+a^{2}}}.

(57)

The physical component where $B$ is the magnetic field strength at the equator ( $\rho=R,\psi=\pi/2$ ) is

B_{\hat{\rho}}=B\,\frac{R}{\rho}\,\sqrt{\frac{R^{2}+a^{2}}{\rho^{2}+a^{2}\sin^{2}\psi}}.

(58)

The physical component at the pole $B_{\hat{\rho}}^{\rm pole}$ is stronger because

B_{\hat{\rho}}^{\rm pole}=B\,\sqrt{1+\frac{a^{2}}{R^{2}}}.

(59)

For small prolateness $a\ll R$ , the potential and the field becomes


$\displaystyle\phi_{M}$	$\displaystyle\approx V\,\frac{R}{\rho}\,\left[1+\frac{1}{6}\,\frac{a^{2}}{R^{2}}\,\left(1-\frac{R^{2}}{\rho^{2}}\right)\right]$	(60a)
$\displaystyle B_{\hat{\rho}}$	$\displaystyle\approx B\,\frac{R^{2}}{\rho^{2}}\,\left[1+\frac{1}{2}\,\frac{a^{2}}{R^{2}}\,\left(1-\frac{R^{2}}{\rho^{2}}\sin^{2}\psi\right)\right]$	(60b)

reducing to the spherical monopole in the limit of a perfect sphere.

At large distances, the physical component reduces to

B_{\hat{\rho}}=B\,\frac{R^{2}}{\rho^{2}}\,\sqrt{1+\frac{a^{2}}{R^{2}}}.

(61)

Here also we observe a correcting factor but of value $\sqrt{1+a^{2}/R^{2}}$ compared to a spherical monopole.

3.2.2 Dipole solution

For the prolate dipole, the potential with parallel $V_{1}^{0}$ and perpendicular $V_{1}^{1}$ contribution is

\phi_{M}(\rho,\psi,\varphi)=V_{1}^{0}\,\frac{Q_{1}^{0}(\zeta)}{Q_{1}^{0}(\zeta_{0})}\,Y_{1}^{0}(\psi,\varphi)+V_{1}^{1}\,\frac{Q_{1}^{1}(\zeta)}{Q_{1}^{1}(\zeta_{0})}\,Y_{1}^{1}(\psi,\varphi).

(62)

or explicitly with parallel $V_{\parallel}$ and perpendicular $V_{\perp}$ contributions (getting rid of spherical harmonic normalization factors)

-\phi_{M}(\rho,\psi,\varphi)=V_{\parallel}\,\frac{\zeta\,\textrm{arccoth }\zeta-1}{\zeta_{0}\,\textrm{arccoth }\zeta_{0}-1}\,\cos\psi+\\ V_{\perp}\,\frac{(\zeta^{2}-1)\,\textrm{arccoth }\zeta-\zeta}{(\zeta_{0}^{2}-1)\,\textrm{arccoth }\zeta_{0}-\zeta_{0}}\,\sqrt{\frac{\zeta_{0}^{2}-1}{\zeta^{2}-1}}\,\sin\psi\,e^{i\,\varphi}.

(63)

As for oblate stars in vacuum, because of linearity, we solve separately for an aligned and an orthogonal rotator. For the aligned rotator, the magnetic field is


$\displaystyle B_{\rho}$	$\displaystyle=-\frac{V_{\parallel}}{\rho}\left[1-\frac{\rho^{2}}{a^{2}}\,\frac{\textrm{arccoth }\zeta}{\zeta}\right]\,\frac{1}{\zeta_{0}\,\textrm{arccoth }\zeta_{0}-1}\,\cos\psi$	(64a)
$\displaystyle B_{\psi}$	$\displaystyle=V_{\parallel}\,\frac{1-\zeta\,\textrm{arccoth }\zeta}{\zeta_{0}\,\textrm{arccoth }\zeta_{0}-1}\,\sin\psi$	(64b)
$\displaystyle B_{\varphi}$	$\displaystyle=0.$	(64c)

By introducing a typical magnetic field strength at the surface, we get


$\displaystyle B_{\rho}$	$\displaystyle=2\,B_{\parallel}\,\frac{R}{\rho}\,\left[1-\xi^{2}\,\frac{\textrm{arccoth }\zeta}{\zeta}\right]\,\left[1-\xi_{0}^{2}\,\frac{\textrm{arccoth }\zeta_{0}}{\zeta_{0}}\right]^{-1}\,\cos\psi$	(65a)
$\displaystyle B_{\psi}$	$\displaystyle=2\,B_{\parallel}\,R\,\left[\zeta\,\textrm{arccoth }\zeta-1\right]\,\left[1-\xi_{0}^{2}\,\frac{\textrm{arccoth }\zeta_{0}}{\zeta_{0}}\right]^{-1}\,\sin\psi$	(65b)
$\displaystyle B_{\varphi}$	$\displaystyle=0.$	(65c)

For the orthogonal rotator, we found


$\displaystyle B_{\rho}$	$\displaystyle=\frac{V_{\perp}}{a}\,\frac{\sqrt{\zeta_{0}^{2}-1}\,((\zeta^{2}-1)\,\zeta\,\textrm{arccoth }\zeta-\zeta^{2}+2)}{((\zeta_{0}^{2}-1)\,\textrm{arccoth }\zeta_{0}-\zeta_{0})\,(\zeta^{2}-1)\,\zeta}\,\sin\psi\,e^{i\,\varphi}$	(66a)
$\displaystyle B_{\psi}$	$\displaystyle=V_{\perp}\,\frac{(\zeta^{2}-1)\,\textrm{arccoth }\zeta-\zeta}{(\zeta_{0}^{2}-1)\,\textrm{arccoth }\zeta_{0}-\zeta_{0}}\,\sqrt{\frac{\zeta_{0}^{2}-1}{\zeta^{2}-1}}\,\cos\psi\,e^{i\,\varphi}$	(66b)
$\displaystyle B_{\varphi}$	$\displaystyle=i\,V_{\perp}\,\frac{(\zeta^{2}-1)\,\textrm{arccoth }\zeta-\zeta}{(\zeta_{0}^{2}-1)\,\textrm{arccoth }\zeta_{0}-\zeta_{0}}\,\sqrt{\frac{\zeta_{0}^{2}-1}{\zeta^{2}-1}}\sin\psi\,e^{i\,\varphi}$	(66c)

or with the typical surface magnetic field strength $B_{\perp}$


$\displaystyle B_{\rho}$	$\displaystyle=2\,B_{\perp}\,\frac{\zeta_{0}}{\zeta}\,\frac{\zeta_{0}^{2}-1}{\zeta^{2}-1}\,\frac{(\zeta^{2}-1)\,\zeta\,\textrm{arccoth }\zeta-\zeta^{2}+2}{(\zeta_{0}^{2}-1)\,\zeta_{0}\,\textrm{arccoth }\zeta_{0}-\zeta_{0}^{2}+2}\,\sin\psi\,e^{i\,\varphi}$	(67a)
$\displaystyle B_{\psi}$	$\displaystyle=2\,a\,B_{\perp}\,\frac{\zeta_{0}\,(\zeta_{0}^{2}-1)}{\sqrt{\zeta^{2}-1}}\frac{(\zeta^{2}-1)\,\textrm{arccoth }\zeta-\zeta}{(\zeta_{0}^{2}-1)\,\zeta_{0}\,\textrm{arccoth }\zeta_{0}-\zeta_{0}^{2}+2}\,\cos\psi\,e^{i\,\varphi}$	(67b)
$\displaystyle B_{\varphi}$	$\displaystyle=2\,i\,a\,B_{\perp}\,\frac{\zeta_{0}\,(\zeta_{0}^{2}-1)}{\sqrt{\zeta^{2}-1}}\frac{(\zeta^{2}-1)\,\textrm{arccoth }\zeta-\zeta}{(\zeta_{0}^{2}-1)\,\zeta_{0}\,\textrm{arccoth }\zeta_{0}-\zeta_{0}^{2}+2}\,\sin\psi\,e^{i\,\varphi.}$	(67c)

At large distances, for the aligned component we get


$\displaystyle B_{\rho}$	$\displaystyle=\frac{4}{3}\,B_{\parallel}\,\frac{a^{3}}{\rho^{3}}\,\frac{\zeta_{0}\,\sqrt{\zeta_{0}^{2}-1}}{\zeta_{0}+(1-\zeta_{0}^{2})\,\textrm{arccoth }\zeta_{0}}\,\cos\psi$	(68a)
$\displaystyle B_{\psi}$	$\displaystyle=\frac{2}{3}\,B_{\parallel}\,\frac{a^{3}}{\rho^{2}}\,\frac{\zeta_{0}\,\sqrt{\zeta_{0}^{2}-1}}{\zeta_{0}+(1-\zeta_{0}^{2})\,\textrm{arccoth }\zeta_{0}}\,\sin\psi$	(68b)
$\displaystyle B_{\varphi}$	$\displaystyle=0.$	(68c)

and for the orthogonal component


$\displaystyle B_{\rho}$	$\displaystyle=\frac{8}{3}\,B_{\perp}\,\frac{a^{3}}{\rho^{3}}\,\frac{\zeta_{0}\,(\zeta_{0}^{2}-1)\,\sin\psi\,e^{i\,\varphi}}{(\zeta_{0}^{2}-1)\,\zeta_{0}\,\textrm{arccoth }\zeta_{0}-\zeta_{0}^{2}+2}$	(69a)
$\displaystyle B_{\psi}$	$\displaystyle=\frac{4}{3}\,B_{\perp}\,\frac{a^{3}}{\rho^{2}}\,\frac{\zeta_{0}\,(\zeta_{0}^{2}-1)\,\cos\psi\,e^{i\,\varphi}}{\zeta_{0}^{2}-2-(\zeta_{0}^{2}-1)\,\zeta_{0}\,\textrm{arccoth }\zeta_{0}}$	(69b)
$\displaystyle B_{\varphi}$	$\displaystyle=\frac{4}{3}\,i\,B_{\perp}\,\frac{a^{3}}{\rho^{2}}\,\frac{\sin\psi\,e^{i\,\varphi}}{\zeta_{0}^{2}-2-(\zeta_{0}^{2}-1)\,\zeta_{0}\,\textrm{arccoth }\zeta_{0}}.$	(69c)

This achieves the implementation of the initial background magnetic field set up to start the numerical simulations in section 5. A last important topic concerns the field normalization convention used to compare results obtained with different neutron star geometries.

3.3 Field normalization

Our main goal is to compare spheroidal stars to perfect spherical stars by computing some relevant quantities like for instance the Poynting flux. An important related question concerns the normalization of the spheroidal field compared to the corresponding spherical case. The impact of the stellar shape will drastically influence these quantities. Therefore we need a reference configuration with an appropriately chosen magnetic field strength at the surface. However there exists obviously several approaches to fix the electromagnetic field at the surface like keeping a fixed value of the magnetic field strength at the pole or at the equator when deforming the stellar surface. Nevertheless this seems not satisfactorily because the spin down luminosity for instance varies not only because of the spheroidal shape but also because of the artificial field strength variation related to the evolving boundary. This problem is reminiscent to the normalization of the magnetic dipole or multipole in a curved space time of general relativity. There the normalization at the surface is chosen in order to keep the asymptotic structure at large distances identical whatever the compactness and frame dragging effects. We decided to use the same techniques to normalize the surface spheroidal field, imposing a dipole magnetic field at large distances tending always to the same perfect spherical dipole keeping a constant asymptotic expression. Would we normalize it differently, the estimate of the spin down would change significantly. With our procedure, we expect to minimize all variations imputed to the field strength value at the surface retaining only the true effect of spheroidal shapes. This normalization must be carefully exposed in order to compare with previous results like for instance Finn & Shapiro (1990) who assumed a star keeping a constant magnetic flux and not a constant asymptotic magnetic dipole moment.

4 Rotating spheroidal stars

The quasi-stationary solution is acceptable for slowly rotating stars or when looking in regions well inside the light-cylinder. However when seeking for the field behaviour at large distances, it switches to a wave nature around the light-cylinder due to the finite speed of propagation of electromagnetic fields $\mathbf{F}$ . In this section, we solve for these fields in the whole space outside the star. Both the electric and the magnetic part satisfy the vector wave equation in vacuum given by

\frac{1}{c^{2}}\,\frac{\partial^{2}\mathbf{F}}{\partial t^{2}}-\Delta\mathbf{F}=\mathbf{0}.

(70)

For periodic motion, the time dependence becomes harmonic and the field $\mathbf{F}$ varies in time according to $\mathbf{F}\propto e^{-i\,\omega\,t}$ . Introducing the wave number $k=\omega/c$ , the vector wave equation reduces to the vector Helmholtz equation

\Delta\mathbf{F}+k^{2}\,\mathbf{F}=\mathbf{0}.

(71)

Next we show how to find exact analytical solutions in spheroidal coordinates at least formally.

4.1 Time harmonic solutions

Before treating the vector Helmholtz equation, it is useful to remind the results about the scalar Helmholtz equation (21). It can be shown that if $W$ is a solution of (21) then $\@vec{\nabla}W$ , $\@vec{\Phi}=\@vec{\nabla}\wedge(W\,\@vec{r})$ and $\@vec{\Psi}=\@vec{\nabla}\wedge\@vec{\Phi}$ are all solutions of the vector Helmholtz equation (71) (Leitner & Spence, 1950; Gumerov & Duraiswami, 2004). Moreover, $\@vec{\Phi}$ and $\@vec{\Psi}$ are solenoidal, meaning that $\@vec{\nabla}\cdot\@vec{\Phi}=\@vec{\nabla}\cdot\@vec{\Psi}=0$ and they satisfy $\@vec{\nabla}\wedge\@vec{\Psi}=k^{2}\,\@vec{\Phi}$ .

The time harmonic vacuum Maxwell equations can then be solved by expanding the electric and magnetic field into (Asano & Yamamoto, 1975)


$\displaystyle\@vec{E}$	$\displaystyle=\sum_{\ell,m}a^{\rm E}_{\ell m}\,\@vec{\Psi}_{\ell m}+b^{\rm E}_{\ell m}\,\@vec{\Phi}_{\ell m}$	(72a)
$\displaystyle c\,\@vec{B}$	$\displaystyle=\sum_{\ell,m}a^{\rm B}_{\ell m}\,\@vec{\Psi}_{\ell m}+b^{\rm B}_{\ell m}\,\@vec{\Phi}_{\ell m}$	(72b)

automatically satisfying $\@vec{\nabla}\cdot\@vec{E}=\@vec{\nabla}\cdot\@vec{B}=0$ . The numbers $\ell$ and $m$ label the multipole expansion modes similarly to the vector spherical harmonics (Pétri, 2013). The coefficients $a^{\rm E/B}_{\ell m}$ and $b^{\rm E/B}_{\ell m}$ are constants of integration depending on the boundary conditions. In order to satisfy the time-harmonic Maxwell equations with temporal behaviour $e^{-i\,\omega\,t}$ , these coefficients must be related by


$\displaystyle c\,b^{\rm E}_{\ell m}$	$\displaystyle=+i\,\omega\,a^{\rm B}_{\ell m}$	(73a)
$\displaystyle c\,b^{\rm B}_{\ell m}$	$\displaystyle=-i\,\omega\,a^{\rm E}_{\ell m}.$	(73b)

Full solutions to Maxwell equations in vacuum are therefore summarized by the expansion


$\displaystyle\@vec{E}$	$\displaystyle=\sum_{\ell,m}a^{\rm E}_{\ell m}\,\@vec{\Psi}_{\ell m}+i\,k\,a^{\rm B}_{\ell m}\,\@vec{\Phi}_{\ell m}$	(74a)
$\displaystyle c\,\@vec{B}$	$\displaystyle=\sum_{\ell,m}a^{\rm B}_{\ell m}\,\@vec{\Psi}_{\ell m}-i\,k\,a^{\rm E}_{\ell m}\,\@vec{\Phi}_{\ell m}.$	(74b)

The only independent coefficients are therefore $a^{\rm E/B}_{\ell m}$ . The central problem is to find explicit expressions for the vector $\@vec{\Psi}_{\ell m}$ and $\@vec{\Phi}_{\ell m}$ in spheroidal coordinates and to adjust the coefficients $a^{\rm E/B}_{\ell m}$ to fit the stellar boundary conditions.

4.2 Electric field

As often done for neutron star interiors, we assume a perfect conductor inside the star leading to an electric field in the observer frame given by $\@vec{E}+\@vec{v}\wedge\@vec{B}=\@vec{0}$ where

\@vec{v}=\@vec{\Omega}\wedge\@vec{r}

(75)

is the rotation speed of the star. Outside the star, the electric field is divergenceless as the magnetic field. Adapted to the spheroidal coordinates, the radial component $E^{\rho}$ is unconstrained, $E^{\varphi}=0$ and

E^{\psi}=-\sqrt{\frac{g_{\rho\rho}\,g_{\varphi\varphi}}{g_{\psi\psi}}}\,\Omega\,B^{\rho}.

(76)

These boundary conditions on the electric field completely and uniquely define the whole electromagnetic field in vacuum.

4.3 Approximate solution for oblique rotators

There are no exact analytical closed forms for the solutions presented above. However, the parameter $\gamma^{2}=\pm k^{2}\,a^{2}$ remains always less than one in modulus because the star must remain within its light cylinder thus the constrain $k\,a=a/r_{\rm L}\ll 1$ . Therefore we can expand the solution into a series in $\gamma$ that converges quickly to the exact expression. In this section we follow this path to get more insight into the impact of oblateness/prolateness on the Poynting flux. Our starting point are the expansion formulae for the angular and radial spheroidal wave functions as given by Abramowitz & Stegun (1965). We summarize the important results useful for our reckoning. The oblate geometry is tightly related to the prolate geometry and the results are often only given for prolate functions. Switching to oblate coordinates only requires a change of variables such that $\rho\rightarrow\pm i\,\rho$ and $\gamma\rightarrow\mp i\,\gamma$ . For brevity we only give results for prolate shapes.

The prolate angular spheroidal wave functions of the first kind are expanded into Legendre functions of the first kind via

\textrm{Ps}_{\ell}^{m}(\gamma,\psi)={\sum_{r=0,1}^{+\infty}}^{\prime}d_{r}^{m\ell}\,P_{m+r}^{m}(\psi)

(77)

where the prime indicates summation from 0/1 over even/odd indices when $(\ell-m)$ is even/odd. The expansion coefficients $d_{r}^{m\ell}$ are determined by solving a three-term recurrence relation with coefficients given in Abramowitz & Stegun (1965).

The prolate radial wave functions are then expanded into

S_{\ell}^{m}(\gamma,\zeta)=\frac{\left(1-\zeta^{-2}\right)^{m/2}}{{\sum_{r=0,1}^{+\infty}}^{\prime}d_{r}^{m\ell}\,\frac{(2\,m+r)!}{r!}}\,{\sum_{r=0,1}^{+\infty}}^{\prime}d_{r}^{m\ell}\,\frac{(2\,m+r)!}{r!}\,i^{r+m-\ell}\,z_{m+r}(\gamma\,\zeta)

(78)

where $z_{\ell}$ is any of the spherical Bessel function $j_{\ell},y_{\ell}$ or spherical Hankel function $h_{\ell}^{(1)}$ or $h_{\ell}^{(2)}$ . For our problem of outgoing wave solutions, we require a radial expansion into spherical Hankel function of type $h_{\ell}^{(1)}$ as for the Deutsch solution.

For analytically tractable purpose, we expand all parameters to second order in $\gamma$ . Actually, the coefficients $d_{r}^{m\ell}$ depend only on even powers of $\gamma$ therefore the expansion of any quantity will also follow an expansion in even powers of $\gamma$ . Consequently, the first correcting term for magnetic field structure, Poynting flux, electromagnetic torque and so on will depend on $\gamma^{2}$ . The key expansion coefficients are remind in appendix B.

4.4 Dipole radiation

Unfortunately the eigenvector expansion in spheroidal coordinates does not allow an identification term by term of each mode $(\ell,m)$ as would be possible in spherical coordinates. To get more analytical insight into the Poynting flux perturbed by the shape, we need to resort to a series expansion of the eigenvectors. We identify various contributions to the Poynting flux, the magnetic dipole being the dominant loss channel. For a spherical star, the rotating magnetic dipole induces a quadrupole electric field that also radiates. Nevertheless, this quadrupole brings in corrections to a point dipole that are of much higher order in spin parameter $w=\Omega\,R/c$ , at least $w^{4}$ compared to $1$ and $w^{2}$ for the dipole. This is due to the fact that the quadrupole is already the result of rotating a magnetic field, thus an $w^{2}$ strength for a $w^{2}$ spin down rate.

We expect this assertion to hold for a spheroidal star, meaning that useful and exact corrections can be found solely by computing the magnetic radiation part to second order in $w$ without contributions from the electric quadrupole. Interesting results are then derived by solving for the coefficient $a^{\rm B}_{1,1}$ only. The Poynting flux in such an approximation then behaves like

L\approx\frac{c}{2\,\mu_{0}}\,|a^{\rm B}_{1,1}|^{2}.

(79)

The dominant term in the magnetic dipole radiation is given by the model $(\ell,m)=(1,1)$ . The prolate radial wave function is therefore to second order in $\gamma$ given by

S_{1}^{1}(\gamma,\zeta)\approx\left[1-\left(\frac{2\,\gamma^{2}}{25}+\frac{a^{2}}{2\,\rho^{2}}\right)\right]\,h_{1}^{(1)}(k\,\rho)-\frac{2}{25}\,\gamma^{2}\,h_{3}^{(1)}(k\,\rho)+O(\gamma^{4})

(80)

For outgoing wave boundary conditions we have set $z_{\ell}=h_{\ell}^{(1)}$ . Following the procedure explained in detail in Pétri (2015), the spin down dependence on $a/R$ and $R/r_{\rm L}$ is approximately proportional to $|S_{1}^{1}(\gamma,\zeta)|^{-2}$ . An expansion to lowest order in these two parameters gives a correction to the luminosity as

L_{\rm prolate}\approx L_{\perp}\,\left[1-\left(\frac{R}{r_{\rm L}}\right)^{2}+\frac{27}{5}\,\left(\frac{a}{R}\right)^{2}-\frac{36}{5}\,\left(\frac{a}{r_{\rm L}}\right)^{2}\right]

(81)

where the spin down $L_{\perp}$ of the vacuum orthogonal point magnetic dipole is

L_{\perp}=\frac{8\,\pi}{3\,\mu_{0}\,c^{3}}\,\Omega^{4}\,B^{2}\,R^{6}.

(82)

The above approximation gives only some important hints about the spin down change due to the spheroidal shape. It does not take into account the normalization of the surface magnetic field strength. Similar analytical investigation can be performed for an oblate star, but contrary to vector spherical harmonics, vector spheroidal harmonics as defined in this work do not permit to impose the stellar surface boundary conditions on the electric field in a closed analytical form because even though the scalar spheroidal harmonics naturally embrace the spheroidal shape, their vector counterpart do not decompose easily into tangential and normal component on a spheroidal object. Therefore, imposing continuity of the normal component of the magnetic field and continuity of the tangential component of the electric field is a non trivial task. Nevertheless, it clearly shows that leading corrections scale as $\gamma^{2}=k^{2}\,a^{2}=a^{2}/r_{\rm L}^{2}$ and $a^{2}/R^{2}$ , thus are of second order in $a$ . We will use this results to fit the numerical simulations performed in the next section.

5 Time-dependent simulations

Analytically solving the Helmholtz equation with separate coordinates helps to get insight into the effect of oblateness or prolateness on the spin-down luminosity and magnetic field structure. However, the boundary conditions on the stellar surface cannot be imposed with a finite number of terms in angular spheroidal wave functions contrary to spherical harmonics for perfect spheres. It is therefore enlightening to compute numerical solutions by performing time-dependent simulations of Maxwell equations in vacuum by properly taking into account the boundary conditions on the surface with high accuracy to catch the effect of the surface electric field. This last section presents the results of such computations, first showing the structure of field lines, then investigating the spin-down luminosity and eventually tracing the shape of the polar caps.

5.1 Numerical setup

Maxwell equations are solved with our pseudo-spectral code developed in Pétri (2012). However in order to better resolve the inner computational domain, we map the usual Chebyshev grid to a truncated rational Chebyshev grid, increasing the resolution in the inner part with respect to the outer part (Boyd, 2001). This allows us to use a coarser grid of only $N_{r}\times N_{\theta}\times N_{\varphi}=129\times 32\times 64$ . The neutron star is a perfect conductor imposing an electric field on its surface given by Eq. (76). The neutron star radius is set to $R/r_{\rm L}=0.3$ which coincides with the inner boundary of the computational domain $R_{1}=R$ . The outer boundary is equal to $R_{2}/r_{\rm L}=7$ . The oblateness or prolateness is controlled by the parameter $a$ defining $\eta$ in spheroidal coordinates. This parameter $a/R$ spans the range $[0,1]$ although it is not bounded by $R$ but by the fact that the equatorial radius of the star cannot exceed the light-cylinder radius. The obliquity $\textstyle\chi$ is taken in the set ${\mathchoice{\raisebox{0.0pt}{$\displaystyle\chi$}}{\raisebox{0.0pt}{$\textstyle\chi$}}{\raisebox{0.0pt}{$\scriptstyle\chi$}}{\raisebox{0.0pt}{$\scriptscriptstyle\chi$}}}\in\{0\degr,30\degr,60\degr,90\degr\}$ .

5.2 Magnetic field lines

Let us first compare the impact of the spheroidal shape onto the magnetic field line structure for an aligned and a perpendicular rotator for ease to plot accurately in 2D. They are shown respectively in Fig. 1 and Fig. 2 for oblate and prolate stars with different boundary conditions, either single multipolar spheroidal field in vacuum or spherical field in vacuum. The spheroidal parameter is chosen as $a/R=\{0,0.5,1\}$ .

For the aligned rotator, Fig. 1, we observe the deformation of the surface as a change in the position of the foot-points of the magnetic field lines. Far from the star, especially outside the light-cylinder, there is hardly a hint about the nature of the spheroidal star. The impact is highest on the surface, and can be quantified by the polar cap rim change as will be shown in the corresponding subsection.

Refer to caption — Figure 1: Magnetic field lines in the meridional plane $xOz$ for an aligned spheroidal rotator with oblateness (first two columns) or prolateness (last two columns) parameter $a/R=\{0,0.5,1\}$ respectively in blue, green and red. The blue disk on the bottom left represents the spherical star.

For the orthogonal rotator, Fig. 2, magnetic field lines are shown in the equatorial plane. For prolate shapes, the stellar deformation is not seen because at the equator its size does not vary. The two-armed spiral pattern typical of Deutsch solution is preserved for spheroidal stars with slight changes.

As for an offset dipole or for a dipole plus multipole components, at large distances, outside the light-cylinder, the field reduces to the magnetic dipole in vacuum, washing the structure at the surface to keep only the leading lowest order term. In our case, the dominant and most relevant multipole component is the dipole, decreasing like $r^{-3}$ , thus even for a spheroidal star, we expect the spin down luminosity to follow expressions very similar to a spherical star with a dependence on inclination $\textstyle\chi$ like $\sin^{2}{\mathchoice{\raisebox{0.0pt}{$\displaystyle\chi$}}{\raisebox{0.0pt}{$\textstyle\chi$}}{\raisebox{0.0pt}{$\scriptstyle\chi$}}{\raisebox{0.0pt}{$\scriptscriptstyle\chi$}}}$ as will be shown in the next paragraph.

5.3 Poynting flux

The spin-down rate of a magnetic dipole is controlled by the Poynting flux. Exact analytical formulae exist for spherical stars but for spheroidal ones, we have to resort to numerical approximations. Fig. 3 shows the evolution of the spin down depending on the parameter $a$ normalized to the stellar radius $R$ for an oblate or a prolate star, in solid and dashed lines respectively. The background magnetic field is set to the static spheroidal solution presented in Section 3. The obliquity is set to ${\mathchoice{\raisebox{0.0pt}{$\displaystyle\chi$}}{\raisebox{0.0pt}{$\textstyle\chi$}}{\raisebox{0.0pt}{$\scriptstyle\chi$}}{\raisebox{0.0pt}{$\scriptscriptstyle\chi$}}}\in\{0\degr,30\degr,60\degr,90\degr\}$ as given in the legend.

Fig. 4 shows the same evolution but for a star keeping a perfect spherical dipole structure.

As a general trend, we observe a decrease in the spin-down luminosity when the normalised asphericity $a/R$ increases. To very good accuracy the vacuum spin-down can be approximated by quadratic corrections in $a/R$ such that

\frac{L}{L_{\perp}}\approx\left[k_{1}-k_{2}\,\left(\frac{a}{R}\right)^{2}\right]\,\sin^{2}{\mathchoice{\raisebox{0.0pt}{$\displaystyle\chi$}}{\raisebox{0.0pt}{$\textstyle\chi$}}{\raisebox{0.0pt}{$\scriptstyle\chi$}}{\raisebox{0.0pt}{$\scriptscriptstyle\chi$}}}.

(83)

The coefficients are listed in table 1.

Model	$k_{1}$	$k_{2}$
oblate	0.921	0.0490
prolate	0.921	0.0186
oblate spherical	0.921	0.0459
prolate spherical	0.921	0.0159

Table 1: Fitted coefficients

k_{1}

and

k_{2}

as given by eq. (83).

Actually the coefficient $k_{1}$ is known analytically and expressed in terms of spherical Hankel functions $h^{(1)}_{\ell}(R/r_{\rm L})$ , see Pétri (2015). For the particular values of our simulations, we should find approximately $k_{1}\approx 0.919$ . However, because the outgoing wave boundary conditions stand at a finite radius, relatively close to the light-cylinder, the numerical flux is impacted by these boundary surface. As carefully shown in Pétri (2014), the accuracy scales as $R_{2}^{-2}$ . However, the error remains very weak, amounting to only 0.2%.

What is the impact of this spin down on for instance stellar magnetic field inference? The accreting millisecond X-ray pulsar IGR J00291+5934 with a spin frequency of 599 Hz is a good example, assuming a mass $M=1.4\,M_{\odot}$ and a radius $R=12$ km it would have $a/R\approx 0.55$ according to the MacLaurin spheroid expression in eq. (4). This represents a large deformation of the stellar surface. The more realistic model of Silva et al. (2021) with an SLy4 equation of state would give a slightly lower value of $a/R\approx 0.48$ but still large. From equation (81) we can then find that the spheroidal formula for the luminosity gives a correction of a factor of order 2, which is significant. However from eq. (83) we expect a much weaker impact due to the removal of the $a/R$ dependence by our normalization procedure at large distances. This discrepancy has to be kept in mind because of the indeterminacy of a relevant normalization.

The spin down luminosity correction can be large depending on the choice of neutron star sequences used to compute the spheroidal electromagnetic field. Indeed, the normalization significantly affects the correcting factor. The key process is to choose a meaningful sequence by keeping some physical parameters constant while deforming the stellar surface from a sphere to a spheroid. Several choices are possible, for instance keeping the equatorial or the polar magnetic field strength constant as done in section 3. But we could keep the magnetic moment or the magnetic flux threading the star constant or the asymptotic field structure at large distance as done in the numerical simulations. Therefore the central question arise: to which spherical star should a spheroidal star be compared? There is no unique answer and the best normalization must be adapted to the problem under scrutiny. The spheroidal corrections to the spin down can be of the same order of magnitude as those arising from the force-free corrections which reach up to a factor 3 for the orthogonal rotator. This leads to a weaker magnetic field estimate compared to the standard vacuum magneto-dipole losses as shown by (Pétri, 2019).

We infer that the combination of spheroidal geometry and force-free magnetosphere will lower even more the dipole magnetic field strength expectation. So how does the above fitting compare with the force-free fitting of a spherical star as found by Spitkovsky (2006)? A quantitative answer would require computation of force-free spheroidal magnetospheres which is out of the scope of the present work. It is difficult to compare eq. (83) to Spitkovsky (2006) formula because in vacuum we observe only a $L_{\rm vac}\approx L_{\perp}\,\sin^{2}{\mathchoice{\raisebox{0.0pt}{$\displaystyle\chi$}}{\raisebox{0.0pt}{$\textstyle\chi$}}{\raisebox{0.0pt}{$\scriptstyle\chi$}}{\raisebox{0.0pt}{$\scriptscriptstyle\chi$}}}$ dependence, which has to be contrasted with the $L_{\rm FFE}\approx 3/2\,L_{\perp}\,(1+\sin^{2}{\mathchoice{\raisebox{0.0pt}{$\displaystyle\chi$}}{\raisebox{0.0pt}{$\textstyle\chi$}}{\raisebox{0.0pt}{$\scriptstyle\chi$}}{\raisebox{0.0pt}{$\scriptscriptstyle\chi$}}})$ dependence of the force-free model. Nevertheless we guess that both constant $\ell_{1}$ and $\ell_{2}$ in the spheroidal force-free fit $L_{\rm FFE}^{\rm spheroid}\approx 3/2\,L_{\perp}\,(\ell_{1}+\ell_{2}\,\sin^{2}{\mathchoice{\raisebox{0.0pt}{$\displaystyle\chi$}}{\raisebox{0.0pt}{$\textstyle\chi$}}{\raisebox{0.0pt}{$\scriptstyle\chi$}}{\raisebox{0.0pt}{$\scriptscriptstyle\chi$}}})$ will no longer remain constant but depend on the ratio $a/R$ at least to second order $(a/R)^{2}$ such that $\ell_{i}=\alpha_{i}-\beta_{i}\,(a/R)^{2}$ for $i=\{1,2\}$ , $\alpha_{i}$ and $\beta_{i}$ being positive numbers with $\alpha_{i}\approx 1$ due to the spherical force-free results.

5.4 Polar caps

The polar cap rims associated to the field lines structure in Fig. 1 and Fig. 2 as well as for some other obliquities are shown in Fig. 5. We used angles ${\mathchoice{\raisebox{0.0pt}{$\displaystyle\chi$}}{\raisebox{0.0pt}{$\textstyle\chi$}}{\raisebox{0.0pt}{$\scriptstyle\chi$}}{\raisebox{0.0pt}{$\scriptscriptstyle\chi$}}}\in\{0\degr,30\degr,60\degr,90\degr\}$ . As a check, the spherical case is compared to the Deutsch solution shown in orange dashed line and marked with a ”D” in the legend. Both curves overlap to high precision and are indistinguishable.

The oblateness tends to elongate the polar caps in the azimuthal direction, that is in the sense of rotation and a slight contraction in the orthogonal direction for an almost perpendicular rotator. We notice an overall substantial increase in surface area for $a=R$ as seen in the first row. The second row corresponds to the same oblateness but for a spherical dipole boundary condition at the stellar surface. In this case the polar cap rim inflates in all directions whatever the obliquity.

Contrary to an oblate star, the prolate star shrinks its cap shape for almost aligned rotators as seen in the third row. While increasing the inclination angle, the polar cap elongates in the meridional direction with a slight squeezing in the sense of rotation. If spherical dipole boundary conditions are applied at the surface, the variation of polar cap shapes becomes irrelevant for an orthogonal rotator, see fourth row. This is simply explained by the fact that the stellar surface does not significantly vary around the equatorial plane for prolate shapes.

Finding realistic polar cap shapes for fast rotating neutron stars is important to model the hot spot emission seen in thermal X-rays. A careful analysis of such pulsed emission in X-ray from the NICER collaboration led to an estimate of the neutron star compactness and equatorial radius (Riley et al., 2019; Bogdanov et al., 2019). The impact of the stellar surface shape on these hot spots is discussed in Silva et al. (2021), taking also the observer orientation into account. Our simulation could help to reckon even more realistic polar caps.

6 Conclusions

Extending well known results from spherical magnetic stars, we showed how to express multipolar vacuum magnetic fields around spheroidal magnetized objects, being oblate or prolate. Exact analytical solutions have been derived in the case of static stars, involving Legendre functions of the third kind for the radial part, the angular part being expanded into spherical harmonics. For rotating stars, the problem cannot be solved exactly in a closed analytical form because of the introduction of radial and angular spheroidal wave functions. We showed however that some approximate solutions can be found for any realistic configuration, by an expansion into the small parameter $|\gamma|=a/r_{\rm L}\ll 1$ . From a practical point of view the lowest order correction is enough to achieve high accuracy.

As a check, we also solved numerically for spheroidal rotating stars, computing the magnetic field line structure as well as the Poynting flux and the polar cap shape. This study is particular relevant for millisecond pulsars for which strong centrifugal forces inflates the equatorial part leading to an oblate shape. The change in the polar cap rim could have a significant impact on the thermal X-ray emission, modifying their light-curve in addition to general-relativistic effects like frame-dragging and light bending.

The neutron star surrounding is seldom vacuum, a pair plasma usually fills its magnetosphere, producing currents and charges modifying the electromagnetic field outside the star. We plan to add such plasma effects in the description of a spheroidal rotating magnetized celestial body.

Acknowledgements.

I am grateful to the referee for helpful comments and suggestions. This work has been supported by the CEFIPRA grant IFC/F5904-B/2018 and ANR-20-CE31-0010.

References

Abramowitz & Stegun (1965) Abramowitz, M. & Stegun, I. A. 1965, Handbook of mathematical functions with formulas, graphs, and mathematical tables (Dover Books on Advanced Mathematics, New York: Dover, |c1965, Corrected edition, edited by Abramowitz, Milton; Stegun, Irene A.)
Arfken & Weber (2005) Arfken, G. B. & Weber, H.-J. 2005, Mathematical methods for physicists, 6th edn. (Boston: Elsevier)
Asano & Yamamoto (1975) Asano, S. & Yamamoto, G. 1975, Appl. Opt., AO, 14, 29, publisher: Optical Society of America
Bogdanov et al. (2019) Bogdanov, S., Lamb, F. K., Mahmoodifar, S., et al. 2019, ApJL, 887, L26, publisher: American Astronomical Society
Bonazzola et al. (2015) Bonazzola, S., Mottez, F., & Heyvaerts, J. 2015, Astronomy and Astrophysics, 573, A51
Boyd (2001) Boyd, J. P. 2001, Chebyshev and Fourier Spectral Methods (Springer-Verlag)
Chandrasekhar (1970) Chandrasekhar, S. 1970, Ellipsoidal Figures of Equilibrium (New Haven: Yale University Press)
Deutsch (1955) Deutsch, A. J. 1955, Annales d’Astrophysique, 18, 1
Finn & Shapiro (1990) Finn, L. S. & Shapiro, S. L. 1990, The Astrophysical Journal, 359, 444
Gumerov & Duraiswami (2004) Gumerov, N. A. & Duraiswami, R. 2004, Fast multipole methods for the Helmholtz equation in three dimensions, 1st edn., Elsevier series in electromagnetism (Amsterdam: Elsevier), oCLC: 249094845
Horedt (2004) Horedt, G. P. 2004, Polytropes: applications in astrophysics and related fields, Astrophysics and space science library No. v. 306 (Dordrecht ; Boston: Kluwer Academic Publishers)
Jackson (2001) Jackson, J. D. 2001, Electrodynamique classique : Cours et exercices d’electromagnétisme (Paris: Dunod)
Leitner & Spence (1950) Leitner, A. & Spence, R. D. 1950, Journal of the Franklin Institute, 249, 299
Li et al. (2001) Li, L.-W., Kang, X.-K., & Leong, M.-S. 2001, Spheroidal Wave Functions in Electromagnetic Theory, 1st edn. (New York: Wiley-Interscience)
Makishima et al. (2016) Makishima, K., Enoto, T., Murakami, H., et al. 2016, Publications of the Astronomical Society of Japan, 68
Makishima et al. (2019) Makishima, K., Murakami, H., Enoto, T., & Nakazawa, K. 2019, Publications of the Astronomical Society of Japan, 71
Meixner & Schäfke (1954) Meixner, J. & Schäfke, F. W. 1954, Mathieusche Funktionen und Sphäroidfunktionen: Mit Anwendungen auf Physikalische und Technische Probleme, ed. J. Meixner & F. W. Schäfke, Die Grundlehren der Mathematischen Wissenschaften (Berlin, Heidelberg: Springer)
Morse & Feshbach (1953) Morse, P. M. & Feshbach, H. 1953, Methods of theoretical physics I, mcgraw-hill book company edn.
Olver (2010) Olver, F. W. J. 2010, NIST handbook of mathematical functions (Cambridge ; New York: Cambridge University Press : National Institute of Standards and Technology (U.S.)), oCLC: ocn502037224
Pétri (2012) Pétri, J. 2012, MNRAS, 424, 605
Pétri (2013) Pétri, J. 2013, MNRAS, 433, 986
Pétri (2014) Pétri, J. 2014, MNRAS, 439, 1071
Pétri (2015) Pétri, J. 2015, MNRAS, 450, 714
Pétri (2019) Pétri, J. 2019, MNRAS, 485, 4573
Riley et al. (2019) Riley, T. E., Watts, A. L., Bogdanov, S., et al. 2019, ApJL, 887, L21
Shapiro & Teukolsky (1983) Shapiro, S. L. & Teukolsky, S. A. 1983, Black holes, white dwarfs, and neutron stars: The physics of compact objects (Research supported by the National Science Foundation. New York, Wiley-Interscience, 1983, 663 p.)
Silva et al. (2021) Silva, H. O., Pappas, G., Yunes, N., & Yagi, K. 2021, Phys. Rev. D, 103, 063038, publisher: American Physical Society
Spitkovsky (2006) Spitkovsky, A. 2006, ApJ, 648, L51
Zanazzi & Lai (2015) Zanazzi, J. J. & Lai, D. 2015, MNRAS, 451, 695
Zeppenfeld (2009) Zeppenfeld, M. 2009, New J. Phys., 11, 073007

Appendix A Legendre functions of third type

Legendre polynomials $P_{\ell}(x)$ and their associated functions $P_{\ell}^{m}(x)$ are usually defined in the interval $x\in[-1,1]$ . In spheroidal coordinates, we often require solutions of the Legendre differential equations outside this range. For instance when studying fields outside an object up to infinity we require $x>1$ . We are particularly interested in solutions decreasing to zero at large distances, meaning an electric and magnetic potential falling to zero when $x\to+\infty$ . These solutions are then called Legendre functions of the third type $Q_{\ell}^{m}(x)$ and related to prolate coordinates for $|x|\leq 1$ by analytical continuation in the complex plane, starting from $x\in[-1,1]$ . They are solutions of the second order linear differential equation

\frac{d}{dx}\left[(1-x^{2})\,\frac{dQ_{\ell}^{m}(x)}{dx}\right]+\left(\ell\,(\ell+1)-\frac{m^{2}}{1-x^{2}}\right)\,Q_{\ell}^{m}(x)=0

(84)

where $\ell$ and $m$ are two integers. It corresponds to the radial part of a multipole of order $(\ell,m)$ in prolate spheroidal coordinates and related to the spherical harmonic $Y_{\ell}^{m}(\theta,\phi)$ .

For practical purposes, we list the first few useful functions for the monopole $\ell=0$ , the dipole $\ell=1$ and the quadrupole $\ell=2$ , which are real-valued and given by


$\displaystyle{Q}_{0}^{0}(x)$	$\displaystyle=\textrm{arccoth }x$	(85a)
$\displaystyle{Q}_{1}^{0}(x)$	$\displaystyle=x\,\textrm{arccoth }x-1$	(85b)
$\displaystyle{Q}_{1}^{1}(x)$	$\displaystyle=\frac{(x^{2}-1)\,\textrm{arccoth }x-x}{\sqrt{x^{2}-1}}$	(85c)
$\displaystyle{Q}_{2}^{0}(x)$	$\displaystyle=\frac{(3\,x^{2}-1)\,\textrm{arccoth }x-3\,x}{2}$	(85d)
$\displaystyle{Q}_{2}^{1}(x)$	$\displaystyle=\frac{3\,x\,(x^{2}-1)\,\textrm{arccoth }x+2-3\,x^{2}}{\sqrt{x^{2}-1}}$	(85e)
$\displaystyle{Q}_{2}^{2}(x)$	$\displaystyle=3\,(x^{2}-1)\,\textrm{arccoth }x+x\,\frac{5-3\,x^{2}}{x^{2}-1}.$	(85f)

These expressions are used to compute the radial profile of the electric or magnetic potential in prolate spheroidal coordinates.

For oblate spheroidal coordinates, we require solutions $\mathrm{Q}_{\ell}^{m}$ for $x>1$ such that

\frac{d}{dx}\left[(1+x^{2})\,\frac{d\mathrm{Q}_{\ell}^{m}(x)}{dx}\right]+\left(\frac{m^{2}}{1+x^{2}}-\ell\,(\ell+1)\right)\,\mathrm{Q}_{\ell}^{m}(x)=0

(86)

with the correspondence $\mathrm{Q}_{\ell}^{m}(x)=i^{\ell}\,Q_{\ell}^{m}(i\,x)$ . Note the change in sign in front of the factor $x^{2}$ of eq. (86) compared to eq. (84).

For practical purposes, here also we list the first few useful functions for the monopole $\ell=0$ , the dipole $\ell=1$ and the quadrupole $\ell=2$ , which are again all real-valued and given by


$\displaystyle\mathrm{Q}_{0}^{0}(x)$	$\displaystyle=\operatorname{arccot}x$	(87a)
$\displaystyle\mathrm{Q}_{1}^{0}(x)$	$\displaystyle=1-x\,\operatorname{arccot}x$	(87b)
$\displaystyle\mathrm{Q}_{1}^{1}(x)$	$\displaystyle=\frac{x-(1+x^{2})\,\operatorname{arccot}x}{\sqrt{1+x^{2}}}$	(87c)
$\displaystyle\mathrm{Q}_{2}^{0}(x)$	$\displaystyle=\frac{(1+3\,x^{2})\,\operatorname{arccot}x-3\,x}{2}$	(87d)
$\displaystyle\mathrm{Q}_{2}^{1}(x)$	$\displaystyle=\frac{3\,x\,(1+x^{2})\,\operatorname{arccot}x-2-3\,x^{2}}{\sqrt{1+x^{2}}}$	(87e)
$\displaystyle\mathrm{Q}_{2}^{2}(x)$	$\displaystyle=3\,(1+x^{2})\,\operatorname{arccot}x-x\,\left(3+\frac{2}{1+x^{2}}\right).$	(87f)

These expressions are used to compute the radial profile of the electric or magnetic potential in oblate spheroidal coordinates.

Appendix B Angular functions expansion

The prolate angular wave functions $\mathrm{Ps}_{\ell}^{m}(\gamma,\eta)$ are expanded into Legendre functions $P_{\ell}^{m}(\eta)$ according to

\mathrm{Ps}_{\ell}^{m}(\gamma,\eta)={\sum_{r=0,1}^{+\infty}}^{\prime}d_{r}^{m\ell}\,P_{m+r}^{m}(\eta).

(88)


$\displaystyle d_{\ell-m}^{m\ell}$	$\displaystyle\approx 1+O(\gamma^{4})$	(89a)
$\displaystyle d_{\ell-m-2}^{m\ell}$	$\displaystyle\approx-\frac{(\ell+m-1)\,(\ell+m)}{2\,(2\,n-1)^{2}\,(2\,n+1)}\,\gamma^{2}+O(\gamma^{4})$	(89b)
$\displaystyle d_{\ell-m+2}^{m\ell}$	$\displaystyle\approx-\frac{(\ell-m+1)\,(\ell-m+2)}{2\,(2\,n+3)^{2}\,(2\,n+1)}\,\gamma^{2}+O(\gamma^{4}).$	(89c)

If the subscript is negative, the coefficient vanishes by convention. To this level of approximation, we neglect corrections being of order $\gamma^{4}$ or higher and therefore no corrections apply to the dominant coefficient $d_{\ell-m}^{m\ell}\approx 1$ .

For the expansion of the first multipoles with $m=1$ , we give the expression of $d_{r}^{m\ell}$ to $\gamma^{2}$ order for $\ell-m=\pm 2$ for $m=1$ .


	$\displaystyle d_{2}^{1,1}\approx\frac{\gamma^{2}}{75}+O(\gamma^{4})$	(90a)
	$\displaystyle d_{3}^{1,2}\approx\frac{3\,\gamma^{2}}{245}+O(\gamma^{4})$	(90b)
$\displaystyle d_{0}^{1,3}\approx-\frac{6\,\gamma^{2}}{175}+O(\gamma^{4})\qquad;\qquad$	$\displaystyle d_{4}^{1,3}\approx\frac{2\,\gamma^{2}}{189}+O(\gamma^{4})$	(90c)
$\displaystyle d_{1}^{1,4}\approx-\frac{10\,\gamma^{2}}{441}+O(\gamma^{4})\qquad;\qquad$	$\displaystyle d_{5}^{1,4}\approx\frac{10\,\gamma^{2}}{1089}+O(\gamma^{4}).$	(90d)

Explicitly, for the first wave functions, this means


$\displaystyle\mathrm{Ps}_{1}^{1}(\gamma,\eta)$	$\displaystyle\approx P_{1}^{1}(\eta)+\frac{\gamma^{2}}{75}\,P_{3}^{1}(\eta)+O(\gamma^{4})$	(91a)
$\displaystyle\mathrm{Ps}_{2}^{1}(\gamma,\eta)$	$\displaystyle\approx P_{2}^{1}(\eta)+\frac{3\,\gamma^{2}}{245}\,P_{4}^{1}(\eta)+O(\gamma^{4})$	(91b)
$\displaystyle\mathrm{Ps}_{3}^{1}(\gamma,\eta)$	$\displaystyle\approx P_{3}^{1}(\eta)-\frac{6\,\gamma^{2}}{175}\,P_{1}^{1}(\eta)+\frac{2\,\gamma^{2}}{189}\,P_{5}^{1}(\eta)+O(\gamma^{4})$	(91c)
$\displaystyle\mathrm{Ps}_{4}^{1}(\gamma,\eta)$	$\displaystyle\approx P_{4}^{1}(\eta)-\frac{10\,\gamma^{2}}{441}\,P_{2}^{1}(\eta)+\frac{10\,\gamma^{2}}{1089}\,P_{6}^{1}(\eta)+O(\gamma^{4})$	(91d)