© J. Pétri 2022

1. Introduction

Self-gravitating astrophysical objects are usually set into rotation due to diverse formation and evolution scenarios. Gravitational attraction is compensated for by pressure, and also by centrifugal forces. For fast rotating astrophysical objects, the centrifugal force substitutes a substantial fraction of the internal pressure to sustain the system in a quasi-stationary equilibrium. If, moreover, the gas is magnetised, it radiates an electromagnetic wave, removing energy and angular momentum from the system. Single stars or clouds are typical configurations where such phenomena occur.

This problem is particularly relevant in the realm of high energy astrophysics. Indeed, fast rotating neutron stars are privileged places where the system significantly deviates from a spherical shape, such as for millisecond pulsars. Magnetars are also expected to show strong magnetic constrictions, leading to oblate or even prolate surfaces not directly connected to a geometry imposed by their rotation. A detailed analysis of such systems evolving in vacuum was recently performed by Pétri (2022). If a vacuum is replaced by an ideal plasma such as in electron and positron pairs, a spheroidal force-free neutron star magnetosphere can be constructed, generalising the extensive literature about the force-free magnetosphere of spherical stars (Spitkovsky 2006; Komissarov 2006; McKinney 2006; Pétri 2012; Cao et al. 2016).

From an observational point of view, thermal X-ray emission from neutron star hot spots was used to constrain the compactness and radius of these compact objects (Riley et al. 2019; Bogdanov et al. 2019). Their estimates rely on fitting the X-ray pulsed profile taking general relativity and oblate surfaces into account. However, this modelling requires accurate knowledge of the hot spot shape, its temperature profile, the stellar oblateness, and the general-relativistic environment. All of these ingredients are still uncertain, especially for millisecond pulsars, significantly deviating from a perfect sphere. Nevertheless, these observational hints tend to reliably constrain the hot spot location on the star and their area. In the present work, we start from a theoretical perspective and construct polar cap areas from the knowledge of the electromagnetic field in the force-free regime of a prolate or an oblate star starting from a theoretical point of view.

Computing light curves with general-relativistic effects was done by Cadeau et al. (2007) and Morsink et al. (2007), who used an oblate Schwarzschild approximation to the exterior gravitational field. AlGendy & Morsink (2014) studied the impact of fast rotation on the gravitational field on the surface of a neutron star. They also provide analytical fits to the equatorial radius. Recently Silva et al. (2021) constructed equilibria configurations using realistic neutron star equations of state. Such investigations give a good estimate of the expected stellar deformation due to centrifugal forces and it helps to constrain the oblateness knowing the neutron star rotation period.

The purpose of this paper is to compute force-free neutron star magnetospheres for a spheroidal star surrounded by an ideal pair plasma. In Sect. 2 we present the simulation set-up, solving the time-dependent Maxwell equations using our fully pseudo-spectral code. In Sect. 3, we plot magnetic field lines, spin-down luminosities, polar cap shapes, and current densities. Some discussion about fast rotating neutron stars is proposed in Sect. 4. Conclusions and perspectives are drawn in Sect. 5.

2. Time-dependent simulation set-up

The present study relies on our previous work (Pétri 2022), adding an electric charge and current density in the force-free approximation. We examined the model before considering the salient features of a force-free spheroidal magnetosphere.

The stellar interior was treated as a perfect conductor of spheroidal shape, being oblate or prolate. The time-dependent Maxwell equations were solved in a spheroidal coordinate system with coordinates (ρ, ψ, ϕ) and an oblateness or prolateness parameter a. The coordinate system belongs to one of the eleven separable coordinate systems and forms an orthogonal basis. The electric and magnetic fields E and B were evolved in time according to

$$\begin{array}{cc}& \mathbf{\nabla }·\mathit{B}=0\hfill \end{array}$$(1a)

$$\begin{array}{cc}& \mathbf{\nabla }\times \mathit{E}=-\frac{\partial \mathit{B}}{\partial t}\hfill \end{array}$$(1b)

$$\begin{array}{cc}& \mathbf{\nabla }·\mathit{E}=\frac{{\rho }_{\mathrm{e}}}{{\epsilon }_0}\hfill \end{array}$$(1c)

$$\begin{array}{cc}& \mathbf{\nabla }\times \mathit{B}={\mu }_0\mathit{j}+\frac{1}{c^2}\frac{\partial \mathit{E}}{\partial t}·\hfill \end{array}$$(1d)

