The pulsar magnetosphere with machine learning: methodology

ABSTRACT

In this study, we introduce a novel approach for deriving the solution of the ideal force-free steady-state pulsar magnetosphere in three dimensions. Our method involves partitioning the magnetosphere into the regions of closed and open field lines, and subsequently training two custom Physics Informed Neural Networks to generate the solution within each region. We periodically modify the shape of the boundary separating the two regions (the separatrix) to ensure pressure balance throughout. Our approach provides an effective way to handle mathematical contact discontinuities in force-free electrodynamics. We present preliminary results in axisymmetry, which underscore the significant potential of our method. Finally, we discuss the challenges and limitations encountered while working with Neural Networks, thus providing valuable insights from our experience.

magnetic fields, pulsars: general, methods: numerical

1 INTRODUCTION

55 yr after the discovery of pulsars (Hewish et al. 1968), there still does not exist a complete physical model for their strongly magnetized plasma-filled magnetospheres. The ground for their modelling was set by Goldreich & Julian (1969), who suggested that one part of the magnetosphere rigidly corotates with the central star (the ‘closed line region’ or ‘dead zone’), while the parts around the magnetic poles (the ‘polar caps’) open up to infinity beyond the light cylinder. An electric field generated by the stellar rotation fills the magnetosphere with electric charges and the open field lines with electric currents. For 30 years, everyone’s focus was on the light cylinder, which was considered to be the source of magnetospheric dissipation. Our understanding changed dramatically when the first solution of the force-free ideal aligned pulsar magnetosphere was obtained (Contopoulos, Kazanas & Fendt 1999; hereafter CKF). It was clearly proposed then that magnetospheric dissipation most probably takes place along the large-scale magnetospheric electric current sheet that the new solution clearly identified, and not along the light cylinder.

CKF introduced a novel iterative numerical method that relaxed to the steady-state solution of the axisymmetric pulsar magnetosphere. The non-axisymmetric problem, however, was solved with a different approach. Most studies of 3D pulsar magnetospheres are performed with time-dependent numerical simulations that start with a dipole magnetic field that is set into rotation at t = 0. The simulations run for long enough times to attain a rotating steady-state configuration which is then presented as the solution of the pulsar magnetosphere. This has been the approach of the pioneering FFE (force-free electrodynamics) simulations of Spitkovsky (2006) and Kalapotharakos & Contopoulos (2009), the MHD (magnetohydrodynamics) simulations of Tchekhovskoy, Spitkovsky & Li (2013), the ‘ab initio’ PIC (particle-in-cell) simulations of Philippov & Spitkovsky (2014), Philippov, Spitkovsky & Cerutti (2015), and Kalapotharakos et al. (2018), and most recently, the hybrid PIC–MHD simulations of Soudais & Cerutti (2022). These solutions have confirmed that the bulk of the magnetosphere operates under ideal force-free conditions, while electromagnetic dissipation, particle acceleration and high-energy radiation take place mostly along the equatorial current sheet. The steady-state problem was also solved with spectral methods (Parfrey, Beloborodov & Hui 2012; Pétri 2012) which, unfortunately, behave very poorly in the equatorial current sheet.

Time-dependent three-dimensional (3D) numerical simulations have evolved dramatically in the past 10 years. Unfortunately, there exists no reference 3D ideal force-free solution with which to compare and validate their results. Thus, several important issues are left unanswered and new problems arose:

In all PIC simulations, the magnetosphere opens-up a significant distance inside the light cylinder (their so-called Y-point may appear down to 85 per cent of the light cylinder radius; e.g. Hu & Beloborodov 2022). If true, this would dramatically modify the pulsar spindown rate, the pulsar braking index, as well as the shape and spectrum of the high-energy pulsar radiation originating around the Y-point. We do not understand the physical origin of this effect in PIC simulations. We do not agree with the common explanation that it is due to the artificially high inertia of the PIC super-particles. Instead, we suspect that it is related to a subtle minimum of the magnetospheric electromagnetic energy (Contopoulos, Ntotsikas & Gourgouliatos 2024). We still need to understand where the Y-point truly lies, and how fast it moves out as the pulsar spins down.
One major complication of pulsar magnetospheres is the appearance of electric current sheets inside which the ideal MHD approximation breaks down. Current sheets exist as mathematical contact-discontinuities across which B² − E² is everywhere continuous (see below). MHD/FFE methods smooth out such discontinuities across several computational grid cells. Unfortunately, ideal FFE conditions are not valid inside current sheets. Moreover, different MHD numerical methods treat the equatorial current sheet differently (e.g. CKF and Timokhin 2006 erroneously place it inside the closed line region), and it is not possible in general to differentiate between artificial (numerical) and true (physical) features in the published solutions (e.g. interchanging multiple positive and negative charge layers as in Kalapotharakos et al. 2012a; Mahlmann & Aloy 2022).Unfortunately, there does not exist a reference solution of the ideal (dissipation-less) force-free 3D problem against which to evaluate our numerical simulations. A reference solution exists only in 2D (CKF; Timokhin 2006).
Recent global PIC simulations (Cerutti, Philippov & Dubus 2020; Hakobyan, Philippov & Spitkovsky 2023) show that the electromagnetic (Poynting) flux remaining in the pulsar magnetosphere decreases as the logarithm of the distance r over the radius of the light cylinder R_LC, namely ∝ 1 − β_recln (r/R_LC). With β_rec ∼ 0.1 as obtained from local PIC simulations of relativistic reconnection layers (e.g. Lyubarsky 2005; Werner et al. 2018), this results in the full pulsar magnetosphere completely disappearing via dissipation in the equatorial current sheet within a few hundred light cylinder radii. According to Hakobyan et al. (2023), this logarithmic gradual electromagnetic energy dissipation applies to the full magnetosphere, not only to its undulating non-axisymmetric component as Cerutti et al. (2020) claim. This leaves no MHD pulsar wind to extend to the termination shock to power the X-ray pulsar nebula. We do not understand this result.
Time-dependent simulations relax to one final steady-state solution. In nature, however, pulsars often exhibit mode-switching which cannot be accounted for by a unique solution. We must look for alternative global pulsar magnetospheric solutions, and we need to understand why and how pulsars may spontaneously jump between them (e.g. Ntotsikas et al. 2024).
The resolution of current PIC simulations is inadequate by several orders of magnitude to properly model the microphysics of the equatorial current sheet. In order to generate pulsar light curves and spectra that may be compared with observations, the simulation results (particle Lorentz factors, magnetic field values, electromagnetic spectra) are extrapolated by several orders of magnitude. Unfortunately, there is no agreement among different research groups on the particular method of extrapolation, and as a result, there is still no generally accepted understanding of the physical origin of the high-energy radiation from pulsars (e.g. synchrotron as in Hakobyan et al. 2023; gamma-ray fundamental plane as in Kalapotharakos et al. 2022; inverse Compton as in Richards & Lyutikov 2018). The remedy would be to investigate the trajectories of accelerated particles for realistic physical parameters (not PIC parameters) in the current sheet obtained in a reference solution.

We believe that it is probably too early to derive safe conclusions about the pulsar magnetosphere from current global (‘ab initio’) PIC simulations, thus we can only trust them qualitatively (not quantitatively) to make meaningful comparisons with observations. Their parameters are orders of magnitude away from realistic values, and except for the aligned (non-pulsar) rotator, there is currently no way to validate their results by comparison with a reference solution. Given the inherent problems of all numerical methods discussed above, we believe that this is a dangerous development that may strongly mislead pulsar research. We thus propose to return to the basics and independently obtain the reference ideal force-free magnetosphere with a novel numerical method. In this paper, we will present our first solutions of the axisymmetric problem, but our results will be directly generalized in the case of the 3D inclined rotator in a forthcoming publication. We present the general mathematical formulation of the problem in Section 2. In Section 3, we propose a novel numerical approach with machine learning that will allow us to obtain the solution in 3D. In Section 4, we formulate the problem in the axisymmetric case of the aligned rotator and present our first solutions of the pulsar equation with a trained Physics Informed Neural Network (PINN). In Section 5, we summarize our results and discuss the challenges and limitations encountered while working with Neural Networks (NNs).

2 GENERAL MATHEMATICAL FORMULATION

The magnetospheres of isolated neutron stars are generally considered to be dominated by the magnetic field. The physical conditions in the magnetosphere, i.e. the low mass of the charge carriers and the low density of the magnetospheric plasma, allow us to neglect gravity, thermal pressure and particle inertia as they are several orders of magnitude smaller than the electromagnetic forces. Particle inertia may only become important at the tip of the closed line region very close to the light cylinder. Other than that, force balance in the bulk of the pulsar magnetosphere is reduced to

