text:
# On the Numerical Stability of the Einstein Equations

## I Introduction

Numerically solving the full 3D nonlinear Einstein equations is, for several reasons, a daunting task. Still, numerical relativity remains the best method for studying astrophysically interesting regions of the solution space of the Einstein equations in sufficient detail and accuracy in order to be used to interpret measurements made by the up and coming gravitational wave detectors. Even though numerical relativity is almost 35 years old, some of the same problems faced by researchers three decades ago are present today. Aside from the computational complexity of implementing a numerical solver for the nonlinear Einstein equations, there exist several unsolved problems, including the well-posedness of certain initial value formulations of the Einstein equations and the proper choice of gauge. Not the least of these problems is numerical stability. A common thread in numerical relativity research over the past three decades is the observation of high frequency (Nyquist frequency) noise growing and dominating the numerical solution. Traditionally, numerical studies have been performed with the initial value formulation of the Einstein equations known as the ADM 3+1 formulation, in which the 3-metric and extrinsic curvature are the dynamically evolved variables.

Lately, a formulation based on variables in which the conformal factor of the 3-metric and the trace of the extrinsic curvature are factored out and evolved separately has been studied. This conformal-traceless (CT) formulation was first introduced by Nakamura and Shibata  and later slightly modified by Baumgarte and Shapiro . The stability properties of the CT formulation were shown in  to be better than those of the ADM formulation for linear waves. The improvement in numerical stability in the CT formulation versus the ADM formulation was demonstrated in strong field dynamical cases in . A step toward understanding the improved stability properties of the CT formulation was taken in  where it was shown by analytically linearizing the ADM and CT equations about flat space that the CT system effectively decouples the gauge modes and constraint violating modes. It was conjectured that giving the constraint violating modes nonzero propagation speed results in a stable evolution.

Here, we take another step towards understanding the improved stability properties of the CT system by performing a von Neumann stability analysis on discretizations of both the ADM and CT systems. We are led to the von Neumann stability analysis by Lax’s equivalence theorem , which states that given a well posed initial value problem and a discretization that is consistent with that initial value problem (i.e., the finite difference equations are faithful to the differential equations), then stability is equivalent to convergence. Here, the words “stability” and “convergence” are taken to mean very specific things. Convergence is taken to mean pointwise convergence of solutions of the finite difference equations to solutions of the differential equations. This is the pièce de résistance of numerical relativity. After all, what we are interested in are solutions to the differential equations. Stability, on the other hand, has a rather technical definition involving the uniform boundedness of the discrete Fourier transform of the finite difference update operator (see  for details). In essence, stability is the statement that there should be a limit to the extent to which any component of an initial discrete function can be amplified during the numerical evolution procedure (note that stability is a statement concerning the finite difference equations, not the differential equations). Fortunately, the technical definition of stability can be shown to be equivalent to the von Neumann stability condition, which will be described in detail in the next section.

While one cannot apply Lax’s equivalence theorem directly in numerical relativity (the initial value problem well-posedness assumption is not valid for the Einstein field equations in that the evolution operator is not, in general, uniformly bounded ), numerical relativists often use it as a “road map”; clearly consistency and stability are important parts of any discretization of the Einstein equations (curiously, convergence is usually implicitly assumed in most numerical studies). Code tests, if done at all, usually center around verifying the consistency of the finite difference equations to the differential equations (as an example of the extents to which some numerical relativists will go to check the consistency of the finite difference equations to the differential equations, see, e.g., ). Stability, on the other hand, is usually assessed postmortem. If the code crashes immediately after a sharp rise in Nyquist frequency noise and/or if the code crashes sooner in coordinate time at higher resolutions, the code is deemed unstable.

We suggest that the stability of the code can be assessed before (and perhaps even more importantly, while) numerical evolutions take place. As will be seen in the next section, the stability properties of any given nonlinear finite difference update operator depend not only on the Courant factor $`\frac{\mathrm{\Delta }t}{\mathrm{\Delta }z}`$, but also on the values of the discrete evolution variables themselves. Therefore, during numerical evolutions of nonlinear problems, as the evolved variables change from discrete timestep to timestep, the stability properties of the finite difference operator change along with them! Ideally, one would want to verify that the finite difference update operator remains stable for each point in the computational domain at each timestep. While the computational expense of this verification would be prohibitive, verification at a reasonably sampled subset of discrete points could be feasible.

The remainder of the paper is outlined as follows. Section II will describe the von Neumann stability analysis for a discretization of a general set of nonlinear partial differential equations in one spatial dimension. Included will be results from a von Neumann stability analysis of the linear wave equation discretized with an iterative Crank-Nicholson scheme. Section III will present the ADM and CT formulations of the Einstein equations restricted to diagonal metrics and further restricted to dependence on only one spatial variable. The von Neumann stability analysis is performed on iterative Crank-Nicholson discretizations of both formulations with flat initial data. Section IV will repeat the von Neumann stability analysis from section III with a nonlinear wave solution for initial data.

## II The von Neumann Stability Analysis

Let a set of partial differential equations be given as

$$\frac{\stackrel{}{u}}{t}=S(\stackrel{}{u},_z\stackrel{}{u},_{z}^{}{}_{}{}^{2}\stackrel{}{u})$$
(1)
where $`\stackrel{}{u}=\stackrel{}{u}(z,t)`$ is a vector whose components consist of the dependent variables that are functions of the independent variables $`0zN`$ and $`t0`$ (while we ignore boundary conditions in this paper, as we consider only interior points of equations with finite propagation speeds, the von Neumann analysis presented here is easily extended to include boundary treatments ).

Now, consider a discretization of the independent variables $`z`$ and $`t`$:

$`z_j=j\mathrm{\Delta }z,`$ $`j=0,1,2,\mathrm{},m`$ (2)
$`t_n=n\mathrm{\Delta }t,`$ $`n=0,1,2,\mathrm{}`$ (3)
along with a discretization of the dependent variables $`\stackrel{}{u}(z,t)`$:

$$\stackrel{}{u}_j^n=\stackrel{}{u}(z_j,t_n).$$
(4)
Furthermore, consider a consistent discretization of the set of differential equations in the form of Eq. 1 that can be written in the form

$$\stackrel{}{u}_j^{n+1}=\stackrel{}{S}_j^n(\mathrm{},\stackrel{}{u}_{j2}^n,\stackrel{}{u}_{j1}^n,\stackrel{}{u}_j^n,\stackrel{}{u}_{j+1}^n,\stackrel{}{u}_{j+2}^n,\mathrm{}).$$
(5)
As we will only be analyzing the iterative Crank-Nicholson discretization scheme, we assume a two-step method as written in Eq. 5. However, a von Neumann analysis does not depend on this, and could easily be performed for three-step methods such as the leapfrog scheme.

Now, assume some initial data

$$\stackrel{}{u}_j^0,\mathrm{\hspace{0.33em}\hspace{0.33em}0}jm$$
(6)
for the discrete variables is given at time $`t=0`$. The amplification matrix at $`z_j`$ is given by