We note that ε₀ and μ₀ are the electric permittivity and magnetic permeability, respectively, ρ_e is the charge density, and j is the current density. For an extensive and detailed discussion of the spheroidal coordinate systems, we refer readers to Pétri (2022) in order to avoid reproducing lengthy formulae here.

As usual, in the relativistic force-free regime, the electric current density was deduced from the force-free prescription and it reads as follows:

$$\begin{array}{c}\hfill \mathit{j}={\rho }_{\mathrm{e}}\frac{\mathit{E}\wedge \mathit{B}}{B^2}+\frac{\mathit{B}·\mathbf{\nabla }\times \mathit{B}/{\mu }_0-{\epsilon }_0\mathit{E}·\mathbf{\nabla }\times \mathit{E}}{B^2}\mathit{B}.\end{array}$$(2)

The system of equation was therefore fully determined.

The numerical set-up for performing simulations is similar to the one used in Pétri (2022). The gravitating fluid’s outer boundary is denoted by ρ_in = R and should not to be confused with the spherical radius because we use spheroidal coordinates. The light cylinder radius is r_L and the spheroidal inner surface of the computational domain was set to ρ_in/r_L = 0.3. The outer boundary of the computation domain is located at ρ_out/r_L = 7. We imposed outgoing wave boundary conditions at the outer edge of the simulation domain and imposed tangential electric field component as well as normal magnetic field component continuity on the stellar surface. A numerical convergence analysis showed that a grid resolution of N_ρ × N_ψ × N_φ = 129 × 32 × 64 was sufficient to reckon accurate results.

3. Results

As physical relevant outputs of our simulations, we draw attention to the magnetic field structure, the spin-down luminosity, the polar cap geometry, and the polar cap current density. These quantities are the baselines underlying a more detailed investigation of observational consequences such as thermal X-ray emission as well as radio and gamma-ray pulse profiles. However, the computation of such emission processes are postponed to future works.

3.1. Magnetic field lines

In the force-free regime, the magnetic field lines are attached to the particles. Because outside the light cylinder these particles cannot corotate with the star any more, field lines are opening up, with the poloidal component tending asymptotically to radial lines wound up by rotation. This wind-up is a general feature of any force-free magnetosphere reaching ultra-relativistic corotating speeds close to the light cylinder region and trying to enforce corotation beyond it, irrespective of the exact inner boundary on the stellar surface. We verified this assertion by plotting some of these field lines for spherical, oblate, and prolate stars.

The asymptotic radial structure of poloidal field lines and the transition between closed and open field lines for an aligned rotator are shown quantitatively in Fig. 1 for an oblateness parameter a/R = {0, 0.5, 1} in blue, green, and red, respectively. The blue quarter disk on the bottom left represents the spherical star. The light cylinder is depicted as a vertical black-dashed line. An oblate star inflates in the equatorial direction, whereas a prolate star inflates along the rotation axis. Field lines of spheroidal stars do, therefore, not start right at the blue disk but at a larger distance.

Fig. 1. Magnetic field lines in the meridional plane xOz for an aligned spheroidal rotator with an 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 black-dashed line represents the light cylinder.

Field lines in the equatorial plane for orthogonal rotators are shown in Fig. 2 with the same spheroidal parameter a/R = {0, 0.5, 1} in blue, green, and red, respectively. The blue disk in the centre represents the spherical star. The two-armed spiral shape, typical of a rotating dipole, develops from the light cylinder up to large distances. Its location in space is identical for all simulations.

Fig. 2. Magnetic field lines in the equatorial plane xOy for a perpendicular spheroidal rotator with an 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 black-dashed circle represents the light cylinder.

3.2. Spin-down

The spin-down depends on the surface shape, and also on the normalisation convention used to compute the luminosity. This important issue is discussed in depth by Pétri (2022). In order to get rid of effects that are not directly related to the change in shape, we used a normalisation where the magnetic dipole tends asymptotically to the same magnetic moment value at infinity, irrespective of the oblateness or prolateness parameter a/R.

As for the vacuum case, the spin-down decreases with the increasing parameter a/R, as seen in Fig. 3 for the single dipole field on the surface and in Fig. 4 for the spherical dipole imposed on the surface.

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

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

All simulation runs are summarised by a spin-down fitted formally with

$$\begin{array}{c}\hfill L_{\mathrm{FFE}}^{\mathrm{spheroid}}\approx L_{\perp }({\ell }_1+{\ell }_2{sin}^2\chi ),\end{array}$$(3)

