{'fraction_non_alphanumeric': 0.027757153463932833, 'fraction_numerical': 0.011232318617092162, 'mean_word_length': 4.655469422911283, 'pattern_counts': {'":': 0, '<': 8, '<?xml version=': 0, '>': 8, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Controlling the Growth of Constraints in Hyperbolic Evolution Systems

Motivated by the need to control the exponential growth of constraint violations in numerical solutions of the Einstein evolution equations, two methods are studied here for controlling this growth in general hyperbolic evolution systems. The first method adjusts the evolution equations dynamically, by adding multiples of the constraints, in a way designed to minimize this growth. The second method imposes special constraint preserving boundary conditions on the incoming components of the dynamical fields. The efficacy of these methods is tested by using them to control the growth of constraints in fully dynamical 3D numerical solutions of a particular representation of the Maxwell equations that is subject to constraint violations. The constraint preserving boundary conditions are found to be much more effective than active constraint control in the case of this Maxwell system.


I. INTRODUCTION
Perhaps the most critical problem facing the numerical relativity community today is the exponential growth of constraints in evolutions of black hole spacetimes. The evolution equations guarantee that constraints that are satisfied exactly on a spacelike surface will be satisfied throughout the domain of dependence of that surface. However, this result does not guarantee that small initial violations of the constraints will remain small, or that constraint violations will not be injected into the computational domain through timelike boundaries. Experience has shown that constraint violations tend to grow exponentially in the numerical evolution of black hole spacetimes (e.g., [1,2,3]). These constraint violating instabilities have been shown to be numerically convergent and thus represent unstable solutions to the partial differential equations. At present these instabilities are the limiting factor that prevents these numerical simulations from running for the needed time with the required accuracy.
A variety of approaches have been explored in a number of attempts to control the growth of the constraints: 1) Fully constrained evolution, in which the constraint equations are re-solved periodically (e.g. at each time step) have been used with great success in spherically symmetric and axisymmetric problems [4,5,6,7,8,9,10]. These methods have not gained widespread use in 3D simulations, however, due in part to the high cost of solving the elliptic constraint equations. Difficult questions also remain unresolved for this method about the proper boundary conditions to impose on the constraint equations at black hole excision boundaries. With the development of more efficient elliptic solvers and the absence of a better alternative however, fully constrained methods are starting to be developed and tested in 3D now as well [11,12,13].
2) Auxiliary dynamical fields have been introduced into the system whose evolution equations are designed to drive the constraints toward zero [14]. This technique has the disadvantage that it requires the size of the dynamical system to be significantly expanded. It has not been tested extensively, but the first numerical results were not uniformly successful [15,16].
3) More sophisticated boundary conditions have been introduced whose purpose is to control the influx of constraint violation through the timelike boundaries of the computational domain [17,18,19,20,21,22,23,24,25,26]. This approach seems very promising, although the current methods may not be fully compatible with the physical requirement that waves pass through the boundaries without reflection. Further these boundary condition methods may not completely solve the constraint violating instability problem in systems like the Einstein evolution equations, where constraint violations are driven both by bulk and by boundary terms in the equations. But this technique can (as we will demonstrate below) significantly improve the influx of constraint violations through the timelike boundaries of the computational domain. 4) Dynamically changing the evolution equations, through the addition of terms proportional to the constraints, has been proposed as a way to minimize constraint growth. This method (developed by Manuel Tiglio and his collaborators [27,28]) has had some success in controlling the growth of constraints in simple numerical solutions of the Einstein evolution equations. We find that this technique when used in combination with standard boundary conditions is not effective however in controlling the influx of constraint violations through the boundaries of the computational domain in fully dynamical situations.
In this paper we explore in some detail two of these methods for controlling the growth of constraint violations in hyperbolic evolution systems. First, we develop a refined version of the dynamical constraint control method being used by Tiglio and collaborators [27,28]. In particular we introduce a more natural norm on the constraints, which has the property that its evolution can be predicted numerically with greater accuracy. We expect that dynamical constraint control based on this new constraint norm should be more stable and robust than the current method. And second, we explore the use of constraint preserving boundary conditions. In this method (explored previously by Calabrese and collaborators [19,25]) the constraints are decomposed into characteristic ingoing and outgoing fields of the constraint evolution equations. Setting the incoming components of the constraint fields to zero provides boundary conditions for the incoming parts of the dynamical fields of the principal evolution system. We test both of these methods by applying them to a non-trivial hyperbolic evolution system (a particular representation of the Maxwell system [29,30]) that is analogous to, but much simpler than, the Einstein evolution system. Our tests-using fully dynamical time dependent solutions on domains with open boundaries-reveal that the constraint preserving boundary conditions are much more effective than active constraint control for this Maxwell system. Some features of this system are rather special, and it is possible that in more generic systems (like the Einstein equations) the active constraint control method may be complimentary to the constraint preserving boundary condition method.
We define and review in Sec. II the particular form of the Maxwell evolution system [29,30] that we use to illustrate and test the constraint control methods studied here. We refer to this system as the "fat" Maxwell system since it replaces the usual representation of the Maxwell system, which has six independent field components, with a representation having twelve. We also present in Sec. II the decomposition of the dynamical fields used in this fat Maxwell system into characteristic parts. In Sec. III we develop the equations needed to perform active constraint control, in particular on the fat Maxwell evolution system. We determine the constraint evolution equations for this system, and derive an improved norm on the constraint fields. We show how the evolution of this new constraint norm should generically be more accurately determined (and hence should provide better constraint control) in numerical solutions. In Sec. IV we develop the particular form of constraint preserving boundary conditions studied here. We present the decomposition of the constraint fields into characteristic parts, and show how these can be used to provide boundary conditions for the principal evolution system. Finally in Sec. V we use these methods to control the growth of constraints in fully dynamical 3D numerical evolutions of the fat Maxwell evolution system. We note that both the active constraint control mechanism and the constraint preserving boundary conditions developed here are applicable to rather general hyperbolic evolution systems. We focus our discussion on the fat Maxwell system in order to make the analysis less abstract, and to provide a simple system on which to perform numerical experiments.

