General-relativistic monopole magnetosphere of neutron stars: a pseudo-spectral discontinuous Galerkin approach

The close vicinity of neutron stars remains poorly constrained by observations. Although plenty of data are available for the peculiar class of pulsars we are still unable to deduce the underlying plasma distribution in their magnetosphere. In the present paper, we try to unravel the magnetospheric structure starting from basic physics principles and reasonable assumptions about the magnetosphere. Beginning with the monopole force-free case, we compute accurate general-relativistic solutions for the electromagnetic field around a slowly rotating magnetized neutron star. Moreover, here we address this problem by including the important effect of plasma screening. This is achieved by solving the time-dependent Maxwell equations in a curved space-time following the 3+1~formalism. We improved our previous numerical code based on pseudo-spectral methods in order to allow for possible discontinuities in the solution. Our algorithm based on a multi-domain decomposition of the simulation box belongs to the discontinuous Galerkin finite element methods. We performed several sets of simulations to look for the general-relativistic force-free monopole and split monopole solutions. Results show that our code is extremely powerful in handling extended domains of hundredth of light-cylinder radii~$\rlight$. The code has been validated against known exact analytical monopole solutions in flat space-time. We also present semi-analytical calculations for the general-relativistic vacuum monopole.


INTRODUCTION
It is well admitted that pulsars are strongly magnetized and rotating neutron stars surrounded by electron-positron pairs filling their magnetosphere. However an accurate description of the interaction between this plasma and the neutron star electromagnetic field remains poorly constrained by observations. Moreover a realistic model should also include some radiative processes. We are still far from a comprehensive and self-consistent picture of the pulsar machinery. Both plasma flows and strong gravity impact on the structure of the magnetosphere. Curvature and frame-dragging effects are indeed important due to the high compactness of neutron stars. For typical models of neutron star interiors, the compactness is about Ξ = Rs/R ≈ 0.5 where Rs = 2 G M/c 2 is the Schwarzschild radius, M is the mass of the neutron star, R its radius, G the gravitational constant and c the speed of light.
It is the purpose of this paper to study the response of the electromagnetic field to the combined effect of plasma ⋆ E-mail: jerome.petri@astro.unistra.fr screening and curved space-time in the vicinity of a neutron star. To this aim we compute general-relativistic solutions in the force-free approximation. The problem we therefore address is similar to the electrodynamics of black hole magnetospheres. Actually, the numerical technique employed are the same expected that in our case we do not have any complication arising from the presence of an event horizon. Komissarov (2004b) was the first to report on numerical simulations of general-relativistic monopole magnetospheres of black holes in the magnetohydrodynamic regime. Komissarov (2004a) also investigated the properties of the magnetospheric plasma in the force-free limit. Since then, several authors followed the effort of modelling generalrelativistic magnetospheres of compact objects. McKinney (2006a) designed a general-relativistic code for force-free magnetospheres and McKinney (2006b) applied it also to neutron stars. Pétri (2013) showed that multipole vacuum solutions in general relativity can be computed semi-analytically via the 3+1 formalism through a vector spherical harmonic expansion method introduced by Pétri (2012). However it is well known that a neutron star cannot be surrounded by vacuum. Indeed, the electric field induced by the rotation of the magnetic field generates huge Lorentz forces able to extract particles from the crust and therefore filling the magnetosphere. For simplicity, as a first step towards more realistic magnetospheres, the force-free assumption is often quoted. In that case the plasma dynamics is completely dominated by the electromagnetic field which is a good approximation for neutron star magnetospheres. The resulting force-free geometry has been investigated by many authors like for instance in the aligned case by Contopoulos et al. (1999); Parfrey et al. (2012) and the general oblique rotator by Spitkovsky (2006); Kalapotharakos & Contopoulos (2009); Kalapotharakos et al. (2012); Pétri (2012). Generalrelativistic force-free neutron star magnetospheres have been less investigated so far. But Beskin (1990) already mentioned that general-relativistic effects can significantly distort the parallel component (with respect to the magnetic field) of the electric field. This can have important implications for particle creation, acceleration and radiation in the polar caps. Indeed, deviation from the corotation charge density leads to a parallel component of the electric field determined by the magnetic field geometry. Therefore, as also claimed by Muslimov & Tsygan (1992), space-time curvature and frame dragging effects are important for the electrodynamics of the gaps. Several numerical techniques have been applied to model such magnetospheres. Usually the schemes are closely related to the finite volume algorithm, a well tested method for computational fluid dynamics due to its conservative properties. High resolution shock capturing techniques enable an increase of the spatial order of the method but at the expense of larger stencils. Such scheme are also useful to solve Maxwell equations. Recently, another arbitrary high order method, the discontinuous Galerkin approach, has been tested in general relativity by Radice & Rezzolla (2011).
Our goal in this paper is to quantify precisely the distortion induced by general-relativistic effects, namely curvature of space-time and frame dragging. To this end, we solve the time-dependent Maxwell equations in curved space-time in spherical coordinates. Nevertheless, as a starting point we restrict the solutions to the monopole field in order to elucidate the consequences of general relativity avoiding complications induced for instance by the presence of a cusp at the light-cylinder in the case of an aligned dipolar magnetic field. Strictly speaking, at this Y-point, the magnetic field strength vanishes and can lead to problems in the force-free approximation due to the electric current density prescription. Nevertheless, in order to show how efficiently the code can handle discontinuities such as current sheets for instance in the equatorial plane, we present Newtonian as well as general-relativistic simulations of the split monopole field. Consequently, we use the 3+1 formalism of electrodynamics as briefly reminded in Section 2. Next we give approximate solutions to the vacuum monopole field in Section 3 which will be useful for benchmarking the code whose algorithm is described in Section 4 and then tested in flat space-time in Section 5. Application of our new code to vacuum and forcefree curved space-time monopoles are presented in Section 6. We extend our study to the split monopole case to demonstrate the ease of handling discontinuities. Conclusions and ongoing work are drawn in Section 7.