where the vacuum luminosity for a perpendicular rotator is

$$\begin{array}{c}\hfill L_{\perp }=\frac{8\pi }{3{\mu }_0{c}^3}{\mathrm{\Omega }}^4{B}^2{R}^6,\end{array}$$(4)

where B is the magnetic field strength at the equator of a spherical star, Ω is its rotation rate, and R is its radius. The coefficients ℓ_i for i = {1, 2} depend on the ellipticity according to

$$\begin{array}{c}\hfill {\ell }_i={\alpha }_i-{\beta }_i{\left(\frac{a}{R}\right)}^2.\end{array}$$(5)

The fitted values are reported in Table 1. The spheroidal shape always decreases the spin-down luminosity with respect to the spherical star. This effect is most pronounced for prolate stars.

3.3. Polar cap shape

Polar caps are important observables. Indeed, the shape of their rims can be indirectly probed by their thermal X-ray pulsation as detected by the Neutron Star Interior Composition Explorer (NICER), see for instance Riley et al. (2019) and Bogdanov et al. (2019). As a starting point for computing thermal X-ray emission, we present realistic polar cap surfaces relying on our force-free simulations of spheroidal stars in this paper. Figure 5 compiles the geometry of the polar cap edges for various neutron star surface deformations, being either oblate or prolate, different inclination angles χ = {0° ,30° ,60° ,90° }, and several normalised spheroidal parameters a/R = {0, 0.5, 1}. For reference, the vacuum Deutsch solution is also shown in orange-dashed lines.

Fig. 5. Polar cap shape for oblate and prolate stars with the oblateness parameter a/R = {0, 0.5, 1} in blue, green, and red, respectively. The orange-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 multipole of order ℓ = 1, the second row is for a spherical dipole magnetic field at the surface, the third row is for a prolate star with one multipole of order ℓ = 1, and the fourth row is for a spherical dipole magnetic field at the surface.

Compared to polar caps deduced from vacuum, in force-free magnetospheres, the cap rims always lie outside the equivalent Deutsch case. For the aligned rotator, the size and area of the cap, being perfectly circular due to the axisymmetry, increases with oblateness or prolateness, except for the third row, as seen in the left column of Fig. 5. For large obliquities, the polar cap surface area also increases with the spheroidal parameter a. Nevertheless, for prolate stars, the deformation of the polar cap remains weak because along the equator, the stellar shape remains almost spherical, tending to the spherical force-free results. For oblate stars, the situation is the opposite because the maximum deformation arises around the equator and the polar cap area increases significantly with oblateness. Contrary to the spin-down comparison performed in the previous section where some fiducial spherical model had to be chosen, the polar cap rims solely rely on geometrical effects of magnetic field lines. Therefore, the polar surfaces discussed in this section are insensitive to the particular normalisation used for the magnetic dipole. The results are robust and can be applied to any magnetic field intensity.

3.4. Polar cap current

Just as for the previous quantities, the current produced within the magnetosphere is significantly altered by the spheroidal shape of the stellar surface. The local current density expelled and returning to the polar caps offers an important view of the electrodynamics impacted by the spheroidal star. This current density is best reckoned in the star corotating frame, removing the convective current produced by the electric drift motion of the charge-separated plasma component. This current is not directly related to the conduction current which is the primordial current to plot in force-free electrodynamics. Endean (1974) presented an interesting derivation of the relativistic force-free equation involving only the magnetic field, therefore generalising the usual non-relativistic equation

$$\begin{array}{c}\hfill (\mathbf{\nabla }\wedge \mathit{B})\wedge \mathit{B}/{\mu }_0=\mathit{j}\wedge \mathit{B}=\mathbf{0}\end{array}$$(6)

to a relativistic version given by

$$\begin{array}{c}\hfill (\mathbf{\nabla }\wedge {\mathit{B}}^{\ast })\wedge \mathit{B}/{\mu }_0=\mathbf{0},\end{array}$$(7)

where he introduced a new vector B^* defined by

$$\begin{array}{c}\hfill {\mathit{B}}^{\ast }=((1-\frac{x^2+y^2}{r_{\mathrm{L}}^2})B^{\rho },(1-\frac{x^2+y^2}{r_{\mathrm{L}}^2})B^{\psi },B^{\varphi }).\end{array}$$(8)

A similar expression was given by Mestel (1973) at the same time and also explained by a simple Lorentz transformation. Therefore, according to Eq. (7), the term