$$𝐆_j^n=\underset{j^{}=0}{\overset{m}{}}\frac{\stackrel{}{S}_j^n}{\stackrel{}{u}_j^{}^n}e^{i\mathrm{\Delta }zk_z(jj^{})}$$
(7)
The condition for numerical stability  is that the spectral radius of the amplification matrix $`𝐆_j^n`$ be 1 or less for each wavenumber $`k_z`$. That is, each eigenvalue $`\lambda _i`$ of the amplification matrix $`𝐆_j^n`$ should have a modulus of 1 or less for each discrete mode $`k_z`$

$$\left|\lambda _i\right|1.$$
(8)
First, notice that for linear finite difference equations (Eq. 5), the amplification matrix $`𝐆_j^n`$ does not depend on either the initial data $`\stackrel{}{u}_j^0`$ nor on the spatial index $`j`$. It only depends on the discretization parameters $`\mathrm{\Delta }z`$ and $`\mathrm{\Delta }t`$ as well as the mode wavenumber $`k_z`$. That is, for the linear case, once one specifies the discretization parameters $`\mathrm{\Delta }z`$ and $`\mathrm{\Delta }t`$, one only needs to verify that the eigenvalues of the amplification matrix $`𝐆_j^n`$ have a modulus of less than or equal to 1 for all wavenumbers $`k_z`$ to insure the numerical stability of a numerical update operator, regardless of initial data $`\stackrel{}{u}_j^0`$. Contrast this with the nonlinear case where the amplification matrix $`𝐆_j^n`$ is not only a function of discretization parameters $`\mathrm{\Delta }z`$, $`\mathrm{\Delta }t`$ and the mode wavenumber $`k_z`$, but also of initial data $`\stackrel{}{u}_j^0`$. In this case, not only must one verify that the amplification matrix $`𝐆_j^n`$ have a spectral radius of 1 or less for each spatial index $`j`$ (assuming the initial data $`\stackrel{}{u}_j^0`$ depends on $`j`$, which it will in general), but one must carry out this verification at every time step $`n`$, as the data $`\stackrel{}{u}_j^n`$ will, in general, change with increasing $`n`$. So, in the nonlinear case, the amplification matrix $`𝐆_j^n`$ will depend on both the spatial index $`j`$ and the temporal index $`n`$. In principle, one must verify that $`𝐆_j^n`$ have a spectral radius of 1 or less for each mode wavenumber $`k_z`$, for each spatial index $`j`$ and for each time step $`n`$, in order to be confident that the solutions of the finite difference equations are converging to solutions of the differential equations (which is, after all, what we are ultimately interested in).

### A Example: Linear Wave Equation

In , a von Neumann analysis of the advection equation was presented for the iterative Crank-Nicholson scheme. Here, as a prelude to computing a von Neumann analysis for discretizations of the equations of general relativity, we present a von Neumann analysis for a 2-iteration iterative Crank-Nicholson discretization of the wave equation in 1 dimension

$$\frac{^2\varphi }{t^2}=\frac{^2\varphi }{z^2}.$$
(9)
We write the wave equation in first-order form as in Eq. 1. Defining $`D\frac{\varphi }{z}`$, $`\mathrm{\Pi }\frac{\varphi }{t}`$, and

$$\stackrel{}{u}\left[\begin{array}{c}\varphi \\ D\\ \mathrm{\Pi }\end{array}\right],$$
(10)
the wave equation in one spatial dimension becomes

$$\frac{\stackrel{}{u}}{t}=\left[\begin{array}{c}\mathrm{\Pi }\\ \frac{\mathrm{\Pi }}{z}\\ \frac{D}{z}\end{array}\right].$$
(11)
The iterative Crank-Nicholson discretization procedure begins by taking a FTCS (Forward Time Centered Space) step which is used to define the intermediate variables $`{}_{}{}^{(1)}\stackrel{~}{\stackrel{}{u}}_{j}^{n+1}`$:

$$\frac{{}_{}{}^{(1)}\stackrel{~}{\stackrel{}{u}}_{j}^{n+1}\stackrel{}{u}_j^n}{\mathrm{\Delta }t}=\left[\begin{array}{c}\mathrm{\Pi }_j^n\\ (\mathrm{\Pi }_{j+1}^n\mathrm{\Pi }_{j1}^n)/(2\mathrm{\Delta }z)\\ (D_{j+1}^nD_{j1}^n)/(2\mathrm{\Delta }z)\end{array}\right].$$
(12)
This intermediate state variable is averaged with the original state variable

$${}_{}{}^{(1)}\stackrel{}{u}_{j}^{n+1/2}=\frac{1}{2}({}_{}{}^{(1)}\stackrel{~}{\stackrel{}{u}}_{j}^{n+1}+\stackrel{}{u}_j^n)$$
(13)
which, in turn, is used to calculate $`{}_{}{}^{(1)}\stackrel{}{u}_{j}^{n+1}`$, the state vector produced from the first iteration of the iterative Crank-Nicholson scheme:

$$\frac{{}_{}{}^{(1)}\stackrel{}{u}_{j}^{n+1}\stackrel{}{u}_j^n}{\mathrm{\Delta }t}=\left[\begin{array}{c}{}_{}{}^{(1)}\mathrm{\Pi }_{j}^{n+1/2}\\ ({}_{}{}^{(1)}\mathrm{\Pi }_{j+1}^{n+1/2}{}_{}{}^{(1)}\mathrm{\Pi }_{j1}^{n+1/2})/(2\mathrm{\Delta }z)\\ ({}_{}{}^{(1)}D_{j+1}^{n+1/2}{}_{}{}^{(1)}D_{j1}^{n+1/2})/(2\mathrm{\Delta }z)\end{array}\right].$$
(14)
The averaged state $`{}_{}{}^{(2)}\stackrel{}{u}_{j}^{n+1/2}`$ is calculated

$${}_{}{}^{(2)}\stackrel{}{u}_{j}^{n+1/2}=\frac{1}{2}({}_{}{}^{(1)}\stackrel{}{u}_{j}^{n+1}+\stackrel{}{u}_j^n)$$
(15)
and used to compute the final state variable $`{}_{}{}^{(2)}\stackrel{}{u}_{j}^{n+1}`$

$$\frac{{}_{}{}^{(2)}\stackrel{}{u}_{j}^{n+1}\stackrel{}{u}_j^n}{\mathrm{\Delta }t}=\left[\begin{array}{c}{}_{}{}^{(2)}\mathrm{\Pi }_{j}^{n+1/2}\\ ({}_{}{}^{(2)}\mathrm{\Pi }_{j+1}^{n+1/2}{}_{}{}^{(2)}\mathrm{\Pi }_{j1}^{n+1/2})/(2\mathrm{\Delta }z)\\ ({}_{}{}^{(2)}D_{j+1}^{n+1/2}{}_{}{}^{(2)}D_{j1}^{n+1/2})/(2\mathrm{\Delta }z)\end{array}\right].$$
(16)
Using this second iteration as the final iteration, we can write Eq. 5 explicitly as the three equations