$$\begin{eqnarray} &&\rho _{\rm e}{\boldsymbol E}+{\boldsymbol J}\times {\boldsymbol B}/c=0\ , \end{eqnarray}$$

(1)

where B and E are the magnetic and the electric field, respectively, ρ_e is the electric charge density, and J is the electric current

$$\begin{eqnarray} {\bf J} &\equiv & \frac{c}{4\pi }\nabla \cdot {\boldsymbol E}\ \frac{{\boldsymbol E}\times {\boldsymbol B}}{B^2}+ \frac{c}{4\pi }\frac{{\boldsymbol B}\cdot \nabla \times {\boldsymbol B}-{\boldsymbol E}\cdot \nabla \times {\boldsymbol E}}{B^2}\ {\boldsymbol B}\ . \end{eqnarray}$$

(2)

The above general expression was obtained by Gruzinov (1999) and Blandford (2002), and simplifies considerably in steady state.

All fields are calculated in the non-rotating inertial lab frame. In that frame, electromagnetic fields need to satisfy Maxwell’s equations, namely

$$\begin{eqnarray} &&\nabla \cdot {\boldsymbol B}=0\ ,\ \nabla \times {\boldsymbol B} = \frac{1}{c}\frac{\partial {\boldsymbol E}}{\partial t}+\frac{4\pi }{c}{\boldsymbol J}\ , \\ && \nabla \cdot {\boldsymbol E}= \frac{4\pi }{c}\rho _{\rm e}\ ,\ \nabla \times {\boldsymbol E} = -\frac{1}{c}\frac{\partial {\boldsymbol B}}{\partial t}\ . \end{eqnarray}$$

(3)

We will be searching for the steady-state solution to the above set of equations. Following the approach of Muslimov & Harding (2009), we will define the transformation of partial time derivatives between our lab frame, and a mathematical (not physical) reference frame with the pulsar angular velocity Ω around the axis θ = 0, namely

$$\begin{eqnarray} \frac{\partial {\boldsymbol B}}{\partial t} & = & \left.\frac{\partial {\boldsymbol B}}{\partial t}\right|_{\rm corot}- \nabla \times (r\sin \theta \ \Omega \ \hat{\phi }\times {\boldsymbol B})\ , \\ \frac{\partial {\boldsymbol E}}{\partial t} & = & \left.\frac{\partial {\boldsymbol E}}{\partial t}\right|_{\rm corot}- \nabla \times (r\sin \theta \ \Omega \ \hat{\phi }\times {\boldsymbol E}) +r\sin \theta \ \Omega \ \hat{\phi }\ \nabla \cdot {\boldsymbol E}\ , \\ \frac{\partial \rho _{\rm e}}{\partial t} & = & \left.\frac{\partial \rho _{\rm e}}{\partial t}\right|_{\rm corot} +r\sin \theta \ \Omega \ \hat{\phi }\cdot \nabla \rho _{\rm e}\ . \end{eqnarray}$$

(4)

Time derivatives in that frame (denoted by the subscript ‘corot’) vanish in steady state, and therefore, the Maxwell’s equations that describe the steady state of the force-free pulsar magnetosphere become

$$\begin{eqnarray} \nabla \times \left({\boldsymbol E}+\frac{r\sin \theta \ \Omega }{c}\ \hat{\phi }\times {\boldsymbol B}\right) &=& 0\ , \end{eqnarray}$$

(5)

$$\begin{eqnarray} \nabla \times \left({\boldsymbol B}-\frac{r\sin \theta \ \Omega }{c}\ \hat{\phi }\times {\boldsymbol E}\right) &=& \frac{4\pi }{c}{\boldsymbol J}-\frac{r\sin \theta \ \Omega }{c}\ \hat{\phi }\ \nabla \cdot {\boldsymbol E}\ . \end{eqnarray}$$

(6)

Let us define here the light cylinder radius R_LC ≡ c/Ω. From equation (5) above we obtain that

$$\begin{eqnarray} &&{\boldsymbol E}\equiv {\boldsymbol E}_{\rm p}=-\frac{r\sin \theta }{R_{\rm LC}}\ \hat{\phi }\times {\boldsymbol B}= \frac{r\sin \theta }{R_{\rm LC}}\ (B_\theta ,-B_r,0) \end{eqnarray}$$

(7)

in component form, where the subscript ‘p’ denotes a poloidal component. Following equation (7), E · ∇ × E = 0, and the expression for J in equation (2) simplifies considerably. Substituting everything back in equation (6) we obtain ∇ × {B_p (1 − (rsin θ/R_LC)²) + B_ϕ} × B = 0 or equivalently,

$$\begin{eqnarray} && \nabla \times \left\lbrace {\boldsymbol B}_{\rm p}\left(1-\left(\frac{r\sin \theta }{R_{\rm LC}}\right)^2\right)+{\boldsymbol B}_\phi \right\rbrace =\alpha {\boldsymbol B}\ . \end{eqnarray}$$

(8)

Here, α is the force-free parameter that obeys the constraint

$$\begin{eqnarray} &&{\boldsymbol B}\cdot \nabla \alpha =0\ , \end{eqnarray}$$

(9)

namely that α is constant along magnetic field lines. In other words, field lines lie along iso-contours of α in 3D. Notice that solving for the steady-state configuration in the co-rotating frame has an important advantage over time-dependent simulations in the non-rotating (lab) frame. In the latter, the final configuration is rotating, i.e. time dependent. All its complex features (current sheets, Y-points, etc.) rotate in the simulation frame of reference; hence, it is difficult to treat them and to guarantee that a steady state is truly reached. On the other hand, solving for the steady-state configuration in the co-rotating frame is a relaxation approach that allows a better treatment of current sheets and more naturally reaches the final steady state. Equation (8) has first been obtained by Endean (1974) and Mestel (1975) but has never been solved before in 3D.

We now need to introduce the vector potential A such that B = ∇ × A. We will work in spherical coordinates (r, θ, ϕ) centred on to the central star. Our notation will be slightly simplified if we introduce the vector magnetic flux components

$$\begin{eqnarray} &&\Psi _r\equiv rA_r\ ,\ \Psi _\theta \equiv rA_\theta \ ,\ \Psi _\phi \equiv r\sin \theta \ A_\phi \ . \end{eqnarray}$$

(10)

In that notation,

$$\begin{eqnarray} B_r &\equiv & \frac{1}{r^2\sin \theta }\left(\frac{\partial \Psi _\phi }{\partial \theta }-\frac{\partial \Psi _\theta }{\partial \phi }\right)\ , \\ B_\theta &\equiv & \frac{1}{r^2\sin \theta }\left(\frac{\partial \Psi _r}{\partial \phi }-r\frac{\partial \Psi _\phi }{\partial r}\right)\ , \\ B_\phi &\equiv & \frac{1}{r}\frac{\partial \Psi _\theta }{\partial r}-\frac{1}{r^2}\frac{\partial \Psi _r}{\partial \theta }\ , \end{eqnarray}$$

(11)

and equation (8) becomes

$$\begin{eqnarray} &&\frac{\partial ^2\Psi _r}{\partial \phi ^2}-r\frac{\partial ^2\Psi _\phi }{\partial r\partial \phi }=\frac{r\sin \theta }{1-(r\sin \theta /R_{\rm LC})^2} \left\lbrace -\alpha \left(\frac{\partial \Psi _\phi }{\partial \theta }-\frac{\partial \Psi _\theta }{\partial \phi }\right)\right. \\ &&+ \cos \theta \ \frac{\partial \Psi _\theta }{\partial r}-\frac{\cos \theta }{r}\frac{\partial \Psi _r}{\partial \theta } +\left. \sin \theta \ \frac{\partial ^2\Psi _\theta }{\partial \theta \partial r}-\frac{\sin \theta }{r}\frac{\partial ^2\Psi _r}{\partial \theta ^2}\right\rbrace \ , \end{eqnarray}$$

(12)

$$\begin{eqnarray} &&\frac{\partial ^2\Psi _\theta }{\partial \phi ^2}- r\frac{\partial ^2\Psi _\phi }{\partial \theta \partial \phi }= \frac{r\sin \theta }{1-(r\sin \theta /R_{\rm LC})^2} \left\lbrace -\alpha \left(\frac{\partial \Psi _r}{\partial \phi }-r\frac{\partial \Psi _\phi }{\partial r}\right)\right. \\ && -\left. r\sin \theta \ \frac{\partial ^2\Psi _\theta }{\partial r^2} -\frac{\sin \theta }{r}\frac{\partial \Psi _r}{\partial \theta }+\sin \theta \ \frac{\partial ^2\Psi _r}{\partial \theta \partial r}\right\rbrace \ , \end{eqnarray}$$

