# A High-Order WENO-based Staggered Godunov-type Scheme with Constrained Transport for Force-free Electrodynamics

###### Abstract

The force-free (or low inertia) limit of magnetohydrodynamics (MHD) can be applied to many astrophysical objects, including black holes, neutron stars, and accretion disks, where the electromagnetic field is so strong that the inertia and pressure of the plasma can be ignored. This is difficult to achieve with the standard MHD numerical methods because they still have to deal with plasma inertial terms even when these terms are much smaller than the electromagnetic terms. Under the force free approximation, the plasma dynamics is entirely determined by the magnetic field. The plasma provides the currents and charge densities required by the dynamics of electromagnetic fields, but these currents carry no inertia. We present a high order Godunov scheme to study such force-free electrodynamics. We have implemented weighted essentially non-oscillatory (WENO) spatial interpolations in our scheme. An exact Riemann solver is implemented, which requires spectral decomposition into characteristic waves. We advance the magnetic field with the constrained transport (CT) scheme to preserve the divergence free condition to machine round-off error. We apply the third order total variation diminishing (TVD) Runge-Kutta scheme for the temporal integration. The mapping from face-centered variables to volume-centered variables is carefully considered. Extensive testing are performed to demonstrate the ability of our scheme to address force-free electrodynamics correctly. We finally apply the scheme to study relativistic magnetically dominated tearing instabilities and neutron star magnetospheres.

###### keywords:

instabilities – magnetic fields – methods: numerical – pulsars: general.^{†}

^{†}pagerange: A High-Order WENO-based Staggered Godunov-type Scheme with Constrained Transport for Force-free Electrodynamics–A

^{†}

^{†}pubyear: 2010

## 1 Introduction

Many high energy astrophysical objects such as black hole magnetospheres, ultrastrongly magnetized neutron stars (magnetars) and probably relativistic jets in active galactic nuclei and gamma ray bursts involve relativistically magnetically dominated plasma. In those situations, the magnetic energy density conspicuously exceeds the thermal and rest mass energy density of particles. Magnetic fields play crucial roles in the dynamics of these astrophysical scenarios. They drive the outflows from astrophysical black holes, neutron stars and the surrounding accretion disks. The magnetic dissipation can also give rise to remarkable non-thermal emissions in the these high energy astrophysical phenomena. Force-free electrodynamics are believed to play an important role in those regimes. As the magnetic energy density greatly exceeds the rest mass energy, conservative magnetohydrodynamic simulation often crashes in such circumstances. However, force-free electrodynamics can behave well in such extreme magnetically dominated situations for the less important terms, such as the inertia and pressure, are entirely ignored in the force-free electrodynamics formulation. In force free electrodynamics the lorentz force disappears, where and are charge and current densities. This approximation allows us to understand the field structure of magnetospheres without solving for the plasma dynamics. However it is still generally difficult to solve FFE to find even stationary solutions. The first solutions of this kind came out for the axisymmetric aligned rotator (Contopoulos, Kazanas & Fendt 1999, hereafter CKF). The force-free constraint can be cast into an elliptic pulsar equation, which specifies the equilibrium configuration of the magnetic flux as a function of the poloidal current. CKF assumed that the last closed field line extends to the light cylinder. This statement was recently relaxed (Goodwin et al. 2004; Contopoulos 2005; Timokhin 2006). These new steady-state solutions differ in spin-down energy loss, structure of current sheets and Y-point demarcating the magnetosphere and have been important for our understanding of pulsar physics. But analytic solutions can not address the stability and variability problems due to the assumption of stationarity. McKinney (2006) performed force-free electrodynamic simulations and found a unique stationary solution for the axisymmetric rotating pulsar magnetospheres. This solution is similar to the solution found by Contopoulos et al. (1999). But the force-free electrodynamic simulations of Komissarov (2006) failed to produce such results due to the high numerical dissipation in his code. Spitkovsky (2006) and Kalapotharakos & Contouplous (2009) solved the force-free electrodynamic equations via the Finite Difference Time Domain (FDTD) method. Due to the wide applications of force-free electrodynamics, such as magnetospheres around black holes, pulsars and accretion disks (Goldreich & Julian 1969; Blandford & Znajek 1977; Contopoulos, Kazanas & Fendt 1999; Lynden-Bell 2003), force-free jets stabilities (Narayan et al. 2009), we are motivated to develop a force-free electrodynamics code, taking advantage of recent progress in computational fluid dynamics (Toro 1997).