II. FAT MAXWELL EVOLUTION SYSTEM
Our primary interest here is to understand how to control the growth of constraints in hyperbolic evolution systems. We will focus our attention on quasi-linear systems of the form, where u α represents the dynamical fields, and A k α β and F α may depend on u α but not its derivatives. We assume that the evolution system described in Eq. (1) is also subject to a set of constraints, c A = 0, which we assume have the general form where K A k β and L A may depend on the dynamical fields u α but not their derivatives. We assume that these constraints are preserved as a consequence of the evolution equations. In particular we assume that the constraints satisfy an evolution equation of the form where A k B D may depend on the dynamical fields u α , while F B D may depend on u α and its spatial derivatives ∂ k u α . When this constraint evolution system is hyperbolic the constraints will remain satisfied within the domain of dependence of the initial surface if they are satisfied initially. We note that multiples of the constraints of the form given in Eq.
(2) may be added to the principal evolution system Eq. (1) without changing the physical (constraint satisfying) solutions of the system or the basic structure of Eq. (1). Systems with this general form include most of the evolution equations of interest in mathematical physics, including for example the Einstein evolution equations, the Maxwell equations, the incompressible fluid equations, etc. In order to explore and test some of the ideas for controlling the growth of constraints in these hyperbolic evolution systems, we adopt a simple example system on which to perform our analysis and to carry out numerical tests. We have chosen to use a form of the vacuum Maxwell evolution equations (introduced independently by Kidder [29] and Reula [30]) that fits nicely into this framework, and that admits constraint violations if nothing is done to control them. The dynamical fields in this formulation are a co-vector that represents the electric field E i , and a second rank tensor D ij that represents the gradient of the spatial parts of the vector potential (i.e. D ij = ∂ i A j , although we impose the relationship between D ij and the vector potential only implicitly as a constraint on this system). We refer to this as the fat Maxwell system, since the usual representation of the Maxwell equations with six dynamical field components is replaced with this larger representation that has twelve. The evolution equations for this system are, where g ab is the Euclidean metric with inverse g ab , and ∂ a is the covariant derivative compatible with this metric, (i.e. just partial derivatives in Cartesian coordinates). The scalar potential φ is a gauge quantity assumed here to be a given function of space and time. This system has the same general form as Eq. (1) with u α = {E i , D ij }. In order to represent the vacuum (i.e. charge and current free) Maxwell system these equations are also subject to the constraints, These constraints have the same general form as those described in Eq. (2) with c A = {C, C ijk }. The second of these constraints is the integrability condition that guarantees that D ij is the gradient of a vector potential. As mentioned above we are free to add multiples of the constraints to the evolution system: where γ 1 and γ 2 are constants. The addition of these constraint terms leaves the physical (constraint preserving) solutions to the system unchanged, and also leaves the system with the same basic structure as Eq. (1). For hyperbolic evolution systems, such as those in Eq. (1), it is often quite useful to decompose the dynamical fields u α into characteristic fields. These characteristic fields are defined with respect to a spatial direction at each point, represented here by the unit normal covector field n k . Given a direction field n k we define the eigenvectors eα α of the characteristic matrix A k α β : where v (α) denotes the eigenvalue (also called the characteristic speed). The indexα labels the various eigenvectors and eigenvalues, and there is no summation over α in Eq. (10). Since we are interested in hyperbolic evolution systems, the space of eigenvectors will have the same dimension as the space of dynamical fields, and the matrix eα β will be invertible. Given these characteristic eigenvectors it is often useful to re-express the dynamical fields in terms of this eigenvector basis. Thus we define the characteristic fields uα (or the characteristic projection of the dynamical fields) as It is straightforward to show that the evolution of the characteristic fields is determined by where the projection operator orthogonal to n i is defined by P ij = g ij − n i n j , and spatial indices are raised and lowered with g ij and g ij respectively.
The characteristic fields for the fat Maxwell evolution system are a collection of fields of the form uα = The characteristic fields Z 1 , Z 2 i and Z 3 ij have characteristic speed 0; the fields U 1± i have speeds ±1, and the fields U 2± have speeds ± √ γ 1 γ 2 . All the characteristic speeds are real, and the characteristic fields are linearly independent (and depend continuously on the unit vector n k ) whenever γ 1 γ 2 > 0. Consequently the fat Maxwell system is strongly hyperbolic when γ 1 γ 2 > 0. We also find that the fat Maxwell evolution system is symmetric hyperbolic when the parameters γ 1 and γ 2 satisfy the more restrictive conditions, 0 < γ 1 < 4, and 1 3 < γ 2 . We note that the characteristic eigenvectors eα α for the fat Maxwell system depend only on the spatial metric g ij and the normal vector n i . It follows that the last term on the right side of Eq. (12) does not depend on any derivatives of u α at all. Thus the right side of Eq. (12) depends only on the transverse (to n i ) derivatives of u α : This important feature of the characteristic evolution equations is satisfied by many systems of interest to us, including the Einstein evolution system. It is also useful to have the inverse transformation u α = e αα uα, where e αα is the inverse of eα α . For the fat Maxwell system this inverse transformation has the form: The characteristic decomposition of the dynamical fields is essential for fixing boundary conditions. We will return to a more complete discussion of boundary conditions in Sec. IV.

