text:
# Search for Extra Dimensions at LHC

## 1 Introduction

Models with extra dimensions (ED) (see intro ) are very attractive extensions of the Standard Model (SM), in particular with respect to the hierarchy problem. While ED have so far escaped detection limits , they could manifest themselves at LHC via a rich and varied phenomenology.

This review focuses on the three main classes of ED models with prediction at the TeV scale. The aim of the studies was two-fold: establishing the sensitivity to ED signals using detector simulation and various physics backgrounds, as well as assessing whether enough information could be extracted in order to distinguish ED signatures from other new physics signals.

The studies are based on fast simulation tools which describe accurately the expected detector performance. The relevant aspects of the simulation have been validated tdr  in full simulation and with test-beam data whenever possible. Except when stated otherwise, all the results and plots presented here are for an integrated luminosity of 100 fb-1 collected by one of the experiments (i.e. one year at the nominal luminosity of LHC, 1034 cm-2s-1).

## 2 Large Extra Dimensions

In this scenario, the SM fields are confined in our 4D world and only gravity propagates in the bulk. The model is characterized by the number of extra dimensions $`\delta `$ and by the new fundamental scale $`M_D`$. The graviton expands in 4D into a tower of Kaluza-Klein (KK) excitations which couple universally to all SM fields. Even though this coupling is small ($`1/M_{Pl}`$), the large number of states and their small mass splitting lead to sizeable cross-sections at the LHC.

The virtual exchange of KK excitations of graviton can lead to deviations in Drell-Yan cross-sections and asymmetries in SM processes. The left plot of Fig. 1 illustrates such deviations in the $`\gamma \gamma `$ invariant mass distribution virtled . This kind of signatures is clear, very sensitive to new physics and could signal the existence of extra dimensions. However the underlying parameters of the model cannot be extracted in this case because the model is sensitive to unknown ultra-violet physics.

The second class of signatures is the direct production of KK excitations of graviton which will escape detection in 4D: $`q\overline{q}gG^{(k)}`$, $`gqqG^{(k)}`$ and $`gggG^{(k)}`$. In this case, the main signature to look for is some missing energy accompanied by a mono-jet (Fig. 1, right plot) led . Within the allowed region for the effective theory ($`\sqrt{\widehat{s}}<M_D`$), those processes can be reliably calculated and the parameters of the model can be constrained from the measurements. Models with up to four extra dimensions could be probed at LHC. For 100 fb-1, the maximum reach in $`M_D`$ is between 9.1 TeV ($`\delta =2`$) and 6.0 TeV ($`\delta =4`$), corresponding to a radius of compactification between 8 $`\mu `$m and 1 pm. The two parameters of the model can in principle be extracted from the absolute cross-section of the processes or more definitely by collecting $`50`$ fb-1 of data at a different center-of-mass energy.

## 3 TeV-1-sized Extra Dimensions

In this case, the ED are small enough to allow the propagation of gauge bosons in the bulk without contradicting the existing electroweak measurements. Only models with one small ED are considered here, with the fermions localized in the brane at different locations (M1 and M2 models). The main feature of this model is the production of KK excitations of the $`Z`$ and $`W`$ bosons. Their mass spectrum is defined by: $`m_{Z^{(k)},W^{(k)}}^2=m_{Z,W}^2+k^2M_C^2`$ where $`M_C`$ is the compactification scale, known to be $`4`$ TeV from precision electroweak data. Hence only the first excitations $`Z^{(1)}`$ and $`W^{(1)}`$ can be seen at LHC. In the electron channel kkz  the experimental resolution is smaller than the natural width of the $`Z^{(1)}`$. The expected signal 111NB: this plot is at parton level. The simulation study gives similar results. is shown on Fig. 2, left plot. The direct observation of a peak is possible if $`M_C5.8`$ TeV. However, this reach can be improved drastically by using all the information with a maximum likelihood, in particular using the region before the peak (Fig. 2, left) with either constructive or destructive interferences between the $`Z^{(1)}/\gamma ^{(1)}`$ and $`Z/\gamma `$. The sensitivity is thus increased to $`M_C^{max}=9.5`$ TeV and could even reach 13.5 TeV by combining the electron and muon channels for 300 fb-1.

Furthermore, the spin-1 $`Z^{(1)}`$ signal can be distinguished from a spin-2 narrow graviton resonance (section 4) using the angular distribution of its decay products. Thanks to the contributions of the higher lying states, the interference terms and to the additional $`\sqrt{2}`$ factor in its coupling to SM fermions, the $`Z^{(1)}`$ can also be distinguished from a $`Z^{}`$ with SM-like couplings: the distribution of the forward-backward asymmetry for the various cases are shown on Fig. 2 for 4 TeV resonances. The $`Z^{(1)}`$ hypothesis can be discriminated for masses up to about 5 TeV with an integrated luminosity of 300 fb-1.

## 4 Warped Extra Dimensions

We consider here one of the Randall-Sundrum models of warped ED, where gravity propagates in a 5D bulk limited by two branes. The SM fields are confined in the first brane. The metric includes an exponential warp factor which curves the space and connects the $`M_{EW}`$ scale in the SM brane to the Planck scale in the other brane, hence reformulating the hierarchy problem. The phenomenology is defined by two parameters: the scale of physical processes in the SM brane, $`\mathrm{\Lambda }_\pi 1`$ TeV, and the curvature scale $`c=k/M_{Pl}`$. The main feature of this model is the production of KK graviton narrow resonances, whose masses are given by $`m_n=x_n\mathrm{\Lambda }_\pi c`$ where $`x_n`$ are the roots of the $`J_1`$ Bessel function. Here again only the first excitation is likely to be seen at LHC, the other modes being suppressed by the falling parton distribution functions. Some of the best channels are $`G^{(1)}e^+e^{}`$ and $`G^{(1)}\gamma \gamma `$, thanks to the energy and angular resolutions of the LHC detectors. The signal expected with 100 fb-1 in the channel $`e^+e^{}`$ is shown on Fig. 3, left plot, for a pessimistic hypothesis of $`c=0.01`$. In this case, masses up to 2 TeV can be probed rs1 . The reach goes up to about 4 TeV for $`c=0.1`$. The sensitivity is summarized on Fig. 4 in the $`(m_{G^{(1)}},c)`$ plane rs2 : if no signal is found, the area left to the curves labeled ”Discovery” can be excluded at 95% CL. The interesting region in this plane is limited by the constraint $`\mathrm{\Lambda }_\pi 10`$ TeV (otherwise the model would no longer be interesting for solving the hierarchy problem) and by the $`c<0.1`$ limit (low-energy consistency). After one year of data-taking the LHC covers completely this region.

The spin-2 nature of $`G^{(1)}`$ can be measured as well as shown on Fig. 3, right plot. It is worth noting that the acceptance at large pseudo-rapidities ($`1.5<|\eta |<2.5`$) is essential for the spin discrimination. The curves labelled ”Spin-2” on Fig. 4 show that for graviton masses up to 2.3 TeV ($`c=0.1`$), there is a 90% chance that the spin-2 nature of the graviton can be determined with a 95% CL. It has also been shown rs1  that this resonance can be seen in many other channels ($`\mu \mu `$,$`\gamma \gamma `$,$`jj`$,$`b\overline{b}`$, $`t\overline{t}`$,$`WW`$,$`ZZ`$), hence allowing to check the universality of its couplings; and that the size $`R`$ of the ED could also be estimated with a 10% accuracy from the mass and cross-section measurements.

Another intriguing signature of this model is the existence of a new massive scalar, the radion, allowing to stabilize the spacing between the two branes at the distance required for solving the hierarchy problem ($`kR12`$). The radion is very similar to the SM Higgs, and can actually mix with it. Its width is smaller though, as shown on Fig. 5, left plot; and its partial widths can be different, with in particular an enhanced coupling to gluons. There are three additional parameters for the radion: its mass, the vacuum expectation value $`\mathrm{\Lambda }_\varphi `$ and the radion-Higgs mixing parameter $`\xi `$. The results for the SM Higgs have been reinterpreted for the radion case by folding in the new branching ratios and the appropriate detector resolutions radion . Fig. 5 (right plot) shows the expected significances for a radion signal in the $`\gamma \gamma `$ and $`ZZ^{()}`$ channels. Discovery is possible over the whole mass range if $`\mathrm{\Lambda }_\varphi `$ TeV. If its mass permits, the radion can also decay into a pair of Higgs scalars. If $`m_\varphi =300`$ GeV, only 2(4) fb-1 are needed for a $`5\sigma `$ discovery in the very clean channel $`\varphi hh\gamma \gamma b\overline{b}`$ for $`\xi =0(1/6)`$. With 30 fb-1 collected at low luminosity (to maintain excellent $`b`$-tagging), scales up to $`\mathrm{\Lambda }_\varphi =2.2`$ TeV can be probed in this channel, and up to 1.0 TeV in the $`\varphi hhb\overline{b}\tau ^+\tau ^{}`$ for $`m_\varphi =600`$ GeV. Discrimination between the radion and a Higgs scalar will require very precise measurements of their couplings.

## 5 Other models and signatures

The existence of ED could also be probed with the di-jet cross-section laforge  or via the polarization of $`\tau `$ leptons in some models with two Higgs doublets and a singlet neutrino in the bulk ketevi . Finally, the production of black holes at LHC is another striking signature of ED greg .

## 6 Conclusion

The LHC will be able to probe the relevant region of the parameter space for most of the models with ED studied so far. Moreover, in most cases it will be possible to discriminate such signals from other new physics scenarios and infer some information about the underlying model.

## Acknowledgements

The studies presented here have been performed within the ATLAS and CMS collaborations and as such have made use of tools which are the results of collaboration-wide efforts. I would like to thank in particular G. Azuelos, F. Gianotti, I. Hinchliffe, L. Pape, L. Poggioli, S. Shmatov, P. Traczyk and G. Wrochna.


format: text

text:
# Dynamical Domain Wall Fermions

## 1 INTRODUCTION

In the quenched approximation Domain Wall Fermions  (DWF) have been found to be an extremely successful approach to simulating QCD on the lattice. Crucial to this success is the fact that, when working at weak couplings ($`a^12\mathrm{G}\mathrm{e}\mathrm{V}`$) and using improved gauge actions such as the DBW2 action, the degree of explicit chiral symmetry breaking is very small for practically useful sizes of the fifth dimension ($`O(10)`$).

Early simulations of dynamical DWF, performed at relatively coarse couplings, suggested that $`L_sO(100)`$ would be needed before the degree of chiral symmetry breaking was small enough to be acceptable. Here we will report on the progress of a preliminary study  of $`N_f=2`$ dynamical DWF in which we both adopt the DBW2 gauge action, and move to weaker coupling, in an attempt to find a region of parameter space where dynamical DWF simulations are practical.

## 2 SIMULATION PARAMETERS

All the results that will be presented were generated using the DWF action with $`L_s=12`$ and $`M_5=1.8`$, and the DBW2 gauge action with $`\beta =0.80`$ on $`16^3\times 32`$ lattices. Using the HMC algorithm we have generated three separate evolutions for bare masses of $`m_f=0.02`$, $`m_f=0.03`$ and $`m_f=0.04`$. Table 1 summarises the total number of trajectories collected so far, together with the acceptance. Each HMC trajectory is of length 0.5 in HMC time and is split up into 50 leapfrog integration steps for mass of 0.02 and 0.03, and 40 integration steps for 0.04.

## 3 ALGORITHMIC DETAILS

The number of degrees of freedom of DWF grows with $`L_s`$, but the number of physical degrees of freedom does not. To cancel off this bulk divergence a set of Pauli-Villars fields,

$$\mathrm{\Phi }^{}D^{}(m_f=1)D(m_f=1)\mathrm{\Phi },$$
(1)
is added to the DWF Lagrangian for dynamical simulations. Previous work has used two sets of pseudo-fermion fields to represent the DWF action: one for the fermion piece of the action and one for Pauli-Villars. The cancellation between these two terms is therefore only apparent after the average over the pseudo-fermion fields. Here we have used the fact that