THE 3+1 FORMALISM
In this section we briefly remind the set of Maxwell equations in curved space-time following the 3+1 formalism for a fixed background metric. We split space-time into an absolute space {x, y, z} and a universal time t, similar to our all day experience. Advantages of such a split have been demonstrated in many numerical simulations about neutron stars and black hole magnetospheres.

The split of the space-time metric
The four dimensional space-time is split into a 3+1 foliation such that the background metric g ik can be expressed as (1) where x i = (c t, x a ), t is the time coordinate or universal time and x a some associated space coordinates. The Landau-Lifschitz convention is used for the metric signature given by (+, −, −, −) (Landau & Lifchitz 1989). α is the lapse function, β a the shift vector and γ ab the spatial metric of absolute space. By convention, latin letters from a to h are used for the components of vectors in absolute space, in the range {1, 2, 3}, whereas latin letters starting from i are used for four dimensional vectors and tensors, in the range {0, 1, 2, 3}. A fiducial observer (FIDO) is defined by its 4-velocity n i such that This vector is orthogonal to the hyper-surface of constant time coordinate Σt. Its proper time τ is measured according to For a slowly rotating neutron star, the lapse function is and the shift vector We use a spherical coordinate system (r, ϑ, ϕ) and an orthonormal spatial basis (er, e ϑ , eϕ). The metric of a slowly rotating neutron star remains close to the usual flat space, except for the radial direction. Indeed the components of the spatial metric are given in Boyer-Lindquist coordinates by For this slow rotation approximation, the spatial metric does not depend on the spin frequency of the massive body but only on M through α. The spin a is related to the angular momentum J by J = M a c. It follows that a has units of a length and should satisfy a Rs/2. Introducing the moment of inertia I, we also have J = I Ω , Ω being the spin frequency. In the special case of a homogeneous and uniform neutron star interior with spherical symmetry, the moment of inertia reads Thus the spin parameter can be expressed as For the remainder of this paper, we will use this expression to relate the spin parameter intervening in the metric to the spin frequency of the neutron star. From the above expression, note that the parameter a/Rs remains smaller than 0.4 because R = 2 Rs and rL 2 R in our set of simulations.

Maxwell equations
Maxwell equations in absolute space take a form very similar to their traditional expression in Newtonian space except that space is curved. The time-dependent Maxwell equations in a prescribed metric (possibly time-dependent) read The source terms (ρ, J ) will be specified by the force-free condition, see paragraph below. The three dimensional vector fields are not independent, they are related by two important constitutive relations, namely ε0 is the vacuum permittivity and µ0 the vacuum permeability. The curvature of absolute space is taken into account by the lapse function factor α in the first term on the righthand side and the frame dragging effect is included in the second term, the cross-product between the shift vector β and the fields. The derivation of the above equations is given in Komissarov (2004a). From the auxiliary fields (E , H ) we get the Poynting flux through a sphere of radius r by computing the two dimensional integral on this sphere by where dΩ is the infinitesimal solid angle and Ω the full sky angle of 4 π sr.

Force-free conditions
The source terms have not yet been specified. They are deduced from the force-free condition that in the 3+1 formalism become which implies E · B = 0 and therefore also D · B = 0. As in the special relativistic case, the current density is found to be, see the derivation for instance in Komissarov (2011) B and D/ε0 can be interpreted as the magnetic and electric field respectively as measured by the FIDO. Moreover its electric current density j is given by Maxwell equations (9a)-(9d), the constitutive relations in equations (10a) and (10b) and the prescription for the source terms in equation (13) set the background system to be solved for any prescribed metric in the force-free approximation.

VACUUM MONOPOLE FIELD
Before dealing with the force-free solution, we recall the exact vacuum electromagnetic field in flat space-time and extend the result to the general-relativistic monopole field, valid up to first order in the spin parameter of the star. Although the monopole assumption is not realistic for the zeroth order magnetic field of the neutron star, it gives us insight into the effects of curved space-time on to force-free magnetospheres. Such solutions will also serve as benchmark for testing and checking current and forthcoming electromagnetic codes in general relativity.

Newtonian solution
We start with a simple monopole magnetic field anchored in a perfectly conducting star of radius R and rotating at a speed Ω around an axis passing through its centre. Let us denote this axis by ez. The strength of the magnetic field at the surface is B. Thus, in Minkowski space-time, the exterior vacuum solution for a rotating magnetic monopole is given by assuming that the electric field in the comoving frame vanishes in the interior of the star. The induced electric field is therefore of dipolar nature. Note that the relation between E and D is simply ε0 E = D. In terms of the "potential", see equation (19) below, we can write it as This means that the only non-vanishing coefficient is All other coefficients of the expansion like f E l,m and g E l,m should be equal to zero. Remember that a divergencelessness vector field E can be expanded according to where Φ l,m are vector spherical harmonics, see for instance Pétri (2013). Later, we will use this expression to check the numerical accuracy of our code, see section 5.

General-relativistic solution
In order to look for the analytical solution to the general-relativistic monopole field in vacuum, we use the formalism developed in depth by Pétri (2013). Closed analytical expressions have only been found for the first order expansion of the electric field D as described in the first part of this section. For higher order approximations, we have to resort to numerical solutions which are exposed in the second part of this section.

First order expansion
The background monopolar magnetic field in eq. (15a) remains exact for the curved space-time geometry. We are looking for a first order approximation to the electric field such that which automatically satisfies ∇ · D = 0. To zeroth order, the magnetic field is not perturbed, we leave it unchanged. For the electric field, the function f D 1,0 satisfies A particular solution vanishing at infinity is The general solution of the homogeneous equation also vanishing at infinity is For the general and particular solutions we have respectively For convenience, we introduce the following constants The index R means that quantities are evaluated on the neutron star surface. Then the constant of integration in eq. (22) reads To summarize, to first order in spin parameter, the electric field satisfies We will use this analytical expressions to check our code in the general-relativistic case to the lowest order in the spin parameter expansion.
In the limit of a weak gravitational field, the solution reduces to equation (17) as it should. To the next leading order in a, we expect a dipolar perturbation of the magnetic field, thus we write The function f B 2,0 will be a solution of Taking into account the boundary conditions, f B 2,0 has to vanish at infinity and at the neutron star surface. This corresponds to the Deutsch approach where the radiative disturbances of the normal component of B are not taken into account. So far we have not found any analytical expression to solve this boundary value problem. We have to resort to numerical integration. This is explained in the next paragraph.