III. ACTIVE CONSTRAINT CONTROL
Unless the constraint evolution system Eq. (3) is hyperbolic, it will not guarantee that the constraints remain satisfied (within the domain of dependence of an initial surface) even if they are satisfied initially. Thus the constraint evolution system must be hyperbolic in any self-consistent and physically reasonable system of constrained evolution equations. We assume that any system considered here has a hyperbolic constraint evolution system. We also assume that the constraint evolution system satisfies the somewhat stronger condition of symmetric hyperbolicity: In particular we assume that there exists a symmetric, positive-definite tensor S AB on the space of constraint fields which symmetrizes the characteristic matrices of the constraint system, for all k. When such a symmetrizer exists, we can define a natural norm on the constraints: The constraint energy E and its associated current E k are defined by This constraint energy can be used to define a norm E on the constraints, since E = 0 if and only if all the constraints are satisfied at each point. It is straightforward to determine the time evolution of E using the constraint evolution equations for any symmetric hyperbolic constraint evolution system: The quantities E k and E AB may depend on the dynamical fields u α and their spatial derivatives ∂ k u α (but not on higher spatial derivatives of u α ). In an evolution system Eq. (1) that includes parameters γ a multiplying constraint terms, such as the system defined by Eqs. (8) and (9), the associated constraint evolution system Eq. (3) and the constraint energy system Eq. (25) will also include terms that depend linearly on these parameters. Integrating Eq. (25) over the spatial slice for such a system, we get an expression for the time evolution of the constraint norm which has the general form, where Q and R a are integrals of quantities that depend on the dynamical fields and their first spatial derivatives. The basic idea of active constraint control then is to adjust the parameters γ a that appear in Eq. (26) to control the evolution of the constraint norm E . For example the growth of E might be prevented by making the right side of Eq. (26) non-positive at the beginning of each time step in the numerical evolution. This control mechanism is a special case of the constraint control method developed by Tiglio and his collaborators [27,28]. It differs from Tiglio's particular implementation [27,28] in that the quantities Q and R a in our expression do not depend on second derivatives of the dynamical fields. Since these higher derivatives are more difficult to evaluate accurately in a numerical simulation, we expect that our constraint control mechanism will be more stable and robust. The constraints associated with the vacuum fat Maxwell system introduced in Sec. II satisfy the following evolution system as a consequence of Eqs. (8) and (9), This system has the same general form as Eq. (3) with In order to define a constraint energy, we need this constraint evolution system to be symmetric hyperbolic. The most general symmetrizer for this system (that can be constructed from the spatial metric g ab ) is given by where and ǫ ijk is the totally antisymmetric tensor volume element. The parameters χ a must be positive χ a > 0, and γ 1 γ 2 must also be positive γ 1 γ 2 > 0 to ensure that S AB is positive definite. We note that these conditions put no additional limits on the allowed ranges of the parameters: every strongly hyperbolic representation of the principal evolution system has a symmetric hyperbolic constraint evolution system. We now evaluate the various quantities that determine the evolution of the constraint energy, Eq. (25), for the fat Maxwell system. We find: Thus the expression for the time derivative of the constraint energy becomes, The right side of Eq. (34) is a divergence, so the integral of this equation over a spatial surface results in an expression that involves only boundary integrals: where n k is the outward directed unit normal to the boundary. Active constraint control for this system consists then of adjusting the sign of the parameter γ 2 to force the constraint norm E to decrease with time whenever it gets unacceptably large. We note that the fat Maxwell system is rather degenerate, since the right side of Eq. (35) contains only a surface integral. Thus constraint violation in the fat Maxwell system arises only from the influx of constraint violations through the timelike boundaries of the computational domain. This property makes this system rather simpler than the Einstein evolution equations where constraint violation can be generated from bulk terms in the equations as well. The simplicity of the fat Maxwell system allows us to study how best to control the influx of constraint violations across boundaries in some detail, but it does not allow us to evaluate how effective these methods are for controlling violations that arise from bulk terms in the equations.