$`\varphi _j^{n+1}`$ $`=`$ $`\varphi _j^n+{\displaystyle \frac{\sigma \mathrm{\Delta }t}{4}}(D_{j+1}^nD_{j1}^n)+`$ (18)
$`{\displaystyle \frac{\sigma ^2\mathrm{\Delta }t}{16}}(\mathrm{\Pi }_{j+2}^n\mathrm{\Pi }_{j2}^n)+\mathrm{\Delta }t(1{\displaystyle \frac{\sigma ^2}{16}})\mathrm{\Pi }_j^n`$
$`D_j^{n+1}`$ $`=`$ $`(1{\displaystyle \frac{\sigma ^2}{4}})D_j^n+{\displaystyle \frac{\sigma ^2}{8}}(D_{j+2}^n+D_{j2}^n)+`$ (20)
$`{\displaystyle \frac{\sigma ^3}{32}}(\mathrm{\Pi }_{j+3}^n\mathrm{\Pi }_{j3}^n)+{\displaystyle \frac{\sigma }{2}}(1{\displaystyle \frac{3}{16}}\sigma ^2)(\mathrm{\Pi }_{j+1}^n\mathrm{\Pi }_{j1}^n)`$
$`\mathrm{\Pi }_j^{n+1}`$ $`=`$ $`(1{\displaystyle \frac{\sigma ^2}{4}})\mathrm{\Pi }_j^n+{\displaystyle \frac{\sigma ^2}{8}}(\mathrm{\Pi }_{j+2}^n+\mathrm{\Pi }_{j2}^n)+`$ (22)
$`{\displaystyle \frac{\sigma ^3}{32}}(D_{j+3}^nD_{j3}^n)+{\displaystyle \frac{\sigma }{2}}(1{\displaystyle \frac{3}{16}}\sigma ^2)(D_{j+1}^nD_{j1}^n),`$
where we have denoted $`\sigma `$ as the Courant factor $`\mathrm{\Delta }t/\mathrm{\Delta }z`$. We compute the amplification matrix Eq. 7, which is independent of the spatial index $`j`$ and the time index $`n`$:

$$𝐆=\left[\begin{array}{ccc}G_{11}& G_{12}& G_{13}\\ G_{21}& G_{22}& G_{23}\\ G_{31}& G_{32}& G_{33}\end{array}\right]$$
(23)
where the components are

$$\begin{array}{c}G_{11}=1\hfill \\ G_{21}=G_{31}=0\hfill \\ G_{12}=\frac{i\sigma \mathrm{\Delta }t}{2}\mathrm{sin}\theta _k\hfill \\ G_{13}=\mathrm{\Delta }t(1\frac{\sigma ^2}{8})+\frac{\mathrm{\Delta }t\sigma ^2}{8}\mathrm{cos}(2\theta _k)\hfill \\ G_{22}=G_{33}=(1\frac{\sigma ^2}{4})+\frac{\sigma ^2}{4}\mathrm{cos}(2\theta _k)\hfill \\ G_{23}=G_{32}=\frac{i\sigma ^3}{16}\mathrm{sin}(3\theta _k)i\sigma (1\frac{3}{16}\sigma ^2)\mathrm{sin}\theta _k,\hfill \end{array}$$
(24)
where $`\theta _kk_z\mathrm{\Delta }z`$. The modulus of the eigenvalues of $`𝐆`$ are easily calculated and found to be (ignoring the eigenvalue that is identically 1)

$$|\lambda _1|=|\lambda _2|=\sqrt{1\frac{1}{4}\sigma ^4\mathrm{sin}^4\theta _k+\frac{1}{16}\sigma ^6\mathrm{sin}^6\theta _k}$$
(25)
By inspection, the von Neumann condition for stability, i.e. $`|\lambda |1`$ for all $`k_z`$, is $`\sigma \frac{\mathrm{\Delta }t}{\mathrm{\Delta }z}2`$. However, it is instructive to look at the $`\theta _k`$ dependence of the stability criterion. Figure 1 shows a plot of Eq. 25 as a function of $`\theta _kk_z\mathrm{\Delta }z`$ for various values of the Courant factor $`\sigma `$. As can be seen, the first eigenmode that goes unstable (i.e., has an eigenvalue whose modulus is greater than 1) as $`\sigma `$ is increased is the mode $`\theta _kk_z\mathrm{\Delta }z=\pi `$. This corresponds to modes of wavelength $`2\mathrm{\Delta }z`$, i.e. Nyquist frequency modes, which are precisely the modes that usually crop up and kill numerical relativity simulations.

## III General Relativity

Here, we present the analytic equations for the ADM and CT formulations that we will use to study the stability properties of numerical methods in general relativity. We will discretize the general relativity evolution equations using the same discretization used in the previous section for discretizing the scalar wave equation, namely, a 2-iteration iterative Crank-Nicholson scheme. Before analyzing the stability of the discretization for non-linear waves, we will perform a von Neumann analysis for the discretizations of both the ADM and CT equations about flat space.

### A ADM Equations in 1-D (plane symmetry)

The form of the metric that we will use to study stability properties of discretizations of the ADM form of the Einstein equations is given by

$$ds^2=\alpha ^2dt^2+g_{xx}dx^2+g_{yy}dy^2+g_{zz}dz^2,$$
(26)
where the lapse $`\alpha `$ and metric functions $`g_{xx}`$, $`g_{yy}`$, and $`g_{zz}`$ are functions of the independent variables $`z`$ and $`t`$. To put the evolution equations in first order form (Eq. 1), we introduce the extrinsic curvature functions $`K_{xx}`$, $`K_{yy}`$, and $`K_{zz}`$, each of which are also function of the independent variables $`z`$ and $`t`$. The evolution equations for this ADM system is given in the form of Eq. 1 as

$`{\displaystyle \frac{g_{ii}}{t}}`$ $`=`$ $`2\alpha K_{ii}`$ (27)
$`{\displaystyle \frac{K_{ii}}{t}}`$ $`=`$ $`\alpha R_{ii}2\alpha K_{i}^{}{}_{}{}^{j}K_{ij}+\alpha KK_{ii}𝒟_i𝒟_i\alpha `$ (28)
$`{\displaystyle \frac{\alpha }{t}}`$ $`=`$ $`\alpha K`$ (29)
where the index pair $`ii`$ takes on the values $`\{xx,yy,zz\}`$, the index $`j`$ is summed over the values $`\{x,y,z\}`$, $`K`$ denotes the trace of the extrinsic curvature $`K=K_{xx}/g_{xx}+K_{yy}/g_{yy}+K_{zz}/g_{zz}`$, $`𝒟`$ denotes the covariant derivative operator compatible with the 3-metric, and $`R_{ii}`$ denotes the three components of the 3-Ricci tensor, which are given explicitly as

