Open Access
Volume 657, January 2022
Article Number A73
Number of page(s) 13
Section Astrophysical processes
Published online 14 January 2022

© J. Pétri 2022

Licence Creative CommonsOpen Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

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 Fcen breaking the isotropy imposed by gravity Fgrav. 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:


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


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


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


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 (ex, ey, ez) 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




where we define and a is a real and positive parameter related to the oblateness. The equatorial radius of the surface is then . The ellipticity is defined by (Shapiro & Teukolsky 1983)


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




The position vector r is then expanded into


with Δo = ρ2 + a2 cos2ψ. The oblate coordinate system is orthogonal, thus leading to a diagonal metric g with the following coefficients:




Its determinant is


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 e1 and e2 such that (with the indices 1 and 2 being ρ, ψ, or φ)


Defined by its covariant components, we get


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 ρ = ρ0 follows as


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




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 . The definition of the ellipticity must be changed accordingly, such that


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




The position vector r is given by


with Δp = ρ2 + a2 sin2ψ. The prolate coordinate system is also orthogonal and the metric is given by the following:




Its determinant is


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


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


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 Φ(φ)∝eimφ, 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



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


with γ2 = −k2a2 < 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




with ξ = ρ/a. The radial (j = 1, 2) and angular 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 with


which represents a generalization of the spherical Hankel functions and . The angular regularity condition imposes s2 = 0. In complex form, a particular solution for the scalar Helmholtz equation Eq. (21) with outgoing wave boundary conditions is


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 γ2 = k2a2 > 0, meaning that γ is real. The solution for outgoing waves now reads as


with . The arguments of 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


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


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


where and are the Legendre functions of first and second kind, respectively, and are the spherical harmonics (see Appendix A). The potential is imposed at ρ = ρ0 with ϕM(ρ0, ψ, φ) = V(ψ, φ) and must vanish at infinity at ρ = +∞. Therefore, the coefficients vanish. Moreover, the coefficients are determined by the decomposition of the surface potential into spherical harmonics


and then these are identified as the surface potential coefficients such that . The solution for any magnetic potential is therefore


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


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


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

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


getting rid of the spherical harmonic normalization factor by setting . This potential actually only depends on the coordinate ρ. The magnetic field therefore only has a ρ component,


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


and the magnetic field becomes


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


The physical component at the equator is stronger because


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


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


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 Br and a Bθ component because of the natural basis expressions in Eq. (7).

At large distances, the field simplifies into


showing that the oblateness still has an imprint far away from the surface depicted by the correcting factor (1 + a2/R2). 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 and perpendicular contribution is


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


The magnetic field is decomposed into an aligned rotator as




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




and an orthogonal rotator as




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




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 for aligned and for perpendicular rotators.

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




for the aligned part and




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


with and .

For small prolateness tending to a spherical shape, a tends to be around zero and the potential simplifies to a multipole such as


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


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


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


such that


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


The physical component at the pole is stronger because


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



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

At large distances, the physical component reduces to


Here we also observe a correcting factor, but with a value of compared to a spherical monopole.

3.2.2. Dipole solution

For the prolate dipole, the potential with a parallel and perpendicular contribution is


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


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




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




For the orthogonal rotator, we find




or with the typical surface magnetic field strength B we find




At large distances, for the aligned component we get




and for the orthogonal component we get




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


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


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, Φ =  ∧ (Wr), 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  ∧ Ψ = k2Φ.

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



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 and are constants of integration depending on the boundary conditions. In order to satisfy the time-harmonic Maxwell equations with temporal behavior eiωt, these coefficients must be related by



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



The only independent coefficients are therefore . The central problem is to find explicit expressions for the vector Ψm and Φm in spheroidal coordinates and to adjust the coefficients to fit the stellar boundary conditions.

4.2. Electric field

As is 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 E + v ∧ B = 0, where


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ρ is unconstrained, Eφ = 0, and


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 γ2 = ±k2a2 always remains less than one in modulus because the star must remain within its light-cylinder thus the constrain ka = a/rL ≪ 1. Therefore we can expand the solution into a series in γ that converges quickly to the exact expression. In this section we follow this path to get more insight as to the impact of oblateness and prolateness on the Poynting flux. Our starting points are the expansion formulae for the angular and radial spheroidal wave functions as given by Abramowitz & Stegun (1965). We summarize the important results that are 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 ρ → ±iρ and γ → ∓iγ. For brevity, we only give results for prolate shapes.

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


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