$$\frac{\mathrm{det}(D^{}(1)D(1))}{\mathrm{det}(D^{}(m_f)D(m_f))}$$
(2)
is equal to

$$\mathrm{det}(D(1)[D^{}D(m_f)]^1D^{}(1))$$
(3)
to represent the fermion and Pauli-Villars pieces of the action with a single pseudo-fermion field. With this approach the cancellation happens step by step in the leapfrog integration. We find the acceptance of the algorithm is increased by $`1020\%`$ while the inversion costs are reduced by $`2030\%`$ when using this modified force term.

We have also implemented the chronological inverter technique of , leading to a performance improvement of a factor of $`1.7`$. After this improvement calculating a single trajectory takes approximately $`1.6\times 10^4`$, $`9\times 10^3`$ and $`6\times 10^3`$ dirac matrix applications for $`m_f=0.02`$, $`0.03`$ and $`0.04`$ respectively.

## 4 RESULTS

While the lattices collected represent part of a larger RBC collaboration project to calculate many hadronic quantities of phenomenological interest, here we will concentrate on a few mesonic observables to determine the basic properties of our simulations such as scale and quark mass. To calculate these quantities we have used every 50th trajectory, leaving out the first $`600`$ trajectories to allow the evolutions to thermalise. All quoted errors will be from a jackknife estimate of the statistical error.

To quantify the chiral symmetry breaking from finite $`L_s`$, we have measured the residual mass, $`m_{\mathrm{res}}`$ as defined from the breaking term in the Ward-Takahashi identity . To extract this we look at

$$R(t)=\frac{_{x,y}J_{5q}^a(y,t)J_5^a(x,0)}{_{x,y}J_5^a(y,t)J_5^a(x,0)},$$
(4)
which for time greater than some $`t_{min}`$ should be time independent and equal to $`am_{\mathrm{res}}`$ .

Figure 1 shows this for the $`m_f=0.02`$ evolution at the dynamical point. As can be seen, a plateau is evident for $`t2`$, with an error weighted average between timeslice 6 and the end giving a value of 0.00137(2). This number is relatively insensitive to the quark mass with linear extrapolation to $`m_f=0`$ giving $`m_{\mathrm{res}}=0.00136(5)`$.

While this value for the residual mass is relatively small compared to the input quark mass, to properly interpret the value we must know the lattice spacing and mass renormalisation. For the purpose of this preliminary analysis we will assign a single lattice spacing for all three evolutions based on a linear extrapolation in the dynamical quark mass of the $`\rho `$ meson mass.

This is shown in Fig 2, and leads to an inverse lattice spacing of $`1.806(60)\mathrm{GeV}`$.

Figure 3 shows $`M_\pi ^2`$ for degenerate quark masses versus dynamical mass. Fitting to the naive, first order chiral perturbation theory, expectation that

$$M_\pi ^2=B_\pi \left(m_1+m_2\right),$$
(5)
where $`B_\pi `$ is a constant and $`m_1`$ and $`m_2`$ are the quark masses, gives a $`\chi ^2`$ per degree of freedom of 0.2 and results that are consistent with the pion mass vanishing at $`m_f=m_{\mathrm{res}}`$ with $`M_\pi ^2(m_f=m_{\mathrm{res}})=1.6(11)\times 10^4\mathrm{MeV}^2`$. The experimental value of $`M_K^2`$ is shown on the figure as a dotted line. Together with Eq 5 this suggests our lightest quark mass is a little above half the strange quark mass. Going further, following the same approach as , we extract preliminary values of the renormalised light and strange quark masses of $`3.94(31)\mathrm{MeV}`$ and $`103(8)\mathrm{MeV}`$ in the $`\overline{MS}`$-scheme at $`2\mathrm{GeV}`$. An alternative way to get a rough idea of the size our input quark mass in physical terms, which is independent of the way we are setting our scale, is to calculate the ratio $`M_\pi /M_\rho `$. This is 0.541(8), 0.598(8) and 0.637(8) for $`m_f=0.02,0.03`$ and $`0.04`$ respectively.

In quenched simulations using the DBW2 action it was noticed that tunneling between different topological sectors is suppressed with respect to more standard gauge actions . As such, it is important to study how the topological charge varies with trajectory number in our dynamical simulations.

Fig 4 shows this for the $`m_f=0.02`$ evolution, determined using a classically $`O(a^4)`$ improved definition of the topological charge calculated on each lattice after applying 20 steps of APE smearing with a coefficient of 0.45. While it is encouraging that the value of the topological charge is changing, it is clear that, with a separation of 50 trajectories between lattices, strong correlations are present.

## 5 CONCLUSIONS

A preliminary study of $`N_f=2`$ dynamical DWF QCD, using the DBW2 gauge action at an inverse lattice spacing of $`1.8\mathrm{GeV}`$, shows that a regime exists for which the explicit chiral symmetry breaking is small for a computationally practical extent of the fifth dimension.


text:
# References

Addendum to “Groups of Ribbon Knots”

by Ka Yi Ng

The purpose of this document is to clarify the inductive step described in the proof of Theorem 3.2 in . In the second last sentence of the proof, it says, ‘it follows from the inductive proof that $`_n`$ is of index two.’ The question of how this assertion is verified was first raised by Dror Bar Natan and his student Ofer Ron . To avoid any confusions that may arise in the future, the author of the paper  would like to fill in the details to show that $`_n`$ is of index two. Readers of this document are assumed to have familiarity with . For your reference, Theorem 3.2 is stated below.

Theorem 3.2 $`_n`$ forms a subgroup of the free abelian group $`𝒱_n`$ of index two. So its rank is the same as the rank of $`𝒱_n`$ and is the number of linearly independent primitive rational invariants of order $`n`$.

To show that $`_n`$ is of index two, we shall establish the following two lemmas.

Lemma A: Let $`D_1,\mathrm{},D_r`$ be all the distinct non-split chord diagrams of order $`n2`$. Let $`\lambda _1,\mathrm{},\lambda _r`$ be $`r`$ arbitrary integers. Then for $`n>2`$ , there exists a ribbon knot $`K_n`$ whose additive invariants of order $`<n`$ are trivial and $`V(K_n)=_{i=1}^r\lambda _iV(D_i)`$ for all additive invariants $`V`$ of order $`n`$. When $`n=2`$, there exists a ribbon knot $`K_2`$ with $`V(K_2)=2\lambda _1V(D_1)`$ for all additive invariants $`V`$ of order 2.

Lemma B: Given a knot $`K`$ whose additive invariants of order $`<n`$ are trivial, there exist n-chord diagrams $`D_1,\mathrm{},D_r`$ and integers $`\lambda _1,\mathrm{},\lambda _r`$ such that $`V(K)=_{i=1}^r\lambda _iV(D_i)`$ for all additive invariants $`V`$ of order $`n`$.

Lemma A is an extended version of Lemma 3.1 in  and can be proved using the same arguments as given in . Lemma B is proved in . We shall now make use of these two lemmas to prove Lemma C.

Lemma C: Let $`X`$ be the right-handed trefoil. Let $`K`$ be an arbitrary knot and let $`n2`$ be an integer. Define $`s`$ to be 1 if the Arf invariant of the knot $`K`$ is non-zero and 0 otherwise. Then one can construct a ribbon knot $`K_n`$ so that $`V(K\mathrm{\#}sX\mathrm{\#}K_n)=0`$ for all additive invariants $`V`$ of order $`n`$.

Proof of Lemma C: For $`n=2`$ , the knot can always be represented by an integral multiple $`\lambda `$ of the non-split 2-chord diagram $`D`$ so that $`V(K)=\lambda V(D)`$ for all additive invariants of order 2. Write $`\lambda `$ as $`2ps`$ where $`p`$ is an integer. By Lemma A, we can find a ribbon knot $`K_2`$ with $`V(K_2)=2pV(D)`$ for all additive invariants of order 2. The right-handed trefoil knot has $`V(X)=V(D)`$ and thus, the case $`n=2`$ is proved.

Suppose a ribbon knot $`K_n`$ is constructed so that $`V(K\mathrm{\#}sX\mathrm{\#}K_n)=0`$ for all additive invariants of order $`n`$. Lemma B tells us that there exist integers $`\lambda _1,\mathrm{},\lambda _r`$ and $`(n+1)`$-chord diagrams $`D_1,\mathrm{},D_r`$ so that $`V(K\mathrm{\#}sX\mathrm{\#}K_n)=_{i=1}^r\lambda _iV(D_i)`$ for all additive invariants of order $`n+1`$ . Apply Lemma A to construct a ribbon knot $`R`$ whose additive inariants of order $`n`$ are trivial and $`V(R)=_{i=1}^r\lambda _iV(D_i)`$ for all additive invariants of order $`n+1`$. Then $`V(K\mathrm{\#}sX\mathrm{\#}K_n\mathrm{\#}R)=0`$ for all additive invariants of order $`n+1`$.

Lemma C implies that $`V(K\mathrm{\#}sX\mathrm{\#}K_n)=0`$ for all rational-valued invariants of order $`n`$. According to the proof of Theorem 3.2 in , $`_n`$ is a subgroup of $`𝒱_n`$. Thus, as an element in $`𝒱_n`$, we have $`[K]_n[sX]_n+_n`$. This completes the proof that $`_n`$ is of index two.

Ka Yi Ng, e-mail: Kayi.Ng@WallStreetSystems.com


text:
# Evolutions in 3D numerical relativity using fixed mesh refinement

## I Introduction

Currently many researchers are involved in efforts to predict the gravitational waveforms emitted by the inspiral and merger of compact binaries. The direct numerical simulation of these binary sources has been regarded as the best, if not only, way to obtain information about the waves emitted during the merger event itself. Many numerical simulation methods are currently in use, and among the most popular of these is the use of finite difference approximations to Einstein’s equations on a uniform mesh or grid.

It is often the case in simulations of physical systems that the most interesting phenomena may occur in only a subset of the computational domain. In the other regions of the domain it may be possible to use a less accurate approximation, thereby reducing the computational resources required, and still obtain results which are essentially similar to those obtained if no such reduction is made. In particular, we may consider using a computational mesh which is non-uniform in space and time, using a finer mesh resolution in the “interesting” regions where we expect it to be necessary, and using a coarser resolution in other areas. This is what we mean by mesh refinement (MR). In the case of large-scale simulations, where memory and computing time are at a premium, the effective use of MR can allow for simulations achieving effective (fine-grid) resolutions which would be, on a practical level, impossible otherwise. This also implies that some systems requiring a large dynamic range of resolution may only be possible to simulate using mesh refinement. One hopes that the use of fewer total grid points in an efficient mesh refinement application, as opposed to a uniform, single grid or unigrid application, may result in a saving of computing time required to arrive at a solution. It is always the case, however, that the reduction of memory requirements afforded by mesh refinement, especially in 3D, could allow for highly accurate simulations which simply may not fit in the memory of a given computer in a unigrid setting.

In many simulations it may be desirable for the computational mesh to adaptively refine in response to some criterion such as the local truncation error. Such adaptive mesh refinement (AMR) systems provide, ideally, a very high degree of computational efficiency for a given problem. AMR schemes have been in use for decades Berger and Oliger (1984), particularly in the fluid dynamics community Berger and Colella (1989); Plewa (2003). Applications to to astrophysical systems have been underway for quite some time (e.g., Winkler et al. (1984a, b); Plewa and Mueller (2001); Khokhlov (1998)). In numerical relativity, Choptuik Choptuik (1989, 1993) introduced many researchers to AMR through his results using a 1D code for studying critical phenomena. Since the early 90’s efforts have been underway in the numerical relativity community to develop and employ AMR applications in two and three dimensions, in studies of waves Choi (1998); Hawley (2000); Lanfermann (1999); Choi et al. (2003); Papadopoulos et al. (1998); New et al. (2000), critical collapse Choptuik (1993); Liebling (2002); Pretorius (2002); Choptuik et al. (2003), the initial data problem Khokhlov (1998); Diener et al. (2000); Lowe and Brown (2003), inhomogeneous cosmologies Hern (1999), Schwarzschild black holes Brügmann (1996); Lanfermann (1999); Imbiriba et al. (2003), characteristic methods Pretorius and Lehner (2003), and a binary black hole system Brügmann (1999). In several of these cases, the use of mesh refinement provided not only a more efficient means of computing a solution, but even allowed the authors to obtain solutions which would not have been possible via standard unigrid methods, given existing computing technology at the time.

