Diffusivity in force-free simulations of global magnetospheres Free

ABSTRACT

1 INTRODUCTION

Energy flows induced into magnetically dominated relativistic magnetospheres of compact objects are commonly modelled by numerical simulations in the force-free electrodynamics (FFE) limit. Fueled by the track record of observations in the era of multimessenger astrophysics, current targets for such simulations include the magnetospheres of rapidly spinning black holes, spiraling neutron stars, magnetars, and pulsars. The tenuous, magnetically dominated atmosphere (magnetosphere) of pulsars is an active field of scientific interest. They fascinate both observers (e.g. Lorimer et al. 1995; Ransom et al. 2005; Abdo et al. 2013; Jankowski et al. 2018) and theorists (e.g. Kennel & Coroniti 1984; Lyubarskii 1996; Contopoulos, Kazanas & Fendt 1999; Goodwin et al. 2004; Timokhin 2006; Timokhin & Arons 2013; Contopoulos 2019; Pétri 2020). With the remarkable progress in scientific computing, their rotating magnetosphere has captured designers of numerical methods that integrate FFE and magnetohydrodynamics (MHD) with ever-improving accuracy (e.g. Komissarov 2006; Spitkovsky 2006; Tchekhovskoy, Spitkovsky & Li 2013; Parfrey, Spitkovsky & Beloborodov 2017; Carrasco & Shibata 2020). Recently, particle-in-cell (PIC) simulations were able to resolve a broad range of scale separations and allow for unprecedented insight into the microphysics of pulsar magnetospheres across the global scale (Cerutti et al. 2015; Philippov, Spitkovsky & Cerutti 2015a; Kalapotharakos et al. 2018; Philippov & Spitkovsky 2018; Guépin, Cerutti & Kotera 2020). In this fascinating flurry of outcomes, only few references scrutinized whether the results from ideal plasma simulations are the best possible model for the pulsar magnetosphere that contains an inherently non-ideal region, namely the equatorial current sheet (ECS) beyond the closed zone (Contopoulos 2016; Contopoulos 2019; Contopoulos, Pétri & Stefanou 2020). Here, we study with rigorous technical depth how this non-ideal region can affect the global dynamics of the force-free aligned rotator magnetosphere – effectively serving as a blueprint for force-free magnetospheres of other compact objects.

Ideal FFE evolve Maxwell’s equations for the electromagnetic fields |${\boldsymbol E}$| and |${\boldsymbol B}$| while rigorously maintaining the force-free conditions |${\boldsymbol E}\cdot {\boldsymbol B}=0$| (equivalent to the no Ohmic heating condition |${\boldsymbol E}\cdot {\boldsymbol j}=0$|⁠, where |${\boldsymbol j}$| is the electric current) and |${\boldsymbol E}^2-{\boldsymbol B}^2\lt 0$| (magnetic dominance). Non-ideal force-free fields are those fields that allow for perturbations to the condition |${\boldsymbol E}\cdot {\boldsymbol B}=0$| by non-negligible electric fields E_∥ along the magnetic field |${\boldsymbol B}$| (e.g. Lyutikov 2003). Such fields were used in the literature to induce the concept of resistivity into FFE, either by specifically designed driving currents (e.g. Alic et al. 2012) or by self-consistently modelled alterations to the current of ideal FFE (Komissarov 2004; Parfrey et al. 2017; Li et al. 2019).

In numerical models of FFE, it is common to enforce the preservation of the |${\boldsymbol E}\cdot {\boldsymbol B}=0$| and |${\boldsymbol E}^2-{\boldsymbol B}^2\lt 0$| conditions algebraically by resetting the electric field instantaneously wherever they are not fulfilled (e.g. Palenzuela et al. 2010; Mahlmann, Levinson & Aloy 2020). Exploiting the specifics of our numerical methodology (Mahlmann et al. 2021a), which combines the instantaneous algebraic extraction of all non-ideal electric fields while evolving the charge continuity equation, we identified another mechanism that adds diffusivity to an ideal FFE scheme (Mahlmann et al. 2021b). Namely, the misalignment of electric fields |${\boldsymbol E}$| and charge density ρ can significantly alter an FFE evolution. In this numerical survey, we further develop the findings obtained with an idealized setup that triggers tearing modes in force-free current sheets (Mahlmann et al. 2021b) to the astrophysical relevant scenario of pulsar magnetospheres.

In the not uncontroversial realm of magnetospheric simulations in the force-free limit, we announced in Mahlmann et al. (2021a) that ambiguities in the standard reference of the force-free aligned rotator required further attention. In fact, Contopoulos (2016) pointed out these ambiguities (trade secrets) that arise when simulating magnetospheres with non-ideal regions, such as the pulsar and Wald magnetospheres, in the ideal plasma limit, as one finds in ideal MHD and FFE. A sensitivity for the dazzling amount of calibration that time-dependent numerical simulations require is rarely transmitted along with the visually appealing results themselves. This manuscript is an effort to bring transparency to the modelling of one astrophysical scenario that crosses the constraints set by FFE. It aims at enabling the reader to ask crucial questions when evaluating results from the simulations of magnetospheres and intends to place some landmarks for the development of future hybrid methods that are not restricted to the ideal regime.

This manuscript is organized as follows. In Section 2, we review the employed numerical methodology as well as the pulsar magnetosphere initial data (Section 2.2). Section 3 presents the outcome of the conducted simulations of a force-free aligned rotator. The results are grouped by different topics. First, we examine the dependency of the luminosity at the light cylinder (LC) on the method employed to preserve the consistency between the charge distribution and the currents (Section 3.1). Section 3.2 analyses the conservation of energy beyond the LC. An array of ancillary high-resolution models yields additional insights into the subtleties of ideal FFE simulations in Section 4. Specifically, we assess the role of the magnetic dominance condition (Section 4.1), compare algebraic corrections of force-free violations to driving currents (Section 4.1.1), and study the effect of resistivity models beyond the LC (Section 4.1.2). The discussion of Section 5 includes views on the propagation of force-free violations (Section 5.1), the diffusive time-scales set by the employed hyperbolic/parabolic cleaning (Section 5.2), and a general picture of diffusivity in force-free magnetospheres (Section 5.4). We conclude this survey by summarizing the main takeaways of the presented results in Section 6.

2 METHODOLOGY

The aligned rotator problem has been studied vastly throughout the last 20 yr and now appears to be a well-established test case for FFE codes. Even if it is likely a simplification of the more complex problem of a magnetic dipole that is misaligned with respect to the rotational axis of the pulsar, it still contains an element of special relevance for the overall problem, namely, the ECS. We approach this investigation by means of time-dependent FFE simulations performed with the numerical code presented in Mahlmann et al. (2021a; 2021b). Our method is an enhanced high-order conservative realization of FFE as introduced by Komissarov (2004) and vastly benefits from the Carpet driver (Goodale et al. 2003; Schnetter, Hawley & Hawke 2004) and its extension to spherical coordinates (Mewes et al. 2018, 2020) supported by the infrastructure of the Einstein Toolkit.¹