Multipole expansion
The most general situation including multipoles to any order is exposed in this paragraph. We look for solutions that can be expanded in the following series Each of the coefficient f D l,0 and f B l,0 has to satisfy the differential equation which is given for the electric and magnetic field respectively by ∂r(α 2 ∂r(r f D l,0 )) − The Kronecker symbol δ l,1 appearing in the differential equation for the electric field represents the contribution from the monopole magnetic field, that cannot be expressed in terms of a curl. We add it explicitly. Let us write down these equations for the three first coefficients in B and D. The system of partial differential equations then reads The associated boundary conditions are where quantities have to be evaluated on the neutron star surface, at r = R. Details on the derivation of these equations can be found in Pétri (2013). We emphasize that the magnetic field at the neutron star surface is exactly matched to the expression for the general-relativistic monopole, equation (15a). All other multipole fields f B l,0 with l 1 vanish at r = R by our definition.

Numerical solution
The above system of boundary value problems is efficiently solved by means of rational Chebyshev polynomials. The technique is presented in detail in Pétri (2013). Here we only report the results for the coefficients f D l,0 and f B l,0 for the monopole. For concreteness, in all the computations, we use the following set of parameters namely R/Rs = {2, 2000} and rL/R = {10, 1000} which should correspond to a compact and a non compact star as well as to a mildly rotating and a slowly rotating star.
First we only consider the dipolar electric field component induced by the rotation of the neutron star. Strictly speaking, we should retrieve the analytical approximation equation (27). This is indeed what we checked. In figure 1 we show on the left panel the absolute value of these expansion coefficients f D 1,0 on a logarithmic scale and on the right panel the relative error. We consider two sets, the first one computed from the analytical exact expression and the second one obtained from the numerical integration of the boundary value problem. The agreement between both solutions is excellent, the error being less than 10 −15 which correspond to the double precision arithmetic of ε = 10 −16 . The coefficients decrease exponentially fast demonstrating the rapid convergence of the series to the exact solution. This exponential convergence to the exact solution is typical for spectral methods when the solution is C ∞ . The relative error increases systematically when the coefficients become of the order ε f D 1,0 . These weak coefficients cannot be computed accurately because of the finite precision of the computer. This is of no concern as in any expansion series, they become irrelevant because not contributing to the summation in a significant way.
After this first test of the solution to the boundary value problem, we switch to the next order of approximation including a perturbation in the magnetic field which will be of quadrupolar order. We thus have to solve simultaneously for f D 1,0 and f B 2,0 . In order to show the rapid convergence of the coefficients, we plot again their absolute values in logarithmic scale, as depicted in figure 2. For the next approximation, we add the multipolar coefficients f D 3,0 and f B 4,0 . Convergence is proven by inspection of figure 3 showing an exponential decay of the coefficients with respect to the index k. Finally, for the most accurate solution we put two other multipolar components, namely f D 5,0 and f B 6,0 . Figure 4 compares the relative importance of each multipolar component with respect to each other. We always observe the characteristic exponential convergence as expected in this smooth boundary value problem. All the coefficients of the electric and magnetic field functions decrease 5,0 , f B 6,0 } are present. Adding higher mulitpole components to the expansion series will not drastically change the lowest order coefficients f D 1,0 and f B 2,0 , at most only starting from the fifth digit. Indeed, for a given couple {R/Rs, rL/R}, all curves almost overlap whatever the number of multipoles added into the expansion. Multipolar fields higher than l = 4, although present are definitely too weak to have an influence on the electric dipole and magnetic quadrupole fields.
To conclude this section, we plot the radial dependence of the functions ,0 } in the four cases corresponding to a slowly or rapidly rotating star, compact or not, with parameters R/Rs = {2, 2000} and rL/R = {10, 1000}, see figure 6. These functions can then directly be compared to the output of our numerical simulations in section 6. The functions are normalized in order to put them on a same graph except for f D 1,0 which is the leading term. The non compact object case with R/Rs = 2000 remains very close to the flat space-time solution. Thus a good approximation to the electric field is given by equation (17). This is clearly seen in the upper left panel of figure 6. In principle, we are able to compute the electromagnetic field to any order to get the solution to any required precision. Actually we stopped with a three terms expansion in the electric and magnetic field respectively because high order multipole moments become negligible compared to the lowest order. We also think that it is largely enough to compare with the numerical code we now describe in details.
Their absolute values are shown on a logarithmic scale. The inset legend shows the couple of ratios {R/Rs, r L /R}. Note the different scales used for each plot.

CODE DESCRIPTION
We now give the general outline of our pseudo-spectral discontinuous Galerkin finite element algorithm. The main ingredients are, the expansion on to vector spherical harmonics for divergencelessness fields in spherical shells, an exact imposition of boundary conditions on the neutron star surface, an explicit time stepping with a fourth-order Runge-Kutta integration scheme, a spectral filtering in the longitudinal and latitudinal directions and a limiting procedure in the radial direction. The radial part is solved with a high-order finite volume scheme whereas the spherical part is solved through a pseudo-spectral approach.