One step in the development of an AMR code is the ability to handle regions of different resolution which are known in advance — fixed mesh refinement (FMR). Thus we view our present work as in some sense a precursor to full 3D AMR simulations. However, the FMR code may be sufficient in its own right, as full AMR may not be necessary for problems where the region of interest is well known beforehand.

An important point should be made regarding the benefits of mesh refinement for reducing storage requirements (and computational costs). It is well known that pseudospectral collocation methods offer exponential convergence as the number of collocation points is increased (for sufficiently smooth functions, which are present in many systems of interest to numerical relativists), and such methods are able to yield extremely high accuracies via simulations which require storage afforded by a typical desktop computer Kidder and Finn (2000); Kidder et al. (2000a); Pfeiffer et al. (2003); Scheel et al. (2002, 2003). Since FMR makes use of finite difference methods, we can obtain at best only polynomial convergence as we increase the resolution of all grids. Our choice of finite difference methods and FMR are based on a choice of physical systems of interest, and on an observation.

One reason we are motivated to use mesh refinement is that we are interested in systems which are “non-smooth” by the standards of pseudospectral methods: hydrodynamics, and gravitational collapse. In hydrodynamics, the formation of discontinuities — shocks — from smooth initial data is a generic phenomena. In gravitational collapse, features may appear on smaller and smaller spatial scales, requiring a means of resolving these features as they appear, either via a truly adaptive algorithm (in which, e.g., the truncation error is used as a refinement criterion Berger and Oliger (1984)) or a “progressive” mesh refinement system in which nested grids are “turned on” at smaller radii and higher resolutions as the simulation proceeds.

The observation mentioned above is that unigrid finite difference codes are generally regarded as being significantly easier to develop than pseudospectral codes. Our interest, from a code-development standpoint, has been to provide an infrastructure whereby existing unigrid codes can be “endowed” with mesh refinement in a way which is somewhat automatic, such that the developer of the original unigrid code is spared the details of implementing mesh refinement. If the introduction of mesh refinement does not significantly alter the dynamics of the system, then one should be able to obtain results comparable to a high resolution unigrid run via an FMR run with appropriately placed fine grids which share the resolution of the unigrid run. In fact, this is the hope of this paper, and the criterion by which we evaluate our results obtained by FMR: ideally, the use of mesh refinement should produce results of comparable quality to a corresponding unigrid run, in terms of stability, accuracy and convergence behaviour. Thus our mesh refinement infrastructure could, in applicable cases, provide a service to the community of researchers who commonly develop unigrid finite difference codes, by providing these researchers a means by which to achieve more accurate results than their current computer allocations allow. It is in this spirit that we are making the FMR system, called Carpet (authored by Erik Schnetter, with refinements offered by several others), freely available as a driver thorn of the open-source Cactus computational infrastructure Allen et al. (2001a, b); Talbot et al. (2000); The Cactus Team ; Goodale et al. (2002).

## II Overview of our FMR method

### II.1 Cactus and mesh refinement

Cactus is an application framework that provides some computational infrastructure such as parallelisation and disk I/O, and lets users specify the modules, called “thorns”, which contain all the “physics”. The main difference between an application framework and a library of subroutines is *control inversion*, meaning that it is Cactus calling the users’ routines, while the main program is part of Cactus. The part of Cactus that controls the creation of initial data and the time stepping is called the *driver*, which is a thorn that interacts with the Cactus scheduler in order to determine which routines are applied to which grid at what time. Control inversion has the important advantage that, by replacing the driver, one can change a Cactus application from unigrid to mesh refinement without rewriting any of the users’ thorns.

In practice, the ability to use the same code with both unigrid and mesh refinement drivers places restrictions on the implementation of a routine. For most evolution thorns these modifications are at most minor. Some thorns, particularly those solving elliptic equations, require more substantial alterations. In most circumstances the restrictions imposed are only technical points, and it should be simple for any new code to be implemented to work with a unigrid or mesh refinement driver interchangeably. Our experience is also that analysis tools for e.g. wave extraction or apparent horizon finding continue to work almost without changes.

### II.2 Mesh refinement method

The mesh refinement driver that we use is called *Carpet* and is available together with the application framework Cactus. It uses the Berger–Oliger approach Berger and Oliger (1984), where the computational domain as well as all refined subdomains consist of a set of rectangular grids. In particular, we base our scheme on the so-called “minimal Berger–Oliger” setup popularised by Choptuik Choptuik (1995). In this simplified version of Berger–Oliger, the grid points are located on a grid with Cartesian topology, and the grid boundaries are aligned with the grid lines. In our version, we also allow fine grid boundaries to occur in between coarse grid points, as shown in Figure 1. (Note that this definition still allows, e.g., spherical coordinate systems.) Furthermore, there is a constant refinement ratio between refinement levels (described below).

We use the following notation. The grids are grouped into refinement levels (or simply “levels”) $`L^k`$, each containing an arbitrary number of grids $`G_{}^{k}{}_{j}{}^{}`$. Each grid on refinement level $`k`$ has the grid spacing (in one dimension) $`\mathrm{\Delta }x^k`$. The grid spacings are related by the relation $`\mathrm{\Delta }x^k=\mathrm{\Delta }x^{k1}/N_{\text{refine}}`$ with the integer *refinement factor* $`N_{\text{refine}}`$. An example is shown in Figure 1. In what follows we will assume that $`N_{\text{refine}}`$ is always set to $`2`$. The base level $`L^0`$ covers the entire domain (typically with a single grid) using a coarse grid spacing. The base level need neither be rectangular nor connected. The refined grids have to be *properly nested*. That is, any grid $`G_{}^{k}{}_{j}{}^{}`$ must be completely contained within the set of grids $`L^{k1}`$ of the next coarser level, except possibly at the outer boundaries.

Although Cactus does allow both vertex and cell centred grids, current relativity thorns only use vertex centring. Hence Carpet currently only supports vertex centred refinement, i.e. coarse grid points coincide with fine grid points (where fine grid points exist).

The times and places where refined grids are created and removed are decided by some *refinement criterion*. The simplest criterion, which is also indispensable for testing, is manually specifying the locations of the refined grids at fixed locations in space at all times. This is called *fixed mesh refinement*. A bit more involved is keeping the same refinement hierarchy, but moving the finer grids according to some knowledge about the simulated system, tracking some feature such as a black hole or a neutron star. This might be called “moving fixed mesh refinement”. Clearly the most desirable strategy is an automatic criterion that estimates the truncation error, and places the refined grids only when and where necessary. This is what is commonly understood by *adaptive mesh refinement*. Carpet supports all of the above in principle, but we will only use fixed mesh refinement in the following. It should be noted that automatic grid placement is a non-trivial problem (see, e.g., Bell et al. (1994); Berger and Rigoutsos (1991)).

### II.3 Time evolution scheme

The time evolution scheme follows that of the Berger and Oliger Berger and Oliger (1984) AMR scheme, in which one evolves coarse grid data forward in time before evolving any data on the finer grids. These evolved coarse grid data can then be used to provide (Dirichlet) boundary conditions for the evolution of data on the finer grids via *prolongation*, i.e. interpolation in time and space. This is illustrated in Figure 2. For hyperbolic systems, where a Courant-like criterion holds, a refinement by a factor of $`N_{\text{refine}}`$ requires time step sizes that are smaller by a factor $`N_{\text{refine}}`$, and hence $`N_{\text{refine}}`$ time steps on level $`k+1`$ are necessary for each time step on level $`k`$. At time steps in which the coarse and fine grids are both defined, the fine grid data are restricted onto the coarse grid (via a simple copy operation) after it has been evolved forward in time. If there are more than two grid levels, then one proceeds recursively from coarsest to finest, evolving data on the coarsest grid first, interpolating this data in time and space along boundaries of finer grids, evolving the finer grid data, and restricting evolved data from finer to coarser grids whenever possible. This is illustrated in Figure 3.

For time evolution schemes that consist only of a single iteration (or step), the fine grid boundary condition needs to be applied only once. Most higher-order time integrations schemes, such as Runge-Kutta or iterative Crank-Nicholson, are actually multi-step schemes and correspondingly require the fine grid boundary condition to be applied multiple times. If this is not done in a consistent manner at each iteration, then the coarse and the fine grid time evolution will not couple correctly, and this can introduce a significant error. We explain this in more detail in Appendix A.

There are several ways to guarantee consistent boundary conditions on fine grids. Our method involves not providing any boundary condition to the individual integration substeps, but instead using a larger fine grid boundary, as demonstrated in Figure 4. That is, each of the integration substeps is formally applied to a progressively smaller domain, and the prolongation operation re-enlarges the domain back to its original size. Note that this “buffering” is done only for prolongation boundaries; outer boundaries are handled in the conventional way. Also, this is done only for the substeps due to the time integration scheme, so that the prolongation *is* applied at fine grid times when there is no corresponding coarse grid time. Note also that the use of buffer zones is potentially more computationally efficient.

We emphasise that the use of these buffer zones is not always necessary. To our knowledge the buffer zones are necessary only when the system of equations contains second spatial derivatives, *and* a multi-step numerical method is used for time integration. This issue arises for the BSSN system discussed below. We also give a simple example using the scalar wave equation in section IV.1 and appendix A.

### II.4 Inter-grid transport operators

As described above, the interaction between the individual refinement levels happens via prolongation and restriction. For prolongation, Carpet currently supports polynomial interpolation, up to quadratic interpolation in time, which requires keeping at least two previous time levels of data. It also supports up to quintic interpolation in space, which requires using at least three ghost zones. We usually use cubic interpolation in space, which requires only two ghost zones. (Quadratic interpolation in space introduces an anisotropy on a vertex centred grid.) For restricting, Carpet currently uses sampling (i.e., a simple copy operation). These transport operators are not conservative. Since our formulation of Einstein’s equation (see below) is not in a conservative form, any use of conservative inter-grid operations offers no benefit. However, the transport operators can easily be changed. (For more discussion of the situations where conservative inter-grid operators are useful or not, see Quirk (1991).)

### II.5 Initial data generation

Initial data generation and time evolution are controlled by the driver Carpet. Initial data are created recursively, starting on the coarsest level $`L^0`$. This happens as follows: On refinement level $`L^k`$, the initial data routines are called. This fills the grids on this level. Then the refinement criterion is evaluated (which might be nothing more than a fixed mesh refinement specification). If necessary, grids on a finer level $`L^{k+1}`$ are created, and initial data are created there, and on all finer levels recursively. Then, the data from level $`L^{k+1}`$ are restricted to level $`L^k`$ to ensure consistency.

In many cases, the initial data specification is only valid for a single time $`t=0`$, such as when using a time-symmetric approach, or when solving an elliptic equation. However, for the time interpolation necessary during prolongation (see above), it may be necessary to have data on several time levels. One solution is to use only lower order interpolation during the first few time steps. We decided instead, according to the Cactus philosophy, that the data that are produced during the initial data creation should in principle be indistinguishable from data produced by a time evolution. Hence we offer the option to evolve coarse grid data *backwards in time* in order to provide sufficient time levels for higher order interpolation in time at fine grid boundaries. This ensures that no special case code is required for the first steps of the time evolution.

