iint

# The Pulsar Magnetosphere with Machine Learning: Methodology

Ioannis Contopoulos1, Ioannis Dimitropoulos2, Vassilis Mpisketzis3, Evangelos
Chaniadakis4
1 Research Center for Astronomy and Applied Mathematics, Academy of Athens,
Athens 11527, Greece
2 Department of Physics, University of Patras, Rio 26504, Greece
3 Department of Physics, National and Kapodistrian University of Athens,
Athens 15784, Greece
4 Department of Informatics and Telecommunications, National and Kapodistrian
University of Athens, Athens 15784, Greece E-mail<EMAIL_ADDRESS>
###### Abstract

We propose a new method for obtaining the general solution of the ideal force-
free steady-state pulsar magnetosphere in 3D. We divide the magnetosphere in
the regions of closed and open field lines and train two custom Physics
Informed Neural Networks (PINNs) to yield the solution in each of these two
regions. We also periodically adjust the shape of the separatrix between the
two regions to satisfy pressure balance everywhere. Our method introduces
several innovations over traditional methods that are based on numerical grids
and finite differences. In particular, it introduces a proper treatment of
mathematical contact discontinuities in FFE. We present preliminary results in
axisymmetry which confirm the significant potential of our method.

###### keywords:

pulsars – magnetic fields – numerical methods – machine learning

##  1 Introduction

55 years 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 (Magneto-Hydro-Dynamics) 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 Cerutti et al. (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 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:

  1. 1.

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% 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 2023). We still need to
understand where the Y-point truly lies, and how fast it moves out as the
pulsar spins down.

  2. 2.

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^{2}-E^{2}$ is everywhere continuous (see below). MHD/FFE methods smooth out
such discontinuities accross several computational grid cells. Unfortunately,
MHD/FFE conditions are not valid inside current sheets. Moreoever, 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.
2012). 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).

  3. 3.

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_{\rm LC}$, namely $\propto
1-\beta_{\rm rec}\ln(r/R_{\rm LC})$. With $\beta_{\rm rec}\sim 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, Philippov & Spitkovsky (2023), this logarithmic gradual
electromagnetic energy dissipation applies to the full magnetosphere, not only
to its undulating non-axisymmetric component as Cerutti, Philippov & Dubus
(2020) claim. This leaves no pulsar wind to extend to the termination shock to
power the X-ray pulsar nebula. We do not understand this result.

  4. 4.

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
spotaneously jump between them (e.g. Ntotsikas, Gourgouliatos & Contopoulos
2023).

  5. 5.

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, Philippov &
Spitkovsky 2023; gamma-ray fundamental plane as in Kalapotharakos, Wadiasingh,
Harding & Kazanas 2022; inverse Compton as in Richards & Lyutikov 2018; etc.).
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 the present 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 the next section. In § 4 we
propose two novel numerical approaches that will allow us to obtain the
solution in 3D. In § 5 we formulate the problem in the axisymmetric case of
the aligned rotater and present our first solutions of the pulsar equation
with a trained Physics Informed Neural Network. We conclude in § 6 with a
discussion of our results and forthcoming investigations.

##  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

$\displaystyle\rho_{\rm e}{\bf E}+{\bf J}\times{\bf B}=0\ ,$ (1)

where ${\bf B}$ and ${\bf E}$ are the magnetic and the electric field
respectively, $\rho_{\rm e}$ is the electric charge density, and ${\bf J}$ is
the electric current