The scheme that we introduced in Mahlmann et al. (2021a) has since been compared to another force-free MHD method (bhac, Ripperda et al. 2019a,b) in the context of Alfvén-wave interactions in the highly magnetized limit (Ripperda et al. 2021). It achieved a striking convergence of results across different methods. Also in the context of the aligned force-free rotator (cf. Section 5.2 in Mahlmann et al. 2021a), our FFE method reproduced the main characteristics of the pulsar magnetosphere: co-rotating magnetic field lines, a force-free closed zone in the wake of the Y-point, and an ECS. However, Mahlmann et al. (2021a) observed a shift of the Y-point away from the LC and a Poynting flux that was above the value which is treated as an established reference throughout the literature (e.g. Spitkovsky 2006; Tchekhovskoy et al. 2013; Etienne et al. 2017). In the following subsections, we briefly introduce the governing equations for the problem at hand, and the methods used to integrate them. We also provide the detailed numerical setup in Section 2.1.

2.1 Numerical method

In a variation of the default numerical scheme, we ignore the charge conservation equation (5) and compute the charge density appearing in the source terms by imposing |$\rho =\nabla \cdot {\boldsymbol E}$|⁠. Specifically, the divergence of the electric field is computed from the cell-centred (volume averaged) electric field values in equations (2) and (6). We maintain the hyperbolic equation for the scalar potential Φ (4) in order to dissipate and transport away any misalignment between charges and currents and preserve the charge conservation equation (5) up to truncation error. We note that the |$\nabla \cdot {\boldsymbol E}$| term appearing in equation (4) is computed as a numerical flux for the temporal update of Φ and, as such, it is obtained from the inter-cell values of the electric field. These interface values are monotonically reconstructed from cell-centred volume-averaged values of the electric field. Hereafter, we will refer to this variation of the base scheme as the local charge reconstruction (LCR) method. We note that we employ a finite volume method where none of the variables is staggered off the cell centers (differently from, e.g. Spitkovsky 2006; Mignone et al. 2019). Thus, the evaluation of |$\nabla \cdot {\boldsymbol E}$| at each numerical cell is performed by employing a fourth-order finite difference approximation based on a suitable number of neighboring cell values of |${\boldsymbol E}$| (with an stencil similar to that used for the evaluation of |${\boldsymbol j}_\parallel$|⁠, Mahlmann et al. 2020, cf. Section 3.4).

Finally, we also considered a second variation of the default numerical scheme that combines the two previous strategies into a hybrid method (HCC hereafter). In the HCC algorithm, we restrict the use of the LCR method to places where the magnetic dominance condition is violated during any sub-step of the time-integration, and the CC method is applied elsewhere. In practice, for the specific context of the aligned rotator magnetosphere, this limits the application of the LCR method to numerical cells affected by the ECS.

2.2 Simulation setup

This rather complicated boundary is needed to prevent spurious numerical behaviors at the stellar surface, which happen if all the magnetic field components are held fixed at r_in ≤ r ≤ r_*. Likewise, this boundary acts equivalently to a conducting boundary in as much as it preserves the continuity of the electric field components parallel to the surface and the magnetic field perpendicular to it (i.e. no jumps in B^r, E^θ or E^ϕ develop at the stellar surface). It allows to relax the initial values of the electromagnetic field at the stellar surface to the approximate equilibrium values attained after a few rotational periods. The electric charge inside the pulsar boundary layer is calculated in every sub-step of the time-integration via |$\rho =\nabla \cdot {\boldsymbol E}$|⁠. To ensure consistency at the inner boundary, we set Φ = 0 for r < r_*, while evolving Ψ freely to evacuate errors to the solenoidal constraint on the magnetic field |${\boldsymbol B}$|⁠.

3 FORCE-FREE ALIGNED ROTATOR

In this section, we evaluate and interpret results from an extensive array of simulations (Table 1). Specifically, Sections 3.2 and 3.1 present results from 41 simulations in the ideal FFE limit, spanning different resolutions as well as a parameter range for |$\kappa _\Phi$| and |$c_\Phi$|⁠. Fig. 1 marks the starting point – and preceding motivation – of our exploration: The FF equilibrium aligned rotator magnetosphere is vastly different when the LCR method is employed (top left-hand panels), as opposed to a conservative evolution of the charge continuity equation in our CC scheme.

Figure 1.

Magnetic field component B^ϕ (normalized to B₀) and parallel force-free current |${\boldsymbol j}_\parallel /j_0$| for different values of the diffusivity parameter α ∈ [0.015, 0.145, 1.447, 4.630, 9.260, 18.52, 36.17] in the CC and in the LCR methods for a resolution of N_r = 64. In all cases, we employ |$c_\Phi =1$|⁠, and for the run with the LCR method, we choose the extreme value α = 36.17. Black fieldlines are seeded at the same latitude on the stellar surface and may serve as a reference points for comparability. The LC position is indicated by a solid black line, and an estimate of the Y-point location by a dashed black line.

Open in new tabDownload slide

3.1 Luminosity at the LC

Figure 2.

Poynting flux as a function of the distance from the central rotator. We present results after 16.7 rotation periods for different damping constants |$\kappa _\Phi$| (and a constant advection parameter |$c_\Phi =1$|⁠) for different resolutions. The tests employing the LCR method are indicated by thick red lines.

Open in new tabDownload slide

Figure 3.

As Fig. 1 but for selected damping constants |$\kappa _\Phi$| and different choices for the advection parameter |$c_\Phi$|⁠.

Open in new tabDownload slide

Fig. 1 lucidly illustrates that the Y-point moves away from the LC with increasing values of α. A Y-point closer to the stellar surface boosts the electromagnetic luminosity of the pulsar (Timokhin 2006). The amount of open magnetic field lines increases for a decreasing dimensionless distance x₀ ≡ r_Y/r_LC of the Y-point from the central object. This is directly related to the angular extension of the polar cap, which can be quantified by measuring the colatitude of the closed zone region. In Fig. 4, we present a large selection of simulated Y-point luminosities (L_Y) versus their estimated Y-point position, excluding results that did not yield equilibrium magnetospheres for very small or very large values of α (see discussion in Section 5). To approximate the Y-point location, we evaluate the drift velocity along the equator and assign x₀ to the position where the velocity is comparable to a small parameter, which equals the grid spacing Δr in magnitude. We employ a similar approximation technique to the evaluation of the closed zone colatitude θ_c (Tables 1 and 2), with a suitably strong decay at the transition to the co-rotation region. We find that, approximately, |$\sin \theta _{\rm c}\sim (r_*/r_{\rm Y})^{1/2} \sim 0.5x_0^{-1/2}$|⁠. Roughly in line with the findings from Timokhin (2006), we measure a correlation between the Y-point location and luminosity. For large values of α, associated to smaller values of x₀, the errors in these measurements – averaged over several pulsar revolutions – become larger. These errors are likely linked to degrading numerical accuracy induced by the increased stiffness of the augmentation equations (3) and (4) for large values of α (see Section 5). We note that the correlation between x₀ and L_Y/L₀ found numerically approaches the theoretical relation of Timokhin (2006), in which |$L_{\rm Y}\propto x_0^{\lambda }$|⁠, with λ = −2.065, more so as the numerical resolution increases. Yet another interesting correlation derived from our results is that L_Y/L₀ ≈ 0.57(sin θ_c)^0.19, which supports the claim that a larger luminosity at the Y-point correlates with a larger polar cap angle. Interpreting our findings, we cautiously suggest a possibility to obtain quasi-stationary magnetospheres in which the Y-point is not located at the LC but inside it. The location of the Y-point depends on the diffusivity at the ECS. Smaller diffusivity moves the Y-point inward and increases the outgoing Poynting flux, as we show in Fig. 4. The root of this dependency is the (magnetic fieldline) coupling between the fieldline footpoints at the stellar surface and part of the ECS (precisely, the part adjacent to the Y-point). The existence of such a region of coupling between the stellar surface and the ECS was suggested, e.g. by Contopoulos (2019) for stationary hybrid solutions combining FFE and a non-ideal ECS. Our time-dependent models reproduce many aspects of these equilibrium states. On the right-hand panel of Fig. 4, we evaluate the specific dependency of the luminosity on the diffusion parameter α. A power law of the form L_Y/L₀ ∝ α^0.1 is found from a fit to our computed models. The parameter α controls the numerical diffusion of the electric field, if we assume that numerical diffusion mimics the physical one to some extent (i.e. α ∝ 1/η, where η is the resistivity, see discussion Section 5 and Mahlmann et al. 2021b). Thus, our results suggest that the luminosity of the aligned rotator is inversely proportional to the resistivity at the ECS. Nevertheless, because of the small power-law index in the aforementioned α–luminosity relation, large changes in the resistivity are required to produce significant variations of the luminosity at the LC.