One dimensional scalar conservation law
To present our new code, we will focus on the one dimensional scalar conservation law which is an archetypal of equations often used to model physical phenomena. Consider therefore the simple conservation law of a scalar field denoted by u with a physical flux function denoted by f such that the conservation of u is expressed as a partial differential equation written as This equation has to be solved for any time t 0 and for all x ∈ [a, b] where [a, b] is the computational domain. Note that in our code x should be interpreted as the radial coordinate r. We subdivide the domain [a, b] in K cells not necessarily of the same length. In each of these cells which we denote by D k with k ∈ [0..K − 1], the solution is expanded on to a basis of spatial functions φ k i such that the approximate solution in the cell k reads valid in the cell k given by the interval [x k l , x k r ]. The basis possesses Np + 1 functions. The spatial method is therefore of order Np. After injecting this expansion into the conservation law equation (34) and projecting on to the basis functions φ k i , performing two successive integrations by part in each cell independently, starting from we arrive at the strong form of the partial differential equation such that We introduced a numerical flux f * which tells to the system how to communicate information between adjacent cells as in classical finite volume schemes. Taking for instance These matrices can be computed analytically and exactly. The semi-discrete system to be solved then becomes or in a pure matrix notation with P k the column vector of the normalized Legendre polynomials. Note that for an orthogonal basis, the mass matrix is diagonal hence very easy to invert. Inverting the mass matrix M k , each coefficient of the expansion evolves according to the first order ordinary differential equation The state of the art in the discontinuous Galerkin methods resides in the choice of the numerical flux f * which has to satisfy several stability and consistence properties. The reader is referred to the excellent book by Hesthaven & Warburton (2008) for a detailed discussion about the implementation of modal and nodal discontinuous Galerkin methods in one dimension and the tricks to deal with non-linear problems, introducing limiting and filtering processes. Here we only give guide lines on the way to implement the techniques for spherical geometries. Let us first discuss the main advantage of the code, namely the flexibility in the choice of the grid.

The grid
Our goal is to look deeply into the light-cylinder with very small ratios of neutron star radius to light-cylinder radius, R/rL ≪ 1, as well as far away from the light-cylinder at distances r much larger than rL, r/rL ≫ 1. In our previous work Pétri (2012), we had some difficulties to achieve such demanding parameters because we used only one radial domain to expand on to Chebyshev polynomials. We thought that the code could greatly benefit from a more advantageous domain decomposition in the radial direction. Indeed, this allows us to zoom into the light-cylinder with very fine grids close to the surface but keeping a coarser grid outside the light-cylinder where we can afford a loss in precision for sufficiently large radii. Due to the flexibility of domain decomposition methods, we are able to use a non-uniform grid when moving from one radial cell to the next one. This technique is called spectral element method (Canuto et al. 2007). It can be seen as a high-order finite volume scheme. To use all the advantages of the conservative form of such finite volume formulation, we prefer to expand the radial direction into normalized Legendre polynomials instead of Chebyshev polynomials. Such expansion makes the algorithm rigorously conservative, meaning that the average value of the unknown quantities are perfectly conserved during the simulation, within numerical round-off errors.
The arbitrary nature of the radial scale is used to fix small volumes close to the neutron star whereas larger shells are sufficient farther away. To be more specific, we employ the usual Fourier transform in the {ϑ, ϕ} directions and expand the radial coordinate into K sub-intervals (which can be seen as finite volume elements), the boundary of each cell is given by [r k g , r k d ] with k ∈ [0..K − 1] dividing the global interval [R1, R2] into non necessarily equal sub-intervals. In each of these volumes, we expand the radial part into normalized Legendre polynomials by rescaling each interval [r k g , r k d ] into [−1, 1] through a scaling function. Let us assume that the computational domain is comprised between the neutron star surface at R1 = R and an arbitrary outer radius R2. The spherical shell is decomposed into K cells but with increasing thickness. We introduce two temporary variables y1 = log(R1/rL) and y2 = log(R2/rL) and a logarithmic thickness by h = (y2 − y1)/K. Each cell, labelled with a superscript k, possesses then two interfaces located at The thickness of the cell labelled k is h k = r k d − r k g . In that way, the ratio between the size of two successive cells is constant and equal to e h . We will show that such variable cell size drastically improves the accuracy in the innermost parts of the simulation box while preserving good accuracy well outside the light-cylinder.

Vector expansion and divergencelessness constraint on B
We use again a clever expansion of the vector fields B and D. Indeed, electric and magnetic fields are expanded onto vector spherical harmonics (VSH) according to Such expansion is done in each cell. However, in order to deal with the divergencelessness of the magnetic field whatever the configuration of the electromagnetic field, loaded or not with plasma it is more appropriate to use an expansion of B into where {f B lm (r, t), g B lm (r, t)} are the expansion coefficients of B. The monopole part, eq. (15a), is added by hand. To impose the divergencelessness constraint, we project the magnetic field on to the subspace subtended by the expansion in equation (44). Actually, because spectral methods for smooth problems are very accurate, the projection is not required at each time step. We perform it only when the divergence becomes larger than a threshold defined by the user.

Numerical flux
As in any other finite volume scheme, communication between cells goes through a numerical flux f * chosen to resolve as accurately as possibly the conservation laws. In the force-free limit, the dynamics reduce to the solution of Maxwell equations with source terms. So we only need to find an appropriate numerical flux for the linear advection problem in one dimension, namely the radial direction. The efficiency of the numerical code will strongly depend on the choice of the numerical flux. For Maxwell equations, we employ a first order upwind scheme as described in Hesthaven & Warburton (2008). Starting from the 3+1 formalism, we consider the one dimensional system of Maxwell equations in spherical geometry and relevant for propagation in the radial direction. Thus only the components (E ϑ , E ϕ , H ϑ , H ϕ ) are meaningful. In this way we get the following equations describing the propagation of the electromagnetic field in the radial direction in general relativity by By a change of variables through the quantity the above system becomes strictly conservative, assuming that the lapse function is time-independent. We can then apply standard discontinuous Galerkin methods to our problem. Introducing the jumps of the electromagnetic field components at the cell interface, denoted by du = u d − ug, the associated numerical upwind flux becomes From these expressions, we deduce the right hand side on the left interface of a cell by and the corresponding right hand side on the right interface of a cell by These numerical fluxes close the overall description of the basic algorithm. We now switch to the delicate problem of non-linearities and how to overcome aliasing effects and related numerical instabilities.