$`R_{xx}`$ $`=`$ $`{\displaystyle \frac{1}{4g_{xx}g_{zz}}}\left({\displaystyle \frac{g_{xx}}{z}}\right)^2{\displaystyle \frac{1}{4g_{yy}g_{zz}}}{\displaystyle \frac{g_{xx}}{z}}{\displaystyle \frac{g_{yy}}{z}}+`$ (31)
$`{\displaystyle \frac{1}{4g_{zz}^{}{}_{}{}^{2}}}{\displaystyle \frac{g_{xx}}{z}}{\displaystyle \frac{g_{zz}}{z}}{\displaystyle \frac{1}{2g_{zz}}}{\displaystyle \frac{^2g_{xx}}{z^2}}`$
$`R_{yy}`$ $`=`$ $`{\displaystyle \frac{1}{4g_{yy}g_{zz}}}\left({\displaystyle \frac{g_{yy}}{z}}\right)^2{\displaystyle \frac{1}{4g_{xx}g_{zz}}}{\displaystyle \frac{g_{xx}}{z}}{\displaystyle \frac{g_{yy}}{z}}+`$ (33)
$`{\displaystyle \frac{1}{4g_{zz}^{}{}_{}{}^{2}}}{\displaystyle \frac{g_{yy}}{z}}{\displaystyle \frac{g_{zz}}{z}}{\displaystyle \frac{1}{2g_{zz}}}{\displaystyle \frac{^2g_{yy}}{z^2}}`$
$`R_{zz}`$ $`=`$ $`{\displaystyle \frac{1}{4g_{xx}^{}{}_{}{}^{2}}}\left({\displaystyle \frac{g_{xx}}{z}}\right)^2+{\displaystyle \frac{1}{4g_{yy}^{}{}_{}{}^{2}}}\left({\displaystyle \frac{g_{yy}}{z}}\right)^2+`$ (36)
$`{\displaystyle \frac{1}{4g_{xx}g_{zz}}}{\displaystyle \frac{g_{xx}}{z}}{\displaystyle \frac{g_{zz}}{z}}+{\displaystyle \frac{1}{4g_{yy}g_{zz}}}{\displaystyle \frac{g_{yy}}{z}}{\displaystyle \frac{g_{zz}}{z}}`$
$`{\displaystyle \frac{1}{2g_{xx}}}{\displaystyle \frac{^2g_{xx}}{z^2}}{\displaystyle \frac{1}{2g_{yy}}}{\displaystyle \frac{^2g_{yy}}{z^2}}`$
Throughout this paper, we use the so-called “1 + log” slicing condition in Eq. 27. This local condition on the lapse has been used successfully in several recent applications . Moreover, it is a local condition, and thus, the von Neumann analysis remains local (this would be in contrast with a global elliptic condition on the lapse, e.g. maximal slicing, in which the calculation of the sum in the definition of the amplification matrix in Eq. 7 would have nonzero global contributions).

The initial value formulation is completed by specifying initial data that satisfies the Hamiltonian and momentum constraints, given respectively by

$``$ $``$ $`{}_{}{}^{3}R+K^2K_{jk}K^{jk}=0`$ (37)
$`_z`$ $``$ $`𝒟_jK_{}^{j}{}_{z}{}^{}𝒟_zK`$ (38)
### B Conformal-Traceless (CT) Equations in 1-D (plane symmetry)

The form of the metric that we will use to study stability properties of the discretizations of the CT form of the Einstein equations, as defined in  is given by

$$ds^2=\alpha ^2dt^2+e^{4\varphi }(\stackrel{~}{g}_{xx}dx^2+\stackrel{~}{g}_{yy}dy^2+\stackrel{~}{g}_{zz}dz^2),$$
(39)
where the lapse $`\alpha `$, the conformal function $`\varphi `$, and the conformal metric components $`\stackrel{~}{g}_{xx}`$, $`\stackrel{~}{g}_{yy}`$, and $`\stackrel{~}{g}_{zz}`$ are functions of the independent variables $`z`$ and $`t`$. The determinant of the conformal 3-metric is identically 1. Instead of evolving the extrinsic curvature components as in the ADM formalism, the extrinsic curvature is split into its trace ($`K`$) and traceless ($`\stackrel{~}{A}_{ii}`$) components:

$$K_{ii}=e^{4\varphi }(\frac{1}{3}\stackrel{~}{g}_{ii}K+\stackrel{~}{A}_{ii}).$$
(40)
In addition, the conformal connection function $`\stackrel{~}{\mathrm{\Gamma }}^z`$, defined by

$$\stackrel{~}{\mathrm{\Gamma }}^z=_j\stackrel{~}{g}^{jz}$$
(41)
is also treated as an evolved variable. There are therefore 10 evolution equations that are in the form of Eq. 1, and given explicitly as

$`{\displaystyle \frac{\varphi }{t}}`$ $`=`$ $`{\displaystyle \frac{1}{6}}\alpha K`$ (42)
$`{\displaystyle \frac{\stackrel{~}{g}_{ii}}{t}}`$ $`=`$ $`2\alpha \stackrel{~}{A}_{ii}`$ (43)
$`{\displaystyle \frac{K}{t}}`$ $`=`$ $`\alpha \stackrel{~}{A}_{jk}\stackrel{~}{A}^{jk}+{\displaystyle \frac{1}{3}}\alpha K^2𝒟_j𝒟^j\alpha `$ (44)
$`{\displaystyle \frac{\stackrel{~}{A}_{ii}}{t}}`$ $`=`$ $`\alpha (e^{4\varphi }R_{ii}{\displaystyle \frac{1}{3}}\stackrel{~}{g}_{ii}\stackrel{~}{A}^{jk}\stackrel{~}{A}^{jk}+{\displaystyle \frac{2}{9}}\stackrel{~}{g}_{ii}K^2+K\stackrel{~}{A}_{ij}2\stackrel{~}{A}_{ij}\stackrel{~}{A}_{i}^{}{}_{}{}^{j})`$ (46)
$`e^{4\varphi }(𝒟_i𝒟_i\alpha {\displaystyle \frac{1}{3}}g_{ii}𝒟_j𝒟^j\alpha )`$
$`{\displaystyle \frac{\stackrel{~}{\mathrm{\Gamma }}^z}{t}}`$ $`=`$ $`\alpha ({\displaystyle \frac{4}{3\stackrel{~}{g}_{zz}}}{\displaystyle \frac{K}{z}}+{\displaystyle \frac{12}{\stackrel{~}{g}_{zz}^{}{}_{}{}^{2}}}\stackrel{~}{A}_{zz}{\displaystyle \frac{\varphi }{z}}+2\stackrel{~}{\mathrm{\Gamma }}_{}^{z}{}_{ij}{}^{}\stackrel{~}{A}^{ij})`$ (47)
$`{\displaystyle \frac{\alpha }{t}}`$ $`=`$ $`\alpha K`$ (48)
where the indices of $`\stackrel{~}{A}_{ij}`$ are raised and lowered with the conformal metric $`\stackrel{~}{g}_{ij}`$, the index on the covariant derivative operator $`𝒟`$ with respect to the physical metric is raised and lowered with the physical metric $`g_{ij}=e^{4\varphi }\stackrel{~}{g}_{ij}`$, and $`\stackrel{~}{\mathrm{\Gamma }}_{}^{k}{}_{ij}{}^{}`$ are the Christoffel symbols related to the conformal metric $`\stackrel{~}{g}_{ij}`$. Note that the Hamiltonian constraint has been substituted in for the 3-Ricci scalar in the equations for $`\frac{K}{t}`$ and $`\frac{\stackrel{~}{A}_{ii}}{t}`$. Also, the momentum constraint has been substituted in the equation for $`\frac{\stackrel{~}{\mathrm{\Gamma }}^z}{t}`$. This corresponds to the “Mom” system from , and the $`\sigma =1`$, $`m=1`$, $`\xi =0`$ system from . Note that the Ricci components $`R_{ii}`$ can be written in terms of the conformal Ricci components $`\stackrel{~}{R}_{ii}`$:

