Solving the Pulsar Equation using Physics-Informed Neural Networks
In this study, Physics-Informed Neural Networks (PINNs) are skilfully applied
to explore a diverse range of pulsar magneto-spheric models, specifically
focusing on axisymmetric cases. The study successfully reproduced various
axisymmetric models found in the literature, including those with non-dipolar
configurations, while effectively characterizing current sheet features. Energy
losses in all studied models were found to exhibit reasonable similarity,
differing by no more than a factor of three from the classical dipole case.
This research lays the groundwork for a reliable elliptic Partial Differential
Equation solver tailored for astrophysical problems. Based on these findings,
we foresee that the utilization of PINNs will become the most efficient
approach in modelling three-dimensional magnetospheres. This methodology shows
significant potential and facilitates an effortless generalization,
contributing to the advancement of our understanding of pulsar magnetospheres.
Petros Stefanou, Jorge F. Urbán, José A. Pons
# Solving the Pulsar Equation using Physics-Informed Neural Networks

###### Abstract

In this study, Physics-Informed Neural Networks (PINNs) are skilfully applied to explore a diverse range of pulsar magnetospheric models, specifically focusing on axisymmetric cases. The study successfully reproduced various axisymmetric models found in the literature, including those with non-dipolar configurations, while effectively characterizing current sheet features. Energy losses in all studied models were found to exhibit reasonable similarity, differing by no more than a factor of three from the classical dipole case. This research lays the groundwork for a reliable elliptic Partial Differential Equation solver tailored for astrophysical problems. Based on these findings, we foresee that the utilization of PINNs will become the most efficient approach in modelling three-dimensional magnetospheres. This methodology shows significant potential and facilitates an effortless generalization, contributing to the advancement of our understanding of pulsar magnetospheres.

keywords: magnetic fields; pulsars; stars: neutron;

## 1 Introduction

Physics-Informed Neural Networks (PINNs) (Lagaris et al. (1997); Raissi et al. (2019)) is a relatively new but very promising family of PDE solvers based on Machine Learning (ML) techniques. This method is suitable to obtain solutions of Partial Differential Equations (PDEs) describing the physical laws of a given system by taking advantage of the very successful modern ML frameworks and incorporating physical knowledge. In recent years, PINN-based solvers have been used to solve problems in a great variety of fields: fluid dynamics (Cai et al. (2021)), turbulence in supernovae (Karpov et al. (2022)), radiative transfer (Korber et al. (2023)), black hole spectroscopy (Luna et al. (2023)), cosmology (Chantada et al. (2023))large scale structure of the universe (Aragon-Calvo (2019)), galaxy model fitting (Aragon-Calvo & Carvajal (2020)), inverse problems (Pakravan et al. (2021)) and many more. In our previous work (Urban et al. (2023), Paper I hereafter) we presented a PINN solver for the Grad-Shafranov equation, which describes the magnetosphere of a slowly rotating neutron star endowed with a strong magnetic field (a magnetar) in the axisymmetric case. In that paper, we demonstrated the ability of the network to be trained for various boundary conditions and source terms simultaneously. In this work, our purpose is to extend our PINN approach to the more general - and more challenging - case of rapidly rotating neutron stars (pulsars). The rotating case presents new interesting challenges related to the presence of current sheets in the magnetosphere. Our implementation is able to deal with these pathological regions sufficiently well, demonstrating its potential for problems where classical methods struggle.

Pulsar magnetospheres have been studied extensively in the last 25 years with various approaches, each with its advantages and limitations. Contopoulos et al. (1999) were the first to solve the axisymmetric, time -independent problem, a result that was later confirmed and improved by Gruzinov (2005) and Timokhin (2006). Spitkovsky (2006) used a full MHD time-dependent code and was able to acquire solutions for aligned and oblique rotators. Solutions for arbitrary inclination were also obtained by Petri (2012) who used a time - dependent pseudo-spectral code. More recent approaches involve large-scale Particle-in-Cell simulations that include the influence of accelerated particles (Cerutti et al. (2015); Philippov & Spitkovsky (2018)). All of these, and many related, works have improved our understanding of the pulsar magnetosphere. However, some important questions remain still without an answer.

The structure of the paper is the following. After a brief summary of the relevant equations and boundary conditions in section SS2, we will describe the PINN method in SS3, with emphasis in the new relevant details with respect to the magnetar problem (paper I). In section SS4 we present our results, showing that we are able to reproduce the well-established axisymmetric results encountered in the literature, but also indicating that new, possibly unexplored solutions can be encountered. We summarize our most relevant conclusions and discuss possible improvements and future extensions in section SS5.

