Pulsar magnetospheres in General Relativity

The main contribution to the pulsar power can be calculated by assuming a rotating magnetically-dominated magnetosphere described by the force-free approximation. Although this simple model has been used thoroughly to study pulsar magnetospheres in the flat spacetime regime, only few works have considered the relativistic corrections introduced by the curvature and frame-dragging effects induced by a rotating neutron star. Here we revisit the problem and describe pulsar magnetospheres within full General Relativity, quantifying the corrections as a function of the angular velocity, the compactness of the star and the misalignment angle between the spin and the magnetic dipole. We provide analytical expressions for the pulsar luminosity by fitting our numerical results. Finally, we also analyze the effect of the relativistic corrections on the braking index, which indicates a slight increment in its value.

The main contribution to the pulsar power can be calculated by assuming a rotating magneticallydominated magnetosphere described by the force-free approximation. Although this simple model has been used thoroughly to study pulsar magnetospheres in the flat spacetime regime, only few works have considered the relativistic corrections introduced by the curvature and frame-dragging effects induced by a rotating neutron star. Here we revisit the problem and describe pulsar magnetospheres within full General Relativity, quantifying the corrections as a function of the angular velocity, the compactness of the star and the misalignment angle between the spin and the magnetic dipole. We provide analytical expressions for the pulsar luminosity by fitting our numerical results. Finally, we also analyze the effect of the relativistic corrections on the braking index, which indicates a slight increment in its value.

I. INTRODUCTION
Pulsars are bright sources of electromagnetic radiation, emitting from radio to gamma-ray frequencies. Even though the main picture -consisting on a rotating magnetized neutron star-is rather simple, a full solution that explain the diverse observational phenomenology still remains elusive due to the extreme conditions found on neutron stars. With typical radius R = 10 − 13 km and stellar masses between M ≈ 1.2 − 2.0 M 1 , they are very compact astrophysical objects, only exceeded by black holes. Moreover, their rotation periods ranges from seconds to millisecond, and their surface magnetic field intensities, inferred by timing properties, vary from ∼ 10 8 to ∼ 10 15 G (see McGill and ATNF pulsar catalogs 2 3 ).
In a seminal work, Goldreich & Julian [1] first demonstrated that the neutron star rotation induces electric fields which are strong enough to strip charged particles off the stellar surface. This mechanism and the pair production induced by the interaction of high-energy photons eventually populate the surrounding magnetosphere with a tenuous plasma. In such rarefied environments, the electromagnetic force dominates over particle inertia and leads to the force-free assumption, which is thought to hold everywhere except in small regions of space called gaps, where the observed radiation is produced by means of particle acceleration. This low-inertia limit of relativistic magnetohydrodynamics (MHD), known as forcefree electrodynamics (FFE), simplifies the problem considerably and has been widely used to study the global properties of pulsar magnetospheres. Assuming a dipo- * fcarrasc@famaf.unc.edu.ar † carlos.palenzuela@uib.es ‡ reula@famaf.unc.edu.ar 1 www3.mpifr-bonn.mpg.de/staff/pfreire/NS masses.html 2 www.physics.mcgill.ca/∼pulsar/magnetar/main.html 3 www.atnf.csiro.au/people/pulsar/psrcat/ lar magnetic field aligned with the rotation axis of the neutron star a canonical model emerged, beginning with a solution found by Contopoulos, Kazanas & Fendt [2] (CKF hereafter). The CKF solution confirmed the basic qualitative picture originally sketched by Goldrich & Julian; namely, the existence of a closed zone that co-rotates with the star, together with a polar outflow along open magnetic field lines that extends to infinity. These two regions are separated by thin return current layers emanating from the poles and meeting at the Ypoint, where the light cylinder intersects the equatorial plane. Beyond the Y-point there is a strong current sheet which extends along the equator up to infinity. After the CKF solution first appearance, other steady solutions were found [3][4][5], displaying a Y-point located at different positions within the light cylinder and thus raising the question about uniqueness of the CKF solution. Time-dependent force-free simulations later revealed that axi-symmetric solutions indeed converge to a CKF-type configuration. These simulations have been developed an extended by several authors using diverse numerical techniques. Both aligned [5][6][7][8][9] and oblique [10][11][12] pulsar magnetospheres were widely studied on flat spacetime. The effects of plasma pressure and inertia were further included (e.g. [13][14][15][16][17][18]); and the assumption of a dipolar magnetic field has been recently relaxed [19][20][21][22] by considering more general multipolar field configurations in electrovacuum. On the other hand, works including curvature effects due to general relativity (GR) are somewhat scarce. The first time-dependent simulations incorporating GR effects appeared in [23], reporting an approximate 20% enhancement of the aligned rotator spin-down luminosity relative to its flat spacetime value. Later, Ruiz et al. [24] presented a more systematic study of the GR aligned rotator, providing fitting formulas that shows the luminosity dependence on the star rotation and compactness. These works, although using different numerical codes, were both performed by using finite-difference schemes and matching the stellar ideal MHD interior to a force-free exterior. Recently, Petri [25] has employed a pseudo-spectral code with discontinuous Galerkin methods to evolve three-dimensional pulsar magnetospheres, modeling the NS by suitable boundary conditions at the stellar surface. This work thoroughly analyzes the spin-down luminosity dependence on the rotation rate and the misalignment angle between the spin and the magnetic dipolar moment, comparing Newtonian and GR results at a fixed compactness. More recently, there has been also an interesting analytical approach to the slowly rotating limit [26,27], using some tools from differential geometry such as the exterior calculus. In the last years there have been significant progress on particle-in-cell (PIC) simulations, which had allowed to reproduce global properties of pulsar magnetospheres, including self-consistently the regions of plasma production and acceleration (see e.g. [28][29][30]). One of the general lessons emerging from these PIC simulations is that pulsar environments are nearly force-free everywhere except for the thin return current layers and the current sheet outside the light cylinder, where particles are produced and part of the radiation is generated. Moreover, luminosities obtained within the PIC framework are consistent with previous force-free simulations; even for the oblique rotator case where the 1 + sin 2 χ dependence on the inclination angle χ has been retrieved [31]. Notice however that the inclusion of curvature -and in particular, the frame dragging effects-has shown to be crucial in allowing pair formation near the polar cap regions [32].
In this paper we aim to elucidate the role of relativistic and curvature effects on the spin-down luminosity of a pulsar, by performing full 3D time-dependent simulations of force-free pulsar magnetospheres within general relativity. The presence of the rotating neutron star is modeled by means of suitable boundary conditions at the stellar surface, while the exterior spacetime is approximated by the Kerr metric. It has been further assumed that the star posses a dipolar magnetic field as in most previous studies. Our numerical code, based on highorder finite differences schemes over a multiple patch infrastructure [33], has been extended from previous efforts in the context of black holes magnetospheres [34], e.g. [33]). The numerical domain is equipped with a Kerr metric (in appropriate coordinates) and accommodates quite naturally to the geometry of the problem. This code has been further developed in the present work, by including novel boundary conditions to represent the perfectly conducting surface of the NS. The method to deal with such boundary conditions relies on the penalty technique [35][36][37], and uses the characteristic decomposition of the force-free equations employed for the evolution. Although this numerical implementation is very different from those of previous studies like Ruiz et al. [24] and Petri [25], the results found on this paper are in good agreement in the regimes on which they overlap. As a final result, we provide a general expression describing the spin-down luminosity in terms of the three adimensional parameters that specify the pulsar: the spin rate of the neutron star, the stellar compactness and the misalignment angle between spin and magnetic dipole axis.
This article is organized as follows: In Section II we present the main aspects of our numerical approach and setup. Especial attention has been devoted to the treatment of boundary conditions at the stellar surface, while some related technical details -as well as other examples of application-were deferred to Appendix A. The main results are presented in Section III, starting with some tests that shows the correct implementation of boundary conditions and constraints behavior. We compare our numerical results with previous studies and then propose a generic formula for the pulsar spin-down luminosity, deduced from fittings of the numerical data. Conclusions and some perspectives for future projects are presented on Section IV. Throughout, we adopt geometrized units in which c = G = 1 and Lorentz-Heaviside units for the electromagnetic field.

