You are an expert at summarizing long articles. Proceed to summarize the following text: 
this is the first in a series of papers about the numerical treatment of collapsing gravitational waves using a new pseudospectral code developed for the problem . in the early  @xmath1s critical phenomena were discovered in gravitational collapse  @xcite , in spherical symmetry , with general relativity minimally coupled to a massless scalar field . one aspect of the finding was that , amazingly , the critical solution dividing the formation of a blackhole from dissipation of the field , was unique , in the sense that if one takes any one parameter family of initial data , with the parameter controlling somehow the strength of the data , and tunes this parameter to the threshold of blackhole formation , one finds that the same solution is always obtained , regardless of the family ! shortly thereafter similar phenomenology was reported in axisymmetric , vacuum general relativity  @xcite , or in other words in the collapse of gravitational waves . since then multiple studies have been performed to reproduce this finding , albeit with different initial data and numerical approaches , but without success . perhaps most strikingly , in  @xcite , numerical evidence of a different critical solution was presented . even if one completely accepts the available evidence for criticality in vacuum collapse , this obviously begs the question whether or not the naive expectation of uniqueness of the critical solution in axisymmetric , rather than spherical , collapse holds . roughly speaking there are two types of code being used used in  3d numerical relativity . the first uses the moving puncture method  @xcite , which consists , in essence , of a clever choice of evolved variables and gauge conditions , normally treated numerically by finite differencing . secondly is the pseudospectral method , most prevalently used with a first order generalized harmonic formulation of general relativity by the spec code  @xcite . recently , we presented a study of the collapse of gravitational waves using the moving puncture method  @xcite , in part to establish how close to the critical regime one can get with this standard approach . the conclusion being ; not very . major difficulties included the formation of coordinate singularities and a lack of accuracy . therefore one would like to tackle the problem using the pseudospectral approach to establish what can be achieved in that setting . we have thus developed a new pseudospectral code along the lines of spec , specializing the continuum and numerical method as much as possible towards the problem of vacuum gravitational collapse . the present paper represents the first outcome of this maneuver . herein we describe the formulation of gr employed , our boundary conditions , the code , calibration of the method , our initial data , our approach to axisymmetric apparent horizons , plus a suite of validation tests for gauge waves , gravitational waves , blackhole and collapse spacetimes . throughout we compare our results carefully with those in the literature . we aim to give a body of evidence for the correctness of the method that the reader will find compelling . with this out of the way , in subsequent papers we turn to the problem of critical collapse . a summary follows before the main text .    in section  [ section : ghg ] we look at a slightly modified version of the first order generalized harmonic formulation of  @xcite . we consider constraint preserving , radiation controlling boundary conditions , paying special attention to the constraint preserving boundaries . by considering the reflection of outgoing waves in the linear approximation we ultimately suggest modified conditions that should reduce spurious reflections caused by the use of constraint damping . we also suggest alternative gauge boundary conditions . next , in section  [ section : bamps ] , we outline the ` bamps ` code , including our carefully constructed cubed - sphere grids , which avoid clustering of grid - points in unfortunate positions of the domain . for the discretization we employ a pure chebyschev approach . we also discuss our ` octant ' symmetry implementation , the crucial patching - penalty approach for communicating data between neighboring coordinate patches , and finally the boundary implementation . in the follow - up section  [ section : energypenalty ] we complete the presentation of the penalty method by computing the penalty parameters appropriate for the semi - discrete system .    given the difficulties in the literature in reproducing the results of  @xcite it seems necessary to solve the problem in axisymmetry before moving to examine the collapse of fully  3d waves without symmetry . in our moving puncture gauge study @xcite a major disadvantage in using the bam code was that 3d grids were employed to evolve axisymmetric data . in section  [ section : axi ] we present our approach to evolving axisymmetric spacetimes with the ` bamps ` code , for which we employ the cartoon method  @xcite to reduce from the standard ` bamps ` 3d domains to a plane , by using the killing vector to evaluate any angular derivatives . we discuss various flavors of axisymmetric moment - of - time - symmetry initial data and their numerical construction . these initial data sets are evolved in a forthcoming study . we also give a detailed description of our formulation of the apparent horizon conditions in axisymmetry . to the best of our knowledge this is the first time that the regularity conditions on the symmetry axis have been carefully taken care of . this is important in later work as the search for apparent horizons will be our key diagnostic tool . the next three sections  ( [ section : experiments]-[section : supercritical ] ) contain a write - up of our development and validation tests . the tests include evolutions with the proposed gauge boundary conditions , which we find are helpful when using large gauge source parameters , as desired . they also include runs comparing the fully 3d , octant symmetry and cartoon evolutions , demonstrating that the various symmetry setups are well - behaved . in the evolution of single blackholes we test different gauges and boundary conditions , and following  @xcite , look at evolutions in which the blackhole is perturbed by a gravitational wave injected through the outer boundary . our results are in good agreement with the earlier studies . we then examine the evolution of supercritical brill waves , where , after the formation of an apparent horizon the run is continued after interpolation onto an excision grid , as used to evolve a single blackhole , which is needed to evolve data with a horizon for long - times . finally we conclude in section  [ section : conclusions ] . [ [ the - evolution - system ] ] the evolution system : + + + + + + + + + + + + + + + + + + + + +    we use the first order reduction of the generalized harmonic formulation with several free parameters . the full reduction from the second order einstein equations is presented in detail elsewhere  @xcite so here we give only a brief overview to establish our notation . throughout the paper in continuum equations we use the latin  @xmath2 for four dimensional indices , but  @xmath3 for spatial indices , with the exception of  @xmath4 and  @xmath5 , whose meaning when used as indices will be described shortly . greek indices are used to refer to the position in a state - vector , grid indices , or where otherwise needed . we start from the vacuum generalized harmonic formulation in second order form , @xmath6\,,\label{eqn : ghg_22}\end{aligned}\ ] ] for the unknown spacetime metric  @xmath7 with christoffels  @xmath8 . the constraints of the system are  @xmath9 , plus the standard hamiltonian and momentum constraints of gr . the gauge source functions  @xmath10 are freely specifiable , provided that they do not include derivatives of the metric , which would affect the principal part of the pde . the terms involving  @xmath11 are included so as to damp away high - frequency constraint violations  @xcite . the parameters  @xmath12 and  @xmath13 control whether or not the constraint addition made in the construction of the formulation is done either with the covariant or the partial derivative , or some combination . the latter choice has the effect of simplifying the constraint subsystem . in the code we use a first order reduction by introducing the variables  @xmath14 and  @xmath15 . the equations of motion are , @xmath16c_c\ , , \label{eqn : ghg11evsys}\end{aligned}\ ] ] with shorthands to be defined momentarily . the formulation here agrees with that of  @xcite except for the inclusion of the  @xmath12 and  @xmath13 parameters . we will either take the new parameters to vanish , or choose  @xmath17 . the lapse and shift are denoted  @xmath18 and  @xmath19 respectively . the unit normal to the spatial slices of constant coordinate time  @xmath20 is written  @xmath21 . when the normal is contracted with a tensor we sometimes use the abbreviation  @xmath22 , and likewise for the arbitrary unit spatial vector  @xmath23 . the induced metric on the slice is written  @xmath24 . in matrix notation this system can be written as @xmath25 with  @xmath26 , and principal matrix , @xmath27 and  @xmath28 containing all source terms . we use the shorthand for the christoffel symbols under the first order reduction , @xmath29 and will frequently use the abbreviation  @xmath30 . the system is symmetric hyperbolic , having the same principal part as a particular first order reduction of the wave equation . the characteristic variables are given by , @xmath31 with the projection operator  @xmath32 , and speeds , @xmath33 respectively . for future reference let us also note that a convenient way to transform to the characteristic variables is to write  @xmath34 , where here the indices represent the position in the state - vector  @xmath35 and where the similarity matrix is , @xmath36 which has left inverse  @xmath37 . but note however that  @xmath38 . the strength of this representation in practical terms is in avoiding special cases in the numerical implementation , like for example  @xmath39 , in the characteristic decomposition . [ [ gauge - source - functions ] ] gauge source functions : + + + + + + + + + + + + + + + + + + + + + + +    for the gauge source functions  @xmath10 we choose , @xmath40 our convention differs from that of both  @xcite and  @xcite in a trivial normalization of the spatial part with respect to the lapse function . writing the resulting gauge conditions in terms of the lapse and shift we get , @xmath41 with  @xmath42 the trace of the extrinsic curvature and  @xmath43 the contracted christoffel symbol of the spatial metric . before blackhole formation for the scalar functions @xmath44 we choose , @xmath45 with  @xmath46 some constants . by default we choose  @xmath47 and  @xmath48 , which naturally maintains the shift damping term even if the lapse function is close to zero , in contrast to the standard condition employed in spec  @xcite , which takes  @xmath49 . since we wish to study near - singular gravitational effects in the computational domain and avoid run - away growth of the shift vector this seems reasonable . we will report in later work on adjustments to these choices when evolving near - critical data . when evolving blackholes by excision we follow  @xcite taking instead  @xmath49 , although so far we have not found it necessary to use the  @xmath50 form of  @xmath51 . [ [ the - constraint - subsystem ] ] the constraint subsystem : + + + + + + + + + + + + + + + + + + + + + + + + +    the first order reduced harmonic constraints are , @xmath52 the terms without  @xmath10 are simply  @xmath53 . in these variables the vacuum adm hamiltonian and momentum constraints can be expressed as @xmath54k}+\tfrac{1}{2}d_{j}\phi_{kin } -\tfrac{1}{2}d_{i}\phi_{jkn } \nonumber\\ & \quad\quad-\tfrac{1}{2}\pi_{j[i}\phi_{j]nn}+\gamma^{lm}\phi_{mk[j}\phi_{i]ln}\nonumber\\ & \quad\quad+2\gamma_{an[i}\,\gamma^a{}_{k]j}\big)\,.\end{aligned}\ ] ] as stated above we use a subscript  @xmath4 to denote contraction with the normal vector  @xmath21 , but with the convention that  @xmath55 stands for the partial derivative , but _ with any such contraction outside of the derivative . _ we can put the hamiltonian and momentum constraints together as a four - vector of constraints , @xmath56 working with the first order system creates the reduction and closely related ordering constraints , @xmath57ab}=0\,.\end{aligned}\ ] ] the constraints  @xmath58 and  @xmath59 evolve according to , @xmath60\nonumber\\ & \quad-\gamma_1\gamma_2\beta^i\big(\tfrac{1}{2}g^{cd}c_{icd}n_a - c_{ina}\big)\,,\nonumber\\ \p_tc_{iab}&=\beta^j\big(\p_jc_{iab}+\gamma_1\p_ic_{jab}\big ) + \alpha\,\big[(1+\gamma_1)d_ig_{jn}\,c^j{}_{ab}\nonumber\\ & \quad-\gamma_2c_{iab}+\phi^j{}_{ab}\,c_{ijn}+\frac{1}{2}c_{inn}\pi_{ab}\big]\,,\end{aligned}\ ] ] where we have introduced the constraint , @xmath61 and where the notation  @xmath62 means take the partial derivative , and afterwards replace all first derivatives of the metric with the reduction variable  @xmath14 . up to lower derivatives in the contraints we find , @xmath63 where  @xmath64 denotes equality up to non - principal terms , the remainder having been suppressed for brevity . the equation of motion for  @xmath65 is readily derived by taking derivatives of that of  @xmath59 . notice that the parameter  @xmath66 serves to damp the reduction constraint . in the description of  @xcite the equivalent reduction variable is called  @xmath67 , with , including  @xmath12 and  @xmath13 in the natural way , @xmath68 in our variables . the difference is not substantial , being only that  @xmath69 appears slightly more naturally in the second order form of the equations . note that in  , the final term contains a piece which is simply the harmonic constraint in the pure harmonic case , but will act as a non - zero coefficient otherwise . [ [ first - order - reduction - of - the - constraint - subsystem ] ] first order reduction of the constraint subsystem : + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +    following  @xcite , a first order reduction of the constraint subsystem is formally introduced by defining the new variable  @xmath70 with , @xmath71 the principal part of this formal reduction is given by , @xmath72 the characteristic variables of the constraint subsystem are then found to be @xmath73 with speeds  @xmath74 and  @xmath75 respectively , where we use upper case latin indices to denote those projected by  @xmath76 . a suitable norm of the constraint violation is given by the constraint monitor which is defined as @xmath77    [ [ the - gravitational - wave - degrees - of - freedom ] ] the gravitational wave degrees of freedom : + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +    in vacuum the weyl scalars  @xmath78 can be expressed as , @xmath79\,,\nonumber\\ \psi_4&=m^am^b[\perp^{(p)bd}\!\!\!{}_{ab}k^ak^cr_{abcd}]\,,\label{eqn : weyl_scalars}\end{aligned}\ ] ] respectively . here we have introduced the null tetrad @xmath80 with  @xmath81 and  @xmath82 mutually orthogonal unit spatial vectors , and the projection operator , @xmath83 in terms of the first order ghg variables we can express the principal part of the riemann tensor as , @xmath84}{}^i -\gamma^j{}_{b}\p_i\phi_{ja[c}\gamma_{d]}{}^i + n_a\p_i\pi_{b[c}\gamma_{d]}{}^i\nonumber\\ & \quad - n_b\p_i\pi_{a[c}\gamma_{d]}{}^i + \gamma^i{}_{a}\p_i\pi_{b[c}n_{d ] } -\gamma^i{}_{b}\p_i\pi_{a[c}n_{d]}\nonumber\\ & \quad- n_a \gamma^{ij}\p_i\phi_{jb[c}n_{d ] }   + n_b \gamma^{ij}\p_i\phi_{ja[c}n_{d]}\nonumber\\ & \quad-\gamma_1\gamma_2 n_an^k\p_kg_{b[c}n_{d ] } + \gamma_1\gamma_2 n_bn^k\p_kg_{a[c}n_{d]}\nonumber\\ & \quad-\gamma_2\gamma^i{}_a\p_ig_{b[c}n_{d ] } + \gamma_2\gamma^i{}_b\p_ig_{a[c}n_{d]}\,.\end{aligned}\ ] ] of course this expression is unique only up to constraint additions . note that upon contraction with  @xmath85 and  @xmath86 to form the weyl scalar  @xmath87 , and after a single addition of  @xmath65 , we naturally form a projection of the incoming characteristic variable  @xmath88 . this is used in the construction of the boundary condition . the spatial vector  @xmath89 is taken to be the unit spatial normal to the boundary . [ [ boundary - conditions ] ] boundary conditions : + + + + + + + + + + + + + + + + + + + +    at the outer boundary we need to control incoming constraint violation , gauge perturbations and physical radiation . by default we initially impose , @xmath90 on the constraint subsystem assuming that the characteristic variable  @xmath91 is always incoming . these conditions are essentially those of  @xcite , with just the additional  @xmath92 term . other conditions for this variable will be motivated and tested in what follows . the remaining constraint subsystem characteristic variables may or may not be incoming , and are dealt with on this basis as described in section  [ subsection : bc_imp ] , but always according to the same prescription . for the gravitational wave degrees of freedom we choose , @xmath93 the lowest order member of a cascade of conditions on incoming radiation  @xcite , with given data  @xmath94 . examining   it is obvious that this is equivalent to setting , @xmath95 which is in practice how the conditions are implemented . for the remaining gauge degrees of freedom we choose either the improved gauge boundary conditions of  @xcite , @xmath96&\,\hat{=}\,0\,,\label{eqn : g_bcs}\end{aligned}\ ] ] or the alternative , @xmath97+\gamma_2\phi_{scd}\nonumber\\ & \qquad\qquad+r^{-1}(u_{cd}^{\hat{+}}-2n_{(c}h_{d)}+\gamma_2g_{cd})\ , \big]\,\hat{=}\,0\,,\label{eqn : g_bcs_freeze}\end{aligned}\ ] ] with given data  @xmath98 , which we will often take to vanish , and where the overbar derivative notation has the same meaning as in equation  . these conditions are similar to the ` freezing ' gauge boundary conditions employed in  @xcite , but taking into consideration the discussion of gauge reflections given in  @xcite , and constructed so that the conditions are naturally applied to metric components ( in adm form ) and their derivatives , but excluding contributions from the gauge sources . we will typically try to choose the given data to be fixed in time , and such that initially the time derivatives vanish for these quantities . here we have introduced the gauge projection operator , @xmath99 the above boundary conditions are implemented in ` bamps ` using the bjrhus method  @xcite as in spec . details of the method are explained in section  [ subsection : bc_imp ] . for completeness here the constraint projection operator  @xmath100 is , @xmath101 and also plays an important role in the implementation of the boundary conditions , as they are again naturally written in the form  @xmath102 . [ [ generalized - harmonic - constraint - subsystem ] ] generalized harmonic constraint subsystem : + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +    we already saw the constraint subsystem of the first order reduction of the ghg system . but to get a better idea of the effect of the different constraint preserving boundary conditions let us consider now the subsystem without the reduction . we have , @xmath103 the shorthand  @xmath104 and the variable  @xmath69 that follows will be related to quantities present in the first order reduction of the ghg formulation shortly . we can equivalently express this as , @xmath105+(n^b\nabla_bn^c)\big[\nabla_cc_a + 2\gamma_4\gamma^d{}_{ac } c_d\nonumber\\ & -(\gamma_4-\gamma_5)g_{ac}\gamma^d c_d - 2 \gamma_0 n_{(a}c_{c)}\big ] + \gamma^c{}_{ab}n^bc_c\nonumber\\ & + r_{ab}c^b\,,\label{eqn : ft2s_cons}\end{aligned}\ ] ] where the variable , @xmath106 is used to allow for the most convenient form of these expressions , and the final term of   is in fact of second polynomial order in the constraints because of the vacuum field equations  . different choices of the constraint addition parameters  @xmath107 result in different behavior in terms of growth of the constraint fields . it is also obvious that different choices of these parameters can simplify the constraint subsystem , the natural choice apparently being  @xmath17 . [ [ linearization ] ] linearization : + + + + + + + + + + + + + +    let us linearize and consider the behavior of a set of fields that satisfies these equations on a fixed constraint satisfying background . we start with equation   and use the tetrad consisting of the null vectors  @xmath108 defined in   to decompose the first index of  @xmath109 . from this we obtain , @xmath110 for the linearization , where we are free to use the notation  @xmath58 for the linearized violation because the constraints are satisfied in the background . [ [ boundary - conditions-1 ] ] boundary conditions : + + + + + + + + + + + + + + + + + + + +    taking the standard setup at the outer boundary so that  @xmath23 , used in the construction of the tetrad , denotes the outward pointing spatial unit vector normal to the boundary . restricting our attention to boundary conditions that contain at most one derivative of the constraints , _ geometrically _ the most natural choice seems to be  @xmath111 in the first order ghg language these conditions are , @xmath112 whereas , discarding the first order reduction , those of   are instead , @xmath113 with either conditions one might guess that the choice  @xmath114 reduces reflections from the boundary , especially when using a non - harmonic  @xmath115 gauge . incidentally this choice also makes the two conditions almost coincident . suppose all derivatives of  @xmath116 tangent to the boundary vanish , and that the background is flat . then we can analyze the solutions in a plane wave approximation .    [ [ mode - solutions - on - flat - space ] ] mode solutions on flat - space : + + + + + + + + + + + + + + + + + + + + + + + + + + + + +    when linearized around flat - space this system takes the form , @xmath117 the right - travelling mode solutions are , @xmath118 with eigenfrequencies , @xmath119 a very desirable property for our boundary conditions would be that they absorb outward going waves perfectly , that is , without reflection . with this motivation high - order derivative boundary conditions on the gravitational wave degrees of freedom have been studied  @xcite , and implemented in the spec code  @xcite in order to absorb higher spherical harmonics of the weyl scalar  @xmath120 . in the current context absorption means that outgoing mode solutions , those associated with an  @xmath121 , lie in the kernel of the boundary conditions . this is only the case if we switch off the damping  @xmath122 . since the low order spherical harmonics might be expected to dominate in the gauge and constraint subsystems , optimizing against this phenomena may be more important than using high - order conditions for the gauge and constraint subsystems whilst neglecting the damping terms . [ [ remainder - of - mode - solutions ] ] remainder of mode solutions : + + + + + + + + + + + + + + + + + + + + + + + + + + + +    substituting these mode solutions into the boundary conditions  , or the natural geometric conditions  , each after appropriate linearization , and expansion at large frequency  @xmath123 gives remainders of order  @xmath124 , indicating that neither is the optimal that can be obtained by adding source terms to the constraint boundary conditions . taking instead , @xmath125 the remainder is rather of order  @xmath126 . there is some freedom in expressing these conditions in the first order ghg language , but we choose , @xmath127 the conditions   can be similarly rewritten . a similar analysis can be performed using the pure gauge subsystem presented in  @xcite , but we currently find that existing gauge boundary conditions are sufficient for our needs , so we do not present these calculations here . tests with the various boundaries are presented in section  [ section : experiments ] . having discussed the continuum system in the previous section , we now discuss details of our numerical implementation of the ghg system . for this we present the ` bamps ` code , which uses a pseudospectral method on cubed - sphere grids . the basic idea of the code is based on spec  @xcite , but in many details , such as the actual grid implementation and the outer boundary treatment , ` bamps ` differences are present . [ [ grid - types ] ] grid types : + + + + + + + + + + +    the numerical domain on which we solve the evolution equations in ` bamps ` is either a cubed - ball or a cubed - shell grid . each type is built up of multiple deformed cubes . each patch is described by two fundamental overlapping charts . in  _ local coordinates _ @xmath128 and  @xmath129 it is a rectangular box @xmath130\times[-1,1]\times[-1,1]$ ] . in _ global cartesian coordinates _ @xmath131 and  @xmath132 the cubes are transformed and rotated in such a way that when added together they build the desired domain . we give a detailed description in the following . the cubed - ball - grid includes the origin and has a spherical outer boundary . it consists of  @xmath133 coordinate patches :    the central cube : : :    is centered around the origin and ranges    from  @xmath134 to  @xmath135 in    the global cartesian coordinate directions . the transition shell : : :    transfers the grid from the inner cube grid to a spherical shell with    radius @xmath136 . it contains six patches . the outer shell : : :    consists of six patches which extends the grid with additional    cubed - shells up to the outer grid boundary    at  @xmath137 . the cubed - shell - grid is an excision grid , meaning that it does not include the origin . it is a special case of the cubed - ball - grid , consisting only of the six outer shell coordinate patches . [ [ cubed - sphere - coordinate - transformation ] ] cubed - sphere coordinate transformation : + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +    the coordinate transformation used in ` bamps ` to construct the grids introduced above relies on the so called `` cubed sphere '' construction . it was introduced in  @xcite and first applied in the context of numerical relativity in  @xcite . since then this idea was implemented in multi - patch approaches  @xcite . in contrast to many of the earlier examples , the numerical method of ` bamps ` does not require overlapping grids , which simplifies the discussion . in  @xcite , the coordinates are constructed by considering great arcs parametrized by equidistant angles . such angle coordinates are used in  @xcite , while  @xcite use an intermediate set of coordinates also given in  @xcite that does not have the equidistant angle property . in ` bamps ` the latter type of coordinates are employed . the concrete coordinate transformation is the following . first , the local coordinates of each patch are transformed to temporary global coordinates @xmath138 this patch , which is orientated in positive @xmath139 direction , will later be referred to as the master patch . from here , cyclic permutation is used to rotate the patches to their location in the sphere . the denominator  @xmath140 depends on where the coordinate transformation happens . for the patches of the outer shell it is @xmath141 in the transition shell its definition includes a transition function @xmath142 @xmath143 this coordinate transformation is constructed to transition from the inner cube to the outer shells . note that this transformation is uniform along the 3d diagonals , where the distance between inner and outer shell boundary is smallest . this significantly improves the time - stepping restriction in the transition shell . [ [ subpatches ] ] subpatches : + + + + + + + + + + +    each coordinate patch can be further divided into subpatches . subpatches are helpful for increasing resolution , and form the backbone of the parallelism of ` bamps ` . each master patch can be split into  @xmath144 subpatches with coordinates @xmath145\,,\end{aligned}\ ] ] with  @xmath146 and  @xmath147 . in practice we ensure subdivisions are made in such a way that subgrids of two neighboring patches match , and that neighboring patches and subpatches share grid - point positions on their respective boundaries . this is necessary because our current penalty - communication method does not deal with interpolating penalties . concretely we split the inner cube into @xmath148 subpatches . the transition and outer shell are divided in  @xmath149 or  @xmath150 subpatches in the radial direction . for the angular direction we choose the number of subpatches to be  @xmath151 . in fig . [ fig : bamps_subpatch ] we show a  2d sketch of the ` bamps ` cubed - ball grid subdivided into subpatches .    [ [ discussion ] ] discussion : + + + + + + + + + + +    it is straightforward to specify a mapping between a rectangular master patch and a cubed sphere , although some book keeping for the different patches and different types of shell transitions is involved . it may be useful to examine different such mappings in terms of a numerical quality criterion , say the size of the jacobian , and to minimize the distortions associated with the coordinate transformation . [ [ spatial - discretization ] ] spatial discretization : + + + + + + + + + + + + + + + + + + + + + + +    ` bamps ` uses the method of lines with a standard ode integrator to integrate in time . the right - hand - sides are approximated using a pseudospectral method . we use a linear transformation to map the local coordinates of each subpatch @xmath152 into a unit cube  @xmath153 ^ 3 $ ] . we discretize the subpatch by choosing gauss - lobatto collocation points in each dimension , for example , @xmath154 with  @xmath155 , and similarly in the other directions . the number of of grid points @xmath156 depends on the patch location in the grid . the central cube is discretized with  @xmath157 points . the radial directions of the transition and outer shell are filled with  @xmath158 and  @xmath159 points respectively . the number of angular points we chose to be the same as in the central cube to assure that we have matching grids . in fig . [ fig : bamps_subpatch ] we show on the right the gauss - lobatto discretization of a subpatch . [ [ basis - expansion ] ] basis expansion : + + + + + + + + + + + + + + + +    on the collocation points we expand all evolution fields  @xmath160 in each dimension in a spectral basis using chebyshev polynomials  @xmath161 @xmath162 and analogously in the remaining two directions . we use the pseudospectral approach and store not the expansion coefficients @xmath163 but the function values @xmath164 at the collocation points @xmath165 . [ [ derivatives ] ] derivatives : + + + + + + + + + + + +    the spatial derivatives of the evolution fields are computed by a matrix multiplication . for example in the  @xmath166-direction we have @xmath167 with the gauss - lobatto derivative matrix , @xmath168 where  @xmath169 at boundary points and  @xmath170 elsewhere . in practice we do not compute diagonal terms of the derivative matrix by the analytic formulas stated above but use the identity @xmath171 this _ negative - sum - trick _ maps a constant function explicitly to zero and is known to give the derivative matrix better stability as regards rounding errors  @xcite . in preliminary experiments we found that this gives slightly more accurate derivatives , but have not studied the influence on the accuracy of the simulations presented later in the paper .    [ [ filtering ] ] filtering : + + + + + + + + + +    we find that a crucial ingredient for numerical stability is the use of a filter against high - frequency growth . for this we follow  @xcite exactly . after every full time - step we apply the filter in each dimension . the filter is easily implemented as a matrix multiplication . for example , in the @xmath166 direction , we filter the function values by , @xmath172 with the filter matrix @xmath173 where  @xmath174 and  @xmath175 and @xmath176 are the chebyshev synthesis and analysis matrices respectively . [ [ time - integration ] ] time integration : + + + + + + + + + + + + + + + + +    we integrate the fields forward in time using a 4th order runge - kutta scheme . unless otherwise stated we fix the time - step ,  @xmath177 , with @xmath178 being the minimal cartesian spatial grid spacing of the whole domain . empirically we find that this choice for the time step always leads to stable numerical evolutions , in the sense that increasing resolution results in smaller errors . we have not not carried out a stability analysis of the fully discrete system . [ [ bamps - octant - grid ] ] bamps octant grid : + + + + + + + + + + + + + + + + + + +    when evolving octant symmetric data in ` bamps ` , it is possible to only evolve one eighth of the cubed ball grid . this saves computational and memory costs . in the ` bamps ` octant mode we choose an odd number of subpatches  @xmath179 and a odd number of grid points  @xmath180 and reduce the numerical domain to  @xmath181 ,   @xmath182 and  @xmath183 . this means that all subpatches containing one of the cartesian axes are cut in half along them . in these patches we use the symmetry conditions to construct special matrices which compute the derivatives and filters .      to glue all subpatches together , we have to impose appropriate conditions at the connecting boundaries of the subpatches . for this we apply the penalty method as is described in  @xcite . the main idea of this method is to add penalty terms for each incoming characteristic variable at the boundary to the right hand side of the evolution equations . we use the characteristic variables of the evolution system to formulate boundary conditions . on the boundary surface we define the outward pointing spatial normal vector @xmath89 . the characteristic variables of the evolution system are given in equation ( [ eqn : ghg11_char_vars ] ) with speeds ( [ eqn : ghg11_char_speeds ] ) . in vector notation we write @xmath184 incoming characteristic variables to the subpatch boundary have positive speeds . on these we want to impose the condition that they are equal to the outgoing characteristic variables of the neighboring patch . table  [ tab : in_out_char ] summarizes all incoming and outgoing characteristic depending on the lapse function  @xmath18 and the shift in @xmath89 direction , @xmath75 . [ cols="^,^,^,^,^,^,^",options="header " , ]     where the numbers in parentheses denote the number of those grids that were cut in half ( at the axis ) once , twice , or three times respectively , for the 3d grids , and once or twice for the cartoon grids . note that our current non - octant cartoon implementation is not optimal because we evolve the whole  @xmath139-@xmath132 plane , wasting effectively a factor of two . currently we use the code most often in cartoon octant mode , so fixing this does not have a high priority . looking at the table the main observation is that the expected reduction factor of eight ( four ) in the total number of grid points is present between the 3d ( cartoon ) and octant grids , but that this number is not so closely reflected in the grid breakdown , where we get only a factor six ( three ) in the total number of grids . this is obviously because there are many grids with fewer points . since our parallelization does not take this fact into account , it is possible that one mpi process is given _ all _ non - cut grids , and so we can expect that the speedup rate is determined to a large extent by ratio in the number of grids . as we make the domain larger the relative number of cut grids decreases , so we might expect that asymptotically the full speedup factors of eight or four can are attained . [ [ basic - dynamics-1 ] ] basic dynamics : + + + + + + + + + + + + + + +    although irrelevant for the octant cartoon comparison , since these data have not been used before , we give a brief description of their evolution . initially the peak of the kretschmann scalar occurs at  @xmath185 with a value  @xmath186 . this profile then oscillates about three times at the origin , attaining a peak value of  @xmath187 before rapidly dispersing . looking at the lapse we see the familiar behavior that at the origin it oscillates slightly before presenting a longer decrease , although at the minimum is only  @xmath188 , having started from  @xmath189 everywhere . afterwards this pulse propagates out , roughly following the disturbance in the kretschmann . looking at the shift component  @xmath190 along the @xmath139-axis we find that early on there is a growth which peaks at  @xmath191 , with value  @xmath192 . the development of the shift looks more like a slowly oscillating standing wave than a localized propagating feature . [ [ d - octant - cartoon - and - octant - cartoon - comparison ] ] 3d , octant , cartoon and octant - cartoon comparison : + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +    taking first the 3d and octant evolutions , we see near perfect agreement throughout the evolution . there are small differences however , starting from the beginning of the simulation at the level of round - off ; differences of  @xmath193 in metric components , which slowly drift as the evolution goes on . this behavior is expected because the derivative approximation differ at this level . similar differences were found between the other setups . these differences are never larger than the constraint violation , in for example  @xmath194 , and we have looked at convergence ( see section  [ section : convergence ] for more discussion ) with each setup , although not for this data , and find no indication of a problem . for the speed comparison we ran the code with each setup on 24 cores ( with hyperthreading ) of our local cluster  _ core12 _ with intel xeon x5650 processors . the octant run was a little more than  @xmath195 times faster than the 3d run , as expected given the foregoing discussion . the octant cartoon run was about  @xmath196 times faster than the pure cartoon test , which is a little disappointing . going from  @xmath197 to  @xmath198 radial subdivisions in the outer shells , this value increases to  @xmath199 , demonstrating the expected dependence . comparing the full 3d and octant cartoon runs , there was a gratifying speed up of nearly a factor  @xmath200 . the ` bamps ` numerical method gives us two options for increasing resolution . the first is to add grid - points in every domain , the second is to subdivide grids further , keeping the number of points inside each subpatch fixed . given fixed finite computational resources it is not obvious what is the optimum strategy to achieve the smallest possible error , because although we would expect that adding points brings spectral convergence , it also comes with a  @xmath201 dependence in the allowed time - step , whereas on the other hand , as we will see , adding more subpatches allows the code to scale up to a large number of processors . probably the optimal strategy relies on a balance between each . to examine the effect of each strategy in the simplest possible way , we evolved gauge wave initial data on the minkowski spacetime , which was setup by choosing  @xmath202 $ ] , @xmath203 , with  @xmath204 as usual , and otherwise the flat spatial cartesian metric and vanishing extrinsic curvature . the results are plotted in the four panels of fig . [ fig : bamps_conv ] and confirm our expectations . +        to demonstrate the necessity of the filter   we evolved a centered  @xmath205 brill wave . the results are plotted in fig . [ fig : bamps_filter ] . in the left panel we see that without filtering the constraint violation starts to grow exponentially in time , whereas with filter the growth is completely absent and the norm of constraints remains steady at a very low value . in the right panel we plot the magnitude of the fourth highest spectral coefficient of  @xmath206 in the transition shell as a function of time . this coefficient is the first that is directly unaffected by the filter . we see that the growth in the constraints seems to be associated with an explosion in the higher spectral coefficients . interestingly we tried the same experiment with gauge wave initial data , and did not see the effect , at least in the same time - frame . we expect that the same behavior would manifest if we were to evolve long enough . the obvious conclusion we draw from this is that it is important to test these methods with several data types to get a reliable picture of their properties . [ [ strong - scaling ] ] strong - scaling : + + + + + + + + + + + + + + +    the current ` bamps ` parallelization strategy is to obtain perfect scaling using many subpatches , and splitting these subpatches across many processors . the key is that , contrast to buffer zones required in the decomposition of a finite differencing grid , only  @xmath207d surfaces of points need be passed by network communication , making the relative time spent there negligible . in a finite differencing approach the relative size of the buffer zones decreases with resolution , but in practice can still be significant in production runs . in fig . [ fig : bamps_scaling ] we present strong scaling plots performed on the _ supermuc _ cluster located in lrz garching , with intel xeon  e@xmath208-@xmath209c processors . we ran the code in  @xmath210d . we took a grid with  @xmath211 total subpatches , and increased the number of cores used until we were computing one patch per core . we find perfect scaling . on the other hand ` bamps ` is currently not parallelized whatsoever at the subpatch level , which means that the maximum number of points per subpatch is in principle determined by the amount of memory available to one core . at least when running the code in cartoon mode however we do not find , in practical terms , that this is problematic . instead the  @xmath201 restriction in the time step makes increasing the number of points infeasible long before we are close to filling the available memory . in @xmath210d this may no longer be the case . we leave such considerations for future work .     on the supermuc cluster . here a grid with @xmath212 sub patches was used . in total this grid consists of @xmath211 patches . [ fig : bamps_scaling ] ] the main thrust of our development has been towards treating collapsing axisymmetric gravitational waves accurately . for super - critical data the cubed - ball grid is unsuitable after the formation of an apparent horizon . therefore the strategy for long - term evolution is to take the data and interpolate onto a cubed - shell grid , with the excision surface suitably positioned , changing the lapse and shift to be sure that the excision surface is a true outflow boundary . a necessary requirement is to treat a single blackhole , which is what we discuss here . [ [ kerr - schild - coordinates ] ] kerr - schild coordinates : + + + + + + + + + + + + + + + + + + + + + + + +    we evolve the schwarzschild solution in kerr - schild coordinates as was done with an earlier version  @xcite of the present code . although the current numerical method is not particularly close to that used previously , some components of the older code were inherited . importantly evolving this data allows a simple comparison with the previous method and results . in spherical polar coordinates the metric and extrinsic curvature take the form , @xmath213 with  @xmath214 the flat metric on the two - sphere , and @xmath215\,,\end{aligned}\ ] ] respectively . inside the code the line - element is written in cartesian coordinates in the standard way . more discussion of kerr - schild coordinates can be found in  @xcite . [ [ harmonic - killing - coordinates ] ] harmonic killing coordinates : + + + + + + + + + + + + + + + + + + + + + + + + + + + + +    we additionally evolve starting from the harmonic killing slicing described in  @xcite , which serves as a convenient starting point when transitioning from one generalized harmonic gauge to another . for this initial data , in spherical polar coordinates , the metric and extrinsic curvature are @xmath216 and @xmath217 with the remaining components vanishing . for this data spatially harmonic coordinates are obtained by building cartesians according to , @xmath218 the resulting metric has a coordinate singularity at  @xmath219 , with  @xmath220 implicitly defined in the obvious way from the new coordinates . the coordinate singularity is not a principle problem as we could just put the excision surface outside this radius . but ` bamps ` relies on standard cartesian coordinates in several places . so in the code we could transform in the standard way but then choose the gauge source function , @xmath221 with  @xmath222 the jacobian between the standard  @xmath223-index cartesians and harmonic cartesian  @xmath224 index coordinates  , the compound object  @xmath225 is defined by , @xmath226 with  @xmath227 and where indices are manipulated in the obvious way with  @xmath7 to obtain  . instead we just choose the gauge source function to be fixed at its initial value , as will momentarily be discussed . in this section we use the code exclusively in cartoon mode , on a cubed sphere grid . we start with the excision surface at  @xmath228 , and the outer boundary at  @xmath229 . in our base setup we take  @xmath230 radial subpatches each with  @xmath231 points per direction . the runs were performed on a desktop machine with an eight - core intel i7 cpu , which was able to compute at about  @xmath232 , the base run requiring about  @xmath233  mb of ram . [ [ killing - gauge - sources ] ] killing gauge sources : + + + + + + + + + + + + + + + + + + + + + +    given initial data which admit a time - like killing vector , we can ensure that the evolution of the system is trivial , at the continuum level , neglecting the effect of outer boundary conditions , by choosing the killing lapse and shift , and taking the gauge source functions  @xmath10 so that  @xmath234 initially . in particular we must choose , @xmath235    [ [ kerr - schild - evolutions - with - spec - ghg ] ] kerr - schild evolutions with spec ghg : + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +    we began by evolving the kerr - schild initial data with the standard formulation parameters of  @xcite , namely  @xmath236 and  @xmath237 on our base grid as just described , using the gauge boundary conditions  . immediately we see that the innermost subpatch has the largest constraint violation , peaked at around  @xmath238 in the  @xmath194 component of the harmonic constraint . this is not surprising because the innermost subpatch contains the part of the solution with the largest derivatives . the evolution successfully continues until the final time  @xmath239 . but after the initial expansion to  @xmath238 , a slow expansion in  @xmath194 is visible , and this growth becomes more rapid as the simulation continues . by the end , the maximum value of  @xmath194 is around  @xmath240 , with peaks appearing at the inner and outer boundary of roughly the same size . we then increased resolution from the base grid to  @xmath241 and  @xmath242 . the  @xmath243 point grid runs at about  @xmath244 , and the initial peak in the  @xmath194 constraint violation is reduced by a factor of about two , with this ratio of improvement slowly declining until the end of the evolution . the  @xmath245 grid runs at  @xmath246 , with both the initial magnitude of the violation and the ` slow expansion ' of the  @xmath194 constraint quashed , the peak being a factor  @xmath247 smaller than in the base run at the end of the simulation . the highest resolution  @xmath242 point grid ran at  @xmath248 , with the final improvement in  @xmath194 against the base run being a factor of  @xmath249 . since the largest constraint violation occurs in the excision subpatch an obvious question is whether or not the excision and outer boundaries would interact badly if they were on the same grid . although the issue is of little practical concern for production runs , for development it deserves a little attention , and therefore we evolved our base grid from before , but cutting the outer two subpatches so that the outer boundary lies at  @xmath250 . this test is not completely fair because the outer boundary conditions are expected to perform better as they are applied further out . we find that the initial peak in the violation of the  @xmath194 constraint is about five times greater than in the base run at  @xmath251 . at the end of the evolution again at  @xmath239 by coincidence the constraint violation in the restricted domain is smaller , but this is just because the slow oscillations in each simulation are out of phase .    [ [ kerr - schild - incoming - wave - evolutions - with - spec - ghg ] ] kerr - schild incoming wave evolutions with spec ghg : + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +    next we evolved the same initial data and gauge , but this time with the same domain as in fig . @xmath210 of  @xcite . to do this we took  @xmath252 radial subpatches , with the same base resolution as before , so that the outer boundary is placed at  @xmath253 . we similarly specify exactly the same given data for an incoming gravitational wave as in that study , taking in particular , @xmath254 with the vectors here the coordinate vectors defined in the obvious way . we take , @xmath255\,,\end{aligned}\ ] ] with  @xmath256 , @xmath257 and  @xmath258 . in fig . [ fig : incoming_wave ] we show the results from these experiments , obtained with a sequence of different resolutions . we plot the weyl scalar  @xmath120  , averaged over the outer boundary , @xmath259 the surface area of the outer boundary is  @xmath260 . fitting the highest resolution data between  @xmath261 and  @xmath262 we find a ring - down frequency of  @xmath263\sim0.372 $ ] as expected  @xcite . in this evolution we found that the apparent horizon oscillates slightly as the gravitational wave is absorbed , increasing the horizon mass   by about  @xmath264 , with  @xmath265 the adm mass of the analytic initial data . note that the gauge boundary condition being employed here is not identical to that used in  @xcite , so the agreement is remarkable . the effect of much larger pulses of gravitational radiation falling onto a blackhole using similar gauge conditions was studied in  @xcite .     in the outer boundary in the evolution of schwarzschild perturbed by a small gravitational wave injected through the boundary . in the left panel we see convergence of the constraints as resolution is increased . at lower resolutions a drift is present in the ring - down . there is good agreement with fig . @xmath210 of  @xcite , and the ring - down frequency agrees well with the analytical computation  @xcite . at the end of the test there is some disagreement with  @xcite , but since square - roots of very small quantities are being taken we expect this is caused by round - off error . it seems that on the cubed - sphere grid more resolution is needed to obtain clean results than with the spherical harmonic discretization used in  @xcite . this is perhaps not surprising , since the latter discretization is well - suited to the given data . [ fig : incoming_wave ] ]    [ [ discussion - of - and - comparison - with ] ] discussion of and comparison with  @xcite : + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +    the prior ` bamps ` study focussed on obtaining numerical stability in the evolution of a single schwarzschild blackhole with the kerr - schild slicing . the numerical method used a chebyschev - fourier - fourier spatial discretization on a single shell with a spin weighted spherical harmonic filter to prevent high frequency growth of the error . in that study the outer boundary condition employed simply fixed the incoming characteristic variables   to some given data , namely their initial values . this approach is possible only when the analytic solution is known , otherwise incoming constraint violations are generated . placing the inner boundary at  @xmath228 and the outer boundary at  @xmath266 , very long evolutions , say until at least  @xmath267 , could be performed with little resolution , in accordance with  @xcite . on the other hand , using this method , the naive boundary conditions rapidly deteriorated as the outer boundary was pushed out , and , crucially resolution did not help but rather made the problem worse . a possible explanation for the latter effect is that no filter was being applied in the radial ( chebyschev discretized ) direction , which have already seen is a crucial ingredient for stability with the current method . the likely cause of the boundary problem is that , as explained in  @xcite , boundary conditions that just freeze the incoming ghg characteristic variables are orders of magnitude more reflecting than the sommerfeld like choice contained in  . evidence for this is obtained in the current code by changing from the gauge boundary condition   to use instead , @xmath268&\,\hat{=}\,0\ , , \label{eqn : g_bcs_ruined}\end{aligned}\ ] ] evolving once more the kerr - schild initial data on the base grid . placing the outer boundary further out then results in greater reflections . however rather than trying to improve a condition only suitable for evolving known data , we immediately moved to the constraint preserving , radiation controlling conditions , with which this issue is completely absent . the first attempted implementation of a regular center in the ` bamps ` code was to use the chebyschev - fourier - fourier discretization with a double covering in the radial direction , similar to that employed in  @xcite . the approach was not successful , as we always eventually found irregularities in the numerical solution at the origin . an exponential filter was applied to the chebyschev coefficients in the radial direction , but to little effect . eventually we settled on the cubed sphere approach , in part because of the expectation that they will later be more convenient for mesh - refinement . other possible solutions to the problems we faced would be to use one - sided jacobi polynomials as in spec  @xcite or to employ a filter that projects the solution in another basis onto these polynomials . [ [ kerr - schild - evolutions - with - simplified - constraint - subsystem ] ] kerr - schild evolutions with simplified constraint subsystem : + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +    using our standard choice for the formulation parameters  @xmath17 , and taking  @xmath269 , returning to our base resolution from the tests with the spec version of ghg , we find that by  @xmath262 the  @xmath194 constraint is about  @xmath208 times larger than that we obtained before , and by the end of the simulation the new run has accrued a  @xmath194 constraint violation with a sharp peak at the outer boundary of order  @xmath270 . this result seems to be in contradiction to those of section  [ section : c_ex ] , until we remember that there the gauge boundary condition   was employed . increasing the constraint damping to  @xmath237 , the initial violation is comparable to the spec ghg evolution previously described throughout the evolution , and the spike at the outer boundary is suppressed by roughly an order of magnitude . at the end of this run the maximum of the  @xmath194 constraint occurs at the excision boundary with a value around  @xmath240 . this experiment thus highlights that the choice of the damping parameters and boundary conditions can be rather subtle .    . in the upper panel we plot the logarithm of the constraint violation  @xmath194 . in the latter case the violation is greatly reduced . in the lower two panels we show the lapse and shift ; the drift present when using   is practically absent with  .[fig : bh_bndy_cmp],title="fig : " ] . in the upper panel we plot the logarithm of the constraint violation  @xmath194 . in the latter case the violation is greatly reduced . in the lower two panels we show the lapse and shift ; the drift present when using   is practically absent with  .[fig : bh_bndy_cmp],title="fig : " ]    [ [ kerr - schild - evolutions - with - alternative - boundary - conditions ] ] kerr - schild evolutions with alternative boundary conditions : + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +    next we returned to the base grid , and switched to the alternative gauge boundary conditions  , with  @xmath17 and  @xmath237 . we find that the aforementioned growth in the constraints is completely eradicated , and the drift in the lapse and shift is also suppressed . evolving the same data with the same formulation and gauge boundary condition , but using the modified constraint preserving boundary condition   gives almost identical results . using instead the reflection reducing conditions   we see a small improvement in the violation throughout the simulation . repeating the experiment with the incoming gravitational wave injected through the boundary with the standard constraint preserving condition   and the gauge boundary conditions  , the growth visible in fig .  [ fig : incoming_wave ] is also completely absent , even on the base resolution  @xmath231 grid . [ [ harmonic - killing - slice - evolutions ] ] harmonic killing slice evolutions : + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +    we now returned to our base grid and resolution , taking the formulation parameters  @xmath17 , and  @xmath237 , evolving the harmonic killing slice with the gauge boundary condition  . the test successfully runs to  @xmath239 . comparing with the equivalent evolution of kerr - schild data , we see that initially near the excision boundary the  @xmath194 constraint violation is significantly greater in the harmonic killing test . by  @xmath251 this difference has accrued to around two orders of magnitude . later however , as the violation in the kerr - schild killing evolution starts to grow , it overtakes that of the harmonic killing evolution . at @xmath239 the peak of the constraint violation in the harmonic killing run is about an order of magnitude smaller than in the earlier test . as remarked before , in the kerr - schild test the inner and outer boundaries have roughly the same magnitude in the  @xmath194 constraint violation . interestingly the twin peaks are not present in the harmonic killing data because the outer boundary is hugely improved . this finding is consistent with the gauge wave tests presented in section  [ section : exp_g_bcs ] , although this test is somewhat easier for the gauge boundary conditions because of the complete lack of dynamics present in the gauge wave test . in the harmonic killing evolution we are evolving with pure harmonic slicing , and some non - zero spatial gauge source functions , which suggests perhaps that the growth at the outer boundary is predominantly caused by the use of a non - trivial gauge source function for the lapse function , as it interacts with the boundary . indeed looking once more at the lapse function towards the end of the kerr - schild evolution we see that it is drifting from its initial value , but that this effect converges away with resolution . in any case , the peak in the constraint violation at the outer boundary in the killing kerr - schild data is suppressed as the outer boundary is placed further out . [ [ harmonic - killing - slice - with - gauge - perturbation ] ] harmonic killing slice with gauge perturbation : + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +    a desirable property for a set of dynamical coordinates is that in the presence of a , perhaps approximate , timelike killing vector they quickly asymptote to a time - independent state . for an arbitrary physical or gauge perturbation there is no hope that this will occur , and nor can any finite set of numerical experiments prove that that there is a basin of attraction to a stationary state . we can however look for some indication of this behavior . to do so we start by taking the initial data for the killing harmonic coordinates , and then perturb the initial lapse function by gaussian as in the previous gauge wave evolutions . in terms of the first order ghg variables this is a slightly fiddly procedure , as compared with the use of lapse , shift and spatial metric , so we give a quick summary :    * set spatial metric and extrinsic curvature from the exact solution . * take the killing lapse and shift . use the conditions  @xmath271 and  @xmath272 to set the gauge source functions  @xmath10 . * add the desired perturbation to the lapse ( or shift ) and then transform to the first order ghg variables . we perturbed the lapse by a gaussian , @xmath273\,,\end{aligned}\ ] ] with  @xmath274 and  @xmath275 . a similar experiment was made in  @xcite , but starting from a maximal slice of schwarzschild to test the gauge driver system . we find that the perturbation in the lapse propagates away , rapidly leaving behind the solution with the harmonic killing data with unperturbed spatial coordinates , or at least negligibly perturbed . the greatest danger to the evolution is probably that the excision boundary fails to be outflow , but at least with this perturbation that does not occur . [ [ harmonic - evolutions - with - incoming - gravitational - wave ] ] harmonic evolutions with incoming gravitational wave : + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +    giving the same gravitational wave data   as previously , evolving with the standard boundary conditions   and   but using the harmonic killing gauge source functions . it is not obvious how , if at all the spacetime computed is related to that considered before , but in any case we find a very similar decay in  @xmath120 . remarkably the growth present in fig . [ fig : incoming_wave ] is absent even in this low resolution  @xmath231 test . [ [ the - transition - function ] ] the transition function : + + + + + + + + + + + + + + + + + + + + + + + +    as elsewhere , we follow  @xcite to transform from one generalized harmonic gauge  @xmath276 to another  @xmath277 . the composite source function is simply , @xmath278\,h^2_a\,.\end{aligned}\ ] ] the transition function is , @xmath279 in the following experiments we choose  @xmath280 and  @xmath281 . note that care must be taken to construct the time and space derivatives of  @xmath10 with the transition function . this choice results in gauge source functions that are only  @xmath282 at  @xmath283 , which could be avoided with a different transition function . it is not clear if this finite differentiability will have a large effect on extracted physical quantities from a simulation . [ [ kerr - schild - initial - slice ] ] kerr - schild initial slice : + + + + + + + + + + + + + + + + + + + + + + + + + +    for our first phase - in test , we started with the kerr - schild slicing of schwarzschild and evolved with  @xmath17 and  @xmath237 , on our base resolution grid . we took the gauge boundary condition   and the constraint preserving condition   ( including a  @xmath92 term ) . we used the wave gauge parameters  @xmath284 and  @xmath285 . the value of  @xmath286 here is much smaller than in our wave collapse evolutions . the reason for this is that when evolving a blackhole it is crucial that the excision boundary is pure outflow in the pdes sense . in other words the characteristic speeds must all have the same outward pointing sign . since the speeds in the  @xmath89 direction are like  @xmath287 this means that the shift can not become too small or else the excision boundary will fail , which in turn means that  @xmath286 can not be chosen too large . we therefore place the excision boundary deeper into the blackhole so that  @xmath288 and carefully monitor the coordinate lightspeeds at the inner boundary . note that this requirement is likely to cause difficulties when computing extreme gravitational waves , because on the one hand large shifts can result in poor resolution of important features , but on the other they may be required in some other region so that we may successfully excise the blackhole region . in the evolution we immediately see significant dynamics and that for example the peak of the  @xmath194 constraint violation along the  @xmath139-axis is two orders of magnitude greater than in our initial kerr - schild base run with killing gauge sources . the reason for this is presumably the presence non - trivial dynamics , plus the fact that we are excising nearer the physical singularity similar to the effect we saw with the harmonic killing slice . regardless , by  @xmath289 the data seem very close to stationary . the simulation then evolves to the target time @xmath239 , and remarkably at the end of simulation the constraint violation in  @xmath194 along the @xmath139-axis has a maximum value which is an order of magnitude smaller than in the base run . at no point does the excision boundary fail to be outflow . as a check of the axisymmetric apparent horizon finder we compare the results obtained with the simpler algebraic condition , @xmath290 which characterizes the position of the apparent horizon in spherical symmetry . we find near perfect agreement throughout . the apparent horizon moves from its initial radius @xmath291 inwards until it reaches @xmath292 around @xmath293 . from there the horizon starts to grow again and seems to settle down to @xmath294 . however in our lowest resolution run , a small drift of the horizon outwards is visible . at late time of the simulation , around @xmath295 , this drift accelerates and we observe that the horizon becomes aspherical . higher resolution runs show that this effect converges away .    , and the radius at which the outward lightspeed vanishes  @xmath296 , computed on our base grid with inner boundary at  @xmath297 . to successfully excise , the speed must be negative at the inner boundary . observe that excision _ exactly _ on the apparent horizon is not possible throughout all of the run.[fig : harm_phase - in ] ]    [ [ harmonic - initial - slice ] ] harmonic initial slice : + + + + + + + + + + + + + + + + + + + + + + +    since the stationary fully harmonic coordinates are singular at  @xmath219 , one might guess that the stationary spatial generalized harmonic coordinates with gauge source functions   are also singular at some radius on the killing slice , at least for some range of the parameters  @xmath44 . given the broad experience in using these coordinates in binary blackhole simulations , the naive expectation would be that , if present , this coordinate singularity is pushed further towards the physical singularity rather than out towards the event horizon for standard choices of the gauge source functions . but this behavior is not clear . to truly resolve the issue one could simply solve for such coordinates along the lines of  @xcite , but this we defer for the future . instead we performed simulations varying the initial excision surface from the base grid excision radius  @xmath298 down to  @xmath299 in steps of  @xmath300 . unsurprisingly we find that initially the constraint violation , is greater in the excision subpatch as the inner boundary is placed closer to the singularity , amounting to about an order of magnitude in the  @xmath194 constraint between the  @xmath288 and  @xmath301 boundary runs by  @xmath302 . besides this there is little to distinguish between the five runs , and at least down to this excision radius no sign of a coordinate singularity forming . by eye the lapse function in the shared part of the domain agrees very well throughout the evolution . although a slight drift between them is present towards the end of the test , this is acceptable since the outer boundary conditions are being imposed at different radii , the solutions need not agree everywhere . there is however a time around  @xmath303 above which the runs with inner boundary  @xmath304 fail to be outflow at the excision surface . assuming that this is not caused by numerical error this means that boundary conditions are required at the surface . it furthermore means that convergence of the numerical scheme as resolution is increased is impossible . the fact that this does not correspond to a catastrophic failure of the code is inconvenient , because it indicates that great care must be taken in monitoring the excision surface . on the other hand , since placing the excision boundary very far in has a large cost in accuracy , a careful balance must be struck . in the spec code this is taken care dynamically of by a control mechanism  @xcite which ` bamps ` does not yet have . in fig . [ fig : harm_phase - in ] the relationship between the character of the excision boundary and the apparent horizon is examined . comparing the initially harmonic and kerr - schild slice evolutions with excision radius  @xmath288 we find that although the lapse functions initially disagree , by about  @xmath305 they have exactly the same profile and lie almost on top of one another . after this time the agreement is maintained . in this section we present the numerical evolution of a centered brill wave , see section  [ subsection : brill ] , with  @xmath306 . this highly supercritical initial data is used as a test case for our excision algorithm for a dynamically forming blackhole . our dynamical excision method currently consists of the following steps :    * 1 . evolve to collapse : * evolve on cubed ball grid , running the apparent horizon finder in ` daemon ' mode . the finder then triggers a ` bamps ` checkpoint once a horizon is found .    * 2 . go - to excision grid : * next interpolate the checkpoint data onto a cubed - sphere grid . in this step we want to place excision boundary just inside the apparent horizon , but as we have already seen in the single blackhole evolutions this may not always be possible , as some wiggle room is needed to allow for dynamical behavior of the horizon . this can require some experimentation , although fine - tuning does not seem necessary . regauge : * adjust the lapse and shift to ensure that the excision boundary is pure outflow . as a particular choice , we take lapse and shift from kerr - schild slicing of schwarzschild , @xmath307 and translate to cartesian components in the obvious way . it is desirable that the radial coordinate light - speeds are close to zero , preferably positive , at the apparent horizon , since this determines the direction of motion of the horizon . therefore we choose the  @xmath308 parameter to satisfy this condition reasonably well , although again without particular fine tuning . safety - net evolution : * we then use single blackhole gauge source parameters like  @xmath309 and  @xmath310 . during the evolution we use a safety net . if any coordinate light - speed on the excision boundary reaches a given threshold , typically  @xmath311 we again regauge to guarantee the outflow character is maintained . we monitor the apparent horizon , and if it falls off of the numerical domain we return to an earlier checkpoint , regauging with a smaller  @xmath308 to avoid this behaviour . as the horizon expands we monitor the position and periodically return to the * go - to * step above , excising further out and regauging with a greater  @xmath308 . [ [ discussion-1 ] ] discussion : + + + + + + + + + + +    as currently implemented this procedure requires that some steps be performed by hand . the numerical results in the following subsection serve to demonstrate ` proof of principle ' of this algorithm . on the other hand it seems at least clear how those steps should be automated . at the regauge step the use of the first order ghg variables is again a little fiddly . much more convenient would be if the lapse and shift were readily available as variables . but the procedure is similar to that described in the gauge perturbation tests in section  [ subsection : frozen_sources ] , so we do not give full details . also at the regauge step , it might be good to choose lapse and shift by abandoning the spherical ansatz and imposing that the coordinate light - speeds at the apparent horizon vanish . the spec approach to controlling the excision surface is much more sophisticated , employing a control mechanism  @xcite , we hope to avoid that investment in the near future . because we are interested in the collapse of waves to form , presumably , a single blackhole it seems reasonable to use a simple approach if at all possible . one aspect of the method that is not very aesthetically appealing , is that by changing the lapse and shift in discrete steps we are computing a spacetime , or patch of spacetime in coordinates that are not globally smooth . another issue associated with this is that of geometric uniqueness , which for the ibvp is an open question . nevertheless one expects that the differences to the computed spacetime with one choice of regauging parameters or another will be rather small in practice , so this does not represent an immediate practical concern .      [ [ initial - data - and - grid - setup ] ] initial data and grid setup : + + + + + + + + + + + + + + + + + + + + + + + + + + + +    we evolved a centered brill wave as described in section  [ subsection : brill ] , with seed function  . we chose a centered  @xmath312 wave with  @xmath306 . the adm mass of this initial data is  @xmath313 . the maximum of the kretschmann scalar in the initial data occurs at the origin , taking the value  @xmath314 . following the algorithm just outlined , we began on a cubed - ball grid with  @xmath315 , and  @xmath316 points per cube , with internal boundaries  @xmath317 ,  @xmath318 and the outer boundary placed at  @xmath319 . we ran the code in cartoon mode on our local cluster  _ quadler _ with  @xmath320 cores . we evolved with the generalized harmonic gauge , as in section  [ section : experiments ] in the evolution of a _ much _ weaker  @xmath321 brill wave , now with the gauge parameters  @xmath322 and  @xmath323 . at coordinate time @xmath324 we first found an apparent horizon with mass  @xmath325 . [ [ continuation - to - code - crash ] ] continuation to code crash : + + + + + + + + + + + + + + + + + + + + + + + + + + +    if we continue this evolution without going to an excision grid after the apparent horizon forms , we find that the constraints inside the apparent horizon rapidly grow along with the kretschmann scalar . the run then crashes at roughly  @xmath326 . this gives the clear signal that if we are to examine the final masses of blackholes formed during collapse , using the ghg formulation , a robust excision algorithm will be essential . in fact at  @xmath327 the horizon has a mass of  @xmath328 on the cubed - ball grid , but at the end of our excision simulation , to be described momentarily , we find that  @xmath329 after apparent horizon formation it has mass  @xmath330 . in the first critical gravitational wave collapse paper  @xcite , the blackhole masses were evaluated roughly  @xmath331 after apparent horizon formation , according to a prescription based on the quasinormal modes of schwarzschild . comparing those values with ours is difficult because we use different time coordinates , but the basic expectation is that the maximal slicing condition is more `` singularity avoiding '' than one of our generalized harmonic gauges , and therefore we might expect to obtain comparable results if we can evolve for a similar coordinate time after the appearance of a horizon . this is , however , not clear and deserves further investigation . in any case without excising the blackhole region , the meager  @xmath332 after collapse is clearly insufficient . we have seen in  @xcite that with the moving - puncture method this type of data also did not result in successful evolutions beyond apparent horizon formation . but here at least a concrete improvement has been made , in that we find an apparent horizon _ before _ the method fails !    [ [ evolution - on - excision - grid ] ] evolution on excision grid : + + + + + + + + + + + + + + + + + + + + + + + + + + +    checkpointing the solution at  @xmath333 we then interpolating , again with barycentric lagrange interpolation as used in the apparent horizon finder , onto a cubed - sphere grid with excision radius at  @xmath334 with the outer boundary position fixed , and with  @xmath335 with  @xmath336 angular patches , now with  @xmath337 points per cube , naturally again evolving in cartoon mode . in the regauge step we choose here  @xmath338 . this step immediately removes most of the constraint violation from the computational domain , and the largest spatial derivatives , so that the constraint monitor is  @xmath339 as compared to  @xmath340 on the original cubed - ball . this difference seems very troublesome until we take into account that , for example the peak of the kretschmann scalar on the cubed ball grid is  @xmath341 whereas on the cubed sphere it is  @xmath342 . so the reduction in the constraints obviously occurs because we are removing the most extreme part of the domain . note also that our definition of the constraint monitor does not include a normalization by the size of the solution , as in for example  @xcite and subsequent papers . in view of this our reduction in resolution is justified . the evolution then proceeded , now on  @xmath343 cores using  @xmath344 . the regauge safety - net was triggered  @xmath210 times up to  @xmath345 , having fixed  @xmath346 , but the apparent horizon remains on the computational domain throughout the calculation . at  @xmath345 we perform the `` * go - to * '' step of our algorithm again , this time excising at  @xmath347 choosing  @xmath348 . after this the regauge safety - net was not called before  @xmath331 , when we changed cubed - sphere grid once more , keeping the same grid parameters but excising at  @xmath349 , and regauging with  @xmath350 . the evolution continued  @xmath351 , at which time we changed grid for the final time , before which the safety - net was again not called . in the last grid we took the excision radius to be  @xmath352 and regauged with  @xmath353 . after this the regauge safety - net was not called , and the evolution was terminated at  @xmath354 after apparent horizon formation . note that in this evolution the `` * go - to * '' step also employed the phase - in for the generalized harmonic gauge , as described in our single blackhole evolutions in section  [ subsection : phasing ] , taking the same parameters employed in those earlier tests , but now with the initial source functions chosen so that the lapse and shift were frozen as the evolution starts on the new grid . other experiments show that this procedure is not strictly necessary . it may be that some refinement is required to this method to allow the evolution of supercritical data _ indefinitely _ after the collapse , but examining the mass of the apparent horizon , we interpret the solution as having mostly settled down , which should be good enough to diagnose a final mass of the blackhole . centered brill wave . the green planes indicate the times at which the `` * go - to * '' step was applied , and what parameter  @xmath308 was chosen in that procedure . the left plot shows a successful choice , and on the right what happens if this parameter is chosen less carefully . in the upper part of the right hand plot one sees that the horizon contracts , and also sees that the ` regauge ' step is frequently applied , resulting in kinks in the horizon . [ fig : ah_dynamics ] ]    [ [ dynamics - of - the - apparent - horizon ] ] dynamics of the apparent horizon : + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +    in the computation described above , as can be seen in in the left panel of fig . [ fig : ah_dynamics ] , the apparent horizon is always present on the computational domain . the horizon mass initially rapidly grows to a value around  @xmath355 where it remains roughly constant . throughout we see that when the regauge safety - net is triggered a slight oscillation in the horizon mass follows . on the other hand when we change grid we see that the horizon mass exhibits a more prominent kink . in the right - hand panel of fig . [ fig : ah_dynamics ] we plot the apparent horizons obtained when , less - wisely , the parameter  @xmath356 is chosen in the last `` * go - to * '' at  @xmath351 . with this choice the apparent horizon rapidly contracts , although the code fails before it leaves the domain . the safety - net is called ever - more frequently as the method insists on forcing the inner boundary to remain pure outflow , until eventually the code crashes at  @xmath357 . the physical interpretation of this is that the excision boundary is falling off of the domain , which starts to drift outside the blackhole region , and that the safety - net then forces the worldline of the excision boundary to be spacelike . this interpretation would be clearer if we had an event horizon finder , but is given credence by performing evolutions of a schwarzschild blackhole with the  @xmath308 gauge parameter similarly poorly chosen . in such tests we see that the areal radius of the excision boundary can indeed fall outside of the event horizon at  @xmath358 . we have developed a pseudospectral numerical relativity code , ` bamps ` , and in so doing have made a series of improvements and investigations into the approach employed in the spec code . we presented a set of experiments carefully performed so that direct comparison with either published work , or independent computations of the bam finite differencing code could be made . these included evolutions of gauge waves , convergence tests , the use of different constraint damping and ghg formulation parameters , scaling tests , perturbed blackhole evolutions and the treatment of supercritical gravitational waves . ultimately we conclude that the ` bamps ` code is working efficiently , scales as desired up to large numbers of processors , and works on sufficiently general grid setups to evolve initial data of interest . particularly surprising to us was the sensitivity of the method to our modifications of the ghg boundary conditions , even within the class of constraint preserving conditions . this was the case even in our simple evolutions of the schwarzschild spacetime , so it would be very interesting to see the extent to which such results carry over to compact binary evolutions , be it in spec , or in the more distant future in ` bamps ` . from the physics point of view , however , our focus is presently on the collapse of axisymmetric gravitational waves . much of the development reflects this fact . most notably the implementation of octant symmetry with the cartoon method gives orders of magnitude speedups over evolving the same data in full 3d . for a recent complimentary approach see  @xcite . we have additionally written a bespoke axisymmetric apparent horizon finder , which already proved a valuable diagnostic tool , crucial in the evolution of supercritical data , where the existence of an apparent horizon was used as the criterion for moving to an excision grid . naturally further developments to the code may be desirable . for physical interpretation , an event horizon finder would complement our apparent horizon finder . a control system like that of spec  @xcite would be useful in controlling the positions of the apparent horizons . but the highest priority will likely be in generalizing available grid setups to enable dynamical mesh - refinement . we have also considered various different types of axisymmetric moment of time - symmetry gravitational wave initial data . in forthcoming work we use ` bamps ` to evolve this initial data , close to the critical amplitude separating dispersion and collapse to a blackhole . we are grateful to sebastiano bernuzzi , david garfinkle , enno harms , sascha husa , nathan kieran johnson - mcdaniel , harald pfeiffer and hannes rter for interesting discussions . dh would like to express special gratitude to helmut friedrich for interesting discussions and for his warm encouragement . this work was supported in part by the deutsche forschungsgemeinschaft ( dfg ) through its transregional center sfb / tr7 `` gravitational wave astronomy '' , by the dfg research training group 1523/1 `` quantum and gravitational fields '' , and the graduierten - akademie jena . computations were performed primarily at the lrz ( munich ) .
we present a new pseudospectral code , ` bamps ` , for numerical relativity written with the evolution of collapsing gravitational waves in mind . 
 we employ the first order generalized harmonic gauge formulation . the relevant theory is reviewed and the numerical method is critically examined and specialized for the task at hand . in particular 
 we investigate formulation parameters , gauge and constraint preserving boundary conditions well - suited to non - vanishing gauge source functions . 
 different types of axisymmetric twist - free moment of time symmetry gravitational wave initial data are discussed . a treatment of the axisymmetric apparent horizon condition is presented with careful attention to regularity on axis . 
 our apparent horizon finder is then evaluated in a number of test cases . moving on to evolutions , 
 we investigate modifications to the generalized harmonic gauge constraint damping scheme to improve conservation in the strong field regime . 
 we demonstrate strong - scaling of our pseudospectral penalty code . 
 we employ the cartoon method to efficiently evolve axisymmetric data in our  @xmath0 dimensional code . 
 we perform test evolutions of schwarzschild perturbed by gravitational waves and by gauge pulses , both to demonstrate the use of our blackhole excision scheme and for comparison with earlier results . 
 finally numerical evolutions of supercritical brill waves are presented to demonstrate durability of the excision scheme for the dynamical formation of a blackhole .