IV. CONSTRAINT PRESERVING BOUNDARY CONDITIONS
A standard boundary condition used for hyperbolic systems is the maximally dissipative condition, which we define here to be the condition that sets the incoming components of the dynamical fields to zero. (More generally the term maximally dissipative has been used to describe a larger class of boundary conditions that guarantee that a certain energy flux at the boundaries is strictly outgoing, e.g. see [18].). To impose such a condition, we first decompose the dynamical fields into characteristic parts, as was done in Eq. (11), and then set to zero at the boundary all those characteristic fields whose characteristic speeds are negative. Let Παβ denote the projection operator that annihilates all the non-incoming characteristic fields: that is, let For a maximally dissipative boundary condition, we set Παβuβ = 0 at the boundaries. We often use a variation on this boundary condition, in which we set to zero the time derivatives of the incoming components of the characteristic fields: For the case of the fat Maxwell system discussed in Sec. II, these "freezing" boundary conditions reduce to where the incoming characteristic fields U 1− i and U 2− are defined in Eqs. (16) and (17). As we shall see in Sec. V, this "freezing" boundary condition does a poor job of preventing the influx of constraint violations through the boundaries.
Calabrese, et al. [19] have proposed an alternate method for constructing boundary conditions that prevent the influx of constraint violations. Their method involves decomposing the constraint fields c A into characteristic parts: where eÂ A represents the eigenvectors of the characteristic matrix of the constraint evolution system, and v (Â) represents the corresponding eigenvalue (or characteristic speed). The idea is to impose what amounts to maximally dissipative boundary conditions on the constraint evolution equations: that is, we set where ΠÂB is the projection operator that annihilates the non-incoming characteristic constraint fields. This condition must now be converted into a boundary condition on the dynamical fields of the principal evolution system u α . This is done through the equation that defines the constraints in terms of u α and its derivatives, Eq. (2). In particular we solve Eq. (41) for the normal derivatives of the incoming characteristic fields, in terms of the outgoing characteristic fields and tangential derivatives of the incoming fields. When this is possible, Eq. (41) becomes a Neumann-like boundary condition on (some of) the incoming characteristic fields. This boundary condition has the following general form (42) We illustrate this procedure below more explicitly (and perhaps more clearly) for the specific case of the fat Maxwell system.
The characteristic fields for the fat Maxwell constraint system are the collection of fields of the form cÂ = The fields Z 4 i and Z 5 ij have characteristic speed 0, while the fields U 3± have speeds ± √ γ 1 γ 2 . The only incoming characteristic field is U 3− . So the constraint preserving boundary condition sets U 3− = 0 on the boundaries of the computational domain. Using the definition of U 3− above, we see that this boundary condition is equivalent to the condition on the boundaries. For a solution that satisfies the constraint preserving boundary condition, Eq. (46), the evolution of the constraint energy norm Eq. (35) becomes Thus the constraint preserving boundary condition ensures that the constraint norm does not grow. Quite generally the constraint preserving boundary conditions of this type will ensure that surface flux terms do not contribute to the growth of the constraint energy. In order to convert the constraint preserving boundary condition into an explicit condition on the dynamical fields, we must express the incoming constraint field U 3− in terms of the characteristic fields uα. Using Eqs. (13)- (17) and (45) we obtain Setting U 3− = 0 we obtain an expression for the normal derivative of U 2− : This has the form of a Neumann-like boundary condition on U 2− , and has the same form as the general expression Eq. (42). The version of our code used to perform the numerical tests described in Sec. V imposes boundary conditions on the time derivatives of the incoming characteristic fields. We therefore convert the Neumann-like boundary condition on U 2− in Eq. (49) into a condition on its time derivative using the characteristic field evolution equation, Eq. (12). We simply replace the normal derivative n k ∂ k U 2− that appears in Eq. (12) with the expression from Eq. (49). Simplifying the results gives the following equation for the time derivative of U 2− at the boundary, This condition together with the freezing boundary conditions ∂ t U 1− i on the remaining incoming characteristic fields constitute our version of constraint preserving boundary conditions on the fat Maxwell system.