## 2 Pulsar Magnetospheres.

Our aim is to extend our previous results from paper I to find numerical solutions of the pulsar equation, that can be written as follows:

\[\mathbf{\nabla}\times\left(\mathbf{B}-\beta^{2}\mathbf{B}_{p}\right)=\alpha\mathbf{B}, \tag{1}\]

where \(\beta=v/c=|\Omega\sigma|/c=\varpi/\kappa_{\rm LC}\) is the co-rotational speed in units of the speed of light, with \(\kappa_{\rm LC}=c/\alpha\) being the light-cylinder (LC)
radius and \(\mathbf{B}_{P}=\mathbf{B}-\mathbf{B}_{\phi}\) is the magnetic field perpendicular to the direction of rotation. Here, \(\alpha\) is a scalar function given by

\[\alpha=\frac{4\pi}{c}\left(\mathbf{J}-\rho_{e}\mathbf{v}\right)\cdot\frac{\mathbf{B}}{B^{2}}. \tag{2}\]

representing the ratio between the field-aligned component of the current in the corotating frame and the local magnetic field strength. We refer to the comprehensive and thorough recent review by Philippov & Kramer (2022) (and references therein) for a complete historical overview of magnetospheric physics and the mathematical derivation of the equations.

As in Paper I, we focus in the axisymmetric case and, for convenience, we use compactified spherical coordinates \((q,\mu,\phi)\) where \(q=1/r,\mu=\cos\theta\). In these coordinates, any axisymmetric magnetic field can be written in terms of a poloidal and a toroidal scalar stream function \(P\) and \(T\) as

\[\mathbf{B}=\frac{q}{\sqrt{1-\mu^{2}}}\left(\mathbf{\nabla}P\times\hat{\phi}+T\hat{\phi }\right). \tag{3}\]

Here \(P\) is related to the magnetic flux and is constant along magnetic field lines. Plugging this expression into Eq. (1) and taking the toroidal component, we get

\[\mathbf{\nabla}P\times\mathbf{\nabla}T=0, \tag{4}\]

