The Pulsar Magnetosphere with Machine Learning: Methodology

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 (PINNs) 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 (FFE). 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.


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 so-⋆ E-mail: icontop@academyofathens.gr lution 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: (i) 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 2024).We still need to understand where the Y-point truly lies, and how fast it moves out as the pulsar spins down.
(ii) 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 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. 2012, Mahlmann et al. 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).
(iii) 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 − β rec ln(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, Philippov & Spitkovsky (2023), this logarithmic gradual electromagnetic energy dissipation applies to the full magnetosphere, not only to its undulating nonaxisymmetric component as Cerutti, Philippov & Dubus (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.
(iv) 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).
(v) 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 con-clusions 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 § 2. In § 3 we propose a novel numerical approach with machine learning that will allow us to obtain the solution in 3D.In § 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.In the final section § 5 we summarize our results and discuss the challenges and limitations encountered while working with Neural Networks.

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 where B and E are the magnetic and the electric field respectively, ρ e is the electric charge density, and J is the electric current 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 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 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 Let us define here the light cylinder radius R LC ≡ c/Ω.From eq. ( 5) above we obtain that in component form, where the subscript 'p' denotes a poloidal component.Following eq. ( 7), E • ∇ × E = 0, and the expression for J in eq. ( 2) simplifies considerably.Substituting everything back in eq. ( 6) we obtain Here, α is the force-free parameter that obeys the constraint 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 nonrotating (lab) frame.In the latter, the final configuration is rotating, i.e. time-dependent.All its complex features (current sheets, Ypoints, 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.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 A such that B = ∇ × A. We will work in spherical coordinates (r, θ, ϕ) centered onto the central star.Our notation will be slightly simplified if we introduce the vector magnetic flux components In that notation, and eq.( 8) becomes From the r.h.s. of eq. ( 14), the reader can check directly that the general regularization condition at the light cylinder r sin θ = R LC becomes Eq. ( 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.

NOVEL NUMERICAL APPROACH
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 subsections 3.2 & 3.3 below 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.

Machine learning
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).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 et al. (2023Stefanou 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: (i) Eqs. ( 12)-( 14) over the computational domain (ii) Eqs. ( 9) over the computational domain (iii) Eq. ( 15) at r sin θ = R LC (it may not be necessary to impose this condition since it is satisfied if (i) is satisfied) (iv) α = 0 along field lines that do not cross the light cylinder (we will assume here an untwisted closed line region) (v) 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 m with a starcentered spherical coordinate system (r, θ m , ϕ m ) along that direction, then From the above star-centered components Ψ r , Ψ θm , Ψ ϕ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 eq. ( 16) are |Ψ r |, |Ψ θm |, and |Ψ ϕm − Ψ max sin 2 θ 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 § 4.5 below).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 § 4 below) we used 1,000 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 5,000 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 figure 5 below).
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 r sin θ < 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 subsection 3.3 below we propose yet another way to open up the magnetosphere.

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, etc.), while in others it extends to only about 85% 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 figure 1).In principle, one may obtain an analytic expression for the shape of the separatrix surface, namely r S = r S (θ, ϕ) . (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 eq. ( 17) that originates around the 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) across current sheets and obtained that 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, Hu & Beloborodov 2022, Mahlmann et al. 2022).We propose here a better way to treat such discontinuities: (i) 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 θ m = θ pc and θ m = π − θ pc on the stellar surface that corresponds to which yields Obviously, for such an ad hoc shape of the separatrix surface, the quantity B 2 −E 2 will in general be discontinuous across that surface.
(ii) Iteratively adjust the shape of the separatrix surface r = r S (θ, ϕ) at all points (θ, ϕ) 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 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 2 − E 2 is continuous everywhere across the separatrix.

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 What we have done here is to express the boundary conditions in terms of the vector magnetic flux components Ψ r , Ψ θm , Ψ ϕm in the inclined spherical system of coordinates (r, θ m , ϕ m ) in which a pure magnetic dipole is described as Ψ r = Ψ θm = 0 and Ψ ϕm = Ψ max sin 2 θ m .Obviously, in the open line region on the surface of the star, 0 ⩽ Ψ ϕm ⩽ 2Ψ S , whereas in the closed line region on the surface of the star, Ψ S ⩽ Ψ ϕm ⩽ Ψ max (Ψ S ≡ Ψ max sin 2 θ 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 eq. ( 11), namely ∂Ψ ϕ /∂θ/(r 2 sin θ), reaches in general a non-zero value that is immediately reversed below the equatorial current sheet.In other words, 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 eq. ( 22).Finally, by choosing the angular size θ pc of the polar cap, this approach allows as to choose the amount of magnetic flux Ψ S ≡ Ψ max sin 2 θ pc that will open to infinity. 1 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.
1 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-centered 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.

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, eqs.(11) yield and eq.( 14) becomes Eq. ( 24) is the well known pulsar equation (Wagoner & Schalermann 1982) in spherical coordinates.Here, I(Ψ) = r sin θB ϕ , and I ′ (Ψ) ≡ dI/dΨ = α.The condition that I = I(Ψ) (eq.9) may be rewritten as 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 (eq.15) becomes As before, eq. ( 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. 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.

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 In those new variables, eq. ( 24) may be rewritten as 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 2 θ Ψ max , we define new functions f, I as the new exits of the NN such that In these new functions, the pulsar equation (eq.27) may be rewritten as while the condition that I = I(Ψ) (eq.25) may be rewritten as The source term in eq. ( 29) is equal to In this reformulation of the problem, the entries of the PINN are the spatial coordinates (q, µ), and the exits are ( f, I).The loss functions must then satisfy the conditions (i) Eq. ( 29) (ii) Eq. ( 30) (iii) f = 1 on the stellar surface (iv) I = 0 in the closed line region Notice that there is no special provision along the light cylinder as expressed by eq. ( 15), since if eq. ( 29) is satisfied, condition ( 15) is also automatically satisfied.This is a welcome advantage of the PINN methodology.

Domain decomposition
The domain decomposition is simpler in the case of the aligned pulsar where the direction m of the stellar magnetic dipole is along θ = 0, and therefore θ m ≡ θ and ϕ m ≡ ϕ (see figure 2).The open line region is r Here, the initial choice of the separatrix radius The shape of r S (θ) will be adjusted iteratively in order to attain continuity of across the separatrix as described in subsection 4.1 above.In the latter expression, I = 0 in the closed line region, and I 0 in the open line region across the separatrix.We want to impose that Ψ(r S , θ) = Ψ S = Ψ max sin 2 θ pc in the separatrix between the two domains.We, therefore, introduce one more loss function in both domains, namely one that requires

Field reversal
The boundary conditions for the magnetic field in the closed line region r * ⩽ r ⩽ r S (θ), θ pc < θ < π − θ pc are As before, θ pc is the opening of the polar cap, and Ψ S ≡ Ψ max sin 2 θ pc .We implement the mathematical field reversal in the open line region outside the closed line region with boundary conditions 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 2 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, Ntotsikas & Gourgouliatos (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).

First results
We implemented three independent NNs with multiple internal (hidden) layers: (i) One PINN with two entries q, µ, three hidden layers with 64 nodes each, and two exits f, 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.
(ii) 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 I = 0.The training of this PINN is stable and fast.
(iii) 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 2 (r S , θ) − E 2 (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 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 highlyinclined 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 eq. ( 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 figure 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 figure 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) 2 over the stellar surface, of (Ψ − Ψ S ) 2 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, eqs.( 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 § 5).As we can see from figure 3, PDE errors (i.e. the residuals of eq.29) at the end of the training in the open line region are around 10 −2 , 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 −4 , and can certainly be improved in the future with further development of our method.
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 −2 ).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%.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 2 θ r * /r; figures 6 & 7).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 figures 4 and 5 for θ pc = 1.176(r * /R LC ) 1/2 and Ψ open ≡ Ψ S = Ψ max sin 2 θ 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 = Ψ max r * /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.Timkokhin 2006) the Y-point lies very close to the light cylinder.
We found that R Y = 0.88R LC which is outside R dipole (Ψ S ) = r * / sin 2 θ 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 eq. ( 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 In other words, in all previous numerical solutions there is no welldefined 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

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 day 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 behavior 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 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 = Ψ max sin 2 θ 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 § 4 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 2024 for a detailed analysis of this effect).This result differs from Gruzinov (2005) that the magnetospheric energy contained around the Ypoint 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 Ypoint 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.

