text:
# An Unsplit, Cell-Centered Godunov Method for Ideal MHD

## 1 Introduction

In this paper we present a new Godunov method for the equations of multidimensional ideal magnetohydrodynamics (MHD). We give results from an implementation of the unsplit, second-order method of Colella  for these equations. The base scheme solves the ideal MHD equations using a second-order predictor-corrector formalism. To the base scheme we add three algorithmic components, whose effects upon accuracy and stability are measured. The first component is a MAC projection  step, which uses a Poisson solver to ensure that the cell-edge centered fields used to calculate fluxes are divergence-free to solver tolerance. This component, though commonplace in the context of incompressible Navier-Stokes simulation, is new to the MHD community. The second component is an approximate projection  that uses another solution of the Poisson equation to ensure that the cell-centered field is divergence-free to second-order. The last component is a filter  that also acts to suppress monopole sources in the cell-centered field.

The section that follows covers recent work and some of the schemes used for ideal MHD simulation. It also introduces methods for enforcing the divergence-free constraint, and introduces some issues surrounding multidimensional MHD. Section 3 introduces our basic algorithm and the extensions we have implemented. A suite of linear and nonlinear test problems will be used to determine which of our algorithmic extensions are best suited to each problem type. These tests and results are covered in Section 4. The overall purpose is to find one combination of these extensions that is well suited to all of the problems considered. This will be done through comparisons to published results and in some cases to an eight-wave MHD algorithm we have implemented. We find that a projection step is indeed required for accuracy and stability of the schemes. For the eight-wave scheme in particular, use of the MAC projection is essential to obtain correct MHD shock jumps.

## 2 Background

The study of numerical algorithms for magnetohydrodynamics simulations remains an active one, with no one method having become the standard. Two generic algorithms are the most widely used at present: the Method of Characteristics/Constrained Transport (MOC/CT; common in the astrophysics community)  and shock-capturing (Godunov) methods . Each has distinct benefits and drawbacks. Codes implementing the MOC/CT algorithm are relatively simple in design, and satisfy the divergence-free constraint to machine precision. However, the method of characteristics used by the ZEUS scheme, as outlined in , is by construction second-order on Alfvén and advective waves, but does not address the two compressive waves of ideal MHD. Moreover, Falle  found that ZEUS exhibits spurious rarefaction shocks in certain 1-D MHD shock tubes for a non-isothermal equation of state. Codes implementing the shock-capturing algorithm on the other hand, while more complex, give highly accurate results even for strong shocks. They suffer from the drawback that the divergence-free constraint is only satisfied to truncation error, which can be large in the region of large jumps. In order to treat this difficulty, a variety of techniques have come into use. One such is the hybrid CT/shock-capturing scheme , for which the constraint is satisfied by design like in the MOC/CT case. The cost to the accuracy of the underlying shock-capturing scheme is unclear. Another approach, originally due to Brackbill and Barnes  and implemented by workers such as Ryu et al , uses a divergence cleaning step on the cell-centered fields to enforce the constraint. In a third approach, Powell and co-workers  use a eight-wave reformulation of the ideal MHD equations originally due to Godunov . Tóth  (hereafter T00) implements all three types of schemes, among others, using them as the basis for a comparison on a variety of 1-D and 2-D tests. More recently, Dedner et al  compare several hyperbolic schemes with additional waves and divergence-damping terms on the T00 tests.

The 2-D tests of T00 serve to underscore the importance of using multidimensional problems in evaluating different algorithms, since it was mainly in this context that differences between them became apparent. This is to be expected, since errors due to non-solenoidal fields will generally only show up for problems in two or more dimensions. Many shock-capturing MHD schemes use an operator-split (or dimensionally-split) formalism to treat multidimensional ideal MHD. This means that, for each spatial dimension of the scheme, the one dimensional MHD equations are applied once. Unsplit schemes, which instead use the full multidimensional version of the equations, have been implemented and shown to give results equivalent to those of split methods for hydrodynamical problems . Unsplit shock-capturing schemes for multidimensional ideal MHD are relatively new, however. One of our main goals is to assess the efficacy of different approaches for enforcing the divergence-free constraint in one such unsplit scheme.

### 2.1 The Divergence-free Constraint

Some means must be employed to ensure that the field satisfies the divergence-free constraint, since this is only guaranteed to within truncation error in shock-capturing schemes. A dramatic example of even small errors in the divergence-free condition leading to instability is given in Section 4.1. There, small amplitude, non-propagating waves become unstable unless non-solenoidal fields are smoothed. This phenomenon can be explained in terms of a modified equation analysis, discussed in Appendix A. The possibility of incorrect field topologies, incorrect dynamics, and numerical instability motivate efforts to formulate and understand different means of supressing non-solenoidal errors in the magnetic field.

The eight-wave MHD algorithm, as implemented by Powell et al  and others, addresses the problem by adding additional terms corresponding to monopoles to the ideal MHD equations. The resultant equations are symmetrizable, so that they are Galilean invariant and transport $`\stackrel{}{B}`$ . The additional terms show up in two ways for shock-capturing schemes. Since they modify the 1-D MHD equations used for characteristic tracing to include an additional eighth wave that travels at the flow speed, monopoles will be advected along with the flow. Such monopoles could be carried out of the domain, or they might build up at a stagnation point. Secondly, the additional terms appear as source terms, making the system non-conservative if the divergence-free condition isn’t already satisfied.

Dedner et al  test a scheme that extends the eight-wave concept to damp and advect monopole sources, even at stagnation points in the flow. This is done through the magnetic analog of an artificial compressibility term, an approach that appeared earlier in the context of the Maxwell equations . To the extent that the computed auxiliary field remains continuous, the scheme will remain divergence free. Any monopole sources will be advected at the fastest speed allowed under the Courant condition, and damped as they are advected. This method is very useful on unstructured grids, where solving the Poisson equation in order to project out the solenoidal component of the field is difficult.

Divergence cleaning, or Hodge projection, in shock-capturing schemes can address the problem of non-solenoidal fields in two ways. The most widely discussed  involves projection of the cell-centered field onto the space of divergence-free fields. Projecting in this manner with a centered difference approximation to the divergence is consistent with the underlying cell-centered scheme. One is left with fields at the advanced time which are divergence-free to machine accuracy. Such a projection has been found to give correct field topologies in shock tube problems . On the other hand, choosing to eliminate the divergence as measured in one metric does not guarantee that unphysical effects are not entering into the dynamics. An illustrative example from incompressible flow  where an analogous constraint on the fluid velocity, $`\stackrel{}{v}=0`$, occurs shows that checkerboard modes in the velocity can cause instabilities when using a centered difference to approximate the divergence. Such modes must be damped by a suitably chosen filter in order to regain stability.

Another option for cleaning of non-solenoidal fields is to project the fields at cell-edges, which are then used to calculate fluxes. In a MAC discretization , the divergence-free constraint is enforced to within solver tolerance, which can be set as low as machine precision. We use a multigrid solver to solve the associated Poisson equation. The MAC projection has the advantage that it does not affect the conservation properties of the scheme.

### 2.2 Multidimensional MHD

The MHD equations in two or more dimensions are decidedly more complex to solve than in one dimension. In MHD simulations with variations along the x-axis alone, there is no change in the field along the x-axis. The divergence-free constraint is therefore trivially satisfied. Obtaining the solution to the Riemann problem upon which Godunov methods are based is also relatively straightforward in 1-D.

In multiple dimensions, we are solving the full equations of ideal MHD, which in conservation form are

$`_t\rho +(\rho \stackrel{}{u})`$ $`=`$ $`0`$ (1)
$`_t(\rho \stackrel{}{u})+\left[\rho \stackrel{}{u}\stackrel{}{u}+\left(P+{\displaystyle \frac{B^2}{8\pi }}\right)𝐈{\displaystyle \frac{1}{4\pi }}\stackrel{}{B}\stackrel{}{B}\right]`$ $`=`$ $`0`$ (2)
$`_t(\stackrel{}{B})+\left[\stackrel{}{u}\stackrel{}{B}\stackrel{}{B}\stackrel{}{u}\right]`$ $`=`$ $`0`$ (3)
$`_t(\rho E)+\left[\left(\rho E+P+{\displaystyle \frac{1}{8\pi }}B^2\right)\stackrel{}{u}{\displaystyle \frac{1}{4\pi }}(\stackrel{}{u}\stackrel{}{B})\stackrel{}{B}\right]`$ $`=`$ $`0,`$ (4)
subject to the constraint $`\stackrel{}{B}=0`$. Here $`\rho `$ is the mass density, $`\rho \stackrel{}{u}`$ the momentum density, $`\stackrel{}{B}`$ the magnetic field, and $`\rho E=\frac{1}{2}\rho |\stackrel{}{u}|^2+\frac{1}{8\pi }|\stackrel{}{B}|^2+\frac{1}{\gamma 1}P`$ the total energy density. The $`_t`$ notation denotes derivatives with respect to time.

## 3 Equations and Algorithms

In general for hyperbolic conservation laws, the conserved variables $`U`$ evolve according to $`_tU+\stackrel{}{F}(U)=0`$. Our code uses a second-order Godunov-type method for hyperbolic conservation laws . Such volume-average schemes follow the flux of conserved quantities such as momentum into and out of each cell comprising the computational domain. Quantities are stored at cell-centers, their value at this point being the volume-average over the entire cell. During the course of a timestep, the flux of the conserved quantities at each edge of each cell is computed. Differencing these fluxes gives the update of the conserved quantity to the next time (see Figure 1).

More specifically, our code is based on the unsplit MHD algorithm of Colella . It ostensibly tracks the evolution of the conserved density $`\rho `$, three components of momentum $`\rho \stackrel{}{u}`$, three components of magnetic field $`\stackrel{}{B}`$, and total energy density $`\rho E`$, in two spatial dimensions. These variables would then be evolved according to the ideal MHD equations (1-4). However, for simplicity and accuracy a set of primitive variables $`W`$ consisting of the density, velocity, field, and pressure are evolved in time:

$$W=[\rho uvwB_xB_yB_zP].$$
(5)
We will switch freely between $`U`$ and $`W`$ here depending on which is most convenient; the indices will indicate the centering of the variables. Thus $`i,j`$ will denote cell-centered states, $`i,j,\pm `$ forward- and backward-interpolations to cell-edges, and $`i+\frac{1}{2},j`$ and $`i,j+\frac{1}{2}`$ (Riemann) states at the edges.

The fundamental aspects of our scheme are as follows. Given the state $`W_{i,j}^n`$ at time $`n`$ and spatial coordinate $`(x_i,y_j)`$, we simultaneously interpolate in space and extrapolate in time to obtain the states $`W_{i,j,\pm ,x}^{n+{\scriptscriptstyle \frac{1}{2}}}`$ on the two x-boundaries of each cell $`(i,j)`$. The same is done for the states at the y-boundaries, $`W_{i,j,\pm ,y}^{n+{\scriptscriptstyle \frac{1}{2}}}`$. All this is done in the normal and transverse predictor steps of the algorithm. Next, the Riemann problem is solved using the $`\widehat{W}_{i,j,\pm ,x}^{n+{\scriptscriptstyle \frac{1}{2}}}`$ states. These Riemann states are used to calculate fluxes at each edge of a cell in the corrector step. The fluxes are then differenced to give the update to the next timestep, $`W_{i,j}^{n+1}`$. We go into more detail on these steps of our scheme in Sections 3.1 through 3.4.