which means that \(T\) is only a function of \(P\) (\(T=T(P)\)) and, therefore, is also constant along magnetic field lines. The poloidal component of Eq. (1) gives the well-known Pulsar Equation (Michel (1973); Scharlemann & Wagoner (1973), which in our coordinates reads

\[\left(1-\beta^{2}\right)\Delta_{\text{GS}}P+2\beta^{2}q^{2}\left(q\partial_{q }P+\mu\partial_{\mu}P\right)+G(P)=0, \tag{5}\]

where \(G(P)=TT^{\prime}\) and \(\Delta_{\text{GS}}\) is the Grad-Shafranov operator, given by

\[\Delta_{\text{GS}}\equiv q^{2}\partial_{q}\left(q^{2}\partial_{q}\right)+ \left(1-\mu^{2}\right)q^{2}\partial_{\mu\mu}. \tag{6}\]

Notice that in the limiting case where \(\beta=0\), we ignore the rotationally induced electric field and Eq. (5) is reduced to the Grad-Shafranov equation. Hereafter we will use the following shorthand notation for the extended operator

\[\Delta_{GS\beta}\equiv\left(1-\beta^{2}\right)\Delta_{\text{GS}}+2\beta^{2}q^ {2}\left(q\partial_{q}+\mu\partial_{\mu}\right). \tag{7}\]

A crucial difference in the modelling of magnetar and pulsar magnetospheres is the source of poloidal currents \(TT^{\prime}\). In magnetars, currents are assumed to be injected into the magnetosphere due to the magnetic field evolution in the NS's crust, which twists the magnetosphere (Akgun et al., 2018). Therefore, the source term in Eq. (5) is given either as a user-parametrised model (as for example in Akgun et al. (2016) for the 2D problem or in Stefanou et al. (2023) for the 3D problem) or in a self-consistent manner by coupling the interior evolution with the magnetosphere. We adopted the latter approach in Paper I (see the astrophysical application in section 5 of that paper), where a series of magnetospheric steady-state solutions were obtained for each time-step of the internal magneto-thermal evolution.

In the pulsar magnetosphere, however, the source of current is the LC. Lines that cross the LC have to open up and bend, giving rise to azimuthal fields and currents. We refer to the region with open field lines as the open region. On the other hand, lines that do not cross the LC, but turn back to the surface, co-rotate rigidly with the the star. We refer to this region as the closed region. The source term \(TT^{\prime}\) must be determined self-consistently to ensure smooth crossing of the lines along the LC. In particular, at the LC where \(\beta=1\), Eq. (5) takes the simple form

\[-q^{2}\left(q\partial_{q}P+\mu\partial_{\mu}P\right)=2B_{z}=TT^{\prime}, \tag{8}\]

which places a constraint on the possible functions \(T(P)\). This makes Eq. (5) complicated to solve, as two functions have to be determined simultaneously. In order to do so, additional physical boundary conditions have to be imposed (see next subsection).

Through this work, we measure distances in units of the radius of the star \(R\) and magnetic fields in units of the surface magnetic field strength at the equator \(B_{0}\) (note that the surface magnetic field strength at the poles is \(2B_{0}\)). In these units, the magnetic flux is measured in units of the total poloidal flux \(P_{0}=B_{0}R^{2}\). It is convenient to additionally define the total magnetic flux carried by field lines that cross the LC in a non-rotating dipole \(P_{1}=P_{0}(R/R_{\text{LC}})\), which will be useful in what follows. Finally, the toroidal stream function is measured in units of \(B_{0}R\).

### Boundary Conditions

We assume that the magnetic field at the surface of the star \(P(q=1)\) is known. For example, in the case of a dipole

\[P(q=1,\mu)=B_{0}(1-\mu^{2}) \tag{9}\]

but we will also explore other options. The last closed field line, called separatrix, will be labelled by \(P=P_{c}\). It marks the border between open and closed (current-free) regions. The value \(P_{c}\) corresponds to the total magnetic flux that crosses the light cylinder. The point where this line meets the equator is the called Y-point and it should lie at a radius \(r_{c}\), somewhere between the surface and the LC (see e.g. Timokhin (2006) for a detailed discussion). Far away from the surface, field lines should become radial, resembling a split monopole configuration (Michel, 1973). Inside the closed region, the magnetosphere is current-free and purely poloidal. No toroidal fields are developed, so that

\[T(P>P_{c})=0 \tag{10}\]

by definition.

In the classical model, a current sheet should develop along the separatrix, supported by the discontinuity of the toroidal magnetic field between closed (\(B_{\phi}=0\)) and open (\(B_{\phi}\neq 0\)) regions. Another current sheet should exist along the equator and beyond the Y-point, supported by the opposed directions of the magnetic field lines between the two hemispheres. The current sheet forms the return current that accounts for the current supported by out-flowing particles along the open field lines and closes the current circuit.

All this motivated the seminal works to impose equatorial symmetry by defining a numerical domain only on one hemisphere and imposing another boundary condition \(P=P_{c}\) at the equator. To reproduce these models, we will impose that along the equator and beyond the Y-point

\[P(q<q_{c},\mu=0)=P_{c}. \tag{11}\]

Here \(q_{c}=1/r_{c}\) is the corresponding location of the Y-point in compactified coordinates. However, later works have relaxed this assumption (Contopoulos et al., 2014), proposing a more general solution without a separatrix and with a transition region close to the equator. We will explore both cases in our results section.

### Total Energy and Radiated Power

The total energy of the electromagnetic field in the magnetosphere is given by

\[\mathcal{E}=\frac{1}{8\pi}\int\left(B^{2}+E^{2}\right)dV. \tag{12}\]
It is convenient to describe the electromagnetic energy content of a given model in terms of the excess energy of a particular magneto-spheric solution with respect to the non-rotating dipole,

\[\Delta\mathcal{E}=\frac{\mathcal{E}-\mathcal{E}_{d}}{\mathcal{E}_{d}}, \tag{13}\]

where

\[\mathcal{E}_{d}=\frac{1}{3}B_{0}^{2}R^{3} \tag{14}\]

is the total magnetic energy of a non-rotating dipole.

The total power radiated away by a rotating magnetosphere can be calculated by integrating the Poynting flux over a sphere far away from the star

\[\hat{\mathcal{E}}=\frac{c}{4\pi}\int\left(\mathbf{E}\times\mathbf{B}\right)\cdot\mathbf{f} d\omega=-\frac{c}{R_{\mathrm{LC}}}\int_{0}^{P_{c}}TdP, \tag{15}\]

where \(\omega\) is the solid angle. Again, it is convenient to measure the relative difference of the radiated power with respect to the classical order of magnitude estimate

\[\hat{\mathcal{E}}_{d}=\frac{B_{0}^{2}R^{6}}{R_{\mathrm{LC}}^{4}}c, \tag{16}\]

and hereafter we will express \(\hat{\mathcal{E}}\) in units of \(\mathcal{E}_{d}\).1

Footnote 1: We should stress that the (unfortunately) often used \(\sin^{2}\theta\) dependence with the inclination angle in \(\hat{\mathcal{E}}_{d}\) is unphysical. It misleads to conclude that an aligned rotator does not emit, but it actually does radiate in a similar amount (but a factor of about 2) than the orthogonal rotator.

## 3 Network Structure and Training Algorithm

In Paper I, we used a PINN to calculate approximate solutions of the axisymmetric Grad-Shafranov equation. We refer the reader to that paper for a more detailed description of the generic implementation. In this section, we will briefly summarise the main points and highlight the novelties introduced to adapt our solver to the pulsar case. The principal changes introduced aim at enforcing the physical constraint (boundary conditions, or the \(T(P)\) requirement) by construction instead of leaving the job to minimize additional terms in the loss function, which usually does not reach the required accuracy.

We consider solutions in points \((q,\mu)\in\mathcal{D}\) in a 2-dimensional domain \(\mathcal{D}\). We denote by \(\partial\mathcal{D}\) the boundary of this domain. To account for the equatorial constraint in Eq. (11), we will consider the equatorial line beyond \(r_{c}\) as part of \(\partial\mathcal{D}\). In order to ensure that \(P\) depends on the coordinates while \(T\) is solely a function of \(P\), we design a network structure consisting of two sub-networks that are trained simultaneously. The output of each sub-network depends only on its corresponding input

The first sub-network takes the coordinates \((q,\mu)\) as input and returns as output a function that we denote by \(N_{P}\). Then, \(P\) is calculated using

\[P(q,\mu)=f_{b}(q,\mu)+h_{b}(q,\mu)N_{P}(q,\mu;\Theta). \tag{17}\]

where \(f_{b}\) can be any function in \(\mathcal{D}\) that satisfies the corresponding boundary conditions at \(\partial\mathcal{D}\). \(h_{b}\) is an arbitrary function representing some measure of the distance to the boundary, that must vanish at \(\partial\mathcal{D}\). Both user-supplied functions \(f_{b}\) and \(h_{b}\) depend only on the coordinates and are unaffected by the PINN. There is some freedom to decide their specific form, as long as they retain certain properties (e.g. they are sufficiently smooth) and they have the desired behaviour at the boundary. The only part that is adapted during training is \(N_{P}\) through its dependence on the trainable parameters \(\Theta\). With this _parametrisation_ (or _hard enforcement_) approach we ensure that boundary conditions are exactly satisfied by construction. This differs from the other usual approach, consisting of adding more terms related to the boundary conditions in the loss function. We will specify the particular form of the functions \(f_{b}(q,\mu)\) and \(h_{b}(q,\mu)\) in the next section when we discuss different cases.

Next, the second sub-network takes \(P\) as input and returns \(N_{T}\) as output. This automatically enforces that \(T=T(P)\), required by Eq. (4). The new network output \(N_{T}(P;\Theta)\) is used to construct the function \(T(P)\) as follows:

\[T(P)=g(P)N_{T}(P;\Theta)\, \tag{18}\]

where \(g\) is another user-supplied function to include other physical restrictions. For example, we can use a step function to treat a possible discontinuity of \(T\) at the border between open and closed regions (Eq. (10)). Another advantage of this approach is that \(T^{\prime}\) is calculated from \(T\) via automatic differentiation in the second sub-network.

A subtle point in the pulsar problem is the determination of the critical value \(P_{c}\) that separates the open and closed regions. We have explored two different approaches: **(a)**\(P_{c}\) is self-determined by the network and **(b)**\(P_{c}\) is fixed to some value. In the first case, \(P_{c}\) is considered an extra output of the first sub-network, closely connected to \(P\), resembling the approach taken by classical solvers for this problem. In contrast, the second method involves fixing \(P_{c}\) as a hyperparameter to a predetermined value, leaving the network to discover a solution that aligns with this fixed value. This introduces the necessity to loosen some of the other imposed constraints, like the position of the Y-point, and encourages the exploration of solution properties from a fresh standpoint.

During the network training process, both \(P\) and \(T\) (and possibly \(P_{c}\), if we follow the case **(a)** above) are obtained simultaneously by minimising the following loss function

\[\mathcal{L}=\omega_{1}\mathcal{L}_{PDE}+\omega_{2}\sigma(P_{c}), \tag{19}\]

where

\[\mathcal{L}_{PDE}=\frac{1}{N}\sum_{(q,\mu)\in\mathcal{D}}\left[\Delta_{GS \beta}P(q,\mu)+G(P(q,\mu))\right]^{2}. \tag{20}\]

Here, \(N\) represents the size of the training set, and \(\omega_{1}\) and \(\omega_{2}\) are adjustable parameters. In approach **(b)**, \(\omega_{2}=0\), whereas in approach **(a)**, \(\omega_{2}\) is adjustable. Additionally, \(\sigma\) denotes the standard deviation of the \(P_{c}\) values for the points in the training set. \(P_{c}\) is a global magnetospheric value that should not depend on the coordinates. However, during training, there is no guarantee that this value will be the same for all the points that are considered. To ensure that the values of \(P_{c}\) at arbitrary points are as close to each other as possible, we minimise the standard deviation (\(\sigma(P_{c})\)) of the values of \(P_{c}\) at all the points of the training set by including an additional term in the loss function (19). This inclusion guarantees that the network's output for \(P_{c}\) remains constant and independent of the coordinates \((q,\mu)\). A schematic representation of our network's structure can be found in Fig. 1.

We consider a fully-connected architecture for both sub-networks consisting of 4 hidden layers for the first one and 2 hidden layers for the second one. All layers have 40 neurons. Our training set consists of \(N=5000\) random points \((q,\mu)\in\mathcal{D}\), which are periodically changed every 500 epochs in order to feed the network with as many distinct points as possible. The total number of epochs is 35000. This number may seem big at first sight, but is necessary because of the
amount of points considered. In order to minimise the loss function, we use the ADAM optimisation algorithm (Kingma & Ba, 2014) with an exponential learning rate decay. We have also considered in this work the idea of introducing trainable activation functions, suggested in PINNs originally by Jagtap et al. 2020. For each hidden layer \(k\), the linear transformation performed at that layer is multiplied by a trainable parameter \(c_{k}\). We found that this practice can accelerate convergence but a more rigorous study is out of the scope of this paper.

## 4 Magnetospheric Models

In this section, we present the results obtained for various cases and under different physical assumptions. Our analysis encompasses the successful reproduction of all distinct axisymmetric models found in literature, which include:

1. The classical solutions (Contopoulos et al., 1999; Gruzinov, 2005),
2. The family of solutions with varying locations of the Y point (Timokhin, 2006),
3. The improved solution where the separatrix current sheet is smoothed out (Contopoulos et al., 2014),
4. The non-dipolar solutions (Gralla et al., 2016).

Each of these models presents unique challenges and characteristics, and we will delve into the outcomes achieved for each one. In particular, the overall magnetospheric configuration, the poloidal flux at the separatrix \(P_{c}\), the functions \(T(P)\) and \(G(P)\), and the energy losses \(\mathcal{E}\) agree with the previous works.

Figure 1: A sketch of the network structure. Two sub-networks are employed to ensure that \(P=P(q,\mu)\) and \(T=T(P)\).

Figure 3: Same as figure 2 but with the Y-point positioned at \(r_{c}=0.4R_{\rm LC}\). The region where the return current flows thrurough open field lines in negligible.

Figure 2: The classical axisymmetric pulsar magnetosphere. Fade black lines: magnetic field lines as contours of \(P\). Thick pink line: separatrix and equatorial current sheet, where \(P=P_{c}\). Vertical green line: light cylinder. Colourmap: the source current \(G(P)\). Labels on the field lines are in units of \(P_{c}\). Colourbar is in symmetrical logarithmic scale. z and x axes are in units of the stellar radius \(R\). The bulk of the return current flows along the current sheet, with a small percentage flowing along open field lines
### Hard enforcement of boundary conditions.

As discussed earlier, we construct \(P\) according to Eq. (17) to fulfill the boundary conditions. For the cases (i) and (ii) we have employed an \(f_{b}\) with the following form:

\[f_{b}(q,\mu)=(1-\mu^{2})\left[P_{c}+q(1-P_{c})\frac{\text{ReLU}^{3}\left(1- \frac{q_{c}}{q}\right)}{(1-q_{c})^{3}}\right], \tag{21}\]

where ReLU (x) is the _Rectified Linear Unit_ function, which returns zero if \(x<0\) and \(x\) if \(x>0\)2, \(q_{c}=\nicefrac{{1}}{{r_{c}}}\) is the position of the Y-point and \(n\) is a free positive parameter.

Footnote 2: Alternative definitions of \(\text{ReLU}(x)\) are \(x\mathcal{H}(x)\) (\(\mathcal{H}\) being the Heaviside step function) and \(\text{max}(0,x)\).

Furthermore, we have chosen the following \(h_{b}\) function:

\[h_{b}(q,\mu)=(1-\mu^{2})(1-q)\sqrt{\text{ReLU}^{3}\left(1-\frac{q_{c}}{q} \right)+\mu^{2}}, \tag{22}\]

where the reader can check that \(h_{b}\) is zero at the boundary. The terms with the ReLU functions enforce that \(P=P_{c}\) at the equator when \(q<q_{C}\), and the third power ensures that \(P\) and its first and second derivatives are all continuous.

The parametrisation for \(T(P)\) requires the additional function \(g\). We use a gaussian3transition beginning at \(P=P_{c}\):

Footnote 3: Any other function that asymptotically approaches the Heaviside step function could be used instead.

\[g(P)=\begin{cases}1&P\leq P_{c}\\ e^{-\frac{(P-P_{c})^{2}}{2(\mathcal{H}P)^{2}}}&P>P_{c}\end{cases} \tag{23}\]

where \(\delta P\) is a small number that controls the width of the current sheet, where the transition of T from a finite value to zero takes place. Note that \(g\) and its first derivative are both continuous at the separatrix, so \(T\) and \(T^{\prime}\) are well defined.

For the cases (iii) and (iv), the condition \(P=P_{c}\) at the equator is lifted and consequently, the parametrisation functions must be modified. Our choice for these cases is:

\[f_{b}(q,\mu)=f_{1}(q)\sum_{l=1}^{l_{\text{max}}}\frac{b_{l}}{l}P _{l}^{\prime}(\mu)+f_{2}(q)(1-|\mu|) \tag{24}\] \[h_{b}(q,\mu)=q(1-q)(1-\mu^{2})\] (25) \[g(P)=-P\ \left(\text{ReLU}\left(1-\frac{|P|}{P_{c}}\right) \right). \tag{26}\]

Adjustable weighting functions, denoted as \(f_{1}\) and \(f_{2}\), are utilized to control the relative importance of the two terms. To be specific, \(f_{1}\) should be significant near the surface but diminish at considerable distances from the star, while \(f_{2}\) should vanish close to the surface but be dominant at \(q\ll 1\). At intermediate distances, approximately around the LC, both \(f_{1}\) and \(f_{2}\) should be much smaller than \(h_{b}\) to ensure that the neural network contribution in Eq. (17) dominates.

In order to allow for more versatile configurations beyond the standard dipole representation, the surface boundary condition is presented as a linear combination of magnetic multipoles. To impose a current-free region, the function \(g\) is employed, and while it does not necessarily indicate a current sheet, it is chosen to be at least quadratic. This ensures that \(T^{\prime}\) and consequently \(G\) experience at least one change of sign, allowing the current circuit to close.

### Classical solutions

The solutions obtained in Contopoulos et al. (1999); Gruzinov (2005) will be referred to as the _classical solutions_. Fig. 2 illustrates contour lines of the poloidal flux function \(P\), with the colormap representing the source current \(G=TT^{\prime}\). All the well-known characteristics of the pulsar magnetosphere are observed in these solutions: The open field lines extend beyond the light cylinder and stretch towards infinity, eventually adopting a split monopole configuration at considerable distances, and a current sheet forms at the equator due to the magnetic field's reversal between the north and south hemispheres.

Figure 4: **(a)** Colourmap of the residuals of the pulsar equation. Colourbar is in log scale. **(b)** Evolution of the loss function with the training epochs. Only the value of the loss every 100 epochs is plotted for clarity. Big spikes correspond to changes in the training set. Small fluctuations can be interpreted as the variance of the PDE residual.
In the region where field lines have \(P>P_{c}\), both the toroidal magnetic field and poloidal current are zero, except for the narrow transition zone \([P_{c},P_{c}+\delta P]\) located just inside the separatrix. Beyond the separatrix, the blue region illustrates the small portion of the return current that flows back to the star along open field lines. In contrast, the nearly white region just inside the separatrix represents the substantial portion of the return current that flows along the current sheet. We summarise in Tab. 1 the model parameters of different solutions. The first line corresponds to the classical model depicted in Fig. 2.

To showcase the accuracy of our findings, we present in Fig. 4a how well our solution aligns with the pulsar equation (5). The color map illustrates the absolute error of the pulsar equation for our model, revealing remarkably low values within the bulk of the domain \((\lesssim 10^{-5})\). As expected, the error is slightly larger in proximity to the separatrix, a region of discontinuity. Nonetheless, this discrepancy does not impede our solver from effectively approximating the solution throughout the rest of the domain, nor does it significantly affect the solution far from these regions.

In general, we anticipate the maximum error to be of the order of \(\sim\sqrt{\mathcal{L}}\) since it corresponds to the error of the PDE for a considerable set of random points. Indeed, as depicted in Fig. 4b, at the conclusion of the training process, the loss reaches values around \(\sim 10^{-7}\), confirming this assumption. The prominent spikes observed in the figure correspond to the periodic changes of the training set of points. However, it is noticeable that as the training epochs progress, the spikes diminish, indicating that the network has learned to generalize to new points without compromising accuracy. The small fluctuations between the spikes, should be interpreted as the variance in the approximation to the solution. As the network adapts its parameters to acquire the solution, it cannot simultaneously reconcile all the training points. Consequently, it fluctuates around a mean instead of finding a single minimum value.

With the same PINN, it is straightforward to produce solutions with different positions of the Y-point, as in Timokhin (2006) simply by varying the parameter \(q_{c}\) in Eq. (21). Fig. 3 shows an example of such a solution, where the Y-point lies at a distance \(r_{c}=0.4R_{\text{LC}}\). In this case, the totality of the return current flows through the current sheet. Interestingly, we observe that the luminosity for \(r_{c}=0.4R_{\text{LC}}\) is an order of magnitude larger than the classical solution (see Table 1).

As discussed in Timokhin (2006), both the luminosity and the total energy stored in the magnetosphere increase with decreasing \(r_{c}\), so the magnetosphere will generally try to achieve the configuration with the minimum energy, that is, with \(r_{c}\) as close as possible to the light cylinder. Therefore, although the configurations with a small \(r_{c}\) are mathematically sound and very interesting from the astrophysical point of view (much larger luminosity), they are probably short-lived and less frequent in nature than the standard configuration.

### Non-fixed boundary condition at the equator

In most previous studies that used classical methods to tackle this problem, it was a common practice to solve it in just one hemisphere and setting boundary conditions at the equator (Eq. (11)), where a kind of "jump" is expected to occur. This approach is reasonable, but it limits the range of magnetospheric solutions that can be obtained. The solutions have equatorial symmetry (in addition to axisymmetry) and end up having a specific configuration: a dipole magnetic field at the surface and a Y-point configuration where the equatorial and separatrix current sheets meet.

However, by utilizing PINNs and taking advantage of their local and flexible nature for imposing constraints and boundaries, we can create solutions with fewer (only physically relevant, rather than mathematical) requirements. In this section, we introduce solutions where boundary conditions are applied only at the surface and at infinity, without restricting the equator as part of the solution domain.

At infinity, we simply demand that the solution approaches a split monopole configuration (last term in Eq. (24)) with a specific value (denoted by \(P=P_{\infty}\)) at the equator:

\[P(\mu,q=0)=P_{\infty}(1-|\mu|)\.\]

Since the equator is free from any boundary conditions, some other constraint must be imposed to distinguish between a possibly infinite class of different solutions. We decide to set the value of \(P_{c}\) beforehand. This value separates the regions with and without electric currents. This approach is just as valid and perhaps more versatile than other constraints, like pinning down the position of the Y-point. As for \(P_{\infty}\), we do not fix its value; instead we let the network figure it out. The value of \(P_{\infty}\) separates the regions with currents of different sign.

In Fig. 5 we present a typical solution for these boundary conditions. This magnetospheric configuration bears resemblance to the one obtained in Contopoulos et al. (2014), where a substantial number of field lines that cross the light cylinder close inside the equatorial current sheet. However, we have arrived to this result by following a different prescription. They enforced specific boundary conditions at the equator to ensure that the perpendicular component of the Lorentz force applied to the equatorial current sheet becomes zero. On the other hand, our approach involved imposing that the solution becomes a split monopole with a certain magnetic flux at infinity while leaving the equator unrestricted.

An interesting generalisation of this set of solutions involves applying non-dipolar surface boundary conditions. Instead of sticking to the basic dipole case, the surface magnetic field can be a combination of various magnetic multipoles (Gralla et al., 2016). This option is not feasible in the classical solution because having contributions of even multipoles breaks the equatorial symmetry. Nevertheless, our solver can handle this situation without any issues, because it is

Figure 5: Same as figure 2 but for the case of non-fixed equatorial boundary condition. Field lines that cross the light cylinder are close through the current sheet. The solid pink line indicates the line \(P=P_{c}=1.4P_{1}\), whereas the dashed pink line corresponds to \(P=P_{\infty}=0.77P_{1}\).
working throughout the entire domain. As we move to significant distances from the star, we anticipate the multipole solution to gradually converge towards the classical dipole solution.

Fig. 6 shows an example of such a case. In particular, we consider a combination of a dipole, a quadrupole, and an octupole, with the corresponding coefficients in Eq. (24) being \(b_{1}=1,b_{2}=-2\) and \(b_{3}=2\). At first glance, when observing the magnetosphere on a larger scale of tens of LC, it appears quite similar to the previous one. However, upon closer inspection of the zoomed-in right panel, the complex and intricate structure in the innermost region of the magnetosphere becomes evident.

To better understand the structure of this last model, in Fig. 7 we show the absolute value of the first six multipole weights computed over spheres at different radii. We see how the even \(l=2\) multipole quickly decreases with distance. Other even multipoles grow in the inner region, but they also become vanishingly small at a distance of a few LC. On the contrary, the odd multipoles approach the asymptotic values corresponding to the split monopole at radial infinity (\(P\propto 1-|\mu|\)), shown with dashed lines in the figure.

## 5 Discussion and final remarks

Table 1 summarises the various models studied in this work. We use the letter \(C\) to denote solutions acquired following the classical approach, the letter \(E\) for solutions without an equatorial constraint and the letter \(M\) for solutions with multipolar content. For the first three models, \(P_{c}\) is determined by the PINN, while for the rest of them is fixed to a predetermined value. \(P_{\infty}\), when present, is always determined by the PINN.

The excess energy for all the models is of the order of a few percent. Models with no equatorial restriction tend to have lower energies, indicating that they are the preferred configuration in nature. With the exception of the model C3, the luminosity for all models does not vary from the classical dipole case by more than a factor of three.

In Fig. 8, we present the function \(T(P)\) for all the models of Table 1. The curves corresponding to models with fixed boundary conditions at the equator exhibit a steep transition to zero, which is modelled by the current sheet (Eq. (23)) and begins at \(P=P_{c}\). In contrast, the models with no restrictions at the equator smoothly reach the value \(P=P_{c}\) through a continuous transition from positive (orange area in Figs. 5, 6) to negative current (blue area in Figs. 5, 6). Importantly, these models do not develop a current sheet to close the current circuit.

The area under each of the curves represents the luminosity (see Eq. (15)), which explains why the model with a lower \(r_{c}\) exhibits a significantly higher energy loss rate, as it has more open lines carrying Poynting flux to infinity. Additionally, the presence or absence of equatorial restrictions leads to distinct behaviours in the models. In the models without equatorial restrictions, not all the Poynting flux carried by lines crossing the light cylinder (LC) escapes to infinity. Instead, up to 50% of the total pulsar spindown energy flux (the area between \(P_{\infty}\) and \(P_{c}\)) remains confined in the equatorial current sheet. This trapped energy could potentially be dissipated in particle acceleration and high-energy electromagnetic radiation within a few times the light cylinder radius, as was first pointed out in Contopoulos et al. (2014). However, the exact fraction of power that can be locally dissipated and reabsorbed, as opposed to the fraction that is genuinely lost and contributes to the spin-down, remains unclear. In the last column of Table 1 we have included the expected range of values for such models.

In this study, extending the results of our previous paper for slowly-rotating magnetar magnetospheres, our success lies in the proficient application of Physics-Informed Neural Networks (PINNs) to obtain numerical results for a diverse range of pulsar magnetospheric models, meticulously exploring various cases under different physical assumptions. Our extensive analysis encompasses the accurate reproduction of several distinct axisymmetric models found in the scientific literature. The decision to focus solely on axisymmetry was a deliberate one, aiming at affirming the remarkable ability of PINNs to effectively capture and characterize the intricate peculiarities exhibited by pulsar magnetospheres, including the features of the current sheet. Throughout our investigation, we purposefully accounted for a rich diversity of models, incorporating non-dipolar configurations.

This work serves as a stepping stone towards the development of a robust and trustworthy general elliptic Partial Differential Equation solver, specifically tailored to address the challenging complexities of this and related astrophysical problems. Looking ahead, our research can be naturally extended to the three-dimensional magnetospheric case, a promising prospect that holds the potential for reaching a deeper understanding of the underlying physics governing pulsar magnetospheres.

## Acknowledgements

We acknowledge the support through the grant PID2021-127495NB-I00 funded by MCIN/AEI/10.13039/501100011033 and by the European Union, the Astrophysics and High Energy Physics programme of the Generalitat Valenciana ASFAE/2022/026 funded by MCIN and the European Union NextGenerationEU (PRTRC-C17.I1) and the Prometeo excellence programme grant CIPROM/2022/13. JFU is supported by the predoctoral fellowship UAFPU21-103 funded by the University of Alicante.

## Data availability

All data produced in this work will be shared on reasonable request to the corresponding author.