The prolate radial wave functions are then expanded to


where z is any of the spherical Bessel functions j, y or spherical Hankel functions, or of order . For our problem of outgoing wave solutions, we require radial expansion into a spherical Hankel function of type as for the Deutsch solution.

For analytically tractable purpose, we expand all parameters to second order in γ. Actually, the coefficients only depend on even powers of γ; therefore, the expansion of any quantity also follows an expansion in even powers of γ. Consequently, the first correcting term for a magnetic field structure, Poynting flux, electromagnetic torque, and so on depends on γ2. The key expansion coefficients are provided 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 (,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 which are of much higher order in spin parameter w = Ω R/c, at least w4 compared to 1 and w2 for the dipole. This is due to the fact that the quadrupole is already the result of rotating a magnetic field, thus an w2 strength for a w2 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 only. The Poynting flux in such an approximation which behaves as follows:


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


For outgoing wave boundary conditions, we set . Following the procedure explained in detail in Pétri (2015), the spin-down dependence on a/R and R/rL is approximately proportional to . An expansion to lowest order in these two parameters gives a correction to the luminosity as


where the spin-down L of the vacuum orthogonal point magnetic dipole is


The above approximation gives only some important hints about the spin-down change due to the spheroidal shape. It does not take the normalization of the surface magnetic field strength into account. 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 one 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 does not decompose easily into a 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 nontrivial task. Nevertheless, it clearly shows that leading corrections scale as and a2/R2, thus they are of second order in a. We use these 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 as to 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 the boundary conditions on the surface into account with a 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ρ × Nψ × Nφ = 129 × 32 × 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/rL = 0.3, which coincides with the inner boundary of the computational domain R1 = R. The outer boundary is equal to R2/rL = 7. The oblateness or prolateness is controlled by the parameter a defining η 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 χ is taken in the set χ = {0°, 30°, 60°, 90°}.

5.2. Magnetic field lines

We 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 in Figs. 1 and 2 for oblate and prolate stars with different boundary conditions, either being a single multipolar spheroidal field in vacuum or a spherical field in vacuum, respectively. The spheroidal parameter is chosen as a/R = {0, 0.5, 1}.

thumbnail Fig. 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} in blue, green, and red, respectively. The blue quarter disk in the bottom left depicts the spherical star. The vertical dashed black line represents the light-cylinder.

thumbnail Fig. 2.

Magnetic field lines in the equatorial plane xOy for a perpendicular spheroidal rotator with oblateness (first two columns) or prolateness (last two columns) parameter a/R = {0, 0.5, 1} in blue, green, and red, respectively. The blue disk in the centre depicts the spherical star. The dashed black circle represents the light-cylinder.

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 as to the nature of the spheroidal star. The impact is highest on the surface, and can be quantified by the polar cap rim change as is shown in the corresponding subsection.

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 its size does not vary at the equator. The two-armed spiral pattern typical of the 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 as 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 χ such as sin2χ as is 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 Sect. 3. The obliquity is set to χ = {0°, 30°, 60°, 90°}, as given in the legend. Figure 4 shows the same evolution, but for a star keeping a perfect spherical dipole structure.

thumbnail Fig. 3.

Spin-down luminosity for oblate and prolate stars in solid and dashed lines, respectively, with single dipole stellar boundary conditions.

thumbnail Fig. 4.

Spin-down luminosity for oblate and prolate stars in solid and dashed lines, respectively, with spherical dipole stellar boundary conditions.

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


The coefficients are listed in Table 1.

Table 1.

Fitted coefficients k1 and k2 as given by Eq. (83).

In fact, the coefficient k1 is known analytically and expressed in terms of spherical Hankel functions (see Pétri 2015). For the particular values of our simulations, we should find approximately k1 ≈ 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 surfaces. As carefully shown in Pétri (2014), the accuracy scales as . However, the error remains very weak, amounting to only 0.2%.

This spin down estimate has an impact on the stellar magnetic field inference of accreting pulsars. For instance, the 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 and a radius R = 12 km, it would have a/R ≈ 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 ≈ 0.48, but it is 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 is done in Sect. 3. But we could keep the magnetic moment or the magnetic flux threading the star constant or the asymptotic field structure at a large distance as is done in the numerical simulations. Therefore the central question becomes how to compare a spherical star to a spheroidal one. There is no unique answer and the best normalization must be adapted for 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 the dipole magnetic field strength expectation even more. The above fitting can be compared 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 LvacL sin2χ dependence, which has to be contrasted with the LFFE ≈ 3/2 L (1+sin2χ) dependence of the force-free model. Nevertheless we guess that both constant 1 and 2 in the spheroidal force-free fit will no longer remain constant, but depend on the ratio a/R at least to second order (a/R)2 such that i = αi − βi (a/R)2 for i = {1, 2}, αi and βi being positive numbers with αi ≈ 1 due to the spherical force-free results.