$$\begin{array}{c}\hfill {\mu }_0{\mathit{j}}_{\mathrm{c}}=\mathbf{\nabla }\wedge {\mathit{B}}^{\ast }\end{array}$$(9)

is directed along B and can be interpreted as the conduction current density j_c in the corotating frame (see also Bai & Spitkovsky 2010 for another derivation of this interpretation). For numerical purposes, we normalised the current density to a fiducial value related to the corotating current density at the pole of an aligned rotator and given by

$$\begin{array}{c}\hfill j_{\mathrm{pol}}=2{\epsilon }_0\mathrm{\Omega }Bc.\end{array}$$(10)

Following the definition given by Eq. (9), Fig. 6 shows the conduction current density on the surface of an aligned oblate star with the oblateness parameter a/R = {0, 0.5, 1} for one pole with ψ ∈ [0° ,90° ]. The other pole is not shown because of the north-south antisymmetry. The conduction current j_c is actually symmetric with respect to the equator located at a colatitude of ψ = 90°. The oblateness decreases the current density by up to a factor of two for a/R = 1. Nevertheless, the polar cap area also increases with oblateness, partially compensating for the current density decrease. The current flowing in the vicinity of the polar cap, therefore, crucially depends on the oblateness. Combining the variation in the area, we expect the associated thermal X-ray emission to be drastically influenced by the stellar surface geometry in addition to the polar cap area and shape.

Fig. 6. Conduction current density on the surface of an aligned oblate star with the oblateness parameter a/R = {0, 0.5, 1}.

Figure 7 shows the same conduction current density on the surface of an aligned prolate star with the oblateness parameter a/R = {0, 0.5, 1}. Here, the current density also decreases, but to a lesser extend compared to the oblate case. The impact on the thermal X-ray emission is consequently less affected than in the previous oblate case.

Fig. 7. Conduction current density on the surface of an aligned prolate star with the oblateness parameter a/R = {0, 0.5, 1}.

To be exhaustive, we also computed the current density for orthogonal rotators in oblate and prolate geometries. Figure 8 shows a cross section, taken around the centre of the polar cap, of the conduction current in oblate stars with the oblateness parameter a/R = {0, 0.5, 1}. The conduction current j_c is actually skew-symmetric with respect to the equator located at a colatitude of ψ = 90°. Eventually, Fig. 9 shows the same cross section for a prolate star with the prolateness parameter a/R = {0, 0.5, 1}. We always notice a monotonic decrease in the peak current density with an increasing spheroidal parameter a.

Fig. 8. Conduction current density on the surface of an orthogonal oblate star with the oblateness parameter a/R = {0, 0.5, 1}.

Fig. 9. Conduction current density on the surface of an orthogonal prolate star with the prolateness parameter a/R = {0, 0.5, 1}.

4. Discussion

Nowadays, thermal X-ray pulsation emanating from the surface of neutron stars is used to probe the location and shape of the underlying hot spots. Recent observations by NICER have clearly demonstrated the strength of such approaches in constraining the mass, the equatorial radius (Riley et al. 2019), and the equation of state (Miller et al. 2019) of these compact objects. However, these investigations rely on observational fittings that are not always quantitatively and accurately linked to neutron star magnetosphere models. For instance, we are not aware of any oblate force-free magnetosphere computation that shows the impact of oblateness on spin-down efficiency and polar cap perturbations induced by the deviation from a spherical star.

Our results unambiguously predict the importance of carefully estimating the oblateness in order to fit hot spot X-ray light curves. Irrespective of the normalisation employed for the magnetic field, we always found an increase in the hot spot surface area compared to a spherical star. The current flowing back to the star, hitting and heating the surface, is therefore also modified by the oblateness or prolateness. As a consequence, the surface brightness temperature and the luminosity of these hot spots are influenced by the stellar boundary, not only because of geometrical effects but also because of electrodynamic effects related to the change in current and charge density existing inside the magnetosphere.

Another important ingredient not included in the present study is general relativity. Because of the large compactness of neutron stars, curved space time also significantly modifies the electrodynamics of its magnetosphere, especially in the vicinity of its surface. Unfortunately, this requires solutions to Einstein equations for an oblate gravitating fluid, and simple analytical formulas for the exterior space time are difficult to find with the exception of Erez & Rosen (1959), corrected by Young & Coulter (1969) and Doroshkevich et al. (1965). A very general solution for the exterior of the star has been proposed by Quevedo (1989). The multipole components are, however, left as free parameters not straightforwardly connected to the stellar oblateness. Nevertheless, a simple and good approximation is obtained by the so-called oblate+Schwarzschild model where the oblate stellar surface is embedded into the Schwarzschild metric. It is simple, but it provides a good impression of these additional effects even if it does not include frame-dragging. Last but not least, it was shown that magnetic multipole components (Bilous et al. 2019) can play a central role in understanding millisecond pulsar light curves (Kalapotharakos et al. 2021).

