© J. Pétri 2022

1. Introduction

Large celestial bodies are approximately spherical because of the preponderance of gravity against other internal forces and stresses. Gravitation does not favor any direction in the sky and therefore a perfect spherical shape is expected for isotropic materials. However, celestial bodies such as stars and molecular clouds are often subject to rotation, producing a centrifugal force F_cen breaking the isotropy imposed by gravity F_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:

$$\begin{array}{c}\hfill \frac{F_{\mathrm{cen}}}{F_{\mathrm{grav}}}\approx \frac{{\mathrm{\Omega }}^2}{{\mathrm{\Omega }}_k^2},\end{array}$$(1)

with Ω being the rotation rate of the star (assumed to be in solid body rotation) and Ω_k being the Keplerian angular frequency at the equator. For a more quantitative description, readers can refer to the discussion from Chandrasekhar (1970) about spheroids and ellipsoids of revolution and to 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 a rotation rate being a fraction of the equatorial Keplerian frequency, we expect their 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

$$\begin{array}{c}\hfill ϵ=\frac{15}{16\pi }\frac{{\mathrm{\Omega }}^2}{G\rho }=\frac{5}{4}\frac{{\mathrm{\Omega }}^2{R}^3}{GM}=\frac{5}{4}\frac{{\mathrm{\Omega }}^2}{{\mathrm{\Omega }}_k^2},\end{array}$$(2)

which is defined by the moment of inertia difference between the 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

$$\begin{array}{c}\hfill e=\sqrt{1-{\left(\frac{R_{\mathrm{pol}}}{R_{\mathrm{eq}}}\right)}^2},\end{array}$$(3)

with R_pol and R_eq being the polar and equatorial radii of the fluid, the frequency Ω given by Chandrasekhar (1970) reads as follows:

$$\begin{array}{c}\hfill \frac{{\mathrm{\Omega }}^2}{2\pi G\rho }=\frac{\sqrt{1-e^2}}{e^3}(3-2e^2)arcsine-3\frac{1-e^2}{e^2}\approx \frac{4}{15}{e}^2.\end{array}$$(4)

The approximation is valid for small ellipticities e ≪ 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 surroundings 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 in 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 aging and spin evolution. So far, exact analytical solutions only exist for a perfect sphere as expressed by Deutsch (1955). Extension to multipolar fields has been thoroughly investigated by Pétri (2015); readers can also refer to 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 11 separable systems (Morse & Feshbach 1953), leading to well-defined solutions for the Laplace and Helmholtz equations. For completeness, we 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 toward more general configurations. Spheroidal coordinate systems possess the interesting and important properties of full separation of a variable for the Laplace and Helmholtz operators. It is therefore possible to solve time harmonic wave emission and propagation in spheroidal coordinates analytically. Moreover for compact objects, anisotropic equation of states for nuclear matter at a very high density as in neutron stars also produce nonspherical astronomical objects.

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

In this paper, we formally compute exact analytical solutions to the electromagnetic wave radiation of a stationary rotating star of a spheroidal shape, oblate, or prolate. In Sect. 2, we recall the useful and important properties of the curvilinear and orthogonal coordinate systems that 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 Sect. 3, we compute static solutions for the multipolar magnetic field sustained by a spheroidal object. We discuss the important question about the normalization of a spheroidal multipole with respect to a spherical multipole being used as the reference solution. Eventually, in Sect. 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 Sect. 5. Conclusions are drawn in Sect. 6.

2. Spheroidal coordinates