Recent years have seen great progress of computational fluid dynamics, especially, the high order Riemann-solver based Godunov schemes. They have many desirable features and are currently believed to be superior to many other numerical schemes for hyperbolic systems. Godunov-type schemes evolve the cell-centered physical variables by incorporating the interactions between neighboring cells. In the simplest case, the first order Godunov scheme (Godunov 1959), the distribution of conserved quantity inside each grid cell is assumed to be a constant. The one-dimensional interaction between two such distributions is the classical Riemann problem (Toro 1997). High order schemes achieve high order accuracy by using high order approximations to compute the interface values and temporal updates. The first second-order Godunov scheme, Monotonic Upwind-centered Scheme for Conservative Laws (MUSCL), was proposed by Van Leer (1979). Piecewise Parabolic Methods (PPM) are high order extension of MUSCL schemes, which introduce parabolae as the basic interpolation functions in a zone allowing for a more accurate representation of smooth spatial gradients, as well as a steeper representation of discontinuities (Colella & Woodward 1984). Essentially Non-Oscillatory (ENO) schemes (Harten et al. 1987) provide a uniformly high order accurate reconstruction, which is total variation bounded (TVB) reconstruction and give robust solutions for flows with discontinuities. They choose only one smoothest stencil out of a set of possible stencils of a fixed length.

Unlike Essentially Non-Oscillatory (ENO) schemes, Weighted Essentially Non-Oscillatory (WENO) schemes (Shu 1997) take a linear combination with coefficients, called weights, of all possible stencils of a given size. The weights in the linear combination add up to unity and are distributed in such a way that stencils that contain discontinuities get extremely small weight. Further, the weights are designed in such a way that when the function is smooth in all stencils, the weights become close to the optimal ones so that the resulting linear combination of the stencils gives a higher order approximation - the same order as the one that the larger, combined, stencil would give. WENO-based schemes are advantageous because they are able to both capture shocks and accurately resolve complex smooth flow structure. Several groups have had success in developing WENO-based schemes in application to relativistic astrophysics (Zhang & MacFadyen 2006) and cosmology (Feng et al. 2004).

The divergence-free evolution of magnetic fields is of vital importance to the correct force-free electrodynamics (Brackbill & Barnes 1980). Unfortunately, straightforward application of Godunov schemes can not produce good results for electrodynamics as the divergence free constraint for the magnetic field can not be automatically satisfied. Komissarov (2004) adopted an augmented system, where the divergence constraint is replaced by an additional evolutionary equation designed to minimize accumulation of error in any one location. In this method may be transported and dissipated like other variables. His method obtained unphysical fast reconnection in the current sheet that led to closed field line far beyond the light cylinder (Komissarov 2006).

Alternatively, we tried a different method, called constrained transport (CT), to satisfy the divergence free constraint. This method is based on a staggered configuration of the magnetic field and the electric field. If the initial magnetic field has zero divergence, then every time step will maintain that to the accuracy of machine round-off error. To construct a numerical scheme based on the CT method and Godunov method, it is important that the volume-centered variables and the face-centered variables be coupled in a consistent manner. Usually the arithmetic average of face-centered values is used and second order accuracy can be achieved (Toth 2000). Recently, Li (2008) proposed a new third order method to map the face-centered variables to the volume-centered variables. In this paper, we would like to make use of both the high-order Godunov scheme and the new mapping method for the CT method and construct a WENO-based Godunov-type scheme with constrained transport for the study of force-free electrodynamics.

This paper is structured as follows: In 2, we summarize the basic force-free electrodynamic equations. We give a basic description of the algorithms in 3. Various numerical tests are given in 4 and applications to physical models are discussed in 5.

## 2 Evolution Equations of Force-Free Electrodynamics

When the plasma inertia and pressure are small, the balance of forces on the plasma is , where and are the charge and current densities. The Maxwell equations together with this force-free constraint read (Gruzinov 1999; Blandford 2002):

(1) |

(2) |

(3) |