5.4. Polar caps

The polar cap rims associated with the field lines structure shown in Fig. 1 and Fig. 2 for χ = 0° and χ = 90°, but also with some other obliquities are shown in Fig. 5. We used angles χ = {0°, 30°, 60°, 90°}. As a check, the spherical case is compared to the Deutsch solution shown in an orange dashed line and marked with a “D” in the legend. Both curves overlap to high precision and are indistinguishable.

thumbnail Fig. 5.

Polar cap shape for oblate and prolate stars with oblateness parameter a/R = {0, 0.5, 1} in blue, green, and red, respectively. The black dashed line shows the reference solution for the Deutsch field as a check. The obliquity from the left column to the right column is χ = {0°, 30°, 60°, 90°}. The first row is for an oblate star with one mode  = 1, the second row is for a spherical dipole magnetic field at the surface, the third row is for a prolate star with one mode  = 1, and the fourth row is for a spherical dipole magnetic field at the surface.

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), also taking 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 have shown 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 have shown however that some approximate solutions can be found for any realistic configuration by an expansion into the small parameter |γ|=a/rL ≪ 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 particularly relevant for millisecond pulsars for which strong centrifugal forces inflate 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 such as frame-dragging and light bending.

The neutron star surrounding is seldom vacuum, and 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 to the description of a spheroidal rotating magnetized celestial body.


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.


  1. Abramowitz, M., & Stegun, I. A. 1965, Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables (Dover Books on Advanced Mathematics, New York: Dover) [Google Scholar]
  2. Arfken, G. B., & Weber, H.-J. 2005, Mathematical Methods for Physicists, 6th edn. (Boston: Elsevier) [Google Scholar]
  3. Asano, S., & Yamamoto, G. 1975, Appl. Opt., 14, 29 [Google Scholar]
  4. Bogdanov, S., Lamb, F. K., Mahmoodifar, S., et al. 2019, ApJ, 887, L26 [Google Scholar]
  5. Bonazzola, S., Mottez, F., & Heyvaerts, J. 2015, A&A, 573, A51 [Google Scholar]
  6. Boyd, J. P. 2001, Chebyshev and Fourier Spectral Methods (Springer-Verlag) [Google Scholar]
  7. Chandrasekhar, S. 1970, Ellipsoidal Figures of Equilibrium (New Haven: Yale University Press) [Google Scholar]
  8. Deutsch, A. J. 1955, Annales d’Astrophysique, 18, 1 [Google Scholar]
  9. Finn, L. S., & Shapiro, S. L. 1990, ApJ, 359, 444 [Google Scholar]
  10. Gumerov, N. A., & Duraiswami, R. 2004, in Fast multipole methods for the Helmholtz equation in three dimensions, 1st edn. (Amsterdam: Elsevier), Elsevier series in electromagnetism, oCLC: 249094845 [Google Scholar]
  11. Horedt, G. P. 2004, in Polytropes: Applications in Astrophysics and Related Fields, (Dordrecht, Boston: Kluwer Academic Publishers), Astrophys. Space Sci. Lib., 306 [Google Scholar]
  12. Jackson, J. D. 2001, Electrodynamique classique: Cours et exercices d’electromagnétisme (Paris: Dunod) [Google Scholar]
  13. Leitner, A., & Spence, R. D. 1950, J. Franklin Inst., 249, 299 [Google Scholar]
  14. Li, L.-W., Kang, X.-K., & Leong, M.-S. 2001, Spheroidal Wave Functions in Electromagnetic Theory, 1st edn. (New York: Wiley-Interscience) [Google Scholar]
  15. Makishima, K., Enoto, T., Murakami, H., et al. 2016, PASJ, 68, S12 [Google Scholar]
  16. Makishima, K., Murakami, H., Enoto, T., & Nakazawa, K. 2019, PASJ, 71, 15 [Google Scholar]
  17. Meixner, J., & Schäfke, F. W. 1954, in Mathieusche Funktionen und Sphäroidfunktionen: Mit Anwendungen auf Physikalische und Technische Probleme, eds. J. Meixner, & F. W. Schäfke (Berlin, Heidelberg: Springer), Die Grundlehren der Mathematischen Wissenschaften [Google Scholar]
  18. Morse, P. M., & Feshbach, H. 1953, Methods of Theoretical Physics I (Mcgraw-Hill Book Company Edn) [Google Scholar]
  19. Olver, F. W. J. 2010, NIST handbook of mathematical functions (National Institute of Standards and Technology (U.S.): Cambridge; New York: Cambridge University Press) oCLC: ocn502037224 [Google Scholar]
  20. Pétri, J. 2012, MNRAS, 424, 605 [Google Scholar]
  21. Pétri, J. 2013, MNRAS, 433, 986 [Google Scholar]
  22. Pétri, J. 2014, MNRAS, 439, 1071 [Google Scholar]
  23. Pétri, J. 2015, MNRAS, 450, 714 [Google Scholar]
  24. Pétri, J. 2019, MNRAS, 485, 4573 [Google Scholar]
  25. Riley, T. E., Watts, A. L., Bogdanov, S., et al. 2019, ApJ, 887, L21 [Google Scholar]
  26. 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), 663 [Google Scholar]
  27. Silva, H. O., Pappas, G., Yunes, N., & Yagi, K. 2021, Phys. Rev. D, 103, 063038 [Google Scholar]
  28. Spitkovsky, A. 2006, ApJ, 648, L51 [Google Scholar]
  29. Zanazzi, J. J., & Lai, D. 2015, MNRAS, 451, 695 [Google Scholar]
  30. Zeppenfeld, M. 2009, New J. Phys., 11, 073007 [Google Scholar]

