particular, this perturbation takes the form of a spher- ical, Gaussian shell of radius r$_{o}$ , width δ , and ampli- tude A added to each component of ⃗ B . We expect the divergence cleaning to propagate this perturbation as a damped wave, and we therefore plot the scaled quantity x 2 [ ⃗ ∇ · ⃗ B ] in Fig. 1. As can be seen by comparing the cases of no cleaning ( top ) and with cleaning ( bottom ), the divergence propagates with damping through the refine- ment boundaries across the grid. As typical with codes dealing with linearly degenerate hyperbolic systems, like those those employed in numeri- cal relativity, a dissipation operator is applied to the met- ric variables. This operation uses a high-order derivative to serve as a low-pass filter and does not affect the accu- racy of the simulation. We have found it useful for keep- ing things smooth. Although this operation is not applied to the fluid variables, we have found it quite important in keeping the magnetic field components smooth. The magnetic field evolution is coupled tightly with that of the velocity, and any nonsmoothness which appears in the velocity can easily affect the smoothness of the mag- netic field. The addition of dissipation to the magnetic field and divergence cleaning field helps control the be- havior of the magnetic field. At the boundaries of the domain, simple outflow boundary conditions are applied to the fluid variables. This is accomplished by copying the values of the conser- vative variables near the boundary outward. Most of the gravitational variables are treated using Sommerfeld-like boundary conditions of the form [46] where η$_{ab}$ is just the Minkowski metric. The rest, which are not so crucial, are set either by maximally dissipa- tive [2] or constraint preserving boundary conditions [31]. A relevant issue when considering the collapse of stars is the formation of a black hole. To deal with such situa- tions, we adopt black hole excision where we dynamically monitor for the appearance of trapped surfaces (which lie inside an event horizon if cosmic censorship holds) and excise cubical region(s) from the computational domain. As discussed in [47, 48], this excision introduces inner boundaries which are of “outflow” type and so no bound- ary condition is required there. However, for a more ro- bust handling of the fluid, we also allow for a modification of the fluid equations inside the trapped surface [7]. The MHD equations are written in balance law form which we modify to include a damping term near the black hole Here the function f ( r ) decreases smoothly with r , from a given value at the excision region to zero at the outermost Generally, for stable stars one sees this quantity oscillate from inherent numerical perturbations about a stable so- lution. For unstable solutions, one expects this quantity to change in a significant way. Similarly, we compute the relative change in baryon mass as a function of time as It is also standard practice to display the maximum of the (primitive) rest mass density, and we compute the fractional change in time as These parameters are computed with respect to the con- servative variable D . in terms of the moment of inertia tensor The initial configuration for the stars is axisymmetric, and we therefore want to be able to measure any change to this structure. For this purpose, we monitor certain distortion parameters [50] defined as extracted at spherical surfaces Σ$_{i}$ located in the wave- zone, far from the sources. We also consider possible corrections required to deal with gauge ambiguities, as discussed in [49]. We refer the reader to that paper for details on the adopted tetrad and required corrections. Finally, gravitational radiation is calculated via the evaluation of the Newman-Penrose scalar Ψ$_{4}$ which is computed by contracting the Weyl tensor, C$_{abcd}$ , with a suitably defined null tetrad { ℓ, n, m, ¯ m } trapped surface (OTS) found, and is zero for r ≥ r$_{OTS}$ , so that the exterior of the BH is causally disconnected from the effect of this extra term. U$_{0}$ is set to zero or to the value of the atmosphere if the corresponding field has one. The coefficient (Δ x ) p ensures that the damping term converges to zero as the grid spacing Δ x is reduced. As long as one chooses p greater than or equal to the order of convergence of the code, this term will not modify the convergence rate. We typically adopt a value of p = 4. 6 Ψ$_{4}$ = C$_{abcd}$n a ¯ m $^{b}$nc ¯ m d (41) ( ∂$_{t}$ + ∂$_{r}$ + 1 r ) ( g$_{ab}$ − η$_{ab}$ ) = 0 , (38) ˙ U + F ( U ) ′ = S, (39) ˙ U + F ( U ) ′ = S − f ( r )(Δ x ) $^{p}$( U − U$_{0}$ ) . (40) Δ M$_{baryon}$ ≡ M$_{baryon}$ ( t ) − M$_{baryon}$ (0) M$_{baryon}$ (0) (47) Δ ρ ≡ max | ρ$_{0}$ ( t ) | − max | ρ$_{0}$ (0) | max | ρ$_{0}$ (0) | . (46) I jk ≡ ∫ D x $^{j}$xk d $^{3}$x. (45) η$_{+}$ = I xx − I yy I xx + I yy (42) η × = 2 I xy I xx + I yy (43) η = √ η 2 + + η 2 × (44) IV. DIAGNOSTIC QUANTITIES