APPENDIX
In order to test the accuracy of our method, we applied it to solve the simple magnetostatic problem ∇ × B = 0, or equivalently in our notation, with boundary condition Ψ(r * , θ) = Ψ max sin 2 θ on the stellar surface.In this problem there is no rotation, thus there is no poloidal electric current I. Nevertheless, we used the formal methodology of our present paper, namely domain decomposition and separatrix re-adjustement.It is also implied that we used a denser distribution of training points in the region of interest of the rotating problem inside the light cylinder.Notice that we chose the value Ψ S = 1.23Ψ dipole LC ad hoc since in this problem there is no light cylinder, thus there is no distinction between open and closed field lines as in the real pulsar problem.The initial separatrix position was taken to be non-dipolar in order to allow the method to obtain its final position (thick black line).In figure 6 we plot the outcome of the PINN Ψ ≡ (1 − µ 2 ) f (solid lines) and compare it with the exact solution Ψ dipole = Ψ max sin 2 θ r * /r (dashed lines).We overplot in color the accuracy of the solution defined as |(Ψ − Ψ dipole )|/Ψ dipole .We see that the accuracy of our method is better than 1.5% in the largest part of the magnetosphere.The two solutions are practically indistinguishable inside the distance that corresponds to the light cylinder of the rotating solution where the density of training points is higher (accuracy around 0.5%).Notice that in figure 6 we tested the accuracy of both the outcome of the PINN and the movement of the separatrix.In order to test the accuracy of the PINN itself we solve eq. ( 38) without a separatrix.This yields the solution shown in figure 7 which has even higher accuracy (better than 0.2%).We conclude that the numerical accuracy of our method is not on par yet with that of standard methods such as finite difference, finite volume or spectral methods.This is also the case in the PINN solution of Stefanou et al. (2023b).On the other hand, the numerical treatment of current sheets by standard methods it too is problematic (numerical current sheets have unphysical thicknesses and the numerical dissipation there is questionable).The purpose of this study has been the determination of the location of the current sheet and the imposition of pressure balance along that surface.This has never been successfully attained before in 3D, thus it is worthwhile to apply our methodology to the oblique rotator too.In order to improve our results in future iterations of this work, we may try to combine standard methods outside current sheets with our methodology of moving the separatrix.This is still work in progress.justment).Plot parameters as in figure 6.The numerical solution (solid lines) is indistinguishable from the analytic solution (red dotted lines).In the largest part of the magnetosphere the accuracy is well below 0.2%.
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 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 figure 1.