Appendix A: Legendre functions of third type

Legendre polynomials P(x) and their associated functions are usually defined in the interval x ∈ [ − 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 → +∞. These solutions are then called Legendre functions of the third type and related to prolate coordinates for |x|≤1 by analytical continuation in the complex plane, starting from x ∈ [ − 1, 1]. They are solutions of the second order linear differential equation


where and m are two integers. It corresponds to the radial part of a multipole of order (,m) in prolate spheroidal coordinates and related to the spherical harmonic .

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







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 for x > 1, such that


with the correspondence . We note the change in sign in front of the factor x2 of eq. (A.3) compared to eq. (A.1).

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







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 are expanded into Legendre functions , according to


where the prime indicates summation from 0/1 over even and odd indices when (m) is even and odd. The expansion coefficients are determined by solving a three-term recurrence relation with coefficients given in Abramowitz & Stegun (1965). The lowest order corrections in γ2 are




If the subscript is negative, the coefficient vanishes by convention. To this level of approximation, we neglect corrections being of the order γ4 or higher and therefore no corrections apply to the dominant coefficient .

For the expansion of the first multipoles with m = 1, we give the expression of to the γ2 order for  − m = ±2 for m = 1. The first non vanishing coeffifients are listed below





Explicitly, for the first wave functions, this means the following:





All Tables

Table 1.

Fitted coefficients k1 and k2 as given by Eq. (83).

All Figures

thumbnail Fig. 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} in blue, green, and red, respectively. The blue quarter disk in the bottom left depicts the spherical star. The vertical dashed black line represents the light-cylinder.

In the text
thumbnail Fig. 2.

Magnetic field lines in the equatorial plane xOy for a perpendicular spheroidal rotator with oblateness (first two columns) or prolateness (last two columns) parameter a/R = {0, 0.5, 1} in blue, green, and red, respectively. The blue disk in the centre depicts the spherical star. The dashed black circle represents the light-cylinder.

In the text
thumbnail Fig. 3.

Spin-down luminosity for oblate and prolate stars in solid and dashed lines, respectively, with single dipole stellar boundary conditions.

In the text
thumbnail Fig. 4.

Spin-down luminosity for oblate and prolate stars in solid and dashed lines, respectively, with spherical dipole stellar boundary conditions.

In the text
thumbnail Fig. 5.

Polar cap shape for oblate and prolate stars with oblateness parameter a/R = {0, 0.5, 1} in blue, green, and red, respectively. The black dashed line shows the reference solution for the Deutsch field as a check. The obliquity from the left column to the right column is χ = {0°, 30°, 60°, 90°}. The first row is for an oblate star with one mode  = 1, the second row is for a spherical dipole magnetic field at the surface, the third row is for a prolate star with one mode  = 1, and the fourth row is for a spherical dipole magnetic field at the surface.

In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.