II. SETUP
Here we are interested on modeling the magnetospheres of neutron stars by solving a particular version of force-free electrodynamics obtained at [38], which has some improved properties in terms of well posedness and involves the full force-free current density 4 . The numerical scheme to solve these equations is based on the multiblock approach [33,[35][36][37], where a particular multiple patch infrastructure has been equipped with the Kerr metric [33]. This provides a domain perfectly adapted to the geometry of the problem, having two global (inner and outer) boundaries with spherical topology. This implementation has been employed to perform accurate studies of force-free jets in black hole spacetimes [34]. Here we introduce however an important modification to the previous scheme: we develop detailed boundary conditions for the inner surface (placed inside the black hole in [34]), as to represent the perfect conducting surface of the star.
We shall start this section by briefly summarizing the generic features of the numerical approach. Then, we introduce the set of equations used to evolve the system. Finally, we describe how to deal with the boundary conditions for the stellar surface through the penalty method and the initial data for the magnetic field.

A. General Scheme
Our numerical domain consists on several touching grids (i.e. there is no overlap among them and only points at the boundaries are sheared), commonly referred as multi-block approach [33]. The equations are discretized at each individual subdomain by using difference operators constructed to satisfy summation by parts (SBP). In particular, we employ sixth-order accurate difference operators on the interior and third-order at the boundaries. Penalty terms [35][36][37] are added to the evolution equations at boundary points. These terms penalize possible mismatches between the different values the characteristic fields take at the interfaces, providing a consistent way of communicate information between the different blocks. More concretely, the penalization terms at the boundary points of a subdomain (say, "A") modify the evolution equations as, where U ν B are the values of the fields on the overlapping points from a neighbor subdomain "B" and the index "a" labels the different characteristic modes, being λ a their eigenvalues and P µ (a)ν (·) the projectors into their characteristic subspaces. Summation is performed only over the incoming modes (respect to subdomain A), λ a > 0. While h is the grid spacing along the normal direction and σ o the coefficient defining a discrete scalar product (see eqn. (2) below) valuated at the boundary. The projectors are build from the eigen-basis U µ For each incoming mode at one side of the interface there is an associated outgoing mode from the other side, and the penalty essentially enforces those values to coincide. At each subdomain, it is possible to find a semidiscrete energy defined by both a symmetrizer of the system at the continuum and a discrete scalar product with respect to which SBP holds, The summation by parts property of the operators implies an energy estimate, up to outer boundary and interface terms. The penalties are constructed such that their contribution to the energy estimate cancel inconvenient interface terms, thus providing an energy estimate which covers the whole integration region across grids. Such semidiscrete energy estimates guarantee the stability of the numerical scheme, provided an appropriate time integrator is chosen [41]. A fourth order Runge-Kutta algorithm is used for time integration in our code. We incorporate numerical dissipation through adapted Kreiss-Oliger operators [42], which are eight-order accurate on the interior and fourth-order at the boundaries.
Each subdomain is handed to a separate processor, while the information required for the interfaces treatment is communicated among them by the message passing interface (MPI) system. The computation for each grid may be, as well, parallelized by means of OpenMP.

B. Evolution Equations
We shall start from the covariant version of force-free electrodynamics for the electromagnetic field F ab and the Faraday tensor F * ab = abcd F cd /2, as presented in ref. [38]:F where the field was defined to extend the system outside of the constraint submanifold, G = 0, being G := F ab F * ab and F := F ab F ab the two electromagnetic invariants. Notice that eq. (3) reduces to the force-free condition, eq. (4) is the Faraday equation with an extra field φ to dynamically control the magnetic divergence-free constraint (see e.g. [43][44][45]), and eq. (5) is just the generalization of the constraint condition ∇ a G = 0.
A covariant hyperbolitation for these equations was found in [38], ensuring well posedness of the symmetrized system (equations (32)- (33) in [38]). The final evolution equations are then obtained after performing a standard 3+1 decomposition of such system. The line element in the adapted coordinates of the foliation is given by, with α,β i and γ ij being the lapse function, the shift vector and the intrinsic metric on the spatial slices, respectively.
We follow the conventions adopted in e.g. [46], where the electric and magnetic field are taken to be E a := F ab n b and B a := F * ab n b , being n a = (1, −β i )/α the normal to the spatial hypersurfaces. Notice the magnetic field differs in a sign respect to the usual definition, which was the one employed on [34,38]. The spatial Poynting vector is defined as, where ijk is the hypersurface induced volume element. The evolution system is finally written as follows: where we have denoted D j (·) ≡ 1 √ γ ∂ j ( √ γ · ) and To have further control on the constraint E · B = 0, we adopt a damping strategy taken from [47], with a moderate coefficient δ ∼ 100 to enjoy the constraint cleaning properties, while avoiding the complications of having stiff terms which would demand implicitexplicit schemes (as pointed out in [40]). In order to deal with current sheets we employ a "standard" approach in which electric field is effectively dissipated in order to maintain the condition that the plasma is magnetically dominated, as discussed in [34].

C. Boundary Conditions
This section is devoted to discuss the physical conditions at the global inner and outer boundaries of our domain and how to implement them numerically via the penalty method. Additionally, we also incorporate an approach introduced in [34] to restrict possible incoming violations of the divergence-free constraint, ∇ · B = 0 5 .
Our implementation of the penalties at the global boundaries is motivated by the interface treatment (1), identifying two main options: either one sets the incoming characteristic fields to a fixed source, regardless of the dynamics in the interior, or one may use the information leaving the system by setting the incoming physical modes to a particular combination of outgoing ones. The first choice has been already employed in [34] at the exterior surface, far away from the source, and we shall use it again here for the outer boundary. However, this approach is not suitable for setting the physical conditions at the stellar surface. Instead, at the inner boundary we adopt the second approach, with a very specific combination of outgoing modes, as we will describe next.

Inner boundary
We shall think the stellar surface as a perfectly conducting layer that separates the force-free exterior from the interior of the star, and let us denote by S its threedimensional world-volume. Our domain extends from the surface exterior S + , in the force-free regime. However, the boundary conditions will be induced by the perfectly conducting fluid at the surface interior S − .
Continuity across S of the Faraday equation, dF = 0, guarantees there are no jumps in the pullback to S of the electromagnetic tensor. Notice however that discontinuities on the remaining components of F may result from the presence of surface charges and/or current densities at the conducting layer. The jump conditions [A] ≡ A + − A − at the interface can be written in terms of the co-moving electric and magnetic field e a := F ab u b and b a := F * ab u b (i.e., where u a is the four-velocity of the co-rotating plasma), namely The definition of the fields at S − arise from interior structure of the magnetic field in the star and the perfect conductor condition for the electric field e a = 0. These values propagate for some of the components to the exterior solution through the junction conditions, namely where f (r, θ, φ) is an arbitrary function defined by the magnetic field structure in the interior of the star. In addition, as noticed in [48], the force-free condition e · b = 0 implies the normal component of the co-rotating electric field must vanish as well (i.e. n · e + = 0, there are no induced charges) whenever the magnetic field is not tangential to the stellar surface.
On the other hand, the ideal MHD condition e a = 0 also allows to easily translate the conditions applied to the co-moving EM fields to the fiducial ones, namely E i = − ijk v j B k . The velocity vector v i is constructed as the projection of the four-velocity u a = W (n a + v a ), being Let us assume then an axi-symmetric star, with a surface located at r = R (i.e., such that the normal vector is radial), that rotates with four-velocity u a = k a + Ω η a , being k a ≡ (∂ t ) a and η a ≡ (∂ φ ) a the two Killing vector fields of the Kerr spacetime. Thus, the resulting boundary conditions for the electric and the magnetic fields can be written as, where E r is to be enforced only if the magnetic field is not completely tangent to the star (i.e. B r = 0). Although written on a slightly different notation, the prescription for the tangential components of the electric field in (16) is equivalent to the given at [25]. Notice also that the remaining components of the magnetic field are free: the junction conditions involve unknown surface currents, so they can not be fixed.
To numerically implement these boundary conditions we proceed as follows. We will keep the normal magnetic field fixed to its interior value, B r = B r o (θ), by enforcing it at each Runge-Kutta substep as done in ref. [9]. The electric field components (16), on the other hand, are going to be prescribed by applying the penalty method to the incoming physical modes. That is, where L (a) (·) are a set of operators that must have a very precise structure in order to guarantee control of the semidiscrete energy by the penalties, namely: Here the co-basis elements Θ a > represents incoming modes, while Θ b ≤ the outgoing and zero modes. The idea is that one prescribes incoming physical modes at the boundary, combining information of the remaining characteristic modes. The coefficients R a b might be interpreted as reflexion coefficients and must be solved according to the particular condition one wishes to impose. A guidance in doing so, is provided by an observation regarding the action of the penalties on the boundary fields. Namely, that they will dynamically enforce: from which it is deduced that the operators L (a) (·) must project into the space of fields components to be specified (in our case, for example, the electric field components involved in (16)). While U o will act as a source, allowing to describe inhomogeneous boundary conditions. In Appendix A, we provide the reader with the specific boundary operators, L (a) (·), to enforce condition (16). We construct these operators, going from the electrovacuum case (i.e. Maxwell theory) to the full forcefree system (8)- (10). In between, we will discuss a simpler version of force-free electrodynamics often used, that might help to understand the transition to the full theory.

Outer boundary
The implementation of the outer boundary condition might be thought as a fictitious interface with an external field, U ext . That is, setting the penalties as: (20) which prescribe the incoming physical modes according to a fixed source we control. Thus, in ref. [34] for instance, U ext represented a uniform magnetic field threatening the BH magnetosphere, sourced by a distant accretion disk. As in the pulsar case there is no expected external electromagnetic sources, we will set U ν ext = 0, corresponding to maximally dissipative boundary conditions. This means no physical signals would enter through the outer surface and all waves reaching it will leave the domain without reflections. In practice, we will bring the exterior field smoothly to zero from its initial data configuration value. This way, we prevent the penalties from introducing spurious perturbations that would propagate inwards.

D. Initial Data
We describe the spacetime around the rotating neutron star by the Kerr metric, which is parametrized by the mass M and spin a. The metric in Kerr-Schild form can be written as, where η ab is the flat metric and a is a null co-vector with respect to the flat and the whole metric. In Cartesian coordinates (t, x, y, z) 6 , the metric function H takes the form, and the co-vector a reads a = 1, Assuming an homogeneous and spherically symmetric neutron star of radius ρ = R, we can fix the dimensionless moment of inertia I := I/M R 2 to the value 2/5 (see e.g. [25,26]). This choice allow us to relate the star angular velocity Ω with the spin parameter in units of the stellar radius, namely Notice that, for the spin range of realistic pulsars a/R 0.15, the Kerr spacetime is rather close to the metric of a neutron star in the slowly rotating limit, which is often used in the literature (e.g. [25][26][27]50]). Indeed, their metric components differ by less than 1% from ours in this range. There are two other important dimensionless quantities that characterize the problem. The first one is the surface rotation velocity, which ranges from 10 −4 to 10 −1 in realistic pulsars. Here R LC = Ω −1 denotes the usual light cylinder radius in flat spacetime. The second quantity is the star compactness, which coincides with the definition in Ruiz et al. [24] for their compaction, and differs in a factor 2 with other definitions like those in refs [25,27,50]. The compactness has a theoretical upper limit given by C < 4/9 [51], and for typical model of neutron star interiors its value is around C ∼ 0.2. We shall assume the magnetosphere is initially populated with a dipolar magnetic field given by the potential [52], where µ is the magnetic dipole moment. The electric and magnetic fields arising in our foliation of Schwarzschild spacetime are just: r + 1 whereas the other components vanish. Notice that this is an exact stationary solution of Maxwell vacuum equations for Schwarzschild but not for Kerr and that the magnetic field will only approximately satisfy the constraint ∇ · B = 0 on Kerr spacetimes. However, its value is dynamically kept very close to zero by using the divergence cleaning approach and constraint preserving boundary conditions. We will show evidence on such dynamic control of the constraint later on (see figure 3).
Notice also that in the Newtonian limit r 2M 1 (where α → 1 and β i → 0) the magnetic field reduces to the usual flat space result in an orthonormal tetradê λ , with vanishing electric field. In order to avoid transient initial currents and sharp profiles of the fields, the star will be smoothly brought from rest to its final angular velocity Ω by using a timedependent function, with t s usually chosen to be around 1/5 of the rotation period. To abandon axial symmetry we shall simply tilt the magnetic dipole axis of the initial configuration on the x − z plane, while keeping the rotation axis of the star (i.e., and the spacetime) unchanged on the zdirection. Tilting the dipole axis corresponds to the following change of Cartesian coordinates, being χ the inclination angle. Whereas for imposing the normal component of the magnetic field at the inner boundary, it will be necessary to keep track of the accumulated rotation phase during the evolution, and apply it to the initial field to get the boundary data for B r .