Slope Limiter
The slope limiting technique is adapted from the classical finite volume community. The idea is to reduce or even kill spurious oscillations that arise from the non-linear evolution or from sharp discontinuities in the solution. The most basic total variation diminishing (TVD) limiters are usually too dissipative for higher-order schemes. Toro (2009) detailed several TVD schemes with application to simple problems and compares the merit of each slope limiter. We refer the reader to this book for more information about the use of TVD method in finite volume algorithms. Indeed, in trouble cells, the polynomial expansion is reduced to at most a linear interpolation and therefore considerably reducing the order of the method around discontinuities. To circumvent such drawbacks, it is necessary to release the TVD property for a less stringent property called total variation bound (TVB) method (Cockburn et al. 1989). The latter does not guaranty strict cancellation of oscillations but only weaken them whereas the former completely avoids oscillations but at the cost of reducing to a low-order scheme. In the simulations shown in this paper, we found that the TVB limiter represents a good compromise between accuracy and spurious oscillations. We implemented both limiters and checked that TVB is preferable to TVD limiters. Unfortunately TVB methods introduce one more parameter, often depicted by the capital letter M . Moreover the value of this parameter is very problem dependent, related to the second spatial derivative of the solution, therefore a priori unknown. So we let the user arbitrarily choose the best limiter parameter M by some trial and error tests. Various examples of limiters can be found in the literature, see Hesthaven & Warburton (2008) for some basic discussion, including the difference between TVD and TVB. In our algorithm, we tried the MUSCL limiter and the less dissipative TVBM limiter. For high enough resolution we did not find any significant difference between both limiters. Thus we will not discuss the influence of these limiters on the solution.

Filtering
The limiter cannot be applied in the latitudinal and longitudinal direction simply because there is no domain decomposition in those directions. We use the classical spherical harmonic expansion. The force-free problem being non-linear due to the electric current in the source terms, we expect the solution to develop sharp gradients or discontinuities also in the spherical directions. It is therefore compulsory to get rid of these high frequencies by some filtering procedure. This is achieved by adding a small damping factor to the high order coefficients of the expansion in Y lm . Filtering is performed at each time step. We use an exponential filter in directions (ϑ, ϕ) given by the general expression where the variable η ranges between 0 and 1. For instance, in the latitudinal direction η = l/(N ϑ − 1) for l ∈ [0..N ϑ − 1], l being the index of the coefficient c l,m in the spherical harmonic expansion f (ϑ, ϕ) = N ϑ −1,Nϕ−1 l,m=0 c l,m Y l,m (ϑ, ϕ) and N ϑ , Nϕ the number of collocation points in the spherical direction (latitude and longitude). The parameter α (not to be confused with the lapse function) is adjusted to values not too large in order to avoid errors in the solution but also not too small in order to sufficiently damp these oscillations.
The above mentioned exponential filter of order β does not strictly satisfy the condition for the smoothing factors as explained in Canuto et al. (2006). However, for numerical purposes we choose α such that e −α is numerically zero i.e. below the machine accuracy ε. In practice, we choose α = 36 assuming double precision computation with ǫ ≈ 10 −15 . The order β of the smoothing influences the dissipation rate in the solution. The low order multipole components are weakly damped and correspond to large scale structures. If the solution shows fine scale structures, the filtering has to be minimized. We will discuss the role of β in the particular case of the split monopole solution in the next sections. We typically tried β ∈ {2, 4, 8}. Actually, because higher order multipoles are almost absent in the solutions, let it be vacuum or force-free, a low order filtering was enough to reach satisfactory accuracy. In all the simulations presented in this work, if not explicitly specified, we systematically used a fourth order filter with β = 4. We also tried a second and eighth order filter without significant variation in the solution. The split monopole is a notable exception for which higher order filtering and a large number of collocation points are necessary to correctly catch the discontinuity induced by the equatorial current sheet.

Exact boundary conditions
As in Pétri (2014) we put exact boundary conditions on the star. In general relativity the correct jump conditions at the stellar surface, continuity of the normal component of the magnetic field Br and continuity of the tangential component of the electric field {Dθ, Dφ} are such that Br(t, R, ϑ, ϕ) = Br 0 (t, ϑ, ϕ) (51a) Dφ(t, R, ϑ, ϕ) = 0 (51c) The continuity of Br automatically implies the correct boundary treatment of the electric field. Br 0 (t, ϑ, ϕ) represents the, possibly time-dependent, radial magnetic field imposed by the star, let it be monopole, split monopole, oblique dipole or multipole. The outer boundary condition cannot be handled exactly. We need to make some approximate assumptions about the outgoing waves we want to enforce in order to prevent reflections from this artificial outer boundary. Using the Characteristic Compatibility Method (CCM) described in Canuto et al. (2007) and neglecting the frame-dragging effect far from the neutron star, the radially propagating characteristics are given to good accuracy by Dθ ± ε0 c Bφ ; Dφ ± ε0 c Bθ In order to forbid ingoing wave we ensure that whereas the other two characteristics are found by the subscript PDE denoting the values of the electromagnetic field obtained by straightforward time advancing without care of any boundary condition. The new corrected values are deduced from the solution of the linear system made of equations (53a)-(54b).

Time integration
One of the strength of pseudo-spectral methods is that they replace a set of partial differential equations (PDE) by a larger set of ordinary differential equations (ODE) for the unknown collocation points or spectral coefficients. Schematically, it can be written as with appropriate initial and boundary conditions. u represents the vector of unknown functions either evaluated at the collocation points or the spectral coefficients. We use a fourth-order Runge-Kutta scheme advancing the unknown functions u in time. See also the discussion in Hesthaven & Warburton (2008) for more details about other time integration schemes especially those called strong stability preserving Runge-Kutta methods including the popular schemes of order two and three (SSPRK2,3).

Initial conditions
The rotation of the neutron star is switched on smoothly as in Pétri (2012). Its spin frequency increase slowly from zero in order to avoid the formation of sharp gradients. This would be especially true at time t = 0 where there is no electric field outside the star but right on its surface. Taking an evolution of the spin frequency as Ω(t) = sin 2 t 8 for t 4 π 1 for t 4 π (56) therefore starting at a null value avoids the initial discontinuity in the electric field. The spin frequency as well as its first derivative are smooth at the initial time of the simulation t = 0. No gradient or sharp features are expected. We next switch to a discussion of the results.