$`R_{ii}`$ $`=`$ $`\stackrel{~}{R}_{ii}+2{\displaystyle \frac{\varphi }{z}}\stackrel{~}{\mathrm{\Gamma }}_{}^{z}{}_{ii}{}^{}2\stackrel{~}{g}_{ii}(\stackrel{~}{𝒟}_j\stackrel{~}{𝒟}^j\varphi +{\displaystyle \frac{2}{\stackrel{~}{g}_{zz}}}({\displaystyle \frac{\varphi }{z}})^2)+`$ (50)
$`\delta _{iz}(4({\displaystyle \frac{\varphi }{z}})^22{\displaystyle \frac{^2\varphi }{z^2}}),`$
where the indices of the covariant derivative operator $`\stackrel{~}{𝒟}`$ with respect to the conformal metric are raised and lowered with the conformal metric. The conformal Ricci components can in turn be written as

$`\stackrel{~}{R}_{ii}`$ $`=`$ $`{\displaystyle \frac{1}{2\stackrel{~}{g}_{zz}}}{\displaystyle \frac{^2\stackrel{~}{g}_{ii}}{z^2}}+\delta _{iz}\stackrel{~}{g}_{zz}{\displaystyle \frac{\stackrel{~}{\mathrm{\Gamma }}^z}{z}}`$ (52)
$`+{\displaystyle \frac{1}{2}}\stackrel{~}{\mathrm{\Gamma }}^z{\displaystyle \frac{\stackrel{~}{g}_{ii}}{z}}+\stackrel{~}{\mathrm{\Gamma }}_{}^{j}{}_{zk}{}^{}(2\stackrel{~}{\mathrm{\Gamma }}_z{}_{}{}^{k}{}_{j}{}^{}+\stackrel{~}{\mathrm{\Gamma }}_j{}_{}{}^{k}{}_{z}{}^{})`$
Notice that, although we have imposed planar symmetry, we have expressed the Ricci tensor in terms of derivatives of the conformal metric $`\stackrel{~}{g}_{ij}`$ and of the conformal connection function $`\stackrel{~}{\mathrm{\Gamma }}^z`$ just as is done for full 3-D numerical relativity .

### C von Neumann stability analysis about flat space

As outlined for the scalar wave equation in Section II A, we discretize both the ADM equations (Eq. 27) and the CT equations (Eq. 42) using a 2-iteration iterative Crank-Nicholson scheme. The resulting finite difference equations, even though planar symmetry and a simplified form of the metric is assumed, are still too complicated to perform a von Neumann analysis by hand. The complication arises due to the recursive nature of the iterative Crank-Nicholson method; a 2-iteration procedure results in the source terms of Eqs. 27 and 42 being computed 3 times in recursive succession.

We have performed the von Neumann analysis of the ADM and CT equations in two independent ways. In the first way, we used the symbolical calculation computer program MATHEMATICA to explicitly calculate the 2-iteration iterative Crank-Nicholson update, and then explicitly calculated the amplification matrix Eq. 7 in terms of the initial data variables. We then took these expressions, substituted the initial data about which we want to compute the von Neumann analysis, and calculated the eigenvalues to arbitrary precision inside MATHEMATICA. The second, independent way of performing the von Neumann analysis was to write an evolution code for the finite difference equations. We then input the initial data about which we want to compute the von Neumann analysis and computed the derivatives in the definition of the amplitude matrix Eq. 7 using finite differencing (this amounts to finite differencing the finite difference equations!). The package EISPACK was then used to compute the eigenvalues of the resulting amplification matrices. To obtain the highest accuracy possible, whenever calculating the amplification matrix Eq. 7 using finite differencing, we finite difference with multiple discretization parameters $`h`$, and use Richardson extrapolation to obtain values of the derivatives.

Both methods were used, and found to produce identical results. All results reported in the paper were produced using both methods as described above. We emphasize the use of two independent methods not only to verify our results, but also for more practical reasons: one may eventually want to perform a von Neumann analysis for a full 3-D numerical relativity code. The analytic method using a symbolical manipulation package such as MATHEMATICA may not be feasible in the near future. It took well over 100 hours on one node of an Origin 2000 running MATHEMATICA to perform the symbolical calculations needed for the von Neumann analysis of a plane symmetric code. It may be several orders of magnitude more expensive to analyze a full 3D evolution update. The finite difference method for computing the amplification matrix is much quicker, and it is reassuring that the finite difference method, used in conjunction with Richardson extrapolation, is accurate enough to reproduce the same detailed structures of the eigenvalues of the amplification matrices as doing the analytic calculation.

In Figure 2 we plot the maximum of the modulus of the eigenvalues (neglecting eigenvalues that are exactly 1) of the amplification matrix of the 2-iteration iterative Crank-Nicholson discretization scheme of the ADM equations from Section III A using flat space as initial data, namely, $`g_{ii}=\alpha =1`$, and $`K_{ii}=0`$. We see that the spectral radius of the amplification matrix is less than or equal to 1 for $`\sigma =\mathrm{\Delta }t/\mathrm{\Delta }x1`$. Notice that for the Nyquist frequency mode ($`\theta _kk_z\mathrm{\Delta }z=\pi `$), when the Courant factor $`\sigma `$ is $`1`$, all eigenvalues of the amplification matrix have a modulus of exactly 1. For $`\sigma >1`$, all Nyquist frequency modes are unstable (i.e. they all have amplification matrices with spectral radii $`>1`$).

In Figure 3 we plot the maximum of the modulus of the eigenvalues (neglecting eigenvalues that are exactly 1) of the amplification matrix of the 2-iteration iterative Crank-Nicholson discretization scheme of the CT equations from Section III B, again using flat space as initial data. The resulting plot is identical to that of the ADM equations for flat space. The stability properties of the CT equations about flat space are exactly the same as the stability properties of the ADM equations about flat space.

## IV Von Neumann analysis for nonlinear plane waves