where . Note that we have set the light speed and throughout this paper. Equation (3) is a prescription for the plasma current to fulfill the force-free constraint. The two terms in the equation have simple physical meanings: the perpendicular current is denoted by the first term. The second term is the parallel component of the current. The current density is essential for nonlinear interactions. In practical numerical simulation, the second parallel term of current is cumbersome to calculate because it needs the interpolation of both fields and field derivatives. As an alternative way, we update the force-free Maxwell equations with only the perpendicular component of the current, and then remove the accumulated parallel component of the electric field. This procedures achieve the same purpose of simulating a perfectly conducting plasma as the full equation (3) in a simpler way.

## 3 The Algorithm of Force-Free Electrodynamics

The two dimensional force-free electrodynamic equations can be cast in the following compact form,

(4) |

where the primitive variables , the fluxes , and the source terms are

The two Jacobians , are

### 3.1 Riemann Solver

Our Godunov scheme requires that the numerical fluxes at the cell interfaces be constructed via the solution to the Riemann problems. The calculation of numerical fluxes needs the eigen-information of the Jacobi matrices and . The eigenvalues and normalized left and right eigenvectors of the matrix are

(5) |

where the prime denotes transposition. Similarly for the matrix , we have

(6) |

The -direction flux of the Roe’s Riemann solver is constructed as follows (e.g., Roe 1981; Toro 1997):

(7) |

where are the coefficients of the projection of onto

(8) |

Here are the th right eigenvector at the position of . and are the WENO spatial reconstruction of the primitive variables at the position of . The superscripts and denote left and right states, respectively. The one dimensional WENO spatial reconstruction will be addressed in detail in the following sections. Note that our multidimensional spatial reconstruction is carried out in a dimension by dimension fashion. Consequently, the -direction flux at the position can be constructed in a similar way. We may define the following operator ,

(9) |

which will be used in the total variation diminishing (TVD) Runge-Kutta temporal update described in section 3.4.

### 3.2 High-Order WENO Spatial Reconstruction

We would apply Weighted Essentially Non-Oscillatory (WENO) reconstruction in our scheme (Shu 1997). Here we give a brief yet self-contained description of the one dimensional reconstruction method. Here we use to represent one arbitrary component of primitive variables . Consider one dimensional grid along the axis, for one cell , the cell center is at . High-order schemes are built upon reconstructing some set of quantities from the cell center to the cell interface. We use the values of in several grid cells around the grid cell in which the reconstruction is being performed; this set of grid cells is called the stencil. Then, we combine the reconstructed profiles inside each of the grid cells to obtain the global reconstruction . Now we can easily obtain the cell-interface representation of by evaluating at the locations of the interfaces.

The interface value of and the interface value of are reconstructed as follows,

(10) |

In the above two equations, (here ) stands for th linear combination of , with combination coefficients being . In parallel, (here ) stands for th linear combination of , with combination coefficients being . Note here that . The values of for the case of are give in Tab 1. The actual value of are weighted sum of with corresponding weight . We may perform similar weighted summations of to get the value of with corresponding weight .

(11) |

For the case of , the weights and are defined as follows:

(12) |

(13) |

where

(14) |

(15) |

and

(16) |

(17) |

Following the above spatial reconstruction procedures, we could get the left and right reconstructed primitive variable values and at the cell interface . Upon the reconstructed primitive variables and , we could build the fluxes , via the Riemann solver, at the cell interface .

To perform reconstructions in more than one dimension, we use the dimension by dimension approach. It uses one dimensional reconstructions, described in this section, as building blocks for a multidimensional reconstruction.

k | r | j=0 | j=1 |
---|---|---|---|

2 | -1 | 3/2 | -1/2 |

2 | 0 | 1/2 | 1/2 |

2 | 1 | -1/2 | 3/2 |

### 3.3 Constrained Transport

Three types of methods are now applied to maintain the divergence free constraint in Godunov schemes. The first is to use a projection to clean the magnetic field of any divergence after each time step (e.g. Balsara & Kim 2003). The second is to add an evolutionary equation for the divergence to minimize the accumulation of error in the computation domain (Komissarov 2006). The third is the constrained transport (CT) method, which designs the magnetic field difference equations to explicitly preserve the divergence free condition and is the method adopted here. We use the staggered mesh technique, which was proposed by Balsara & Spicer (1999), to preserve the divergence free condition of the magnetic field. We define the magnetic field components on the face centers and all the other quantities are defined at cell centers. The method for constructing fluxes of the face-centered field (the fluxes are located at grid cell corners) based upon the fluxes of the cell-centered field (the fluxes are located at grid cell faces and are returned by a Riemann solver) proposed by Toth (2000) is adopted in our scheme.