This initial data generation proceeds in two stages. First the data are evolved *both forwards and backwards in time* one step, leading to the “hourglass” structure illustrated Figure 5. This evolution proceeds recursively from coarsest to finest, so that all data necessary for time interpolation are present. Note that this would not be the case if we evolved two steps backwards in time, as there would not be enough data for the time interpolation for the restriction operation between these two steps.

In the end we must provide initial data only at times *preceding* the initial time $`t=0`$; i.e., the hourglass structure of Figure 5 is invalid as an initial data specification in Cactus. Therefore we perform in the second stage of this scheme one additional step backwards in time on each level, leading to initial data at the times $`t_0`$, $`t_0\mathrm{\Delta }t^k`$, and $`t_02\mathrm{\Delta }t^k`$ on each level $`L^k`$.

## III Physical System

The set of equations we solve are described in detail in Alcubierre et al. (2003a), and although we briefly review the material here, we suggest interested readers refer to the prior publication. The evolution system is that of Shibata-Nakamura Shibata and Nakamura (1995) and Baumgarte-Shapiro Baumgarte and Shapiro (1999), the so-called BSSN formulation. The physical quantities present in a typical ADM York (1979) evolution are the 3-metric $`\gamma _{ij}`$ and the extrinsic curvature $`K_{ij}`$. In the BSSN formulation, one instead evolves a different set of variables: $`K_{ij}`$ is decomposed into its trace $`K`$ and its trace-free part

$$A_{ij}K_{ij}\frac{1}{3}\gamma _{ij}K,$$
(1)
and one applies a conformal transformation,

$$\stackrel{~}{\gamma }_{ij}=e^{4\varphi }\gamma _{ij}.$$
(2)
We choose $`\varphi `$ such that the determinant of $`\stackrel{~}{\gamma }_{ij}`$, denoted by $`\stackrel{~}{\gamma }`$, is $`1`$.

Thus our evolved quantities are related to the ADM and physical quantities by

$`\varphi `$ $`=`$ $`{\displaystyle \frac{1}{12}}\mathrm{log}(\gamma )`$ (3)
$`\stackrel{~}{\gamma }_{ij}`$ $`=`$ $`e^{4\varphi }\gamma _{ij}`$ (4)
$`K`$ $`=`$ $`\gamma ^{ij}K_{ij}`$ (5)
$`\stackrel{~}{A}_{ij}`$ $`=`$ $`e^{4\varphi }(K_{ij}{\displaystyle \frac{1}{3}}\gamma _{ij}K)`$ (6)
One also creates a new evolved variable $`\stackrel{~}{\mathrm{\Gamma }}^i`$, defined as

$$\stackrel{~}{\mathrm{\Gamma }}^i:=\stackrel{~}{\mathrm{\Gamma }}_{jk}^i\stackrel{~}{\gamma }^{jk}.$$
(7)
The evolution equations for these variables are given by

$`_t\varphi `$ $`=`$ $`{\displaystyle \frac{1}{6}}\alpha K+\beta ^k_k\varphi +{\displaystyle \frac{1}{6}}_k\beta ^k`$ (8)
$`_t\stackrel{~}{\gamma }_{ij}`$ $`=`$ $`2\alpha \stackrel{~}{A}_{ij}+\beta ^k_k\stackrel{~}{\gamma }_{ij}+\stackrel{~}{\gamma }_{ik}_j\beta ^k+\stackrel{~}{\gamma }_{jk}_i\beta ^k{\displaystyle \frac{2}{3}}\stackrel{~}{\gamma }_{ij}_k\beta ^k`$ (9)
$`_tK`$ $`=`$ $`D^iD_i\alpha +\alpha (\stackrel{~}{A}_{ij}\stackrel{~}{A}^{ij}+{\displaystyle \frac{1}{3}}K^2)+\beta ^i_iK`$ (10)
$`_t\stackrel{~}{A}_{ij}`$ $`=`$ $`e^{4\varphi }[D_iD_j\alpha +\alpha R_{ij}]^{TF}+\alpha (K\stackrel{~}{A}_{ij}2\stackrel{~}{A}_{ik}\stackrel{~}{A}_j^k)`$ (11)
$`+\stackrel{~}{A}_{kj}_i\beta ^k+\stackrel{~}{A}_{ki}_j\beta ^k{\displaystyle \frac{2}{3}}\stackrel{~}{A}_{ij}_k\beta ^k`$
$`_t\mathrm{\Gamma }^i`$ $`=`$ $`2_j\alpha \stackrel{~}{A}^{ij}+2\alpha (\stackrel{~}{\mathrm{\Gamma }}_{jk}^i\stackrel{~}{A}^{kj}{\displaystyle \frac{2}{3}}\stackrel{~}{\gamma }^{ij}_jK+6\stackrel{~}{A}^{ij}_j\varphi )`$ (12)
$`_j(\beta ^k_k\stackrel{~}{\gamma }^{ij}\stackrel{~}{\gamma }^{kj}_k\beta ^i\stackrel{~}{\gamma }^{ki}_k\beta ^j+{\displaystyle \frac{2}{3}}\stackrel{~}{\gamma }^{ij}_k\beta ^k),`$
where $`D_i`$ is the covariant derivative corresponding to the 3-metric $`\gamma _{ij}`$, $`R_{ij}`$ is the three dimensional Ricci tensor, and the “$`TF`$” superscript denotes the trace-free part of the enclosed expression.

The gauge conditions are given as follows. The lapse $`\alpha `$ is chosen as one of the Bona-Massó Bona et al. (1995) family of slicings,

$$_t\alpha =\alpha ^2f(\alpha )(KK_0)$$
(13)
where $`K_0K(t=0)`$, and $`f(\alpha )`$ is chosen to give us either harmonic slicing ($`f(\alpha )=1`$) or a “$`1+\mathrm{log}`$” slicing ($`f(\alpha )=2/\alpha `$). In this paper the shift will be held constant at the analytic value in all cases.

We use the code described in Alcubierre and Brügmann (2001); Alcubierre et al. (2001) and refer there for details of the finite differencing scheme used.

For certain problems a small amount of artificial dissipation is useful. We use a dissipation of the form

$`_tu`$ $`=`$ $`ϵh^4(_{xxxx}+_{yyyy}+_{zzzz})u`$ (14)
with the grid spacing $`h`$, which is described in Kreiss and Oliger (1973). Although this dissipation operator is typically employed only in systems with first order derivatives in time and space, we find that its use in the BSSN system (which has first order derivatives in time, and second order derivatives in space) is effective at reducing high-frequency oscillations (i.e., noise) in the simulations, but has little effect on the overall convergence behaviour.

## IV Tests

For simplicity we will present tests using only two levels containing one grid each, which we will refer to as the “coarse grid” and the “fine grid”. The fine grid is a box contained in the larger coarse grid box, with the fine grid having a mesh spacing (in space and time) of half that of the coarse grid. The only limitation on the number of grids is the available computational resources, and we have successfully performed tests with up to 24 levels of refinement.

One of the principal criteria we use to evaluate the effectiveness of the FMR scheme is the requirement of second-order convergence in the limit as the mesh spacing goes to zero. Thus we run a given simulation many times at different resolutions. In our examples we compare against the exact solution and check the convergence of the solution error. We show such tests only for the data on the coarsest grid, because the restriction operator ensures that the coarse grid and fine grid data are identical.

### IV.1 Wave Equation: Periodic Boundaries

We tested our code with a simple wave equation in flat space in Cartesian coordinates using several different kinds of initial data and boundary conditions. The first such test was that of sinusoidal plane waves in a 3D box with periodic boundary conditions. From a code-development standpoint, we simply took an existing set of subroutines for solving the wave equation in parallel and ran using the FMR driver instead of the usual unigrid driver. The formulation of the wave equation we used was a single equation with second order derivatives in both time and space, i.e.

$$_{tt}u=^i_iu,$$
(15)
which we solved using a leapfrog-like scheme. Second order convergence was found, as shown in Figure 6.

An alternative formulation of the wave equation uses second order derivatives in space but only first order in time. We write it in the form

$`_tu`$ $`=`$ $`v`$ (16)
$`_tv`$ $`=`$ $`^i_iu`$ (17)
This formulation is comparable to the ADM (and BSSN) formulations of the Einstein equations. When this formulation of the scalar wave equation is evolved using ICN integration without buffer zones, the result is only first-order convergent, as shown using a one dimensional example in Figure 7. The same formulation evolved with buffer zones converges even at extremely high resolution as shown in Figure 8.

Having demonstrated the need for careful handling of refinement boundaries, and having introduced buffer zones as an effective approach to that, we use buffer zones in all the remaining tests discussed in this paper.

### IV.2 Wave Equation: Gaussian pulse

We consider a Gaussian pulse that crosses a mesh refinement boundary, travelling from the fine into the coarse region. This is supposed to mimic the case of gravitational waves propagating from fine, inner grids and radiating out into coarser grids.

We use an effectively one dimensional domain (planar symmetry in 3D) with $`x[0.5,+0.5]`$, and a coarse grid resolution of $`\mathrm{\Delta }x^0=1/100`$. The region $`x>0`$ is refined by a factor of $`2`$. The Gaussian pulse starts in the refined region and travels to the left. Figure 9 shows the pulse after it has crossed the interface, and compares the result to two unigrid simulations. The errors that are introduced at the refinement boundary are very small and converge to second order. In particular, at this resolution about $`10^3`$ of the original pulse is reflected by the refinement boundary.

### IV.3 Wave Equation: $`1/r`$, with excision

Another test we performed was the solution of the wave equation (15) using initial data which scaled inversely with $`r`$, the distance from the origin. This was conceived as a useful step before moving on to black holes Brügmann (Oct. 2001), as the “puncture data” Brandt and Brügmann (1997) we use for the black holes has elements which scale as $`1/r`$ in the vicinity of the puncture. Other systems with isolated central masses may also be expected to have elements which scale in a similar fashion. The $`1/r`$ data are a static solution of the wave equation in 3D, and are compatible with the standard Sommerfeld outgoing boundary condition. To handle the singularity at the centre, we “excise” the centre of the computational domain by choosing some inner boundary at a finite excision radius from the centre and filling the interior region with prescribed data. Thus in this test we also have a test of the use of FMR in the presence of excision, which is commonly used for the interiors of black holes in analogous simulations. We perform no evolution within this excised region.

We use a full 3D grid with $`x_i[1,+1]`$ and coarse grid resolutions of $`\mathrm{\Delta }x^0=1/32`$ and $`\mathrm{\Delta }x^0=1/64`$. The region $`x_i[0.5,+0.5]`$ is refined by a factor 2, and the region $`|x_i|0.125`$ is excised. Graphs of the error at two different times in the evolution are shown in Figure 10, which also shows corresponding unigrid runs for comparison. We see that the solution is fully convergent, and similar to the corresponding unigrid results in the region of refinement. It is interesting to note that even outside this region, the FMR and unigrid results are very similar for the “transient” shown in the left frame, however the late time (“stationary”) behaviour shown in the right frame reveals a notable difference between the FMR and unigrid results outside the refined grid.

Having demonstrated the existence of convergent solutions of the wave equation for oblique angles of incidence to refinement boundaries (in the $`1/r`$ case), and convergent solutions in which the amplitude of the reflected wave is significantly (roughly three orders of magnitude) less than the amplitude of the transmitted wave, we now move on to full solutions of Einstein’s equations. For more detailed calculations of reflection and transmission effects at refinement boundaries see Lanfermann (1999); Choi et al. (2003).

### IV.4 Robust Stability

We have applied a robust stability test Szilágyi et al. (2000, 2002); Alcubierre et al. (2003b) to the BSSN formulation. This is a test that is meant to supplement and extend an analytic stability analysis, especially in cases where such an analysis is difficult or impossible (because the equations contain too many terms). A numerical test has the advantage that it tests the complete combination of evolution equations, gauge conditions, boundary conditions, as well as the discretisation and the implementation, and (in our case) the mesh refinement scheme. Thus, while a numerical test is not as reliable as an analytically obtained statement, it is able to cover more general cases.