NON RELATIVISTIC TESTS
For the remaining of the paper, we adopt the following normalization: the magnetic moment of the star is equal to unity, therefore B R 2 = 1, as well as the stellar angular velocity and the speed of light, Ω = c = ε0 = µ0 = 1, therefore the light-cylinder radius is rL = 1. We start with a discussion about the non-relativistic monopole solutions in vacuum but also in the force-free limit. Interestingly analytical closed expressions do exist in these cases. They are very valuable solutions to check the correctness and accuracy of our code. The generalrelativistic rotator will be treated in section 6.
Although a discontinuous Galerkin method is intended to do better than second order in space, in this paper we only show results with Np = 1 i.e. use linear polynomial interpolation of the unknown fields. Indeed, so far we only implemented TVD and TVBM limiters which fall down to first order at shocks or when a limiting procedure is applied. We plan to add higher order slope limiters in the near future such as the moment limiter described in Biswas et al. (1994). Fortunately we already get accurate solution with linear polynomials.

Vacuum monopole solution
We tested our code against some well known analytical solutions. The starting point is the vacuum monopole field for which the Poynting flux is equal to zero. The solution has been presented in section 3. The analytical solution is exact and easy to compare with the output of our simulations.
We start our computation with a non rotating monopole magnetic field, Ω = 0, and zero electric field outside the star, except for the crust where we enforce the inner boundary condition, see equation (51). Note however that due to our special profile of Ω(t), the electric field at the surface of the star is initially equal to zero. It will slowly increase to its maximal value reached at a normalized time t = 4 π.
We performed simulations with different spin frequencies of the neutron star corresponding to several ratio between stellar radius R and light cylinder radius rL such that rL/R = {2, 10} and between the artificial outer boundary and the light-cylinder Rout/rL = {10, 100, 1000}. Obviously the resolution of the grid should be highest for the largest domain in radius with r/rL ∈ [0.1, 1000]. A minimum resolution of K × Np × N ϑ = 128 × 1 × 4 was necessary. Actually, throughout the paper, we will show results with a higher resolution of K × Np × N ϑ = 256 × 1 × 8. Because of the axisymmetry of the problem, a Fourier transform in the azimuthal direction is not necessary, so we simply put Nϕ = 1. We let the system evolve until it reaches a stationary state inside the simulation box. Thus the final time strongly depends on the location of the outer boundary, it can be as high as t final = 300 π for Rout = 1000 rL.
In the non-relativistic monopole solution, the magnetic field remains unchanged. The only relevant quantity to check is the coefficient f D 1,0 (r) for the electric field. It is understood that all other coefficients should be equal to zero. In figure 7 we show this coefficient f D 1,0 (r) on the left panel and its relative error on the right panel for several sets of parameters. Note that it is plotted on a log-log scale in order to make more visible the outer part of the function. A careful investigation of this outer part shows a slight deviation of the computed solution with respect to the analytical solution. Let us assume that the solution is accurate if the relative error is less than the one reached close to the neutron star surface. Then if Rout = 10 rL the computed solution becomes inaccurate above ≈ 5 rL but if Rout = 100 rL then the discrepancy starts at ≈ 50 rL and finally for Rout = 1000 rL the inaccuracy starts at ≈ 500 rL. This behaviour clearly indicates an influence of the location of the outer boundary on the numerical solution. Such artifact can only be removed by moving away the artificial outer boundary. Using a smaller time step will not help to improve the accuracy or to remove the outer boundary influence. Indeed, we run the same sim-  ulations with a time step 2.5 or 5 times smaller than the one presented here for relative error. We have not noticed any changes in this error so these plots are not shown to avoid congesting the figures. The corresponding Poynting flux is shown in figure 8. As expected it is very close to zero as it should be. The accuracy is better than 10 −3 in the whole simulation box whatever its size. Note that even if the solution is inaccurate at large distances, the associated Poynting flux, although having large errors, remains close to zero. This is explained by the fact that the electromagnetic field in those region is weak. It is impossible to compute the relative error in the Poynting flux because the exact value should be zero.
To conclude with the vacuum case, note that Maxwell equations become linear. Therefore we do not need to apply a strong limiting in the radial direction. In that case, we can use higher order spatial expansions of the unknown fields without destroying the high order of the method. This has been done for instance with a quadratic Np = 2 and a fourth order Np = 4 polynomial expansion. Results of such simulations are shown in figure 9 where the relative error in the function f D 1,0 is plotted and has to be compared with the corresponding plot in fig. 7 with Np = 1. We used the same number of cells in each computation. It is clear that higher order methods are much more accurate. This demonstrates the need for limiters that do preserve the high order accuracy of discontinuous Galerkin schemes.
The above results demonstrate that the code is able to catch accurate solutions of the vacuum electromagnetic field with appropriate boundary conditions on the perfectly conducting star and at large distances. As we now discuss, in the force-free limit the code also gives accurate solutions.

Force-free monopole solution
Next we tackle the problem of an axisymmetric force-free flow known as the monopole field introduced by Michel (1973). We recall that this monopole solution is given by In terms of a vector spherical harmonic (VSH) expansion, this magnetic field is expressed as where g B(exact) 1,0 all other coefficients being equal to zero. The associated Poynting flux is The initial set up is the same as in the previous paragraph. We only add a source term represented by the force-free current given by equation (13). During the evolution of the electromagnetic field, it is easy to show that the component Br remains constant in time and that only the Bϕ component is present with the coefficient g B(exact) 1,0 (r). The numerical value of this coefficient is shown in the left panel of figure 10. Moreover, in order to prove the accuracy of our code, we plot the ratio g B 1,0 /g B(exact) 1,0 (r) and compare it to unity as depicted in figure 10, right panel. The accuracy is better than 6 digits in the whole computational domain. For completeness we also plot the Poynting flux obtained from the simulations as shown in figure 11. From the analytical solution, we known that the Poynting flux is a constant, irrespective of the size of the neutron star. This is indeed what we found. In normalized units, the Poynting flux is equal to unity whatever the ratio rL/R and whatever the location of the outer boundary. The result is very accurate, better than 7 significant digits. Interestingly, contrary to the vacuum monopole field, the force-free solution does not suffer from the location of the outer boundary. We always found the exact analytical expression (to high numerical accuracy) in the whole simulation box. Thus a small digression about these outer boundaries is in order as exposed in the next paragraph. Flat force free monopole Figure 11. The Poynting flux relative error for the force-free monopole solution for Rout/r L = {10, 100, 1000} and a ratio r L /R = {2, 10}. As expected it is equal to Lmono to very high precision, better than 7 digits.