In the CT method, the integral form of the induction equation is based on area rather than volume averages. Starting from the differential form of the induction equation (1), the magnetic fields can be updated via

(18) |

(19) |

where and are defined at the face centers, the flux is located at the cell corners. One can check that the above integration can maintain interior to a grid cell at time if it was zero at time .

To calculate the fluxes for cell-centered variables, the Godunov methods still need the magnetic field components that are defined at the cell centers. Usually the cell-centered magnetic fields are defined to be the average of the face-centered values, which is sufficient for second order accuracy. To achieve third order accuracy, we adopt the method proposed by Li (2008) to construct the cell-centered values from the face-centered values. The detail of the mapping from the face-centered values to the cell-centered values is described in the Appendix.

### 3.4 Runge-Kutta TVD Temporal Integration

In order to achieve high order accuracy in time, we choose to update the variables using the high order TVD Runge-Kutta time integration (Shu & Osher 1988) to get from , which combines the first order Euler steps and involves prediction and correction. For example, the third order accuracy in time is achieved by

(20) |

(21) |

(22) |

where is the operator defined in the equation (9) of the section 3.1. Our explicit scheme is subject to the Courant-Friedrich-Levy (CFL) condition. For two dimensional problem, the time step is determined by

We usually choose a CFL number between 0.5-0.8.

## 4 Numerical Tests

In this section, we present some one dimensional test problems
(Komissarov 2002, 2004) to validate the
implementation of this scheme for force-free electrodynamics.

Fast wave : The initial solution is

Fig. 1 shows the results for a fast wave propagating in the positive direction. Both the head and tail of the wave are more resolved than in Komissarov (2002).

Non-degenerate Alfvén wave : In this problem, we consider the following initial solution,

where and are measured in the wave frame that is moving relative to the grid with speed .

Degenerate Alfvén wave : the initial solution is

where

where and are measured in the wave frame that is moving relative to the grid with speed .

Figs. 2 and 3 show the simulation results of the non-degenerate and degenerate Alfvén wave. The scheme captures both fast and Alfvén waves well.

Three wave problem : The initial condition is

Fig. 4 shows the results of the three-wave problem. The initial discontinuity at evolves into two fast discontinuities and one stationary Alfvén wave.

A problem that evolves to : The initial condition is

In between the range , a linear transition layer is set up to connect left and right state.

Fig 5 shows that, in the transition layer, decreases with time and finally the condition is violated.

Stationary Alfvén wave : The initial condition is

Fig 6 shows that the wave is much more resolved than in Komissarov (2004), indicating the effective diffusion coefficient is quite low.

Current sheet : The initial condition is

## 5 Physical Models

### 5.1 Neutron Star Magnetosphere Simulation

In this section, we aims to get the structure of the neutron star magnetosphere via numerical simulations. Similar problems have been investigated by different methods (Komissarov 2006; McKinney 2006; Spitkovsky 2006; Kalapotharakos & Contopoulos 2009). We revisit this problem in order to see the performance of the scheme for more complicated field structures. We solve force-free electrodynamic equations in spherical coordinates on a uniform grid. The initial magnetic field is assumed to be an axisymmetric dipole. The poloidal magnetic field is specified by the poloidal magnetic flux function via

(23) |

The poloidal magnetic flux function is

(24) |

where is the magnetic dipole moment. The electric field on the star is set to to simulate a rotating conducting sphere. We then follow its evolution as the neutron star rotation is switched on. The radial inner and outer boundaries are at and , where . In our simulation, is set to 0.1 and the light cylinder is at . The inner boundary conditions are the same as those proposed by McKinney (2006). The outer radial boundary conditions proposed by Godon (1996) are adopted as our outer radial boundary. Fig. 9 shows the poloidal magnetic field lines of neutron star magnetosphere. We can see that the solution consists of a closed field line region extending to the light cylinder and an open field line region. This solution is quite similar to that found by Contopoulos et al. (1999). Our simulation result also shows that the outer boundary condition performs quite well. Due to this boundary condition, we can put the outer boundary (at ) relatively close to the light cylinder. Fig. 9 shows that the closed magnetic field line beyond the light cylinder in Komissarov (2006) is avoided in our simulations, which means that our scheme has a lower dissipation.