Figure 4.

Dependence of the luminosity on the Y-point location and diffusion parameter α. For the Y-point location (left-hand panel), measurements from simulations at five different moments in time are averaged (standard deviation indicated in error bars). Fit functions are derived from for all models that find a reliable equilibrium (solid dots) and are represented by coloured dashed lines. As a reference, we compare these results to the corresponding relation in Timokhin (2006, black line). Errors to these fit functions are represented by the lightly shaded regions of the respective colour. Outliers that do not find stable equilibria are denoted by coloured circumferences. The right-hand panel shows the normalized luminosity measured at the LC, evaluated against the diffusion parameter α. A fit function to the data is shown with a thick solid line. Uncertainties to the fit function are given as a shaded region around the central line.

Open in new tabDownload slide

3.2 Energy dissipation beyond the Y-point

Tchekhovskoy et al. (2013) hold a reliable benchmark for the convergence of pulsar magnetosphere modelling and its comparison between FFE and MHD. We emphasize the following property of their results: MHD models show a decay of the Poynting flux beyond the LC that diminishes with higher resolution. However, their FFE models of the aligned rotator magnetosphere show |$43{-}50{{\ \rm per\ cent}}$| dissipation that converge to a stable amount of dissipation (⁠|$43{{\ \rm per\ cent}}$|⁠) for the highest resolution. We observe such behaviour for the case of LCR, as it is marked by the thick red lines in Fig. 2.

In all other cases using the CC method, the dissipation of Poynting flux along the radial direction is significantly less. This observation is also supported by the field line images of Fig. 1. There, field lines notably reconnect beyond the LC in the limit of low α, while the cleaning of numerical errors shapes the ECS over longer distances. Stationary solutions of the force-free aligned rotator magnetosphere involving an ECS (where the FFE approximation does not hold) have been obtained (e.g. Contopoulos 2007; Contopoulos, Kalapotharakos & Kazanas 2014; Contopoulos et al. 2020). The configurations built by Contopoulos et al. (2020) show that there is a region of the ECS that may extend beyond the LC, where magnetic field lines can still be closed. Our time-dependent models certainly reveal such a region, and its properties depend on the model parameters, which ultimately determine the level of (numerical) dissipation of the method. Small-scale structures in the equatorial plane (such as plasmoid-like formations) do not automatically appear by increasing resolution. Rather, it seems that an efficient damping of charge conservation errors allows them to emerge.

Remarkably, not only the smallest luminosity L_LC is attained by the LCR method, but, most importantly, also the largest relative decrease of the luminosity, ΔL/L_LC = [L(r = 5r_LC) − L_LC]/L_LC between the LC and r = 5r_LC. Taking ΔL/L_LC as a measure of the diffusivity of the algorithm, we conclude that the LCR method is more diffusive than the CC method. This conclusion is numerically robust, since duplicating the spatial resolution does not notably change the dissipation beyond the Y-point observed in the individual models.

Figs 2 and 3 show a significant variability of the Poynting flux beyond the LC for small values of |$\kappa _\Phi$|⁠. For most values of |$\kappa _\Phi$|⁠, the Poynting flux beyond the LC is not monotonically decaying, but shows spikes associated to the ejection of plasmoid-like structures that move outwards along the ECS (see Fig. 1). The fact that these blobs of strong currents move outwards even if they are produced inside the LC is relevant: They do not contribute to the growth of the closed magnetospheric region. This finding is in contrast to Komissarov (2006), who claims that part of the plasmoids will move inwards, increasing the size of the closed magnetosphere with time and, hence, pushing the Y-point towards the LC. This assertion has also led Spitkovsky (2005) and McKinney (2005) to suggest that all configurations with x₀ < 1 would be unstable or transitory (see also Kalapotharakos & Contopoulos 2009; Yu 2011; Parfrey, Beloborodov & Hui 2012a). We do not notice such behaviour here and, thus, our results suggest that closed magnetospheric configurations with a Y-point inside the LC may survive for many dynamical times.

In order to compute ΔL/L_LC when it is spatially (and temporarily) variable (i.e. using the CC method with small values of α), it is necessary to smooth out the data taking a suitable moving average. Larger values of |$\kappa _\Phi$| or smaller advection speeds |$c_\Phi$| (i.e. larger values of α) yield magnetospheres with smaller Poynting flux dissipation beyond the LC. We note that the dissipation takes place at the ECS, where the force-free approximation is not strictly valid, specifically because the electric field strength becomes larger than the magnetic one and non-ideal electric fields with |${\boldsymbol E}\cdot {\boldsymbol B}\ne 0$| become dynamically important. Hence, the different amounts of relative dissipation are closely connected with the numerical handling of the force-free constraints and the propagation of errors from the regions where FFE is breached.

Extreme values of α suppress the hyperbolic/parabolic cleaning, letting the charge and the corresponding electric field divergence become misaligned over time. Nevertheless, very strong damping cleans errors very rapidly but effectively decouples the evolution equations of Φ from the underlying system of physical balance laws (in this case, Maxwell’s equations), and renders the scalar potential dynamically negligible. As cases of intermediate α do not show this variability, we empirically demonstrate that there is an optimal range of values of |$\kappa _\Phi$|⁠, corresponding to α ∼ 1. In this range, the CC method works very efficiently, effectively minimizing the dissipation at the ECS. We find strong evidence for a much better conservation of energy flux beyond the LC than what is quoted throughout the literature, where results seem to correspond to the most diffusive case of our parameter exploration, as to say the LRC method. As a matter of fact, for the intermediate values of α ∈ [1.45, 4.63], the current |${\boldsymbol j}_\parallel$| is very efficiently suppressed in the equatorial region beyond the Y-point (Fig. 1). Indeed, this is the reason which explains the significantly reduced dissipation beyond the LC in the models using the CC method and optimal values of α ∼ 1. For the smallest values of α, the variability time-scales are of the order of the polar cap light-crossing time, namely ∼r_*θ_c/c ∼ r_*/2c (see also the discussion in Section 5.4).

4 SOURCES OF DISSIPATION

4.1 The role of the violations of FFE constraints