Influence of the location of the outer boundary
Imposing exact outgoing wave boundary conditions on a sphere of finite radius is a tedious work. Indeed Novak & Bonazzola (2004) showed that the Sommerfeld radiation condition is only valid for the monopole field. For dipolar or even multipolar structures, restricting the infinite domain to a sphere of radius Rout will lead to some deviation from a perfect outgoing wave. To elucidate the influence of the location of this outer sphere, we looked at the error of the Poynting flux with respect to the location of Rout defined by where Lana and Lnum are the analytical and numerical Poynting fluxes respectively. We report our results in this brief paragraph for the flat space-time, choosing a radius of the neutron star equal to rL/R = 2 and Rout/rL = {10, 100, 100}. For the vacuum or force-free field we know exact solutions. As we already showed in figure 7 there is a slight influence for the vacuum field. Nevertheless we did not found any influence on the force-free solution.
We demonstrated in this section that our pseudospectral discontinuous Galerkin code is mature and able to compute accurately vacuum as well as force-free electromagnetic fields in flat space-time. Boundary conditions have been implemented in an efficient way avoiding spurious reflections and artificial inner boundaries as usually required for finite difference/volume methods. Before looking at the general-relativistic solution we finish the test in flat spacetime by a discussion about the important situation where a current sheet is present in the solution.

Split monopole solution
Our first intention to implement the discontinuous Galerkin method was to handle multi-domain computational boxes, allowing for non-uniform grids and therefore larger scales. However, this method is also well suited for the study of solutions presenting discontinuities. So we decided to test our code against a magnetic field structure showing a current sheet in the equatorial plane as for instance in the split monopole field. It is well known that analytically the solution is made of two half monopole fields of opposite "magnetic charge" separating the space into two hemispheres where the above force-free monopole applies separately. We have not met any particular problem to deal with this discontinuous solution. Let us investigate in more details the split monopole.
At the surface of the star, the radial component of the magnetic field reverses polarity at the equator. It therefore represents a step function in the ϑ variable on which we perform a series expansion. This jump will introduce the well-known Gibbs phenomenon and decrease the convergence rate to the worst case: first order. The Gibbs phenomenon produces an associated overshoot in Br that do not decrease by increasing the number of terms in the expansion, i.e. N ϑ . This is proved rigorously mathematically. The filtering explained in the code description section will help to enforce a lowering of these spurious oscillations. In any case, the current sheet does not pollute or even destroy the solution in the simulation domain. The Poynting flux is shown in fig. 12 for the ratio Rout/rL = {10, 100, 1000} and rL/R = {2, 10}. Theoretically, we know that this flux should be equal to the forcefree monopole luminosity, so in normalized units it should equal to unity. But the filtering and limiting procedures, useful to prevent strong numerical oscillations and possible non-linear instabilities, introduce some nonphysical dissipation. This is clearly recognized in fig. 12 where the computed Poynting flux decreases with radius. The rate of dissipation can be controlled by the resolution of the simulation and the filtering. This is shown in fig. 13 where the azimuthal component Bϕ is plotted against the colatitude ϑ at three different radii, namely at the neutron star surface, at some point inside the simulation box and at the outer boundary. We recognize the Gibbs phenomenon through its oscillatory nature in the vicinity of the discontinuity. The solution becomes more accurate when we increase the number of coefficients in the ϑ expansion and/or if we reduce the influence of the filtering on the lowest multipole coefficients.
The dissipation outside the light-cylinder is close to 25%. We plan to reduce this strong dissipation by replacing the fist order TVBM limiter by higher order filtering and increasing the number of discretization points in both directions. Nevertheless, this improvement of our code is left for future work.

GENERAL-RELATIVISTIC MONOPOLE SOLUTIONS
We now present new results about the monopole force-free solution in general relativity. We adopt the fixed background metric for a slowly rotating neutron star in Boyer-Lindquist coordinates as described in section 2. The same spin frequencies than those for the nonrelativistic solutions are used, corresponding to Rout/rL = {10, 100, 1000} whereas the spin frequency is such that Flat force free split monopole Figure 13. Azimuthal component of the magnetic field Bϕ for the split monopole solution at three different radii, at the neutron star surface r = R (red curve), at some point inside the simulation box r = r L (blue curve) and at the outer boundary r = 10 r L (green curve), using different filtering orders, to the left, β = 4 and to the right β = 8. The parameters are N ϑ = 32, Rout/r L = 10 and r L /R = 10. Bϕ is multiplied by r to get ride of the radial dependence. In the exact analytical solution, all three curves should overlap.