(13)

$$\begin{eqnarray} &&\frac{\partial ^2\Psi _\phi }{\partial r^2}+\frac{1}{r^2}\frac{\partial ^2\Psi _\phi }{\partial \theta ^2}- \frac{1}{r^2}\frac{\partial ^2\Psi _\theta }{\partial \theta \partial \phi } -\frac{1}{r}\frac{\partial ^2\Psi _r}{\partial r\partial \phi } \\ && +\frac{1}{r^2}\frac{\partial \Psi _r}{\partial \phi } - \frac{\cos \theta }{r^2 \sin \theta }\left(\frac{\partial \Psi _\phi }{\partial \theta }- \frac{\partial \Psi _\theta }{\partial \phi }\right)= \\ && \frac{\sin \theta }{1-(r\sin \theta /R_{\rm LC})^2}\left\lbrace -\alpha \left(\frac{\partial \Psi _\theta }{\partial r}-\frac{1}{r}\frac{\partial \Psi _r}{\partial \theta }\right)\right. \\ && +\left. \frac{2\cos \theta }{R^2_{\rm LC}}\left(\frac{\partial \Psi _\phi }{\partial \theta } -\frac{\partial \Psi _\theta }{\partial \theta }\right) -\frac{2\sin \theta }{R^2_{\rm LC}}\left(\frac{\partial \Psi _r}{\partial \phi } -r\frac{\partial \Psi _\phi }{\partial r}\right) \right\rbrace . \end{eqnarray}$$

(14)

From the rhs of equation (14), the reader can check directly that the general regularization condition at the light cylinder rsin θ = R_LC becomes

$$\begin{eqnarray} \alpha =\left. 2\ \frac{\cos \theta \ B_r-\sin \theta \ B_\theta }{B_\phi R_{\rm LC}}\right|_{\rm LC}\equiv \left.\frac{2B_z}{B_\phi R_{\rm LC}}\right|_{\rm LC}\ . \end{eqnarray}$$

(15)

Equation (15) determines the value of the function α along all field lines that cross the light cylinder and never return to the star. In an untwisted magnetosphere, all other field lines that do not cross the light cylinder and form the closed line region have α = 0.

3 NOVEL NUMERICAL APPROACH

We propose to solve the steady-state pulsar magnetosphere problem using Machine Learning techniques, and in particular by training a custom NN. In Sections 3.2 and 3.3, we will introduce two innovations that will lead to a better solution than what was obtained before both in 2D and 3D. Our formulation is general, but our method will also be implemented in axisymmetry in the next section.

3.1 Machine learning

A NN that is used to solve the partial differential equations that describe a physical problem is called a PINN. Quoting Wikipedia, ‘PINNs are a type of universal function approximators that can embed the knowledge of any physical laws that govern a given data-set in the learning process, and can be described by partial differential equations’. The entries of such NN are basically the spatial coordinates r, θ, ϕ, and the exits are Ψ_r, Ψ_θ, Ψ_ϕ, and α. The NN consists of several internal layers (such a NN is termed a deep NN) with several tens of nodes each. Obviously, the exits of the NN are functions of the entries r, θ, ϕ; hence, all order partial derivatives of the exits with respect to the entries are known. The fact that the derivatives of the fields are known locally is an important property of PINNs that differentiates them from finite difference grid methods.

The first implementation of PINNs in pulsar magnetospheres were obtained in the pioneering works of Stefanou, Pons & Cerdá-Durán (2023a) and Stefanou et al. (2023b), where one can find a detailed description of their construction and training. PINNs are trained with several loss functions that implement the constraints of the physical problem and are added together into a single loss function. The objective of the training is to adjust the internal parameters (weights) of the NN that will minimize the total loss calculated over random points in the computational domain, i.e. bring it as close to zero as possible everywhere. In other words, we choose random points over the computational domain r_* ≤ r ≤ r_max, 0 ≤ θ ≤ π, 0 ≤ ϕ ≤ 2π (here, r_* is the radius of the central neutron star) and try to satisfy all the following constraints on them:

Equations (12)–(14) over the computational domain
Equations (9) over the computational domain
Equation (15) at rsin θ = R_LC (it may not be necessary to impose this condition since it is satisfied if (i) is satisfied)
α = 0 along field lines that do not cross the light cylinder (we will assume here an untwisted closed line region)
Boundary magnetic field conditions on the stellar surface. If we consider a central dipole magnetic field on a star of radius r_* with the axis of the dipole along some direction |$\hat{m}$| with a star-centred spherical coordinate system (r, θ_m, ϕ_m) along that direction, then
$$\begin{eqnarray} \Psi _r(r_*,\theta _{\rm m},\phi _{\rm m}) &=& \Psi _{\theta _{\rm m}}(r_*,\theta _{\rm m},\phi _{\rm m}) = 0, \\ \Psi _{\phi _{\rm m}}(r_*,\theta _{\rm m},\phi _{\rm m}) &=& \Psi _{\rm max}\sin ^2\theta _{\rm m}. \end{eqnarray}$$
(16)
From the above star-centred components |$\Psi _r,\Psi _{\theta _m},\Psi _{\phi _m}$| of the vector potential on the stellar surface, we will obtain the actual components Ψ_r, Ψ_θ, Ψ_ϕ of the vector potential on the stellar surface as boundary conditions of our problem.

One nice thing about NNs is that we can consider any number of conditions that the physical problem requires and, given enough internal layers and nodes, training points and training steps, the NN will manage to satisfy all of them by minimizing a global loss function which is just the sum of the various constraints. Note that boundary conditions are also introduced in the form of extra constraints. For example, the loss functions that need to be minimized to zero and describe equation (16) are |$|\Psi _r|, |\Psi _{\theta _m}|,$| and |$|\Psi _{\phi _m}-\Psi _{\rm max}\sin ^2\theta _m|$| summed over the (r_*, θ_m, ϕ_m) training points chosen randomly along the stellar surface. One other condition that needs to be satisfied everywhere is E < B. We have checked that this condition is automatically satisfied when magnetic field lines that cross the light cylinder open up to infinity, thus it is not necessary to impose it through a loss function term.

In our present PINN implementation, we used the PyTorch Deep Learning library. Code development was implemented in python on Anaconda Jupyter Notebooks running on local GPUs, while production runs were performed in the Cloud in Google’s Colab. What are the optimal number of internal nodes and training steps, the optimal activation functions and NN optimization parameters are not known a priori, and are determined after several trial trainings of the problem (see one particular implementation in Section 4.5). We implemented Adam optimizers, SiLU activation functions, and ReduceLROnPlateau schedulers from the PyTorch library. Another nice thing about NNs is that they are mesh-less, and the distribution of training points in the computational domain may be randomly chosen. This allows us to train the NN more precisely around the critical regions of interest (e.g. the separatrix between closed and open field lines, the Y-point at the tip of the closed line region, the equatorial current sheet, etc.). It also allows us to easily change the computational domain when the separatrix changes shape as we will see below. Both of the above benefits are very hard to implement in classical finite-difference numerical integration techniques with adaptive mesh refinement in a moving and deformable numerical grid. In our 2D axisymmetric runs (see Section 4), we used 1000 random training points in the open field lines and 450 in the closed ones, 200 random training points along the boundaries r = r_* (the stellar surface), q = 0 (infinity), θ = π/2 (the equator), and 200 random training points along the separatrix. These random training points were updated every 5000 training steps. We trained our PINNs for 50 000 training steps before we were satisfied with the loss convergence. One particular encouraging result was the reproduction of important known features of the solution, such as the distribution of the poloidal electric current in the open line region (see Fig. 5).

The general 3D PINN numerical procedure described above should relax to a solution where the closed line region extends up to some distance inside the light cylinder. The main problem with all numerical methods is that the solution of the pulsar magnetosphere contains electric current sheets along the separatrix and along the equator where the force-free ideal MHD condition breaks down. Unfortunately, the pulsar equation is not informed of their presence, and therefore, without special attention from the part of the programmer, it is incapable of supporting the mathematical contact discontinuities along them. This problem is addressed in the following subsection. One other problem is that, in past numerical methods, the imposition of the constraint E < B was used to open up magnetic field lines that cross the light cylinder (E < B is automatically satisfied inside the light cylinder where rsin θ < R_LC). Therefore, if we want to study solutions where the closed line region does not extend up to the light cylinder, we must find other ways to open up the magnetosphere outside the closed line region. One way may be to impose the constraint that E < νB with ν < 1. In Section 3.3, we propose yet another way to open up the magnetosphere.

3.2 Decomposition of the computational domain