In this section, we study the stability properties of the ADM and CT equations about nonlinear plane waves. We require initial data that corresponds to nonlinear plane waves that satisfy the constraints and takes on the form of our simplified metric Eqs. 26 and 39. We choose an exact plane wave solution first given by . The metric is assumed to take the form

$$ds^2=dt^2+L^2(e^{2\beta }dx^2+e^{2\beta }dy^2)+dz^2,$$
(53)
where $`L=L(v)`$, $`\beta =\beta (v)`$, and $`v=tz`$. Given an arbitrary function $`\beta (v)`$, the Einstein equations reduce to the following ordinary differential equation for L(v):

$$\frac{d^2L(v)}{dv^2}+\left(\frac{d\beta (v)}{dv}\right)^2L(v)=0$$
(54)
In this paper, we take $`\beta (v)`$ to be given as

$$\beta (v)=\frac{3}{10}e^{v^2/16}$$
(55)
Eq. 54 is then solved with a 4th order Runge-Kutta solver. Initial data is obtained by setting $`t=0`$ (and thus, $`v=z`$), shown in Figure 4. We present results of a von Neumann analysis about the point $`z=3.0`$ with $`\mathrm{\Delta }z=0.15`$. We have investigated different values of $`z`$ and $`\mathrm{\Delta }z`$, and find the results presented to be generic for different values of $`z`$ and $`\mathrm{\Delta }z`$, as long as we remain in the nonlinear regime, $`z<8`$.

The main difference between the stability properties of discretizations of the CT and ADM systems is shown in Figures 5 and 6. These show, respectively, plots of the maximum modulus of the eigenvalues (ignoring eigenvalues that are exactly 1) of the amplification matrices for discretizations of the ADM and CT systems at $`z=3.0`$. These plots show the range of the Courant factor $`0.3\sigma \mathrm{\Delta }t/\mathrm{\Delta }z0.9`$ and the range of modes $`0.7\pi \theta _kk_z\mathrm{\Delta }z\pi `$. Notice how, for these ranges of Courant factor $`\sigma `$ and mode wavenumber $`k_z`$, all of the eigenvalues for the CT system are less than or equal to 1, where as for the ADM system, there is always at least one eigenvalue that has a modulus that is greater than 1. For these ranges of $`\sigma `$ and $`k_z`$, we conclude that the CT system is stable, while the ADM system is unstable.

By looking at Courant factors of values $`\sigma >1`$, we see the typical Nyquist frequency instability, shown in Figures 7 and 8 for the ADM system and CT system, respectively. Notice in both Figures 7 and 8 that the largest modulus of the eigenvalues of the amplification matrices for long wavelength modes ($`\theta _kk_z\mathrm{\Delta }z<0.3\pi `$) are all greater than $`1`$. In fact, this is the case for all values of Courant factor $`\sigma `$. This is due to the fact that there is an exponentially growing gauge mode in the analytic solution to the analytic equations. Recall that we are using a different gauge choice than that given by the exact solution in Eqs. 53 \- 55. Using the gauge choice given in Eq. 27 and Eq. 42, there exists an exponentially growing gauge mode in $`K_{zz}`$. To take into account equations that admit exponentially growing solutions, the von Neumann condition, Eq. 8, must be modified (see  for details). In order that finite difference discretizations of equations that admit solutions that have exponentially growing modes remain stable, we must have

$$\left|\lambda _i\right|1+𝒪(\mathrm{\Delta }t).$$
(56)
To verify that the long wavelength ($`\theta _kk_z\mathrm{\Delta }z<0.3\pi `$) phenomena observed in Figures 7 and 8 is simply due to the existence of (long wavelength) exponentially growing modes in the analytic solution to the analytic equations, we repeat the calculations leading to Figures 7 and 8, but decrease the discretization parameter $`\mathrm{\Delta }z`$ by a factor of 2 ($`\mathrm{\Delta }z=0.15/2=0.075`$). The results for both the ADM and CT systems are similar, so we only show the results for the CT system in Figure 9. Notice that by decreasing $`\mathrm{\Delta }z`$ by a factor of 2 and holding the Courant factor $`\sigma `$ constant, we also decrease $`\mathrm{\Delta }t`$ by a factor of 2. By comparing the long wavelength ($`\theta _kk_z\mathrm{\Delta }z<0.3\pi `$) sections of Figures 8 and 9, we indeed see that the difference between the maximum modulus of the eigenvalues of the amplification matrices and 1 decreases by a factor of 2 when $`\mathrm{\Delta }t`$ is decreased by a factor of 2. Therefore, the discretizations of both the ADM and CT systems are stable for the long wavelength, exponentially growing gauge mode. Of course, the maximum eigenvalues for the high frequency sections of Figures 8 and 9 for $`\sigma >1`$ do not approach 1 as $`𝒪(\mathrm{\Delta }t)`$, which signifies a true von Neumann instability for $`\sigma >1`$.

## V Discussion and Conclusions

First, while the results presented in this paper are specific to the discretization method and initial data chosen here, one of the main points of this paper is to show that a von Neumann analysis can indeed be applied directly to discretizations of the Einstein equations. In the past, particular discretization methods were only analyzed with the von Neumann method through simple equations, such as the linear wave equation. Here we show that it is possible to carry out a von Neumann analysis on discretizations of equations as complicated as the Einstein equations. We would like to point out that a von Neumann analysis could also be used to test and/or formulate boundary conditions. The stability of outer boundary conditions, as well as the stability of inner boundary conditions for black hole evolutions, could be tested first through a von Neumann analysis instead of the traditional (and painful) method of coding up an implementation and looking for (and usually finding) numerical instabilities during the numerical evolution.

In this paper, we have shown that the stability properties, as determined by a von Neumann stability analysis, of a common discretization (a 2-iteration iterative Crank-Nicholson scheme) of the ADM and CT systems about flat space are similar to the stability properties of the scalar wave equation. However, as we would like to emphasize again, the stability properties of a nonlinear finite difference update operator depend on the values of the discrete evolution variables. Therefore, it is not enough to study the stability of numerical relativity codes about flat space. In principle, one must verify the stability of a nonlinear finite difference update operator about every discrete state encountered during the entire discrete evolution process, as argued in Section II. As a first step in this direction, we studied the stability properties of the ADM and CT systems about highly nonlinear plane waves by performing a von Neumann analysis for these scenarios. Several interesting features presented themselves.

The main difference between the stability of the ADM and CT systems about the nonlinear wave solution is seen in Figures 5 and 6. There, we see that, for a wide range of Courant factors $`\sigma =\mathrm{\Delta }t/\mathrm{\Delta }z`$ and wavenumbers $`k_z`$ (which includes the mode that is usually the most troublesome in numerical relativity, the Nyquist frequency mode whose wavelength is $`2\mathrm{\Delta }z`$), the CT system is stable (i.e., the amplification matrix has a spectral radius of 1 or less) whereas the ADM system is unstable (i.e., the amplification matrix has a spectral radius greater than 1). Note that the ADM system, in this region of Courant factor $`\sigma `$ and wavenumber $`k_z`$, has an amplification matrix whose spectral radius is approximately $`1+10^5`$. This is a very small departure from unity. Thus, instabilities arising from these modes could take a long time to develop during numerical evolutions. For example, the 2-iteration iterative Crank-Nicholson update operator for the scalar wave equation from Section II A has a spectral radius of $`1+10^5`$ for a Courant factor $`\sigma =\mathrm{\Delta }t/\mathrm{\Delta }z=2.000004`$, and this update operator can be used to evolve initial data sets many wavelengths before the Nyquist frequency instability sets in.