E. Analysis Quantities
We shall compute the spin-down luminosity as the integrated Poynting flux on spherical shells, where p a := −T ab k b is the conserved electromagnetic four-momentum (being T ab the electromagnetic stressenergy tensor and k a = (∂ t ) a the Killing vector field), p r ≡ p a (dr) a its radial component and d 2 w, the differentials of the angular coordinates (see e.g. [34,53]) 7 . We will consider the spin-down luminosity to be the one computed at the light cylinder, as in Petri [25]; while Ruiz et al. [24] has taken instead the asymptotic value, which might depend on the numerical resistivity and dissipation of the numerical scheme at the current sheet. The spin-down power of the aligned rotator has been estimated analytically in flat spacetime, yielding a dependence proportional to µ 2 Ω 4 (in geometrized units). This relation turned-out to be quite accurate, agreeing within a factor close to unity with most numerical simulations in flat spacetime (see e.g. [54] and references therein). Therefore, it is common to set 7 Notice this definition is consistent with the one in [25] written on an slightly different language, namely: L := S L (E ∧ H)r dS.
as a reference luminosity in order to rescale the results. The power of Ω appearing on (33) is directly related with another interesting quantity, known as the braking index n, defined asΩ = −K Ω n , or, equivalently, n = ΩΩ/Ω 2 .
Since the spin-down luminosity is associated to the star rotational energy loss by L = −I ΩΩ, one can formally derive n in terms of the luminosity and the surface rotation velocity v s , as done in ref. [21], The braking index measures how the spin period and its derivatives change with time, and it can be inferred from the physical mechanism regulating the spin-down.
For the aligned point-dipole magnetic field, this quantity is known to give exactly n = 3. However, considering neutron stars of finite size -as well as allowing curvature effects and misalignment-induces deviations from this value, as we shall see on Section III C.