$\displaystyle{\bf J}$ $\displaystyle\equiv$
$\displaystyle\frac{c}{4\pi}\nabla\cdot{\bf E}\ \frac{{\bf E}\times{\bf
B}}{B^{2}}+\frac{c}{4\pi}\frac{{\bf B}\cdot\nabla\times{\bf B}-{\bf
E}\cdot\nabla\times{\bf E}}{B^{2}}\ {\bf B}\ .$ (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

$\displaystyle\nabla\cdot{\bf B}=0\ ,\ \nabla\times{\bf
B}=\frac{1}{c}\frac{\partial{\bf E}}{\partial t}+\frac{4\pi}{c}{\bf J}\ ,$
$\displaystyle\nabla\cdot{\bf E}=\frac{4\pi}{c}\rho_{\rm e}\ ,\
\nabla\times{\bf E}=-\frac{1}{c}\frac{\partial{\bf B}}{\partial t}\ .$ (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
$\Omega$ around the axis $\theta=0$, namely

$\displaystyle\frac{\partial{\bf B}}{\partial t}$ $\displaystyle=$
$\displaystyle\left.\frac{\partial{\bf B}}{\partial t}\right|_{\rm
corot}-\nabla\times(r\sin\theta\ \Omega\ \hat{\phi}\times{\bf B})\ ,$
$\displaystyle\frac{\partial{\bf E}}{\partial t}$ $\displaystyle=$
$\displaystyle\left.\frac{\partial{\bf E}}{\partial t}\right|_{\rm
corot}-\nabla\times(r\sin\theta\ \Omega\ \hat{\phi}\times{\bf E})+r\sin\theta\
\Omega\ \hat{\phi}\ \nabla\cdot{\bf E}\ ,$
$\displaystyle\frac{\partial\rho_{\rm e}}{\partial t}$ $\displaystyle=$
$\displaystyle\left.\frac{\partial\rho_{\rm e}}{\partial t}\right|_{\rm
corot}+r\sin\theta\ \Omega\ \hat{\phi}\cdot\nabla\rho_{\rm e}\ .$ (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

$\displaystyle\nabla\times\left({\bf E}+\frac{r\sin\theta\
\Omega}{c}\hat{\phi}\times{\bf B}\right)$ $\displaystyle=$ $\displaystyle 0\
,$ (5) $\displaystyle\nabla\times\left({\bf B}-\frac{r\sin\theta\
\Omega}{c}\hat{\phi}\times{\bf E}\right)$ $\displaystyle=$
$\displaystyle\frac{4\pi}{c}{\bf J}-\frac{r\sin\theta\ \Omega}{c}\hat{\phi}\
\nabla\cdot{\bf E}\ .$ (6)

Let us define here the light cylinder radius $R_{\rm LC}\equiv c/\Omega$. From
eq. (5) above we obtain that

$\displaystyle{\bf E}\equiv{\bf E}_{p}=-\frac{r\sin\theta}{R_{\rm
LC}}\hat{\phi}\times{\bf B}=\frac{r\sin\theta}{R_{\rm
LC}}(B_{\theta},-B_{r},0)$ (7)

in component form, where the subscript ‘p’ denotes a poloidal component.
Following eq. (7), ${\bf E}\cdot\nabla\times{\bf E}=0$, and the expression for
${\bf J}$ in eq. (2) simplifies considerably. Substituting everything back in
eq. (6) we obtain $\nabla\times\\{{\bf B}_{p}(1-(r\sin\theta/R_{\rm
LC})^{2})+{\bf B}_{\phi}\\}\times{\bf B}=0$ or equivalently,

$\displaystyle\nabla\times\left\\{{\bf
B}_{p}\left(1-\left(\frac{r\sin\theta}{R_{\rm LC}}\right)^{2}\right)+{\bf
B}_{\phi}\right\\}=\alpha{\bf B}\ .$ (8)

Here, $\alpha$ is the force-free parameter that obeys the constraint

$\displaystyle{\bf B}\cdot\nabla\alpha=0\ ,$ (9)

namely that $\alpha$ is constant along magnetic field lines. In other words,
field lines lie along isocontours of $\alpha$ 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 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 more naturally reaches
the final steady state. Eq. (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 ${\bf A}$ such that ${\bf
B}=\nabla\times{\bf A}$. We will work in spherical coordinates
$(r,\theta,\phi)$ centered onto the central star. Our notation will be
slightly simplified if we introduce the vector magnetic flux components

$\displaystyle\Psi_{r}\equiv rA_{r}\ ,\ \Psi_{\theta}\equiv rA_{\theta}\ ,\
\Psi_{\phi}\equiv r\sin\theta\ A_{\phi}\ .$ (10)

In that notation,

$\displaystyle B_{r}$ $\displaystyle\equiv$
$\displaystyle\frac{1}{r^{2}\sin\theta}\left(\frac{\partial\Psi_{\phi}}{\partial\theta}-\frac{\partial\Psi_{\theta}}{\partial\phi}\right)\
,$ $\displaystyle B_{\theta}$ $\displaystyle\equiv$
$\displaystyle\frac{1}{r^{2}\sin\theta}\left(\frac{\partial\Psi_{r}}{\partial\phi}-r\frac{\partial\Psi_{\phi}}{\partial
r}\right)\ ,$ $\displaystyle B_{\phi}$ $\displaystyle\equiv$
$\displaystyle\frac{1}{r}\frac{\partial\Psi_{\theta}}{\partial
r}-\frac{1}{r^{2}}\frac{\partial\Psi_{r}}{\partial\theta}\ ,$ (11)

and eq. (8) becomes

$\displaystyle\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\\{-\alpha\left(\frac{\partial\Psi_{\phi}}{\partial\theta}-\frac{\partial\Psi_{\theta}}{\partial\phi}\right)\right.$
$\displaystyle+\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\\}\
,$ (12)
$\displaystyle\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\\{-\alpha\left(\frac{\partial\Psi_{r}}{\partial\phi}-r\frac{\partial\Psi_{\phi}}{\partial
r}\right)\right.$ $\displaystyle-\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\\}\ ,$ (13)
$\displaystyle\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}$
$\displaystyle+\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)=$
$\displaystyle\frac{\sin\theta}{1-(r\sin\theta/R_{\rm
LC})^{2}}\left\\{-\alpha\left(\frac{\partial\Psi_{\theta}}{\partial
r}-\frac{1}{r}\frac{\partial\Psi_{r}}{\partial\theta}\right)\right.$
$\displaystyle+\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\\}$ (14)

From the r.h.s. of eq. (14), the reader can check directly that the general
regularization condition at the light cylinder $r\sin\theta=R_{\rm LC}$
becomes

$\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}\ .$ (15)

Eq. (15) determines the value of the function $\alpha$ 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 $\alpha=0$.

##  3 Machine Learning

We propose to solve the steady-state pulsar magnetosphere problem using
Machine Learning techniques, and in particular by training a custom Neural
Network (hereafter NN). In § 4 below we will introduce several innovations
that will lead to a better solution than what was obtained before in both 2D
and 3D. Our present formulation is general, but our method will also be
implemented in axisymmetry in § 5.

A NN that is used to solve the partial differential equations that describe a
physical problem is called a Physics Informed Neural Network (hereafter PINN).
The entries of such NN are basically the spatial coordinates $r,\theta,\phi$,
and the exits are $\Psi_{r}$, $\Psi_{\theta}$, $\Psi_{\phi}$ and $\alpha$. 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,\theta,\phi$, 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.

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_{*}\leq r\leq r_{\rm max}$, $0\leq\theta\leq\pi$, $0\leq\phi\leq
2\pi$ (here, $r_{*}$ is the radius of the central neutron star) and try to
satisfy all the following constraints on them:

  1. 1.

Eqs. (12)-(14) over the computational domain

  2. 2.

Eqs. (9) over the computational domain

  3. 3.

Eq. (15) at $r\sin\theta=R_{\rm LC}$ (it may not be necessary to impose this
condition since it is satisfied if (i) is satisfied)

  4. 4.

$\alpha=0$ along field lines that do not cross the light cylinder (we will
assume here an untwisted closed-line region)

  5. 5.

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-centered spherical
coordinate system $(r,\theta_{m},\phi_{m})$ along that direction, then

$\displaystyle\Psi_{r}(r_{*},\theta_{m},\phi_{m})$ $\displaystyle=$
$\displaystyle\Psi_{\theta_{m}}(r_{*},\theta_{m},\phi_{m})=0$
$\displaystyle\Psi_{\phi_{m}}(r_{*},\theta_{m},\phi_{m})$ $\displaystyle=$
$\displaystyle\Psi_{\rm max}\sin^{2}\theta_{m}$ (16)

From the above star-centered components
$\Psi_{r},\Psi_{\theta_{m}},\Psi_{\phi_{m}}$ of the vector potential on the
stellar surface, we will obtain the actual components
$\Psi_{r},\Psi_{\theta},\Psi_{\phi}$ 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, the
NN will manage to satisfy all of them by minimizing a global loss function. As
we will see below, it is not necessary to impose boundary conditions at the
outer radial boundary if we use rescaled variables (e.g. $q\equiv 1/r$) such
that the outer radial boundary lies at infinity. 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.

What are the optimal number of internal nodes and training steps, the optimal
activation functions and NN optimization parameters are not known, and are
determined after several trial trainings of the problem. 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.

##  4 Innovations

The general 3D PINN numerical procedure described in § 3 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 $r\sin\theta<R_{\rm 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. This problem will be addressed in subsection 4.2 below.

###  4.1 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. The two domains are separated by
the separatrix surface characterized by the line where it originates on the
surface of the star (see figure 1). In principle, one may obtain an analytic
expression for the shape of the separatrix surface, namely

$r_{\rm S}=r_{\rm S}(\theta,\phi)\ .$ (17)

Obviously, this shape is unknown and must be determined. Let’s assume for the
moment that we choose ad hoc some initial form of $r_{\rm S}$ in eq. (17) that
originates around the magnetic poles of the central star at inclined polar
angles $\theta_{m}=\theta_{\rm pc}$, $0\leq\phi_{m}\leq 2\pi$ and
$\pi-\theta_{m}=\theta_{\rm pc}$, $0\leq\phi_{m}\leq 2\pi$. Here, as before,
$\theta_{m},\phi_{m}$ are spherical coordinates centered around the inclined
magnetic axis $\hat{m}$ of the central star. Notice that for angles
$\theta_{m}<\theta_{\rm pc}$ and $\theta_{m}>\pi-\theta_{\rm pc}$, the radius
connecting the point $(r,\theta,\phi)$ with the origin does not cross the
separatrix. Thus, in the open line region, the domain of integration becomes
$r\geq r_{*}$, $0\leq\theta_{m}\leq\theta_{\rm pc}$ and $\pi-\theta_{\rm
pc}\leq\theta_{m}\leq\pi$, and $r\geq r_{\rm S}$, $\theta_{\rm
pc}\leq\theta_{m}\leq\pi-\theta_{\rm pc}$. In the closed line region, the
domain of integration becomes $r_{*}\leq r\leq r_{\rm S}$, $\theta_{\rm
pc}\leq\theta_{m}\leq\pi-\theta_{\rm pc}$.

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 eq. (1)
accross current sheets and obtained that

$(B^{2}-E^{2})_{\rm below}=(B^{2}-E^{2})_{\rm above}\ .$ (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. 2012). We propose here a
better way to treat such discontinuities:

  1. 1.

Solve eqs. (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 $\theta_{m}=\theta_{\rm pc}$ and $\theta_{m}=\pi-\theta_{\rm
pc}$ on the stellar surface that corresponds to

$\Psi_{\rm dipole}\equiv\Psi_{\rm
max}\frac{r_{*}}{r}\sin^{2}\theta_{m}\equiv\Psi_{\rm max}\sin^{2}\theta_{\rm
pc}\ ,$ (19)

which yields

$r_{\rm S}\equiv r_{*}\frac{\sin^{2}\theta_{m}}{\sin^{2}\theta_{\rm pc}}\
\mbox{if}\ \theta_{\rm pc}\leq\theta_{m}\leq\pi-\theta_{\rm pc}$ (20)

Obviously, for such an ad hoc shape of the separatrix surface, the quantity
$B^{2}-E^{2}$ will in general be discontinuous accross that surface.

  2. 2.

Iteratively adjust the shape of the separatrix surface $r=r_{\rm
S}(\theta,\phi)$ at all points $(\theta,\phi)$ so as to satisfy the condition
that $B^{2}-E^{2}$ 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^{2}-E^{2}$ at each point of the separatrix in the adjacent closed
and open-line regions, and displace $r_{\rm S}(\theta_{m},\phi_{m})$
proportionally to the difference between the two corresponding values in the
direction of the smaller value. 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^{2}-E^{2}$ is continuous
everywhere accross the separatrix.

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_{\rm S}(\theta,\phi)$ (thick red line).
Open-line region: outside the separatrix ($r>r_{\rm S}(\theta,\phi)$). 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).

###  4.2 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 that is contained within the separatrix surface inside the
light cylinder. Since now the open-line region is independent from the closed-
line region, there is no problem to numerically invert the direction of the
field lines that leave the star from the northern polar cap and return from
infinity on the southern polar cap. With this numerical trick, all field lines
in the open-field region outflow from both polar caps, the field discontinuity
along the equatorial current sheet numerically ‘disappears’, and the position
of the equatorial current sheet coincides with the surface of flux direction
inversion in the open 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.

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

$\displaystyle\Psi_{\phi_{m}}(r_{*},\theta_{m},\phi_{m})$ $\displaystyle=$
$\displaystyle\Psi_{\rm max}\sin^{2}\theta_{m}\ \mbox{for}\
0\leq\theta_{m}\leq\theta_{\rm pc}\ ,$
$\displaystyle\Psi_{\phi_{m}}(r_{*},\theta_{m},\phi_{m})$ $\displaystyle=$
$\displaystyle 2\Psi_{\rm S}-\Psi_{\rm max}\sin^{2}\theta_{m}$ (21)
$\displaystyle=$ $\displaystyle\Psi_{\rm max}(2\sin^{2}\theta_{\rm
pc}-\sin^{2}\theta_{m})$ $\displaystyle\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \
\ \mbox{for}\ \pi-\theta_{\rm pc}\leq\theta_{m}\leq\pi$

Here $\theta_{\rm pc}$ is the opening of the polar cap. Obviously, in the
open-line region on the surface of the star, $0\leq\Psi_{\phi_{m}}\leq
2\Psi_{\rm S}$, whereas in the closed-line region on the surface of the star,
$\Psi_{\rm S}\leq\Psi_{\phi_{m}}\leq\Psi_{\rm max}$ ($\Psi_{\rm
S}\equiv\Psi_{\rm max}\sin^{2}\theta_{\rm pc}$). This mathematical trick
simplifies tremendously the numerical treatment of the large-scale open-line
region, and there is no problem with the closed-line region because the two
regions are in contact along the separatrix where all flux functions defined
in eqs. (10) are the same from both regions. After the solution is obtained,
the field will be reversed in the open-line region originating in the southern
polar cap, and the true magnetic field configuration will be obtained with the
equatorial current sheet.

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 $\theta$ approaches the equatorial current sheet at
$\theta=\theta_{\rm eq.c.s.}$ from above, the first term in the expression for
the field component $B_{r}$ in eq. (11), namely
$\partial\Psi_{\phi}/\partial\theta/(r^{2}\sin\theta)$, reaches in general a
non-zero value that is immediately reversed below the equatorial current
sheet. In other words,

$\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.}^{+})\neq 0\ .$ (22)