### 5.2 Tearing Instability in Relativistic Magnetically Dominated Plasma

As an illustrative application of our scheme, we study the tearing instability in relativistic magnetically dominated plasma (Lyutikov 2006). In this simulation, to account for the non-ideal effects of resistivity, we adopt the same Ohm’s law as Komissarov, Barkov & Lyutikov (2007). We carry out two dimensional simulations in the plane (all quantities have only dependence on and , but no dependence on ). The computational domain is , with and . We impose periodic boundary conditions at and zero gradient conditions at . The initial force-free equilibrium current sheet is:

(25) |

(26) |

We adopt the units such that and . The initial equilibrium current sheet is normal to the y-axis and its symmetry plane is . In our simulations, the resistivity is set to . The initial equilibrium state is perturbed by adding the following perturbations to the initial equilibrium states :

To get accurate solutions, we must have the resistive sub-layer well resolved. Far from the resistive sub-layer the resolution constraint is not that stringent. We adopt a variable resolution in the direction in the same manner as Komissarov et al. (2007). In Fig. 10, we show the time evolution of the current sheet with . The current sheet bocomes thinner in the middle of the computational domain and becomes thicker at the boundaries. Two large magnetic islands develop subsequently. By the end of the simulation a third smaller island forms in the middle. The simulation shows that our scheme well captures all the detailed physical features inherent in this problem. Our simulation results also confirm the results by Lyutikov (2006) that the shortest growth time is equal to the geometric mean of the Alfvén and resistive time-scales, which is a well-known result in the solar physics (Priest & Forbes 2000).

## 6 Conclusions and Future Work

We have implemented the WENO-based staggered Godunov scheme to study the force-free electrodynamics. We advance the magnetic field with the constrained transport (CT) scheme to preserve the divergence free condition to machine round-off error. We apply the third order TVD Runge-Kutta scheme for the temporal integration. The mapping from face-centered variables to volume-centered variables is carefully considered. A good variety of testings are performed to demonstrate the ability of our scheme to address force-free electrodynamics correctly. We also apply the scheme to study the neutron star magnetosphere. A force-free solution is found for the neutron star magnetosphere with a dipole surface field, the unphysical fast reconnection encountered by Komissarov (2006) is avoided. We finally investigate the tearing instability in the relativistic magnetically dominated plasma. Both primary and secondary magnetic island are well resolved in our simulations. Numerical results of the astrophysical models show that the scheme is robust and accurate.

Strong field electrodynamics is proposed by Gruzinov (2008) recently. It is interesting to study different prescriptions of resistivity using our current scheme, which is left for future work.

This paper mainly concerns with the two dimensional problems. Three dimensional extension of the scheme is straightforward. Another potential advantage is in its extension to the adaptive mesh refinement (AMR) grid, which would allow us to study relevant problems with much higher spatial resolution and larger spatial extent.

When the drift velocity exceeds the speed of light, the force-free electrodynamic approximation breaks down. Under such circumstances, the dynamical effects of matter becomes important and can no longer be ignored. To enter this regime, we need to use the model of relativistic MHD proposed by Lyutikov & Uzdensky (2003) and Lyubarsky (2005).

## Acknowledgments

I thank Dr. Jun Fang for some helpful discussions. This work was funded by the Natural Science Foundation of China (10703012 and 10778702) and Western Light Young Scholar Program. The computation is performed at HPC Center, Kunming Institute of Botany, CAS, China.

## References