The additional parts of our MHD scheme consist of additional terms in the normal predictor step to ensure correct multidimensional behavior, and several steps that address the divergence-free constraint. These are also discussed in Sections 3.1 through 3.4.

### 3.1 Normal Predictor

The predictor computes $`W_{i,j,\pm ,x}^{n+{\scriptscriptstyle \frac{1}{2}}}`$ and $`W_{i,j,\pm ,y}^{n+{\scriptscriptstyle \frac{1}{2}}}`$ using a Taylor series expansion:

$$W_{i,j,\pm ,x}^{n+{\scriptscriptstyle \frac{1}{2}}}=W_{i,j}^n\pm \frac{\mathrm{\Delta }x}{2}_xW\frac{\mathrm{\Delta }t}{2}𝐀^x_xW\frac{\mathrm{\Delta }t}{2}\frac{U}{W}_yF^y.$$
(6)
Here the matrix $`𝐀`$ is related to the flux $`F^x`$ of of the conserved variable $`U`$ by $`𝐀=(_UW)(_WF^x)`$. We give formulae for $`W_{i,j,\pm ,x}^{n+{\scriptscriptstyle \frac{1}{2}}}`$, those for $`W_{i,j,\pm ,y}^{n+{\scriptscriptstyle \frac{1}{2}}}`$ being similar. The first three terms on the right-hand side are computed in the normal predictor, and we label this intermediate result $`\widehat{W}_{i,j,\pm ,x}^{n+{\scriptscriptstyle \frac{1}{2}}}`$. We now separate out the evolution of the normal field $`B_n=B_x`$ through the following notation:

$$\widehat{W}=\left[\begin{array}{c}\stackrel{~}{W}\\ B_x\end{array}\right],𝐀=\left[\begin{array}{cc}\stackrel{~}{𝐀}& a_B\\ 0& 0\end{array}\right].$$
(7)
The matrix $`\stackrel{~}{𝐀}`$ corresponds to the usual 1-D MHD equations, with its seven characteristics: forward- and backward-propagating fast, slow, and Alfvén, plus the advective wave. Second-order accuracy is achieved in part through the use of characteristic analysis to calculate second-order accurate derivatives $`_xW`$ in the spatial interpolation. This characteristic interpolation is based on calculation of the eigenvalues $`\lambda _k`$ and left- and right-eigenvectors $`l_k`$ and $`r_k`$ of the matrix $`\stackrel{~}{𝐀}`$, giving the following expression for the interpolation of the $`\stackrel{~}{W}`$ variables to cell-edges:

$`\stackrel{~}{W}_{i,j,\pm ,x}^{n+{\scriptscriptstyle \frac{1}{2}}}=\stackrel{~}{W}_{i,j}^n+{\displaystyle \frac{1}{2}}{\displaystyle \underset{k:\lambda _k0}{}}\left(\pm 1{\displaystyle \frac{\mathrm{\Delta }t}{\mathrm{\Delta }x}}\lambda _k\right)\alpha _kr_k`$ (8)
$`\alpha _k=\{\begin{array}{cc}\mathrm{Min}(\alpha ^0,\alpha ^+,\alpha ^{})\hfill & \mathrm{if}\alpha ^+\alpha ^{}0\hfill \\ 0\hfill & \mathrm{otherwise}\hfill \end{array}`$ (11)
$`\alpha ^0=\frac{1}{2}l_k(W_{i+1,j}^nW_{i1,j}^n)`$ (12)
$`\alpha ^+=2l_k(W_{i+1,j}^nW_{i,j}^n)`$ (13)
$`\alpha ^{}=2l_k(W_{i,j}^nW_{i1,j}^n)`$ (14)
The $`\alpha _k`$ represent the strength of the kth wave in the interpolant. The sum over $`\lambda _k<0`$ would correspond to backward-propagating waves used in the interpolation to $`\stackrel{~}{W}_{i,j,,x}^{n+{\scriptscriptstyle \frac{1}{2}}}`$, and similarly for $`\lambda _k>0`$ and $`\stackrel{~}{W}_{i,j,+,x}^{n+{\scriptscriptstyle \frac{1}{2}}}`$.

A full accounting for all x-derivative terms in the 2-D MHD equations shows that $`a_B`$ is given by

$$a_B=[0,\frac{B_x}{4\pi \rho },\frac{B_y}{4\pi \rho },\frac{B_z}{4\pi \rho },v,w,\frac{\stackrel{}{u}\stackrel{}{B}}{4\pi }]^T.$$
(15)
These terms are essential to the second-order accuracy of the scheme, in particular on multidimensional problems such as waves not propagating along the coordinate axes. They are incorporated into our algorithm through a simple finite differencing of the normal derivative of the normal field. The terms are added to those already present due to the characteristics-based interpolation:

$$\stackrel{~}{W}_{i,j,\pm ,x}^{n+{\scriptscriptstyle \frac{1}{2}}}:=\stackrel{~}{W}_{i,j,\pm ,x}^{n+{\scriptscriptstyle \frac{1}{2}}}\frac{\mathrm{\Delta }t}{2}a_B(𝒟_x^0B_x^n)_{i,j}.$$
(16)
The $`(𝒟_x^0B_x^n)_{i,j}=((B_x)_{i+1,j}^n(B_x)_{i1,j}^n)/(2\mathrm{\Delta }x)`$ correction term is the centered-difference approximation to $`_xB_x`$. Note that this approximation to the derivative is not limited. The need for such correction terms in an unsplit scheme for multidimensional MHD was not addressed in C90. It was first noted by Stone  during an examination of the accuracy of an unsplit Godunov scheme on advected flux rings.

### 3.2 Transverse Predictor

The last term in the evolution equation (6) is included via the transverse predictor. The basic idea is to approximate transverse derivatives (in this case in the y-direction) using a 1-D Godunov method. We take the states calculated in the normal predictor, $`\widehat{W}_{i,j,\pm ,y}^{n+{\scriptscriptstyle \frac{1}{2}}}`$, and first use them to solve the Riemann problem at each y-boundary in the domain. The resulting Riemann states $`U_{i,j+{\scriptscriptstyle \frac{1}{2}}}^{n+{\scriptscriptstyle \frac{1}{2}}}`$ are subsequently used to calculate the fluxes needed for the last term in equation 6. In more formal terms,

$`U_{i,j+{\scriptscriptstyle \frac{1}{2}}}^{n+{\scriptscriptstyle \frac{1}{2}}}=(\widehat{W}_{i,j,+,y}^{n+{\scriptscriptstyle \frac{1}{2}}},\widehat{W}_{i,j+1,,y}^{n+{\scriptscriptstyle \frac{1}{2}}})`$ (17)
$`F_{i,j+{\scriptscriptstyle \frac{1}{2}}}^{y,n+{\scriptscriptstyle \frac{1}{2}}}=F^y\left(U_{i,j+{\scriptscriptstyle \frac{1}{2}}}^{n+{\scriptscriptstyle \frac{1}{2}}}\right).`$ (18)
Here $`(,)`$ denotes the Riemann problem solution using the two states on either side of an edge as input; see Section 3.4. It is then a straightforward matter to use these fluxes to calculate a finite-difference approximation to $`_yF^y`$, and thereby complete the calculation of the edge-centered states.

### 3.3 Corrector

The corrector first calculates fluxes at all cell-edges using another Riemann problem solve. At this stage, we have the first-order accurate approximation to the interpolated states $`W_{i,j,\pm ,}^{n+{\scriptscriptstyle \frac{1}{2}}}`$ from the predictor in hand. The Riemann solver takes these states and returns a single state for each cell-edge, $`U_{i+{\scriptscriptstyle \frac{1}{2}},j}^{n+{\scriptscriptstyle \frac{1}{2}}}`$ and $`U_{i,j+{\scriptscriptstyle \frac{1}{2}}}^{n+{\scriptscriptstyle \frac{1}{2}}}`$. For instance,

$$U_{i+{\scriptscriptstyle \frac{1}{2}},j}^{n+{\scriptscriptstyle \frac{1}{2}}}=(W_{i,j,+,x}^{n+{\scriptscriptstyle \frac{1}{2}}},W_{i+1,j,,x}^{n+{\scriptscriptstyle \frac{1}{2}}}).$$
(19)
The formula for $`U_{i,j+{\scriptscriptstyle \frac{1}{2}}}^{n+{\scriptscriptstyle \frac{1}{2}}}`$ is similar. See Section 3.4 for more details on the Riemann problem solution.

The states at cell-edges are not guaranteed to be divergence-free. We modify the C90 algorithm to enforce the divergence-free condition for the Riemann problem states. These Riemann states have a non-solenoidal component that we treat through a MAC projection, earlier used in the context of incompressible fluid computations . The edge-centered fields $`(\stackrel{}{B}^{})_{i+{\scriptscriptstyle \frac{1}{2}},j}^{n+{\scriptscriptstyle \frac{1}{2}}}`$ and $`(\stackrel{}{B}^{})_{i,j+{\scriptscriptstyle \frac{1}{2}}}^{n+{\scriptscriptstyle \frac{1}{2}}}`$ are used to calculate a cell-centered monopole charge density $`q_M=\stackrel{}{B}^{}`$, and a Poisson solver is in turn used to find the scalar field $`\varphi `$ implied by this monopole charge density distribution. The scalar field satisfies the following relations, in which $`𝒟_x^\pm `$ correspond to the forward- and backward-difference approximations to the derivative $`\frac{}{x}`$, and similarly for $`𝒟_y^\pm `$:

$`(q_M)_{i,j}`$ $`=`$ $`𝒟_x^{}(B_x^{})_{i+{\scriptscriptstyle \frac{1}{2}},j}+𝒟_y^{}(B_y^{})_{i,j+{\scriptscriptstyle \frac{1}{2}}}`$ (20)
$`\left[𝒟_x^+𝒟_x^{}+𝒟_y^+𝒟_y^{}\right]\varphi _{i,j}`$ $`=`$ $`(q_M)_{i,j}.`$ (21)
The correction to the field is calculated from $`\varphi `$ as follows:

$`(B_x)_{i+{\scriptscriptstyle \frac{1}{2}},j}`$ $`=`$ $`(B_x^{})_{i+{\scriptscriptstyle \frac{1}{2}},j}𝒟_x^+\varphi _{i,j}`$ (23)
$`(B_y)_{i+{\scriptscriptstyle \frac{1}{2}},j}`$ $`=`$ $`(B_y^{})_{i+{\scriptscriptstyle \frac{1}{2}},j}\frac{1}{2}\left[𝒟_y^0\varphi _{i+1,j}+𝒟_y^0\varphi _{i,j}\right]`$ (24)
$`(B_x)_{i,j+{\scriptscriptstyle \frac{1}{2}}}`$ $`=`$ $`(B_x^{})_{i,j+{\scriptscriptstyle \frac{1}{2}}}\frac{1}{2}\left[𝒟_x^0\varphi _{i,j+1}+𝒟_x^0\varphi _{i,j}\right]`$ (25)
$`(B_y)_{i,j+{\scriptscriptstyle \frac{1}{2}}}`$ $`=`$ $`(B_y^{})_{i,j+{\scriptscriptstyle \frac{1}{2}}}𝒟_y^+\varphi _{i,j}`$ (26)
With this correction to the magnetic field of the Riemann states, the $`L^1`$ norm of the MAC monopole density is reduced from its initial value by a user-settable multiplicative factor, in our case $`10^{12}`$.