The oblateness parameter a can be deduced from numerical models of neutron star interiors. According to Table 1 of Silva et al. (2021), in using realistic equations of state for the construction of stars in stationary equilibrium for fast rotating neutron stars with a frequency of about 600 Hz, the ratio between the polar radius R_p and equatorial radius R_e is about r = R_p/R_e ≈ 0.9, implying an oblateness parameter of $a=\sqrt{r^{-2}-1}\approx 0.48$. The deformation of the stellar surface is therefore substantial, leading to a slight increase in the polar cap surface area.

The impact on the spin-down luminosity depends on the internal structure of the magnetic field compared to the equivalent spherical star. This issue is, however, not solved. For some neutron stars, a birth period around 1 ms is expected, even for strongly magnetised ones, coined millisecond magnetars and suspected to be the central engine of hypernovae or superluminous supernovae. For such high periods, the aspect ratio is about r ≈ 0.7, implying an oblateness of a ≈ 1. If we assume an ideal plasma compression due to the star slowing down, the magnetic flux is approximately conserved and the magnetic field strength at the equator increases by a factor of 1/r ≈ 1.42. The associated spin-down augments by a factor of 1/r² ≈ 2.0. Therefore, at least two mechanisms contribute to a modification of the neutron star luminosity, first the change in the stellar surface and second the variation in the surface magnetic field strength. We believe that these effects must be taken into account for the evolution of the early stage of a neutron star relaxing to a spherical shape as it slows down.

5. Conclusions

Although neutron stars are subject to enormous gravitational binding energy, millisecond rotation periods can lead to significant centrifugal forces deforming their surface into an oblate shape or even into a prolate shape due to magnetic compression. We computed force-free magnetospheres for spheroidal neutron stars and showed the repercussion on magnetic topology, spin-down efficiency, polar cap geometry, and the current density. We observed significant changes in the polar cap rims for reasonable spheroidal parameters. Compared to a spherical star, the impact on spin-down luminosity remains small with our normalisation convention. However, the likelihood of the magnetic field geometry having repercussions on the polar cap shape is significant, irrespective of the magnetic field strength. Such shapes can indirectly be probed via thermal X-ray emission from hot spots, as recently demonstrated by the NICER Collaboration.

From our study, we can compute such X-ray light curves, knowing the size and shape of the hot spot. Nevertheless, in order to better stick to realistic neutron stars, we need to take the compactness of the star into account and, therefore, we are obliged to include some general-relativistic effects using, for instance, an oblate Schwarzschild metric or more accurate models of spheroidal neutron star gravitational fields. These investigations are the natural steps that we pursue in forthcoming studies.

Acknowledgments

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

References

All Tables

All Figures

	Fig. 1. Magnetic field lines in the meridional plane xOz for an aligned spheroidal rotator with an 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 black-dashed line represents the light cylinder.
In the text

	Fig. 2. Magnetic field lines in the equatorial plane xOy for a perpendicular spheroidal rotator with an 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 black-dashed circle represents the light cylinder.
In the text

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

	Fig. 4. Spin-down luminosity for oblate and prolate stars and the associated fits 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 the oblateness parameter a/R = {0, 0.5, 1} in blue, green, and red, respectively. The orange-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 multipole of order ℓ = 1, the second row is for a spherical dipole magnetic field at the surface, the third row is for a prolate star with one multipole of order ℓ = 1, and the fourth row is for a spherical dipole magnetic field at the surface.

In the text

	Fig. 6. Conduction current density on the surface of an aligned oblate star with the oblateness parameter a/R = {0, 0.5, 1}.
In the text

	Fig. 7. Conduction current density on the surface of an aligned prolate star with the oblateness parameter a/R = {0, 0.5, 1}.
In the text

	Fig. 8. Conduction current density on the surface of an orthogonal oblate star with the oblateness parameter a/R = {0, 0.5, 1}.
In the text

	Fig. 9. Conduction current density on the surface of an orthogonal prolate star with the prolateness parameter a/R = {0, 0.5, 1}.
In the text