The robust stability test proceeds in three stages of increasing difficulty:

Minkowski (or some other simple) initial data with small random perturbations. The simulation domain is a three dimensional box with periodic boundary conditions. The perturbations should be linear, so we chose a maximum amplitude of $`10^{10}`$. The periodicity means that there is effectively no boundary, so that this stage is independent of the boundary condition. A code is deemed to be robustly stable if it shows at most polynomial growth and if the growth rate is independent of the grid resolution. This is different from other definitions of stability, where exponential growth is often deemed to be stable if the rate of exponential growth is independent of the grid resolution.

The same as stage I, except that the boundaries in the $`x`$ direction are now Dirichlet boundaries. In addition to the noise in the initial data, noise is also applied to these Dirichlet boundaries. This tests the consistency of the formulation with the boundary conditions, but without the complications of edges and corners.

The same as stage II, except that there are now Dirichlet boundary conditions in all directions. This tests whether edges and corners are handled correctly by the combination of the formulation and boundary conditions.

In accordance with Alcubierre et al. (2003b) we chose an effectively one dimensional, planar symmetric domain that extends almost only in the $`x`$-direction with $`x[0.5,+0.5]`$. The domain has only 3 grid points in the $`y`$ and $`z`$ directions. Thus we have to omit stage III of the test here. We used a resolution of $`1/50`$, and refined the centre of the domain. Although the domain was thus essentially one dimensional, the simulation was performed with the full three dimensional code.

Figure 11 shows the $`L_{\mathrm{}}`$-norm of the Hamiltonian constraint versus time for $`1000`$ crossing times for the stages I and II for three different resolutions. All test were run without artificial dissipation. The results show that our implementation is robustly stable.

### IV.5 Gauge Wave

As another test of the BSSN system we implemented the gauge wave of Alcubierre et al. (2003b). From the typical Minkowski coordinates $`\{\widehat{t},\widehat{x},\widehat{y},\widehat{z}\}`$, one defines new coordinates $`\{t,x,y,z\}`$ via the coordinate transformation

$`\widehat{t}`$ $`=`$ $`t{\displaystyle \frac{Ad}{4\pi }}\mathrm{cos}\left({\displaystyle \frac{2\pi (xt)}{d}}\right)`$ (18)
$`\widehat{x}`$ $`=`$ $`x+{\displaystyle \frac{Ad}{4\pi }}\mathrm{cos}\left({\displaystyle \frac{2\pi (xt)}{d}}\right)`$ (19)
$`\widehat{y}`$ $`=`$ $`y`$ (20)
$`\widehat{z}`$ $`=`$ $`z,`$ (21)
where $`d`$ is the size of the simulation domain. In these new variables, the 4-metric is

$$ds^2=Hdt^2+Hdx^2+dy^2dz^2,$$
(22)
where

$$H=H(xt)=1A\mathrm{sin}\left(\frac{2\pi (xt)}{d}\right).$$
(23)
This test provides us with an exact solution to which we can compare our numerical results. In addition to the exact values $`\alpha =\sqrt{H}`$ and $`g_{xx}=H`$, we will compare the extrinsic curvature $`K_{ij}`$, for which the only nonzero component is

$$K_{xx}=\frac{\pi A}{d\sqrt{H}}\mathrm{cos}\left(\frac{2\pi (xt)}{d}\right).$$
(24)
Since $`\beta ^i=0`$ in the analytic solution we do not evolve the shift but keep it set to zero at all times.

For this simulation we find it useful to add dissipation to the evolution equations to suppress high frequency noise at very high resolutions. The reason is evident from unigrid simulations at high resolutions, as demonstrated in Figure 12.

The simulation domain was set up in the same way as in Alcubierre et al. (2003b), which is also the same as is used in Section IV.4. That is, the simulation domain extended almost only in the $`x`$ direction with $`x[0.5,+0.5]`$. The domain has only 3 grid points in the $`y`$ and $`z`$ directions. All boundaries are periodic. Although the domain is thus essentially one dimensional, the simulation was performed with the full three dimensional code. Figure 13 shows the results after $`5`$ crossing times for the metric component $`g_{xx}`$, comparing refinement and unigrid runs. We see perfect second order convergence.

### IV.6 Schwarzschild Black Hole with Excision

Next, we evolved a static Schwarzschild black hole in Kerr-Schild coordinates which are manifestly static. The line element is

$$ds^2=\left(1\frac{2M}{r}\right)dt^2+\left(\frac{4M}{r}\right)dtdr+\left(1+\frac{2M}{r}\right)dr^2+r^2d\mathrm{\Omega }^2$$
(25)
where $`M`$ is the mass of the black hole. We choose to use excision to handle the singularity. We use the excision method together with one of the tests described in Alcubierre and Brügmann (2001). The shift is held static at the analytic solution and the lapse is evolved using $`1+\mathrm{log}`$ slicing. For simplicity we do not search for an apparent horizon but merely excise those points within a cube with corners at $`1M`$. This small excision region removes the divergences due to the singularity while retaining some of the steep gradients. However, as noted by Kidder et al. (2000b); Calabrese et al. (2003), this does not guarantee that the future light cone at the excision boundary is contained within the excised region.

We evolve only one octant of the grid to take advantage of the symmetries present. We use a coarse grid spacing of $`\mathrm{\Delta }x^0=0.4M`$ and $`\mathrm{\Delta }t^0=0.1M`$ with $`29^3`$ points, giving an outer boundary at $`10M`$ (ignoring two symmetry points and two outer boundary points). The fine grid also contains $`29^3`$ points. The points are not staggered about the origin, so there is always a grid point at $`r=0`$ which must be excised.

We compare the runs using mesh refinement with simulations using the unigrid code as described in Alcubierre and Brügmann (2001); Alcubierre et al. (2001). The coarse unigrid test is identical to the coarse grid in the simulation using refinement. That is, $`29^3`$ points are used with a grid spacing of $`\mathrm{\Delta }x=0.4M`$ giving an outer boundary at $`10M`$. For the medium unigrid run the resolution is doubled whilst the outer boundary location is held fixed, giving the same effective resolution near the excision region as the simulation using refinement. That is, $`55^3`$ points are used with a grid spacing of $`\mathrm{\Delta }x=0.2M`$. To compare with the high resolution run using refinement we also perform a unigrid run with the same effective resolution. That is, $`105^3`$ points are used with a grid spacing of $`\mathrm{\Delta }x=0.1M`$.

Results showing the norm of the change of the lapse in time are shown in Figure 14. Firstly we note that the use of refinement combined with an excision boundary has no qualitative effect on the simulation. As in Alcubierre and Brügmann (2001) the change in the lapse shows an exponentially damped oscillation in the absence of instabilities.

However, in this case we do see an exponentially growing instability which sets in only at the highest resolutions. By tracing back the magnitude of this mode we see that it appears to come from floating point round-off error at the initial time. It appears in runs both with and without mesh refinement and the growth rate is the same in both cases. The origin of this instability is unclear, especially as very similar simulations with the same resolutions were shown to be stable in Alcubierre and Brügmann (2001). However, the important point for this paper is that the stability of the simulations are independent of the mesh refinement.

The convergence of the Hamiltonian constraint at $`t=500M`$ is shown in figure 15, which should be compared with Figure 3 in Alcubierre and Brügmann (2001). Second order convergence for the unigrid and mesh refinement simulations are clear away from the excision boundary, whilst at the excision boundary the convergence is not so obvious due to the low resolution. The error in the mesh refinement runs is comparable to the unigrid runs with the same effective resolution, although only near the excision boundary are the full benefits seen. Given the computational resources required, illustrated in Table 1, the benefits of mesh refinement are clear.

## V Conclusions

We have shown comparisons of results obtained via fixed mesh refinement (FMR) techniques and via standard unigrid methods. Ideally, the appropriate use of FMR should produce results very similar in accuracy to the corresponding unigrid results, while retaining the stability and convergence properties of the unigrid simulations. “Corresponding” in this case implies that the unigrid simulation has the same resolution as the finest grid of the FMR simulation. This ideal is only expected to hold for simulations where most of the “interesting behaviour” occurs in a limited region which is covered by refined grids, such as in simulations of compact objects.

We find that particular care must be taken at boundaries between coarse and fine grids, in order to retain stability and convergence. In cases where the system contains second spatial derivatives and a multistep time integration method is used, we introduce buffer zones to guarantee the convergence of the mesh refinement implementation. The need for these buffer zones can be clearly seen in very fine resolution simulations of the wave equation when written in “mixed” form (with second spatial derivatives and first derivatives in time). The same effects are seen in simulations of gauge waves even at intermediate resolutions when using the BSSN formulation in a similar mixed form.

In order to use parabolic interpolation in time, even from the initial time, we designed a scheme by which data is initially evolved forwards and backwards in time. This means that the fine grid evolution can begin with three time levels of coarse grid data with which to do the time interpolation.

We are able to obtain results from FMR simulations which compare favourably with corresponding unigrid simulations for the following systems: the wave equation with $`1/r`$ data and excision near the origin, robust stability tests, a nonlinear gauge wave, and a Schwarzschild black hole with excision. The robust stability tests indicate that numerical noise combined with inter-grid transport operations will not lead to exponential blow ups, and that the introduction of refinement boundaries does not produce instabilities. Both the wave equation test and the gauge wave test show that proper second order convergence can be obtained irrespective of waves propagating through the interfaces between coarse and fine grids. For the case of the Schwarzschild black hole, although we do see an instability at high resolutions, this instability is present in the unigrid system as well. At lower resolutions, we see stable behaviour showing close agreement between FMR and unigrid results.

There are few other implementations of mesh refinement in more than one dimension in numerical relativity with which to compare. The method of Choptuik et al. (2003) is very close to ours. However, it is an axisymmetric code designed to evolve relativistic scalar fields. The code presented in Brügmann (1999) uses the ADM formalism (which is known to be unstable) and exhibits problems at the refinement boundaries. In the implementation of Choi et al. (2003) which uses Paramesh (MacNeice et al. (2000)), all refinement levels are evolved using the same timestep size. This avoids the difficulties that we observe in obtaining convergence at refinement boundaries. However it may be less efficient, especially where a large range of refinement levels is employed. Our implementation provides both the full efficiency of the Berger–Oliger method and is a generic interface allowing the simple implementation of other formulations and systems such as relativistic hydrodynamics.

Finally, we note that the evolutions described in this paper were performed by taking routines for unigrid simulations of each physical system of interest, and after only slight modifications to these routines, the original unigrid driver was exchanged for the “Carpet” FMR driver. The application of FMR techniques to existing unigrid systems is thus something which, from a code development viewpoint, can be performed generally. A clear advantage of this approach is that most existing analysis tools (such as for wave extraction or apparent horizon finding) will continue to work. Thus we invite other researchers to make use of our freely available code to perform their finite-difference-based simulations in an FMR setting, and thus achieve higher accuracies with less computing resources than their current unigrid simulations may require.

###### Acknowledgements.

We thank Bernd Brügmann and Nils Dorband for their early contributions to the work on the wave equation in FMR, Tom Goodale, Gabrielle Allen and the Cactus team for their assistance regarding the infrastructure development of Cactus which makes FMR possible, and Edward Seidel for much encouragement over the past years. We also thank our collaborators in Wai-Mo Suen’s WuGrav group, with which we shared an early version of the code. Finally we thank Richard Matzner, Bernd Brügmann, and Denis Pollney for their helpful comments regarding the preparation of this manuscript. Work on this project was funded by the AEI, the DFG’s Collaborative Research Centre SFB 382, and by NSF grant PHY 0102204.

## Appendix A Coupling coarse and fine grid time evolution