V. NUMERICAL RESULTS
In this section we present numerical experiments that illustrate the effectiveness of the various constraint control strategies discussed in this paper. All of these tests use the fat Maxwell evolution system, with a variety of topologies for the computational domain and with a variety of initial data for the dynamical fields. In Sec. V A we illustrate what happens when the equations are solved without any constraint control. These tests show that significant constraint violations (and in some cases constraint violating instabilities) occur in dynamical solutions of the fat Maxwell system on computational domains with open boundaries. In Sec. V B we study the use of the active constraint control mechanism described in Sec. III. Our tests show that this method is not numerically convergent, and is not very effective in controlling the growth of constraints in this system. And finally in Sec. V C we describe the results of using the constraint preserving boundary conditions described in Sec. IV. This method is shown to be numerically convergent and quite effective in controlling the growth of constraints in the symmetric hyperbolic subset of the fat Maxwell system.
All numerical computations presented here are performed using a pseudospectral collocation method. Our numerical methods are essentially the same as those we have applied to the evolution problem in full general relativity [1,2,3,31] and in studies of scalar fields in Kerr spacetime [32]. Given a system of partial differential equations where u α is a vector of dynamical fields, the solution u α (x, t) is expressed as a time-dependent linear combination of N spatial basis functions φ k (x): Spatial derivatives are evaluated analytically using the known derivatives of the basis functions: The coefficientsũ α k (t) are chosen so that Eq. (51) is satisfied exactly at N c collocation points x i selected from the spatial domain. The values of the coefficients are obtained by the inverse transform where w i are weights specific to the choice of basis functions and collocation points. It is straightforward to transform between the spectral coefficientsũ α k (t) and the function values at the collocation points u α N (x i , t) using Eqs. (52) and (54). The partial differential equations, Eq. (51), are now rewritten using Eqs. (52)-(54) as a set of ordinary differential equations for the function values at the collocation points, where G α i depends on u α N (x j , t) for all j. This system of ordinary differential equations, Eq. (55), is integrated in time using a fourth-order Runge-Kutta method. Boundary conditions are incorporated into the right-hand side of Eq. (55) using the technique of Bjørhus [33]. The time step is typically chosen to be half the distance between the closest collocation points, which ensures that the Courant condition is satisfied.
In order to provide a quantitative measure of convergence and the amount of constraint violation of our numerical solutions, we have defined several norms on the constraints c A and the dynamical fields u α . We have already defined the constraint energy E in Eq. (24), which provides a norm on the constraint fields. In computing E for these numerical studies we fix χ 1 = χ 2 = χ 3 = 1. We also define norms on the dynamical fields themselves: We compute the volume integrals in these norms, e.g. in Eq. (24) or (56), exactly using the appropriate form of Gaussian quadrature, and the maximum in Eq. (57) is taken over the appropriate set of collocation points at a particular instant of time. These norms are most useful for comparing solutions evaluated with different numerical resolutions. Thus we define where δu α = {δE i , δD ij } is the difference between u α at a given resolution and u α at the best (highest) resolution we computed. Differences between quantities at different resolutions are computed by evaluating and then subtracting the spectral series for each resolution at the points on the finer grid. In order to provide meaningful scales for these normed quantities we typically plot dimensionless ratios of expressions such as ||δu|| 2 L 2 /||u best || 2 L 2 and ||δu|| 2 L ∞ /||u best || 2 L ∞ . In the case of the constraint energy we typically plot E /||∂u|| 2 , where is a norm on the gradients of the fields. We are interested in seeing how these ratios behave as the resolution of the numerical solution is increased: Order unity values of these ratios, ||δu|| 2 /||u best || 2 or E /||∂u|| 2 , indicate a complete lack of numerical convergence or solutions that are dominated by constraint violations, respectively. Values of these ratios of order 10 −34 correspond to double precision roundoff error.