III. NUMERICAL RESULTS
In this section we present a detailed numerical study of the pulsar magnetosphere, mainly focusing on how the general relativistic effects impact on the resulting spindown luminosity. First we study the aligned case, when the magnetic dipole coincides with the spin of the star and the solution is axially symmetric. Our code reproduces the well known properties of the aligned rotator magnetosphere on flat spacetime. Basic tests of our numerical implementation are included here, as well as a convergence analysis of the luminosity for different numerical resolutions and evidences of the correct behavior of the solenoidal constraint. Later, we incorporate the full general relativistic effects and compare our results with previous studies. We have carefully investigate the Poynting luminosity dependence on stellar surface velocity v s and compactness C (see (35) and Fig. 7); and then generalize it to the case where the magnetic dipole and rotation axes are not aligned, thus adding the dependence with the misalignment angle χ (see (36) and Fig. 10). Our main result is summarized in a fitted formula which is a good approximation for any star, depending only on three adimensional parameters: the surface velocity v s = R/R LC , the stellar compactness C = M/R and the misalignment angle χ between the spin and the magnetic dipolar moment. Finally, we have used this generic formula to estimate corrections to the braking index in terms of these three parameters.
A. Aligned Rotator

Magnetic field topology and tests
Our simulations reproduce all the well known features of an aligned rotator in flat spacetime (see e.g.  7,9,13,15]). The late time solutions exhibit a closed zone that co-rotates with the star, extending up to the light cylinder. This zone ends on a Y-point at the equatorial plane from where a strong equatorial current sheet begins. Our maximally dissipative boundary conditions allow us to place the outer numerical boundary fairly close to the central region, at around 4R LC . The evolution reaches steady state after ∼ 2 rotational periods of the star. A typical equilibrium solution is depicted in Fig. 1, showing the global magnetic field structure of the magnetosphere: lines represents the poloidal magnetic field and color illustrates the toroidal component. Notice that some of the field lines close outside the light cylinder due to effective resistivity of the current sheet, but the co-rotating region with vanishing toroidal magnetic field remains within the light cylinder.
The luminosity is constant between the stellar surface and the light cylinder, as expected from conservation of the electromagnetic energy-momentum tensor. Beyond the light cylinder the energy flux slowly decreases with radius due to dissipation at the equatorial current sheet. In Fig. 2 it is displayed the luminosity, integrated on spherical surfaces, as a function of the radius for a typical late-time solution of the aligned rotator. The curves corresponds to various numerical resolutions N θ ×N φ ×N r with N r = 2N φ = 4N θ . Clearly, the flux is already converging for a moderate resolution N θ = 40. We generally employ a higher resolution, with N θ = 80, for the results presented in this paper. We confirmed numerically that the resulting luminosities, normalized conve- niently with L o = µ 2 Ω 4 , do not depend on the given angular velocity Ω but only on the adimensional relation v s = R/R LC ≡ R Ω.
The behavior of the solenoidal constraint ∇ · B = 0 on a typical run is displayed in Fig. 3, where the quotient ||∇ · B|| 2 /||B|| 2 is plotted as a function of time for Kerr, Schwarzschild and flat spacetimes. As it can be seen, its value remains always smaller than 10 −3 , thus showing a decent control of the constraint during the evolution. Even on Kerr spacetime -where the initial data slightly violates the solenoidal constraint-we observe the divergence cleaning mechanism act to keep its value small, and moreover, it ends-up driving it to a similar value to that of Schwarzschild after an initial transient.