Figure 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).

Figure 2 .
Figure 2. Same as figure 1 in axisymmetry.Note that we have implemented artificial field-direction reversal in the open line region.
B 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 − Y , θ = π/2) that is required to satisfy pressure balance accross the T-point (the superscript − is used to denote 'right inside the Y-point'; figure 5 middle; see also figure 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; figure 5 bottom).The above observations merit further investigation in a future publication.

Figure 3 .
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': eq.(29), the pulsar equation (note that I = 0 in the closed line region); 'sep in'/ 'sep out': separatrix boundary condition; 'paral': eq.(30), I = I(Ψ).Abrupt loss changes correspond to updates in the training set.

Figure 4 .
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 = Ψ max sin 2 θ 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.

Figure 5 .
Figure 5. Top: electric current distribution I(Ψ)R LC /Ψ S across open magnetic field lines 0 ⩽ Ψ/Ψ S ⩽ 1 for the solution of figure 4. Dotted line: analytic expression for split monopole.Our result is in excellent agreement with figure 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 3 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 inwardx squeezing of the field right behind the Y-point that is required to satisfy pressure balance (in agreement with figure 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).

Figure 6 .
Figure 6.Solution to test the PINN methodology together with separatrix adjustment.Hyperparameters as in § 4.5.Ψ and separatrix iso-contours as in figure 4. Red dotted lines: magnetostatic dipole with known analytic solution.Color map: accuracy defined as |(Ψ − Ψ dipole )|/Ψ dipole .In order to test the convergence of the method, we initialized the surface Ψ = Ψ S = 1.23Ψ dipole LC at a position different from its final dipole position.Accuracy is better than 1.5% in the largest part of the magnetosphere and around 0.5% inside the distance that corresponds to the light cylinder of the rotating solution where the density of training points is higher.

Figure 7 .
Figure 7. Solution to test the PINN methodology alone (no separatrix ad-justment).Plot parameters as in figure6.The numerical solution (solid lines) is indistinguishable from the analytic solution (red dotted lines).In the largest part of the magnetosphere the accuracy is well below 0.2%.