Again, it must be pointed out that these results are specific to plane wave spacetimes, simplified (diagonal) forms of the metric, choice of gauge, choice of initial data used here, and the 2-iteration iterative Crank-Nicholson scheme. It remains to be seen whether or not these results are generic to other, more general discretizations of the Einstein equations and for more general data, such as black hole and neutron star discrete evolutions. The one thing that is certain is that the von Neumann analysis can be used as a diagnostic tool for determining the stability of discretizations of nonlinear differential equations as complicated as the Einstein equations.

###### Acknowledgements.

We would like to thank Josh Goldberg, Alyssa Miller, David Rideout, Peter Saulson, Rafael Sorkin, and Wai-Mo Suen for encouragement and interesting discussions. This research is supported by NSF (PHY 96-00507 and PHY 99-79985) and NASA (NCCS5-153).


format: text

text:
# 1 Introduction

## 1 Introduction

Diffusion in a disordered system is a well-studied problem . Different approaches have been taken – such as models which consider inter-site barriers of varying heights, or models which picture sites for diffusing particles as wells with varying depths, also combinations of barrier and well models . Problems where the diffusing particle may encounter randomly distributed traps are also studied . However, one class of problems involve dynamic disorder, i.e. wells, barriers or traps whose distribution in space changes with time. In other words the location of the traps get redistributed at certain time intervals. Such models describe systems such as ion-conducting polymers, or glasses above the glass transition temperature .

Models with dynamic disorder have also been studied for some time. The ‘dynamic bond percolation model’ (DBPM) by the Northwestern University group, and its many variations study one aspect of the problem quite extensively, effective medium approaches have also been used .

A different formulation of the dynamic disorder problem, considering a ‘well’ model, rather than a ‘barrier’ model was suggested by Bhattacharyya and Tarafdar . Mandal et al  applied model  to explain experimental results on the PEO-NH4ClO4 system .The model  considers a random distribution of two types of sites on a square lattice, both having finite dwell times, i.e. non-zero jump rates. It is not a ‘percolation’ model in the sense that it never has zero diffusivity. The diffusion coefficient for different ratios of the two types of sites, and for different rearrangement times was calculated by computer simulation.

In the present work we take the rearranging lattice model with finite traps  to the limit where one type of site is an infinite trap. In this limit the diffusion coefficient as function of the trap concentration $`c`$ and rearrangement time $`\tau _r`$ are calculated in three-dimensions. The results are compared with the finite trap model as well as DBPM, and are supported by computer simulations.

## 2 The rearranging trap model

We consider a cubic lattice with a fraction $`c`$ of sites occupied by traps. A particular configuration of traps remains constant from $`t=0`$ to a time interval $`t=\tau _r`$. After that, the traps get rearranged, though the trap concentration $`c`$ is constant. If $`N_0`$ walkers started a random walk on the lattice from different sites at $`t=0`$, during time $`\tau _r`$ a certain number get trapped. The total square distance that would have been travelled by the $`N_0`$ walkers in the absence of traps is

$$(r_{0}^{}{}_{}{}^{2})_{tot}=2.d.N_0.D_0\tau _r$$
(1)
If there are traps this value gets reduced to say $`(r_{tr}^{}{}_{}{}^{2})_{tot}`$. Here $`D_0`$ is the diffusion coefficient with no traps, and $`d`$ is the dimension of the system.

After the interval $`\tau _r`$, since traps are rearranged, previously trapped walkers may be set free and resume their walk. In the next interval from $`t=\tau _r`$ to $`t=2\tau _r`$ they again cover a total distance$`(r_{tr}^{}{}_{}{}^{2})_{tot}`$. Fig 1 shows how the square distance covered increases with time in our simulation described in section 3. So for a time interval $`>>\tau _r`$, the diffusion coefficient on a rearranging lattice is given by

$$D_{tr}=(r_{tr}^{}{}_{}{}^{2})_{tot}/2.d.N_0.\tau _r$$
(2)
This argument is the same as in . The problem is thus reduced to calculating $`(r_{tr}^{}{}_{}{}^{2})_{tot}`$.With traps having infinite depth (i.e. infinite dwell time), this is estimated as follows, without involving any complicated mathematics.

### 2.1 Calculating the diffusion coefficient $`D(c,\tau _r)`$

Suppose $`N_0`$ particles start a random walk at $`t=0`$ on a cubic lattice occupied randomly by $`c`$ fraction of traps. The remaining $`(1c)`$ sites are ‘allowed’ sites. If a particle encounters a trap it gets stuck until the next rearrangement time , $`r_{tr}^{}{}_{}{}^{2}`$ is calculated as the mean distance travelled by each of the $`N_0`$ particles $`(r_{}^{2}{}_{tot}{}^{})/N_0`$. If there were no traps $`(r^2)_{tot}`$ would be given by equation 1. This is the total distance travelled by all the particles in time $`\tau _r`$. In presence of traps more and more particles get stuck and the total square distance travelled does not reach the value in eqn (1), but falls short by an amount $`\mathrm{\Delta }r`$.

$$\mathrm{\Delta }r=2.d.N_0.D_0.\tau _r(r_{tr}^2)_{tot}$$
(3)
knowing $`\mathrm{\Delta }r`$ one can calculate $`D(c,\tau _r)`$ from the relation

$$D(c,\tau _r)=(r_{tr})_{tot}^2/2.d.N_0.\tau _r$$
(4)
### 2.2 Calculation of $`\mathrm{\Delta }r`$ : the trapping law

To calculate $`\mathrm{\Delta }r`$ we assume a simple law for trapping.
We assume that the number of particles $`dN`$ trapped in the time interval $`t`$ to $`t+dt`$ is proportional to $`N`$ – the number of particles at time $`t`$, $`c`$ – the trap concentration and $`dt`$. So,

$$dN=a.c.N.dt$$
(5)
$`adS(t)/dt`$, where $`S(t)`$ is the number of distinct sites visited in time $`t`$ . For three or higher dimensions

$$S(t)t$$

, so $`a`$ is a constant depending on the dimension of the system but independent of $`t`$. We have thus the simple trapping rule

$$N=N_0.exp(a.c.t)$$
(6)
It may be mentioned here, that problems of absorption by traps, is well studied, and exact analytical results have been worked out for the asymptotic limit $`t\mathrm{}`$ . In the asymptotic limit the number of surviving particles at time t is given by