The algorithm now proceeds to calculate the fluxes associated with the Riemann states. These fluxes are then differenced to give the update to the next time $`U_{i,j}^{n+1}`$:

$`F_{i+{\scriptscriptstyle \frac{1}{2}},j}^{x,n+{\scriptscriptstyle \frac{1}{2}}}=F^x\left(U_{i+{\scriptscriptstyle \frac{1}{2}},j}^{n+{\scriptscriptstyle \frac{1}{2}}}\right)`$ (27)
$`U_{i,j}^{n+1}=U_{i,j}^n\mathrm{\Delta }t𝒟_x^{}F_{i+{\scriptscriptstyle \frac{1}{2}},j}^{x,n+{\scriptscriptstyle \frac{1}{2}}}\mathrm{\Delta }t𝒟_y^{}F_{i,j+{\scriptscriptstyle \frac{1}{2}}}^{y,n+{\scriptscriptstyle \frac{1}{2}}}`$ (28)
After the update to $`t^{n+1}`$, we are left with a cell-centered field $`\stackrel{}{B}_{i,j}^{,n+1}`$ that is no longer divergence-free by a centered-difference divergence metric. To what extent, if any, this is a problem depends on the physical problem being considered, and will be addressed later.

For those cases where a reduction of the divergence is required, two algorithmic extensions have been implemented. The first follows from noting that it is desirable to have a diffusive term of the form (see ):

$$\frac{(\stackrel{}{B})}{t}=\eta ^2(\stackrel{}{B})$$
(29)
act on the divergence of $`\stackrel{}{B}`$. This may be rewritten to eliminate a spatial derivative by pulling out a divergence operator, giving

$$\frac{\stackrel{}{B}}{t}=\eta (\stackrel{}{B})$$
(30)
A simple, single-step filter may be derived  as a finite-difference approximation of equation 30:

$`B_x:=B_x+\eta \mathrm{\Delta }t(𝒟_x^+𝒟_x^{}B_x+𝒟_x^0𝒟_y^0B_y)`$ (31)
$`B_y:=B_y+\eta \mathrm{\Delta }t(𝒟_x^0𝒟_y^0B_x+𝒟_y^+𝒟_y^{}B_y)`$ (32)
The advantage of this choice of discretization of Equation 30 lies in its effective damping of checkerboard modes; see  for analysis and comparison of different filtering schemes.

In order to choose a value for $`\eta `$ we use Fourier stability analysis, giving a stability condition for the scheme in equations 31 and 32 of $`\mathrm{\Delta }t\frac{2(\mathrm{\Delta }x)^2}{5\eta }`$. Since $`\mathrm{\Delta }t`$ is set by the Courant condition, we are able to derive a condition on $`\eta `$, giving the maximum amount of diffusion of the monopole sources possible, given our timestep and grid spacing:

$$\eta =C\frac{(\mathrm{\Delta }x)^2}{\mathrm{\Delta }t}$$
(33)
with $`C\frac{2}{5}`$. A stronger condition, $`C\frac{1}{5}`$, will always damp monopole modes. This formulation will both decrease the cell-centered divergence and damp checkerboard modes. The filtering is always used when an approximate projection is performed.

The numerical effect of the filter can be modified through the parameter $`C`$. We found, through the nonlinear tests outlined in Section 4, that values in the range $`10^210^1`$ worked best. We note here that we chose to apply the filter to the conserved variables $`U^{n+1}`$, so that the changed magnetic field causes no change in the total energy. Any addition (subtraction) of magnetic energy shows up as a decrease (increase) in the internal energy.

The second means for treating non-solenoidal cell-centered fields is an approximate projection, in which the cell-centered divergence is used to calculate a monopole charge density $`q_M=𝒟_x^0B_x^{}+𝒟_y^0B_y^{}`$. Note that the multigrid Poisson solver uses the finite-difference operator $`𝒟^+𝒟^{}`$, not $`𝒟^0𝒟^0`$, for the Laplacian $`^2`$. The resulting solution is not discretely divergence-free, as it would be if we used the operator $`𝒟^0𝒟^0`$ , but instead second-order accurate. However, the extension of exact solvers to adaptive meshes is extremely complicated , whereas the approximate projection used here is straightforward to extend to AMR . We note moreover that use of a centered-difference measure of divergence leaves unchanged unphysical checkerboard modes . An extensive analysis of approximate projections is given in .

The Poisson solve yields the scalar field $`\varphi `$ satisfying equation 21. We difference $`\varphi `$ as follows to give the corrected field:

$`B_x^{n+1}=B_x^{,n+1}𝒟_x^0\varphi _{i,j}`$ (34)
$`B_y^{n+1}=B_y^{,n+1}𝒟_y^0\varphi _{i,j}`$ (35)
We find in our tests that an approximate projection does not have a significant impact on the quality of our results when used in conjunction with a MAC projection. Given the significant computational expense of using both, we advise against such a scheme.

### 3.4 Riemann Solver

We solve the Riemann problem for ideal MHD using a linearized solver . It employs characteristic analysis, like that of the normal predictor of Section 3.1, to solve for the states at cell edges. In this case, eigenvalues $`\lambda _k`$ and eigenvectors $`l_k`$ and $`r_k`$ are calculated using the arithmetic average of the states to the left and right of the edge. (Call this the base state.) For example, the solution to the Riemann problem at $`x_{i+{\scriptscriptstyle \frac{1}{2}}}`$, $`W_{i+{\scriptscriptstyle \frac{1}{2}}}^{RP}`$, is built from the appropriately directed jumps across waves as follows:

$$W^{RP}=\{\begin{array}{cc}\hfill \frac{1}{2}\left(W^R_{k,\lambda _k>0}\alpha _kr_k+W^L+_{k,\lambda _k<0}\alpha _kr_k\right),& \text{ if }|u|ϵ_v\hfill \\ \hfill W^L+_{k,\lambda _k<0}\alpha _kr_k,& \text{ if }u<ϵ_v\hfill \\ \hfill W^R_{k,\lambda _k>0}\alpha _kr_k,& \text{ if }u>ϵ_v\hfill \end{array}$$
(36)
Here $`|u|`$ is the advective velocity (along $`\widehat{e}_x`$ in this case) of the base state. The constant $`ϵ_v`$ is chosen to be $`10^{14}`$ times the largest characteristic speed $`\lambda _k`$. Note the averaged solution used for small advective speeds. This ensures that advective velocities on the order of machine precision, whose sign is random, do not cause asymmetries in the Riemann problem solution. Finally, note that in order to maintain consistency, the longitudinal component of the magnetic field in the solution (in this case $`B_x^{RP}`$) is assigned the value of the longitudinal magnetic field of the base state, $`\frac{1}{2}(B_x^L+B_x^R)`$.

For strongly nonlinear problems, it sometimes happens that the CFL condition is not sufficient to keep a scheme stable. Ordinarily, we use a simple minimum of timesteps implied by the one-dimensional CFL condition. If accelerations are large, velocities can grow such that in one timestep the pressure or density becomes negative. In order to dynamically adjust to such situations, we implement a scheme that checks for negative pressures or densities in the cell-centered states, $`W_{i,j}^{n+1}`$, after the corrector step. If one is encountered, that timestep is restarted with the CFL number lowered by a factor of two from its nominal value, down to a minimum of one-eighth of the initial CFL number. Once the code has proceeded for several timesteps without again encountering negative values, the CFL number is raised by a factor of two, eventually reaching its initial value. In practice, negative values are encountered only rarely, so that the average CFL number over the course of a run is close to the nominal value. Note that even if the CFL ramps quickly back up to the nominal value, the timestep is never allowed to increase by more than 10% in one iteration. In this way, the effect of an increase in CFL number is in fact spread over several iterations. Note that the average CFL number remains above 0.4 in all tests of Section 4.

## 4 Numerical Tests

In this section, we compare the behavior of the code on a variety of linear and nonlinear problems. Both seven- and eight-wave MHD codes were used. The code implementing the eight-wave MHD algorithm also uses a predictor-corrector formalism. The characteristic analysis performed in the predictor steps uses an eighth wave carrying changes in the normal field at the advection velocity, in addition to the seven waves of ideal MHD. As a result, in equation 7 of Section 3.1, the $`(2,2)`$ entry of $`𝐀`$ is equal to the advective speed $`u`$ and not zero. The terms needed for second-order accuracy in ideal MHD (equation 15) are accounted for in this case, so that $`a_B=0`$. The eight-wave algorithm also adds a source term,

$$S=[0,\frac{B_x}{4\pi },\frac{B_y}{4\pi },\frac{B_z}{4\pi },u,v,w,\frac{\stackrel{}{u}\stackrel{}{B}}{4\pi }]^T(\stackrel{}{B}),$$
(37)
to the right-hand side of the MHD equations 1-4. The source term is calculated in the transverse predictor and corrector steps. It is added to the updated states along with the differenced fluxes. We note, however, that performing a MAC projection guarantees that $`\stackrel{}{B}=0`$ for the edge-centered states, and as a result the source term calculated in the corrector is numerically zero in this case. This implies that only with a MAC projection is the scheme conservative and satisfies the jump relations. The source term in the transverse predictor, while generally small, is not numerically zero.

For both MHD implementations, different variations on the base algorithm were tested. In what follows, codes with conservative filtering are labeled with ’CF’. An ’AP’ denotes codes with an approximate projection; it is always accompanied by conservative filtering, in order to suppress checkerboard modes. Those codes utilizing a MAC projection are labeled with ’MAC’. So, for example, a MAC projection code with filter is labeled ’MAC+CF’.

### 4.1 Simple Linearized MHD Waves

We tested the performance of the code on all four varieties (advective, fast, slow, and Alfvén) of linearized MHD waves. These tests comprised waves propagating along x- and y-axes (wavenumbers $`\stackrel{}{n}=(1,0)`$ and $`\stackrel{}{n}=(0,1)`$), along with waves at slopes of 1:1 ($`\stackrel{}{n}=(1,1)`$) and 2:1 ($`\stackrel{}{n}=(2,1)`$). The simulation domain had length $`L=1`$ in both dimensions, and the boundaries in all cases were periodic. The Alfvén waves used $`\rho _0=1`$, $`\stackrel{}{u}_0=0`$, $`\stackrel{}{B}_0=B_0\widehat{b}=\sqrt{4\pi }\widehat{b}`$ with unit vector $`\widehat{b}=(\frac{1}{\sqrt{2}},\frac{1}{\sqrt{2}})`$, and $`P_0=1`$. The perturbation $`\delta W`$ is

$$\delta W=\left[\begin{array}{c}0\\ 0\\ 0\\ c_A\\ 0\\ 0\\ B_0\\ 0\end{array}\right]\delta _{\mathrm{pert}}\mathrm{sin}\left(\stackrel{}{k}\stackrel{}{x}\right),\text{with}\stackrel{}{k}=2\pi \stackrel{}{n},$$
(38)
where $`\stackrel{}{n}`$ is the aforementioned vector of integers chosen so as to be consistent with the periodic boundaries. The Alfvén speed $`c_A=B_0/\sqrt{4\pi \rho _0}=1`$.

Fast and slow wave expressions are somewhat more complicated. In this case, $`\widehat{b}`$ lies at 45 degrees from the unit wavevector $`\widehat{k}`$, all other aspects of the unperturbed state remaining unchanged from the Alfvén case. The perturbation is

$`\delta W=\left[\begin{array}{c}\rho _0\\ \frac{\sqrt{2}c_{F/S}^2\widehat{b}_ya^2\widehat{n}_y}{c_{F/S}}\\ \frac{a^2\widehat{n}_x\sqrt{2}c_{F/S}^2\widehat{b}_x}{c_{F/S}}\\ 0\\ \sqrt{2}B_0\frac{c_{F/S}^2a^2}{c_A^2}\widehat{n}_y\\ \sqrt{2}B_0\frac{c_{F/S}^2a^2}{c_A^2}\widehat{n}_x\\ \rho _0a^2\end{array}\right]\delta _{\mathrm{pert}}\mathrm{sin}\left(\stackrel{}{k}\stackrel{}{x}\right)`$ (46)
$`a^2={\displaystyle \frac{\gamma P_0}{\rho _0}}`$ (47)
$`c_{F/S}^2={\displaystyle \frac{1}{2}}\left(a^2+c_A^2\pm \sqrt{c_A^4+a^4}\right)\text{corresponding to fast/slow speeds}`$ (48)
$`\widehat{n}=(\widehat{n}_x,\widehat{n}_y)={\displaystyle \frac{\stackrel{}{n}}{\stackrel{}{n}}}`$ (49)
All tests gave results consistent with second order accuracy, although waves at 2:1 showed slightly smaller ($`1.8`$, versus 2.0) rates of convergence in some components. The modifications suggested by Stone were essential in obtaining second-order convergence for the 1:1 and 2:1 slope tests, increasing the convergence rates from first- to second-order for waves propagating along the diagonal. An example of the convergence rates for small amplitude fast waves is given in Table 1.

A more stringent test of the code is to advect the linearized waves so that their profile remains stationary. In these tests, the background state has a non-zero velocity equal to the wavespeed:

$`\text{Alfvén:}\stackrel{}{u}_0=c_A\widehat{n}`$ (50)
$`\text{Fast/Slow:}\stackrel{}{u}_0=c_{F/S}\widehat{n}`$ (51)
An analysis of the eigenstructure of MHD shows (Appendix A) that such waves should cause trouble for the seven-wave MHD codes that do not suppress non-solenoidal fields. Errors generated in this case are not advected away, causing difficulty in the case that they are not diffused or otherwise dealt with. We find that the seven-wave MHD algorithm without the application of any projection or filter is unstable for low amplitude ($`\delta _{\mathrm{pert}}=10^4`$) fast waves with $`\widehat{n}`$ at 2:1 slope. The instability starts as a high frequency oscillation in the field parallel to $`\stackrel{}{n}`$. It grows with time and spreads to the other components, causing a low order of convergence even at early times. Once the oscillations reach a certain level, the solution becomes unstable. The oscillations are almost completely absent, and the code stable, using a MAC projection. Adding filtering of the field further reduces the errors. A code including filtering but no MAC projection, and one with approximate projection and filtering both stabilize the scheme. We conclude that either projection or filtering is essential to stability for this class of problems.

### 4.2 Decay of Linearized MHD Waves

The equations of ideal MHD neglect the effects of viscosity and electrical resistivity. Numerical dissipation, however, can affect the solution in ways that mirror these physical effects. Following Ryu, Jones, and Frank  (hereafter RJF95), we measure the decay of Alfvén, fast, and slow waves, and use the implied physical resistivity as a measure of the numerical resistivity of our ideal MHD scheme.

We use exactly the same set-up as RJF95, with Alfvén, fast, and slow waves propagating at a 1:1 slope and wavelength of $`\sqrt{2}`$ times the length of one side of the computational domain. The decay of these waves was measured by fitting a decaying exponential to a measure of the wave strength,

$$\epsilon =\underset{i,j}{}l_k(\delta W)_{i,j},$$
(52)
where $`(\delta W)_{i,j}=W_{i,j}W_0`$ is the perturbation in the primitive variables. $`l_k`$ is the left-eigenvector associated with the mode in question, evaluated at $`W_0`$, the unperturbed state. Figure 2 shows the dependence of the decay rate on resolution for runs with number of cells per dimension $`N=`$16, 32, 64, 128, 256, and 512. Comparison with Figure 4 of RJF95 reveals a much smaller decay rate for our unsplit method. Furthermore, we find that the decay rate varies according to the power law $`\mathrm{\Gamma }N^3`$, not $`N^2`$. The former is consistent with second-order accuracy, since the truncation errors have the following form:

$$\tau =C_2(\mathrm{\Delta }x)^2\frac{^3U}{x^3}+C_3(\mathrm{\Delta }x)^3\frac{^4U}{x^4}+𝒪(\mathrm{\Delta }x^4).$$
(53)
We see here that the second-order error term is proportional to $`\frac{^3U}{x^3}`$, which is dispersive in nature. The third-order term is the first error term that is dissipative. As a result, we expect the dissipation to decrease as $`\mathrm{\Delta }x^3`$, or $`N^3`$. We have no good explanation for the $`𝒪(\mathrm{\Delta }x^2)`$ decay law observed in RJF95, aside from the possible effects of artificial viscosity.

### 4.3 MHD Shock Tube

The second test is the MHD shock tube problem from Ryu, Jones, and Frank . The solution consists of two fast shocks, one slow shock, one slow rarefaction, and a contact discontinuity. We have run the problem in two orientations: (1) with the shock velocity aligned with the x-axis of the computational domain (referred to below as 1-D), and (2) with this velocity inclined at a 2:1 slope. The latter configuration follows the 2-D shock tube test case from T00. It was chosen in order to test the multidimensional behavior of the code, meanwhile ensuring there were no serendipitous cancellations of errors, as might be the case in a $`45^{}`$ inclined case. In runs of a 1-D tube at $`R_{512}`$ (ie. 512 cells per linear dimension) on a domain of size $`L=1`$ to a time of $`t=0.08`$, we are able to reproduce the results given in Table VI of Dai & Woodward  (hereafter DW94) with errors of $`0.12\%`$. These results are independent of the type of filter used and whether an approximate projection was performed. Both the seven- and eight-wave codes give exactly the same result. When combined with an observed first-order convergence rate, they give us great confidence in the performance of both codes on 1-D shock tube problems.

The inclined version of this shock tube problem was run on grids of size $`2N\times N`$, with $`N=64`$, 128, and 256. Unlike in the aligned shock tube, in this case divergence cleaning is required. Since the natural boundary conditions for the problem would be tedious to implement in our multigrid solver, we chose instead to embed the physical domain within a larger one, four times its size. The data in the region outside of the original one is filled using the natural symmetry of the problem. The Poisson equation is then solved on this enlarged domain, using homogeneous Neumann boundary conditions. We have verified that the field in the original domain is clean to within the tolerances used previously.

We compared the results from the inclined shock tube with coarsened versions of a $`R_{4096}`$ 1-D shock tube results by first taking cuts of the data that included all cells lying along a line at 2:1 slope ($`\alpha =26.57^{}`$). Then, all velocities and fields in the cut were rotated by $`\alpha `$. A check using the initial conditions showed that such a cut of the inclined initial conditions matched perfectly with the 1-D shock tube ICs outside of the jump region. (In the jump region, the inclined ICs were slightly different due to the volume averaging performed in producing them. This leads the initial jump to be spread over two zones instead of one.) Note that we chose the grid spacing for the inclined runs to be a factor of $`\sqrt{5}`$ smaller than the 1-D grid spacing, so that the shock covers the same distance in physical space in a given time. The domain was therefore of length $`L=1`$ in the direction of shock propagation, and the final time $`t=0.08`$, in both the aligned and inclined shock runs.

Figure 3 shows results for the inclined shock tube overlaid on the (coarsened) aligned result. They compare a $`R_{256}`$ 1-D run with a $`N=256`$ 2-D run that has equivalent effective grid spacing. The ideal MHD codes shown give fractional errors with an $`L^1`$ norm of 3% or less, averaged over all components. Errors in the normal component of the field are relatively large, however. As illustrated here and by Figure 11 in T00, shock-capturing codes generally produce unphysical variation in the normal component of the magnetic field inside the shock transition layer in non-grid-aligned shock tube problems such as this. More generally, such codes exhibit non-monotonic behavior inside the shock structure, so that Riemann invariants are not exactly preserved there . Importantly, despite errors in $`B_n`$ inside the jump region, the jump relations are still satisfied outside of it. This is clearly illustrated in Figure 4, where we see the normal field for seven- and eight-wave schemes with filtering. Note in particular that while the eight-wave scheme converges to an answer, it does not give the correct shock structure.

Figures 6 and 7 show plots of the $`L^1`$ norm of the error versus resolution for the 2-D shock tube runs. We first note that neither the seven- nor eight-wave base code converges at first order. This is most evident in the errors in the transverse velocity and transverse field. The same plots indicate filtering alone is not sufficient to produce reasonable convergence rates. A MAC projection step added to either scheme produces much improved results (label ’MAC’ in the figures), and the further addition of a filtering step gives first-order convergence in all components (label ’MAC+CF’). The use of an approximate projection (labels ’CF+AP’ and ’MAC+CF+AP’) produces decreased error in the normal field, indicating a more accurate solution at given resolution. However, no improvement in the convergence rates or errors in the other five components is observed. When averaged over all components, the decrease in error observed with an approximate projection is still noticeable, though less so than indicated by considering the normal field alone.

While the results for seven- and eight-wave codes with a MAC projection step are very similar, this picture changes somewhat when the MAC projection is replaced with an approximate projection (label ’CF+AP’). The seven-wave code gives identical results whether or not a MAC projection is used with the approximate projection and filter. The eight-wave code still exhibits smaller errors when an approximate projection is utilized. However, convergence of the errors for the eight-wave code with approximate projection alone stalls at higher resolution. The reason for this behavior becomes more clear upon closer inspection of the normal field, plotted in Figure 5. The normal field converges to an answer that is 0.1% above the expected value. This error, while masked by the larger errors at the jumps at lower resolution, becomes significant at higher resolution and consequently hampers convergence. This error is observed to grow steadily over time.

### 4.4 Magnetized Flux Tube

This problem involves a high-field, low gas pressure region bounded on both sides by a high-gas pressure, zero-field region. The physical domain is of length $`L_x=L_y=1`$ on both sides. The base state has the entire 2-D domain in pressure balance. The boundaries between magnetized and unmagnetized regions are discontinuous and lie along $`x=\pm 0.2`$. In both regions, $`\rho =1`$, $`\stackrel{}{u}=0`$ and $`B_x=B_z=0`$ initially. In the magnetized region, $`x[0.2,0.2]`$, $`B_y=\sqrt{80\pi }`$ and $`P=1`$. Outside of this, $`B_y=0`$ and $`P=11`$. To the background state is added a sinusoidal perturbation upon the x-velocity whose amplitude is $`\delta _{\mathrm{pert}}=0.01`$ times the Alfvén speed, and covers the entire domain:

$$\delta u=\delta _{\mathrm{pert}}c_A\mathrm{sin}\left(2\pi y\right)\text{and}\delta \stackrel{}{B}=\delta \rho =\delta P=\delta v=\delta w=0$$
(54)
The Alfvén speed in this case is $`c_A=B_0/\sqrt{4\pi \rho }=\sqrt{20}`$. The strong discontinuity in the field is expected to cause problems for algorithms that do not suppress non-solenoidal fields. In particular, at the stagnation points where the perturbation velocity is zero, truncation errors leading to non-solenoidal fields can build up and cause numerical schemes to go unstable. This is in fact what we find for base schemes run without a projection step. Non-solenoidal fields build up at these stagnation points, causing the production of spurious velocities and triggering numerical instabilities.

Figure 8 shows the initial conditions for the flux tube problem. The perturbation plucks the field. After an initial transient at start-up, the perturbation develops into a standing wave in the magnetized region. The unmagnetized region sloshes back and forth from one side of the flux tube to the other, owing to the periodic boundary conditions. We expect the initial velocity perturbation to cause standing Alfvénic waves in the magnetized region and, because of the slightly different oscillation frequencies of the two regions, compressive waves at the boundaries of the tube. The change in the internal energy tracks the compressive waves. The problem was run to a time $`t=6.0`$, corresponding to about 60 Alfvén crossings of the short dimension of the tube. This was enough time for the problem to exhibit stable periodic behavior and subsequently evolve for many periods of this oscillation.

In this case we tested both seven- and eight-wave codes, but only with either a MAC projection and filter, or an approximate projection and filter. This is justified by the sub-par performance of the other variants on the inclined shock tube problem. These four versions of the code performed almost identically on this test, as evidenced by Figure 9. There, we plot global quantities such as kinetic and magnetic energy and $`L^1`$ norm of the monopole density. Note that the cell-centered divergence of the field is smaller for the approximate projection, though the dynamics remain the same.

### 4.5 Inclined Flux Tube

A second version of the flux tube problem, in which it is inclined at an angle of 45 degrees with respect to the original, is better at differentiating between the algorithms. It constitutes a strong test of the robustness and stability of the codes.

The size of the domain was $`L_x=L_y=\sqrt{2}`$, and both magnetized and unmagnetized regions have the same physical extent as in the aligned case. In both regions, we have $`\rho =1`$, $`\stackrel{}{u}=0`$ and $`B_z=0`$ initially. In the magnetized region, $`B_x=\sqrt{40\pi }`$, $`B_y=\sqrt{40\pi }`$, and $`P=1`$. In the unmagnetized region, $`P=11`$ and $`B_x=B_y=0`$. The perturbation is again applied to the entire domain, and has strength $`\delta _{\mathrm{pert}}=0.01`$:

$`\delta u=\delta _{\mathrm{pert}}{\displaystyle \frac{c_A}{\sqrt{2}}}\mathrm{sin}\left(2\pi {\displaystyle \frac{x+y}{\sqrt{2}}}\right)`$ (55)
$`\delta v=\delta _{\mathrm{pert}}{\displaystyle \frac{c_A}{\sqrt{2}}}\mathrm{sin}\left(2\pi {\displaystyle \frac{x+y}{\sqrt{2}}}\right)`$ (56)
$`\delta \stackrel{}{B}=\delta \rho =\delta P=\delta w=0`$ (57)
One effect of rotating the flux tube is the addition of a numerical perturbation, or gridding effect. The volume-averaging required for producing the initial conditions creates a region of intermediate pressure and field between the fully magnetized and unmagnetized regions. The total pressure in this intermediate region was kept the same as elsewhere in the domain. The physical extent, and presumably physical effects, of this transition region become smaller as the resolution is increased. The second effect of our rotation is that the physical domain was somewhat larger. As mentioned above, this was done in order to keep the physical width of the magnetized and unmagnetized regions the same as for the aligned flux tube. The result is a factor of $`\sqrt{2}`$ mismatch between the grid spacing in the aligned and inclined flux tube runs.

We expect many qualitative similarities between the aligned and inclined flux tube results. The characteristics of the start-up transients are somewhat different, due in part to the gridding effects. The slightly different grid spacing also means that slightly higher resolution is required to obtain comparable results. Indeed, these tests were run at resolutions ranging by factors of two from $`R_{128}`$ to $`R_{1024}`$, a factor of two higher than for the aligned case.

We again chose to run only those codes that showed close to first-order convergence for the inclined shock tube on this problem: seven- and eight-wave codes with either MAC projection and filter, or approximate projection and conservative filter. There was very little difference between the two seven-wave codes. The difference between seven- and eight-wave codes was more noticeable. The eight-wave codes both experienced large increases in kinetic energy, indicating that the magnetized tube was no longer in equilibrium. The size of this increase was largest ($`10^3\mathrm{KE}_{\mathrm{initial}}`$) for the eight-wave code with approximate projection, which also experienced substantial deviations from energy conservation. The size of these kinetic energy jumps decreased with resolution.

Figure 11 shows a resolution study of the inclined flux tube, with global quantities plotted out to $`t=1`$. (This is about 5 Alfvén crossing times; a shorter time was chosen in order to show more detail. All codes that were stable to this time were also stable out to $`t=6`$.) The plots include kinetic, internal, and magnetic energies, all scaled to their initial values. We also show the $`L^1`$ norm of cell-centered measure of the divergence, $`D^0\stackrel{}{B}`$, along with the magnetic energy change due to the MAC projection and filtering steps. First, note that the internal and magnetic energies appear to be converging, while convergence of the kinetic energy is not so clear. This result can be attributed to the large difference in the size of these quantities; the initial kinetic energy is $`10^5`$ times smaller than the initial internal energy, and $`10^4`$ times smaller than the initial magnetic energy. Thus, small errors in either of these quantities can appear as sizeable changes in the kinetic energy. We note in particular that increases in the kinetic energy seem to correlate with increases in the fraction of the magnetic energy subtracted in either the filter or the MAC projection.

The expected qualitative similarities between aligned and inclined flux tube results are present, such as oscillations in kinetic and magnetic energy of similar magnitudes. The size of the oscillations in the $`L^1`$ norm of $`v_t`$ and $`v_n`$ were also found to be very similar. However, substantial differences are also present. Comparing Figures 10 and 9, we note that the transient rise in kinetic energy early in the simulation begins later for the inclined tube than for the 1-D tube. In fact, the rise begins later and has smaller amplitude as resolution is increased, a fact that we attribute in part to the decrease in physical extent of the transition region.

As a final test of the nonlinear behavior of the code, we reproduce the first MHD rotor problem outlined in T00, earlier performed by Balsara and Spicer . This problem constitutes a central, high-density ($`\rho =10.0`$) region surrounded by a low-density medium. The central region is given a constant angular velocity. The entire domain is threaded by an initially constant field, $`\stackrel{}{B}=5.0\widehat{e}_x`$. The rotation winds the magnetic field, sending Alfvén waves propagating into the surrounding medium.

This test was run using both seven-wave codes (MAC+CF and AP+CF) and one eight-wave variant (MAC+CF). We see no negative pressures in any of the code variants, and obtain first-order convergence (errors decreasing as the grid spacing to the 0.95 power, as averaged over all variables) using the seven-wave codes, and the eight-wave code with a MAC projection. A reproduction of Figure 18 from T00 is given in Figure 12. Note the steeper gradient in Mach number in comparison to the T00 result; we do not see the peaks therein attributed to pressure undershoots.

## 5 Conclusions

We have presented an unsplit method for ideal MHD which, when combined with projection and filtering steps, shows no effects of non-solenoidal fields, while retaining the co-location of all physical quantites at cell-centers. The latter point is important because such a uniform centering makes it easier to extend the scheme to adaptive meshes, and has the advantage of using a well-understood Godunov method for time integration. These are in contrast to the Constrained Transport (CT) approach, where staggered grids add additional software complexity. Furthermore, it is not obvious how to discretize diffusion operators to include non-ideal MHD effects on staggered grids. Lastly, we add that the effective advection scheme for the magnetic field components is not of the standard type, and it is unclear what its properties are in the presence of under-resolved gradients.

The complexity and cost of the scheme relative to others is an important consideration. Our scheme requires three characteristic analysis steps plus twelve Riemann solution steps in 3-D. This is comparable to the six characteristic analysis and six Riemann solution steps required in a two-step Runge-Kutta scheme. The computational cost of projection is not insignificant, and adds additional software complexity. On a single grid and a single processor, the cost of the projection when computed using FFTs is a small fraction of the overall cost of the computation (less than 200 floating-point operations per cell in 3D). Current research on analysis-based solvers for Poisson’s equation  have the potential to make the cost of this remain less than that of the hyperbolic calculation, on both adaptive grids and on parallel processors.

The Hodge projection of fields either at cell-edges (MAC projection) or cell-centers (in our case via an approximate projection) was found to be essential to accuracy and stability. The MAC projection was essential for accuracy of the eight-wave code in the presence of discontinuities. The seven-wave code results on the inclined shock tube were less sensitive to the type of projection used. The inclined flux tube results painted essentially the same picture, though the seven-wave MAC projected result was somewhat more robust.

Use of a filter alone on linear problems reduces computational cost and does not affect accuracy. We find that both filter and projection are required on strongly nonlinear problems such as the inclined flux tube. Also essential to the accuracy of the scheme are the modifications suggested by Stone – linear waves not propagating along coordinate-axis directions are not second-order accurate without it.

We found that in determining the accuracy of our code, both measuring the rate of convergence on nonlinear problems and comparing absolute magnitude of the error were required. Significant differences were found between the convergence rates even of individual primitive variables. These data can be crucial in choosing between different algorithmic variants of a base code. We encourage future authors to incorporate convergence testing, taking care to calculate errors in each variable separately, of simple nonlinear problems such as the shock tubes run here into their test suites.

Overall the ideal MHD code with a projection, either MAC or approximate, and a filter performed best on the suite of tests presented. Filtering of the magnetic field was important to code stability in the nonlinear problems, and accuracy in the deficient wave problems.

The magnetized, perturbed flux tube constitutes a strong test of the stability of our schemes. The combination of a stationary discontinuity, low ratio of thermal to magnetic pressure (beta), and imposed perturbation caused significant problems for the base unsplit code without projection or filter. With a suitably chosen value of the filtering coefficient, the filter helped stability and worked to decrease the magnitude of a centered-difference measure of the divergence of the field.

## 6 Acknowledgements

The authors would like to thank Jim Stone for pointing out the missing multidimensional MHD terms.

The research of Robert K. Crockett was supported in part by A Division at Lawrence Livermore National Laboratory.

The research of Christopher F. McKee was supported in part by NSF grant AST-0098365.

Work at the Lawrence Berkeley National Laboratory is sponsored by the US Department of Energy Applied Mathematical Sciences program under contract DE-AC03-76SF00098 and by the NASA Earth and Space Sciences Computational Technologies Program under interagency agreement number S-44830-X.

The research of Richard I. Klein was supported in part by a NASA ATP grant NAG5-12042. Both Richard I. Klein and Robert T. Fisher are supported under the auspices of the US Department of Energy at the Lawrence Livermore National Laboratory under contract W-7405-ENG-48.

## Appendix A Appendix A: Divergence Constraints, Modified Equation Analysis, and Eigenvector Deficiencies

In this section, we will present a heuristic analysis of the effect of numerical errors in the divergence-free constraint on the stability of finite-difference methods for the ideal MHD equations. Our starting point will be the modified equation approach to analyzing the effect of truncation error on solution error. For any finite difference method, the modified equation is given by the original system of PDEs, with forcing terms given by the truncation error. In the present setting, the modified equation takes the following form:

$`_tU^{Mod}+\stackrel{}{F}(U^{Mod})=\tau _U(U^{Mod})`$ (58)
$`\stackrel{}{B}^{Mod}=\tau _D(U^{Mod}).`$ (59)
Here $`\tau _U`$ is the usual truncation error for the numerical method obtained from applying the difference operator to a solution to the differential equation evaluated on the grid. The equation for the evolution of $`\tau _D`$ is obtained by taking the divergence of the modified equation for the evolution of $`U^{Mod}`$,

$$_t\tau _D=\tau _B.$$
(60)
The truncation error forcing terms mimic the effect of numerical error on the computed solution. Specifically, we expect $`U^{Mod}`$, the solution to the modified equation, to satisfy $`U^{\mathrm{\Delta }x}U^{Mod}=𝒪(\mathrm{\Delta }x^{p+1})`$, where $`U^{\mathrm{\Delta }x}`$, the solution obtained from the $`p^{\mathrm{th}}`$-order scheme on a grid with spacing $`\mathrm{\Delta }x`$, satisfies $`U^{\mathrm{\Delta }x}U=𝒪(\mathrm{\Delta }x^p)`$.

In particular, for MHD, the effect of numerical error can be understood in terms of the truncation-error forcing in the modified equation causing the solution to violate the divergence-free constraint. Without that constraint being satisfied, the remaining ideal MHD equations can exhibit eigenvector deficiencies in the linearized-coefficient matrix $`𝐀`$, leading to anomalous loss of regularity and ill-posedness. In the numerical simulation, this translates into loss of accuracy and possibly instability of the underlying difference method.

To see this, we consider the case of a small-amplitude wave corresponding to one of the eigenmodes of $`𝐀`$ (see Equation 8):

$`W(\stackrel{}{x},t)=W_0+\alpha (x\lambda _kt)r_k`$ (61)
$`𝐀_0r_k=\lambda _kr_kr_k=(\stackrel{~}{r}_k,\mathrm{\hspace{0.17em}0})^T.`$ (62)
Then $`W(\stackrel{}{x},t)`$ satisfies the MHD equations up to terms of $`𝒪(\alpha ^2)`$.

Without loss of generality, we take the direction of propagation to be in the x-direction in 3-D. However, we allow our computational spatial grid to have an arbitrary orientation in space. In that case, the modified equation corresponding to our numerical solution to the PDE in primitive variables $`W`$ is given by (see Equation 7)

$`_tW^{Mod}+𝐀_0_xW^{Mod}=\tau _W`$ (63)
$`W^{Mod}=(\stackrel{~}{W}^{Mod},B_x^{Mod}).`$ (64)
If we define the new variable $`\alpha ^{Mod}=l_k(\stackrel{~}{W}^{Mod}\stackrel{~}{W}_0)`$, the modified equation dynamics can be reduced to the following system of two equations:

$`_t\alpha ^{Mod}+\lambda _k_x\alpha ^{Mod}+(l_ka_B)_xB_x^{Mod}=l_k\stackrel{~}{\tau }`$ (65)
$`_tB_x^{Mod}=\tau _B`$ (66)
In the system here, if one of the computational spatial grid coordinate axes is aligned with $`\widehat{x}`$, $`\tau _B0`$. However, if the direction of propagation is not aligned with one of the computational coordinates, then in general $`\tau _B0`$. If $`l_ka_B0`$ and $`\lambda =0`$, then the left-hand side of 65-66 is an example of a first-order system with an eigenvector deficiency. Such systems have an obvious loss of spatial regularity: $`\alpha ^{Mod}`$ grows like the derivative of $`\tau _B`$. This is in contrast to the behavior of well-posed hyperbolic systems, in which the solution has the same spatial regularity as the forcing. In other words, since regularity implies that there are as many derivatives in the solution as in the forcing term, discontinuous forcing of a hyperbolic system leaves the problem ill-posed. Discontinuity in $`\tau _B`$ implies that $`\alpha ^{Mod}`$ can grow without bound, owing to its dependence on the derivative of $`\tau _B`$.

In terms of a numerical method, we expect that the presence of such terms would lead to either an anomalous loss of accuracy or numerical instability. In the latter case, the forcing of $`\alpha ^{Mod}`$ in 65 takes the form of a finite difference operator applied to $`B_x^{Mod}`$, whose spatial variation is due entirely to $`\tau _B`$. If $`\tau _B`$ fails to be smooth, either because of lack of smoothness in the initial data or in the finite difference formulae (eg. limiters), such lack of smoothness is immediately amplified.

This discussion also provides an explanation for the behavior of the method described here. The use of the MAC projection and the filter does not eliminate the truncation error terms that lead to the eigenvector deficiency, but regularizes it by smoothing. For example, the application of the filter in the plane-wave example corresponded to adding a diffusion term to the equation for $`B_x^{Mod}`$,

$$_tB_x^{Mod}=\tau _B+\eta _x^2B_x^{Mod},$$
(67)
with $`\eta =𝒪(\mathrm{\Delta }x)`$. The use of the filter alone is sufficient to stabilize the small-amplitude plane wave solution in Section 4.1, and in that case leads to a second-order accurate result. The MAC projection performs a more drastic smoothing, but only on the intermediate form of $`B_x^{Mod}`$ used to compute $`_xB_x^{Mod}`$ in equation 65.


format: text

text:
# Hubble Space Telescope Observations of UV Oscillations in WZ Sagittae During the Decline from Outburst1footnote 11footnote 1 Based on observations made with the NASA/ESA Hubble Space Telescope, obtained at the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-26555.

## 1 Introduction

WZ Sge is the prototype of a class of cataclysmic variables that exhibit extreme characteristics: short orbital periods, extremely large dwarf nova outbursts, long outburst recurrence times, and low–mass companion stars. WZ Sge displays $``$7.5 mag outbursts separated by $``$2–3 decades and has a stellar mass ratio of 13–25 : 1 with a companion star that is less than 0.11 M and perhaps as low as 0.03 M (Steeghs et al. 2001, Ciardi et al. 1998). These extreme characteristics make WZ Sge an ideal proving ground for accretion disk theory.

Among the many challenges WZ Sge offers is the origin of the rapid 27.87 s and 28.96 s oscillations. First seen in the optical by Robinson, Nather & Patterson (1978), these oscillations have been detected in the UV (Welsh et al. 1997, Skidmore et al. 1999), infrared (Skidmore et al. 2002) and in the X-ray (Patterson et al. 1998), though weakly. The oscillations are complex, with large changes in amplitude, phase jitter, and transient signals at nearby periods. Sometimes both periodicities are present in the light curve. Robinson et al. (1978) proposed that the oscillations are due to non–radial $`g`$–mode pulsations of the white dwarf, a hypothesis supported by the simultaneous presence of the incommensurate 27.87 s and 28.96 s periodicities. Further support comes from analogy with GW Lib, a cataclysmic variable that almost certainly is a bona fide pulsating white dwarf (van Zyl et al. 2000, Szkody et al. 2002).

An alternative hypothesis, proposed by Patterson (1980), attributed the oscillations to a magnetized white dwarf channelling the accretion flow. The white dwarf’s rotation provides the stable clock driving the 27.87 s fundamental periodicity. The other complex, transient oscillations are due to reprocessing of radiation in the accretion disk (see Patterson et al. 1998 for more details). WZ Sge would then be a member of the DQ Herculis (or intermediate polar) class of cataclysmic variables. The detection of a 27.87 s periodicity in ASCA observations added much support to this magnetic accretor model. However, despite the success of this “oblique rotator” model in other cataclysmic variables, the presence of the simultaneous incommensurate periodicities in WZ Sge remains an obstacle for this interpretation.

In their investigation of quasi–periodic oscillations in cataclysmic variables, Warner & Woudt (2002) developed the “low inertia magnetic accretor” (LIMA) model and this yields a third possibility for the origin of the oscillations. The model contains two elements: (i) magnetically–controlled accretion onto a rapidly rotating belt on the white dwarf and (ii) a prograde travelling wave at the inner edge of the disk that produces a vertical thickening which acts as a site for reprocessing and can also occult part of the disk and/or white dwarf. In this scenario, the 27.87 s periodicity arises from magnetic accretion onto the white dwarf, but unlike the DQ Her model, the belt can slip on the white dwarf’s surface and hence the oscillation can vary in phase, amplitude and even period quite naturally. The 28.96 s signal comes from the reprocessing of the 27.87 s signal by the travelling wave whose orbital period is $``$740 s. While the model can satisfactorally explain the periodicities, and perhaps even the enigmatic $``$15 s quasi–periodic oscillation (QPO) discovered by Knigge et al. (2002), no other dwarf nova shows persistent QPOs in quiescence, requiring WZ Sge to be unique in this regard.

Understanding the nature of the periodicities in WZ Sge is important not only for its own sake, but because of the implications for accretion disk physics in general. Due to its extreme outburst characteristics, WZ Sge is believed to have a low mass–transfer rate from its companion star (e.g. Osaki 1996). Yet this alone is insufficient to allow accretion disk models to match the outburst size and timescale. One solution requires that the accretion disk have a “hole” at its center (Lasota et al. 1995). This inner disk region may be evacuated by (i) a magnetic field, as in the DQ Her stars where the white dwarf’s magnetic field truncates the inner disk (Warner et al. 1996; Lasota et al. 1999); or (ii) evaporation via a “siphon” into a hot corona (Meyer & Meyer–Hofmeister 1994). This coronal siphon mechanism has properties related to ADAF mechanisms in black-hole accreting systems (e.g. Mineshige et al. 1998, Meyer–Hofmeister & Meyer 2001). Another solution requires the quiescent viscosity in WZ Sge’s accretion disk to be very low: $`\alpha _{cold}\stackrel{<}{}0.001`$, compared to a more typical value of $`\alpha _{cold}0.03`$ in “normal” accretion disks (e.g. Smak 1993; Osaki 1996; Meyer–Hofmeister, Meyer & Liu 1998). If WZ Sge can be shown to be a DQ Her–like system, then the “universality” of the viscosity parameter is preserved and no seemingly ad hoc low value for $`\alpha _{cold}`$ is necessary. On the other hand, if WZ Sge is a white dwarf pulsator, then either the coronal siphon must be operating and we have a white dwarf analog of a black hole accretion disk, or an anomalously low viscosity is required. Either way, understanding the nature of the oscillations in WZ Sge can have an important impact on our understanding of accretion disks.

In 2001 July WZ Sge went into outburst, 23 years after its previous outburst (see Patterson et al. 2002 for extensive coverage of the 2001 outburst). Sion et al. (2003) obtained Hubble Space Telescope (HST) observations of WZ Sge on the decline from outburst and in this paper we present a time series analysis of those same data. These observations provide compelling evidence against a simple white dwarf non–radial $`g`$–mode pulsation interpretation of the oscillations.

## 2 The HST Observations

As part of the Director’s Discretionary time, WZ Sge was observed with the HST at four epochs: 2001 September 11, October 10, November 10, and December 11. Each observation lasted approximately 39 min, just under half of the 81.6 min binary orbital period. The UV light was probably dominated by emission from the heated white dwarf from October to December; in September the disk contribution may have been substantial and spectral decomposition into disk and white dwarf emission is much less certain (Sion et al. 2003). See figure 1 in Long et al. (2003) for the outburst light curve and the placement of the September, October and November HST observations.

The HST STIS FUV MAMA detector was used in TIME-TAG mode with a wavelength coverage of 1140 – 1735 Å (see Sion et al. 2003 for additional information regarding the observations). Data reduction was carried out with the standard STSDAS pipeline software. The data were cast into 2 s or 5 s time bins then binned over wavelength into a single UV light curve at each epoch.

The September observation took place just prior to the minimum before the final re–brightening (or “echo outburst”) of WZ Sge’s slow and erratic decline. The light curve shows large amplitude rapid flickering; a rms variation of 10.5% of the mean flux (21.8 mJy) was present. In contrast, the October, November and December light curves show much less flickering, with variations only slightly above the noise level. Figure 1 shows the four light curves, with respect to the start time of each observation. The mean UV flux dropped by a factor of 2.4 between the September and December observations.

## 3 Power Density Spectra

Figure 2 shows the power spectrum at each epoch, after prewhitening the light curves with a sliding boxcar average of width $``$146 s. Since the duration of the light curves were short and therefore the frequency resolution poor ($`4.292\times 10^4`$ Hz, or formally equivalent to a period resolution of only 0.36 s at 28.96 s), the power spectra were oversampled by a factor of 2, so only every other bin is independent. The power density is normalized such that a sinusoid of amplitude $`A`$ gives a power of $`\frac{1}{2}A^2`$. In all four cases, a spike consistent with the previously known 28.96 s periodicity is present. Henceforth we shall refer to this feature as the P29 signal. The period and amplitude of the P29 signal were determined via least–squares sine fitting. Table 1 lists the observed characteristics of the P29 signal.

The luminosity of the P29 signal dropped by a factor $``$4.5 between September and October, but thereafter remained fairly constant. In flux units, the amplitude changed from $``$0.42 mJy in September to $``$0.09 mJy in October–December, while the fractional amplitude dropped from $``$1.9% to $``$0.9%.

The P29 luminosity experienced a larger decrease than the mean UV luminosity, suggesting that P29 responds more quickly to changes in the global mass accretion rate than the mean system. One possible explanation is that the long cooling time of the white dwarf keeps the mean system luminosity well above its quiescent value; indeed in December WZ Sge was still much brighter in the UV than when in quiescence.

The 27.87 s periodicity was not detected, nor was the $``$15 s quasi–periodicity seen by Knigge et al. (2002) in HST UV observations earlier in the outburst. However, a broad feature near 18 s is present in the September observations which we identify as a quasi–periodic oscillation (QPO).

In Fig. 3 we show the September observations in more detail. The power spectrum was recomputed using only a linear detrending, thus the rapid rise in power at low-frequencies due to flickering remains present. Fig. 3 also shows the mean light curve and the amplitude (from the square root of the power) of the P29 and QPO periodicities. The P29 and QPO amplitudes were computed from short segments of the light curve, each 118 s in duration. Fifty such segments were used; they are not independent, having 61% overlap between them. The frequency resolution of such short light curves is very poor so the power spectrum was oversampled by a factor of 4. After removing the white noise background, the power was summed over the intervals 27.8–33.7 s for the P29 signal and 16.9–19.7 s for the QPO signal. The flickering red noise background remains present in these bandpasses, so these amplitudes are still somewhat biased high; the uncertainties in the amplitudes are roughly half the values of amplitudes themselves. The most noticeable features of the amplitudes are their rapid variations and lack of correlation with the mean light curve or each other. No periodicity in P29 or QPO amplitude is apparent. It is interesting to note that despite the P29 signal being present at 4 epochs spanning 3 months, the signal can vary from being quite strong to undetectable in a span of less than a few hundred seconds. Long–term stability is present while short–term stability is not. Although hard to quantify, the amplitude of the QPO is roughly 0.31 mJy or 1.4% of the mean system brightness. If only a single sinusoid is fit to the broad QPO signal, an amplitude of 0.27 mJy (1.3%) at a period of 17.95 s is measured. The relationship, if any, between this weak QPO and the much stronger $``$15 s and $``$6.5 s QPOs seen 20 days earlier in the outburst (Knigge at al. 2002) is unclear.

## 4 Discussion

### 4.1 The White Dwarf on the Decline from Outburst

Based on the same HST UV spectra as used in this work, Sion et al. (2003) measured the white dwarf’s temperature decline from 29,000 K in 2001 September to 18,000 K in December. Although some of the lines used in the analysis seem unusually narrow, the interpretation of a cooling white dwarf is fairly secure based on the success of the models in matching the observations, in particular the continuum and broad Ly $`\alpha `$ absorption. Further evidence of an elevated white dwarf temperature comes from the lack of strong quasi–molecular H2 or C I absorption lines that are prominent in quiescence. Nevertheless, the interpretation that the white dwarf completely dominates the UV spectrum at these epochs is not fully established. If we omit the September observations on the presumption that the accretion disk contributes significantly to this UV spectrum there still remains a 5000 K drop in white dwarf temperature between October and December. The uncertainties in these temperature estimates are roughly 10%, though the errors are largely systematic rather than statistical.

A similar change in white dwarf temperature was deduced from analysis of Far Ultraviolet Spectroscopic Explorer (FUSE) observation by Long et al. (2003), with temperatures dropping from $``$25,000 K to $``$22,600 K between September 29 and November 8. Part of the difference in estimated white dwarf temperature between the non–simultaneous HST and FUSE spectra comes from the uncertainty in the contribution of the disk to the total spectrum. Nevertheless the trend — a substantial cooling of the white dwarf — seems secure. Post–outburst cooling of the white dwarf has been well-documented in other dwarf novae, e.g., in VW Hyi (Gänsicke & Beuermann 1996, Sion et al. 1995), U Gem (Sion et al. 1998), and AL Com (Szkody et al. 1998).

Because of its critical importance in the interpretation to come, we repeat the statement in Godon et al. (2003) regarding the temperature of the white dwarf in WZ Sge. The temperature increase of the white dwarf, due to the enhanced accretion during outburst, is not confined to a thin outermost layer. Such a layer would very quickly radiate away heat gained during the outburst and therefore would remain at an elevated temperature only while substantial accretion is present. However, compressional heating takes place deeper in the star (though still in the outer layers) and takes months to cool, as is the case in WZ Sge (Sparks et al. 1993). So while irradiation from the boundary layer and inner disk heats the surface of the white dwarf, compression gets under the white dwarf’s skin (Sion 1995, 1999). Thus the dwarf nova outburst accretion episode, which dumps on the order of $`10^9\text{M}\text{}`$ of matter onto the white dwarf (Godon et al. 2003), produces a deep and significant heating of the white dwarf’s envelope.

### 4.2 The White Dwarf Non–radial $`g`$–mode Pulsation Hypothesis

Before we present the primary conclusion of our investigation, we briefly summarize the plausibility of the pulsation hypothesis.

The original, and in some ways still best, argument in favor of the white dwarf pulsation hypothesis is the observed simultaneous presence of the incommensurate P28 and P29 periodicities (Robinson et al. 1978). Incommensurate periodicities are common for the ZZ Ceti class of stars (isolated non–radial $`g`$–mode pulsating white dwarfs; see Kepler & Bradley 1995 for a good review of white dwarfs and pulsation). Also, the white dwarf in GW Lib is known to exhibit ZZ Ceti–type pulsations (van Zyl et al. 2000; Szkody et al. 2002). Additional support, albeit weak, comes from the UV spectrum of the oscillations that are consistent with a white dwarf photospheric origin (Skidmore et al. 1999, Welsh et al. 1997).

Yet the ZZ Ceti stars (and GW Lib) have much longer periods (hundreds of seconds) than WZ Sge. Furthermore, the quiescent temperature of the white dwarf in WZ Sge, $``$14,600 K (Cheng et al. 1997, Godon et al. 2003), exceeds the nominal upper temperature limit (the “blue edge”) of $``$13,500 K of the instability strip for ZZ Ceti stars. However, more massive white dwarfs have a hotter “blue edge” than the canonical 0.6 M white dwarf (Bradley & Winget 1994; Bergeron et al. 1995; Giovannini et al. 1998) and recent studies do in fact suggest WZ Sge has a massive white dwarf, exceeding $``$0.8 M and perhaps even as high as 1.2 M (Steeghs et al. 2001; Skidmore et al. 2000; Knigge et al. 2002; Long et al. 2003; Godon et al. 2003). By linearly extrapolating the blue edge temperature versus mass presented in Giovanninni et al. 1998, a $`g`$–mode pulsating white dwarf with a temperature of 14,600 K would require a white dwarf mass $`\stackrel{>}{}1.2\text{M}\text{}`$, barely within acceptable mass estimates. However, the extrapolation to values above $`1\text{M}\text{}`$ is highly uncertain so this crude mass limit should be taken as nothing more than a statement that $`g`$–mode pulsations cannot be ruled out by quiescent white dwarf temperature arguments. Furthermore, as was suggested by Szkody et al. (2002) for GW Lib, it could be that the temperature determined for WZ Sge is a global average, and a cooler region could be responsible for the pulsations — for GW Lib, a two–temperature white dwarf model fits substantially better then a one temperature model. Note that the thickness of the surface convective zone responsible for the pulsations in WZ Sge should be very thin compared to other ZZ Ceti stars since the timescale of the periodicity depends on the depth of the base of the convective zone (e.g. see Bergeron et al. 1995) — fast periodicities require a shallow convective zone.

Skidmore et al. (1999) extrapolated the observed temperature–period relation of ZZ Ceti stars using the data in Clemens (1993) and found the extrapolation falls tantalizingly close to the quiescent temperature and period of WZ Sge. This bolstered the viability of the white dwarf pulsation hypothesis, and mitigated the objection that the 28 s and 29 s periods are too short to be related to the ZZ Ceti phenomenon. We re-examined this extrapolation of figure 2 of Clemens (1993) and while the observed temperature–period relation for ZZ Ceti stars does project close to WZ Sge, the theoretical relation makes a dramatic drop to very fast periodicities at temperatures just higher than the observed ZZ Ceti star temperatures. From figure 3 of Clemens (1993), for a temperature of 14,600 K the thermal timescale at the bottom of the ionization zone where the pulsations are generated (which should be approximately equal to the observed periodicity) is on the order of only 0.01 s. There are large uncertainties in both the theory and observed temperatures, so this 0.01 s timescale should not be overinterpreted, but the conclusion remains thus: extrapolation of the observed temperature–period relation of ZZ Ceti stars to the WZ Sge regime cannot be trusted. However, we can be confident that if WZ Sge is a non–radial $`g`$–mode pulsator, it should follow, at least crudely, the trend of the temperature–pulse period relation for isolated white dwarf pulsators: hotter surface temperatures stars have shorter periods and smaller pulsation amplitudes (Clemens 1993).

Given the significant and deep heating of the white dwarf’s envelope discussed in §4.1, it would be astounding if the periods of any non-radial $`g`$–mode pulsations in WZ Sge did not change substantially as the white dwarf’s temperature changed. Yet in this work we have shown the continued presence of the 29 s periodicity while the white dwarf cools by over 10,000 K. During the outburst decline, WZ Sge no longer lies near the extrapolated empirical ZZ Ceti period–temperature relation and, more importantly, the observed period of the P29 oscillation is independent of the white dwarf’s temperature. These observations therefore rule out non–radial $`g`$–mode pulsations of the white dwarf as the sole origin of both the 28 s and 29 s pulsations.

### 4.3 Implications for Magnetic Accretor Hypotheses

Having claimed the observations of the P29 oscillation herald the end of the simple non–radial $`g`$–mode pulsation hypothesis, we are now left with a challenge: Can the observations be interpreted within the framework of the magnetically–channelled accretion flow models?

Patterson (1980) proposed the magnetic oblique–rotator model for WZ Sge (i.e. the DQ Her–type scenario) based on the stability of the 27.87 s oscillation and the ability to match some of the complex sideband structure in the power spectrum (the “satellite periods”) with simple integer multiples of the orbital frequency. In this model the rotation of the white dwarf provides a stable clock that produces the primary 27.87 s oscillation while the accretion disk provides (transient) structures that reprocess the primary signal. Cheng et al. (1997) measured a rapidly rotating white dwarf ($`V_{rot}\mathrm{sin}i=1200_{400}^{+300}`$ km s-1 in quiescence) supporting the magnetic rotator model, though this rotational velocity has not been confirmed with post–outburst observations. As Patterson et al. 1998 showed, this fast rotation requires a high–mass (0.8–1.1M) white dwarf if the period of rotation of the white dwarf is 28 s. Finally, the presence of a periodicity at 27.86 s in ASCA 2–6 keV light curves, along with a weaker signal in the 0.4–2 keV band (one half cycle out of phase), lends considerable support for the DQ Her interpretation (Patterson et al. 1998). Lasota et al. (1999) expand upon the DQ Her model and suggest that the rapidly spinning white dwarf actually ejects most of the matter transferred from the secondary star via a magnetic propeller. The P29 signal is caused by reprocessing of the P28 signal by blobs at the outer edge of the disk. A prediction of the model is that the white dwarf cannot be massive ($`\stackrel{>}{}1`$ M), potentially a serious problem if recent estimates suggesting a massive white dwarf are correct.

The LIMA model of Warner & Woudt (2002) is similar in spirit to the DQ Her model of Patterson (1980), but has added flexibility: the periodicities are not tied to the white dwarf’s rotation, but instead to a thin belt on the white dwarf. Since the belt contains much less inertia than the white dwarf, significant short timescale changes in rotation rate are possible. In this model the belt itself amplifies the white dwarf’s magnetic field, allowing a white dwarf with very low magnetic field strength to channel the accretion flow. One particularly attractive feature of the model is that the vertical thickening at the inner disk can also explain the “dips” seen in WZ Sge’s light curve. While the model seems very promising for the dwarf nova oscillations (DNOs) and QPOs seen in other cataclysmic variables, WZ Sge presents more of a challenge. For all other systems the DNOs are present only in outburst,333at least in the optical; Warner & Woudt (2002) point out an X-ray QPO seen during quiescence in HT Cas. while in WZ Sge the DNO signal, P28, is often present in quiescence and seems to have disappeared during the outburst.

#### 4.3.1 The 28s Periodicity

The P28 signal is not present in outburst (e.g. Patterson et al. 1981; Knigge et al. 2002; this work), and remained hidden for nearly 18 years after the 1978 outburst. In the standard DQ Her scenario, the lack of the P28 oscillation during outburst can be explained by the greatly enhanced mass accretion rate (see Patterson et al. 1998): the ram pressure of the accreting gas overwhelms the magnetic pressure and effectively crushes the magnetosphere back onto the surface of the white dwarf and accretion takes place along the equator, not at the magnetic poles. The magnetic field is too weak to channel the accretion flow during outburst, hence the photometric modulations vanish.

Even when deep in quiescence, observations do not always show the P28 signal; this is a puzzle. When present, the P28 signal maintains long–timescale clock–like stability, though short–term phase and amplitude changes are frequent. Inhomogeneities in the disk and or accretion flow must be such that the P28 signal is either not produced, hidden, or canceled for extended durations of time. Patterson et al. (1998) claim that sometimes when the P28 signal was absent other nearby periodicities that could be identified as sidebands of the P28 signal were observed, suggesting the P28 signal was present, just not detected. While the prolonged absence of the P28 signal after the 1978 outburst is not quite fully understood, a plausible reason exists, and the stability, sidelobes, and X–ray modulation can be interpreted, at least consistently, within the context of the DQ Her model.

#### 4.3.2 The 29s Periodicity

As seen here and many times previously, the P29 signal can exhibit rapid changes in amplitude and phase on timescales as short as minutes (e.g. Skidmore et al. 1997), yet the mean period is constant and on occasion the signal remains quite stable — perhaps even retaining phase coherence for up to a week (Patterson et al. 1998). The P29 signal exists throughout most of the dwarf nova eruption cycle: before, during and after the superoutburst — though perhaps not close to the peak of outburst (Knigge et al. 2002). The P29 signal does not correlate with the instantaneous luminosity of the system (and hence the inferred instantaneous $`\dot{M}`$ — see Fig. 3), though its amplitude is somewhat tied with the average brightness of the system on the decay from outburst. As we have shown in this work, P29 is de-coupled from the white dwarf temperature.

As Patterson (1998) stated, the P29 signal does not naturally appear as any known simple ratio or sum of periodicities in the system (orbital or spin). In the DQ Her model, the P29 signal is created by the reprocessing of the P28 signal by some structure in the disk: P29 is the “beat” between P28 and a blob orbiting in the disk with a period of $``$740 s, most likely at the inner edge of the disk. Similarly, in the LIMA model of Warner & Woudt (2002), the P28 signal corresponds to a DNO while the P29 signal corresponds to a QPO, created by the reprocessing of the DNO signal. While a periodicity of 741 s has been detected in WZ Sge’s light curve (Warner & Woudt 2002), it is very rare, and the DNO driving signal (i.e. P28) was not present at that time.

In the DQ Her model, the P28 period remains constant during outburst since it is maintained by the rotation of the entire white dwarf. In contrast, on the rise to outburst the period of the P29 signal should increase as the radius of the reprocessing inner disk shrinks, and the signal disappear entirely when the disk reaches the white dwarf. As the disk hole (i.e. magnetosphere) re–opens during the outburst decline, the P29 oscillation should re–appear and decrease back to its quiescent period.

In the LIMA model, the DNO period should decrease as more and more angular momentum is accreted into the belt during the rise to outburst. The P29 (QPO) signal would also experience a period change, for two reasons: it is being driven by the faster DNO and the inner radius of the disk where reprocessing is occurring is shrinking. On the decline from outburst, the periods should increase. A strong period–luminosity relationship through the outburst should exist, as is the case for DNOs of other systems (see Woudt & Warner 2002 and references therein).

In both magnetic accretor scenarios, the P29 signal should vary during the outburst. Yet as we have shown here, this is not the case in WZ Sge: the P29 signal remains constant in period on the decline from outburst, and with the same period it has in quiescence. In fact, just the presence of the P29 signal in the absence of the P28 signal during the outburst decay presents a significant problem for the magnetic accretor models. In the DQ Her model, if the magnetosphere is crushed onto the white dwarf thereby quenching the P28 signal, the inner edge of the disk is now coincident with the white dwarf’s surface. So not only is there no inner disk region to reprocess any modulation, there is no modulation to reprocess. And yet the P29 signal is quite strong, stronger in fact than in quiescence, at least in the 2001 September observations. As in the DQ Her model, the LIMA model requires the presence of the DNO at 28 s to drive the QPO at 29 s. The complete absence of the P28 signal in the observations presented here is vexing. We clearly do not understand the mechanism producing the P29 signal.

The continued presence of the P29 signal during the outburst decay, as the disk reforms and the mass accretion rate and white dwarf temperature dramatically change, suggests that a clock of exceptional (long–term) stability is involved. The most immediate candidate is the rotation of the white dwarf, but how this could work is not at all apparent. We conclude that while the rotating magnetic accretor models remain viable, they are sorely incomplete.

## 5 Summary

We present ultraviolet light curves of the dwarf nova WZ Sge as it declined from its 2001 July superoutburst. The observations were made with the HST STIS and were obtained on 2001 September, October, November and December. The well–known 28.96 s periodicity was present at all four epochs, but the 27.87 s periodicity was absent. We also note the appearance of a QPO with period $``$18 s in the September light curve.

Between 2001 September and December the temperature of the white dwarf dropped from $``$30,000 K to $``$18,000 K (Sion et al. 2003). We argue that by analogy with the ZZ Ceti stars, the white dwarf non–radial $`g`$–mode pulsation model for the origin of the 27.87 s and 28.96 s periodicities (Robinson et al. 1976) demands a large change in the pulsation period as the white dwarf’s temperature changes. However, these observations do not support this prediction: the 29 s oscillation period is observed to be independent of the white dwarf’s temperature. We consider this the coup de grace of the non–radial $`g`$–mode pulsation hypothesis as the sole origin of both the 28 s and 29 s pulsations.

In light of these observations, the alternative hypotheses for the origin of the periodicities, the venerable DQ Her magnetic rotator model (Patterson 1978) and the LIMA model of Warner & Woudt (2002), do not fare particularly well either. In these models the 28 s periodicity is related to rotation of either the entire white dwarf or an accretion belt on the white dwarf. The 28 s modulation is ultimately responsible for driving all the other periodicities including the 29 s modulation. The presence of a relatively strong 29 s periodicity in the total absence of the 28 s periodicity is a problem for the models. While the models remain viable, they must be considered incomplete. Finally, we remark that while non–radial $`g`$–mode pulsations cannot be responsible for the oscillations in WZ Sge, this does not necessarily imply that other modes of white dwarf pulsation are ruled out. In particular, the periods of non–radial $`r`$–modes (Papaloizou & Pringle 1978) can be near those of WZ Sge and are not very temperature sensitive so they may play a role in WZ Sge’s puzzling behavior.

This work was supported by NASA through grants GO-09304 and GO-09459 from the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA Contract NAS 5-26555. Support was also provided in part by NSF grant 99-01195, NASA ADP grant NAG5-8388 (EMS), and from NASA through STScI grant HST-GO-08156 (WFW). BTG acknowledges support from a PPARC Advanced Fellowship.