When certain formulations of the equations are evolved with multi-step time integration methods, we find that the standard Berger–Oliger approach leads to a loss of convergence at the boundaries of fine grids. We will consider only second order methods here. We will argue that the use of Berger–Oliger style prolongation to find boundary conditions for refined grids at the intermediate steps causes a second order truncation error after a fixed number of time steps. Therefore the global error will be first order for a hyperbolic system where to integrate to a fixed time the number of timesteps increases with the grid resolution.

We first give an example, considering a simple one dimensional example on an infinite domain. The wave equation

$`_t^2u`$ $`=`$ $`_x^2u`$ (26)
can be reformulated into a system that is first order in time, but still second order in space:

$`_tu`$ $`=`$ $`v`$ (27)
$`_tv`$ $`=`$ $`_x^2u\text{.}`$ (28)
Discretising this system on a uniform grid in space and time, we represent the function $`u`$ and $`v`$ through

$`u(t,x)`$ $`=`$ $`U_i^n`$ (29)
$`v(t,x)`$ $`=`$ $`V_i^n`$ (30)
at the collocation points $`t=n\mathrm{\Delta }t`$, $`x=i\mathrm{\Delta }x`$, where $`n`$ and $`i`$ number the grid points in the temporal and spatial direction, and $`\mathrm{\Delta }t`$ and $`\mathrm{\Delta }x`$ are the temporal and spatial grid spacings, respectively. The ratio $`\alpha =\mathrm{\Delta }t/\mathrm{\Delta }x`$ is the CFL factor. $`h`$ will be used to refer to terms that are order $`\mathrm{\Delta }x`$ or $`\mathrm{\Delta }t`$.

For the spatial discretisation of the operator $`_x^2`$ we will consider the second order centred differencing stencil

$`D_x^2F(x)`$ $`:=`$ $`{\displaystyle \frac{1}{(\mathrm{\Delta }x)^2}}\left(F(x\mathrm{\Delta }x)2F(x)+F(x+\mathrm{\Delta }x)\right)\text{.}`$ (31)
For time integration we consider one form of the ICN (iterative Crank-Nicholson) method. It consists first of an Euler step

$`T_E[F_0]`$ $`:=`$ $`F_0+\mathrm{\Delta }t_tF_0`$ (32)
and then several (in our case, two) ICN iterations

$`T_I[F_0,F]`$ $`:=`$ $`F_0+\mathrm{\Delta }t_t{\displaystyle \frac{1}{2}}\left(F_0+F\right)`$ (33)
so that the overall ICN integration with $`2`$ iterations is given by

$`T_{\mathrm{ICN}(2)}[F_0]`$ $`:=`$ $`T_I[F_0,T_I[F_0,T_E[F_0]]].`$ (34)
Note that there exist several slightly different variants of the ICN method.

This discretisation is correct up to second order, which means that for a fixed CFL factor $`\alpha `$, the discretisation error is $`O(h^2)`$. This means also that functions $`f(t,x)`$, which depend on $`t`$ and $`x`$ in a polynomial manner with an order that is not higher than quadratic, are represented exactly, i.e. without discretisation error. (Such functions $`f`$ can be written as

$`f(t,x)`$ $`:=`$ $`{\displaystyle \underset{p=0}{\overset{2}{}}}{\displaystyle \underset{q=0}{\overset{2}{}}}C_{pq}t^px^q`$ (35)
with constant coefficients $`C_{pq}`$.)

### A.1 Performing a single step

It is illustrating to perform a time integration step by step. We will use for that the solution of the wave equation

$`u(t,x)`$ $`=`$ $`{\displaystyle \frac{1}{2}}t^2+{\displaystyle \frac{1}{2}}x^2`$ (36)
$`v(t,x)`$ $`=`$ $`t`$ (37)
which will test whether the formulation is indeed capable of representing all functions of the form (35).

Starting from the exact solution at time $`t_0`$,

$`U(t_0,x)`$ $`=`$ $`{\displaystyle \frac{1}{2}}t_0^2+{\displaystyle \frac{1}{2}}x^2`$ (38)
$`V(t_0,x)`$ $`=`$ $`t_0,`$ (39)
the result of the Euler step of the ICN integration is

$`U^{(0)}(t_0+\mathrm{\Delta }t,x)`$ $`=`$ $`U(t_0,x)+\mathrm{\Delta }tV(t_0,x)`$ (40)
$`=`$ $`{\displaystyle \frac{1}{2}}t_0^2+{\displaystyle \frac{1}{2}}x^2+(\mathrm{\Delta }t)t_0`$ (41)
$`V^{(0)}(t_0+\mathrm{\Delta }t,x)`$ $`=`$ $`V(t_0,x)+\mathrm{\Delta }tD_x^2U(t_0,x)`$ (42)
$`=`$ $`t_0+\mathrm{\Delta }t`$ (43)
where $`F^{(k)}`$ denotes the result after $`k`$ ICN iterations. In these expressions, $`U^{(0)}(t_0+\mathrm{\Delta }t,x)`$ differs from the true solution $`U(t_0+h,x)`$ by a term $`O(h^2)`$.

The first ICN iteration then leads to

$`U^{(1)}(t_0+\mathrm{\Delta }t,x)`$ $`=`$ $`U(t_0,x)+\mathrm{\Delta }t{\displaystyle \frac{1}{2}}\left(V(t_0,x)+V^{(0)}(t_0+\mathrm{\Delta }t,x)\right)`$ (44)
$`=`$ $`{\displaystyle \frac{1}{2}}t_0^2+{\displaystyle \frac{1}{2}}x^2+\mathrm{\Delta }tt_0+{\displaystyle \frac{1}{2}}(\mathrm{\Delta }t)^2`$ (45)
$`=`$ $`{\displaystyle \frac{1}{2}}(t_0+\mathrm{\Delta }t)^2+{\displaystyle \frac{1}{2}}x^2`$ (46)
$`V^{(1)}(t_0+\mathrm{\Delta }t,x)`$ $`=`$ $`t_0+\mathrm{\Delta }t{\displaystyle \frac{1}{2}}\left(D_x^2U(t_0,x)+D_x^2U^{(0)}(t_0+\mathrm{\Delta }t,x)\right)`$ (47)
$`=`$ $`t_0+\mathrm{\Delta }t`$ (48)
which is already the correct result. The second ICN iteration does not change the above expressions.

### A.2 Mesh refinement

Let us now introduce mesh refinement, so that the spatial and temporal resolution is not uniform any more. Let the grid points be staggered about the origin, so that the grid points with $`x<0`$ have a coarse spatial grid spacing $`2\mathrm{\Delta }x`$, and the grid points with $`x>0`$ have a fine spatial grid spacing $`\mathrm{\Delta }x`$. Let the CFL factor $`\alpha `$ remain uniform, so that the temporal grid spacing is correspondingly $`2\mathrm{\Delta }t`$ for $`x<0`$ and $`\mathrm{\Delta }t`$ for $`x>0`$. We will use the Berger–Oliger time stepping scheme, which is explained in Section II.3 and illustrated in Figure 3 of the main text.

The discretisation of $`_x^2`$ needs so-called ghost points on the fine grid, which are filled by prolongation. Our third order polynomial spatial prolongation operator is given by

$`P_S[F](x)`$ $`:=`$ $`\{\begin{array}{cc}F(x)\hfill & \text{when aligned}\hfill \\ \frac{1}{16}[1,+9,+9,1]\hfill & \\ F(x+[3,1,+1,+3]\mathrm{\Delta }x)\hfill & \text{otherwise}\hfill \end{array}`$ (52)
and the second order polynomial temporal prolongation operator by

$`P_T[F](x)`$ $`:=`$ $`\{\begin{array}{cc}F(t)\hfill & \text{when aligned}\hfill \\ \frac{1}{8}[1,+6,+3]\hfill & \\ F(t+[3,1,+1]\mathrm{\Delta }t)\hfill & \text{otherwise}\hfill \end{array}`$ (56)
where $`[\mathrm{}]`$ denotes a vector, and the operator $`()`$ a dot product.

Let us now consider a time evolution of the above solution (3637). According to the Berger–Oliger time stepping scheme, the coarse grid evolution happens as on a uniform grid. The fine grid can now be evolved in two different ways: (a) with a Dirichlet boundary condition from prolongating from the coarse grid, or (b) without a boundary condition, i.e. only as an IVP (as opposed to an IBVP). In case (b), the “lost” points have to be filled by prolongation after the time step; in that case, prolongation therefore has to fill not $`1`$, but $`k+1`$ grid points for an ICN integration with $`k`$ iterations. This is illustrated in Figure 4 of the main text.

### A.3 ICN integration with prolongation

Let us examine case (a) in more detail, in which case the fine grid boundary values at $`x<0`$ are given through prolongation. For the solution (3637), time integration and prolongation are exact, hence it is for $`x<0`$

$`U^{(k)}(t_0+\mathrm{\Delta }t,x)`$ $`=`$ $`{\displaystyle \frac{1}{2}}(t_0+\mathrm{\Delta }t)^2+{\displaystyle \frac{1}{2}}x^2`$ (57)
$`V^{(k)}(t_0+\mathrm{\Delta }t,x)`$ $`=`$ $`t_0+\mathrm{\Delta }t`$ (58)
for all ICN iterations $`k`$. (This assumes that all ICN iterations end up the final time. There are also different ways of expressing ICN.)

The Euler step then leads to, for $`x>0`$,

$`U^{(0)}(t_0+\mathrm{\Delta }t,x)`$ $`=`$ $`{\displaystyle \frac{1}{2}}t_0^2+{\displaystyle \frac{1}{2}}x^2+\mathrm{\Delta }tt_0`$ (59)
$`V^{(0)}(t_0+\mathrm{\Delta }t,x)`$ $`=`$ $`t_0+\mathrm{\Delta }t`$ (60)
as in (4043) above. For all values of $`x`$, this can be written as

$`U^{(0)}(t_0+\mathrm{\Delta }t,x)`$ $`=`$ $`{\displaystyle \frac{1}{2}}(t_0+\mathrm{\Delta }t)^2+{\displaystyle \frac{1}{2}}x^2\mathrm{\Theta }(x>0){\displaystyle \frac{1}{2}}(\mathrm{\Delta }x)^2`$ (61)
$`V^{(0)}(t_0+\mathrm{\Delta }t,x)`$ $`=`$ $`t_0+\mathrm{\Delta }t`$ (62)
by using the $`\mathrm{\Theta }`$ function, which is defined as

$`\mathrm{\Theta }(L)`$ $`=`$ $`\{\begin{array}{cc}0\hfill & \text{when }L\text{ is false}\hfill \\ 1\hfill & \text{when }L\text{ is true}\text{.}\hfill \end{array}`$ (65)
The first ICN iteration leads to

$`U^{(1)}(t_0+\mathrm{\Delta }t,x)`$ $`=`$ $`{\displaystyle \frac{1}{2}}(t_0+\mathrm{\Delta }t)^2+{\displaystyle \frac{1}{2}}x^2`$ (66)
$`V^{(1)}(t_0+\mathrm{\Delta }t,x)`$ $`=`$ $`t_0+\mathrm{\Delta }t+\mathrm{\Theta }(0<x<H){\displaystyle \frac{1}{2}}\alpha ^2\mathrm{\Delta }t\text{.}`$ (67)
On a uniform grid, and for our solution $`u`$, the first ICN iteration is already exact. In general, the first ICN iteration should lead to an error which is $`O(h^2)`$. This is here not the case; clearly $`V^{(1)}`$ has an $`O(h)`$ error. However, the error is confined to a region of length $`\mathrm{\Delta }x`$, so an integral norm of the error can be $`O(h^2)`$; it is not clear what this means in practice.

The result of the second ICN iteration is