$$N=N_0.exp(a.c^{2/(d+2)}.t^{d/(d+2)})$$
(7)
where the exponents of $`c`$ and $`t`$ as well as $`a`$ are dimension dependent. However this result is valid only when the probability of survival becomes very small $`10^{13}`$ . In the regime we are interested in the simple expression eq(6) agrees much better with simulation results than eq (7) as we show in fig 2 for trapping in three-dimensions. We may now proceed to calculate $`D`$. The particles trapped at time $`t`$ to $`t=t+dt`$ contribute an amount $`\delta r`$ less to $`(r_{tr}^2)_{tot}`$ than they would if not trapped.

$$\delta r=2.d.(\tau _rt)$$
(8)
Henceforth we consider the case $`d=3`$ only.

Integrating this over the interval $`\tau _r`$ for all particles trapped during this time, we have

$$\mathrm{\Delta }r=_{t=0}^{\tau _r}\delta rdN=6_{0}^{}{}_{}{}^{\tau _r}a.c.N_0.exp(a.c.t)dt.(\tau _rt).D_0$$
(9)
or,

$$=6.D_0.N_0(\tau _r+exp(a.c.\tau _r)/a.c1/a.c)$$
(10)
Since

$$(r_{tr}^2)_{tot}=6.N_0.D_0\tau _r\mathrm{\Delta }r$$
(11)
the diffusion coefficient in the presence of traps with a rearrangement time $`\tau _r`$ is

$$D=(D_0/ac\tau _r)(1e^{ac\tau _r})(1c)$$
(12)
the factor $`(1c)`$ accounts for the particles which are trapped at the outset i.e. at $`t=0`$.

### 2.3 The limiting cases

For $`\tau _r0`$ that is for very rapid rearrangement

$$D/D_0=1c$$
(13)
as expected. For $`\tau _r\mathrm{}`$ that is when the system is effectively quenched, with no rearrangement

$$D/D_0=1$$
(14)
for $`c=0`$ and

$$D/D_0=0$$
(15)
for $`c>0`$. So there is a percolation threshold at $`c=0`$. Any nonzero trap concentration, however small gives $`D=0`$ in the limit $`\tau _r0`$. Fig 3 shows $`D`$ vs $`c`$ for different $`\tau _r`$. The value of $`a`$ is taken as $`a=.68`$, as given by our simulation results for the survival of random walkers in a three-dimensional lattice with traps.

## 3 Simulating diffusion on the rearranging trap model

We have simulated the random walk on a three-dimensional rearranging lattice, and calculated the diffusion coefficient for different $`c`$ and $`\tau _r`$. The algorithm is somewhat similar to the finite trap model in .

We work on an effectively infinite lattice, thus avoiding finite size effects. A random walker starts to walk on a cubic lattice of unit spacing with a concentration $`c`$ of traps. The location of the traps are not preassigned, the walker decides whether the current site is a trap from a random number, as it goes on. However, once assigned, the trap location remains fixed for time $`\tau _r`$ until the next rearrangement. Each trap is of infinite depth, so the walker encountering a trap has to remain there until the next rearrangement. To incorporate this in the algorithm, as soon as the walker falls in a trap, we freeze that walk upto the next rearrangement. To ensure that a site once assigned as ‘allowed’, is not seen as a trap on a subsequent visit, we maintain a list of all allowed sites visited during $`\tau _r`$, so that if these are revisited they are still assigned as ‘allowed’.

We average over $`10^5`$ walks to get $`r^2`$ for different time intervals $`t>>\tau _r`$.$`r^2`$ vs. $`t`$ are shown in fig 1, the curve is exactly like the results shown by . $`D(c,\tau _r)`$ calculated from the average slope of $`r^2`$ vs. $`t`$ are plotted against $`\tau _r`$ for different $`c`$ and shown in fig. 4. In fig. 5 we compare the simulated value of $`D`$ with the value calculated from equation 12.

## 4 Discussion

Figure 5 shows that the calculated results are always lower than the value obtained from our simulations. The reason is probably that our simplified calculation ignores the distribution of distances travelled by different walkers and assumes an equal average distance covered by all the walkers. However the agreement is quite satisfactory inspite of this.

The present model though closely related to the DBPM  and the finite trap model  shows features distinct from either of them. We discuss here the points of similarity and dissimilarity between the models.

Physically the difference between the finite trap model  and the present one is only that here the traps are infinite. At this limit however, a qualitative difference enters which is the appearance of the percolation threshold. Now, in the quenched limit $`\tau _r\mathrm{}`$, when the trap positions are frozen, the diffusion coefficient is zero for nonzero $`c`$, so the threshold is at $`c=0`$. The position of the threshold is at $`c=0`$, whatever be the dimension of the system. So, the difference with the model  is that there is a threshold. On the other hand the DBPM does have a percolation threshold $`p_c`$, which is at $`c=0`$ for 1-dimension, but for higher dimensions it is at the appropriate bond percolation threshold for that dimension.

Our simulation results shown in figure 1 resemble the DBPM results for $`p<p_c`$. So, the trapping effect is stronger here than in DBPM. This is expected since with sites as infinite traps, a random walker must fall in it given sufficient time, even if the trap concentration is very small. In the bond model ‘forbidden’ sites are inaccessible for the walker and may be avoided for $`c<p_c`$.

The finite trap model was shown to account successfully for the dynamic disorder in the polymer electrolyte $`PEONH_4Cl_4`$ . There the crystalline regions in the polymer are assumed to be the sites with longer dwell time, and the amorphous regions the highly conducting sites with smaller dwell time. This model is suited for systems with stronger trapping, and also systems where reacting species are diffusing in a viscous medium. For example – formation for excimers in fluorescence experiments, where fluorescing monomers combine with each other to form an ‘excited dimer’ or excimer . The excimer has a certain lifetime after which it may dissociate, i.e. the walker is released from the trap. However, in this case $`dNc^2`$ rather than $`c`$ in equation (5) if $`c`$ is the monomer concentration.

The primary interest of this model is that it is possible to calculate the dynamic disorder effect in a very simple way, compared to the DBPM which requires very involved mathematics, and the results are qualitatively very similar.

## 5 Acknowledgement

We thank UGC for financial assistance. S.Mandal and R.Dasgupta acknowledge CSIR for the award of Senior Research Fellowship and Research Associateship respectively. Authors also thank Dr. P. Nandy for useful discussions.

Figure Captions

Figure 1 : Plot of average squared distance versus time for c = .007 from computer simulation.

Figure 2 : Plot of fraction of walkers surviving versus time steps from simulation, equation (6) & equation (7) for c = .001 & c = .01

Figure 3 : Plot of calculated diffusion coefficient versus trap concentration for renewal time = 1, 10, 50, 250.

Figure 4 : Plot of diffusion coefficient versus renewal time for c = .001, c = .005, c = .01 & c = .05 from simulation.

Figure 5 : Comparison of diffusion coefficient versus trap concentration from simulation and theory \[equation (12)\] for renewal times = 250, 2000.