The first thing we want to do is to decompose the computational domain into the region of closed and open field lines. This is a central element of our methodology, namely that we choose from the very beginning of our simulation which field lines will be open and which ones will be closed. This is equivalent to an ad hoc determination from the very beginning of the extent of the polar cap (see below). In general, the solution that we will obtain will contain a closed line region that does not extend all the way to the light cylinder. It will, however, be a valid solution as we know since Contopoulos (2005) and Timokhin (2006). Up to this day it is not yet clear why in some time-dependent simulations the closed line region extends all the way to the light cylinder (e.g. Spitkovsky 2006; Kalapotharakos & Contopoulos 2009), while in others it extends to only about |$85{{\ \rm per\ cent}}$| of the light cylinder (e.g. Hu & Beloborodov 2022; Hakobyan et al. 2023). Which size of the polar cap corresponds to the physical (true) solution is a very important question that we want to elucidate in our next paper.

The two domains are separated by the separatrix surface characterized by the line where it originates on the surface of the star (see Fig. 1). In principle, one may obtain an analytic expression for the shape of the separatrix surface, namely

$$\begin{eqnarray} r_{\rm S}=r_{\rm S}(\theta ,\phi)\ . \end{eqnarray}$$

(17)

Obviously, the shape of the separatrix surface is unknown and must be determined. Let’s assume for the moment that we choose ad hoc some initial form of r_S in equation (17) that originates around the magnetic poles of the central star at inclined polar angles θ_m = θ_pc, 0 ≤ ϕ_m ≤ 2π and π − θ_m = θ_pc, 0 ≤ ϕ_m ≤ 2π. Here, θ_pc is the angular opening of the polar cap, and θ_m, ϕ_m are spherical coordinates centered around the inclined magnetic axis |$\hat{m}$| of the central star. Notice that for angles θ_m < θ_pc and θ_m > π − θ_pc, the radius connecting the point (r, θ, ϕ) with the origin does not cross the separatrix. Thus, in the open line region, the domain of integration becomes r ≥ r_*, 0 ≤ θ_m ≤ θ_pc and π − θ_pc ≤ θ_m ≤ π, and r ≥ r_S, θ_pc ≤ θ_m ≤ π − θ_pc. In the closed line region, the domain of integration becomes r_* ≤ r ≤ r_S, θ_pc ≤ θ_m ≤ π − θ_pc. The above r, θ_m, ϕ_m intervals correspond to r, θ, ϕ intervals and express mathematically what is obvious geometrically as closed line and open line regions in the schematic of Fig. 1.

Figure 1.

Schematic of the decomposition of the computational domain into a region of closed and open field lines. Closed line region: Between the central star r = r_* and the separatrix r = r_S(θ, ϕ) (thick red line). Open line region: Outside the separatrix (r > r_S(θ, ϕ)). In the open line region, all field lines are artificially outflowing and there is no equatorial current sheet during the PINN training. When the PINN is trained, the open field lines in the Southern hemisphere are reversed and the equatorial current sheet in the open line region appears along the surface of field-direction reversal (dotted red line).

Open in new tab Download slide

The reason we decided to separate the region of open and closed field lines is that, due to the electric current sheet that flows along the separatrix between the two regions, a strong contact discontinuity (of practically infinitesimal thickness) in the magnetic and electric fields exists along the separatrix. Uzdensky (2003) (see also Lyubarskii 1990) integrated equation (1) across current sheets and obtained that

$$\begin{eqnarray} \left(B^2-E^2\right)_{\rm below}=\left(B^2-E^2\right)_{\rm above}\ . \end{eqnarray}$$

(18)

Treating such discontinuities inside any computational domain is very problematic and any computation method would generate spurious Gibbs oscillations around the separatrix. This is seen in all previous solutions of the pulsar magnetosphere (e.g. Kalapotharakos et al. 2012b; Hu & Beloborodov 2022; Mahlmann et al. 2022). We propose here a better way to treat such discontinuities:

Solve equations (12)–(14) in the two regions independently for an initial arbitrary choice of the separatrix surface between them. An informed but not necessary unique initial choice would be the dipole magnetic field surface that originates at θ_m = θ_pc and θ_m = π − θ_pc on the stellar surface that corresponds to
$$\begin{eqnarray} \Psi _{\rm dipole}\equiv \Psi _{\rm max}\frac{r_*}{r}\sin ^2\theta _{\rm m}\equiv \Psi _{\rm max}\sin ^2\theta _{\rm pc}\ , \end{eqnarray}$$
(19)
which yields
$$\begin{eqnarray} r_{\rm S} \equiv r_* \frac{\sin ^2\theta _{\rm m}}{\sin ^2\theta _{\rm pc}}\ \mbox{if}\ \theta _{\rm pc}\le \theta _{\rm m}\le \pi -\theta _{\rm pc}. \end{eqnarray}$$
(20)
Obviously, for such an ad hoc shape of the separatrix surface, the quantity B² − E² will in general be discontinuous across that surface.
Iteratively adjust the shape of the separatrix surface r = r_S(θ, ϕ) at all points (θ, ϕ) so as to satisfy the condition that B² − E² is continuous at all points of the separatrix.
In other words, we choose an initial dipolar shape for the separatrix surface, train the two PINNs in the open and closed line regions for a number of steps, calculate B² − E² at each point of the separatrix in the adjacent closed and open line regions, and displace that point proportionally to the difference between the two corresponding values in the direction of the smaller value. We must acknowledge that there is a high probability that the shape of the separatrix r_S(θ, ϕ) may not be a single-valued function of θ and ϕ for highly inclined rotators. We will thus take particular care in 3D to displace each separatrix point along the direction perpendicular to its surface. After each re-adjustment of the separatrix surface continue the training of the PINN for a few more steps, and then repeat the re-adjustment. This procedure has never been tried before and, as we will see below, yields the final solution where B² − E² is continuous everywhere across the separatrix.

3.3 The numerical ‘disappearance’ of the equatorial current sheet

One second improvement to the solution method has to do with the numerical treatment of the equatorial current sheet that originates at the tip of the closed line region. The reason there exists an equatorial current sheet is that magnetic field lines leave one pole of the star, open up to infinity, and return from infinity to the other pole of the star. In doing so, they also carry the same amount of poloidal electric current in each hemisphere from each pole of the star to infinity. These two electric currents return to the star through the equatorial current sheet. In other words, the equatorial current sheet is there to close the global poloidal electric current circuit (CKF). It is thus obvious that, if we artificially (numerically) invert the direction of the field lines that leave the star from the southern pole, the electric current direction in the southern hemisphere will be inverted, and therefore there will be no need to close the global poloidal electric current circuit through an equatorial current sheet. This configuration is clearly artificial (it is equivalent to a magnetic monopole), but it is mathematically and dynamically equivalent to the configuration that we are investigating in the open line region, only without the mathematical discontinuity of the equatorial current sheet! We are able to implement this trick because we are treating the open line region independently from the closed line region. This is the same trick assumed by Bogovalov (1999) when he obtained the solution for the tilted split monopole. We here generalize his approach and show that it is also valid (and very helpful) in the numerical treatment of the open line region in the more general dipole magnetosphere. We have thus found a way to make the equatorial current sheet discontinuity ‘numerically disappear’.

A numerical implementation of this trick is to assume that, in the open line region outside the closed line region, the boundary conditions for the magnetic flux function on the stellar polar caps will be

$$\begin{eqnarray} \Psi _{\phi _{\rm m}}(r_*,\theta _{\rm m},\phi _{m})&=&\Psi _{\rm max}\sin ^2\theta _{\rm m}\ \mbox{for}\ 0\le \theta _{\rm m}\le \theta _{\rm pc}\ , \\ \Psi _{\phi _{\rm m}}(r_*,\theta _{\rm m},\phi _{m})&=&2\Psi _{\rm S}-\Psi _{\rm max}\sin ^2\theta _{\rm m} \\ &=&\Psi _{\rm max}(2\sin ^2\theta _{\rm pc}-\sin ^2\theta _{\rm m}) \\ &&\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \mbox{for}\ \pi -\theta _{\rm pc}\le \theta _{\rm m}\le \pi \\ \Psi _{r}(r_*,\theta _{\rm m},\phi _{\rm m})&=&0\ , \\ \Psi _{\theta _{\rm m}}(r_*,\theta _{\rm m},\phi _{\rm m})&=&0\ . \end{eqnarray}$$

(21)