- (2004) Balsara D. S., & Kim J., 2004, ApJ, 602, 1079
- (1999) Balsara D. S., & Spicer D. S., 1999, J. Comput. Phys., 149, 270
- (2002) Blandford R. D., 2002, in Gilfanov M., Sumyaev R., & Churazov E., eds, Proc. MPA/ESO, Lighthouse of the Universe: the most Luminous Celestial Objects and Their Use for Cosmology. ESO, Garching, p. 381
- (1977) Blandford R. D., & Znajek R. L., 1977, MNRAS, 179, 433
- (1980) Brackbill J. U., Barnes D. C., 1980, J. Comput. Phys., 35, 426
- (1984) Colella P., Woodward P. R., 1984, J. Comput. Phys., 54, 174
- (2005) Contopoulos I., 2005, A&A, 442, 579
- (1999) Contopoulos I., Kazanas D., & Fendt C., 1999, ApJ, 511, 351
- (2004) Feng L. L., Shu C. W., Zhang M., 2004, 612, 1
- (1996) Godon P., 1996, MNRAS, 282, 1107
- (1959) Godunov S. K., 1959, Math. Sb., 47, 271
- (1969) Goldreich P., Julian W. H., 1969, ApJ, 157, 869
- (2004) Goodwin S. P., Mestel, J., Mestel L., Wright G. A. E., 2004, MNRAS, 349, 213
- (1999) Gruzinov A., 1999, [arXiv:9902288]
- (2008) Gruzinov A., 2008, [arXiv:0802.1716]
- (1987) Harten A., Engquist B., Osher S., Chakravarthy S., 1987, J. Comput. Phys., 71, 231
- (2009) Kalapotharakos C., & Contopoulos I., 2009, 496, 495
- (2002) Komissarov S. S., 2002, MNRAS, 336, 759
- (2004) Komissarov S. S., 2004, MNRAS, 350, 427
- (2006) Komissarov S. S., 2006, MNRAS, 367, 19
- (2007) Komissarov S. S., Barkov M., Lyutikov M., 2007, MNRAS, 374, 415
- (2008) Li S., 2008, J. Comput. Phys., 227, 7368
- (2003) Lynden-Bell D., 2003, MNRAS, 341, 1360
- (2005) Lyubarsky Y. E., 2005, MNRAS, 358, 113
- (2006) Lyutikov M., 2006, MNRAS, 367, 1594
- (2003) Lyutikov M., Uzdensky D., 2003, ApJ, 589, 893
- (2006) McKinney J. C., 2006, MNRAS, 368, L30
- (2009) Narayan R., Li J., Tchekhovskoy A., 2009, ApJ, 697, 1681
- (2000) Priest E., Forbes T., 2000, Magnetic Reconnection. MHD Theory and Applications. Cambridge Univ. Press, Cambridge
- (1981) Roe P. L., 1981, J. Comput. Phys., 43, 357
- (1988) Shu C. W., Osher S., 1988, J. Comput. Phys., 77, 439
- (1997) Shu C. W., 1997, Technical Report, Essenentially Non-Oscillatory and Weighted Essentially Non-Oscillatory Schemes for Hyperbolic Conservation Laws, ICASE Rep. 97-65, NASA Langley Research Center, VA
- (2006) Spitkovsky A., 2006, ApJ, 648, L51
- (2006) Timokhin A. N., 2006, MNRAS, 368, 1055
- (1997) Toro E. F. 1997, Riemann Solvers and Numerical Methods for Fluid Dynamics, Springer-Verlag, Berlin
- (2000) Toth G., 2000, J. Comput. Phys., 161, 605
- (1979) Van Leer B., 1979, J. Comput. Phys., 32, 101
- (2006) Zhang W., MacFadyen A. I., 2006, ApJS, 164, 255

## Appendix A Mapping from face-centered field to volume-centered field

We assume the magnetic field at the cell faces has the following parabolic form,

(27) |

(28) |

where the superscript designates that the coefficients are at the cell faces. To match the above two equations at the cell’s four boundaries, the reconstructed polynomials interior to one cell must take the form

(29) |

(30) |

Matching the interior magnetic fields with the magnetic fields at cell faces, we have (Li 2008)

(31) |

(32) |

(33) |

(34) |

(35) |

(36) |

(37) |

(38) |

(39) |

(40) |

where the superscript L, R, T, B denote the values at the left, right, top, and bottom faces for a particular cell. After some manipulations, the volume-averaged cell-centered magnetic field is obtained by

(41) |

(42) |

The coefficients at the left and right cell faces can be calculated as follows,

(43) |

(44) |

(45) |

(46) |

(47) |

(48) |

The coefficients at the top and bottom cell faces can be obtained similarly. In our simulations the equations (A15) and (A16) are used to calculate the cell-centered magnetic fields from the face-centered magnetic fields.