We compare the final (quasi-)equilibrium of three different numerical treatments of the electric charge in the magnetosphere in models |$(\mathbf {Ba})-(\mathbf {Bf})$|⁠. Since, we have introduced changes in the boundary conditions and in the numerical grid, we confirm the results from the previous section (incorporating the aforementioned modifications with respect to models in Section 3). First, we use the LCR method throughout the entire magnetosphere for two distinct values of α (models |$\mathbf {Ba}$| and |$\mathbf {Bd}$|⁠). Second, we use the CC method in models |$\mathbf {Bb}$| and |$\mathbf {Be}$|⁠. Finally, we apply the HCC method in models |$\mathbf {Bc}$| and |$\mathbf {Bf}$| using different values of α.

Fig. 5 demonstrates that the cases where charge is not conservatively evolved in the entire domain (models |$\mathbf {Ba}$| and |$\mathbf {Bd}$|⁠) produce congruent results independent of the cleaning parameter α. Such solutions have a large amount of reconnecting field lines beyond the LC and a luminosity at r_LC that approaches L_LC/L₀ ∼ 1. The small relative differences observed in the position of the Y-point (⁠|$\lesssim 2{{\ \rm per\ cent}}$|⁠) and the insignificant change in L_LC/L₀ compared with the |$\simeq 10{{\ \rm per\ cent}}$| relative difference in ΔL/L_LC highlight the fact that the variation of α (by a factor of ≈6.6) mostly affects the dissipation dynamics of the ECS and its neighborhood. When charge is evolved by a separate continuity equation, any misalignment of the electric field divergence and the charge density is cleaned out by the scalar potential Φ. As we presented above (Section 3), altering the cleaning parameter α shifts the position of the Y-point and changes the Poynting flux at r_LC. While most of the cases employing the CC method have a well-maintained current sheet beyond the Y-point (contrasting the reconnection in models |$\mathbf {Ba}$| and |$\mathbf {Bd}$|⁠), excessive cleaning (driven by large values of α) can induce variability of the current sheet along the vertical direction. This is obvious when comparing models |$\mathbf {Bb}$| and |$\mathbf {Be}$|⁠, but can also be observed for large values of α in Fig. 1. As we see when comparing models |$\mathbf {Bb}$| and |$\mathbf {Bc}$|⁠, using the HCC scheme stabilises the current sheet on the equator, damping vertical displacements in the ECS. Inside the LC, models |$\mathbf {Bc}$| and |$\mathbf {Bf}$| only slightly differ from their CC counterparts. The location of the Y-point, and the luminosity differ less than |$\simeq 2{{\ \rm per\ cent}}$| among HCC and CC models for α = 0.7, while for α = 4.6, HCC models show values of x₀ and L_LC/L₀ in between of the ones of CC and LCR models (Table 2). Beyond the LC, the differences between HCC and CC models are driven by insufficient damping, especially in the model with larger α (⁠|$\mathbf {Bb}$|⁠), where the ECS is distorted and the violations of the magnetic dominance condition are more patchy and intermittent along it. The fact that HCC models |$\mathbf {Bc}/\mathbf {Bf}$| look more like |$\mathbf {Bb}/\mathbf {Be}$| than like |$\mathbf {Ba}/\mathbf {Bd}$|⁠, respectively, results from a subtle combination of two effects. On the one hand, the corrections to the electric fields after violations of the magnetic dominance condition are larger when using a local reconstruction of charge, either globally (LCR) or locally (HCC), than in CC models. This seems natural as the largest charges are created as a result of the sharp discontinuities of the electric field across the ECS, where LCR and HCC models undergo the same corrections in the charge density. On the other hand, the commonality in the (hyperbolic) evacuation of the errors triggered by violations of the FFE constraints at the ECS results in a greater similarity between HCC and CC models than to LCR models (the feedback of |$\nabla \cdot {\boldsymbol E}$| on Ψ is significantly reduced because of the explicit form in which ρ is constrained in the LCR method; Section 2.1). We can, thus, conclude that the large diffusivity in case of LCR method is primarily induced by corrections to the electric fields after violations of the magnetic dominance condition.

Figure 5.

Comparison of different treatments of the charge ρ for some configurations specified in Table 2. We display the toroidal magnetic field in the left-hand panel of each respective model. The corresponding right-hand panels visualize the force-free violations accumulating during one full time-step of the FFE integration. We show the growth of non-ideal electric fields (absolute value) emerging transiently due to the violation of the |${\boldsymbol E}\cdot {\boldsymbol B}=0$| constraint by blue colours, and locations where the magnetic dominance condition is breached by yellow contours. The vertical solid and dashed black lines denote the position of the LC and of the Y-point, respectively.

Open in new tabDownload slide

The (algebraic) correction of violations to the force-free conditions has different consequences for each of the constraints. Deviations from the |${\boldsymbol E}\cdot {\boldsymbol B}=0$| condition build up continuously and are distributed throughout the domain, though they are larger close to current sheets. We visualize this in the panels of Fig. 5 that show the |${\boldsymbol E}\cdot {\boldsymbol B}$| errors normalized to the local magnetic field strength in a blue gradient. Violations include both very small inaccuracies that result from truncation errors of the algorithm, and strong non-ideal electric fields emerging, for example, at current sheets. In models employing the LCR method (⁠|$\mathbf {Ba}$| and |$\mathbf {Bd}$|⁠), the hyperbolic part of the cleaning still operates to drive violations of the |${\boldsymbol E}\cdot {\boldsymbol B}$| constraint away from the ECS (we note a wave pattern – concentric structures – apparently emerging from the equatorial plane at about r ≈ 1.2r_LC in the respective panels of Fig. 5). That pattern persists in models |$\mathbf {Bb}$| and |$\mathbf {Be}$|⁠, as well as in |$\mathbf {Bc}$| and |$\mathbf {Bf}$|⁠, where we employed the CC and the HCC method. However, in these cases, the violations of |${\boldsymbol E}\cdot {\boldsymbol B}$| can affect a larger region, especially for small values of α. The mid-panel of Fig. 6 shows that models employing the CC method systematically dissipate less energy by Ohmic processes (here approximated by the amount of current parallel to the electric field added up over the whole domain on a single time-step) than models using the LCR method for α ≲ 5.

Figure 6.

Dissipation via different channels of numerical diffusion as a function of the dimensionless parameter α (equation 12) in a quasi-equilibrium state at t ≈ 7.5t_p. An estimate of the favorable range for α that minimizes the dissipation by force-free violations is indicated as a light blue region. Different colours of the markers indicate different resolutions and treatment of charge density evolution. The shapes of the markers vary with different |$c_\Phi$|⁠. We provide reference values of the respective dissipation channel for the LCR method of highest α as star-shaped data-points.

Open in new tabDownload slide