Comparison with previous results
For concreteness, we shall consider here full GR simulations at a high but still realistic star compactness C = 0.25. As previously observed in Ref. [25], we find that including curvature effects does not significantly change the field topology. It does, however, compress the poloidal magnetic field lines towards the neutron star equator, as it can be seen in Fig. 4. Gravitational effects have been also shown to considerably affect the polar cap structure (see e.g. [26,32,50]). As a result, there is an enhancement of the spin-down luminosity, as we shall see below. Figure [24] and Petri [25]. It can be seen that our numerical data is in very good agreement with these previous studies, despite using different formulation of force-free electrodynamics and numerical methods. In particular, we notice that in the region comprised between v s 0.1 and v s 0.3 the three models coincides almost perfectly. For this particular case with C = 0.25, we find an enhancement of the total emitted power due to relativistic effects of about 15% for v s = 0.1, 30% for v s = 0.2, and even larger for higher spin rates. It it worth mentioning that Petri reported a discrepancy with Ruiz et al. respect to the aforementioned increments. We do not observe such discrepancies and thus believe that it may be related with a misinterpretation of Ruiz numerical data, which used a rather different setting and notation. Poynting luminosities (normalized with Lo = µ 2 Ω 4 ) that include the effects of general relativity at a fixed compactness, C = 0.25. Our results (blue triangles/solid-lines) are compared with previous numerical studies at several rotation rates vs. Green asterisk/dottedlines represents the data taken from Figure 22 of Petri [25], while that the red dashed curve was constructed by fitting the values from Table II of Ruiz et al. [24] (evaluating at C = 0.25). Figure 6 compares the luminosities resulting from full GR simulations with those on the Newtonian regime. In the slow rotation limit, both GR and flat spacetime solutions seems to approach a similar value close to L o . This is consistent with the fact that relativistic corrections tends to disappear for small values of v s , as pointed out recently by ref. [26] for comparisons made at fixed µ in the regime R/R LC 1. For spin rates v s 0.05 deviations become apparent: Newtonian luminosities slightly decreases with rotation, while GR solutions increase their values significantly. Although one may distinguish among curvature effects associated with the spacetime mass and spin, they might be difficult to disentangle. We try to isolate the influence of stellar compactness from frame dragging effects by setting a = 0 (i.e. Schwarzschild metric). The results are shown in the top panel of Fig. 6, where it can be noticed that frame dragging only produces a modest enhancement for v s 0.1. It has been argued that the frame dragging effect is being compensated by an effective reduction of the angular velocity, which happens at the Kerr metric (see e.g. [32]) but not in Schwarzschild. Bottom panel of Fig. 6 presents the non-normalized luminosity as a function of the surface velocity on logarithmic scale, from where an effective braking index might be obtained. Our results compares well with those found by Petri [25], although we get slightly different values for the braking indexes: we estimated n = 2.94 for flat spacetime and n = 3.17 for GR (Fig. 6, bottom  and n = 3.12, respectively 8 .

Dependence on compactness and surface velocity
In order to analyze quantitatively the general dependence of the spin-down luminosity on gravitational effects, we have considered different stellar compactness at various rotation rates. By using our numerical results, we have been able to derive a fitting formula for the spin-down luminosity of an aligned rotator as a function of star compactness C = M/R and surface velocity v s = R/R LC . An excellent fit of the numerical data was 8 A plausible explanation for this small discrepancy is that Petri has considered smaller rotation rates (like vs = 0.01) where the trend is very close to n = 3. Thus, when fitting the curves he gets effective braking indexes closer to this value. Finally, notice that we could find stable solutions for star's compactness up to C ≈ 1/3. Beyond that limit, simulations become unstable near the stellar surface, probably due to the presence of a light ring.

B. Misaligned Rotator
We study now the pulsar magnetospheres for the misaligned cases, occurring when the dipolar magnetic field of the star is inclined an angle χ with respect to its rotation axis. Again, our simulations reproduce previous results in the literature of a misaligned rotator in flat spacetime. Our numerical evolutions settle to steady state after about two stellar rotation periods. As it can be seen in Fig. 8, the magnetic field topology on the µ − Ω plane is reminiscent of the aligned solution with closed and open zones. The current sheet starts at the intersection of the closed zone with the light cylinder and oscillates around the rotational equator. The Poynting fluxes are again constant between the stellar surface and the light cylinder, where they are measured. The energy dissipation at the current sheet decreases as the misalignment angle increases, as reported in reference [18] for the Newtonian case.
Including the effects of general relativity does not change dramatically this qualitative picture. In Figure  9, the equilibrium solution of the orthogonal rotator in GR is compared to the one obtained in the Newtonian regime. The equatorial field lines in both cases resemble the Deutsch solution inside the light cylinder and exhibit an spiral structure outside, being the GR (blue) lines more curved than the Newtonian (red) ones due to the gravity pull. Quantitatively, the presence of spacetime curvature enhances the spin-down luminosity in the misaligned case, although it depends softly on the compactness and the surface velocity. We have evolved the magnetospheres varying the inclination angles χ for three different rotation rates, comparing the luminosities between the Newtonian and the general relativistic regime of a star with compactness C = 0.25. The values of the luminosity obtained from our simulations were fitted by the following expression, L(C, v s , χ) L(C, v s , χ=0) = 1 + (1.24 + 8.14 v s C) sin 2 χ (36) Figure 10 summarizes these results. Numerical data for v s = {0.05, 0.1, 0.2} (at fixed compactness C = 0.25) is displayed together with their corresponding curves taken from (36). Notice for flat spacetime solutions, the normalized luminosity L/L(χ=0) does not depend on the surface velocity; and thus, there is just a single curve describing the zero compactness cases. Interestingly, the relation L ∝ k 1 +k 2 sin 2 χ proposed by Spitkovsky in [10] fits very well the numerical data, although the coefficients k 1 and k 2 now depend on the compactness and the spin rate. A careful comparison shows that our results are in good agreement with those of Spitkovsky for the particular case C = 0 and v s = 0.2; which were later reproduced by many other authors within different frameworks (see fig. 6 of [54] and references therein). Moreover, our re- sults are also consistent with [25], which was the first to incorporate GR effects to the misaligned rotator. Indeed, when fitting our solutions with L/L o = a+b sin 2 χ, we obtain the parameters listed on Table I for the cases v s = {0.1, 0.2}, which are only slightly larger than to those found by Petri (

C. Braking Index
One may rewrite the two fitting formulas (35)-(36) from last sections into a single general expression, where L o = µ 2 Ω 4 ≡ µ/R 2 2 v 4 s and, k(C, v s ) = 1.24 + 8.14 v s C From this approximate formula of the pulsar spin-down luminosity and equation (34), one can readily obtain an estimated braking index n. The calculation is rather straightforward and gives n = n(C, v s , χ) in a closed form. However, the resulting expression is not very enlightening by itself and so, instead of writing it down, we illustrate the results by plotting some representative curves in Fig. 11 at constant C and v s (top and bottom panels, respectively). These plots show a braking index ranging from n 2.8 for a rapidly rotating star on flat spacetime, up to n 4 for a rapidly rotating and highly compact NS. Actually, for realistic pulsars, C 0.1 and v s 0.2, the braking index increases both with the compactness and the surface velocity and it would be in the range n 3.0 − 3.5. Notice that, in the Newtonian regime, the braking index does not vary with the inclination angle χ, whereas in the relativistic case its value increases with misalignment. We also note that, as expected, deviations from the standard value n = 3 disappear for v s 0.1, which is the most common case for the observed pulsars. Unfortunately, the handful of pulsars with a reliable estimation of n have periods of the order of 0.1 − 1 s, for which the relativistic effects on the spin-down formula are negligible 10 .

IV. CONCLUSIONS
We have presented a formalism and a numerical code to carefully analyze the force-free magnetospheres of neutron stars, by extending previous studies on black hole magnetospheres. We have performed several tests to show the correct implementation of the boundary conditions at the stellar surface, which is the main difference with respect to the previous code. A careful and detailed decomposition of the eigenvectors of the evolution equations is required to apply suitable inner boundary conditions of our domain, which corresponds to a perfectly conducting neutron star.
We have studied, through three-dimensional timedependent numerical simulations, general-relativistic neutron star magnetosphere within the force-free approximation. Our results confirm other recent numerical investigations, in particular regarding the total electromagnetic power radiated by the neutron star. By performing suitable fits of our results, we provide a quite generic formula for the luminosity of a rotating dipolar magnetic field as a function of the compactness of the neutron star, its angular velocity, and the misalignment angle between the spin and the magnetic dipolar moment. From this formula, we have estimated deviations from the standard braking index value which will generally depend on these adimensional parameters. These deviations are rather modest, leading to n = 3.2 ± 0.2 for realistic millisecond pulsars. Unfortunately, current astrophysical observations are not accurate enough to distinguish such small differences on the braking index. Even for the pulsars displaying an almost constant braking index which could be measured with high precision, the observed values are systematically below the standard n = 3 [55,56]. As it has been confirmed in this work, these deviations can not be associated to relativistic effects. Therefore, we can only conclude that, if the measurements are correct, there must be other physical mechanisms modifying the luminosity of a pulsar magnetosphere which are not captured by the simple force-free model.

V. ACKNOWLEDGMENTS
We would like to thank Daniele Viganò for several helpful discussions throughout the realization of this work. F.C. and O.R. acknowledge financial support from Conicet, SeCyT-UNC and MinCyT-Argentina. F.C. and C.P. acknowledge support from the Spanish Ministry of Economy and Competitiveness grants AYA2016-80289-P and AYA2017-82089-ERC (AEI/FEDER, UE). C.P. also acknowledges support from the Spanish Ministry of Education and Science through a Ramon y Cajal grant. F.C. thanks Conicet financiamiento parcial de estadias breves en el exterior para becarios posdoctorales program and the Departament de Física & IAC3 for hosting his stay at the Universitat de les Illes Balears, where parts of this work were completed. This work used computational resources from Mendieta Cluster (CCAD, Universidad Nacional de Cordoba) and Pirayu Cluster (supported by the ASACTEI, Gobierno de la Provincia de Santa Fe, Proyecto AC-00010-18), which are both part of the SNCAD -MinCyT-Argentina.

Appendix A: Characteristic Decomposition and Boundary Operators
In this appendix, we build appropriate boundary operators L (a) (·) (see eqn. (18)) to deal with a perfectly conducting surface via penalty terms. The construction strongly relies on the characteristic structure of the particular evolution system under consideration. Thus, before going to the main result concerning the full forcefree system used in this paper, we shall first illustrate the method by applying it to vacuum electrodynamics and to a simpler version of the force-free system. The first one is the usual linear theory with transversal propagation modes, and serves as example for the whole construction. While the second system, already nonlinear, posses a much simpler characteristic structure that helps to better understand the transition to the full theory. This simpler approach to force-free electrodynamics has been widely used in the literature, as in e.g. [53], and so it might be also interesting in its own right. These systems can be described by the same covariant equations, where an extra dynamical field, φ, has been added to handle the magnetic divergence-free constraint (e.g. [43][44][45]). If I a = 0 one recovers usual electrodynamics in vacuum. On the other hand, the presence of the plasma can be captured by an effective force-free current density, which includes the drift current only. In this approach the force-free condition, E · B = 0, needs to be enforced separately: either by damping strategy like the one used on [47], or by cutting electric field at each Runge-Kutta substep as in e.g. [44,53]. Notice, ρ : represents here the electric charge density. The resulting evolution equations are nonlinear and, as it can be shown, they constitute a strongly hyperbolic system.
Let us introduce some important notational conventions for the remaining of this section. We shall denote the set of dynamical fields U µ as, (where the index i is used to denote spatial vectors) and its "dual" element Θ µ , Characteristic decompositions will be taken respect to the wave front propagation direction, given by a unit vector m i that will be identified in this context with the boundary surface normal. Some useful abbreviations are defined: for any vector A i . Notice A i p and A i q are orthogonal to each other, and both tangent to the boundary surface.

Vacuum electrodynamics
The physical modes in Maxwell electrodynamics are transversal to the propagation direction m i , and propagate at the speed of light with eigenvalues λ ± = β m ± α. The associated characteristic basis (and co-basis) elements can be written as, being (e i p , e i q ) any two orthogonal transversal directions. There are in this case three extra modes associated with constraints, playing no role in the construction of the boundary operators. The operators L (a) (·) will only involve here the physical modes, and shall be written as: The idea is now to find the coefficients, R ij , as to remove from both L i the action on any free (i.e. non-restricted) component of the boundary fields. That means, in our case, to remove the components of Θ + 1 and Θ + 2 acting on the tangential magnetic field. An immediate consequence of such requirement is that, The action on the difference dU µ = U µ o − U µ then reads, This way, we are able to control both tangential components of the electric field at the boundary, i.e.: Notice the boundary is acting as a mirror, reflecting back part of the outgoing waves while keeping fixed the tangential components of the electric field.

Alternative force-free system
We look now at the evolution system which arises from a standard (3+1)-decomposition of equations (A1)-(A3). In this system, the transversal subspace (with eigenvalues λ ± = β m ± α) appear slightly modified and there is an extra physical mode 11 , which might be incoming or outgoing according to the sign of λ E := β m +ασ E (where σ E = Sm B 2 ). The complete eigen-system is here, Outgoing case, λ E ≤ 0. Here there are two incoming physical modes, namely: U + 1 and U + 2 . And thus, the 11 The presence of the plasma turns the characteristic mode associated to the constraint ∇ · E = 0 in vacuum electrodynamics, into a new physical propagation linked to electric charge density.
general boundary operators 12 reads: The idea is, as before, to find the coefficients which removes from L i the action on any free component of the boundary fields. That means, to remove the components of Θ + 1 and Θ + 2 acting on the tangential magnetic field. As a consequence, while it seems there is some freedom on the other coefficients R 13 and R 23 . One possible choice being, which results on no explicit enforcing of the normal electric field. The action would then read, Another natural election would be to set R 13 = 0 and R 23 = 0. Obtaining, We shall use (A11)-(A12) for the general case in which the magnetic field is not tangential to the stellar surface, so that the normal component of the electric field is enforced. While for the special case where B m = 0, we shall adopt (A9)-(A10) instead, fixing solely the tangential components.
Incoming case, λ E > 0. Now there are three incoming physical modes, namely: U + 1 , U + 2 and U 3 . The general boundary operators are: In the present situation, the restriction to the operators completely determine all the coefficients: from which one obtains, We observe here, as opposed to the case λ E ≤ 0, is not possible to enforce the tangential electric field alone. But recall this was required only for the case in which B m = 0, and it can be shown that λ E = 0 when the magnetic field is purely tangential. Hence, there are no inconsistencies in the scheme and everything fits together perfectly. It seems the extra condition for E m , arising from the forcefree constraint, is relevant to the characteristic structure of the theory and plays a crucial role in the construction of the boundary operators.
3. Full force-free system a. Characteristic Structure The fully nonlinear force-free system (8)-(10) has two transversal modes (λ ± = β m ± α) and two Alfvén modes ). These are modes associated with physical propagation. The remaining modes correspond to constraints: two of them are related with divergence cleaning propagations (of eigenvalues λ ± 0 = β m ±α) and the third one represents the algebraic force-free constraint with λ Ψ = β m + α Sm B 2 . We write here the full eigen-system: with its associated co-basis, where we have defined: , ∆ := B 2 − E 2 and the normalization, N ± A := (σ ± A ) 2 − 1 .

b. Eigenvalues Analysis
To go further into the construction of the boundary operators it seems necessary to better understand the Alfvén eigenvalues, since they depend on the values that the background fields take at each point during the evolution. Hence, it is not obvious which of these modes will turn-out incoming or outgoing across the stellar surface.
Assuming the perfect conducting condition (16) at the boundary, we were able to perform a very general (but rather straightforward) calculation which gives: where v i := β i + Ω η i . Under very mild assumptions like e.g. ΩR 0.5, a 0.5 and C 0.8, it is possible to infer: This means that there will be always an incoming and an outgoing Alfvén mode at the stellar boundary, except in the special case when the magnetic field is tangential to the surface (i.e. B m = 0). Indeed, in such situations one has λ ψ = λ + A = λ − A = 0, which relates with the well known fact that Alfvén waves cannot propagate orthogonal to the magnetic field.

c. Boundary Operators
According to the previous analysis, we would always have one incoming and one outgoing Alfvén modes (i.e. λ ± A > 0 and λ ∓ A < 0); except in the special case where B m = 0, where one has λ + A = λ − A = 0. On the other hand, we also have incoming and outgoing transversal modes with eigenvalues λ ± = β m ±α. So there are always (provided B m = 0) two incoming physical modes that can be prescribed via penalties, and an extra unphysical Ψ-mode which may be incoming or outgoing 13 .
There are two incoming modes: U + T and U ± A ("±" depending on the signs of λ ± A ). The general boundary operators are: The idea is to apply the operators to a generic field 14 ,Û , and find the coefficients that cancel the terms containing 13 The remaining modes, linked to dynamical enforcement of the constraint ∇ · B = 0, will not be treated on the same footing. This differential constraint is handed with a different method to restrict possible incoming violations from the boundary [34]. 14 It is very important to keep these fields separate from those background fields present in the co-basis elements Θ i . φ andB i . For convenience, we shall considerΨ ≡ Θ Ψ (Û ) instead ofB m on the calculations, by means of the relation: E mBm = ∆ 2Ψ − B mÊm − (E p ·B) − (B p ·Ê).
We shall solve first for L T . After some redefinitions of the coefficients, we obtain the following system: x Ψ which can be solved to get, and finally, Analogously, we solve now for the operator L ± A .
and we obtain the following system: x Ψ which can be solved to get, Special case, B m = 0. In this spacial situation there is just one incoming mode, namely: U + T . The general boundary operator is: Recall we only can prescribe the tangential electric field here, so we need to further remove theÊ m component in the construction. Notice it will be now possible, due to the presence of an extra Alfvén mode in the right hand side. We shall define again the coefficients, but now we obtain instead the system: which leads to x Ψ = x 0 = 0 and, The operator then takes the form: and we see it only penalizes tangential components of the electric field at the boundary.
Incoming case, λ Ψ > 0. Now there are three incoming modes, namely: U + T , U ± A and U Ψ . The general boundary operators are: It is important to notice that there is not enough terms at the right hand sides of the first two expressions above to remove all the components we need; in particular,Ψ. However, we have an extra incoming (constraint) mode that will allow us to enforceΨ = 0 through the penalty method. That allow us to solve the system, and we get: There are only three cases to consider separately from the general one we have described so far. The first two happen when one of the Alfvén subspace degenerates with one of the magnetosonic (i.e. λ A = λ ± ), which is possible when: E m = 0, E i p ⊥ B i p , E 2 p = B 2 p and ∆ ≡ |B m |. The third situation arises when this happens for both Alfvén modes at the same time, only possible if E i = 0 and B i p = 0.
The relevant co-basis elements are, The boundary operators are fairly easy to solve; the only non-vanishing coefficients are R T T = −1 and R AA = 1. It leads to, where the corresponding incoming directions are, The relevant co-basis elements are, The boundary operators are fairly easy to solve; the only non-vanishing coefficients are R T T = −1 and R AA = 1. It leads again to, but now the Alfvén direction in the penalty is given instead by, with λ + A = β m + α(B 2 m − B 2 p )/B 2 , being the positive eigenvalue. Some remarks follows: • The special situation in which B m = 0 does not need any treatment, since is not allowed under these degenerate conditions.
E i = 0 and B i p = 0: In such case, the characteristic system reduces to vacuum electrodynamics. And thus, we adopt the boundary operators discussed in Sec. A 1.