What we have done here is to express the boundary conditions in terms of the vector magnetic flux components |$\Psi _r,\Psi _{\theta _{\rm m}},\Psi _{\phi _{\rm m}}$| in the inclined spherical system of coordinates (r, θ_m, ϕ_m) in which a pure magnetic dipole is described as |$\Psi _r=\Psi _{\theta _{\rm m}}=0$| and |$\Psi _{\phi _{\rm m}}=\Psi _{\rm max}\sin ^2\theta _m$|⁠. Obviously, in the open line region on the surface of the star, |$0\le \Psi _{\phi _{\rm m}}\le 2\Psi _{\rm S}$|⁠, whereas in the closed line region on the surface of the star, |$\Psi _{\rm S}\le \Psi _{\phi _{m}}\le \Psi _{\rm max}$| (Ψ_S ≡ Ψ_maxsin ²θ_pc). This mathematical trick simplifies tremendously the numerical treatment of the large-scale open line region. Notice that, with or without the artificial flux reversal in the southern hemisphere, the pressure balance condition along the separatrix remains the same. After the solution is obtained, the field will be reversed in the southern hemisphere, and the true magnetic field configuration will be obtained with an equatorial current sheet along the surface of flux direction reversal.

There is an important advantage of this method over standard numerical methods that include the equatorial current sheet in their domain of numerical integration. As θ approaches the equatorial current sheet at θ = θ_eq.c.s. from above, the first term in the expression for the field component B_r in equation (11), namely ∂Ψ_ϕ/∂θ/(r²sin θ), reaches in general a non-zero value that is immediately reversed below the equatorial current sheet. In other words,

$$\begin{eqnarray} \frac{\Psi _\phi }{\partial \theta }(r,\theta \rightarrow \theta _{\rm eq.c.s.}^-)=-\frac{\Psi _\phi }{\partial \theta }(r,\theta \rightarrow \theta _{\rm eq.c.s.}^+)\ne 0\ . \end{eqnarray}$$

(22)

Unfortunately, from the symmetry of Ψ_ϕ and its derivatives above and below the equatorial current sheet, in a numerical scheme without artificial field reversal in the southern hemisphere, the latitudinal derivative ∂Ψ_ϕ/∂θ will be equal to zero in the middle of the current sheet, which is not true in general. The numerical trick that we propose in this subsection naturally allows for ∂Ψ_ϕ/∂θ(r, θ → θ_eq.c.s.) ≠ 0, which, after the final field reversal, naturally satisfies equation (22).

Finally, by choosing the angular size θ_pc of the polar cap, this approach allows as to choose the amount of magnetic flux Ψ_S ≡ Ψ_maxsin ²θ_pc that will open to infinity.¹ Due to the artificial field reversal, this amount is ‘forced’ to extend to infinity, since it has no way to return to the star in the open line region. As we acknowledged above, we checked that the condition E < B is automatically satisfied in the open line region beyond the light cylinder.

4 FIRST RESULTS IN AXISYMMETRY

4.1 Mathematical formulation

Before we proceed to a solution of the full 3D problem in a forthcoming publication, we would like to consider first the simpler axisymmetric problem (∂/∂ϕ = 0) where a steady-state solution is known (CKF). In the case of the aligned rotator, it is customary to use the simplified notation Ψ_ϕ ≡ Ψ. In units where R_LC ≡ 1, equations (11) yield

$$\begin{eqnarray} B_r &\equiv & \frac{1}{r^2\sin \theta }\frac{\partial \Psi }{\partial \theta }\ , \\ B_\theta &\equiv & -\frac{1}{r\sin \theta }\frac{\partial \Psi }{\partial r}\ , \\ B_\phi &\equiv & \frac{I(\Psi)}{r\sin \theta }\ , \end{eqnarray}$$

(23)

and equation (14) becomes

$$\begin{eqnarray} && \left(1-r^2 \sin ^2\theta \right) \left[ \frac{\partial ^2 \Psi }{\partial r^2}-\frac{\partial \Psi }{\partial \theta }\frac{\cos \theta }{r^2 \sin \theta }+\frac{1}{r^2}\frac{\partial ^2 \Psi }{\partial \theta ^2}\right] \\ && \quad -2r\sin \theta \left[\frac{\partial \Psi }{\partial \theta }\frac{\cos \theta }{r}+\frac{\partial \Psi }{\partial r} \sin \theta \right]+II^{\prime }(\Psi)=0. \end{eqnarray}$$

(24)

Equation (24) is the well-known pulsar equation (Scharlemann & Wagoner 1973) in spherical coordinates. Here, I(Ψ) = rsin θ B_ϕ, and I′(Ψ) ≡ dI/dΨ = α. The condition that I = I(Ψ) (equation 9) may be rewritten as

$$\begin{eqnarray} \frac{\partial \Psi }{\partial r}\frac{\partial I}{\partial \theta }- \frac{\partial \Psi }{\partial \theta }\frac{\partial I}{\partial r}=0 \end{eqnarray}$$

(25)

I(Ψ) ≠ 0 only in the open line part of the magnetosphere. We consider here an untwisted magnetosphere with I(Ψ) = 0 in its closed line part. Our results may also be generalized in the case of a twisted magnetosphere (see Ntotsikas et al. 2024). In axisymmetry, the regularization condition at the light cylinder (equation 15) becomes

$$\begin{eqnarray} II^{\prime }&=&2(\cos \theta \ B_r-\sin \theta \ B_\theta)|_{\rm LC}\equiv 2B_z|_{\rm LC}\ . \end{eqnarray}$$

(26)

As before, equation (26) determines the value of the function II′(Ψ) along all field lines that cross the light cylinder and never return to the star. In an untwisted magnetosphere, all other field lines that do not cross the light cylinder and form the closed line region have II′(Ψ) = 0.

Equation (24) was first solved by Contopoulos et al. (1999) with a standard elliptic relaxation procedure based on a finite-difference grid in cylindrical coordinates. We will now solve the pulsar equation with machine learning.

4.2 Machine learning

In the aligned rotator, we need to obtain the solution Ψ(r, θ) and I(Ψ). It is helpful to reformulate the problem in new coordinates

$$\begin{eqnarray} q &\equiv & \frac{1}{r}, \\ \mu &\equiv & \cos \theta. \end{eqnarray}$$

In those new variables, equation (24) may be rewritten as

$$\begin{eqnarray} && \left(q^2-\left(1-\mu ^2\right)\right)\left(q^2\frac{\partial ^2\Psi }{\partial q^2} +\left(1-\mu ^2\right)\frac{\partial ^2\Psi }{\partial \mu ^2}\right) \\ &&+2q^3\frac{\partial \Psi }{\partial q}+2\left(1-\mu ^2\right)\mu \frac{\partial \Psi }{\partial \mu } +II^{\prime }=0 . \end{eqnarray}$$

(27)

In our rescaled radial variable q, the outer radial boundary r → ∞ lies at q = 0. We tried to impose various boundary conditions at q = 0 (e.g. B_θ = 0, or Ψ = 0, or Ψ = Ψ_S = const.), and in all cases the resulting poloidal field configuration became asymptotically radial. From this we concluded that, as long as the outer radial boundary of our calculation lies at infinity, the inner magnetospheric solution does not depend on the particular boundary conditions that we impose there. In order to help the solution satisfy the boundary conditions along the axis and the stellar surface, namely Ψ(r, θ = 0, π) = I(r, θ = 0, π) = 0 and Ψ(r_*, θ) = sin ²θ Ψ_max, we define new functions |$f,{\cal I}$| as the new exits of the NN such that

$$\begin{eqnarray} \Psi &\equiv & \left(1-\mu ^2\right)f, \\ I &\equiv & \left(1-\mu ^2\right){\cal I}. \end{eqnarray}$$

(28)

In these new functions, the pulsar equation (equation 27) may be rewritten as