$`U^{(2)}(t_0+\mathrm{\Delta }t,x)`$ $`=`$ $`{\displaystyle \frac{1}{2}}(t_0+\mathrm{\Delta }t)^2+{\displaystyle \frac{1}{2}}x^2+\mathrm{\Theta }(0<x<H){\displaystyle \frac{1}{4}}\alpha ^2(\mathrm{\Delta }t)^2`$ (68)
$`V^{(2)}(t_0+\mathrm{\Delta }t,x)`$ $`=`$ $`t_0+\mathrm{\Delta }t\text{.}`$ (69)
This expression is the final result of the first fine grid time step. It has an error that is $`O(h^2)`$, or $`O(h^3)`$ in an integral norm. As noted above a local error that is $`O(h^2)`$ will lead to a global error of $`O(h)`$ when considered at a fixed time. That is, the method is only first order convergent.

A third ICN step — that we don’t perform, but it would be possible to do so — would result in

$`U^{(3)}(t_0+\mathrm{\Delta }t,x)`$ $`=`$ $`{\displaystyle \frac{1}{2}}(t_0+\mathrm{\Delta }t)^2+{\displaystyle \frac{1}{2}}x^2`$ (70)
$`V^{(3)}(t_0+2\mathrm{\Delta }t,x)`$ $`=`$ $`t_0+\mathrm{\Delta }t`$
$`+\left[2\mathrm{\Theta }(0<x<\mathrm{\Delta }x)+\mathrm{\Theta }(\mathrm{\Delta }x<x<2\mathrm{\Delta }x)\right]{\displaystyle \frac{1}{8}}\alpha ^4\mathrm{\Delta }t`$
which is worse than the result of the previous iteration. The error is again $`O(h)`$, but has now “infected” two grid points. The error is smaller by a factor $`\alpha ^2`$, but the error at the two grid points has a different sign, indicating the start of an oscillation.

It should be noted that using higher order derivatives, e.g. using a five-point stencil for a fourth-order second derivative, does not remove this error. Similarly, using a higher order prolongation operator does not help. (We assume that it would be possible to adapt the prolongation scheme directly to the time integration scheme and arrive at a consistent discretisation in this way. We do not think that this is worthwhile in practice.)

The main problem seems to be caused by taking a second derivative, which has formally an $`O(1)`$ error. The usual arguments why the error should be smaller in practice Choptuik (1991) do not hold near the discontinuity that is introduced at the refinement boundary.

The result of simulating case (a) is shown in Figure 7 of the main text. For high resolutions, convergence seems to be only of first order. Figure 8 shows case (b), which shows second order convergence for all resolutions.

### A.4 The general case

In the general case the coarse grid evolution and the prolongation will not provide exact data at the boundary on the fine grid. Instead we argue heuristically as follows. The initial data $`U^{(0)},V^{(0)}`$ will be correct to order $`h^2`$,

$$U^n=U(t^n,x)+h^2e^nf^n(x)+O(h^3)$$
(72)
where $`e^n`$ is a constant and $`f^n`$ a smooth function of $`x`$. This data will be evolved forward on the coarse grid using ICN to time $`t^{n+1}`$ where it will also be correct to order $`h^2`$,

$$U^{n+1}=U(t^{n+1},x)+h^2e^{n+1}f^{n+1}(x)+O(h^3)$$
(73)
The fine grid data will then be evolved. The first ICN step is an Euler step to $`t^{n+1/2}`$ (as the fine grid timestep is one half the coarse grid timestep). On the interior of the refined grid the result will be first order accurate in time and second order in space,

$$U^{n+1/2}=U(t^{n+1/2},x)+(\mathrm{\Delta }te_t^{n+1/2}+(\mathrm{\Delta }t)^2e_x^{n+1/2})f^{n+1/2}(x)+O(h^3).$$
(74)
The boundary data for the fine grid is then found by prolongation. Interpolation in space and time, assumed to be second order accurate, gives

$$U^{n+1/2}(x<0)=U(t^{n+1/2},x)+h^2e_p^{n+1/2}f^{n+1/2}(x)+O(h^3).$$
(75)
Therefore at the refinement boundary there is a discontinuous jump in the function,

$$\left[U^{n+1/2}\right]_{x=0}^{x=0+}=\left(he_t^{n+1/2}+(\mathrm{\Delta }t)^2\left(e_x^{n+1/2}e_p^{n+1/2}\right)\right)f^{n+1/2}(x).$$
(76)
As shown in the case of simple polynomial data in appendix A.3, the next step in the ICN loop will lead to a second order local error when the second derivative of the function is taken at the discontinuity. Thus the error at a fixed time will only be first order convergent. It is clear that this discontinuity will also lead to first order errors with other multi-step time integration methods such as the Runge-Kutta methods.

We have performed a symbolic calculation of two complete coarse grid time steps with generic initial data using Maple Maple . We find that the errors in $`u`$ and $`v`$ scale with $`\mathrm{\Delta }x^2`$ and $`\mathrm{\Delta }x`$, respectively, both after the first and the second coarse grid steps. This is consistent with the calculation above as well as our numerical results. This means that, after a fixed time $`t`$, the error in $`u`$ is $`O(\mathrm{\Delta }x)`$, so that we expect first order convergence only.


text:
# STM characterization of the Si-P heterodimer

## I Introduction

Currently there are several proposals to use phosphine gas (PH3) to create atomic-scale devices by positioning P atoms on a surface using H-lithography. Tucker and Shen (2000); O’Brien et al. (2001); Schofield et al. (2003) One such idea O’Brien et al. (2001) involves the fabrication of a 2-dimensional array of P atoms in Si for the realization of a silicon based quantum computer. Kane (1998) A requirement of each of these schemes is that the P atoms, relinquished from the PH3 molecules, are incorporated into substitutional lattice sites. Since PH3 is used for doping Si in the semiconductor industry, the interaction of PH3 with the Si(001) surface has been studied for over two decades using a wide variety of surface sensitive techniques (e.g., Ref. Yu and Meyerson, 1984). It is well established that the adsorption of PH3 on the Si(001) surface is partially dissociative, resulting in the adsorption of both PH3 and PH$`_2.`$ Colaianni et al. (1994); Yoo et al. (1995); Hamers et al. (1996); Shan et al. (1996); Lin et al. (1999); Schofield et al. (2003) Moreover, it is known that annealing PH3 dosed surfaces causes the complete dissociation of the adsorbed PH3 and PH2 molecules (estimates of the temperature at which this occurs vary within the range 300C - 425Shan et al. (1996); Hamers et al. (1996); Lin et al. (2000); Colaianni et al. (1994)) resulting in the incorporation of P into the surface to form Si-P heterodimers. Wang et al. (1994a); Hamers et al. (1996); Oberbeck et al. (2002); Schofield et al. (2003) The majority of the investigations mentioned above have studied the interaction of phosphine with the Si(001) surface at high surface coverages. However, at low coverages there are only two papers reporting detailed atomic-resolution scanning tunneling microscopy (STM) investigations of the Si(001):PH3 system, the first by Wang et al. Wang et al. (1994b) and the second by Kipp et al. Kipp et al. (1995) While the experimental results presented in these two papers are in good agreement, their interpretations are contradictory.

The experimental results from Wang et al. Wang et al. (1994b) and Kipp et al. Kipp et al. (1995) showed that after the Si(001) surface was exposed to a low dose of PH3, two types of bright features appeared on the surface. Both features were centered on top of the Si dimer rows, with one feature larger than the other. Wang et al. Wang et al. (1994b) also observed an increase in the number of dimer vacancies on the surface after PH3 dosing. In both studies, after the dosed surfaces were annealed (to 230C, Ref Wang et al., 1994b; and to 350C, Ref Kipp et al., 1995), the two types of bright feature disappeared from the surface and were replaced by 1-dimensional (1-D) chains running perpendicular to the underlying Si dimer rows and a few small 2-dimensional (2-D) islands. The above data was interpreted by Wang et al. Wang et al. (1994b) in the following way: the larger bright feature was a single Si atom (a monomer) ejected from the surface and the smaller bright feature was a PH3 molecule adsorbed on top of a Si dimer. Upon heating the surface to 230C the Si monomers combined to form 1-D rows and small 2-D islands. This annealing also induced complete PH3 dissociation, with the P atoms being incorporated into the surface and the H adsorbed onto the surface to form monohydride dimers. Kipp et al. Kipp et al. (1995) presented a different interpretation of the above data. They concluded that the large bright features were PHn (n=1-3) molecules adsorbed at defects and the small bright features were P-P dimers, resulting from the pairing and subsequent dissociation of two PH3 molecules (with the dissociated hydrogen atoms also bonding to the surface). Upon heating the surface to 350C Kipp et al. Kipp et al. (1995) concluded that the P-P dimers coalesce to form small chains and islands. They suggested that the formation of large islands was hindered by the presence of adsorbed H on the surface.

In contrast to the results of Wang et al. and Kipp et al., separate Fourier transform infrared spectroscopy studies Shan et al. (1996) and a combined electron energy loss spectroscopy and thermal programmed desorption study Colaianni et al. (1994) have indicated that molecular PH3 adsorption is accompanied by a partial dissociation of PH3 to form adsorbed PH2 and H, during phosphine adsorption at room temperature. Therefore, further study of this system is required to consolidate the different interpretations. We have used high resolution STM and Auger electron spectroscopy to study the evolution of a PH3 dosed Si(001) surface, as a function of temperature and coverage. Particular emphasis was placed on the process of P incorporation into the Si(001) surface and the formation of the Si-P heterodimer. To this end we have also performed a study of surface segregated P from a Si encapsulated $`\delta `$-doped layer. Pivotal to these experiments is the ability to resolve the presence of a single P atom in the Si(001) surface using STM. We show that the isolated Si-P heterodimer appears as a zig-zag feature in filled state STM images, but that the heterodimer is not visible if there is a moderate coverage of monohydrides on the surface. However, empty state imaging can be used to identify the heterodimers at all coverages.

## II Experimental

Experiments were performed in ultra-high vacuum (base pressure $`<5\times 10^{11}`$ mbar) using an Omicron variable temperature STM and electrochemically etched tungsten tips. The n-type Si(001) sample had a resistivity of $``$ 1 $`\mathrm{\Omega }`$cm and was prepared in vacuum by out-gassing overnight at 575C using a resistive heater element, followed by flashing to 1175C by passing a direct current through the sample. After flashing, the samples were cooled slowly ($``$3C/s) from 900C to room temperature. For all anneal steps the sample temperature was maintained at the stated anneal temperature for 10 seconds unless stated otherwise. The PH3 dosing was performed by opening a precision leak valve, which faced the sample, between the UHV chamber and a PH3 micro-dosing system. All doses were performed by back-filling the chamber with a 1 $`\times `$ 10-9 Torr pressure of PH3 and controlling the dose in Langmuirs (1 L = 106 $`\times `$ pressure (torr) $`\times `$ dose (seconds)) via the dose time. The ionization gauge was not in line-of-sight of the sample. All STM images were taken at room temperature and are filled state images (negative sample bias) unless stated otherwise. Dual bias STM images, where filled and empty state images are obtained simultaneously, were taken with a negative (positive) sample bias while scanning the tip in forward (reverse) direction.

## III Results and discussion

### III.1 Identification of surface features for room temperature adsorption of PH3 on Si(001)

Figure 1 shows STM images of the Si(001) surface after exposure to low coverages of phosphine (PH3) gas, at room temperature. The clean Si(001) surface consists of rows of Si dimers which are seen as the parallel lines running diagonally across the images.