The condition |${\boldsymbol E}^2-{\boldsymbol B}^2\le 0$| is only relevant when significant non-ideal electric fields have built up, in the setup at hand at the ECS, as it is lucidly illustrated by the yellow contours in Fig. 5. It is a natural impulse to associate the high diffusivity observed for the setups using the LCR method to the violations of the magnetic dominance. Looking at the time evolution of the violation of the magnetic dominance constraint, we observe that in models implementing the LCR method everywhere, one finds regions of |${\boldsymbol E}^2\gt {\boldsymbol B}^2$| that consistently cover the ECS in the range x₀ ≲ r/r_LC ≲ 1.7 (and likely beyond). In contrast, models employing the CC and HCC methods display an intermittent set of spots of smaller extension, where |${\boldsymbol E}^2\gt {\boldsymbol B}^2$| in the equatorial region beyond the Y-point. The top panel of Fig. 6 displays the electric energy subtracted from the whole domain during one iteration of the time-integrator. Contributions to this channel of diffusion result from mesh cells where the magnetic dominance constraint is breached. In practice, we calculate averages from several snapshots of the results to ensure that the quoted dissipation estimates have (more or less) stabilized. LCR models computed at the highest resolution (N_r = 128; magenta squares) systematically dissipate more (the least about the same) energy by restoring the magnetic dominance condition than models computed with the CC method and the same resolution (dark green symbols) for α ≲ 3. Similarly to the electric energy dissipation in the correction of the |${\boldsymbol E}\cdot {\boldsymbol B}=0$| constraint, there is an interval 0.5 ≲ α ≲ 3, where dissipative losses induced by the restoration of the magnetic dominance condition are minimized.

In contrast to the trend found for the CC method, the dissipation triggered by violations of the FFE constraints is quite insensitive to α when evaluated for the LCC method (see magenta symbols in the upper and mid panels of Fig. 6). In order to understand this behaviour, we have computed the electric energy dissipated by the hyperbolic/parabolic cleaning algorithm, which is given by |$|c_\Phi ^2{\boldsymbol E}\cdot \nabla \Phi |$| in the lower panel of Fig. 6. While the gradients of the cleaning potential Φ are large enough, increasing α (i.e. damping the divergence errors faster than shifting them away) reduces the diffusion through this channel, independent of the method used to evolve or reconstruct the charge density (LCR or CC). The dissipation through this channel is significantly larger than that driven by violations of the FFE constraints for α ≲ 5 and dominates the overall diffusion of electromagnetic energy in the magnetosphere. Above that value of α, the total dissipated energy in the |$|c_\Phi ^2{\boldsymbol E}\cdot \nabla \Phi |$| channel is smaller than in the other two channels (for this, we compare the numerical values in the top and mid panels of Fig. 6 to those in the bottom panel). From that point on, the cleaning algorithm in the CC method cannot efficiently evacuate and damp the errors induced in |$\nabla \cdot {\boldsymbol E}$| by the restoration of the FFE constraints. Thus, the energy dissipation in models employing the CC method, becomes dominated by the violation of the FFE constraints above a certain value of α. This reasoning substantiates our claim about the existence of an optimal interval for 0.7 ≲ α ≲ 3, where the overall dissipation in the magnetosphere is smaller using the CC method than the LCR one. The insensitivity of the dissipation due to violations of the FFE constrains on α for LCR models (for a fixed numerical resolution) can be explained by the limited evacuation of FFE violations, e.g. from current sheets. As |$\rho \approx \nabla \cdot {\boldsymbol E}$| (up to discretization errors), the right-hand side of equation (4) is significantly smaller than in the models equipped with the CC method; the two mechanisms (cleaning of divergence errors and enforcing FFE violations) operate independently.

4.1.1 Algebraic corrections versus driver currents