A. No Constraint Control
In this section we illustrate the results of finding numerical solutions to the fat Maxwell evolution system Eqs. (8)-(9) using no constraint control at all. We examine three separate cases: First we look at evolutions on a computational domain with topology T 3 , a 3-torus. The differential equations governing the fat Maxwell system allow no constraint growth on domains without boundaries. So this first test is to verify that our code accurately reproduces "exact" constraint conservation in this case. Next we examine the evolution of a representation of the static Coulomb solution on a computational domain with topology S 2 × R, a spherical shell. Finally we study a highly dynamical solution on a computational domain with topology S 2 × R using freezing boundary conditions and no constraint control.
The evolution of the constraint energy norm E for the fat Maxwell system is driven entirely by a boundary term, Eq. (35). Thus we expect the constraints to be satisfied exactly for evolutions on a computational domain without boundary. To confirm that our numerical code correctly models this, we solve Eqs. (8)-(9) on a computational domain with topology T 3 , i.e. within a 3-torus. In particular we choose coordinates x, y, and z in the interval [0, 2π], and impose periodic boundary conditions. The basis functions used in our pseudospectral method are sines and cosines. We use initial data for this test in which the electric field is set to zero, E i = 0, and each component of the vector potential is set to be a cylindrical Gaussian pulse: where the width of the pulse is set to w = 0.5, and the center is placed in the middle of the computational domain, c y = c z = 3.14. The shape of this pulse is selected so that the value of the Gaussian falls below double precision roundoff, 10 −17 , at the periodicity "boundaries" of the domain, y = 0 and y = 2π etc. This ensures that these data are smooth on T 3 to roundoff accuracy. The initial data for D ij are set to the numerically determined values of ∂ i A j . We use the gauge choice φ = 0 throughout this evolution. Because these initial data are effectively two-dimensional, we can place as few as two collocation points in the x direction for computational efficiency. Figure 1 shows a convergence plot for this case that was run with evolution parameter values γ 1 = 1/γ 2 = −0.1, and resolutions N y = N z = 10, 20, 30, 40, and 50 collocation points. We see from Fig. 1 that the differences converge to zero as the resolution is increased. Figure 2 illustrates the amount of constraint violation in these runs. These curves, which increase approximately linearly with time, have magnitudes that are roughly proportional to the number of numerical operations performed multiplied by double precision roundoff error. Thus, the finer resolutions have larger errors than the coarser ones, since the finer resolutions require a larger number of timesteps and a larger number of numerical operations at each step. As expected from Eq. (35), we see that the constraints are satisfied essentially exactly when the domain has no boundaries. We have also computed evolutions for these initial data on T 3 using other values of the evolution parameters. In particular we have computed evolutions with γ 1 = 1/γ 2 = 0.1, and also evolutions that switch back and forth between these positive and negative values at each time step. In all of these cases, we find the evolutions to be convergent with roundoff level constraint violation.
Next we turn our attention to solving the evolution equations on a computational domain with topology S 2 × R, i.e. within a spherical shell. For our basis functions we choose Chebyshev polynomials for the radial coordinate and spherical harmonics for the angular coordinates. Although our basis functions are written in (r, θ, ϕ) coordinates, our fundamental variables are the Cartesian components of E i and D ij . To eliminate highfrequency numerical instabilities that sometimes develop during our simulations in S 2 × R, we apply a filter to the right-hand sides of Eqs. (55) before and after incorporating boundary conditions via the Bjørhus algorithm. The filter consists of simply setting high-frequency spherical harmonic coefficients to zero: If ℓ max is the largest index used in the basis functions Y ℓm at a particular resolution, then the largest ℓ retained in the right-hand side of the E i equations after filtering is 2ℓ max /3 − 1, and the largest ℓ retained in the right-hand sides of the D ij equations is 2ℓ max /3. This filter is a variation on the standard 2/3 rule used to remove the inevitable effects of aliasing whenever functions are multiplied using spectral methods [34]. For the first of these tests we choose initial data that corresponds to a static point charge that is located in the hole in the center of the computational domain. Thus we choose initial data with E i = −∂ i φ = r −2 ∂ i r and D ij = 0, appropriate for a unit point charge at rest at the origin. We then solve Eqs. (8)-(9) with γ 1 = 1/γ 2 = 0.1 on a computational domain with topology S 2 × R, de-fined by 1.9 ≤ r ≤ 11.9. (This is the same computational domain that we typically use to evolve single black hole spacetimes.) At both the inner and outer spherical boundaries we set the time derivatives of the incoming characteristic fields to zero, i.e. we impose freezing boundary conditions, Eq. (38). The scalar potential φ is held constant in time. We find that these numerical evolutions are stable and convergent and the constraints are preserved, as shown in Figs. 3 and 4. These computations were performed with radial resolutions N r = 11, 21, 31, 41, 51, 61, and 73 collocation points, and a fixed angular resolution with spherical harmonic index ℓ max = 5 (or equivalently, N θ = 6 and N ϕ = 12 angular collocation points). For ℓ max = 9 the results are indistinguishable on the scale of Figs. 3 and 4 except at the highest radial resolutions, indicating that the radial and temporal truncation errors dominate, as expected for a solution with little angular structure. This is a case (as we shall see) in which a time-independent solution is not always the best test problem to investigate the constraint-preserving properties of a system of evolution equations.
Finally, we examine a highly dynamical solution of the fat Maxwell system on the computational domain S 2 ×R, defined by 1.9 ≤ r ≤ 11.9. For this solution we choose initial data with E i = 0 and where r 0 = 6.5 and w = 1.0. The initial values of D ij are set to the numerically determined values of ∂ i A j . These initial data correspond to a pulse of radiation centered at r = r 0 . This pulse is neither spherically symmetric nor even axially symmetric, because the Cartesian components of the vector potential are set to spherically symmetric functions in Eq. (62); however, only a small number of spherical harmonics are sufficient to represent the full solution. The scalar potential is set to φ = 0 for these solutions, and we impose freezing boundary conditions, Eq. (38), on the incoming characteristic fields. Figure 5 shows a convergence plot for this case, using evolution parameters γ 1 = 1/γ 2 = 0.1 and ℓ max = 5, which confirms that the numerical solution is convergent. For ℓ max = 9 the results are indistinguishable on the scale of the figure. Figure 6 shows that significant constraint violations do exist in these solutions with seven different radial resolutions: the data from the three highest resolutions coincide on the scale of this diagram. Thus the constraints are violated, but the constraint energy is still convergent in these solutions. This indicates that the constraint violation is a property of the true solution of the continuum equations with freezing boundary conditions, rather than an effect caused by a defective numerical method. The constraint violation appears to be generated as the initially constraint-satisfying waves interact with the boundaries, starting at about t ≈ 4.
We have also computed evolutions for these dynamical initial data using negative values of the evolution parameters: γ 1 = 1/γ 2 = −0.1; the results are depicted in Figs. 7 and 8. Since the product γ 1 γ 2 is unchanged from the previous runs, the characteristic speeds of the system remain the same. And the definition of the constraint energy E (which depends on the ratio γ 1 /γ 2 ) also remains unchanged; so this allows us to meaningfully compare E for the two cases. For γ a < 0 the fundamental evolution system, Eqs. (8)-(9), is strongly but no longer symmetric hyperbolic as in the γ a > 0 case. Figures 7 and 8 show that these evolutions with negative values of γ a appear to be convergent, and have fractional constraint violations that are comparable with those in the positive γ a case. However as illustrated in Fig. 9, the evolutions in the negative γ a case have constraint violating instabilities. We note that we also find a numerical instability for the γ a < 0 case, and apparently for all cases in which the fundamental evolution system is strongly but not symmetric hyperbolic. This numerical instability appears to be associated with the angular discretization. It grows exponentially in time, and becomes worse at higher angular resolutions. However, for the angular resolutions we use, and for choices of γ a near the symmetric hyperbolic range, such as the case shown here, γ 1 = 1/γ 2 = −0.1, the numerical instability is negligible compared to the constraint-violating instability shown in Fig. 9. Only by going to higher angular resolution can one see any nonconvergent growth at all on the time scales we consider here. For ℓ max = 9 the instability is visible only at late times (t ≈ 200) for the highest radial resolutions in the quantities ||δu|| 2 L 2 and ||δu|| 2 L ∞ , and is not visible in plots of E . To construct a quantity that is sensitive to this instability, we repeated the runs shown in Figs. 5-9 at angular resolutions ℓ max = 5, 7, 9, and 11 and computed the norms of differences of the fundamental fields at adjacent angular resolutions. These norms are plotted in Fig. 10. We see no indication of this numerical instability for values of γ a in which the fundamental evolution system is symmetric hyperbolic (for example, the solid curves in Fig. 10 are convergent). For choices of γ a very far from the symmetric hyperbolic range (such as γ 1 = 1/γ 2 = 10), the instability grows much more rapidly and dominates the results. Although it is possible that the numerical instability can be cured by modifying our angular filtering algorithm, for the purpose of this paper we simply consider only values of γ a and angular resolutions for which the timescale of this instability is longer than that of other effects we wish to study.