Figure 1(a) and (b) show a 21 nm $`\times `$ 17 nm area of a Si(001) surface before and after dosing with 0.002 L of PH3, respectively. The circled features in (a) are single dimer vacancy defects which are inherent to the clean surface and which we identify to distinguish them from species resulting from PH3 dosing. After PH3 dosing \[Fig. 1(b)\] we find that additional features have appeared on the surface, which are: (i) adsorbed PH3, (ii) adsorbed PH2 \+ H, (iii) dimer vacancies (in addition to those present before dosing) and (iv) Si ad-dimers. Figure 1(c) shows an image of a PH3 molecule which has a height/diameter of 0.9 Å/11 Å in filled state images and is adsorbed upon the center of a dimer row, as previously identified by Wang, Shan and coworkers. Wang et al. (1994b); Shan et al. (1996) It can be seen in Fig. 1(c) that this molecule is adsorbed adjacent to a defect, which we find is commonly the case at low coverage. We believe that defects stabilize the adsorption of the PH3 molecules centered over Si dimers. However we cannot rule out the possibility that the feature is a PH2 radical centered on the dimer row, as suggested by Lin et al. Lin et al. (1999) Figure 1(d) shows a PH3 molecule after it has dissociated into PH2 \+ H on a single Si dimer, as described. Schofield et al. (2003) The observation of both adsorbed PH3 in its molecular form, and as the dissociation products PH2 \+ H is in agreement with FTIR Shan et al. (1996) and photoemission Lin et al. (1999) results.

The third type of features in Figure 1(b) are single and multiple dimer vacancies, resulting directly from the removal of Si from the surface layer during PH3 dosing. While dimer vacancies are difficult to distinguish from monohydride dimers (formed when each Si atom of a dimer is terminated with an H atom) the formation of monohydride dimers is not expected at room temperature since single H atoms have insufficient mobility to diffuse and pair-up into monohydrides. Bowler et al. (2000) We assign the fourth feature to be a Si ad-dimer, centered on the dimer rows with a height/diameter of 1.3 Å/13 Å in the filled state image and a height/diameter of 1.9 Å /13 Å in the empty state, in agreement with the observation of Si ad-dimers by Swartzentruber Swartzentruber (1997) and shown in Fig. 1(e) and (f) respectively. The filled state height of the feature (1.3 Å) is of the same order as the distance between  planes of the Si crystal (1.36 Å) as expected for a Si ad-dimer. We discount the suggestion of Wang et al. Wang et al. (1994b) that the feature is a Si monomer, since Si monomers are not observed in filled state images. Swartzentruber (1998) We also rule out the fourth feature being either adsorbed P or PH species because Fourier transform infrared spectroscopy Shan et al. (1996) and photoemission Lin et al. (1999) results show that these species are not formed after room temperature PH3 adsorption on Si(001).

### III.2 P identification using STM and AES

We now describe the incorporation of P atoms into the surface, to form Si-P heterodimers, upon annealing. Fig. 2(a) shows a Si(001) surface that has been exposed to 0.002 L of PH3 and annealed to 350C. Three distinct types of feature are present on the surface, Si-P heterodimers, ejected Si dimer chains and monohydride dimers. The ejected Si chains appear as short bright lines running perpendicular to the dimer rows. These features have the same height above the surface (1.36 Å) as a monolayer Si step Da̧browski and Müssig (2000) and have a characteristic minima running lengthways along their center in empty state images (see Fig. 4(d)), indicating that they are indeed Si dimer chains as suggested by Wang et al. Wang et al. (1994b) and not P dimer chains as proposed by Kipp et al. Kipp et al. (1995) The monohydride dimers (H-Si-Si-H) appear as dark features, similar in appearance to missing dimer defects.

Fig. 2(b) shows a high resolution image of two Si-P heterodimers, which appear as zig-zags along the dimer row, with the highest part $``$ 0.3 Å above the surface dimers. Oberbeck et al. (2002) This zig-zag appearance indicates there is static buckling of the surface dimers on either side of the Si-P heterodimer. Figure 2(b) shows that this buckling decays in amplitude as a function of distance from the Si-P heterodimer and extends 2-3 dimers away from it.

To confirm the interpretation of Fig. 2(a) presented above, we compare the number of Si-P heterodimers in STM images with the fraction of P at the surface measured by Auger electron spectroscopy (AES). In order to control the number of Si-P heterodimers at the surface without the complication of ejected Si dimer chains and adsorbed H (Fig. 2(a)), we have created a buried P $`\delta `$-doped layer and then diffused varying amounts of the P atoms to the surface using a range of anneal temperatures. The $`\delta `$-layer was produced by saturation dosing a Si(001) surface with PH3, annealing to 600C, and encapsulating with 5 monolayers of Si at room temperature. Oberbeck et al. (2003) Upon the annealing of this sample the buried P atoms diffuse to the surface at a rate determined by the anneal temperature. Figures 2(c) and 2(d) show STM images obtained after annealing the P $`\delta `$-doped sample for 60 seconds at 450C and 600C respectively. Figure 2(e) shows the 120 eV P AES peak obtained from the two surfaces. There is a direct correlation between the number of zig-zag features in the two STM images and the intensity of the 120 eV P Auger peak, confirming that the zig-zag feature contain P atoms. From the zig-zag appearance we conclude that each of these features contain only one P atom since P-P dimers are known to have have mirror symmetry about a line running along their center Wang et al. (1994a) and a feature consisting of neighboring incorporated P atoms along would show neighboring atoms with similar intensities.

While the results presented here are the first to characterize the isolated Si-P heterodimer, the observation that P atoms incorporate into the Si(001) surface upon annealing fits well with the Wang et al. Wang et al. (1994b) interpretation of the annealing of the low dosed PH3/Si(001), and contradicts the interpretation of Kipp et al. Kipp et al. (1995) We therefore conclude, following Wang et al., that the processes that occur as a result of the 350C anneal are as follows; the PHx (x=2,3) molecules firstly undergo complete dissociation. Each P atom then undergoes a substitution reaction with one atom of the Si surface dimer to form a Si-P heterodimer in the plane of the surface, thereby ejecting the displaced Si atom onto the surface. At 350C the displaced Si atoms are mobile enough to diffuse short distances Mo et al. (1991) and subsequently form 1-D dimer chains on the surface. The H atoms liberated from PHx dissociation are highly mobile at 350C and pair-up to form the energetically favorable monohydride dimers. Shan et al. (1996)

### III.3 The changing appearance of the Si-P heterodimer with coverage

We have found that the appearance of the Si-P heterodimer varies greatly as a function of initial PH3 coverage. Figure 3 shows a Si(001) surface after exposure to 2 different doses (0.002 L and 0.01 L) of PH3 and subsequent annealing of the dosed surfaces to 350C and then 500C. If we compare the two 350C anneals, Fig. 3(a) and (b), we see that the surface that was exposed to the higher PH3 dose has significantly more Si chains and the onset of 2D Si island formation (Fig. 3(b)), as expected. After both the lower and higher dosed surfaces have been subject to a further anneal of 500C the Si chains have now disappeared, see Fig. 3(c) and (d). This can be explained since at 500C the diffusion of Si atoms on the surface is sufficient that all the Si ad-dimers migrate to step edges to which they preferentially bond.

We now consider what information we can determine about the Si-P heterodimers, from Fig. 3. If we just consider Fig. 3(c) and 3(d), which show the lower and higher dosed surfaces after annealing to 500C, the data fits our explanation of the evolution of the PH3 dosed Si(001) surface i.e. the more PHx (x=2,3) initially adsorbed on the surface, the more Si-P heterodimers are observed on the surface after the 500C anneal. However, the 350C anneals of the lower and higher dosed surfaces are not so easy to understand. We repeatedly observe that lower dosed, 350C annealed surfaces such as that of Fig. 3(a) show moderate densities of Si-P heterodimers and monohydrides, as expected. Conversely, the higher dosed, 350C annealed surfaces such as that of Fig. 3(b) appear to show a high density of monohydrides and no Si-P heterodimers. We know from the 500C anneal of the higher dosed surfaces that Si-P heterodimers are abundant, so we can conclude that a high density of monohydrides inhibits the observation of Si-P heterodimers with the STM. It is known that monohydride dimers exhibit very little buckling Smith (2003) and thus pin neighboring Si dimers in an unbuckled geometry. Hence, if the density of monohydrides is sufficiently high it is reasonable to assume that almost all the surface Si dimers are pinned, even those in the vicinity of Si-P heterodimers. As a consequence, the characteristic zig-zag appearance of the Si-P heterodimer will be suppressed. It is now interesting to look again at Fig. 2(d) which was obtained by annealing a buried P layer (with no H present) to 600C. Whilst a sufficient number of P atoms have diffused to the surface to yield a high density of Si-P heterodimers, the surface is free of H. In this image the Si-P heterodimers are clearly visible, which is consistent with the observation that it is a high density of monohydrides on the surface that inhibits the STM observation of Si-P heterodimers.

Our final anneal of the PH3/Si(001) low dosed surface (not shown) was performed at 700C and resulted in the disappearance of the Si-P heterodimers. By comparing our results with that from the annealing of a saturation dosed PH3/Si(001) surface, Lin et al. (1999) we conclude that the P atoms have desorbed from the surface, leaving only features inherent to the clean Si(001) surface.

### III.4 Identification of the Si-P heterodimer using empty state imaging

In the previous section we showed that for high coverages of monohydrides, the identification of Si-P heterodimers in filled state images was not possible. However, we will show that we can clearly identify Si-P heterodimers, at any coverage, using empty state imaging. Figure 4 shows (a) filled- and (b) empty-state images of the same area of a Si(001) sample after a 0.002 L PH3 dose and a 500C anneal. We have chosen a low PH3 dose and an anneal temperature sufficiently high that ejected Si dimer chains and monohydrides have been completely removed from the surface, so that we can image Si-P heterodimers in isolation from other features. Figure 4(b) shows that in empty state images the Si-P heterodimer appears as a bright circular protrusion on one side of a dimer row with a slight depression on the other side of the row. The height of the protrusion, measured with respect to the surface dimers, varies linearly with bias between $`0.50\pm 0.05`$ Å at +0.8 V to $`0.10\pm 0.02`$ Å at +2.0 V. Dimer buckling in empty state images of Si(001) is not observed for sample biases greater than +0.8 V. Qin and Lagally (1999)

With the Si-P heterodimer well characterized in both filled and empty state imaging we can now also interpret images obtained after higher PH3 doses but a lower anneal temperature, where H is still present. Figures 4(c) and (d) are dual bias images of a Si(001) surface that has been dosed with 0.013 L of PH3 and annealed to 350C. The filled state image shows a number of Si chains and monohydrides but no Si-P heterodimers are observed. We see that the ejected Si chains have a characteristic minima running lengthways along their center in the empty state image, as expected. Wang et al. (1994b); Swartzentruber (1997) However, despite the difficulty in resolving features due to the extreme brightness of the ejected Si, in the empty state image we still see circular protrusions characteristic of the Si-P heterodimers, as well as the Si chains and the monohydrides. This ability to identify incorporated P atoms regardless of the density of other features on the surface is a useful tool for locating P atoms in Si, towards the goal of realizing atomic-scale devices. Tucker and Shen (2000); Schofield et al. (2003)

## IV Conclusion

We have used scanning tunneling microscopy (STM) and Auger electron spectroscopy (AES) to study the incorporation of P into Si(001) by dosing with small quantities of PH3 and annealing. The Si-P heterodimer appears as a zig-zag feature when imaged in filled state STM images due to the static dimer buckling induced by the heterodimer. The Si-P assignment was confirmed from the direct correlation between the number of zig-zag features imaged by STM and the fraction of P at the surface measured by AES. The presence of moderate numbers of monohydrides on the surface removes the buckling induced by the Si-P heterodimer and renders the heterodimer invisible in filled state images. We characterized the Si-P heterodimer using empty state imaging and show that, in the empty state the Si-P heterodimers are visible even with a high coverage of monohydrides. The demonstration of the incorporation of individual P atoms into the Si(001) surface, and subsequent characterization using STM, is an important step towards fabrication of nanoscale electronic devices based on controlled atomic-scale doping of Si.

###### Acknowledgements.

This work was supported by the Australian Research Council, the Australian Government and by the US National Security Agency (NSA), Advanced Research and Development Activity (ARDA) and the Army Research Office (ARO) under contract number DAAD19-01-1-0653. LO acknowledges a Hewlett-Packard Fellowship.