Unfortunately, from the symmetry of $\Psi_{\phi}$ 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 $\partial\Psi_{\phi}/\partial\theta$ 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
$\partial\Psi_{\phi}/\partial\theta(r,\theta\rightarrow\theta_{\rm
eq.c.s.})\neq 0$, which, after the final field reversal, naturally satisfies
eq. (22).

Finally, by chosing the angular size $\theta_{\rm pc}$ of the polar cap, this
approach allows as to choose the amount of magnetic flux $\Psi_{\rm
S}\equiv\Psi_{\rm max}\sin^{2}\theta_{\rm pc}$ that will open to
infinity.111Notice that we can even choose a more general non-circular polar
cap by specifying a non-uniform polar cap angular distribution $\theta_{\rm
pc}(\phi_{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 $\phi$. 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 $\theta_{\rm pc\ south}(\phi_{m})=\pi-\theta_{\rm
pc\ north}(\phi_{m}$). This effect is caused because field lines in the
corotating closed-line region start and end along the same star-centered
meridional angle $\phi_{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. 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. We checked that the condition $E<B$ is
automatically satisfied in the open-line region beyond the light cylinder.

##  5 First results in axisymmetry

###  5.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
($\partial/\partial\phi=0$) where a steady-state solution is known (CKF). In
the case of the aligned rotator it is customary to use the simpified notation
$\Psi_{\phi}\equiv\Psi$. In units where $R_{\rm LC}\equiv 1$, eqs. (11) yield

$\displaystyle B_{r}$ $\displaystyle\equiv$
$\displaystyle\frac{1}{r^{2}\sin\theta}\frac{\partial\Psi}{\partial\theta}\ ,$
$\displaystyle B_{\theta}$ $\displaystyle\equiv$
$\displaystyle-\frac{1}{r\sin\theta}\frac{\partial\Psi}{\partial r}\ ,$
$\displaystyle B_{\phi}$ $\displaystyle\equiv$
$\displaystyle\frac{I(\Psi)}{r\sin\theta}\ ,$ (23)

and eq. (14) becomes

$\displaystyle\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]$
$-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$ (24)

Eq. (24) is the well known pulsar equation (Wagoner & Schalermann 1982) in
spherical coordinates. Here, $I(\Psi)=r\sin\theta B_{\phi}$, and
$I^{\prime}(\Psi)\equiv{\rm d}I/{\rm d}\Psi=\alpha$. The condition that
$I=I(\Psi)$ (eq. 9) may be rewritten as

$\frac{\partial\Psi}{\partial r}\frac{\partial
I}{\partial\theta}-\frac{\partial\Psi}{\partial\theta}\frac{\partial
I}{\partial r}=0$ (25)

$I(\Psi)\neq 0$ only in the open-line part of the magnetosphere. We consider
here an untwisted magnetosphere with $I(\Psi)=0$ in its closed-line part. Our
results may also be generalized in the case of a twisted magnetosphere (see
Ntotsikas, Gourgouliatos & Contopoulos 2023). In axisymmetry, the
regularization condition at the light cylinder (eq. 15) becomes

$\displaystyle II^{\prime}$ $\displaystyle=$ $\displaystyle 2(\cos\theta\
B_{r}-\sin\theta\ B_{\theta})|_{\rm LC}\equiv 2B_{z}|_{\rm LC}\ .$ (26)

As before, eq. (26) determines the value of the function $II^{\prime}(\Psi)$
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^{\prime}(\Psi)=0$.

Eq. (24) was first solved by Contopoulos, Kazanas & Fendt (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.

###  5.2 Machine learning

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

$\displaystyle q$ $\displaystyle\equiv$ $\displaystyle\frac{1}{r}$
$\displaystyle\mu$ $\displaystyle\equiv$ $\displaystyle\cos\theta$

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

$\displaystyle(q^{2}-(1-\mu^{2}))\left(q^{2}\frac{\partial^{2}\Psi}{\partial
q^{2}}+(1-\mu^{2})\frac{\partial^{2}\Psi}{\partial\mu^{2}}\right)$
$\displaystyle+2q^{3}\frac{\partial\Psi}{\partial
q}+2(1-\mu^{2})\mu\frac{\partial\Psi}{\partial\mu}+II^{\prime}=0$ (27)

In order to help the solution satisfy the boundary conditions along the axis
and the stellar surface, namely $\Psi(r,\theta=0,\pi)=I(r,\theta=0,\pi)=0$ and
$\Psi(r_{*},\theta)=\sin^{2}\theta\ \Psi_{\rm max}$, we define new functions
$f,{\cal I}$ as the new exits of the NN such that

$\displaystyle\Psi$ $\displaystyle\equiv$ $\displaystyle(1-\mu^{2})f$
$\displaystyle I$ $\displaystyle\equiv$ $\displaystyle(1-\mu^{2}){\cal I}$
(28)

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

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

while the condition that $I=I(\Psi)$ (eq. 25) may be rewritten as

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

The source term in eq. (29) is equal to

$\displaystyle{\cal S}$ $\displaystyle\equiv$ $\displaystyle{\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)$ (31) $\displaystyle=$
$\displaystyle{\cal I}\frac{\sqrt{q^{2}(1-\mu^{2})\left(\frac{\partial{\cal
I}}{\partial q}\right)^{2}+\left(2\mu{\cal I}-(1-\mu^{2})\frac{\partial{\cal
I}}{\partial\mu}\right)^{2}}}{\sqrt{q^{2}(1-\mu^{2})\left(\frac{\partial
f}{\partial q}\right)^{2}+\left(2\mu f-(1-\mu^{2})\frac{\partial
f}{\partial\mu}\right)^{2}}}$ $\displaystyle\cdot{\rm
sgn}\left(q^{2}(1-\mu^{2})\frac{\partial{\cal I}}{\partial q}\frac{\partial
f}{\partial q}\right.$ $\displaystyle\ \ \ \ \ \ +\left.\left(2\mu{\cal
I}-(1-\mu^{2})\frac{\partial{\cal I}}{\partial\mu}\right)\left(2\mu
f-(1-\mu^{2})\frac{\partial f}{\partial\mu}\right)\right)$

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

  1. 1.

Eq. (29)

  2. 2.

Eq. (30)

  3. 3.

$f=1$ on the stellar surface

  4. 4.

${\cal I}=0$ along the equator inside the light cylinder

###  5.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 $\theta=0$,
and therefore $\theta_{m}\equiv\theta$ and $\phi_{m}\equiv\phi$ (see figure
2). The open line region is $r\geq r_{\rm S}$, $0\leq\theta\leq\pi$. The
closed line region is $r_{*}\leq r\leq r_{\rm S}$, $\theta_{\rm
pc}\leq\theta\leq\pi-\theta_{\rm pc}$. Here, the initial choice of the
separatrix radius $r_{\rm S}$ is

$r_{\rm S}(\theta)=r_{*}\frac{\sin^{2}\theta}{\sin^{2}\theta_{\rm pc}}\ \ \ \
\ \mbox{if}\ \theta_{\rm pc}<\theta<\pi-\theta_{\rm pc}$ (32)

The shape of $r_{\rm S}(\theta)$ 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

$\displaystyle(q^{2}-(1-\mu^{2}))\left(q^{2}(1-\mu^{2})\left(\frac{\partial
f}{\partial q}\right)^{2}+4\mu^{2}f^{2}\right.$
$\displaystyle+\left.(1-\mu^{2})^{2}\left(\frac{\partial
f}{\partial\mu}\right)^{2}-4\mu(1-\mu^{2})f\frac{\partial
f}{\partial\mu}\right)+(1-\mu^{2}){\cal I}^{2}$ (33)

accross the separatrix as described in subsection 4.1 above. In the latter
expression, ${\cal I}=0$ in the closed-line region, and ${\cal I}\neq 0$ in
the open-line region accross the separatrix.

We want to impose that $\Psi(r_{\rm S},\theta)=\Psi_{\rm S}=\Psi_{\rm
max}\sin^{2}\theta_{\rm pc}$ in the separatrix between the two domains. We,
therefore, introduce one more loss function in both domains, namely one that
requires

$f\left(q=\frac{1}{r_{\rm S}},\mu\right)=\Psi_{\rm
max}\frac{\sin^{2}\theta_{\rm pc}}{1-\mu^{2}}\ .$ (34)

###  5.4 Field reversal

The boundary conditions for the magnetic field in the closed-line region
$r_{*}\leq r\leq r_{\rm S}(\theta),\theta_{\rm pc}<\theta<\pi-\theta_{\rm pc}$
are

$\displaystyle\Psi(r_{*},\theta)$ $\displaystyle=$ $\displaystyle\Psi_{\rm
max}\sin^{2}\theta\ ,$ $\displaystyle\Psi(r_{\rm S}(\theta),\theta)$
$\displaystyle=$ $\displaystyle\Psi_{\rm S}\ .$ (35)

As before, $\theta_{\rm pc}$ is the opening of the polar cap, and $\Psi_{\rm
S}\equiv\Psi_{\rm max}\sin^{2}\theta_{\rm pc}$. We implement the mathematical
field reversal in the open-line region outside the closed-line region with
boundary conditions

$\displaystyle\Psi(r_{*},\theta)$ $\displaystyle=$ $\displaystyle\Psi_{\rm
max}\sin^{2}\theta\ \ \ \ \ \ \ \ \ \ \ \mbox{for}\ 0\leq\theta\leq\theta_{\rm
pc}\ ,$ $\displaystyle\Psi(r_{*},\theta)$ $\displaystyle=$ $\displaystyle
2\Psi_{\rm S}-\Psi_{\rm max}\sin^{2}\theta$ $\displaystyle=$
$\displaystyle\Psi_{\rm max}(2\sin^{2}\theta_{\rm pc}-\sin^{2}\theta)\ \
\mbox{for}\ \pi-\theta_{\rm pc}\leq\theta\leq\pi\ \ $
$\displaystyle\Psi(r_{\rm S}(\theta),\theta)$ $\displaystyle=$
$\displaystyle\Psi_{\rm S}\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \
\mbox{for}\ \theta_{\rm pc}<\theta<\pi-\theta_{\rm pc}\ .$ (36)

Obviously, in the open-line region, $0\leq\Psi\leq 2\Psi_{\rm S}$, whereas in
the closed-line region, $\Psi_{\rm S}\leq\Psi\leq\Psi_{\rm 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 $\Psi(r,\theta)=\Psi_{\rm S}$.
After the solution is obtained, the field will be reversed in the open-line
region where $\Psi_{\rm S}\leq\Psi\leq 2\Psi_{\rm S}$, and the true magnetic
field configuration will be presented with the equatorial current sheet along
$\Psi=\Psi_{\rm S}$.

As we discussed above, one extra benefit of this approach is that
$B_{r}(r>1,\theta=\pi/2^{-})$ is naturally allowed to be non-zero and equal to
$-B_{r}(r>1,\theta=\pi/2^{+})$ accross the equatorial current sheet. However,
without the mathematical field reversal that we propose in the open-line
region, due to symmetry,
$B_{r}=(\partial\Psi/\partial\theta)/(r^{2}\sin\theta)=0$ along the equator
$\theta=\pi/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, Ntotsikas & Gourgouliatos (2023) 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).

###  5.5 First results

We present below our first preliminary results for $\theta_{\rm
pc}=1.13(r_{*}/R_{\rm LC})^{1/2}$ and $\Psi_{\rm open}\equiv\Psi_{\rm
S}=\Psi_{\rm max}\sin^{2}\theta_{\rm pc}$ (see figure 3). This particular
value of $\theta_{\rm 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 $\Psi_{\rm dipole\ LC}=\Psi_{\rm max}r_{*}/R_{\rm LC}$, and
$\theta_{\rm dipole\ pc}\approx(r_{*}/R_{\rm LC})^{1/2}$. In previous high-
resolution solutions of the pulsar equation (e.g. Timokhin 2006), the magnetic
field line that crosses the light cylinder corresponds to $\Psi_{\rm LC\
Timokhin}=1.23\Psi_{\rm dipole\ LC}$ and $\theta_{\rm pc\ Timokhin}\approx
1,10(r_{*}/R_{\rm LC})^{1/2}$. Therefore, this particular value of
$\theta_{\rm pc}$ was chosen so that the Y-point lies very close to the light
cylinder. Instead, it was found that $R_{\rm Y}\approx 0.76R_{\rm LC}$ which
is closer to $R_{\rm dipole}(\Psi_{\rm S})=r_{*}/\sin^{2}\theta_{\rm
pc}\approx 0.82R_{\rm LC}$. This unexpected result differs from all previous
numerical solutions and needs further investigation. 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, and $B_{r}(r>R_{\rm Y},\theta=\pi/2^{-})=-B_{r}(r>R_{\rm
Y},\theta=\pi/2^{+})\neq 0$ (superscripts $-/+$ are used to denote ‘right
above/below the equator’ respectively). The above observations merit further
investigation in a future publication.

We implemented three independent PINNs:

  1. 1.

One with two entries $q,\mu$, three hidden layers with 256 nodes each, and two
exits $f,{\cal I}$. This PINN solves the pulsar equation in the open-line
region outside the separatrix. It encounters no problems at the light cylinder
where the functional form $I|_{\rm LC}=I(\Psi|_{\rm LC})$ is determined. Its
training requires tens of thousands of steps, mainly because the condition
that $I=I(\Psi)$ 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.

  2. 2.

One with two entries $q,\mu$, two hidden layers with 256 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.

  3. 3.

The most crucial part of our method is the readjustment of the shape of the
separatrix. We choose random angles $\theta_{\rm pc}<\theta<\pi-\theta_{\rm
pc}$ where we know the values of $B^{2}(r_{\rm S},\theta)-E^{2}(r_{\rm
S},\theta)$ 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 separatric must be adjusted at those
angles. We implemented the adjustment

$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}^{+}}}$ (37)

at those angles $\theta$. 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 $\beta$ is an adjustable positive parameter. After
the new position of the separatrix is defined on these randomly chosen angles
$\theta_{\rm pc}<\theta<\pi-\theta_{\rm pc}$, a third NN with one entry $\mu$,
two hidden layers with 128 nodes each, and one exit $r_{\rm S}$ is trained to
yield the shape of the separatrix at all angles $\theta$ required by the first
two main PINNs.

In practice, the two main PINNs are trained for 5,000 steps, then 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 whole
process is repeated 20 times, which correspond to about 100,000 training steps
for the two main PINNs. In the resulting configuration, pressure balance is
achieved accross the separatrix to within 1%.

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

Figure 3: Representative axisymmetric solution for $\theta_{\rm
pc}=1.13(r_{*}/R_{\rm LC})^{1/2}$ with domain decomposition and field reversal
in the southern hemisphere. Top: poloidal magnetic field lines ($\Psi$ and
$I(\Psi)$ isocontours). Thick line: separatrix $\Psi_{\rm open}\equiv\Psi_{\rm
S}=\Psi_{\rm max}\sin^{2}\theta_{\rm pc}$. 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. Bottom left: electric current
distribution $I(\Psi)R_{\rm LC}/\Psi_{\rm S}$ accross open magnetic field
lines $0\leq\Psi/\Psi_{\rm S}\leq 1$. Dotted line: analytic expression for
split monopole $\Psi/\Psi_{\rm S}(1-\Psi/\Psi_{\rm S})$. Notice that
$I^{\prime}(\Psi_{\rm S})=0$ in agreement with eq. (15) at $r=R_{\rm LC}$ and
$\theta=\pi/2$. Bottom right: the distribution of $B_{r}(r,\theta)/B_{r}(r,0)$
for $r=R_{\rm LC}$. Notice that
$B_{r}(r,\theta=\pi/2^{-})=-B_{r}(r,\theta=\pi/2^{+})\neq 0$ (without field
reversal, $B_{r}(r,\theta=\pi/2)$ would erroneously be equal to zero).

##  6 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
separatric 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.

We have obtained a preliminary solution of the pulsar equation for a
particular value of the open magnetic flux, namely $\Psi_{\rm open}=\Psi_{\rm
max}\sin^{2}\theta_{\rm pc}$ with $\theta_{\rm pc}=1.13(r_{*}/R_{\rm
LC})^{1/2}$, in which the closed-line region is found to extend up to
$0.76R_{\rm LC}$. We would like to correct here an important misunderstanding
found in the literature. With the solution obtained in § 5 above, 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, Ntotsikas &
Gourgouliatos 2023 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 explains the placement of the Y-point inside the light
cylinder in all high-resolution global PIC simulations of the last decade.

In summary, although Stefanou, Pons et al. (2023b) were the first to solve the
pulsar equation with a PINN and their work was an inspiration for our present
work, our implementation differs. 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 followup 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 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.

## Data availability statement

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

## References

  * Blandford (2002) Blandford, R. 2002, ESO Astroph. Symp., M. Gilfanov, R. Sunyaev, E. Churazov Eds., Springer-Verlag, 381
  * Bogovalov (1999) Bogovalov, S. 1999, Astron. Astrophys. , 349, 1017
  * Cerutti et al. (2020) Cerutti, B., Philippov, A., A. & Dubus, G. 2020, Astron. Astrophys. , 642, 204
  * Cerutti et al. (2022) Soudais, A. & Cerutti, B. 2022, An. Meet. French Soc. of Astr. & Astroph., 53
  * Contopoulos, Kazanas & Fendt (1999) Contopoulos, I., Kazanas, D. & Fendt, C. 1999, ApJ, 511, 351 (CKF)
  * Contopoulos, Nathanail & Pugliese (2014) Contopoulos, I., Nathanail, A. & Pugliese, D. 2014, ApJ, 780L, 5
  * Contopoulos, Ntotsikas & Gourgouliatos (2023) Contopoulos, I., Ntotsikas, D., & Gourgouliatos, K. 2023, Mon. Not. Roy. Astron. Soc. (submitted)
  * Endean (1974) Endean 1974, Astrophys. J., 187, 359
  * Goldreich & Julian (1969) Goldreich & Julian 1969, Astrophys. J., 157, 869
  * Gourgouliatos & Lynden-Bell (2019) Gourgouliatos, K. & Lynden-Bell, D. 2019, Mon. Not. Roy. Astron. Soc. , 482, 1942
  * Gruzinov (1999) Gruzinov, A. 1999, arXiv:astro-ph/9902288
  * Gruzinov (2005) Gruzinov, A. 2005, Phys. Rev. Let., 94, 021101
  * Hakobyan et al. (2023) Hakobyan, H., Philippov, A. & Spitkovsky, A. 2023, Astrophys. J., 943, 105
  * Hewish et al. (1968) Hewish, A. et al. 1968, Nature, 217, 709
  * Hu & Beloborodov (2022) Hu, R. & Beloborodov, A. M. 2022, Astrophys. J., 939, 42
  * Kalapotharakos & Contopoulos (2009) Kalapotharakos, C. & Contopoulos, I. 2009, Astron. Astrophys. , 496, 495
  * Kalapotharakos et al. (2012) Kalapotharakos, C., Harding, A. K., Kazanas, D. & Contopoulos, I. 2012, Astrophys. J., 754, 1
  * Kalapotharakos et al. (2018) Kalapotharakos, C., Brambilla, G., Timokhin, A., Harding, A. K. & Kazanas, D. 2018, ApJ, 857, 44
  * Kalapotharakos et al. (2022) Kalapotharakos, C., Wadiasingh, Harding, A. K. & Kazanas, D. 2022, Astrophys. J., 934, 65
  * Lyubarsky (1990) Lyubarskii, Yu. E. 1990, Soviet Astron. Lett., 16, 16
  * Lyubarsky (2005) Lyubarsky Y. E. 2005, Mon. Not. Roy. Astron. Soc. , 358, 113
  * Michel (1982) Michel, F. C. 1982, Rev. Mod. Phys., 54, 1
  * Mestel (1975) Mestel, L. 1975, Soc. Roy. des Sciences de Leège, 8, 79
  * Muslimov & Harding (2009) Muslimov, A. & Harding, A. K. 2009, Astrophys. J., 692, 140
  * Ntotsikas et al. (2023) Ntotsikas, D., Gourgouliatos, C. N. & Contopoulos, I. 2023 (submitted)(
  * Parfrey et al. (2012) Parfrey, K., Beloborodov, A. & Hui, L. 2012, Mon. Not. Roy. Astron. Soc. , 423, 1416
  * Pétri (2012) Pétri, J. 2012, Mon. Not. Roy. Astron. Soc. , 424, 605
  * Pétri (2018) Pétri, J. 2018, Mon. Not. Roy. Astron. Soc. , 477, 1035
  * Philippov & Spitkovsky (2014) Philippov, A. & Spitkovskyt, A. 2014, Astrophys. J., 785, L33
  * Philippov et al. (2015) Philippov, A., Spitkovskyt, A. & Cerutti, B. 2015, Astrophys. J., 801, L19
  * Richards & Lyutikov (2018) Richards , G.& Lyutikov, M. 2018, arXiv:1801.05839
  * Spitkovsky (2006) Spitkovsky, A. 2006, Astrophys. J., 648, L51
  * Stefanou et al. (2023) Stefanou, P., Pons, J. A. & Cerdá-Durán, P. 2023a, Mon. Not. Roy. Astron. Soc. , 518, 6390
  * Stefanou et al. (2023) Stefanou, P., Pons, J. A. et al. 2023b, submitted
  * Tchekhovskoy et al. (2013) Tchekhovskoy, A., Spitkovsky, A. & Li, J. G. 2013, Mon. Not. Roy. Astron. Soc. , 435, L1
  * Timokhin (2006) Timokhin, A. 2006, Mon. Not. Roy. Astron. Soc. , 368, 1055
  * Uzdensky (2003) Uzdensky, D. 2003, Astrophys. J., 598, 446
  * Werner et al. (2018) Werner, G. R., Uzdensky, D. A., Begelman, M. C., Cerutti, B., & Nalewajko, K. 2018, Mon. Not. Roy. Astron. Soc. , 473, 4840