We start with a brief overview of both spheroidal coordinate systems, namely oblate and prolate. We 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 a self-gravitating gas well which is 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 (e_x, e_y, e_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 (ρ, ψ, φ), such that the Cartesian coordinates (x, y, z) are given by

$$\begin{array}{cc}& x=\eta sin\psi cos\phi \hfill \end{array}$$(5a)

$$\begin{array}{cc}\hfill & y=\eta sin\psi sin\phi \hfill \end{array}$$(5b)

$$\begin{array}{cc}\hfill & z=\rho cos\psi ,\hfill \end{array}$$(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_{\mathrm{eq}}=\sqrt{R^2+a^2}$. The ellipticity is defined by (Shapiro & Teukolsky 1983)

$$\begin{array}{c}\hfill ϵ=\frac{R_{\mathrm{eq}}-R}{(R_{\mathrm{eq}}+R)/2}.\end{array}$$(6)

The natural basis vectors derived from Eq. (5) are expressed as follows:

$$\begin{array}{cc}& \eta {\mathit{e}}_{\rho }=\rho sin\psi cos\phi {\mathbf{e}}_x+\rho sin\psi sin\phi {\mathbf{e}}_y+\eta cos\psi {\mathbf{e}}_z\hfill \end{array}$$(7a)

$$\begin{array}{cc}\hfill & {\mathit{e}}_{\psi }=\eta cos\psi cos\phi {\mathbf{e}}_x+\eta cos\psi sin\phi {\mathbf{e}}_y-\rho sin\psi {\mathbf{e}}_z\hfill \end{array}$$(7b)

$$\begin{array}{cc}\hfill & {\mathit{e}}_{\phi }=-\eta sin\psi sin\phi {\mathbf{e}}_x+\eta sin\psi cos\phi {\mathbf{e}}_y.\hfill \end{array}$$(7c)

The position vector r is then expanded into

$$\begin{array}{c}\hfill \mathit{r}=\frac{\rho {\eta }^2}{{\mathrm{\Delta }}_o}{\mathit{e}}_{\rho }+\frac{a^{2}cos\psi sin\psi }{{\mathrm{\Delta }}_o}{\mathit{e}}_{\psi },\end{array}$$(8)

with Δ_o = ρ² + a² cos²ψ. The oblate coordinate system is orthogonal, thus leading to a diagonal metric g with the following coefficients:

$$\begin{array}{cc}& g_{\rho \rho }=\frac{{\rho }^2+a^2{cos}^2\psi }{{\eta }^2}=\frac{{\mathrm{\Delta }}_o}{{\eta }^2}\hfill \end{array}$$(9a)

$$\begin{array}{cc}\hfill & g_{\psi \psi }={\rho }^2+a^2{cos}^2\psi ={\mathrm{\Delta }}_o\hfill \end{array}$$(9b)

$$\begin{array}{cc}\hfill & g_{\phi \phi }={\eta }^2{sin}^2\psi .\hfill \end{array}$$(9c)

Its determinant is

$$\begin{array}{c}\hfill \gamma =\det (g)={({\rho }^2+a^2{cos}^2\psi )}^2{sin}^2\psi ={\mathrm{\Delta }}_o^2{sin}^2\psi .\end{array}$$(10)

For the computation of fluxes along the coordinate ρ, it is useful to have the surface element vector dΣ defined for an orthogonal coordinate system such as the oblate one expressed in terms of the variation of the position vector r along two directions e₁ and e₂ such that (with the indices 1 and 2 being ρ, ψ, or φ)

$$\begin{array}{c}\hfill \mathrm{d}\mathbf{\Sigma }=\mathrm{d}{\mathit{r}}_1\wedge \mathrm{d}{\mathit{r}}_2=\frac{\partial \mathit{r}}{\partial x_1}\wedge \frac{\partial \mathit{r}}{\partial x_2}\mathrm{d}{x}_1\mathrm{d}{x}_2.\end{array}$$(11)

Defined by its covariant components, we get

$$\begin{array}{c}\hfill \mathrm{d}{\mathrm{\Sigma }}_i=\sqrt{\gamma }{ϵ}_{i12}\mathrm{d}{x}_1\mathrm{d}{x}_2,\end{array}$$(12)

with ϵ_ijk being the spatial Levi-Civita tensor (Arfken & Weber 2005). For instance, the radial Poynting flux across a closed surface Σ defined by the spheroid ρ = ρ₀ follows as

$$\begin{array}{c}\hfill L={\displaystyle {∯}_{\mathrm{\Sigma }}\mathit{S}·\mathrm{d}\mathbf{\Sigma }={\displaystyle {∯}_{\mathrm{\Sigma }}}{S}^{\rho }d{\mathrm{\Sigma }}_{\rho },}\end{array}$$(13)

with dΣ_ρ = g_ρρ dΣ^ρ = Δ_o sin ψ dψ dφ and S being the Poynting vector. Exactly similar reckoning is performed for prolate coordinates, as shown in the next paragraph.

2.2. Prolate spheroidal coordinates

The prolate spheroidal coordinate system (ρ, ψ, φ) is given by

$$\begin{array}{cc}& x=\rho sin\psi cos\phi \hfill \end{array}$$(14a)

$$\begin{array}{cc}\hfill & y=\rho sin\psi sin\phi \hfill \end{array}$$(14b)

$$\begin{array}{cc}\hfill & z=\eta cos\psi .\hfill \end{array}$$(14c)

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

$$\begin{array}{c}\hfill ϵ=\frac{R_{\mathrm{pol}}-R}{(R_{\mathrm{pol}}+R)/2}.\end{array}$$(15)

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

$$\begin{array}{cc}& \eta {\mathit{e}}_{\rho }=\eta sin\psi cos\phi {\mathbf{e}}_x+\eta sin\psi sin\phi {\mathbf{e}}_y+\rho cos\psi {\mathbf{e}}_z\hfill \end{array}$$(16a)

$$\begin{array}{cc}\hfill & {\mathit{e}}_{\psi }=\rho cos\psi cos\phi {\mathbf{e}}_x+\rho cos\psi sin\phi {\mathbf{e}}_y-\eta sin\psi {\mathbf{e}}_z\hfill \end{array}$$(16b)

$$\begin{array}{cc}\hfill & {\mathbf{e}}_{\phi }=-\rho sin\psi sin\phi {\mathbf{e}}_x+\rho sin\psi cos\phi {\mathbf{e}}_y.\hfill \end{array}$$(16c)

The position vector r is given by

$$\begin{array}{c}\hfill \mathit{r}=\frac{\rho {\eta }^2}{{\mathrm{\Delta }}_p}{\mathit{e}}_{\rho }-\frac{a^{2}cos\psi sin\psi }{{\mathrm{\Delta }}_p}{\mathit{e}}_{\psi },\end{array}$$(17)

with Δ_p = ρ² + a² sin²ψ. The prolate coordinate system is also orthogonal and the metric is given by the following:

$$\begin{array}{cc}& g_{\rho \rho }=\frac{{\rho }^2+a^2{sin}^2\psi }{{\eta }^2}=\frac{{\mathrm{\Delta }}_p}{{\eta }^2}\hfill \end{array}$$(18a)

$$\begin{array}{cc}\hfill & g_{\psi \psi }={\rho }^2+a^2{sin}^2\psi ={\mathrm{\Delta }}_p\hfill \end{array}$$(18b)

$$\begin{array}{cc}\hfill & g_{\phi \phi }={\rho }^2{sin}^2\psi .\hfill \end{array}$$(18c)

Its determinant is

$$\begin{array}{c}\hfill \gamma =det(g)=\frac{{\rho }^2}{{\eta }^2}{({\rho }^2+a^2{sin}^2\psi )}^2{sin}^2\psi =\frac{{\rho }^2{\mathrm{\Delta }}_p^2}{{\eta }^2}{sin}^2\psi .\end{array}$$(19)

The surface element on a surface ρ = ρ₀ is in contravariant components

$$\begin{array}{c}\hfill \mathrm{d}{\mathrm{\Sigma }}^{\rho }=\eta \rho sin\psi \mathrm{d}\psi \mathrm{d}\varphi .\end{array}$$(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 as

$$\begin{array}{c}\hfill \mathrm{\Delta }W+k^{2}W=0,\end{array}$$(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 into three functions of one independent variable each, including a radial part P(ρ), an angular part Ψ(ψ), and an azimuthal part Φ(φ). The second order linear differential equations derived from this separation generate the radial and angular spheroidal wave functions (Olver 2010; Abramowitz & Stegun 1965).

Explicitly writing the solution with the ansatz W(ρ, ψ, φ) = P(ρ) Ψ(ψ) Φ(φ), the separation of the variable leads to an azimuthal dependence Φ(φ)∝e^i m φ, where m is an integer because Φ(φ) must be single-valued in φ ∈ [0, 2 π] and to two second order linear differential equations for ρ ≥ 0 and ψ ∈ [0, π], given by

$$\begin{array}{cc}& \frac{\mathrm{d}}{\mathrm{d}\rho }\left[({\rho }^2+a^2)\frac{\mathrm{d}P}{\mathrm{d}\rho }\right]+[k^2({\rho }^2+a^2)+\frac{m^2{a}^2}{{\rho }^2+a^2}]P=+\lambda P\hfill \end{array}$$(22a)

$$\begin{array}{cc}\hfill & \frac{1}{sin\psi }\frac{\mathrm{d}}{\mathrm{d}\psi }[sin\psi \frac{\mathrm{d}\mathrm{\Psi }}{\mathrm{d}\psi }]-[k^2{a}^2{sin}^2\psi +\frac{m^2}{{sin}^2\psi }]\mathrm{\Psi }=-\lambda \mathrm{\Psi },\hfill \end{array}$$(22b)

respectively, where λ is a separation constant, being an eigenvalue determined by the boundary conditions. By a change to a new independent variable z, letting ρ = ±i a z for Eq. (22a) and z = cos ψ for Eq. (22b), both equations reduce to the same Sturm-Liouville problem

$$\begin{array}{c}\hfill \frac{\mathrm{d}}{\mathrm{d}z}\left[(1-z^2)\frac{\mathrm{d}f}{\mathrm{d}z}\right]+[\lambda +{\gamma }^2(1-z^2)-\frac{m^2}{1-z^2}]f=0,\end{array}$$(23)

with γ² = −k² a² < 0 meaning that γ is purely imaginary. We 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

$$\begin{array}{cc}& P(\xi )=p_1{S_{\ell }^m}^{(1)}(i\xi ,\gamma )+p_2{S_{\ell }^m}^{(2)}(i\xi ,\gamma )\hfill \end{array}$$(24a)

$$\begin{array}{cc}\hfill & \mathrm{\Psi }(\psi )=s_1{\text{Ps}}_{\ell }^m(\psi ,{\gamma }^2)+s_2{\text{Qs}}_{\ell }^m(\psi ,{\gamma }^2)\hfill \end{array}$$(24b)

$$\begin{array}{cc}\hfill & \mathrm{\Phi }(\phi )=h_{1}cos(m\phi )+h_{2}sin(m\phi ),\hfill \end{array}$$(24c)

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

$$\begin{array}{c}\hfill {S_{\ell }^m}^{(3)}(i\xi ,\gamma )={S_{\ell }^m}^{(1)}(i\xi ,\gamma )+i{S_{\ell }^m}^{(2)}(i\xi ,\gamma ),\end{array}$$(25)

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

$$\begin{array}{c}\hfill W={S_{\ell }^m}^{(3)}(i\xi ,\gamma ){\text{Ps}}_{\ell }^m(\psi ,{\gamma }^2)e^{im\phi }.\end{array}$$(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 γ² = k² a² > 0, meaning that γ is real. The solution for outgoing waves now reads as

$$\begin{array}{c}\hfill W={S_{\ell }^m}^{(3)}(\zeta ,\gamma ){\text{Ps}}_{\ell }^m(\psi ,{\gamma }^2)e^{im\phi },\end{array}$$(27)

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

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

3. Static spheroidal star

For nonrotating stars, the Helmholtz equation reduces to the Poisson equation. By introducing a magnetic scalar potential ϕ_M, the solution to the magnetic field structure in vacuum can be solved as 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 B in vacuum outside a static star is curl-free and divergenceless. It can be written as the gradient of a magnetic scalar potential ϕ_M, such that

$$\begin{array}{c}\hfill \mathit{B}=-\mathbf{\nabla }{\varphi }_M.\end{array}$$(28)

The condition ∇ ⋅ B = 0 implies that the scalar potential satisfies the Laplace equation

$$\begin{array}{c}\hfill \mathrm{\Delta }{\varphi }_M=0.\end{array}$$(29)

This equation is fully separable in oblate spheroidal coordinates. Assuming that the inner boundary is located on the surface ρ = ρ₀, the general solution is expanded with ξ = ρ/a into

$$\begin{array}{c}\hfill {\varphi }_M(\rho ,\psi ,\phi )={\displaystyle \sum _{\ell ,|m|\le \ell }(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)})}{Y}_{\ell }^m(\psi ,\phi ),\end{array}$$(30)

where $P_{\mathcal{l}}^m$ and $Q_{\mathcal{l}}^m$ are the Legendre functions of first and second kind, respectively, and $Y_{\mathcal{l}}^m$ are the spherical harmonics (see Appendix A). The potential is imposed at ρ = ρ₀ with ϕ_M(ρ₀, ψ, φ) = V(ψ, φ) and must vanish at infinity at ρ = +∞. Therefore, the coefficients $a_{\mathcal{l}}^m$ vanish. Moreover, the coefficients $b_{\mathcal{l}}^m$ are determined by the decomposition of the surface potential into spherical harmonics

$$\begin{array}{c}\hfill V(\psi ,\phi )={\displaystyle \sum _{\ell ,|m|\le \ell }{V}_{\ell }^m{Y}_{\ell }^m(\psi ,\phi )}\end{array}$$(31)

and then these $b_{\mathcal{l}}^m$ are identified as the surface potential coefficients such that $b_{\mathcal{l}}^m=V_{\mathcal{l}}^m$. The solution for any magnetic potential is therefore

$$\begin{array}{c}\hfill {\varphi }_M(\rho ,\psi ,\phi )={\displaystyle \sum _{\ell ,|m|\le \ell }{V}_{\ell }^m\frac{Q_{\ell }^m(i\xi )}{Q_{\ell }^m(i{\xi }_0)}{Y}_{\ell }^m(\psi ,\phi ).}\end{array}$$(32)

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

$$\begin{array}{c}\hfill B^{\widehat{i}}=-\frac{1}{\sqrt{g_{\mathit{ii}}}}{\partial }_i{\varphi }_M.\end{array}$$(33)

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

$$\begin{array}{c}\hfill \underset{a\to 0}{lim}\frac{Q_{\ell }^m(i\xi )}{Q_{\ell }^m(i{\xi }_0)}={\left(\frac{R}{r}\right)}^{\ell +1}\end{array}$$(34)

and we retrieve the expressions for standard magnetic multipoles (Pétri 2015). We note that in this limit, ρ₀ = R = R_eq.

3.1.1. Monopole solution

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

$$\begin{array}{c}\hfill {\varphi }_M(\rho ,\psi ,\phi )=V_0^0\frac{Q_0^0(i\xi )}{Q_0^0(i{\xi }_0)}{Y}_0^0(\psi ,\phi )=V\frac{\text{arccot}\xi }{\text{arccot}{\xi }_0},\end{array}$$(35)

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

$$\begin{array}{c}\hfill B_{\rho }=-{\partial }_{\rho }{\varphi }_M=V\frac{a}{({\rho }^2+a^2)\text{arccot}{\xi }_0}.\end{array}$$(36)

Introducing the constant magnetic field strength B at the surface ρ = R by the definition B = B_ρ(R), the magnetic surface potential reads as

$$\begin{array}{c}\hfill V=\frac{R^2+a^2}{a}B\text{arccot}{\xi }_0\end{array}$$(37)

and the magnetic field becomes

$$\begin{array}{c}\hfill B_{\rho }=B\frac{R^2+a^2}{{\rho }^2+a^2}.\end{array}$$(38)

The physical component where B is the magnetic field strength at the pole (ρ = R, ψ = 0) is

$$\begin{array}{c}\hfill B_{\widehat{\rho }}=B\frac{R^2+a^2}{\sqrt{({\rho }^2+a^2)({\rho }^2+a^2{cos}^2\psi )}}.\end{array}$$(39)

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

$$\begin{array}{c}\hfill B_{\widehat{\rho }}^{\mathrm{eq}}=B\sqrt{1+\frac{a^2}{R^2}}.\end{array}$$(40)

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

$$\begin{array}{c}\hfill B_r=B\frac{R^2}{r^2}.\end{array}$$(41)

In the case of a small deformation with a ≪ R, the field expands into

$$\begin{array}{c}\hfill B_{\widehat{\rho }}\approx B\frac{R^2}{{\rho }^2}[1+\frac{a^2}{R^2}(1-\frac{R^2}{2{\rho }^2}(1+{cos}^2\psi ))].\end{array}$$(42)

We stress that the ρ component does not correspond to the spherical radial component and as such the monopolar magnetic field actually has two components that are both contained in the meridional plane. In the spherical coordinate system (r, θ, φ), B_ρ decomposes into a B_r and a B_θ component because of the natural basis expressions in Eq. (7).

At large distances, the field simplifies into

$$\begin{array}{c}\hfill B_{\widehat{\rho }}\approx B\frac{R^2}{{\rho }^2}[1+\frac{a^2}{R^2}]\approx B\frac{R^2}{r^2}[1+\frac{a^2}{R^2}],\end{array}$$(43)

showing that the oblateness still has an imprint far away from the surface depicted by the correcting factor (1 + a²/R²). 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 for an aligned and an orthogonal dipole separately. 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 a parallel $V_1^0$ and perpendicular $V_1^1$ contribution is

$$\begin{array}{c}\hfill {\varphi }_M(\rho ,\psi ,\phi )=V_1^0\frac{Q_1^0(i\xi )}{Q_1^0(i{\xi }_0)}{Y}_1^0(\psi ,\phi )+V_1^1\frac{Q_1^1(i\xi )}{Q_1^1(i{\xi }_0)}{Y}_1^1(\psi ,\phi ),\end{array}$$(44)

or again in getting rid of spherical harmonic normalization factors with parallel V_∥ and perpendicular V_⊥ contributions and a possible negative sign for the magnetic field components, we write

$$\begin{array}{cc}& -{\varphi }_M(\rho ,\psi ,\phi )=V_{\Vert }\frac{\xi \text{arccot}\xi -1}{{\xi }_0\text{arccot}{\xi }_0-1}cos\psi \hfill \\ \hfill & +V_{\perp }\frac{({\xi }^2+1)\text{arccot}\xi -\xi }{({\xi }_0^2+1)\text{arccot}{\xi }_0-{\xi }_0}\sqrt{\frac{{\xi }_0^2+1}{{\xi }^2+1}}sin\psi e^{i\phi }.\hfill \end{array}$$(45)

The magnetic field is decomposed into an aligned rotator as

$$\begin{array}{cc}& B_{\rho }=\frac{V_{\Vert }}{a}\frac{({\xi }^2+1)\text{arccot}\xi -\xi }{({\xi }^2+1)({\xi }_0\text{arccot}{\xi }_0-1)}cos\psi \hfill \end{array}$$(46a)

$$\begin{array}{cc}\hfill & B_{\psi }=-V_{\Vert }\frac{\xi \text{arccot}\xi -1}{{\xi }_0\text{arccot}{\xi }_0-1}sin\psi \hfill \end{array}$$(46b)

$$\begin{array}{cc}\hfill & B_{\phi }=0,\hfill \end{array}$$(46c)

or by introducing a typical surface magnetic field strength B_∥ = B_ρ(R, ψ = 0)

$$\begin{array}{cc}& B_{\rho }=2B_{\Vert }\frac{{\xi }_0^2+1}{{\xi }^2+1}\frac{({\xi }^2+1)\text{arccot}\xi -\xi }{({\xi }_0^2+1)\text{arccot}{\xi }_0-{\xi }_0}cos\psi \hfill \end{array}$$(47a)

$$\begin{array}{cc}\hfill & B_{\psi }=2a{B}_{\Vert }\frac{({\xi }_0^2+1)(1-\xi \text{arccot}\xi )}{({\xi }_0^2+1)\text{arccot}{\xi }_0-{\xi }_0}sin\psi \hfill \end{array}$$(47b)

$$\begin{array}{cc}\hfill & B_{\phi }=0\hfill \end{array}$$(47c)

and an orthogonal rotator as

$$\begin{array}{cc}& B_{\rho }=\frac{V_{\perp }}{a}\frac{({\xi }^2+1)\xi \text{arccot}\xi -{\xi }^2-2}{({\xi }_0^2+1)\text{arccot}{\xi }_0-{\xi }_0}\frac{\sqrt{{\xi }_0^2+1}}{{({\xi }^2+1)}^{3/2}}sin\psi e^{i\phi }\hfill \end{array}$$(48a)

$$\begin{array}{cc}\hfill & B_{\psi }=V_{\perp }\frac{({\xi }^2+1)\text{arccot}\xi -\xi }{({\xi }_0^2+1)\text{arccot}{\xi }_0-{\xi }_0}\sqrt{\frac{{\xi }_0^2+1}{{\xi }^2+1}}cos\psi e^{i\phi }\hfill \end{array}$$(48b)

$$\begin{array}{cc}\hfill & B_{\phi }=i{V}_{\perp }\frac{({\xi }^2+1)\text{arccot}\xi -\xi }{({\xi }_0^2+1)\text{arccot}{\xi }_0-{\xi }_0}\sqrt{\frac{{\xi }_0^2+1}{{\xi }^2+1}}sin\psi e^{i\phi },\hfill \end{array}$$(48c)

or by introducing a typical surface magnetic field strength B_⊥ = B_ρ(R, ψ = π/2, φ = 0)

$$\begin{array}{cc}& B_{\rho }=2B_{\perp }\frac{{\xi }^2+2-({\xi }^2+1)\xi \text{arccot}\xi }{{\xi }_0^2+2-({\xi }_0^2+1){\xi }_0\text{arccot}{\xi }_0}{\left(\frac{{\xi }_0^2+1}{{\xi }^2+1}\right)}^{3/2}sin\psi e^{i\phi }\hfill \end{array}$$(49a)

$$\begin{array}{cc}\hfill & B_{\psi }=2a{B}_{\perp }\frac{({\xi }^2+1)\text{arccot}\xi -\xi }{({\xi }_0^2+1){\xi }_0\text{arccot}{\xi }_0-{\xi }_0^2-2}\frac{{({\xi }_0^2+1)}^{3/2}}{\sqrt{{\xi }^2+1}}cos\psi e^{i\phi }\hfill \end{array}$$(49b)

$$\begin{array}{cc}\hfill & B_{\phi }=2ia{B}_{\perp }\frac{({\xi }^2+1)\text{arccot}\xi -\xi }{({\xi }_0^2+1){\xi }_0\text{arccot}{\xi }_0-{\xi }_0^2-2}\frac{{({\xi }_0^2+1)}^{3/2}}{\sqrt{{\xi }^2+1}}sin\psi e^{i\phi }.\hfill \end{array}$$(49c)

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

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

$$\begin{array}{cc}& B_{\rho }=\frac{4}{3}{B}_{\Vert }\frac{a^3}{{\rho }^3}\frac{{\xi }_0^2+1}{({\xi }_0^2+1)\text{arccot}{\xi }_0-{\xi }_0}cos\psi \hfill \end{array}$$(50a)

$$\begin{array}{cc}\hfill & B_{\psi }=\frac{2}{3}{B}_{\Vert }\frac{a^3}{{\rho }^2}\frac{{\xi }_0^2+1}{({\xi }_0^2+1)\text{arccot}{\xi }_0-{\xi }_0}sin\psi \hfill \end{array}$$(50b)

$$\begin{array}{cc}\hfill & B_{\phi }=0\hfill \end{array}$$(50c)

for the aligned part and

$$\begin{array}{cc}& B_{\rho }=\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\text{arccot}{\xi }_0}sin\psi e^{i\phi }\hfill \end{array}$$(51a)

$$\begin{array}{cc}\hfill & B_{\psi }=\frac{4}{3}{B}_{\perp }\frac{a^3}{{\rho }^2}\frac{{({\xi }_0^2+1)}^{3/2}}{({\xi }_0^2+1){\xi }_0\text{arccot}{\xi }_0-{\xi }_0^2-2}cos\psi e^{i\phi }\hfill \end{array}$$(51b)

$$\begin{array}{cc}\hfill & B_{\phi }=\frac{4}{3}i{B}_{\perp }\frac{a^3}{{\rho }^2}\frac{{({\xi }_0^2+1)}^{3/2}}{({\xi }_0^2+1){\xi }_0\text{arccot}{\xi }_0-{\xi }_0^2-2}sin\psi e^{i\phi }\hfill \end{array}$$(51c)

for the orthogonal part, respectively.

3.2. Prolate magnetic star

Following the same procedure as for the oblate magnetic star, we find the magnetic potential ϕ_M for prolate star as

$$\begin{array}{c}\hfill {\varphi }_M(\rho ,\psi ,\phi )={\displaystyle \sum _{\ell ,|m|\le \ell }{V}_{\ell }^m\frac{Q_{\ell }^m(\zeta )}{Q_{\ell }^m({\zeta }_0)}{Y}_{\ell }^m(\psi ,\phi )}\end{array}$$(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 be around zero and the potential simplifies to a multipole such as

$$\begin{array}{c}\hfill \underset{a\to 0}{lim}\frac{Q_{\ell }^m(\zeta )}{Q_{\ell }^m({\zeta }_0)}={\left(\frac{R}{r}\right)}^{\ell +1},\end{array}$$(53)

which is the same expression as in the previous subsection because oblate and prolate both tend to have a spherical shape when a → 0.

3.2.1. Monopole solution

For the monopole, the potential reads explicitly as

$$\begin{array}{c}\hfill {\varphi }_M(\rho ,\psi ,\phi )=V_0^0\frac{Q_0^0(\zeta )}{Q_0^0({\zeta }_0)}{Y}_0^0(\psi ,\phi )=V\frac{\text{arccoth}\zeta }{\text{arccoth}{\zeta }_0}.\end{array}$$(54)

It actually only depends on the ρ coordinate. The magnetic field therefore only has a ρ component

$$\begin{array}{c}\hfill B_{\rho }=-{\partial }_{\rho }{\varphi }_M=V\frac{a}{\rho \sqrt{{\rho }^2+a^2}\text{arccoth}{\zeta }_0}.\end{array}$$(55)

At the stellar surface, ρ = R and B_ρ = B, thus

$$\begin{array}{c}\hfill V=\frac{R\sqrt{R^2+a^2}}{a}B\text{arccoth}{\zeta }_0,\end{array}$$(56)

such that

$$\begin{array}{c}\hfill B_{\rho }=B\frac{R}{\rho }\sqrt{\frac{R^2+a^2}{{\rho }^2+a^2}}.\end{array}$$(57)

The physical component where B is the magnetic field strength at the equator (ρ = R, ψ = π/2) is

$$\begin{array}{c}\hfill B_{\widehat{\rho }}=B\frac{R}{\rho }\sqrt{\frac{R^2+a^2}{{\rho }^2+a^2{sin}^2\psi }}.\end{array}$$(58)

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

$$\begin{array}{c}\hfill B_{\widehat{\rho }}^{\mathrm{pole}}=B\sqrt{1+\frac{a^2}{R^2}}.\end{array}$$(59)

For small prolateness a ≪ R, the potential and the field become

$$\begin{array}{cc}& {\varphi }_M\approx V\frac{R}{\rho }[1+\frac{1}{6}\frac{a^2}{R^2}(1-\frac{R^2}{{\rho }^2})]\hfill \end{array}$$(60a)

$$\begin{array}{cc}\hfill & B_{\widehat{\rho }}\approx B\frac{R^2}{{\rho }^2}[1+\frac{1}{2}\frac{a^2}{R^2}(1-\frac{R^2}{{\rho }^2}{sin}^2\psi )],\hfill \end{array}$$(60b)

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

At large distances, the physical component reduces to

$$\begin{array}{c}\hfill B_{\widehat{\rho }}=B\frac{R^2}{{\rho }^2}\sqrt{1+\frac{a^2}{R^2}}.\end{array}$$(61)

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

3.2.2. Dipole solution

For the prolate dipole, the potential with a parallel $V_1^0$ and perpendicular $V_1^1$ contribution is

$$\begin{array}{c}\hfill {\varphi }_M(\rho ,\psi ,\phi )=V_1^0\frac{Q_1^0(\zeta )}{Q_1^0({\zeta }_0)}{Y}_1^0(\psi ,\phi )+V_1^1\frac{Q_1^1(\zeta )}{Q_1^1({\zeta }_0)}{Y}_1^1(\psi ,\phi ),\end{array}$$(62)

or explicitly with parallel V_∥ and perpendicular V_⊥ contributions (getting rid of spherical harmonic normalization factors)

$$\begin{array}{cc}& -{\varphi }_M(\rho ,\psi ,\phi )=V_{\Vert }\frac{\zeta \text{arccoth}\zeta -1}{{\zeta }_0\text{arccoth}{\zeta }_0-1}cos\psi \hfill \\ \hfill & +V_{\perp }\frac{({\zeta }^2-1)\text{arccoth}\zeta -\zeta }{({\zeta }_0^2-1)\text{arccoth}{\zeta }_0-{\zeta }_0}\sqrt{\frac{{\zeta }_0^2-1}{{\zeta }^2-1}}sin\psi e^{i\phi }.\hfill \end{array}$$(63)

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

$$\begin{array}{cc}& B_{\rho }=-\frac{V_{\Vert }}{\rho }[1-\frac{{\rho }^2}{a^2}\frac{\text{arccoth}\zeta }{\zeta }]\frac{1}{{\zeta }_0\text{arccoth}{\zeta }_0-1}cos\psi \hfill \end{array}$$(64a)

$$\begin{array}{cc}\hfill & B_{\psi }=V_{\Vert }\frac{1-\zeta \text{arccoth}\zeta }{{\zeta }_0\text{arccoth}{\zeta }_0-1}sin\psi \hfill \end{array}$$(64b)

$$\begin{array}{cc}\hfill & B_{\phi }=0.\hfill \end{array}$$(64c)

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

$$\begin{array}{cc}& B_{\rho }=2B_{\Vert }\frac{R}{\rho }[1-{\xi }^2\frac{\text{arccoth}\zeta }{\zeta }]{[1-{\xi }_0^2\frac{\text{arccoth}{\zeta }_0}{{\zeta }_0}]}^{-1}cos\psi \hfill \end{array}$$(65a)

$$\begin{array}{cc}\hfill & B_{\psi }=2B_{\Vert }R[\zeta \text{arccoth}\zeta -1]{[1-{\xi }_0^2\frac{\text{arccoth}{\zeta }_0}{{\zeta }_0}]}^{-1}sin\psi \hfill \end{array}$$(65b)

$$\begin{array}{cc}\hfill & B_{\phi }=0.\hfill \end{array}$$(65c)

For the orthogonal rotator, we find

$$\begin{array}{cc}& B_{\rho }=\frac{V_{\perp }}{a}\frac{\sqrt{{\zeta }_0^2-1}(({\zeta }^2-1)\zeta \text{arccoth}\zeta -{\zeta }^2+2)}{(({\zeta }_0^2-1)\text{arccoth}{\zeta }_0-{\zeta }_0)({\zeta }^2-1)\zeta }sin\psi e^{i\phi }\hfill \end{array}$$(66a)

$$\begin{array}{cc}\hfill & B_{\psi }=V_{\perp }\frac{({\zeta }^2-1)\text{arccoth}\zeta -\zeta }{({\zeta }_0^2-1)\text{arccoth}{\zeta }_0-{\zeta }_0}\sqrt{\frac{{\zeta }_0^2-1}{{\zeta }^2-1}}cos\psi e^{i\phi }\hfill \end{array}$$(66b)

$$\begin{array}{cc}\hfill & B_{\phi }=i{V}_{\perp }\frac{({\zeta }^2-1)\text{arccoth}\zeta -\zeta }{({\zeta }_0^2-1)\text{arccoth}{\zeta }_0-{\zeta }_0}\sqrt{\frac{{\zeta }_0^2-1}{{\zeta }^2-1}}sin\psi e^{i\phi }\hfill \end{array}$$(66c)

or with the typical surface magnetic field strength B_⊥ we find

$$\begin{array}{cc}& B_{\rho }=2B_{\perp }\frac{{\zeta }_0}{\zeta }\frac{{\zeta }_0^2-1}{{\zeta }^2-1}\frac{({\zeta }^2-1)\zeta \text{arccoth}\zeta -{\zeta }^2+2}{({\zeta }_0^2-1){\zeta }_0\text{arccoth}{\zeta }_0-{\zeta }_0^2+2}sin\psi e^{i\phi }\hfill \end{array}$$(67a)

$$\begin{array}{cc}\hfill & B_{\psi }=2a{B}_{\perp }\frac{{\zeta }_0({\zeta }_0^2-1)}{\sqrt{{\zeta }^2-1}}\frac{({\zeta }^2-1)\text{arccoth}\zeta -\zeta }{({\zeta }_0^2-1){\zeta }_0\text{arccoth}{\zeta }_0-{\zeta }_0^2+2}cos\psi e^{i\phi }\hfill \end{array}$$(67b)

$$\begin{array}{cc}\hfill & B_{\phi }=2ia{B}_{\perp }\frac{{\zeta }_0({\zeta }_0^2-1)}{\sqrt{{\zeta }^2-1}}\frac{({\zeta }^2-1)\text{arccoth}\zeta -\zeta }{({\zeta }_0^2-1){\zeta }_0\text{arccoth}{\zeta }_0-{\zeta }_0^2+2}sin\psi e^{i\phi .}\hfill \end{array}$$(67c)

At large distances, for the aligned component we get

$$\begin{array}{cc}& B_{\rho }=\frac{4}{3}{B}_{\Vert }\frac{a^3}{{\rho }^3}\frac{{\zeta }_0\sqrt{{\zeta }_0^2-1}}{{\zeta }_0+(1-{\zeta }_0^2)\text{arccoth}{\zeta }_0}cos\psi \hfill \end{array}$$(68a)

$$\begin{array}{cc}\hfill & B_{\psi }=\frac{2}{3}{B}_{\Vert }\frac{a^3}{{\rho }^2}\frac{{\zeta }_0\sqrt{{\zeta }_0^2-1}}{{\zeta }_0+(1-{\zeta }_0^2)\text{arccoth}{\zeta }_0}sin\psi \hfill \end{array}$$(68b)

$$\begin{array}{cc}\hfill & B_{\phi }=0,\hfill \end{array}$$(68c)

and for the orthogonal component we get

$$\begin{array}{cc}& B_{\rho }=\frac{8}{3}{B}_{\perp }\frac{a^3}{{\rho }^3}\frac{{\zeta }_0({\zeta }_0^2-1)sin\psi e^{i\phi }}{({\zeta }_0^2-1){\zeta }_0\text{arccoth}{\zeta }_0-{\zeta }_0^2+2}\hfill \end{array}$$(69a)

$$\begin{array}{cc}\hfill & B_{\psi }=\frac{4}{3}{B}_{\perp }\frac{a^3}{{\rho }^2}\frac{{\zeta }_0({\zeta }_0^2-1)cos\psi e^{i\phi }}{{\zeta }_0^2-2-({\zeta }_0^2-1){\zeta }_0\text{arccoth}{\zeta }_0}\hfill \end{array}$$(69b)

$$\begin{array}{cc}\hfill & B_{\phi }=\frac{4}{3}i{B}_{\perp }\frac{a^3}{{\rho }^2}\frac{sin\psi e^{i\phi }}{{\zeta }_0^2-2-({\zeta }_0^2-1){\zeta }_0\text{arccoth}{\zeta }_0}.\hfill \end{array}$$(69c)

This achieves the implementation of the initial background magnetic field setup to start the numerical simulations in Sect. 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 such as 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 obviously exist several approaches to fix the electromagnetic field at the surface such as keeping a fixed value of the magnetic field strength at the pole or at the equator when deforming the stellar surface. Nevertheless, this does not seem satisfactory because the spin-down luminosity, for instance, not only varies because of the spheroidal shape but also because of the artificial field strength variation related to the evolving boundary. This problem is reminiscent of 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, always tending to be the same perfect spherical dipole keeping a constant asymptotic expression. If we were to 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 it with previous results such as from 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 behavior at large distances, it switches to a wave nature around the light-cylinder due to the finite speed of propagation of electromagnetic fields F. In this section, we solve in the whole space outside the star for these fields. Both the electric and the magnetic part satisfy the vector wave equation in vacuum given by

$$\begin{array}{c}\hfill \frac{1}{c^2}\frac{{\partial }^2\mathbf{F}}{\partial t^2}-\mathrm{\Delta }\mathbf{F}=\mathbf{0}.\end{array}$$(70)

For periodic motion, the time dependence becomes harmonic and the field F varies in time according to F ∝ e^{−i ω t}. Introducing the wave number k = ω/c, the vector wave equation reduces to the vector Helmholtz equation

$$\begin{array}{c}\hfill \mathrm{\Delta }\mathbf{F}+k^2\mathbf{F}=\mathbf{0}.\end{array}$$(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, we summarise important results about the scalar Helmholtz equation (21). It can be shown that if W is a solution of (21), then ∇W, Φ = ∇ ∧ (W r), and Ψ = ∇ ∧ Φ are all solutions of the vector Helmholtz equation (71) (Leitner & Spence 1950; Gumerov & Duraiswami 2004). Moreover, Φ and Ψ are solenoidal, meaning that ∇ ⋅ Φ = ∇ ⋅ Ψ = 0 and they satisfy ∇ ∧ Ψ = k² Φ.

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

$$\begin{array}{cc}& \mathit{E}={\displaystyle \sum _{\ell ,m}}{a}_{\ell m}^{\mathrm{E}}{\mathbf{\Psi }}_{\ell m}+b_{\ell m}^{\mathrm{E}}{\mathbf{\Phi }}_{\ell m}\hfill \end{array}$$(72a)

$$\begin{array}{cc}\hfill & c\mathit{B}={\displaystyle \sum _{\ell ,m}}{a}_{\ell m}^{\mathrm{B}}{\mathbf{\Psi }}_{\ell m}+b_{\ell m}^{\mathrm{B}}{\mathbf{\Phi }}_{\ell m,}\hfill \end{array}$$(72b)

automatically satisfying ∇ ⋅ E = ∇ ⋅ B = 0. The numbers ℓ and m label the multipole expansion modes similarly to the vector spherical harmonics (Pétri 2013). The coefficients $a_{\ell m}^{\mathrm{E}/\mathrm{B}}$ and $b_{\ell m}^{\mathrm{E}/\mathrm{B}}$ are constants of integration depending on the boundary conditions. In order to satisfy the time-harmonic Maxwell equations with temporal behavior e^{−i ω t}, these coefficients must be related by