Vacuum monopole
Approximate expressions for the vacuum monopole field in general relativity are given as outlined in section 3. No outgoing electromagnetic wave propagating into vacuum space exists except for a transient regime relaxing to the stationary state. The Poynting flux as seen by an observer at infinity therefore vanishes. We checked this assertion by plotting the Poynting flux in figure 14 according to equation (11). Different runs are shown corresponding to increasing size of the simulation box, namely for the set of ratio rL/R = {10, 100, 1000}. The Poynting flux vanishes everywhere to very good accuracy. Moreover the electromagnetic field evolved to a steady state without reflection at the outer boundary. Our characteristics compatibility method used in flat space-time does also give good results in a curved spacetime, when the outer boundary is kept far from the light cylinder, justifying its numerical use. Note however that the outer edge of the box is not rigorously transparent to electromagnetic waves as was already the case with the flat vacuum monopole. The only remedy to this inaccuracy is to enlarge the box size at the expense of computational time due to the propagation delay between the star and the outer boundary and due to the requirement of higher grid resolutions. In order to give an estimate of the accuracy of our computed solution, the first order approximation of f D 1,0 given by the analytical expression equation (27) is compared with the output of the pseudo-spectral discontinuous Galerkin code. Results are shown in figure 15 for the function f D 1,0 itself, on the left panel, and its relative error on the right panel. We find good agreement between both functions. Although the time-dependent simulations contain multipolar electromagnetic fields with l > 1, the computed solution do not differ much from the analytical expression containing only the dipolar electric field l = 1. As expected, the corrections induced by the mulitpolar components remain negligible. General-relativistic effects stay on a low level. The  (11) and Lmono is given by equation (60). The computed flux vanishes as expected, within the numerical precision of the algorithm. The solution settled down to a stationary state. The inset legend corresponds to the ratio Rout/r L = {10, 100, 1000}. Note the logarithm scale in radius.
flat vacuum function almost overlaps the curved space-time counterpart.
In the general-relativistic case too, a higher order spatial expansion remains more accurate than a low order one. To demonstrate it, we performed here again simulations with Np = 2 or Np = 4. Results are shown in figure 16 and should be compared to the linear approximation in fig. 15 with Np = 1.
To sum up, we proved that our code is able to computed accurate solutions of the electromagnetic field in a fixed curved geometry. We eventually switch to the most interesting case, the general-relativistic force-free monopole field.

Force-free monopole
Current wisdom assumes that pulsars are neutron stars surrounded by relativistic plasmas of electron/positron pairs. If the pulsed radio emission comes from the magnetic poles, then we should look for accurate configurations of the magnetic field in the vicinity of the neutron star including curved space-time and plasma screening effects. This last paragraph is intended to bring us one step closer to this difficult task. As a starting point, we envisage a monopolar magnetic field instead of the more traditional dipolar structure. A detailed investigation of the general-relativistic dipole force-free magnetosphere is left for upcoming work. A simple prescription including the plasma current is based on the force-free approximation as explained in section 2. The simulation set up is the same as in vacuum except that the force-free current is switched on. We summarize the results by showing the Poynting flux for the different runs as presented in figure 17. We found that the power radiated does not significantly deviate from its Minkowski version. For the case rL/R = 10, the normalized Poynting flux is equal to unity with 0.1%. It is the same as the Michel monopole solution given in Michel (1973). However, for the case rL/R = 2, we observe a deviation from the flat space-time monopole Poynting flux around 2%. This decrease of the luminosity is a direct consequence of the frame dragging effect, being more pronounced in that case. If we artificially switch off the frame dragging effect by setting β = 0, we would retrieve to good accuracy the Newtonian case, independently of the ratio rL/R. For completeness the coefficient g B 1,0 is also compared to its flat spacetime version in figure 18. General relativity distorts the field sensitively close to the neutron star. As expected, generalrelativistic effects are important only close to the neutron star surface where curvature and frame-dragging are significant. The lowest order azimuthal magnetic field geometry is distorted with respect to its flat counterpart. Nevertheless the Poynting flux as measured by a distant observer is not significantly affected by such perturbations. Even if multipolar components are present, they remain at a low level compared to the dominant multipole electric and magnetic field. Curved force free monopole Figure 18. Deviation of the general-relativistic magnetic field coefficient g B 1,0 from its Newtonian version in the force-free monopole with Rout/r L = {10, 100, 1000} and a ratio r L /R = {2, 10}. It is compared to the flat monopole through the relative error g B 1,0 /g B(f lat) 1,0 − 1.

Split monopole
The study of the split monopole in flat space-time can be repeated in general relativity. A typical set of runs is shown in fig. 19. The same remarks as for the Newtonian split monopole hold here. Dissipation is again introduced in order to minimize the effect of the Gibbs phenomenon. At the stellar surface the Poynting flux is maximal and close to the true value but as soon as we depart from the stellar surface, energy is dissipated and diminishes the measured outgoing Poynting flux. Energy is dissipated up to 30%. Nevertheless the solution settled down to a stationary state. We also show the azimuthal magnetic field component Bϕ at three different radii, namely at the neutron star surface, at some point inside the simulation box and at the outer boundary, fig. 20. The Gibbs phenomenon is apparent through its oscillations in the vicinity of the discontinuity. We gain accuracy by increasing the number of coefficients in the ϑ expansion and/or by increasing the order of the filtering. The discontinuity is better resolved by keeping higher orders but at the expense of introducing stronger oscillations. Compare the right panel (β = 8) of fig. 20 to its left panel (β = 4).

Performances
As a last point, we show some performances of our code by providing the time needed for the code to produce the presented solutions for different resolutions and time steps.
As an example, we looked at the time spend for computing the force-free monopole solution in curved space-time for the parameters rL/R = 2 and Rout/rL = 10. Computations have been done on a single core processor with clock around 2.2 GHz. Results are summarized in table 1. For reasonable accuracy we need around one hour and for high precision a few days a needed. Note that for larger simulation boxes, such as the one presented in this paper, the computation time has to be multiplied by 10 for Rout = 100 rL or 100 for Rout = 1000 rL. Another factor 5 is required if rL/R = 10. Such parameter space is at the edge of our current computational capability. We plane to write a new version employing the Message Passing Interface (MPI) in the near future to improve these performances and most importantly to be able to compute accurate three dimensional simulation of an oblique pulsar magnetosphere including the Fourier transform in the azimuthal direction.

CONCLUSION
General-relativistic force-free pulsar magnetospheres are the simplest approach to a self-consistent accurate investigation of the electromagnetic field configuration and plasma distribution around compact objects. In this paper, in order to quantify the effects of a curved background metric,  fig. 13, r = R for the red curve, r = r L for the blue curve and r = 10 r L for the green curve.