B. Active Constraint Control
In this section we investigate the use of the active constraint control methods described in Sec. III. For the case of the fat Maxwell system this active control consists of switching the sign of the evolution parameters γ 1 and γ 2 to ensure that the constraint energy E does not increase. In the previous section we presented two sets of numerical evolutions without constraint control differing only by the signs of γ 1 and γ 2 . The characteristic speeds and the definition of the energy E were the same for these evolutions. Both evolutions were convergent on the time scale considered here (t ≤ 100), and the fractional constraint violation was convergent and similar in these cases on the same time scale. We now investigate the possibility of switching between these two cases during a single evolution as a means of reducing the constraint violations. The strategy is to monitor the quantity on the right side of Eq. (35), and to change the signs of γ 1 and γ 2 whenever necessary to keep the right side negative. This should ensure that the constraint energy norm E is always decreasing, so the constraints should remain satisfied. Note that this method should work only as long as our numerical solution satisfies the equation governing the evolution of the constraint energy, Eq. (35). Figure 11 illustrates for the case γ a > 0 that this equation does remain satisfied to truncation error level for the runs discussed in Sec. V A; the plot for γ a < 0 is similar. Consequently we expected good results from this active control method. Figure 12 shows the constraint violation for a case in which the γ a are allowed to change sign at every time step in order to control the constraints. The signs are changed only if the right side of Eq. (35) becomes positive, and the current value of E exceeds the value it had after the first timestep. The latter condition is intended to prevent sign changes that attempt to reduce the constraint violation to less than truncation error. Since this constraint control method depends on Eq. (35) being satisfied, control should only be possible to truncation error level at best. Figure 12 shows that the maximum value of the constraint is smaller than for the uncontrolled case (plotted as a dotted line in Fig. 12). However, the improvement is only an order of magnitude at best, even for the highest resolution run. Even more disturbing is the lack of convergence of the constraint norm. The fundamental fields do not converge very rapidly (if at all) either, as can be seen from Fig. 13. This lack of effective constraint control is confirmed in Fig. 14, which shows that Eq. (35) is not satisfied very well for this case. It appears that the active constraint control mechanism severely degrades the convergence of our numerical simulations in such a way that Eq. (35) no longer holds to the needed or expected accuracy. Consequently the active constraint control method is able to reduce the constraint violations only by a small amount over the uncontrolled case. One effect that can destroy convergence in these tests is the fact that the control mechanism is applied independently for each resolution (as pointed out by Tiglio [27]). Therefore, at a given value of t, the evolution equations used for one resolution can be different (because of the sign of γ 1 and γ 2 ) than the equations used for another res- olution at that time. This effect should have significant consequences only on quantities computed using different resolutions, such as the differences plotted in Fig. 13. But quantities computed using a single resolution, such as E , should not be affected. When these latter quantities are compared for different resolutions on the same plot, as in Figs. 12 and 14, the graph will not look like the "classic" convergence test in which all curves have the same shape. But the curves should (if everything else in the method is convergent) decrease at roughly the correct rate. Because we lose a great deal of accuracy in Figs. 12 and 14 compared to their uncontrolled counterparts, Figs. 6 and 11, we believe that the fact that the control mechanism is applied independently for each res-olution is not the primary cause of the problem. In order to eliminate this effect on convergence, we repeated our simulations, but this time switched the sign of γ a only once at the time t = 4 for each resolution, regardless of the sign of the right side of Eq. (35) or the magnitude of E . In this case, exactly the same evolution equations are being solved at each resolution. As shown in Fig. 15, the convergence rate is severely reduced even in this case when the signs of γ a are switched at t = 4. Furthermore as shown in Fig. 16, Eq. (35) is violated after the signs are switched.
We now believe that this nonconvergence is caused by a lack of smoothness of the fields that is introduced by the discontinuous change in the evolution equations: Suppose the signs of γ a are switched at a time t 0 , and suppose that at time just before this switch (t = t 0 − ǫ) some outgoing characteristic fields at the boundary are nonzero. When the signs of γ a are switched, some of the outgoing and zero-speed characteristic fields will be converted to incoming characteristic fields, and vice versa, as can be seen from Eqs. (13)- (17). (For example, switching the signs of the γ a in Eq. [17] while keeping E i and D ij fixed yields U 2± after = −U 2∓ before .) Therefore, at a time just after this switch (t = t 0 +ǫ) the solution will contain some nonvanishing incoming characteristic fields near the boundary. However, our freezing boundary condition requires the incoming characteristic fields to be constant in time at the boundary. Since the incoming characteristic fields propagate inward, at times after the switch (t > t 0 + ǫ) a kink will appear in the profile of these incoming fields. This type of boundary-condition-induced kink is illustrated in Fig. 17. The sketch on the left in Fig. 17 illustrates an incoming characteristic field at time just after the switch (t 0 + ǫ), and the sketch on the right shows the kink in this field that develops from the boundary condition. The existence of such a kink in the evolution fields greatly reduces the convergence rate of our spectral evolution method. And even for finite-difference methods such a kink is likely to reduce the convergence as well, since a kink in the fundamental quantities translates into a discontinuity in the constraints. Unless great care is taken to ensure that discontinuous solutions are treated properly (a standard problem in computational fluid dynamics but quite foreign to vacuum numerical relativity because the gravitational field is not expected to have physical shocks), Eq. (35) is likely to be violated and the constraint preserving mechanism will fail. We have made several attempts to replace the freezing boundary condition with a condition that smoothly adjusts the value of the incoming fields at the boundary. Unfortunately none of these attempts have been very successful.