$$\begin{eqnarray} && \left(q^2-\left(1-\mu ^2\right)\right)\left(q^2\frac{\partial ^2 f}{\partial ^2 q}+\left(1-\mu ^2\right) \frac{\partial ^2 f}{\partial \mu ^2}\right) \\ &&+2q^3 \frac{\partial f}{\partial q} - (4\mu ^2+2\left(q^2-\left(1-\mu ^2\right)\right)f \\ &&+(2\mu \left(1-\mu ^2\right)-4\mu \left(q^2-\left(1-\mu ^2\right)\right)\frac{\partial f}{\partial \mu }+{\cal S} = 0\ , \end{eqnarray}$$

(29)

while the condition that I = I(Ψ) (equation 25) may be rewritten as

$$\begin{eqnarray} \frac{\partial f}{\partial q}\left(-2\mu {\cal I}+\left(1-\mu ^2\right)\frac{\partial {\cal I}}{\partial \mu }\right)- \frac{\partial {\cal I}}{\partial q}\left(-2\mu f+\left(1-\mu ^2\right) \frac{\partial f}{\partial \mu }\right)=0. \\ \end{eqnarray}$$

(30)

The source term in equation (29) is equal to

$$\begin{eqnarray} {\cal S} & \equiv & {\cal I}\frac{{\rm d}I}{{\rm d}\Psi }= {\cal I} \frac{|\nabla I|}{|\nabla \Psi |}\cdot {\rm sgn}\left(\nabla I \cdot \nabla \Psi \right) \\ &=& {\cal I}\frac{\sqrt{q^2\left(1-\mu ^2\right)\left(\frac{\partial {\cal I}}{\partial q}\right)^2+\left(2\mu {\cal I}-\left(1-\mu ^2\right)\frac{\partial {\cal I}}{\partial \mu }\right)^2}}{\sqrt{q^2\left(1-\mu ^2\right)\left(\frac{\partial f}{\partial q}\right)^2+\left(2\mu f-\left(1-\mu ^2\right)\frac{\partial f}{\partial \mu }\right)^2}} \\ && \cdot {\rm sgn}\left(q^2\left(1-\mu ^2\right)\frac{\partial {\cal I}}{\partial q} \frac{\partial f}{\partial q}\right. \\ && +\left.\left(2\mu {\cal I}-\left(1-\mu ^2\right)\frac{\partial {\cal I}}{\partial \mu }\right)\left(2\mu f-\left(1-\mu ^2\right)\frac{\partial f}{\partial \mu }\right)\right). \end{eqnarray}$$

(31)

In this reformulation of the problem, the entries of the PINN are the spatial coordinates (q, μ), and the exits are |$(f,{\cal I})$|⁠. The loss functions must then satisfy the conditions

Equation (29)
Equation (30)
f = 1 on the stellar surface
|${\cal I}=0$| in the closed line region.

Notice that there is no special provision along the light cylinder as expressed by equation (15), since if equation (29) is satisfied, condition (15) is also automatically satisfied. This is a welcome advantage of the PINN methodology.

4.3 Domain decomposition

The domain decomposition is simpler in the case of the aligned pulsar where the direction |$\hat{m}$| of the stellar magnetic dipole is along θ = 0, and therefore θ_m ≡ θ and ϕ_m ≡ ϕ (see Fig. 2). The open line region is r ≥ r_S, 0 ≤ θ ≤ π. The closed line region is r_* ≤ r ≤ r_S, θ_pc ≤ θ ≤ π − θ_pc. Here, the initial choice of the separatrix radius r_S is

$$\begin{eqnarray} r_{\rm S}(\theta) = r_*\frac{\sin ^2\theta }{\sin ^2\theta _{\rm pc}}\ \ \ \ \ \mbox{if}\ \theta _{\rm pc}\lt \theta \lt \pi -\theta _{\rm pc}. \end{eqnarray}$$

(32)

The shape of r_S(θ) will be adjusted iteratively in order to attain continuity of |$B^2-E^2\equiv B_p^2(1-(r\sin \theta)^2)+B_\phi ^2$|⁠, or equivalently continuity of

$$\begin{eqnarray} &&\left(q^2-\left(1-\mu ^2\right)\right)\left(q^2 \left(1-\mu ^2\right)\left(\frac{\partial f}{\partial q}\right)^2 + 4\mu ^2 f^2 \right. \\ &&+\left. \left(1-\mu ^2\right)^2 \left(\frac{\partial f}{\partial \mu }\right)^2 -4\mu \left(1-\mu ^2\right)f \frac{\partial f}{\partial \mu }\right) +\left(1-\mu ^2\right){\cal I}^2 \end{eqnarray}$$

(33)

across the separatrix as described in Section 4.1. In the latter expression, |${\cal I}=0$| in the closed line region, and |${\cal I}\ne 0$| in the open line region across the separatrix.

Figure 2.

Same as Fig. 1 in axisymmetry. Note that we have implemented artificial field-direction reversal in the open line region.

Open in new tab Download slide

We want to impose that Ψ(r_S, θ) = Ψ_S = Ψ_max sin ²θ_pc in the separatrix between the two domains. We, therefore, introduce one more loss function in both domains, namely one that requires

$$\begin{eqnarray} f\left(q=\frac{1}{r_{\rm S}},\mu \right)=\Psi _{\rm max}\frac{\sin ^2\theta _{\rm pc}}{1-\mu ^2}\ . \end{eqnarray}$$

(34)

4.4 Field reversal

The boundary conditions for the magnetic field in the closed line region r_* ≤ r ≤ r_S(θ), θ_pc < θ < π − θ_pc are

$$\begin{eqnarray} \Psi (r_*,\theta)&=&\Psi _{\rm max}\sin ^2\theta \ , \\ \Psi (r_{\rm S}(\theta),\theta)&=&\Psi _{\rm S}\ . \end{eqnarray}$$

(35)

As before, θ_pc is the opening of the polar cap, and Ψ_S ≡ Ψ_max sin ²θ_pc. We implement the mathematical field reversal in the open line region outside the closed line region with boundary conditions

$$\begin{eqnarray} \Psi (r_*,\theta)&=&\Psi _{\rm max}\sin ^2\theta \ \ \ \ \ \ \ \ \ \ \ \mbox{for}\ 0\le \theta \le \theta _{\rm pc}\ , \\ \Psi (r_*,\theta)&=&2\Psi _{\rm S}-\Psi _{\rm max}\sin ^2\theta \\ &=&\Psi _{\rm max}(2\sin ^2\theta _{\rm pc}-\sin ^2\theta)\ \ \mbox{for}\ \pi -\theta _{\rm pc}\le \theta \le \pi \ \ \\ \Psi (r_{\rm S}(\theta),\theta)&=&\Psi _{\rm S}\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \mbox{for}\ \theta _{\rm pc}\lt \theta \lt \pi -\theta _{\rm pc}\ . \end{eqnarray}$$

(36)

Obviously, in the open line region, 0 ≤ Ψ ≤ 2Ψ_S, whereas in the closed line region, Ψ_S ≤ Ψ ≤ Ψ_max. This mathematical trick simplifies tremendously the numerical treatment of the large-scale open line region. There is no problem with the closed line region because the two regions are in contact along Ψ(r, θ) = Ψ_S. After the solution is obtained, the field will be reversed in the open line region where Ψ_S ≤ Ψ ≤ 2Ψ_S, and the true magnetic field configuration will be presented with the equatorial current sheet along Ψ = Ψ_S.

As we discussed above, one extra benefit of this approach is that B_r(r > 1, θ = π/2⁻) is naturally allowed to be non-zero and equal to −B_r(r > 1, θ = π/2⁺) across the equatorial current sheet. However, without the mathematical field reversal that we propose in the open line region, due to symmetry, B_r = (∂Ψ/∂θ)/(r²sin θ) = 0 along the equator θ = π/2. One might think that in a typical MHD/FFE grid simulation, the equatorial field reversal takes place within one grid cell. Unfortunately, in all known solutions in the literature the equatorial current sheet extends over several grid cells in thickness, which is not correct. This unfortunate result has the following repercussion: all components of the magnetic field are zero in a region of finite thickness immediately outside the Y-point. As a result, the magnetic (and electric) field pressure from immediately inside the tip of the closed line region pushes and protrudes outwards, forming a true Y-point as seen in all MHD/FFE numerical simulations. As is shown in detail in Contopoulos et al. (2024), however, this is not true, and the Y-point is in fact a T-point as first proposed by Uzdensky (2003). This is captured in our present PINN solution with a proper treatment of the separatrix and equatorial current sheet via domain decomposition and field direction reversal in the open line region (see below).

4.5 First results

We implemented three independent NNs with multiple internal (hidden) layers:

One PINN with two entries q, μ, three hidden layers with 64 nodes each, and two exits |$f,{\cal I}$|⁠. This PINN solves the pulsar equation in the open line region outside the separatrix and is the most complex among the three. It encounters no problems at the light cylinder where the functional form I|_LC = I(Ψ|_LC) is determined. Its training requires tens of thousands of steps, mainly because the condition that I = I(Ψ) everywhere is very slowly enforced. The types of the individual activation functions and the weights of the individual losses are carefully chosen so as to yield consistent and stable convergence of the training.
One PINN with two entries q, μ, three hidden layers with 64 nodes each, and one exit f. This PINN solves the pulsar equation in the closed line region inside the separatrix where |${\cal I}=0$|⁠. The training of this PINN is stable and fast.
The most crucial part of our method is the readjustment of the shape of the separatrix. We choose random angles θ_pc < θ < π − θ_pc where we know the values of B²(r_S, θ) − E²(r_S, θ) both inside the outside the separatrix from the two PINNS trained in the closed and open line regions, respectively. In general, these two values are not equal, so the radius of the separatrix must be adjusted at those angles. We implemented the adjustment
$$\begin{eqnarray} r_{\rm S\ new} = r_{\rm S}+2\beta \ \left(\frac{r_{\rm S}-r_*}{1-r_*}\right)^2\ \frac{(B^2-E^2)_{r=r_{\rm S}^-} - (B^2-E^2)_{r=r_{\rm S}^+}}{(B^2-E^2)_{r=r_{\rm S}^-} + (B^2-E^2)_{r=r_{\rm S}^+}} \\ \end{eqnarray}$$
(37)
at those angles θ. In the above expression, the separatrix is displaced due to the pressure imbalance in the third term of the right hand side product. The first term is introduced to avoid moving the separatrix footpoint on the stellar surface, and β is an adjustable positive parameter. After the new position of the separatrix is defined on these randomly chosen angles θ_pc < θ < π − θ_pc, a third NN with one entry μ, two hidden layers with 128 nodes each, and one exit r_S is trained to yield the shape of the separatrix at all angles θ required by the first two main PINNs. We have assumed here that the separatrix shape r_S(θ) is a single-valued function of θ. This is indeed true in axisymmetry, but is not in general true in highly inclined oblique rotators. In that case, we will displace each point of the separatrix in the direction perpendicular to its surface, not in the radial direction as in equation (37).

The number of internal nodes and layers in the two PINNs that solve the pulsar equation was chosen by trial and error to be able to reproduce known features of the solution (e.g. the functional form of the electric current distribution I = I(Ψ) shown in Fig. 5). Our number of internal nodes per layer is comparable to that in the PINN implementation by Stefanou et al. (2023b). The number of internal nodes in the third NN was chosen so that it describes in detail the deformed shape of the separatrix.

In Fig. 3, we show the evolution of our training losses as a function of the number of training steps. The various losses that we implemented to be minimized to zero are the averages of (f − 1)² over the stellar surface, of (Ψ − Ψ_S)² along the separatrix and the equator in the open line region, and of course the average squares of the residuals of the central PDEs of the problem (equations 29 and 30). The total loss is the sum of the above losses, each one weighted by a corresponding adjustable weight factor (see discussion in Section 5). As we can see from Fig. 3, PDE errors (i.e. the residuals of equation 29) at the end of the training in the open line region are around 10⁻², which, for magnetic flux values on the order of Ψ ∼ Ψ_max = 1 and (r_*/R_LC)Ψ_max = 0.25, respectively, correspond to a conservative accuracy of about two significant digits in the value of Ψ. This is better than CKF, but worse than Timokhin (2006) who mentions relative residuals on the order of 10⁻⁴, and can certainly be improved in the future with further development of our method.

$Evolution of the various losses (average square residuals) with the number of training steps for the two PINNs in the closed and open line regions. Initial PINN training: 50 000 steps. Subsequent separatrix shape adjustments and PINN retraining: 10 × 50 000 steps. Line notation for the closed/ open line regions, respectively: ‘pde in’/‘pde out’: equation (29), the pulsar equation (note that ${\cal I}=0$ in the closed line region); ‘sep in’/‘sep out’: separatrix boundary condition; ‘paral’: equation (30), I = I(Ψ). Abrupt loss changes correspond to updates in the training set.$

Figure 3.

Evolution of the various losses (average square residuals) with the number of training steps for the two PINNs in the closed and open line regions. Initial PINN training: 50 000 steps. Subsequent separatrix shape adjustments and PINN retraining: 10 × 50 000 steps. Line notation for the closed/ open line regions, respectively: ‘pde in’/‘pde out’: equation (29), the pulsar equation (note that |${\cal I}=0$| in the closed line region); ‘sep in’/‘sep out’: separatrix boundary condition; ‘paral’: equation (30), I = I(Ψ). Abrupt loss changes correspond to updates in the training set.

Open in new tab Download slide

In practice, the two main PINNs are initially trained for 50 000 steps after which, satisfactory convergence is achieved in the two regions (all residual losses fall below 10⁻²). After that initial training stage, the separatrix is displaced based on the resulting pressure differences between the closed and open line regions, and then the third NN is trained. The two main PINNs are re-trained for another 50 000 steps, and the process is repeated 10 times. In total, we run 500 000 training steps for the two main PINNs. In the resulting configuration, pressure balance is achieved across the separatrix to less than 1 per cent. We further tested the accuracy of our PINN method in the Appendix by applying it to reproduce the standard vacuum magnetostatic dipole whose solution we know analytically (Ψ(r, θ) = Ψ_max sin ²θ r_*/r; Figs A1 and A2). We can see there that the agreement with the analytic solution is acceptable, at least in the region of interest inside the light cylinder where the density of training points is higher.

We present our preliminary results in Figs 4 and 5 for θ_pc = 1.176(r_*/R_LC)^1/2 and Ψ_open ≡ Ψ_S = Ψ_maxsin ²θ_pc for r_* = 0.25R_LC. This particular value of θ_pc was specially chosen with the following in mind: It is straightforward to calculate that, in a dipolar magnetic field configuration, the magnetic field line that crosses the light cylinder corresponds to Ψ_dipole LC = Ψ_maxr_*/R_LC, and θ_dipole pc = arcsin[(r_*/R_LC)^1/2] ≈ (r_*/R_LC)^1/2. In previous high-resolution solutions of the pulsar equation (e.g. Timokhin 2006), the magnetic field line that reaches the light cylinder corresponds to Ψ_{LC Timokhin} = 1.23Ψ_dipole LC and θ_{pc Timokhin} = arcsin[(1.23r_*/R_LC)^1/2] = 1.176(r_*/R_LC)^1/2 for r_* = 0.25R_LC. Therefore, this particular value of θ_pc (and correspondingly Ψ_open) was chosen because in the standard solution (e.g. Timokhin 2006 ) the Y-point lies very close to the light cylinder.

Figure 4.

Representative axisymmetric solution with domain decomposition and field reversal in the southern hemisphere. Solution for θ_pc = 1.176(r_*/R_LC)^1/2, or equivalently, Ψ_open ≡ Ψ_S = 1.23Ψ_dipole LC. Thin lines: Poloidal magnetic field lines (Ψ iso-contours along integer multiples of 0.05Ψ_dipole LC, with Ψ = 0 along the z-axis). Thin red dotted line: Initial dipole position of separatrix surface. Thick lines: Separatrix Ψ_open ≡ Ψ_S = Ψ_maxsin ²θ_pc and equatorial current sheet. Vertical line: Light cylinder. Notice the shape of the tip of the closed line region. In agreement with Uzdensky (2003), its true shape is a T-point.

Open in new tab Download slide

Figure 5.

Top: Electric current distribution I(Ψ)R_LC/Ψ_S across open magnetic field lines 0 ≤ Ψ/Ψ_S ≤ 1 for the solution of Fig. 4. Dotted line: Analytic expression for split monopole. Our result is in excellent agreement with fig. 3 of Timokhin (2006). We consider this an important success of the PINN methodology. Middle: the distribution of B_z(r, θ = π/2) inside (black line) and outside (red line) the Y-point. The dotted line represents the standard 1/r³ dipolar dependence. We clearly infer both the outward stretching of the initially dipolar field (at small r distances the solid line lies below the dotted one), and the inward squeezing of the field right behind the Y-point that is required to satisfy pressure balance (in agreement with fig. 11 of Timokhin 2006). Bottom: The distribution of B_r(r, θ)/B_r(r, 0) for r = R_Y/R_LC/2R_LC/3R_LC (black/green/red/magenta lines, respectively). We see the gradual transition from a dipole field (B_r(r, θ) = B_r(r, 0)cos θ) to a split-monopole field (B_r(r, θ) = B_r(r, 0)sgn(cos θ)) with distance from the star. Notice that B_r(r, θ = π/2⁻) = −B_r(r, θ = π/2⁺) ≠ 0 (without field reversal, B_r(r, θ = π/2) would erroneously be equal to zero).

Open in new tab Download slide

We found that R_Y = 0.88R_LC which is outside R_dipole(Ψ_S) = r_*/sin ²θ_pc = 0.81R_LC, as expected from previous solutions, but also inside R_LC. We suspect that this is due to our treatment of the separatrix which differs from that in all previous MHD/FFE/PIC solutions. In the present work, the separatrix is treated as a perfect (zero thickness) mathematical discontinuity, and is found to lie at a position where pressure balance is satisfied as expressed by equation (18). In all previous solutions, however, the separatrix return current is distributed among a bundle of magnetic field lines that contains a finite amount of poloidal magnetic flux δΨ ≠ 0 (e.g. δΨ = 0.01−0.06Ψ_S in Timokhin 2006). The outer part of this bundle represents the last closed field line that extends very close to the light cylinder and reaches the equator at a non-vertical angle. This position is what is usually characterized as the Y-point. The separatrix return current, however, is distributed till the inner part of this bundle which represents the last closed field line without poloidal electric current. This reaches the equator vertically inside the light cylinder (see fig. 4 in Timokhin 2006, figs 2 and 3 in Contopoulos et al. 2024, for details; see also fig. 1 in Hakobyan et al. 2023). In other words, in all previous numerical solutions there is no well-defined Y-point, only an extended region of distributed separatrix return current. We do not know whether the separatrix return current is indeed distributed. This merits further investigation which is, however, outside the scope of the present work. Under our present assumption that the separatrix return current is not distributed over Ψ, R_Y(Ψ_S) is found to lie slightly inside the light cylinder.

Another interesting point is the squeezing of the separatrix surface in the latitudinal direction right above the stellar surface. This is due to the pressure of the nonzero toroidal field |$B_\phi ^2$| in the open line region. Furthermore, as expected, the Y-point is indeed a T-point, i.e. the separatrix crosses the equator vertically. In fact, the closed line region is squeezed inwards at its tip on the equator, which results in the slight increase of |$B_z(r=R^-_{\rm Y},\theta =\pi /2)$| that is required to satisfy pressure balance accross the T-point (the superscript − is used to denote ‘right inside the Y-point’; Fig. 5 middle; see also fig. 11 in Timokhin 2006). Moreover, B_r(r > R_Y, θ = π/2⁻) = −B_r(r > R_Y, θ = π/2⁺) ≠ 0 (superscripts −/ + are used to denote ‘right above/below the equator’, respectively; Fig. 5 bottom). The above observations merit further investigation in a future publication.

5 DISCUSSION AND CONCLUSIONS

In this work, we have presented a novel numerical method for solving the generalized pulsar equation with PINNs. We have observed that PINNs behave satisfactorily around the light cylinder, but fail to properly treat the separatrix current sheet. For this reason, we have also introduced two innovations, namely the decomposition of the computation domain into two regions (one inside the separatrix current sheet and one outside), and the mathematical field reversal in the open line region. The latter allows us to completely avoid the equatorial current sheet discontinuity outside the light cylinder. The former allows PINNs to relax to a continuous and smooth solution that does not contain any mathematical singularities. The domain decomposition and field reversal that we introduced may also be implemented in more general situations with MHD current sheets, as for example in active regions of the solar corona. Regarding computational cost, axisymmetric PINN training is completed in a few hours on one GPU, whereas standard high-resolution solutions of the pulsar equation require about 1 d on a standard PC. Standard HPC time-dependent 3D simulations require days, and we can only hope that extrapolation of our method to an oblique rotator will be faster. Although the PINN methodology has better behaviour around the light cylinder and slightly faster convergence, we do not see any significant advantage over more traditional methods. One could also use standard numerical techniques to treat the deformable regions of open and closed field lines by constructing curvilinear grids that adapt to the evolving geometry of the boundary between the two regions.

While NNs possess considerable power, it’s important to recognize their inherent challenges and limitations. Without careful supervision, there’s no assurance that they will converge to the accurate solution, even after hundreds of thousands of iterations. Quoting from Moschou et al. (2023), ‘PINNs suffer from competing losses during gradient descent that can lead to poor performance especially in physical setups involving multiple scales’. Our experience has shown that NNs can easily become ‘trapped’ in a configuration with substantial losses, in particular when the set of training points is not balanced, i.e. if one loss is trained over many more training points than others. As such, users of NNs must be vigilant, particularly in increasing the weights of those losses that appear to be inactive during their training. It is crucial that the training process is never left unmonitored. Users should track the evolution of the losses and may need to re-calibrate their weights. With sufficient experience working with NNs, users can develop an understanding of when the training process is failing and when a loss weight readjustment is necessary. This hands-on approach ensures the effective application of NNs and mitigates potential issues.

We have obtained a preliminary solution of the pulsar equation for a particular value of the open magnetic flux Ψ_open = Ψ_maxsin ²θ_pc with θ_pc = 1.176(r_*/R_LC)^1/2, in which the closed line region is found to extend radially up to 0.88R_LC. We would like to correct here an important misunderstanding found in the literature. With the solution obtained in Section 4, we confirm Uzdensky (2003) that in the presence of the separatrix return current sheet, the Y-point is in fact a T-point at its tip. This means that at its tip it has a finite height and it is not at a non-vertical angle with respect to the equator as Gruzinov (2005) surmises. This realization makes a huge difference in the calculation of the electromagnetic field energy contained in this region of finite height at the tip of the closed line region which diverges as the Y-point approaches the light cylinder (see Contopoulos et al. 2024, for a detailed analysis of this effect). This result differs from Gruzinov (2005) that the magnetospheric energy contained around the Y-point is finite (that result was based on his conclusion that the separatrix arrives on the equator at a non-vertical angle; in that case, the electromagnetic energy contained in the tip of the closed line region indeed remains finite). This could be another reason why the Y-point is located well inside the light cylinder in all high-resolution global PIC simulations conducted in the last decade.

In summary, our implementation differs from Stefanou et al. (2023b) who were the first to solve the pulsar equation with a PINN and their work was an inspiration for our present work. We were the first to introduce a proper treatment of mathematical contact discontinuities in FFE. This yielded important details that differ from previous solutions of the pulsar magnetosphere (e.g. the shape and the extent of the tip of the closed line region). A more detailed description of our results in the axisymmetric magnetosphere will be presented in a follow-up publication. A detailed solution of the 3D magnetosphere where the real power of our method will be manifested will also be presented in the near future.

ACKNOWLEDGEMENTS

We would like to thank Petros Stefanou and Jose Pons for interesting and inspiring discussions. We also thank the anonymous referee for insisting that we should improve the accuracy of our preliminary results. This led us to acknowledge the challenges and limitations of NNs. We would also like to thank the International Space Science Institute (ISSI) for providing financial support for the organization of the meeting of ISSI Team No. 459 led by I. Contopoulos and D. Kazanas where the ideas presented in the paper originated. This research work was supported by the Hellenic Foundation for Research and Innovation (HFRI) under the 4th Call for HFRI PhD Fellowships (Fellowship Number: 9239).

DATA AVAILABILITY

The data underlying this article will be shared on reasonable request to the corresponding author.

Footnotes

Notice that we can even choose a more general non-circular polar cap by specifying a non-uniform polar cap angular distribution θ_pc(ϕ_m), as e.g. in Pétri (2018). This will be implemented when we will adjust the azimuthal shape of the polar cap so that the separatrix surface touches the light cylinder at all azimuthal angles ϕ. Notice also that, because we assume no twisting in the closed line region, the angular sizes of the north and south polar caps are the same in azimuth, i.e. the shape of the south polar cap is equal to θ_pc south(ϕ_m) = π − θ_pc north(ϕ_m). This effect is caused because field lines in the corotating closed line region start and end along the same star-centred meridional angle ϕ_m and at the same latitudinal distances from each magnetic pole. This is also true for the last closed field lines immediately interior to the separatrix surface.

References

Blandford

2002

, in

Gilfanov

Sunyaev

Churazov

, eds,

ESO Astroph. Symp

Springer-Verlag

, Berlin, p.

381

Google Scholar

Google Preview

OpenURL Placeholder Text

WorldCat

Bogovalov

1999

A&A

349

1017

Cerutti

Philippov

A. A.

Dubus

2020

A&A

642

204

Month:	Total Views:
January 2024	8
February 2024	93
March 2024	41
April 2024	49
May 2024	45
June 2024	31
July 2024	25
August 2024	16
September 2024	27
October 2024	41
November 2024	17
December 2024	14
January 2025	16
February 2025	20
March 2025	26
April 2025	57
May 2025	2

Article Contents

The pulsar magnetosphere with machine learning: methodology

ABSTRACT

1 INTRODUCTION

2 GENERAL MATHEMATICAL FORMULATION

3 NOVEL NUMERICAL APPROACH

3.1 Machine learning

3.2 Decomposition of the computational domain

3.3 The numerical ‘disappearance’ of the equatorial current sheet

4 FIRST RESULTS IN AXISYMMETRY

4.1 Mathematical formulation

4.2 Machine learning

4.3 Domain decomposition

4.4 Field reversal

4.5 First results

5 DISCUSSION AND CONCLUSIONS

ACKNOWLEDGEMENTS

DATA AVAILABILITY

Footnotes

References

APPENDIX A

Citations

Views

Altmetric

Email alerts

Astrophysics Data System

Citing articles via

Latest

Most Read

Most Cited

This Feature Is Available To Subscribers Only