Issue 
A&A
Volume 657, January 2022



Article Number  A73  
Number of page(s)  13  
Section  Astrophysical processes  
DOI  https://doi.org/10.1051/00046361/202141911  
Published online  14 January 2022 
Spheroidal magnetic stars rotating in vacuum
Université de Strasbourg, CNRS, Observatoire astronomique de Strasbourg, UMR 7550, 67000 Strasbourg, France
email: jerome.petri@astro.unistra.fr
Received:
30
July
2021
Accepted:
8
October
2021
Context. Gravity shapes stars to become almost spherical because of the isotropic nature of gravitational attraction in Newton’s theory. However, several mechanisms break this isotropy, such as their rotation generating a centrifugal force, magnetic pressure, or anisotropic equations of state. The stellar surface therefore slightly or significantly deviates from a sphere depending on the strength of these anisotropic perturbations.
Aims. In this paper, we compute analytical and numerical solutions of the electromagnetic field produced by a rotating spheroidal star of oblate or prolate nature. This study is particularly relevant for millisecond pulsars for which strong deformations are produced by rotation or a strong magnetic field, leading to indirect observational signatures of the polar cap thermal Xray emission.
Methods. First we solve the time harmonic Maxwell equations in vacuum by using oblate and prolate spheroidal coordinates adapted to the stellar boundary conditions. The solutions are expanded in series of radial and angular spheroidal wave functions. Particular emphasis is put on the magnetic dipole radiation. Second, we compute approximate solutions by integrating the timedependent Maxwell equations in spheroidal coordinates numerically.
Results. We show that the spindown luminosity corrections compared to a perfect sphere are, to leading order, given by terms involving (a/r_{L})^{2} and (a/R)^{2} where a is the stellar oblateness or prolateness, R the smallest star radius, and r_{L} the lightcylinder radius. The corresponding perturbations in the electromagnetic field are only perceptible close to the surface, deforming the polar cap rims. At large distances r ≫ a, the solution tends asymptotically to the perfect spherical case of a rotating dipole.
Key words: methods: analytical / methods: numerical / magnetic fields / stars: general / stars: rotation / pulsars: general
© J. Pétri 2022
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), 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 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:
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 protoneutron 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 rotationinduced 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 R_{pol} and R_{eq} 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 welldefined 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 torquefree 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 pseudospectral timedependent 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 selfgravitating 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
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} + a^{2} cos^{2}ψ. 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 e_{1} and e_{2} such that (with the indices 1 and 2 being ρ, ψ, or φ)
Defined by its covariant components, we get
with ϵ_{ijk} being the spatial LeviCivita 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} + a^{2} sin^{2}ψ. 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 Φ(φ)∝e^{i m φ}, where m is an integer because Φ(φ) must be singlevalued 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 ρ = ±i a z for Eq. (22a) and z = cos ψ for Eq. (22b), both equations reduce to the same SturmLiouville problem
with γ^{2} = −k^{2} a^{2} < 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 s_{2} = 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} = k^{2} a^{2} > 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 curlfree 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 = 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
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 B_{r} 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 + a^{2}/R^{2}). This factor corresponds to an increase in the magnetic field strength, compared to a spherical dipole with surface field B as measured by a distant observer.
3.1.2. Dipole solution
The procedure to follow for the dipole or for any multipole is very similar to the monopole case. Because of the linearity of the problem, we solve 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 spindown 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 quasistationary solution is acceptable for slowly rotating stars or when looking in regions well inside the lightcylinder. However, when seeking for the field behavior at large distances, it switches to a wave nature around the lightcylinder 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 ∝ e^{−i ω 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, Φ = ∇ ∧ (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^{2} Φ.
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 timeharmonic Maxwell equations with temporal behavior e^{−i ω 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} = ±k^{2} a^{2} always remains less than one in modulus because the star must remain within its lightcylinder thus the constrain k a = a/r_{L} ≪ 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 threeterm 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 w^{4} compared to 1 and w^{2} for the dipole. This is due to the fact that the quadrupole is already the result of rotating a magnetic field, thus an w^{2} strength for a w^{2} spindown 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 spindown dependence on a/R and R/r_{L} is approximately proportional to . An expansion to lowest order in these two parameters gives a correction to the luminosity as
where the spindown L_{⊥} of the vacuum orthogonal point magnetic dipole is
The above approximation gives only some important hints about the spindown 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 a^{2}/R^{2}, thus they are of second order in a. We use these results to fit the numerical simulations performed in the next section.
5. Timedependent simulations
Analytically solving the Helmholtz equation with separate coordinates helps to get insight as to the effect of oblateness or prolateness on the spindown 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 timedependent 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 spindown luminosity and eventually tracing the shape of the polar caps.
5.1. Numerical setup
Maxwell equations are solved with our pseudospectral 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/r_{L} = 0.3, which coincides with the inner boundary of the computational domain R_{1} = R. The outer boundary is equal to R_{2}/r_{L} = 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 lightcylinder 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}.
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 lightcylinder. 
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 lightcylinder. 
For the aligned rotator, Fig. 1, we observe the deformation of the surface as a change in the position of the footpoints of the magnetic field lines. Far from the star, especially outside the lightcylinder, 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 twoarmed 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 lightcylinder, 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 spindown luminosity to follow expressions very similar to a spherical star with a dependence on inclination χ such as sin^{2}χ as is shown in the next paragraph.
5.3. Poynting flux
The spindown 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.
Fig. 3.
Spindown luminosity for oblate and prolate stars in solid and dashed lines, respectively, with single dipole stellar boundary conditions. 
Fig. 4.
Spindown 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 spindown 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.
In fact, the coefficient k_{1} 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 k_{1} ≈ 0.919. However, because the outgoing wave boundary conditions stand at a finite radius, relatively close to the lightcylinder, 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 Xray 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 spindown 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 forcefree 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 magnetodipole losses as shown by Pétri (2019).
We infer that the combination of spheroidal geometry and forcefree magnetosphere will lower the dipole magnetic field strength expectation even more. The above fitting can be compared with the forcefree fitting of a spherical star as found by Spitkovsky (2006). A quantitative answer would require computation of forcefree spheroidal magnetospheres, which is out of the scope of the present work. It is difficult to compare Eq. (83) to Spitkovsky (2006) formula because in vacuum we observe only a L_{vac} ≈ L_{⊥} sin^{2}χ dependence, which has to be contrasted with the L_{FFE} ≈ 3/2 L_{⊥} (1+sin^{2}χ) dependence of the forcefree model. Nevertheless we guess that both constant ℓ_{1} and ℓ_{2} in the spheroidal forcefree 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 forcefree 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.
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 Xrays. A careful analysis of such pulsed emission in Xray 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 wellknown 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/r_{L} ≪ 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 Xray emission, modifying their lightcurve in addition to generalrelativistic effects such as framedragging 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.
Acknowledgments
I am grateful to the referee for helpful comments and suggestions. This work has been supported by the CEFIPRA grant IFC/F5904B/2018 and ANR20CE310010.
References
 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]
 Arfken, G. B., & Weber, H.J. 2005, Mathematical Methods for Physicists, 6th edn. (Boston: Elsevier) [Google Scholar]
 Asano, S., & Yamamoto, G. 1975, Appl. Opt., 14, 29 [Google Scholar]
 Bogdanov, S., Lamb, F. K., Mahmoodifar, S., et al. 2019, ApJ, 887, L26 [Google Scholar]
 Bonazzola, S., Mottez, F., & Heyvaerts, J. 2015, A&A, 573, A51 [Google Scholar]
 Boyd, J. P. 2001, Chebyshev and Fourier Spectral Methods (SpringerVerlag) [Google Scholar]
 Chandrasekhar, S. 1970, Ellipsoidal Figures of Equilibrium (New Haven: Yale University Press) [Google Scholar]
 Deutsch, A. J. 1955, Annales d’Astrophysique, 18, 1 [Google Scholar]
 Finn, L. S., & Shapiro, S. L. 1990, ApJ, 359, 444 [Google Scholar]
 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]
 Horedt, G. P. 2004, in Polytropes: Applications in Astrophysics and Related Fields, (Dordrecht, Boston: Kluwer Academic Publishers), Astrophys. Space Sci. Lib., 306 [Google Scholar]
 Jackson, J. D. 2001, Electrodynamique classique: Cours et exercices d’electromagnétisme (Paris: Dunod) [Google Scholar]
 Leitner, A., & Spence, R. D. 1950, J. Franklin Inst., 249, 299 [Google Scholar]
 Li, L.W., Kang, X.K., & Leong, M.S. 2001, Spheroidal Wave Functions in Electromagnetic Theory, 1st edn. (New York: WileyInterscience) [Google Scholar]
 Makishima, K., Enoto, T., Murakami, H., et al. 2016, PASJ, 68, S12 [Google Scholar]
 Makishima, K., Murakami, H., Enoto, T., & Nakazawa, K. 2019, PASJ, 71, 15 [Google Scholar]
 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]
 Morse, P. M., & Feshbach, H. 1953, Methods of Theoretical Physics I (McgrawHill Book Company Edn) [Google Scholar]
 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]
 Pétri, J. 2012, MNRAS, 424, 605 [Google Scholar]
 Pétri, J. 2013, MNRAS, 433, 986 [Google Scholar]
 Pétri, J. 2014, MNRAS, 439, 1071 [Google Scholar]
 Pétri, J. 2015, MNRAS, 450, 714 [Google Scholar]
 Pétri, J. 2019, MNRAS, 485, 4573 [Google Scholar]
 Riley, T. E., Watts, A. L., Bogdanov, S., et al. 2019, ApJ, 887, L21 [Google Scholar]
 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: WileyInterscience), 663 [Google Scholar]
 Silva, H. O., Pappas, G., Yunes, N., & Yagi, K. 2021, Phys. Rev. D, 103, 063038 [Google Scholar]
 Spitkovsky, A. 2006, ApJ, 648, L51 [Google Scholar]
 Zanazzi, J. J., & Lai, D. 2015, MNRAS, 451, 695 [Google Scholar]
 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 realvalued 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 x^{2} 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 realvalued 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 threeterm 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
All Figures
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 lightcylinder. 

In the text 
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 lightcylinder. 

In the text 
Fig. 3.
Spindown luminosity for oblate and prolate stars in solid and dashed lines, respectively, with single dipole stellar boundary conditions. 

In the text 
Fig. 4.
Spindown luminosity for oblate and prolate stars in solid and dashed lines, respectively, with spherical dipole stellar boundary conditions. 

In the text 
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 (fulltext 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 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.