C. Constraint-Preserving Boundary Conditions
Finally, we have performed a series of tests on the constraint preserving boundary conditions described in Sec. IV. Figures 18 and 19 show evolutions of our dy- namical initial data on S 2 × R (analogous to that used in Figs. 7 and 8) with the boundary condition on U 2− now set according to the constraint preserving condition Eq. (50). The γ a are negative for the plots in Figs. 18 and 19. The constraints are satisfied, and the simulation appears to be convergent (except for a late-time angular numerical instability, not visible on the plots, that appears identical to the numerical instability discussed at the end of Section V A). We have also performed these evolutions using positive values of γ a , and the simulations appear to be convergent, completely stable, and constraint-preserving in this case. Figure 20 compares the unnormalized constraint energy for these evolutions with those run with freezing boundary conditions. Thus adopting constraint preserving boundary conditions clearly does improve the constraint preserving properties of these evolutions much more than the active constraint control method.
But this is not the entire story. Figure 21 shows the norm of the fundamental dynamical fields ||u|| 2 for evolutions using constraint preserving boundary conditions with γ 1 = −0.1 (solid curves) and γ 1 = 0.1 (dotted curves). This plot shows that while the positive γ 1 evolutions are stable, those with negative γ 1 are not. A more extensive sampling of the parameter space reveals that evolutions preformed with γ 1 = 1/γ 2 = {0.1, 1.0, 2.5, 2.9} (for which the principal evolution system is symmetric hyperbolic) appear to all be convergent, constraint preserving, and stable. Conversely, we find that evolutions performed with γ 1 = 1/γ 2 = {−1.0, −0.1, 3.5, 4.1} (for which the principal evolution system is strongly but not symmetric hyperbolic) are all constraint preserving but unstable. These evolutions are numerically convergent for the resolutions and time scales we have tested (except for a slowly-growing angular numerical instability, not visible on Figure 21, that appears at late times or for high angular resolutions). Therefore, the type of growth seen in Figure 21 appears to represent a solution of the partial differential equations. Since these solutions do satisfy the constraints, the driving force for these instabilities must be an excess of incoming radiation that is reflected back into the computational domain by the boundary condition. We refer to this type of instability as a boundary condition driven instability. Thus the constraint preserving boundary conditions are a great improvement over the other methods studied here, but they do not completely eliminate all the instabilities in these strongly hyperbolic cases. Further study will be needed to determine whether these boundary conditions can be improved.

VI. DISCUSSION
This paper explores the effectiveness of two methods for controlling the growth of constraints in hyperbolic evolution systems. Using an expanded version of the Maxwell evolution system-which we call the fat Maxwell system-we showed that significant constraint violations and in some cases even constraint violating instabilities occur when the evolutions are performed using "standard" numerical methods and boundary conditions. We show that the active constraint control method (which has been studied by Tiglio and his collaborators [27,28]) is not very effective in controlling the growth of constraints in the fat Maxwell system when spectral numer- Norm of the fundamental variables for fat Maxwell on S 2 × R with γa < 0 (solid) and γa > 0 (dotted), using constraint-preserving boundary conditions. Even though the constraints are satisfied for γa < 0, the fundamental quantities increase exponentially, but in a convergent manner. Seven different resolutions are depicted for each case, but only the lowest resolution curves are distinct at the scales shown.
ical methods are used. This lack of effectiveness appears to be caused by the non-smooth nature of the control mechanism for this system. We also show that constraint preserving boundary conditions are very effective in suppressing the constraint violations in this system. Unfortunately these constraint preserving boundary conditions did not eliminate the instabilities for the strongly (but not symmetric) hyperbolic evolution systems. In these cases these boundary conditions merely converted a constraint violating instability into a boundary condi-tion driven instability. Generalizing these methods to more complicated systems like the Einstein evolution equations should be straightforward. For more general systems the analogue of the constraint energy evolution equation, Eq. (50), will contain both boundary terms like the fat Maxwell system and also volume terms, so that constraint violations can be generated both at boundaries and in the bulk. We expect that constraint preserving boundary conditions alone will not be sufficient to control the constraint violating instabilities that occur in these more general systems. Instead we expect that some combination of methods will be needed. The disappointing results obtained here with the active constraint control mechanism suggest that significant improvements will be needed to make this method useful for helping to control the growth of constraints in the Einstein system for evolutions based on spectral methods.