We have seen in the previous section that violations of the FFE constraints are fundamentally connected to the magnetospheric structure. It is, therefore, natural to assess whether different methods to restore the FFE constraints impact on the overall dissipation and topology of the magnetosphere. Throughout the literature, various techniques are used to correct any deviations from the |${\boldsymbol E}\cdot {\boldsymbol B}=0$| condition. These techniques split up into the ones that apply algebraic resets of the electric field at each violation (as done in the models shown so far), and those that modify the Ohm’s law encoded in a suitable current to drive the system into a force-free state (cf. Mahlmann et al. 2021a, and references therein). It is very much justified to suspect that all the variability in the luminosity identified in Section 3 stems from this critical ingredient to force-free codes, rather than from their connection to charge conservation. To dissect this subtle issue, we conduct simulations of the models |$(\mathbf {Bg})-(\mathbf {Bh})$|⁠. Without algebraic resets of the |${\boldsymbol E}\cdot {\boldsymbol B}=0$| condition, we employ κ_I = 8 in the force-free current presented in equation (7). In this way, our method is comparable to the ones that employ driving currents, namely, a formulation of Ohm’s law with a finite resistivity σ_∥ that acts along the direction of the magnetic field (e.g. Alic et al. (e.g. Komissarov 2004; Alic et al. 2012).

In Fig. 5, we contrast the magnetospheric equilibrium found when combining a suitable current to drive |${\boldsymbol E}\cdot {\boldsymbol B}\rightarrow 0$| for the LCR method (models |$\mathbf {Bg}$| and |$\mathbf {Bi}$|⁠), and for the HCC method (models |$\mathbf {Bh}$| and |$\mathbf {Bj}$|⁠). One straightforward observation is that a local reconstruction of charge remains the most diffusive limit in view of the much larger decrease of luminosity with distance measured by ΔL/L_LC (Table 2). For the HCC models |$\mathbf {Bh}$| and |$\mathbf {Bj}$|⁠, the luminosity decrease beyond the LC is as small as in the case in which algebraic cutbacks of the electric field enforce |${\boldsymbol E}\cdot {\boldsymbol B}=0$| (models |$\mathbf {Bc}$| and |$\mathbf {Bf}$|⁠). The luminosity of the HCC models that implement a driving current to enforce |${\boldsymbol E}\cdot {\boldsymbol B}\rightarrow 0$| decreases by |$\sim 18{-}25{{\ \rm per\ cent}}$| with respect to models using algebraic resets of the electric field (models |$\mathbf {c}$| and |$\mathbf {f}$|⁠). This is a direct consequence of the larger value of x₀ in models |$\mathbf {Bh}$| and |$\mathbf {Bj}$|⁠, which place their Y-points closer to the LC. In contrast, models using the LCR method display only a small change of the luminosity (⁠|$\lesssim 1{{\ \rm per\ cent}}$|⁠) when comparing models with different strategies to enforce |${\boldsymbol E}\cdot {\boldsymbol B}=0$| (i.e. models |$\mathbf {Ba}/\mathbf {Bd}$| versus models |$\mathbf {Bg}/\mathbf {Bi}$|⁠). Thus, the smaller luminosity of models using the LCR method is not fully accounted for by the algebraic – rather harsh – corrections we apply to non-ideal electric fields.

4.1.2 Diffusivity models beyond the LC

In Fig. 7, we compare the magnetospheric evolution of different values of η_d for the HCC and CC models (Bk)–(Bn) and (Bo)–(Br), respectively. As we also denote in Table 2, we observe – very much like in the previous section - that the Y-point is closer to the LC, and the luminosity at the LC is lower in case of models using the HCC method. Regarding the dissipation of Poynting flux beyond the LC, the whole series of models (Bk)–(Br) display similar and relatively low values of ΔL/L_LC, especially if we compare these values to the analogous ones obtained with the LCR method and η = 0 (models Bg and Bi). For most of the phenomenological resistivities considered here (say η ≤ 10⁻²), HCC models are more diffusive than CC models and this global measurement (of luminosity) manifests itself at the local level in a smoother velocity profile along the ECS. Contrasting this, local reconnection events with large drift velocities (normal to the ECS) become visible in Fig. 7 (bottom row of panels) for the cases using the CC method (especially in models q and r). At the same time, the width of the resistive layer, in which we can identify an inflow (drift) velocity into the ECS, increases with larger η_d.

Figure 7.

Comparison of magnetospheric dynamics for different values of η in the current given by equation (7). From the left- to right-hand panels, the phenomenological resistivity is increasing as η ∈ {10⁻³, 10⁻², 10⁻¹, 1.0}. The background colour denotes the (drift) velocity component |$v$|^y, as to say the velocity pointing perpendicular to the ECS. The separatrices where |$v$|^y = 0 are highlighted by magenta contours. The top row comprises models implementing the HCC method, the bottom row those of fully conservative evolution of the charge density ρ (CC models).

Open in new tabDownload slide

5 DISCUSSION

5.1 Action of non-ideal electric fields in FFE

As there is no discrepancy between the charge density ρ and |$\nabla \cdot {\boldsymbol E}$|⁠, the cleaning potential Φ reduces its role to a true scalar correction of (very small) numerical truncation errors. Non-ideal fields still alter the system of conservation laws by the loss of both, energy conservation, and charge conservation. The charge density is one constituent of the current that acts as a source of Ampère’s law. A source or sink of it, due to the addition of charge when fixing the violation of the |${\boldsymbol E}\cdot {\boldsymbol B}=0$| constraint, is dynamically relevant and may have a global impact as it can be transported away from the numerical cells where it is initially generated.

5.2 Time-scales of the hyperbolic/parabolic cleaning

As has been noted throughout the literature (e.g. Mignone & Tzeferacos 2010; Mahlmann et al. 2019, 2021b), mostly in the context of cleaning errors to the |$\nabla \cdot {\boldsymbol B}=0$| constraint, a careful calibration of the parameters controlling hyperbolic/parabolic cleaning is necessary. The same holds true for the cleaning of errors to the |$\rho =\nabla \cdot {\boldsymbol E}$| constraint, and we conducted a thorough calibration for the simulations shown in this paper. One way of optimizing the cleaning parameter α is to evaluate the dissipation induced by the source terms to the energy evolution equation, as it was presented in Mahlmann (2020, equation 2.73). The relevant channels are dissipation by cleaning (⁠|$\propto c_\Phi ^2{\boldsymbol E}\cdot \nabla \Phi$|⁠), as well as Ohmic heating fueled by non-ideal electric fields (⁠|$\propto {\boldsymbol E}\cdot {\boldsymbol j}$|⁠).

With the results in Fig. 6 we find that employing the Φ cleaning with α ≈ 1 yields optimal results. This means that an optimal regime is reached when the time-scales of dissipation (τ_d) and advection (τ_a) of numerical errors induced by the violation of FFE constraints (in the aligned force-free rotator, predominantly at the ECS) are of the same order. First, the dissipation by corrections of non-ideal electric fields is minimized across all models in this regime (light blue stripe). Second, the dissipation by the cleaning potential stabilizes at a steady value across resolutions. This plateau roughly coincides with the blue shaded region and can be interpreted as a trade-off between excessively weak cleaning (inducing larger dissipation for small values of α because of the oscillatory behaviour of the electric field after applying corrections to enforce the force-free regime) and the extreme case of over-damping that drives the scalar function Φ to zero very rapidly.

5.3 Magnetospheric structure

In Fig. 8, we display the structure of currents and charges in two-selected regions of the magnetosphere. Specifically, we examine a location close to the polar cap and another one around the Y-point for models using the LCR method (model Bd) and CC method (model Hc). Both models have the same value of α = 0.7 and the same numerical resolution. The overall charge distribution looks qualitatively like the one obtained by previous force-free models (e.g. Parfrey et al. 2012a, in axial symmetry or Kalapotharakos & Contopoulos 2009; Kalapotharakos et al. 2012 in three dimensions) and in a number of PIC simulations of aligned rotators (e.g. Chen & Beloborodov 2014; Philippov & Spitkovsky 2014; Cerutti & Beloborodov 2016; Brambilla et al. 2018). Negative charges fill the regions above the polar caps, while positive ones fill the closed magnetosphere up to the Y-point. The structure of the Y-point and of the ECS adjacent to it consists of a set of charge layers of alternating sign stacked vertically in model Bd (Fig. 8, bottom left-hand panel). This layered structure is also observed in, e.g. Parfrey et al. (2012a, cf. their Fig. 16). Along the equatorial region, a positively charged layer with a thickness ∼0.027r_LC emerges, reproducing the results of a positive surface charge density along the ECS shown in Timokhin (2006). In model Hc, the layered charge structure is destroyed due to time-variable episodes of reconnection along the ECS. Still, the positively charged central layer is intertwined with regions of negative charge. The (more) variable structure of the ECS in models using the CC method arises because of the existence of a finite time induced by restoring the stationary condition |$\nabla \cdot {\boldsymbol E}=\rho$|⁠. This finite time is brought about by the (finite) propagation speed (⁠|$c_\Phi \sim c$|⁠) of the hyperbolic part of the equation controlling Φ, and modulated by the (finite) diffusion time-scale τ_d (see Section 5.2). In the LCR method, this time is zero, as to say the restoration of |$\nabla \cdot {\boldsymbol E}=\rho$| is instantaneous. However, the current does not follow the change in the charge density distribution instantaneously.

Figure 8.

Charge and current structure in time-dependent aligned rotator models. Upper and mid panels: charge (ρ; coloured background, normalized to the absolute value of the Goldreich–Julian charge density, |ρ₀|) and current (⁠|${\boldsymbol j}_\parallel$|⁠; arrows, normalized to ρ₀c) structure of selected models (see panel legends). We visualize two different zooms into the magnetosphere, namely the polar cap region (top panels), and the Y-point vicinity (bottom panels). The direction of the current flow is indicated by arrows for regions where j_∥/j₀ > 0.01. The cyan dashed line denotes the location where ρ = ρ₀ (specifically, in between that line and the stellar surface, the charge density is larger than the Goldreich–Julian charge density). Bottom panels: difference between the current |${\boldsymbol j}$| in our time-dependent models and the current in the steady state case |${\boldsymbol j}_{\rm steady}=\nabla \cdot {\boldsymbol B}$|⁠, normalized to the Goldreich–Julian current density j₀. We display the poloidal magnetic field lines; the vertical black line denotes the position of the LC.

Open in new tabDownload slide

The directions of the poloidal current (displayed by the coloured arrows in Fig. 8) shows a return current flowing along the ECS and continuing over the closed zone separatrix converging on the Y-point and extending up to the stellar surface. However, the structure of the current is not simple. Among currents flowing towards the stellar surface, we also find antiparallel currents flowing away from the surface. We point out the interesting observation that above the polar caps a region with super Goldreigh-Julian charge density, i.e. ρ > ρ₀ (enclosed by the cyan dashed line in the upper panels of Fig. 8), emerges.

Under the assumption of stationarity, in a small layer around the ECS with a thickness 2h ≪ r, Contopoulos et al. (2014) constructed a simple analytic model for the electromagnetic structure, where the electric field toroidal component is zero. However, in our models, we find a time-dependent |E_ϕ(r, |$z$|⁠, t)| ≠ 0 in places along the ECS (and x > 1) where the magnetic dominance condition is breached. Indeed, the model using the CC method shows a quasi-periodic pattern of the form E_ϕ(r, |$z$|⁠, t) ∼ E_ϕ(r, h)sin [8|$\pi$|(x − t)]. For the model using the LCR method, the spatial frequency is about two times larger. The other two components of the electric field roughly follow the analytic model of Contopoulos et al. (2014), namely, they follow the relations E_r(x, h) ≈ −xB_|$z$|(x, |$z$|⁠) and E_|$z$|(x, h) ≈ xB_r(x, |$z$|⁠). Because of the presence of displacement currents and the (comparatively much smaller) contributions of the term |$c_\Phi ^2\nabla \Phi$| in equation (2), the current density in the ECS is not given by |${\boldsymbol j}_{\rm steady}=\nabla \times {\boldsymbol B}$| as in the stationary analytic model of Contopoulos et al. (2014). This is lucidly represented in the bottom panels of Fig. 8, where the ratio |$\Delta J/j_0=(|{\boldsymbol j}|-|\nabla \times {\boldsymbol B}|)/j_{0}$| differs from zero (indicating that the current density includes contributions other than |${\boldsymbol j}_{\rm steady}$|⁠) mostly in the vicinity of the ECS and also along the current sheets surrounding the closed magnetospheric region. At this point, there is an interesting difference between models using LCR and CC methods, namely, models with a LCR have a current smaller than |${\boldsymbol j}_{\rm steady}$| (ΔJ/j₀ < 0) along most of the current sheet. In contrast, the model Hc shows smaller average deviations along the ECS between |${\boldsymbol j}$| and |${\boldsymbol j}_{\rm steady}$|⁠, and surrounding it in intermediate patches beyond the LC. That finding is unexpected given the larger variability along the ECS encountered, in general, for CC models.

5.4 Diffusivity in force-free magnetospheres

Resetting the charge density ρ to the instantaneous value of |$\nabla \cdot {\boldsymbol E}$| assures compatibility of electric fields and the corresponding charge distribution. However, removing significant fractions of the electric field from the domain for ensuring FFE conditions alters the charge distribution in the domain by an amount that is not necessarily of the order of the truncation error. In our default (CC) scheme (Mahlmann et al. 2021a), the charge is evolved with a separate continuity equation. Changes to the electric field by the algebraic reset of force-free errors introduce a misalignment between the charge density and the divergence of electric fields. The degree of localization of this misalignment to inherently non force-free regions, such as current sheets, can be interpreted as a diffusive length-scale. We control this localization in our hyperbolic/parabolic cleaning procedure with the parameter α, that expresses the ratio of diffusion to advection time-scales for the cleaning errors.

The LCR method turns out to be the most diffusive one, using as a measurement of the diffusivity the dissipation of Poynting flux beyond the LC. In the models resorting to the LCR method, the effect of violations to ideal FFE conditions spreads throughout the domain, with notable consequences for the global energetics. While the Y-point relaxes to a location close to the LC (cf. Fig. 5), field lines reconnect across the ECS and induce a notable dissipation of the outward transported energy (⁠|$\sim 30{{\ \rm per\ cent}}$|⁠). When augmented by a conservative evolution of the charge density with a suitable calibration of the divergence cleaning parameter α, reconnection beyond the LC is greatly reduced – as is the dissipation of energy. At the same time, the Y-point is pushed away from the LC (cf. Fig. 5), and the overall luminosity increases. Though the action of the cleaning potential Φ can become excessively large for α ≫ 1, the increase of the luminosity with the corresponding change of the Y-point location roughly follows the trend that is theoretically expected for small and intermediate values of α (cf. Fig. 4).

We explore an interpretation of our results in light of the parameterization of the diffusivity beyond the LC in terms of the pair formation multiplicity κ suggested by Contopoulos et al. (2020). These authors show that smaller values of κ eventuate in more dissipation at the ECS and, conversely, κ ≫ 1 yields very little or no dissipation at all. When drawing a direct comparison to the results presented throughout this paper, the models endowed with the LCR method would correspond to values of κ ∼ 1. Models employing our the CC scheme (our default), would compare to values κ ≫ 1. The physical conditions in a typical pulsar magnetosphere tend to produce many pairs per Goldreich–Julian charge particle in the polar cap (e.g. Timokhin & Harding 2015; Contopoulos 2019, and references therein). Hence, in this regard, our CC models potentially reproduce the usual conditions met in actual rotating pulsars more closely, though the limitations set by the force-free regime naturally persist.

The panels of Fig. 5 that show the location of violations of the FFE constraints illustrate two notable aspects: (i) deviations from |${\boldsymbol E}\cdot {\boldsymbol B}=0$| are not a localized phenomenon. Emerging from genuinely non-ideal regions, such as current sheets, non-ideal fields can spread throughout the domain. (ii) A change in the treatment of charge conservation not only minimizes the dissipation induced by corrections of the |${\boldsymbol E}\cdot {\boldsymbol B}=0$| condition, but also changes the electromagnetic structure of the current sheet itself. Cases |$\mathbf {Ba}/\mathbf {Bd}$| show a region with |${\boldsymbol E}^2-{\boldsymbol B}^2\gt 0$| along the length of the current sheet. Contrasting this, such violations only occur at X-points in the current sheet of models |$\mathbf {Bb}/\mathbf {Be}$|⁠. The acceleration of particles and the production of radiation demand the existence of an electric field component parallel to the magnetic field, namely that |${\boldsymbol E}\cdot {\boldsymbol B}\ne 0$|⁠. Hence, the structure of the regions of the magnetosphere, where the strongest violations of the |${\boldsymbol E}\cdot {\boldsymbol B}=0$| condition occur, are likely closely related to the production of pulsar radiation (e.g. Timokhin & Arons 2013). Looking at Fig. 5, these violations, although extended in the magnetosphere as stated above, are maximized in the vicinity of the Y-point, suggesting that the Y-point is the most important site for particle acceleration in the magnetosphere. This result is backed up by PIC simulations of, e.g. Chen & Beloborodov (2014) and gives support to the theoretical ‘ring-of-fire’ model of Contopoulos & Stefanou (2019).

Violations to the ideal force-free conditions in combination with the transport of charge conservation errors originating during their correction (Section 5.1) are the main driver of diffusivity in force-free aligned pulsar magnetospheres. The relative importance of the dissipation triggered by the enforcement of either the |${\boldsymbol E}\cdot {\boldsymbol B}=0$| or the |${\boldsymbol E}^2-{\boldsymbol B}^2\lt 0$| conditions is similar across these two particular channels, as we can observe in the magnitudes of the electric energy lost within a given time-step in Fig. 6. However, the magnitude of the dissipation triggered by each of them may depend on the order and frequency with which they are applied, as well as on the mesh and time-integrator (cf. Spitkovsky 2006). Nevertheless, the dominant contribution to the dissipation along these channels stems from the ECS, where the magnetic dominance condition is chronically breached. We suggest that the magnetic dominance condition is really the origin of the differences between the various methods of dealing with the charge treatment as presented in this paper. In the HCC models using LCR only in grid zones where |${\boldsymbol E}^2-{\boldsymbol B}^2\gt 0$|⁠, the dissipation through this channel is minimized (see light green coloured squares in the upper panel of Fig. 6), and the Ohmic dissipation is maximized (mid-panel of Fig. 6). Since |${\boldsymbol E}^2-{\boldsymbol B}^2\gt 0$| is only reached at the ECS, it is standing to reason that the specific restoration of the magnetic dominance constraint may induce global changes in the magnetospheric structure, its luminosity, and, certainly, on the amount of dissipation of Poynting flux beyond the LC.

The resistivity models used beyond the LC in Section 4.1.2 provide twofold insight. First, they allow us to estimate the numerically induced diffusivity η₀ across the ECS. The effect of η_d will only become noticeable when the phenomenological resistivity is larger than the numerical diffusivity of the method. From Fig. 7 we can estimate that η₀ ≲ 10⁻¹. Second, and in line with the results presented in Mahlmann et al. (2021b), a choice of η_d ≳ 10⁻¹ in the current presented in equation (7) allows to properly model dynamics of the resistive layer of the ECS. Increasing η_d gradually drives relatively large-scale inflows into the ECS (as traced by the drift velocity component perpendicular to it), effectively mimicking physical dynamics in the non-ideal region. As we established very competitive convergence of our high-order FFE method (Mahlmann et al. 2021b), we suggest that this relatively large value of η₀ is, indeed, induced by the non-ideal fields emerging in the ECS. In this context, it becomes clear why several FFE methods need to employ special treatments of the ECS in the pulsar magnetosphere (McKinney 2006; Etienne et al. 2017), namely, to reduce the extent of the diffusive regions by an ad hoc prescription.

As stated in, e.g. Timokhin & Arons (2013), a necessary criterion imposed by observational constraints on the pulsar magnetosphere is the stationarity in a statistical sense. In other words, any local fluctuations on time-scales smaller than the LC light-crossing time τ_LC = r_LC/c (or more likely, over the rotational period of the pulsar t_p) that average to a stationary state, may also account for the stability of the pulsar mean profiles and sharpness of the peaks in the spectra of gamma-ray pulsars. A very salient feature related to the treatment of the charge conservation equation in FFE is the obvious time-dependence of the magnetosphere, driven by episodes of magnetic reconnection along the ECS. However, the magnetospheres resulting from the CC method (including a suitable conservative treatment of the charge) are stationary if we average them out over time-scales comparable to τ_LC. In order to support this statement, we show the time-averaged map of the toroidal magnetic field over an interval Δt slightly larger than one pulsar rotational period in Fig. 9 (left-hand panels). We find even stronger evidence than in the polar distribution of the toroidal field when evaluating the Poynting flux as a function of distance for the averaged models. It displays a radial dependence which radically smooths the spatial variability (Fig. 9, right-hand panel). We further probed the long-term stability of the representative models |$\mathbf {Ba}$|⁠, |$\mathbf {Be}$|⁠, and |$\mathbf {Bb}$| by tracking their evolution during ≳ 30 rotational periods (Mahlmann & Aloy 2021). The models show stability over such time-scales in all the characteristic properties discussed in the previous section, especially regarding the Y-point location and pulsar luminosity.

Figure 9.

Averaging of magnetospheric fields and Poynting fluxes during a time of Δt = (7.5 − 6.8)t_p = 1.3t_p for selected models (⁠|$\mathbf {Ba}$|⁠, |$\mathbf {Be}$|⁠, and |$\mathbf {Bb}$|⁠). Small variations beyond the LC, especially those observed in Fig. 2 (shown again as thick, transparent lines in the background of the right-hand panel) are compensated across time-scale Δt ≳ 1t_p.

Open in new tabDownload slide

6 CONCLUSIONS

In a deep exploration with our recently developed force-free code (Mahlmann et al. 2021a,b), we exploit the diffusion time-scale induced by hyperbolic/parabolic cleaning of charge conservation errors (Section 5.2) to quantify an aspect that has not been systematically assessed so far in many FFE simulations, including our own. Namely, the global imprint of local violations to the force-free constraints. At the example of the force-free aligned rotator magnetosphere, we demonstrate that balancing the amount of damping and advection of charge conservation errors (encoded in the parameter α) can alter the global structure of the simulated magnetosphere. Specifically, by decreasing the amount of numerical diffusion arising from violations to the FFE constraints, the Y-point moves away from the LC while the outgoing Poynting flux increases by a factor of a few (Section 3).

In summary, our exploration clarifies several technical aspects that should become central for the assessment of (global) FFE simulations. First, the localization of force-free violations to small regions in resistive layers, such as current sheets, are key to reduce the diffusivity that is induced into FFE by non-ideal electric fields or breaches of magnetic dominance (Section 5.1). In our method, we achieve and control such a localization by combining a conservative evolution of charge density with a hyperbolic/parabolic cleaning of errors to Gauss’ law. We suggest that α ≲ 1 is an optimal parameter for the minimization of the combined channels of numerical diffusion (Section 4.1). Second, in the inherently non-ideal aligned rotator magnetosphere, the ECS is the main source of numerical diffusion by inducing strong field gradients and violations to the FFE constraints. We identify a strong dependence of the dynamics on the specific treatment of violations to the magnetic dominance condition (Section 5.4). Third, the extreme nature of algebraic corrections to the force-free conditions are not the main reason for the luminosity dependence of the global magnetosphere. Driving currents yield very similar results (Section 4.1.1). Finally, different treatments of force-free violations, especially at Y-points and current sheets, are likely to change the resistive time-scales of the evolution and to have a notable impact on the equilibrium magnetospheres. We extend our analysis to so-called phenomenological resistivity models, where adapted driving currents mimic the development of resistive layers around current sheets (Section 4.1.2).

FFE is a robust way to model energy flows in highly magnetized plasma, as we find in the magnetospheres of many astrophysical objects. This statement can safely be extended to situations in which the global dynamics of field lines drives the transient appearance of inherently non-ideal regions, such as current sheets. Specifically, we argued this in the context of accretion of magnetic loops on to rapidly spinning BHs (Mahlmann et al. 2020), and the shearing of fields driven by interacting Alfvén waves (Ripperda et al. 2021). However, in situations where areas of genuine (physical) resistivity drive the global field line dynamics, employing FFE methods has to be carefully benchmarked. This work is, to our knowledge, the first extensive calibration of an FFE method for the specific application of astrophysical magnetospheres. We find the global flows of energy to be extremely sensitive to the treatment of FFE violations. Our results agree with comparable numerical surveys throughout the literature only in the limit of strong numerical diffusion induced by the ECS.

Finally, we suggest that scenarios like the aligned rotator should be handled with care when used as a standard test for FFE methods. The stability and magnitude of the Poynting flux beyond the LC can be used as a primer for the diffusivity of the respective method. It could be argued that operating in the well established, but ultimately limited, regimes of ideal fluid approximations for the modelling of global scenarios that are affected or even driven by genuinely resistive effects will require more and more care in the future. The desire to overcome many orders of scale separation has to go hand-in-hand with a deep understanding of the diffusive properties of the employed numerical methods. Until we can transition into an era of multiregime astroplasma codes, we find it reassuring to have the limits of our FFE method laid out transparently.

ACKNOWLEDGEMENTS

This work has been supported by the Spanish Ministry of Science, Education and Universities (PGC2018-095984-B-I00) and the Valencian Community (PROMETEU/2019/071). JM was partially supported by NASA grant 80NSSC18K1099 and NSF grant AST-1909458. We furthermore thank for support from the COST Actions PHAROS CA16214 and GWverse CA16104. This manuscript relies vastly on high performance computing resources. They were provided with the LluisVives machine at the Servei d’Informàtica of the Universitat de València (financed by the FEDER funds for Scientific Infrastructures; IDIFEDER-2018-063) and extensively supplemented by allocations on the MareNostrum and Tirant supercomputers of the Spanish Supercomputing Network (AECT-2021-1-0006, AECT-2021-1-0007).

DATA AVAILABILITY

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

Footnotes

REFERENCES