text:
# Tree–Particle–Mesh: an adaptive, efficient, and parallel code for collisionless cosmological simulation

## 1 Introduction

Numerical simulation of the nonlinear evolution of collisionless dark matter, usually by following a set of point masses as they move under their mutual gravitational influence, has, over the last few decades, played an important role in increasing our understanding of cosmological structure formation. Rather than attempt to make a summary of the extensive literature available on this subject here, we direct the reader to the recent reviews by Bertschinger (1998) and Klypin (2000).

A search of the internet will reveal a diverse set of N-body codes now available to be downloaded by those interested in carrying out cosmological simulations. One of the first to be developed is the direct code of Aarseth (1999). The operation count for implementing any such particle–particle code scales as $`N^2`$, where $`N`$ is the number of particles. One can reach larger $`N`$ and thus higher mass resolution (but spatial resolution limited by the cell size on a grid) in the same amount of CPU time using the PM (Particle–Mesh) method, for example the code by Klypin & Holtzman (1997). The operation count here scales as $`N`$log$`N`$, with a small prefactor; PM is implemented on regular grids, where efficient fast Fourier transform (FFT) algorithms are available. To attain higher, subgrid resolution one can add to a PM code particle–particle interactions and refined subgrids, such as in the Hydra code (Couchman, Thomas & Pearce, 1995). A gridless approach is used in tree codes, for example the ZENO package of Barnes (1998). The operation count in the Barnes–Hut tree algorithm (Barnes & Hut, 1986) also scales as $`N`$log$`N`$ but the prefactor is much larger— by a factor of 10–50 (Hernquist, 1987)— than for a PM code. The tree code of Hernquist (1987) has been made parallel (Becciani & Antonuccio-Delogu, 2001). A tree code with individual particle time steps called GADGET has recently been released in both serial and distributed memory parallel versions (Springel, Yoshida & White, 2001). One promising new development is the use of multigrid methods involving mesh refinements of arbitrary shape (Kravtsov, Klypin & Khokhlov, 1997), such as in MLAPM (Knebe, Green & Binney, 2001). Another cell-based approach has recently been developed which shows improved performance over standard tree codes (Dehnen, 2000, 2002) and is included as part of the NEMO package (Teuben, 1995).

For a given choice of algorithm, a related problem is making it run efficiently on parallel computers (Yahagi, Mori & Yoshii, 1999; Ricker, Dodelson & Lamb, 2000; Viturro & Carpintero, 2000; Lia & Carraro, 2001; Springel, Yoshida & White, 2001; Yahagi & Yoshii, 2001; Dorband, Hemsendorf & Merritt, 2002; Jing & Suto, 2002; Miocchi & Capuzzo-Dolcetta, 2002; Stadel, 2002). Gravity is long-range, so to compute the force on a given particle, one needs some information from every other particle, and communication is thus required if particles are distributed among many processors. Xu (1995) presented a new algorithm, called Tree–Particle–Mesh, based on using domain decomposition in a manner designed to overcome this difficulty. TPM uses the efficient PM method for long-range forces and a tree code for sub-grid resolution. Isolated, overdense regions are each treated with a separate tree, thus ensuring coarse–grained parallelism.

Bode, Ostriker, & Xu (2000), hereafter BOX, made several improvements in the code of Xu (1995). In this paper we present further refinements which improve the accuracy and efficiency of the algorithm. These include the following: allowing individual particle time steps within trees, an improved treatment of tidal forces, new criteria for higher force resolution and for the choice of time step, and parallel treatment of large trees. In §2 we present an overview of latest version of TPM, including time stepping; §3 discusses implementation of domain decomposition and the particle time step criterion, with an overview of how the code works in practice. The accuracy of the code is explored in §4 by comparing results with two other algorithms: P3M and the GADGET tree code. Concluding remarks and details of the public release of the TPM code are given in §5.

Calder et al. (2002) discuss the distinction between verifying a simulation (that is, being sure that the equations are being solved correctly) and validating it (having confidence that the equations and their solution actually resemble a real world problem). Obviously, in this paper only the former is done. Various aspects of the validity of cosmological simulations are addressed in Klypin (2000), Knebe et al. (2000), van Kampen (2000), Hamana, Yoshida & Suto (2002), Knebe (2002), Power et al. (2002), Baertschiger, Joyce & Labini (2002), and references therein. When using an N-body code, the simulator will need to be cognizant of these issues.

## 2 Overview of TPM

TPM begins with a standard PM code (in fact, if there are no particles selected to be in trees, then TPM defaults to a PM code); this portion is similar to the PM code described by Efstathiou et al. (1985). The density on a regular grid is found by interpolating particle positions using the cloud–in–cell (CIC) method, and Poisson’s equation is solved on this grid using Fast Fourier Transforms. As noted in the introduction, the FFT technique is highly efficient and scales as $`N`$log$`N`$ (Hockney & Eastwood, 1981). In order to obtain subgrid resolution a tree code is used. One possible manner of doing this is to compute shorter range forces for every particle with a tree code (Bagla, 1999), in a manner analogous to the P3M algorithm. Instead, TPM uses domain decomposition, taking advantage of the fact that potentials due to matter outside and inside a given domain can be added linearly, and it uses different solvers to find these two potentials. This section begins with a very broad description of the TPM algorithm, and then considers various aspects in more detail.

A basic outline of TPM can be summarized as follows:

1. Identify tree regions and the particles in each such region; particles outside of tree regions are evolved with PM only.
2. Push PM particles to midstep
3. Find the tidal potential in each tree region— the potential due to all mass outside of that region.
4. Integrate each tree region forward to the middle of the PM step.
5. Find the PM potential and update the acceleration for PM particles.
6. Integrate each tree region forward to the end of the PM step.
7. Update PM particle velocities and positions to the end of the time step.

The motivation behind this algorithm is that the motions of particles within a given tree region can be integrated independently using high temporal and spatial resolution: knowledge of the rest of the simulation is not required, except for the time–averaged, externally produced tidal forces.

The stepping in time of positions $`𝐱`$ and velocities $`𝐯`$ is accomplished with a standard second order leapfrog integration in comoving coordinates, the cosmological model determining the scale factor $`a(t)`$:

$`𝐱(t+\frac{1}{2}\mathrm{\Delta }t)`$ $`=`$ $`𝐱(t)+\frac{1}{2}𝐯(t)\mathrm{\Delta }t,`$ (1a)
$`𝐠(t+\frac{1}{2}\mathrm{\Delta }t)`$ $`=`$ $`\mathrm{\Phi }(𝐱),`$ (1b)
$`𝐯(t+\mathrm{\Delta }t)`$ $`=`$ $`{\displaystyle \frac{1H\mathrm{\Delta }t}{1+H\mathrm{\Delta }t}}𝐯(t)+{\displaystyle \frac{a^3}{1+H\mathrm{\Delta }t}}𝐠(t+\frac{1}{2}\mathrm{\Delta }t)\mathrm{\Delta }t,`$ (1c)
$`𝐱(t+\mathrm{\Delta }t)`$ $`=`$ $`𝐱(t+\frac{1}{2}\mathrm{\Delta }t)+\frac{1}{2}𝐯(t+\mathrm{\Delta }t)\mathrm{\Delta }t.`$ (1d)
Here $`a`$, $`H\dot{a}/a`$, and the gravitational potential $`\mathrm{\Phi }`$ are determined at time $`t+\frac{1}{2}\mathrm{\Delta }t`$.

One advantage of the TPM approach is that it allows the use of multiple time steps. Tree particles are required to take at least two steps per PM step, but each particle has an individual time step so that finer time resolution can be used if required. Particles are arranged in an hierarchy of time bins differing by a factor of two, in the manner of Hernquist & Katz (1989):

$$\mathrm{\Delta }t_{tree}=\frac{\mathrm{\Delta }t_{PM}}{2^s},$$
(2)
where $`\mathrm{\Delta }t_{PM}`$ is the PM time step and the integer $`s1`$. In a sense TPM operates along the same lines as a tree code, except that the particles in the $`s=0`$ time step bin are handled differently from the rest. A diagrammatic representation of the time stepping is shown in Fig.1 for the case when all tree particles are in the longest time step bin, $`s=1`$. Beginning at time $`t`$, the PM particle positions are moved forward to midstep (eq. \[1a\] with $`\mathrm{\Delta }t=\mathrm{\Delta }t_{PM}`$). The tidal potential is then found, and the tree particles are evolved forward one full step with $`\mathrm{\Delta }t=\mathrm{\Delta }t_{PM}/2`$. The PM potential is updated (eq. \[1b\]), and for a second time the tree particles are evolved for one full step, to time $`t+\mathrm{\Delta }t_{PM}`$. Finally, the PM particle velocities and positions are updated to the end of the step (eq. \[1c-d\]).

To repeat the TPM algorithm in more detail:

#### 1.

Identify tree regions and the particles in each region. This is done by identifying PM cells above a given density threshold, described in §3.1. Adjoining cells are then grouped into isolated tree regions, as described in BOX.

#### 2.

Push PM particles to midstep (eq. \[1a\]).

#### 3.

Find the tidal potential $`\mathrm{\Phi }_{ext}`$ for each tree region. First the total potential is computed on the grid in the standard PM manner. Then for each of the tree regions, a small portion of the grid containing the PM potential is saved, and the contribution to this potential made by the tree itself is subtracted out, leaving the tidal potential due to all external mass. See BOX for details.

Since $`\mathrm{\Phi }_{ext}`$ can only be calculated once per PM step, it will be calculated at midstep. PM particles are thus already in the proper location for this, but tree particles, with positions at the beginning of the PM time step, are not. However, since tree regions are spatially separated and the halo density profiles are evolving on slower timescales than any individual particle’s orbital period, an approximate position is sufficiently accurate. (Note that since a given tree’s own contribution to $`\mathrm{\Phi }_{ext}`$ is exactly subtracted back out, only the effect on other trees needs to be considered). Thus, an approximation of the tree particles’ positions at the middle of the PM step is taken to be the position one full particle time step ahead; since tree particle time steps are half the PM time step or less, this is no more than the PM midstep. With these advanced positions, the potential on the grid is found in the standard PM manner.

For each tree the PM potential in a cubical subvolume is saved. This volume is slightly larger than the active cell region (in order to allow for finding the gradient of the potential through finite differencing) plus one extra cell on a side in case particles migrate out of the active cell region during the time integration. Since the PM time step is limited by a Courant condition (see §3.2) such that PM particles cannot move more than a fraction of a cell per step, one extra cell is sufficient— particles near the edge of a tree region are at only a slightly higher overdensity than the densest PM regions, and hence have similar velocities. At the beginning of a tree step the portion of the PM potential due to the tree particles themselves is subtracted out, leaving the tidal potential. With $`\mathrm{\Phi }_{ext}`$ saved in this manner, every time a particle’s acceleration is updated, the tidal force is calculated from the grid in the same manner that PM forces are found. This is an improvement over the method used in BOX.

#### 4.

Integrate each tree region forward to the middle of the PM step. This is done with a tree code, adding in the tidal forces. The tree code we use was written by Lars Hernquist (Hernquist, 1987; Hernquist & Katz, 1989; Hernquist, 1990). Any other potential solver could be used, but a tree code is well suited for this type of problem, in that it can efficiently handle a wide variety of particle distributions, can include individual particle time steps, and scales as $`N`$log$`N`$.

Note that for each tree this portion of the code is self–contained; that is, once the particle data and potential mesh for a given tree have been collected together no further information is required to evolve that tree forward. This makes TPM well suited for parallel processing on distributed memory systems. Once the tree data is received by a given processor this step— which is the most computationally expensive part of the algorithm— can be performed without any further communication. Given this coarse–grained parallelism, one can use widely distributed and heterogeneously configured processors, with more capable processors (or groups of processors) reserved for the largest trees.

#### 5.

Find the density and potential on the PM grid again, and update the acceleration for PM particles (eq. \[1b\]). All tree particles are at the PM midstep, the proper location for computing the force on PM particles.

#### 6.

Integrate each tree region forward to the end of the PM step.

#### 7.

Update PM particle velocities and positions to the end of the time step (eqs. \[1c-d\]).

## 3 Other Aspects of the Implementation

### 3.1 Domain Decomposition

TPM exploits the fact that gravitational instability in an expanding universe creates spatially isolated high density peaks; for typical CDM models these peaks will contain a significant fraction of the mass ($`1/4`$) but occupy a tiny fraction of the volume ($`10^3`$; see for example Fig. 2 of Mo & White (2002)). As described in detail in BOX, this is accomplished by tagging active mesh cells based on the cell density (which has already been calculated in order to perform the PM step). Adjacent cells are linked together, dividing the cells into groups in a friends–of–friends manner. Thus isolated high–density regions of space are identified, separated by at least one cell width.

In BOX, the criterion used to decide if a cell $`i`$ was active was if its density $`\rho _i`$ exceeded a global threshold density

$$\rho _i\rho _{\mathrm{thr}}=A\overline{\rho }+B\sigma ,$$
(3)
where $`\overline{\rho }`$ is the mean density and $`\sigma `$ is the dispersion of cell densities. While temporally adaptive (in the sense that the threshold is lower at early times when peaks are rare), this criterion has two main drawbacks. First, small peaks in isolated regions or voids will not be picked up, especially when $`\sigma `$ is large, because the total mass in such a halo would not put a PM cell above $`\rho _{\mathrm{thr}}`$. Secondly, in the case of a region with peak density near $`\rho _{\mathrm{thr}}`$, whether or not a given PM cell is above $`\rho _{\mathrm{thr}}`$ can depend on the offset of the grid with respect to the particles— the halo mass may be divided into two or more cells.

We have developed an improved criterion for selecting tree regions involving the contrast of a cell with its surroundings (rather than a single global threshold value) found by comparing the density smoothed on two scales. For each cell three densities are computed. First $`\rho _s`$, found by boxcar smoothing over a length of 5 cells along each dimension. (Another type of smoothing, e.g. with a Gaussian, would be possible, but the boxcar can be done with minimal interprocessor communication). Second, $`\rho _c`$, the mean density in the eight–cell cube of which the cell under consideration is the lower octant (when using CIC, a particle in this cell would also assign mass to the other cells in this cube). Third, $`\rho _{\mathrm{ext}}`$, a measure of the density surrounding this cube, defined as

$$\rho _{\mathrm{ext}}=\frac{5^3\rho _s2^3\rho _c}{5^32^3}.$$
(4)
All eight cells used to find $`\rho _c`$ are marked as active if

$$\rho _c\rho _{\mathrm{thr}}=\left(A\overline{\rho }+B\sigma \right)\frac{\rho _{\mathrm{ext}}}{\overline{\rho }+\sigma +\rho _{\mathrm{ext}}}.$$
(5)
The reasons for such a choice of form can be seen by examining the behavior of eq.  in various limits. At early times in a cosmological simulation $`\rho \overline{\rho }`$ and $`\sigma \overline{\rho }`$, so $`\rho _{\mathrm{thr}}\frac{1}{2}A\rho _{\mathrm{ext}}\frac{1}{2}A\overline{\rho }`$. Thus with a choice of $`A=2`$ those regions which are only slightly overdense will be selected. At late times $`\sigma \overline{\rho }`$. In void regions $`\rho _{\mathrm{ext}}\sigma `$, which means $`\rho _{\mathrm{thr}}B\rho _{\mathrm{ext}}`$; thus an isolated cell will be chosen if it is a fixed multiple above its surroundings. If on the other hand $`\rho _{\mathrm{ext}}\sigma `$, then $`\rho _{\mathrm{thr}}B\sigma `$, meaning that at each step there is a given density above which all cells will be chosen; this limit will increase with time as $`\sigma `$ increases. Keep in mind that the distribution of cell densities will be approximately lognormal (Kayo, Taruya & Suto (2001) and references therein) when interpreting the value of $`\sigma `$.

After cell selection, the domain decomposition proceeds in the manner described in BOX. Cells are linked by friends–of–friends, yielding regions of space separated by at least one cell length. Any particle which contributes some mass to an active cell when finding the PM density is assigned to the corresponding tree.

Because the volume of all tree regions is less than one percent of the total volume, the amount of memory allocated to hold the tidal potential data is generally not significant compared to that already needed for the PM mesh. However, one complication of this scheme is that at early times, when $`\sigma `$ is small, the resulting trees can be very long filaments. While not containing a large amount of matter, being not very overdense and nearly one–dimensional, the size of a cubical volume enclosing such a tree will sometimes be a significant fraction of the entire simulation volume. This can cause difficulties because of the amount of computer memory required to compute and save the tidal potential in such subvolumes. The value of $`B`$ in eq.  is thus allowed to increase at early times when the spatial extent of trees tends to become larger, and decrease at later times when trees are more compact. We have settled on the following method, which has worked well in a variety of simulations. Space is allocated for a maximum subvolume length one quarter of the PM mesh size. The size of the largest sub-box used is checked at the end of each PM step. At early times (when $`\sigma <1.6`$) $`B`$ is increased by 0.5% if this box has a length more than a third of the allocated value. Once $`\sigma `$ exceeds 1.6 in a typical simulation, the spatial extent of trees has stopped growing, so in this case $`B`$ is decreased by 0.25% if the largest sub-box has a length less than half the maximum allocated. A minimum value is set for $`B`$ to prevent too many tree particles being chosen; a choice of $`B=10`$ will place roughly half of the particles in trees at $`z=0`$. An example of this in practice is discussed in §3.3.

The number of PM grid cells can be more or less than the number of particles, though of course the finer the grid, the greater the accuracy of the PM and tidal forces. Generally, the number of cells should at least equal the number of particles; we have found that eight grid cells per particle works well.

### 3.2 Time Step Criterion

A variety of methods for determining the particle time step have been proposed. This diversity of criteria reflects the fact that the appropriate time step can depend on a number of quite different considerations. These include the local density or the local dynamical time, nearby substructure, the nearest neighbor, and the softening length (Knebe et al., 2000; Power et al., 2002). In addition, for any single criterion, one can imagine a special circumstance where it breaks down. This suggests the possibility of basing $`\mathrm{\Delta }t`$ on more than one criterion. Thus in the TPM code we have decided to combine a number of possible criteria, drawn from variables which do not require any significant computation to find, in the following manner:

$`\mathrm{\Delta }t`$ $`=`$ $`\eta \mathrm{MAX}(\mathrm{\Delta }t_1,\mathrm{\Delta }t_2),\mathrm{where}`$ (6a)
$`\mathrm{\Delta }t_1`$ $`=`$ $`\mathrm{MIN}(ϵ/v,\sqrt{a^3ϵ/g}),\mathrm{and}`$ (6b)
$`\mathrm{\Delta }t_2`$ $`=`$ $`\mathrm{MIN}(v/g,g/\dot{g},{\displaystyle \frac{ϵ^2+r^2}{rv}}).`$ (6c)
Here $`\eta `$ is an adjustable dimensionless parameter, $`ϵ`$ is the spline kernel softening length, and $`r`$ is the distance to the nearest neighboring particle. The latter is found as a byproduct of the force calculation; if at the beginning of a tree step $`r`$ has not been calculated, one can as a proxy find the minimum $`r`$ that is consistent with the particle’s current $`\mathrm{\Delta }t`$. The time derivative of $`g`$ is approximated by $`\dot{g}=|g(t)g(t\mathrm{\Delta }t_s)|/\mathrm{\Delta }t_s`$. This criterion for $`\mathrm{\Delta }t`$ is quite conservative, and is similar in principle to that adopted by Aarseth (1999).

A comparison of the time step returned by eq.  and the often used criterion $`\eta \sqrt{ϵ/g}`$ is given in Fig.2, which shows how these two criteria compare at $`z=0`$ for particles in the largest tree of the simulation discussed in §3.3 and §4.2. This tree is made up of over $`1.4\times 10^5`$ particles, and contains two large halos plus a number of smaller satellites and infalling matter. Each contour level in Fig.2 encloses a tenth of the particles, and the remaining tenth are plotted as points. It can be seen that eq.  tends to yield a smaller timestep than $`\eta \sqrt{ϵ/g}`$; this is the case for 74% of the particles. Generally the differences are not large— less than a factor of two for 60% of the particles. Another trend seen in the figure is that as $`\eta \sqrt{ϵ/g}`$ becomes smaller it is more likely for eq.  to give an even smaller value, and conversely eq.  tends to give a longer time step than $`\eta \sqrt{ϵ/g}`$ as the latter becomes larger. In most cases the value of $`\mathrm{\Delta }t`$ is set by the factor $`(ϵ^2+r^2)/(rv)=(r/v)(1+ϵ^2/r^2)`$. In other words, the distance to the nearest neighbor is not allowed to change by a large factor, unless this distance is small compared to the softening length. Even in cases when this factor is larger than $`\sqrt{ϵ/g}`$, it is still smaller than $`v/g`$ and $`g/\dot{g}`$; since the equations being integrated have the form $`\mathrm{\Delta }x=v\mathrm{\Delta }t`$ and $`\mathrm{\Delta }v=g\mathrm{\Delta }t`$ this limits integration errors.

The PM time step is limited by a Courant condition, such that no PM particle moves more than one quarter of a cell size per step. Also, $`a`$ is not allowed to change by more than 1% per step. Initially this latter criterion is the most restrictive, but over time it allows a longer $`\mathrm{\Delta }t_{PM}`$ and becomes unimportant. The $`\mathrm{\Delta }t_{PM}`$ allowed by the Courant condition also tends to increase over time, although more slowly. This is because particles falling into dense regions are placed into trees, and the velocities of the remaining particles are redshifted (eq. \[1c\]). The largest time step for any particle, $`\mathrm{\Delta }t_{PM}`$, is kept constant until it can safely be increased by a factor of 2, and it is then doubled. The time step for tree particles is unchanged, which means they take twice as many steps per PM step as before; during tree evolution particles will move into a lower time step bin if allowed by eq. . Whenever any particle changes time step, its position is updated as described in Hernquist & Katz (1989) to preserve second order accuracy.

### 3.3 A typical run

Various aspects of a typical TPM run can be demonstrated with a standard cosmological simulation. This test case contains $`N`$=1283 particles in a box $`L`$=40$`h^1`$Mpc on a side. The number of grid points for the PM and domain decomposition portions of the code is eight times the number of particles, or $`256^3`$. The initial conditions were generated using the publicly available codes GRAFIC1 and LINGERS111These codes are available at http://arcturus.mit.edu/grafic/ (Bertschinger, 2001; Ma & Bertschinger, 1995). A spatially flat LCDM model was chosen, with cosmological parameters close to the concordance model of Ostriker & Steinhardt (1995): $`\mathrm{\Omega }_m=0.3,\mathrm{\Omega }_\mathrm{\Lambda }=0.7,h=0.70,\mathrm{\Omega }_bh^2=0.20,n=1`$, and $`\sigma _8=0.95`$. The particle mass is thus $`2.54\times 10^9h^1M_{}`$. Within the tree portion of the code, the opening parameter $`\theta `$ in the standard Barnes–Hut algorithm (Barnes & Hut, 1986) is set to $`\theta =0.5771/\sqrt{3}`$, and the time step parameter $`\eta `$=0.3 (see eq. ). The cubic spline softening length was chosen to be 3.2$`h^1`$kpc so that the spatial dynamic range $`L/ϵ10^4`$. With such a small softening length, halos with fewer than 150 particles are likely to undergo some 2-body relaxation in their cores over a Hubble time (assuming the particle distribution follows an NFW profile with concentration $`c`$=12).

Fig.3 displays the evolution of several important quantities during this run as a function of expansion parameter $`a`$. The top panel shows the dispersion of PM cell densities $`\sigma `$ (where the mean cell density is unity), and the second panel shows the value of $`B`$; both of these factor in eq. , which in turn plays a major part in determining the tree distribution. Various aspects of this distribution are shown in the remaining curves. The third panel shows the size of the largest cubical subvolume required (in units of PM grid cells). Initially small, this grows extremely rapidly as $`\sigma `$ rises from its initial value of 0.3 to $`1`$; in response $`B`$ is increased. When this size is at its greatest (at $`\sigma 1.6`$), the percentages of total volume and mass in trees (shown in the next two panels) are still quite small. The trees at this time tend to follow caustics— they are only slightly overdense and not very massive, but because of their filamentary nature they can have a large spatial extent. These caustics then fragment and collapse, so— even while the total mass in trees and the number of trees (shown in the third panel from the bottom) increase— the maximum subvolume size decreases. The changes in $`B`$ affect mainly the maximum subvolume size; the number of trees and volume contained in tree regions both increase at a steady rate up to $`z`$=1. The number of tree particles increases monotonically throughout the simulation.

After $`z=1`$ (by far the bulk of the computational time), the characteristics of the simulation change much more slowly. Roughly half the particles are in $``$2700 trees ranging from a few to $`10^5`$ particles but occupying only 0.4% of the simulation volume. The penultimate panel shows the number of particles in the first and third largest trees; by $`z=1`$ the growth in mass of these objects has slowed to a fairly constant, low rate. Trees are distributed in mass roughly as a power law, with 67% of the trees having fewer than 100 particles and 95% less than 1000. Most tree particles take two or four steps per PM step, but some are taking up to 64; there were 1075 PM steps total in this run.

As $`\sigma `$ increases, so does the limiting density above which all cells are treated at full resolution. It is important to keep track of this limit when interpreting the results of a TPM run. The bottom panel shows the highest value of the density among those cells evolved only with PM. By $`z=1`$ this density is 120 times the mean, corresponding to 15 particles inside a cell. By the end of the computation, the densest PM–only cell contains 26 particles; so trying to make statements about objects smaller than this will be complicated by the varying spatial resolution of TPM. Note that most objects with $`N25`$ will in fact be followed at full resolution— substructures inside larger objects because they are in regions of higher density, and small isolated halos because they present a density contrast to their surroundings. To be cautious we will limit our analysis to objects 50% larger than this limit, i.e. 40 particles or $`10^{11}h^1M_{}`$.

### 3.4 Load balancing

As seen in the previous section, during later epochs— which take up most of the computational time needed for a run— the mass distribution of trees is generally well fit by a power law ranging from a few particles up to of order 0.01$`N`$. The actual mass of the largest tree will depend on the ratio of the largest non-linear scale in the box to the box size; as this ratio becomes larger so does the mass of the largest tree. Furthermore, the amount of computation required for a tree with $`N_t`$ particles will scale as $`N_t`$log$`N_t`$. Thus, there can be a substantial variation in the computational time required between different trees, and evolving the largest tree can comprise a significant fraction of the total computational load. An efficient parallel code must handle this situation well when dividing work up among NCPU processors.

Load balancing of trees is achieved in three ways. First, trees are sorted by $`N_t`$ and then divided into “copses” of roughly equal amounts of work using the “greedy” algorithm. That is, starting with the largest tree, each one in turn is given to the copse which up to that point has been assigned the least amount of total work. Usually there is one copse per CPU, but there can be two or more per CPU if required by space constraints. There is a communication step, when all data associated with a copse is sent to one CPU. A CPU then evolves each tree in its local copse in turn, starting with the most massive. Additionally, the largest trees can be done in parallel. If a few large trees dominate the work load, then it is impossible for all copses to contain equal amounts of work— ideally, one would want NCPU copses, each comprising 1/NCPU of the total amount of work, but this clearly can’t happen if the largest tree takes a larger fraction just by itself. In this case, the force calculation for these large trees is done in parallel by a small number of CPUs. This is currently done quite crudely, with only the tree walk and force calculation actually carried out in parallel; in theory a fully parallel tree code could be used for every tree, allocating more processors to those copses containing the most work. Only a few nodes are usually required to reduce the time spent on the largest tree to the level required for load balancing. As a final means of balancing the load, when a CPU finishes evolving all the trees in its local copse, it then sends a signal to the other CPUs. A CPU with work still left will send an unevolved tree to the idle CPU for it to carry out the evolution.

The scaling of the current implementation of TPM with number of processors is shown in Fig.4. The test case for these timing runs is a standard LCDM model at redshift $`z`$=0.16 with $`N`$=2563 particles in a 320$`h^1`$Mpc cube. This particular set of runs was carried out on an IA-64 Linux cluster at NCSA named “Titan”. This machine consists of 128 nodes, each with dual Intel 800MHz Itanium processors, and a Myrinet network interconnect. While the total time required depends on processor performance, similar scaling with NCPU has been found on a number of machines with various types of processors and interconnects. The topmost line in Fig.4 is total time, calculated as the number of seconds wallclock time per PM step multiplied by NCPU; perfect scaling would be a horizontal line. TPM performs quite well— at NCPU=64 the efficiency is still 91% as compared to NCPU=4. Beyond this point scaling begins to degrade, with the efficiency dropping to 75% for NCPU=128.

The reason TPM scales well can be seen in the second line from the top, which shows the total time spent in tree evolution. This part of the code takes up most of the CPU time, but it requires no communication and there are enough trees so that just coarse–grained parallelization works reasonably well. The next two curves shown indicate the amount of time in the PM portion of the code and overhead related to trees (identifying tree regions and particles, etc.). These take a small fraction of the total time and scale well since the grid is distributed across all processors. The final curve in Fig.4 shows the time related to communication and imbalance in the tree part of the code. As NCPU increases it becomes more difficult to divide the tree work evenly and processors spend more time either recruiting work from others or waiting for them to finish. This overhead could likely be reduced by incorporating a fully parallel tree code and also by computing the nonperiodic FFT (required to find the tidal potential— see BOX) in parallel. A run with a larger number of particles and grid points, and thus a larger number of trees, would show efficient scaling beyond NCPU=128.

## 4 Performance and Accuracy

### 4.1 Spherical Overdensity Test

As a test of new code we have carried out a simulation of the secondary infall and accretion onto an initially uniform overdensity. This can be compared both to other codes and to the analytic, self-similar solution of Bertschinger (1985). To create the initial conditions, $`64^3`$ particles were placed on a uniform grid with zero velocity. The eight particles on the grid corners were then removed and placed at the center of the volume with a spacing one half that of the regular grid. Thus, this is actually more of a cubic overdensity than spherical, but it quickly collapses and the subsequent infall is independent of the details of the initial state. Also, there is a void located half the box distance away from the overdensity, but the evolution is not carried out for a long enough time for this to be significant. The initial condition is integrated from expansion factor $`a=0.01`$ to $`a=1`$. The values $`A=2`$ and $`B=10`$ were used for the constants in eq. . Since only one halo is forming in the box, $`\sigma `$ remains small, rising to only 0.5 by the end of the run. At the end, there are roughly 500 particles within the turnaround radius. This test can be seen as exploring how well the initial collapse of small objects is followed, which is the beginning stage of halo formation in an hierarchical scenario.

The final state of this run is shown in Fig.5. The top panel shows density as a function of radius. Filled points are the TPM model; error bars are the square root of the number of particles in a given radial bin. The inner edge of the innermost bin is twice the softening length. Also shown is the same run carried out with a P3M code (as open circles; details of this code are discussed in §4.2), and the solution of Bertschinger (1985). The agreement is quite good— the main limitation of this run is likely the mass resolution. As a measure of the phase space density, the lower panel shows $`\rho /v^3`$, where $`v`$ is the velocity dispersion of the particles in each radial bin. Following Taylor & Navarro (2001), we compare this to the Bertschinger (1985) solution for a $`\gamma =5/3`$ gas. Again the agreement is quite reasonable. Differences between TPM and P3M arise from different types of softening (P3M uses Plummer softening; here $`ϵ`$ was set to half the TPM value) and from time stepping (for P3M, $`\mathrm{\Delta }t`$ for all particles was set by the minimum $`\sqrt{0.05ϵ/g}`$ ).

### 4.2 Comparison with other codes

As a test of TPM in a less idealized situation, the initial conditions described in §3.3 were again evolved, but using two other N-body codes. The first is the P3M code of Ferrell & Bertschinger (1994) in a version made parallel by Frederic (1997)— this code was also used in §4.1. As in the TPM run, a mesh of $`256^3`$ grid points was used. This code uses Plummer softening; the softening length was set to half the value used in the TPM run. The time variable in this code is $`\tau `$, defined by $`d\tau =dt/a^2(t)`$; all particles have the same time step, set by the minimum $`\sqrt{0.05ϵ/g}`$. The other N-body code used is the tree code named GADGET222Publically available at http://www.MPA-Garching.MPG.DE/gadget/ of Springel, Yoshida & White (2001). In this code, periodic boundary conditions are handled with Ewald summation, the time variable is the expansion factor $`a`$, and each particle has an individual time step. Following Power et al. (2002), the conservative tree node opening criterion flagged by -DBMAX was used, and the time step was set by $`\sqrt{0.03ϵ/g}`$. In addition, the maximum allowed time step was set to 1% of the initial expansion factor; whenever $`a`$ doubled this maximum was also increased by a factor of two, by checkpointing the run and restarting with the new value. The softening length was set to be the same as the TPM run (allowing for different notational conventions).

As a first comparison of the codes, the two point correlation function $`\xi (r)`$, found by counting the number of particle pairs in bins of separation $`r`$, was calculated. To compute $`\xi (r)`$, 51 logarithmically spaced bins were used, with the minimum pair separation considered being $`2ϵ=6.4h^1`$kpc and the maximum one third of the box size. The results at various redshifts are shown for all three runs in Fig.6. Clearly there is little difference to be seen in the three codes. At larger scales this is to be expected; at smaller scales force approximations, smoothing, and relaxation may become important. The fact that P3M uses Plummer rather than spline softening explains why it shows a lower $`\xi (r)`$ at scales less than a few times the softening length. If high values of $`\xi (r)`$ provide a measure of accuracy, the TPM and GADGET are slightly more accurate than P3M.

The statistics of mass peaks— the dark halos surrounding galaxies, clusters, and so on— are an important product of N-body codes (White, 2002). One commonly used halo finder is friends–of–friends, or FOF333We employed the University of Washington NASA HPCC ESS group’s FOF code, available at http://www-hpcc.astro.washington.edu/ (Davis et al., 1985). The cumulative mass function, found with FOF using a linking length of 0.2 times the mean interparticle separation, is shown in Fig.7. At higher redshifts there is little discernible difference between the three codes. At later times, TPM seems to have fewer halos containing $``$50 particles. It is unclear how to interpret this, because halos with less than 150 particles (i.e. mass $`<3.8\times 10^{11}h^1M_{}`$) are affected by two-body relaxation. Knebe, Green & Binney (2001) found that an adaptive P3M and GADGET both produced more small FOF halos than the adaptive multigrid MLAPM and ART codes, so in this case TPM may agree more closely with the latter.

To investigate further a different halo finder was used, namely the Bound Density Maxima (BDM) algorithm444Publically available at http://astro.nmsu.edu/$$aklypin/ (Klypin et al., 1999). This algorithm is significantly different from FOF, so the resulting mass function is probing different qualities of the dark matter distribution. Density maxima inside spheres of radius 100 $`h^1`$kpc were found, and halo centers were then calculated using spheres of radius 40 $`h^1`$kpc. The mass for each halo is taken to be that inside a sphere containing the overdensity expected for a virialized halo just collapsed from a spherical top–hat; to calculate this overdensity the fit of Bryan & Norman (1998) was used. In BDM (unlike FOF) particles moving faster than the escape velocity are removed from the halo. The BDM cumulative mass functions for the three N-body codes are shown in Fig.8. As with FOF, the agreement between codes is quite good. In this case it appears to be TPM which produces more small halos, but again this effect is for low mass halos likely affected by relaxation.

Using the BDM halos with more than 40 particles, the halo–halo correlation function is shown in Fig.9. The number of pairs of halos in 21 logarithmically spaced radial bins, with the smallest separation being 200 $`h^1`$kpc and the largest one third of the box size, were tallied and compared with the expectation for a random distribution. All three codes give the same result. From the results presented so far in this section it is clear that TPM yields a quite similar evolved matter distribution as compared to other codes.

Finally we turn to the internal properties of halos. One probe of the mass distribution in a halo is the circular velocity $`v_c\sqrt{GM(<r)/r}`$. For each simulation we divide halos up into six mass bins. To avoid relaxation effects, the lowest mass considered is $`3.81\times 10^{11}h^1M_{}`$. The width of each bin is a factor of $`\sqrt{10}`$, so the largest mass bin contains two halos with mass above $`1.20\times 10^{14}h^1M_{}`$. For each set of halos the average circular velocity as a function of radius was calculated; these average velocity curves are shown in Fig.10. It can be seen that these curves are quite similar for all three codes. Two main trends are noticeable, the main one being that P3M shows lower circular velocities at small radii. This is likely due to the use of Plummer softening and less strict time step criterion in the P3M run, both of which would lead to lower densities in the innermost parts of halos. The other main difference is that for the lowest masses it appears that TPM yields higher circular velocities than the other two codes.

Instead of averaging over a number of halos of similar mass, it is also possible to compare halos on a one-by-one basis. To find the appropriate pairs of halos, the list of halos with more than 150 particles found by BDM was sorted by mass, and for each target GADGET or P3M halo a TPM halo was selected as a match if its mass was within 25% and its position within 1 $`h^1`$Mpc of the target halo; if more than one halo passed this test then the nearest in position was selected. The selected TPM halo was removed from further consideration, and the process repeated for next target halo. In this manner a match was found for 588 out of 609 GADGET halos and 573 out of 582 P3M halos. For various measured halo properties, the percentage difference of the TPM halo from the GADGET or P3M halo was calculated. Since the dispersion of these differences increases as less massive halos are considered, halo pairs are split into two groups, with the higher mass group containing all target halos with 394 or more particles (i.e. mass $`10^{12}h^1M_{}`$).

Results are shown in Table 1, which for each property gives the mean percentage difference and one standard deviation, as well as the first quartile, median percentage difference, and third quartile. The first property shown is halo mass $`M`$. The difference here is constrained to be less than 25%, but in most cases the TPM value is within 10% of the other code’s. It appears that P3M yields slightly lower masses than the other two codes, but the latter two agree quite well. The agreement in masses shows that the simple scheme used to find matching halo pairs works well, as does the fact that the difference in position is less than 165$`h^1`$kpc in 95% of the cases. This is reinforced by the agreement between the halo center of mass velocities $`v_{cm}`$ shown in Table 1; the velocity vectors are closely aligned, differing by less than $`7.5\mathrm{°}`$ in 95% of the pairs. Thus only a few percent of the halo pairs are mismatches; since these pairs still have similar masses and are likely in similar environments, they are kept in the comparisons.

In terms of the 3-D root mean square velocity $`v_{rms}`$ and maximum circular velocity $`v_c`$, the more massive halos are very similar in all three codes, with no offset between codes. However, at the lower mass end TPM gives values that tend to be systematically higher by a few percent. This can also be seen for the average circular velocity in the bottom curve of Fig.10; this is curve is the average for halos below roughly $`10^{12}h^1M_{}`$.

The difference between TPM and the other two codes is most pronounced near halo centers. The last property shown in Table 1 is a comparison of central density $`\rho _o`$, which is computed by measuring the amount of mass within $`4ϵ`$=12.8$`h^1`$kpc of the halo center. As would be expected from Fig.10, TPM yields higher central densities than GADGET, with P3M giving lower $`\rho _o`$ than either of these.

One possible source of the differences seen in TPM halos as compared to a pure tree code could simply be the choice of time step. The P3M run took $`2.1\times 10^4`$ steps and the GADGET run $`2.8\times 10^4`$; TPM on the other hand took $`4.4\times 10^4`$, or 60% more than GADGET. Of course this comparison is not entirely straightforward because different particles determine the smallest time step at different times. Part of the difference between GADGET and TPM may be due to the fact that TPM bins time steps by factors of two, whereas GADGET allows time steps to vary more gradually. A second TPM run was carried out to separate other code differences from the time step criterion. This run was identical to the first except that the time step was set not by eq. , but rather by $`\sqrt{0.05a^3ϵ/g}`$ (similar to the GADGET and P3M codes); as a result it took fewer steps than any of the other runs. There is no significant difference in the mass functions and halo–halo correlation functions between this run and the original TPM run. The particle-particle correlation function is different, however. This can be seen in Fig.11, which shows the ratio of the original TPM run’s $`\xi (r)`$ to that of the new run, at redshift $`z`$=0. The new TPM run has a lower $`\xi (r)`$ for $`r<20h^1`$kpc, roughly 5-10% lower than the original TPM and similar to the P3M run. This indicates the internal structure of halos has been affected by the longer time steps, which is confirmed by repeating the comparison of individual halos. The difference between halos in the new TPM run and the GADGET run are given in Table 1. For the higher mass halos the new run, unlike the original one, tends to give lower maximum circular velocities and central densities than the tree code. This indicates that the longer time step is in fact leading to inaccuracies. The differences between TPM and the tree code for low mass halos seen in the original run persist in the new run, though they are not quite as pronounced.

Another point of comparison between codes is efficiency with which they use computing resources. Fig.12 shows the wall-clock time consumed by each of the codes, as a function of expansion parameter, in carrying out the test simulation. All the runs used four 300 MHZ IP27 processors of an SGI Origin 2000; the codes were compiled with the SGI compilers and MPI library. All three codes used roughly the same amount of memory (TPM requires at least 20 reals per particle plus 3 reals per mesh point, divided evenly among processors). At the earliest times both P3M and TPM spend most of their time in the PM FFT, and so they behave similarly. However, as objects begin to collapse P3M begins to consume more time computing particle-particle interactions. The fact that the accelerations of all particles are updated every step also makes P3M use more time than do the two multiple time step codes. The tree code takes more time when the particle distribution is nearly homogeneous, demonstrating that a PM code (which is what the gridded codes basically are in this situation) is very efficient. However, the tree code timing is roughly independent of the particle distribution, and once inhomogeneity develops it does not require more time per update, whereas the other codes do. TPM does well compared to the tree code for a couple of reasons. By $`a`$=1 roughly half the particles in the TPM run are still being handled solely by PM. Also, imagine breaking $`N`$ particles up into $`t`$ trees each with $`N/t`$ particles; the time required to solve all the trees will scale as $`N\mathrm{log}N/t`$, lower by a logarithmic factor than using the same tree solver on all the particles. Overall, despite taking more timesteps, in this test case TPM required less CPU time than the other codes, by a factor of $`4`$ relative to the tree code and $`8`$ relative to P3M.

## 5 Conclusions

This paper has presented a parallel implementation of the TPM algorithm. Several improvements over the implementation of BOX have been made. Particles in tree regions each have an individual time step, half of the PM time step or less, making the tree integration more efficient. The treatment of tidal forces on trees is also improved by saving the tidal potential on a grid and evaluating the force as a function of particle position at each smaller particle time step; thus a greater ratio of tree to PM time steps is allowed. A new, more stringent, time step criterion has been implemented. A new criterion for locating regions for treatment with higher resolution is given; by finding cells with higher density than their surroundings, small halos in lower density regions are located and followed at full resolution. These changes significantly increase the speed and accuracy of TPM: a computation of the test case discussed in §4.2 using the code described in BOX required more CPU time (by a factor of over three) than the current version, but was only accurate for halos with more than 315 particles.

Comparisons with other widely used algorithms were made for a typical cosmological structure formation simulation. These show excellent agreement. The particle-particle correlation functions, the halo mass functions, and the halo-halo correlation functions from TPM agree quite well with those from P3M and tree codes. The internal properties of halos also agree; the main difference being that, for lower mass halos, TPM yields higher $`rms`$ and maximum circular velocities (by a few percent over a tree code). TPM halos also show higher central densities than those of the other two codes, though the mean difference is smaller than the dispersion. This difference disappeared, at least among the more massive halos, when the time step criterion of eq.  was replaced with one similar to that employed by the other two codes. Thus we conclude that a choice of a relatively conservative time step criterion contributed to a slightly improved accuracy.

TPM yielded results of similar accuracy to the other codes used here while using significantly less computational time (about a quarter of that needed by a tree code and an eighth of that needed by P3M). It also scales well on distributed memory parallel machines, such as networked PCs, because this parallelism is built in as part of the design of TPM. However, in committing to an algorithm which accentuates the coarse-grained parallelism inherent in a typical cosmological simulation, a large degree of flexibility is sacrificed. A basic presumption of TPM is that the largest nonlinear structure inside the simulation box is a small fraction of the total mass and volume. To simulate a situation where this is not the case (e.g. two colliding galaxies) another code would be preferred.

The TPM source code can be obtained at http://astro.princeton.edu/$$bode/TPM or by contacting the authors. The code is written in Fortran 77 and uses MPI for message passing; thus it is very portable and can be used on clustered PC’s or other distributed memory systems.

Many thanks are due to Lars Hernquist for generously supplying a copy of his tree code; also Edmund Bertschinger for use of his P3M code, Joe Henawi for help with FOF, and Scott Tremaine for useful discussions. This research was supported by the National Computational Science Alliance under NSF Cooperative Agreement ASC97-40300, PACI Subaward 766. Computer time was provided by NCSA and the Pittsburgh Supercomputing Center.


format: text

text:
# Power-law Bounds on Transfer Matrices and Quantum Dynamics in One Dimension II.

## 1. Introduction

Consider a discrete one-dimensional Schrödinger operator

(1)
$$H\psi (n)=\psi (n1)+\psi (n+1)+V(n)\psi (n),$$

in $`\mathrm{}^2()`$ or $`\mathrm{}^2()`$ (with a Dirichlet boundary condition). We are interested in proving lower bounds on the spreading of an initially localized wavepacket under the dynamics governed by $`H`$. That is, if we consider the initial state $`\psi `$, we ask how fast $`\psi (t)=\mathrm{exp}(itH)\psi `$ spreads out. One is normally interested in initial states that are well localized. In the present paper we shall limit our attention to the case $`\psi =\delta _1`$.

A typical quantity that is considered to measure the spreading of $`\psi (t)`$ is the following: Define

(2)
$$|X|_\psi ^p(T)=\underset{n}{}|n|^pa(n,T),$$

where

(3)
$$a(n,T)=\frac{1}{T}_0^+\mathrm{}e^{2t/T}|\delta _n,\psi (t)|^2𝑑t.$$

Clearly, the faster $`|X|_\psi ^p(T)`$ grows, the faster $`\psi (t)`$ spreads out, at least averaged in time. One typically wants to prove power-law lower bounds on $`|X|_\psi ^p(T)`$ and hence it is natural to define the following quantity: For $`p>0`$, define the lower growth exponent $`\beta _\psi ^{}(p)`$ by

$$\beta _\psi ^{}(p)=\underset{T+\mathrm{}}{lim\; inf}\frac{\mathrm{log}|X|_\psi ^p(T)}{\mathrm{log}T}.$$

There are presently two distinct approaches to proving lower bounds for $`\beta _\psi ^{}(p)`$. The first goes back to works of Guarneri , Combes , and Last  and is based on a study of the Hausdorff dimension of the spectral measure $`\mu _\psi `$ associated with the pair $`(H,\psi )`$. Namely, we have the following bound:

(4)
$$\beta _\psi ^{}(p)p\mathrm{dim}_H(\mu _\psi ).$$

The Jitomirskaya-Last extension  of Gilbert-Pearson theory  allows for a convenient way of investigating $`\mathrm{dim}_H(\mu _\psi )`$ and hence this approach has enjoyed some popularity (see, e.g.,  for applications).

On the other hand, this bound clearly gives nothing in the case of a zero-dimensional spectral measure, for example, in the case of a pure point measure. There are a number of models where one expects (or can prove) pure point spectrum with strictly positive values for $`\beta _\psi ^{}(p)`$. An example is given by the random dimer model; studied, for example, in . It is therefore desirable to have a way of proving lower bounds on the transport exponents which works for such models and, of course, whose input is easy to verify in concrete cases. Such an approach was developed in  (and employed in  to prove the conjectured dynamical lower bound for the random dimer model), and the present article is a continuation of that paper. The necessary input are power-law upper bounds on transfer matrices for certain energies. It may come as a surprise that dynamical bounds can be obtained if there is only one energy where one can exhibit a power-law bound for the transfer matrix. This is indeed necessary for models such as the random dimer model and related ones , where there are only a finite number of such energies.

Another advantage of the approach from  over the bound (4) is the stability of its input with respect to perturbations of the potential $`V`$. It was noted in  that if its approach can be applied to a given model, then it can also be applied to all finitely supported perturbations of the given potential—and it gives the same dynamical bounds for the perturbed models. Such a stability is not true in general for bounds derived using (4). For example, it may happen that the addition of a finitely supported perturbation turns a given singular continuous spectral measure into a pure point measure; see  for many examples illustrating this phenomenon.

In , the general criterion was applied to three prominent models from one-dimensional quasicrystal theory, namely, the Fibonacci model, the period doubling model, and the Thue-Morse model. All these models can be generated by a substitution process. This allows one to study the growth of transfer matrix norms with the help of an associated dynamical system—the trace map—and this provides in particular a very convenient way of verifying the input to the general dynamical criterion.

In the present paper we will prove a more general version of the dynamical result from , involving also the weight assigned by the spectral measure to the set of energies with power-law bounded transfer matrices. This gives stronger dynamical results in cases where such bounds hold for all energies in the spectrum, for example, models with Sturmian potentials. We shall also prove a stronger stability result. Namely, we will show that, for a fixed energy, the power-law bound is stable with respect to power-decaying perturbations. Here, the power-decay of the perturbation that we can allow depends on the transfer matrix power-law bound we start out with. Finally, we shall study a large number of examples and derive dynamical results for them by applying our main theorem, Theorem 1 below. The examples discussed in this paper include, in particular, generalizations of each of the three prominent substitution models studied in .

The organization of the paper is as follows. In Section 2 we prove our main theorem which derives quantum dynamical lower bounds from power-law bounds on transfer matrices. Section 3 discusses the stability of such power-law bounds on transfer matrices with respect to power-decaying perturbations of the potential. Section 4 deals with a class of models that are “sparse” in a certain sense and which includes a variety of substitution models (in particular, generalizations of Fibonacci, period doubling, and Thue-Morse), the prime Schrödinger operator, and moderately sparse models which were studied by Zlatoš in . The hierarchical model, which was studied in detail by Kunz et al. in  from a spectral point of view, will then be considered in Section 5. Finally, we present results for Sturmian models (studied, e.g., in ; see also the reviews ) in Section 6.

## 2. A Quantum Dynamical Lower Bound Derived From Power-Law Transfer Matrix Bounds

In this section, we prove a more general version of the main result from . The general idea of proof is the same and the result derives lower bounds on the dynamical quantity $`\beta _{\delta _1}^{}(p)`$ from power-law bounds on transfer matrices. However, the result established in this section gives improved bounds in many cases, in particular in the case of Sturmian potentials discussed later in the paper.

Recall the notion of a transfer matrix. Consider for some $`E`$, a solution $`\varphi `$ of the difference equation

(5)
$$\varphi (n+1)+\varphi (n1)+V(n)\varphi (n)=E\varphi (n).$$

Denote $`\mathrm{\Phi }(n)=(\varphi (n+1),\varphi (n))^T`$. The transfer matrix $`T(n,m;E)`$ is defined by requiring

$$\mathrm{\Phi }(n)=T(n,m;E)\mathrm{\Phi }(m)$$

for every solution $`\varphi `$ of (5). It is straightforward to verify that for $`n>m`$,

$$T(n,m;E)=T(V(n);E)\times \mathrm{}\times T(V(m+1);E),$$

where

$$T(x;E)=\left(\begin{array}{cc}Ex& \hfill 1\\ 1& \hfill 0\end{array}\right),$$

and similarly for $`n<m`$.

With this notation at hand we can now state:

###### Theorem 1.

The following statements hold:
(a) Suppose that for some $`K>0`$, $`C>0`$, $`\alpha >0`$, the following condition holds: For any $`N>0`$ large enough, there exists a nonempty Borel set $`A(N)`$ such that $`A(N)[K,K]`$ and

(6)
$$T(n,m;E)CN^\alpha EA(N),n,m:|n|N,|m|N$$

(resp., with $`1nN,1mN`$ in the case of $`\mathrm{}^2()`$). Let $`N(T)=T^{1/(1+\alpha )}`$ and let, for $`j=1,2`$, $`B_j(T)`$ be the $`j/T`$-neighborhood of the set $`A(N(T))`$:

$$B_j(T)=\{E:E^{}A(N(T)),|EE^{}|j/T\}.$$

Denote by $`F(z)`$ the Borel transform of the spectral measure of the state $`\psi =\delta _1`$:

$$F(E+i\epsilon )=_{}\frac{d\mu (x)}{x(E+i\epsilon )}.$$

Then for the initial state $`\psi =\delta _1`$ and all $`T>1`$ large enough, the following bound holds:

(7)
$$P(T)\underset{n:|n|N(T)}{}a(n,T)\frac{C}{T}N^{12\alpha }(T)_{B_2(T)}𝑑E(1+\mathrm{Im}^2F(E+i\epsilon )).$$

In particular,

(8)
$$P(T)\frac{C}{T}N^{12\alpha }(T)(|B_1(T)|+\mu (B_1(T))),$$

where $`|B|`$ denotes the Lebesgue measure. This gives the following bound for the time-averaged moments:

(9)
$$|X|_{\delta _1}^p(T)\frac{C}{T}N^{p+12\alpha }(T)(|B_1(T)|+\mu (B_1(T))).$$

(b) Suppose that there exists a set $`A[K,K]`$ of positive measure $`\mu (A)>0`$ such that

$$T(n,m;E)C(|n|^\alpha +|m|^\alpha )$$

for all $`EA,n,m`$. Then

(10)
$$\beta _{\delta _1}^{}(p)\frac{p3\alpha }{1+\alpha }.$$

(c) Assume that

$$T(n,m;E_0)C(E_0)(|n|^\alpha +|m|^\alpha )$$

for some $`E_0`$, uniformly in $`n,m`$, then

(11)
$$|X|_{\delta _1}^p(T)CT^{\frac{p3\alpha }{1+\alpha }}\left(T^1+\mu \left([E_0T^1,E_0+T^1]\right)\right).$$

Assume moreover that $`E_0`$ is an eigenvalue (possible only if $`\alpha >1/2`$), so that there exists $`\psi \mathrm{}^2,\psi 0`$ such that $`H\psi =E_0\psi `$. Suppose that $`\psi (1)0`$ (this is always true in the case of $`\mathrm{}^2()`$). Then

(12)
$$\beta _{\delta _1}^{}(p)\frac{p+12\alpha }{1+\alpha }.$$

###### Proof.

As in  we shall consider the case of $`\mathrm{}^2()`$, because for $`\mathrm{}^2()`$, the proof is similar but simpler. The main part of the proof is virtually identical with that of . For the sake of completeness we shall briefly recall the main lines.

The starting point is the Parseval equality:

$$a(n,T)\frac{1}{T}_0^{\mathrm{}}e^{2t/T}|\delta _n,\mathrm{exp}(itH)\delta _1|^2𝑑t=\frac{\epsilon }{2\pi }_{}|\delta _n,R(E+i\epsilon )\delta _1|^2𝑑E,$$

where $`R(z)=(HzI)^1`$ and $`\epsilon =1/T`$. For $`z=E+i\epsilon ,\epsilon >0`$, we define $`\varphi =R(z)\delta _1`$, $`\mathrm{\Phi }(n)=(\varphi (n+1),\varphi (n))^T`$. For each $`n>1`$, one has the inequality

(13)
$$\mathrm{\Phi }(n)T(n,1;z)^1\mathrm{\Phi }(1)$$

and for each $`n<0`$,

(14)
$$\mathrm{\Phi }(n)T(n,0;z)^1\mathrm{\Phi }(0).$$

An upper bound for the norm of the transfer matrix with complex $`z`$ is obtained using condition (6) and \[8, Lemma 2.1\]. Namely, let us fix some $`T>1,\epsilon =1/T`$ and define $`NN(T)=T^{1/(1+\alpha )}`$. Then for every $`EB_2(T)`$ and $`1nN`$

(15)
$$T(n,1;E+i\epsilon )DN^\alpha ,$$

where $`D=C\mathrm{exp}(3C)`$, and $`C`$ is the constant from (6). A similar bound holds for negative values of $`n`$. Using the bounds (13)–(15), one shows that for every $`EB_2(T)`$,

(16)
$$\underset{n:|n|N/2}{}|\delta _n,R(E+i\epsilon )\delta _1|^2cN^{12\alpha }(|\varphi (0)|^2+|\varphi (1)|^2+|\varphi (2)|^2)$$

with uniform constant $`c>0`$. It was shown in  that under the conditions of the theorem one always has

$$|\varphi (0)|+|\varphi (1)|+|\varphi (2)|c>0$$

with uniform constant. What one can also observe (and this is a new point) is the fact that

$$\varphi (1)=R(z)\delta _1,\delta _1=F(z),$$

where $`F(z)`$ is the Borel transform of the spectral measure corresponding to the pair $`(H,\delta _1)`$. Therefore, it follows from (16) that

$$\underset{n:|n|N/2}{}|\delta _n,R(E+i\epsilon )\delta _1|^2cN^{12\alpha }(1+\mathrm{Im}^2F(E+i\epsilon )).$$

Integrating this bound over $`EB_2(T)`$, one proves (7). Next, one observes that $`1+\mathrm{Im}^2F(z)2\mathrm{I}\mathrm{m}F(z)`$. For any set $`S`$, denote by $`S_\epsilon `$ the $`\epsilon `$-neighborhood of $`S`$. Following , one can see that

$`{\displaystyle _{S_\epsilon }}\mathrm{Im}F(E+i\epsilon )𝑑E`$ $`={\displaystyle _{}}𝑑\mu (x){\displaystyle _{S_\epsilon }}{\displaystyle \frac{\epsilon dE}{(xE)^2+\epsilon ^2}}`$
$`{\displaystyle _S}𝑑\mu (x){\displaystyle _\epsilon ^\epsilon }{\displaystyle \frac{\epsilon du}{u^2+\epsilon ^2}}`$
$`={\displaystyle \frac{\pi }{2}}\mu (S).`$
Taking $`S=B_1(T)`$, we prove (8). The bound (9) immediately follows.

To prove part (b), one just takes $`A(N)=A`$ for every $`N`$. Since $`\mu (B_1(T))\mu (A(N(T))=\mu (A)>0`$, the result follows from the bound (9).

The bound (11) of part (c) follows directly from (9), taking $`A(N)=\{E_0\}`$ for every $`N`$. Finally, to prove the second part of (c), we go back to (7) to obtain

$$|X|_{\delta _1}^p(T)\frac{C}{T}N^{p+12\alpha }(T)_{B_2(T)}\mathrm{Im}^2F(E+i\epsilon )𝑑E,$$

where $`B_2(T)=[E_02\epsilon ,E_0+2\epsilon ]`$. Under the condition $`\psi (1)0`$, one has $`\mu (\{E_0\})>0`$. Thus,

$$\mathrm{Im}F(E+i\epsilon )\frac{c\epsilon }{(EE_0)^2+\epsilon ^2}.$$

Integration over $`B_2(T)`$ yields (12). ∎

Remark. Part (b) of Theorem 1 remains true if

(17)
$$T(n,m;E)C(E)(|n|^\alpha +|m|^\alpha )$$

for all $`n,m`$ and $`EA`$ with $`C(E)<\mathrm{}`$ for $`\mu `$-almost every $`E`$. To prove this, it is sufficient to take a smaller set $`A^{}A`$ of positive measure where $`C(E)C<\mathrm{}`$. The bound (10) should be compared with the well-known result of : If (17) holds for some $`\alpha [0,1/2)`$ on a set $`A`$ of positive $`\mu `$-measure, then the restriction of $`\mu `$ to $`A`$ is $`12\alpha `$-continuous. In particular,

$$\beta _{\delta _1}^{}(p)p(12\alpha ).$$

This bound is better than (10) for small $`p`$, but for $`p`$ large enough, (10) is always better. Moreover, (10) holds also if $`\alpha 1/2`$.

## 3. Stability With Respect to Power-Decaying Perturbations

In this section, we discuss the stability of the crucial input to our dynamical bounds, power-law bounds on transfer matrices, with respect to perturbations of the potential. It is easy to see, and was noted in \[8, Corollary 1.3\], that finitely supported perturbations of the potential cannot destroy such a power-law bound. Here we strengthen this to stability with respect to power-decaying perturbations, where the allowed power depends on the bound we can prove for the unperturbed problem.

###### Theorem 2.

Assume that for some energy $`E`$ and some constant $`C_1`$, the transfer matrices $`T`$ associated with $`\mathrm{\Delta }+V`$ satisfy

(18)
$$T(n,m;E)C_1|nm|^\alpha \text{ for every }n,m\text{ with }nm0.$$

Assume further that, for some $`\epsilon >0`$, the perturbation $`W`$ satisfies

(19)
$$|W(n)|C_2(1+|n|)^{12\alpha \epsilon }\text{ for every }n.$$

Then the transfer matrices $`T^{}`$ associated with $`\mathrm{\Delta }+V+W`$ satisfy

(20)
$$T^{}(n,m;E)C_3|nm|^\alpha \text{ for every }n,m\text{ with }nm0.$$

###### Proof.

We present the proof in the special case where we assume (18) only for $`n0`$ and $`m=0`$ and then prove (20) for $`n0`$ and $`m=0`$. A slight variation of the argument below works for general $`n,m`$ with $`nm0`$ (with a uniform constant $`C_3`$ in (20)).

Our strategy will be to work with solutions and employ a general perturbation method developed by Kiselev et al. in .

Consider the unperturbed equation (5) and the perturbed equation

(21)
$$\psi (n+1)+\psi (n1)+[V(n)+W(n)]\psi (n)=E\psi (n).$$

Note that the transfer matrix $`T^{}(n,0;E)`$ is given by

$$T^{}(n,0;E)=\left(\begin{array}{cc}\psi _\mathrm{D}(n+1)& \psi _\mathrm{N}(n+1)\\ \psi _\mathrm{D}(n)& \psi _\mathrm{N}(n)\end{array}\right),$$

where $`\psi _\mathrm{D},\psi _\mathrm{N}`$ solve (21) and obey

$$\left(\begin{array}{cc}\psi _\mathrm{D}(1)& \psi _\mathrm{N}(1)\\ \psi _\mathrm{D}(0)& \psi _\mathrm{N}(0)\end{array}\right)=I.$$

Fix a complex reference solution $`\varphi `$ of (5). For example, we could set $`\varphi =\varphi _\mathrm{D}+i\varphi _\mathrm{N}`$, where $`\varphi _\mathrm{D},\varphi _\mathrm{N}`$ solve (5) and have the same initial conditions as $`\psi _\mathrm{D},\psi _\mathrm{N}`$. By (18) we have

(22)
$$|\varphi (n)|C|n|^\alpha .$$

Let $`\psi `$ be one of the basic solutions $`\psi _\mathrm{D},\psi _\mathrm{N}`$ of (21). Define $`\rho (n)`$ by

$`\left(\begin{array}{c}\psi (n)\\ \psi (n1)\end{array}\right)`$ $`={\displaystyle \frac{1}{2i}}\left[\rho (n)\left(\begin{array}{c}\varphi (n)\\ \varphi (n1)\end{array}\right)\overline{\rho (n)}\left(\begin{array}{c}\overline{\varphi (n)}\\ \overline{\varphi (n1)}\end{array}\right)\right]`$
$`=\mathrm{Im}\left[\rho (n)\left(\begin{array}{c}\varphi (n)\\ \varphi (n1)\end{array}\right)\right].`$
Write $`\varphi (n)`$ and $`\rho (n)`$ in polar coordinates,

$$\varphi (n)=|\varphi (n)|e^{i\gamma (n)}\rho (n)=R(n)e^{i\eta (n)},$$

and define

$$\theta (n)=\eta (n)+\gamma (n)\text{ and }U(n)=\frac{2W(n)}{\omega }|\varphi (n)|^2,$$

where $`i\omega `$ is the Wronskian of $`\overline{\varphi }`$ and $`\varphi `$, that is,

$$2i\mathrm{Im}(\varphi (n+1)\overline{\varphi (n)})=i\omega \text{ for every }n.$$

Clearly, the assertion of the theorem follows if we can show that $`R(n)`$ remains bounded as $`|n|\mathrm{}`$. The key identity (equation (45) in ) is the following:

(23)
$$R(n+1)^2=R(n)^2[1+U(n)\mathrm{sin}(2\theta (n))+U(n)^2\mathrm{sin}^2(\theta (n))].$$

It follows from (19) and (22) that $`U(n)`$ is summable. Thus, boundedness of $`R(n)`$ follows from this and (23) (cf., e.g., \[20, Lemma 3.5\]). This concludes the proof. ∎

The theorem above implies the stability of the number $`\alpha `$ and of the sets $`A(N)`$, $`B_1(T)`$, $`A`$ under suitable power-decaying perturbations of the potential. On the other hand, the measure of the sets $`\mu (B_1(T)),\mu (A)`$ and the Borel transform $`F(z)`$ may change after such a perturbation. In particular, it is possible that $`\mu (A)=0`$ for the perturbed operator in part (b) of Theorem 1. Thus, the bounds (10) and (12) are in general not stable. Of course, we still get a dynamical bound for the perturbed model. For example, we have the following consequence of Theorem 1 and Theorem 2:

###### Corollary 3.1.

Assume that for some energy $`E_0`$ and some constant $`C_1`$, the transfer matrices $`T`$ associated with $`\mathrm{\Delta }+V`$ satisfy $`T(n,m;E_0)C_1|nm|^\alpha `$ for every $`n,m`$ with $`nm0`$. Assume further that, for some $`\epsilon >0`$, the perturbation $`W`$ satisfies $`|W(n)|C_2|n|^{12\alpha \epsilon }`$ for every $`n`$. Then we have for the operator $`\mathrm{\Delta }+V+W`$,

$$\beta _{\delta _1}^{}(p)\frac{p14\alpha }{1+\alpha }$$

for every $`p>0`$.

###### Proof.

By Theorem 2, we have that the transfer matrices $`T^{}`$ associated with $`\mathrm{\Delta }+V+W`$ satisfy $`T^{}(n,m;E_0)C|nm|^\alpha `$ for every $`n,m`$ with $`nm0`$. Then, an inspection of the proof of Theorem 1 shows that this suffices to prove the bound (11) which yields

$$<|X|_{\delta _1}^p>(T)CT^{\frac{p3\alpha }{1+\alpha }1}$$

and the assertion of the corollary follows. More precisely, one can work independently on the two half-lines and hence needs bounds on $`T^{}(n,m;E)`$ only for the case where $`n,m`$ have the same sign. ∎

## 4. A Class of Pseudo-Sparse Potentials

In this section, we study a class of “sparse” potentials which includes various substitution models and the prime model. These potentials are not all sparse in the standard sense, but the point is that the class we discuss contains sparse potentials, and also a number of other potentials that have been considered before and which can be studied within the same framework.

Let us consider the case where the potential $`V`$ is defined on the half-line $``$ and takes on two values $`a,b`$. We assume the following for $`n`$ large enough, that is, for $`nN`$:

* Occurrences of $`b`$ are always isolated, that is, if $`V(n)=b`$ for some $`n`$, then $`V(n1)=V(n+1)=a`$.
* The value $`a`$ always occurs with odd multiplicity, that is, if $`V(n)=V(n+k+1)=b`$ and $`V(n+j)=a`$, $`1jk`$, then $`k`$ is odd.

Sparseness in this context refers to the $`b`$’s being isolated and the results below holding for arbitrarily long gaps between consecutive $`b`$’s. However, some of the concrete applications—for example the applications to substitution models—will not be sparse in a traditional sense.

We can prove the following:

###### Theorem 3.

Suppose $`V:\{a,b\}`$ is a potential satisfying (S1) and (S2) above. We have for every $`p>0`$,

$$\beta _{\delta _1}^{}(p)\frac{p5}{2}.$$

###### Proof.

Up to an initial piece, the transfer matrices are given by products of matrices of the following form:

$$T(a,E)^{2l+1}\text{ and }T(b,E).$$

Let $`E_0=a`$. Then

$$T(a,E_0)^{2l+1}=\left(T(a,E_0)^2\right)^lT(a,E_0)=\left(I\right)^lT(a,E_0)=\pm T(a,E_0).$$

Up to sign, this gives rise to powers of

$$T(a,E_0)T(b,E_0)=\left(\begin{array}{cc}\hfill 0& \hfill 1\\ \hfill 1& \hfill 0\end{array}\right)\left(\begin{array}{cc}ab& \hfill 1\\ 1& \hfill 0\end{array}\right)=\left(\begin{array}{cc}1& \hfill 0\\ ab& \hfill 1\end{array}\right).$$

Clearly, such powers satisfy a bound which is linear in the number of factors. Thus, the claim follows from (11). ∎

Remark. We can apply Corollary 3.1 and obtain that the dynamical bound in Theorem 3 is stable with respect to perturbations $`W`$ obeying $`|W(n)|C_2n^{3\epsilon }`$ for some fixed $`\epsilon >0`$ and every $`n`$. Similarly, we have stability with respect to power-decaying perturbations for all the dynamical bounds that will be shown in this section and we will not make this explicit for each one of them.

Let us now discuss the case where the $`a`$’s occur with even multiplicities. That is, we assume for $`n`$ large enough,

* The value $`a`$ always occurs with even multiplicity, that is, if $`V(n)=V(n+k+1)=b`$ and $`V(n+j)=a`$, $`1jk`$, then $`k`$ is even.

In this case we can prove a dynamical bound even without assuming the sparseness condition (S1). However, we need that $`|ab|`$ is not too large. Namely, we have the following result:

###### Theorem 4.

Suppose $`V:\{a,b\}`$ is a potential satisfying (S3) above.
(a) If $`|ab|<2`$, then for every $`p>0`$,

$$\beta _{\delta _1}^{}(p)p1.$$

(b) If $`|ab|=2`$, then for every $`p>0`$,

$$\beta _{\delta _1}^{}(p)\frac{p5}{2}.$$

###### Proof.

The argument proceeds in a way similar to the proof above. Again, up to an initial piece, the transfer matrices are given by products of matrices of the following form:

$$T(a,E)^{2l}\text{ and }T(b,E).$$

Again, let $`E_0=a`$. Then

$$T(a,E_0)=\left(\begin{array}{cc}\hfill 0& \hfill 1\\ \hfill 1& \hfill 0\end{array}\right)$$

and hence

$$T(a,E_0)^{2l}=\left(T(a,E_0)^2\right)^l=\left(I\right)^l=\pm I.$$

On the other hand, $`T(b,E_0)`$ is elliptic when $`|ab|<2`$ and parabolic when $`|ab|=2`$. Thus, in the former case, products of matrices of the form $`T(a,E)^{2l}`$ or $`T(b,E)`$ remain bounded, while in the latter case such products satisfy a bound which is linear in the number of factors. The claim thus follows from (11). ∎

Let us note that a result like part (a) of Theorem 4 is implicitly contained in , where mainly random polymer models are studied.

It is clear that whole-line analogs of the above theorems hold. In this case, we need (S1) and (S2) or (S3) to hold for $`|n|`$ large enough.

More importantly, these results cover a variety of seemingly very different cases: First consider the period doubling Hamiltonian, which was already discussed in . On the alphabet $`A=\{a,b\}`$, consider the period doubling substitution $`S(a)=ab`$, $`S(b)=aa`$. Iterating on $`a`$, we obtain a one-sided sequence

$$u=abaaabababaaabaaab\mathrm{}$$

which is invariant under the substitution process. Define the associated subshift $`\mathrm{\Omega }_{\mathrm{pd}}`$ to be the set of all sequences over $`A`$ which have all their finite subwords occurring in $`u`$. Here, we can consider either one-sided or two-sided sequences. This does not matter for the results in this paper, but we remark that for substitution models, one generally considers the two-sided case. For $`\omega \mathrm{\Omega }_{\mathrm{pd}}`$, we define the potential $`V_\omega `$ by $`V_\omega (n)=\omega _n`$. It is easy to check that each $`V_\omega `$ satisfies (S1) and (S2) (even for every $`n`$) and hence an application of Theorem 3 allows us to recover \[8, Theorem 3\]. However, we can prove a more general result. Consider, for example, substitutions of the form

(24)
$$S(a)=a^{2k1}b,S(b)=a^{2l},k,l1.$$

The case $`k=1,l=1`$ corresponds to the period doubling case. The potentials generated by a substitution of the form (24) (by generating a one-sided fixed point and passing to the associated subshift, as in the period doubling case above) are easily seen to obey (S1) and (S2). On the other hand, substitutions of the form

(25)
$$S(a)=a^{2k}b,S(b)=a^{2l},k,l1$$

give rise to potentials satisfying (S3) and hence Theorem 4 applies in these cases. Thus we may state the following:

###### Corollary 4.1.

(a) Let $`S`$ be a substitution of the form (24), $`\mathrm{\Omega }`$ the associated subshift, and for $`\omega \mathrm{\Omega }`$, let $`V_\omega (n)=\omega _n`$, $`n`$. Then, for every $`\omega \mathrm{\Omega }`$, the potential $`V_\omega `$ gives rise to an operator satisfying

$$\beta _{\delta _1}^{}(p)\frac{p5}{2}\text{ for every }p>0.$$

(b) Let $`S`$ be a substitution of the form (25), $`\mathrm{\Omega }`$ the associated subshift, and for $`\omega \mathrm{\Omega }`$, let $`V_\omega (n)=\omega _n`$, $`n`$. Then, for every $`\omega \mathrm{\Omega }`$, the potential $`V_\omega `$ gives rise to an operator satisfying

$$\beta _{\delta _1}^{}(p)p1\text{ for every }p>0\text{ if }|ab|<2$$

and

$$\beta _{\delta _1}^{}(p)\frac{p5}{2}\text{ for every }p>0\text{ if }|ab|=2.$$

Consider the following class of substitutions:

(26)
$$S(a)=a^mb^n,S(b)=a.$$

The case $`m=n=1`$ gives rise to the Fibonacci substitution. Hence, the substitutions in (26) are usually called generalized Fibonacci substitutions. If $`n=1`$, the resulting potentials are Sturmian and will be discussed in this more general context in a later section. Here, we restrict our attention to the case $`n2`$. These substitutions and the associated Schrödinger operators were studied, for example, in .

If $`n`$ is even, it is easily seen that each $`V_\omega `$ satisfies (S3) with the roles of $`a`$ and $`b`$ interchanged, that is, $`b`$’s always occur with even multiplicity. Thus, we can derive a dynamical bound for the associated operators by applying Theorem 4.

If $`n`$ is odd, the model satisfies neither (S2) nor (S3) but we can nevertheless employ a similar argument. As a warmup, let us consider the case $`n=3`$ (the special case $`m=1`$, $`n=3`$ is usually called the nickel mean substitution). Then the transfer matrices are given by products of matrices of the following form:

$$T(a,E)\text{ and }T(b,E)^3.$$

Let $`E_0=b+1`$. Then

$$T(b,E_0)=\left(\begin{array}{cc}\hfill 1& \hfill 1\\ \hfill 1& \hfill 0\end{array}\right)$$

and hence

$$T(b,E_0)^3=I.$$

This would allow us to prove bounds on $`\beta _{\delta _1}^{}(p)`$ in the same way as in the proof of Theorem 4.

Let us now turn to the case of a general odd $`n3`$. Here we can extend the above idea and prove a result which applies to the substitutions in (26) with $`n`$ odd but which is much more general. Denote

* There is some odd $`k3`$ such that the value $`b`$ always occurs with a multiplicity which is a multiple of $`k`$, that is, if $`V(n)=V(n+l+1)=a`$ and $`V(n+j)=b`$, $`1jl`$, then $`l=mk`$ for some $`m`$,

Then we can prove the following:

###### Theorem 5.

Suppose $`V:\{a,b\}`$ is a potential satisfying (S4). Then there is a set $``$ of cardinality $`k1`$ such that for every $`E`$, we have
(a) If $`|aE|<2`$, then for every $`p>0`$,

$$\beta _{\delta _1}^{}(p)p1.$$

(b) If $`|aE|=2`$, then for every $`p>0`$,

$$\beta _{\delta _1}^{}(p)\frac{p5}{2}.$$

###### Proof.

In this case, the transfer matrices are given by products of matrices of the following form:

$$T(a,E)\text{ and }T(b,E)^k.$$

It suffices to exhibit $`k1`$ energies $`E_0`$ with

(27)
$$T(b,E_0)^k=\pm I.$$

This can be seen as follows: The matrix $`T(b,E)^k`$ is the monodromy matrix of the constant potential $`V(n)=b`$, regarded as a $`k`$-periodic potential. This gives rise to an operator with $`k1`$ gaps. However, since the operator with this potential has spectrum $`[b2,b+2]`$, all these gaps are degenerate. Every degenerate gap corresponds to an energy where the monodromy matrix is equal to $`\pm I`$, hence there are exactly $`k1`$ energies $`E_0`$ for which we have (27). ∎

Putting everything together, we obtain the following result for the models generated by substitutions from (26):

###### Corollary 4.2.

Let $`S`$ be a substitution of the form (26), $`\mathrm{\Omega }`$ and the $`V_\omega `$’s as above.
(a) If $`n2`$ is even, then for every $`\omega \mathrm{\Omega }`$, the potential $`V_\omega `$ gives rise to an operator satisfying

$$\beta _{\delta _1}^{}(p)p1\text{ for every }p>0\text{ if }|ab|<2$$

and

$$\beta _{\delta _1}^{}(p)\frac{p5}{2}\text{ for every }p>0\text{ if }|ab|=2.$$

(b) If $`n3`$ odd, then $`T(b,E_0)^n=\pm I`$ has $`n1`$ solutions $`E_0`$ and for each such solution $`E_0`$, we have that for every $`\omega \mathrm{\Omega }`$, the potential $`V_\omega `$ gives rise to an operator satisfying

$$\beta _{\delta _1}^{}(p)p1\text{ for every }p>0\text{ if }|aE_0|<2$$

and

$$\beta _{\delta _1}^{}(p)\frac{p5}{2}\text{ for every }p>0\text{ if }|aE_0|=2.$$

The final substitution model we consider is the following:

(28)
$$S(a)=a^mb^n,S(b)=b^na^m.$$

The case $`m=n=1`$ gives rise to the Thue-Morse substitution. Hence, the substitutions in (28) are usually called generalized Thue-Morse substitutions. They were considered, for example, in . If at least one of $`m,n`$ is even, (S3) holds and we can apply Theorem 4. In the remaining case, where both $`m`$ and $`n`$ are odd (and at least one is $`3`$), (S4) holds and we can apply Theorem 5. Thus, for models generated by generalized Thue-Morse substitutions, we obtain the following dynamical bounds:

###### Corollary 4.3.

Let $`S`$ be a substitution of the form (28), $`\mathrm{\Omega }`$ and the $`V_\omega `$’s as above.
(a) If at least one of $`m,n`$ is even, then for every $`\omega \mathrm{\Omega }`$, the potential $`V_\omega `$ gives rise to an operator satisfying

$$\beta _{\delta _1}^{}(p)p1\text{ for every }p>0\text{ if }|ab|<2$$

and

$$\beta _{\delta _1}^{}(p)\frac{p5}{2}\text{ for every }p>0\text{ if }|ab|=2.$$

(b) If we have $`m3`$ odd, then $`T(b,E_0)^m=\pm I`$ has $`m1`$ solutions $`E_0`$ and for each such solution $`E_0`$, we have that for every $`\omega \mathrm{\Omega }`$, the potential $`V_\omega `$ gives rise to an operator satisfying

$$\beta _{\delta _1}^{}(p)p1\text{ for every }p>0\text{ if }|bE_0|<2$$

and

$$\beta _{\delta _1}^{}(p)\frac{p5}{2}\text{ for every }p>0\text{ if }|bE_0|=2.$$

An analogous result holds if we have $`n3`$ odd.
(c) If $`m=n=1`$, then

$$\beta _{\delta _1}^{}(p)p1\text{ for every }p>0.$$

Part (c) was proved in  and is stated for completeness. One might expect the bound $`\beta _{\delta _1}^{}(p)p1`$ to hold always. In fact, the paper  claims, for every choice of $`m,n,a,b`$, the existence of an energy, where the transfer matrices remain bounded. However, the argument given in that paper is incomplete and it would be interesting to prove or disprove this claim.

Next, we consider the prime Schrödinger operator $`H_{\mathrm{prime}}`$ on $`\mathrm{}^2()`$ whose potential is given by

$$V_{\mathrm{prime}}(n)=\{\begin{array}{cc}a& \text{ if }n\text{ is not prime,}\hfill \\ b& \text{ if }n\text{ is prime.}\hfill \end{array}$$

This operator was studied, for example, in . Based on numerics and heuristics contained in these two papers, one may expect the following: On the one hand, for almost every energy $`E`$, there is an $`\mathrm{}^2`$ solution to $`H_{\mathrm{prime}}\varphi =E\varphi `$, that is, when one varies the boundary condition at the origin, one gets pure point spectrum for almost every boundary condition. On the other hand, the model displays non-trivial transport for every boundary condition. We will confirm the latter below (the proof discusses only the case of a Dirichlet boundary condition, but it readily extends to every other boundary condition). Let us briefly discuss the first point. It is natural to view $`V_{\mathrm{prime}}`$ as a sparse potential. In fact, this point of view was proposed in . However, the current methods in the spectral analysis of models with sparse potentials (see, in particular, ) are clearly insufficient to conclude anything for the prime model. We regard this as an interesting problem and refer the reader also to  for further motivation to consider models of moderate sparseness.

Let us now turn to a dynamical result for the prime model. Clearly, (S1) and (S2) are satisfied for $`n`$ large enough. Hence we get:

###### Corollary 4.4.

For every $`a,b`$, the operator $`H_{\mathrm{prime}}`$ satisfies

$$\beta _{\delta _1}^{}(p)\frac{p5}{2}\text{ for every }p>0.$$

Finally, we discuss a model which is sparse in the standard sense. Namely, pick some integer $`\gamma 2`$ and define $`n_k=\gamma ^k`$ for $`k`$. Let $`V_{\mathrm{sparse}}(n)=b`$ if $`n=n_k`$ for some $`k`$ and $`V_{\mathrm{sparse}}(n)=a`$ otherwise. Schrödinger operators with potentials of this kind were studied in . Clearly, when $`\gamma `$ is even, all $`n_k`$’s are even, and when $`\gamma `$ is odd, all $`n_k`$’s are odd, so we have (S1) and (S2). Thus, Theorem 3 applies and we get

###### Corollary 4.5.

For every $`a,b`$ and $`\gamma \{1\}`$, the potential $`V_{\mathrm{sparse}}`$ gives rise to an operator satisfying

$$\beta _{\delta _1}^{}(p)\frac{p5}{2}\text{ for every }p>0.$$

This can be improved if $`\gamma |ab|`$:

###### Proposition 4.6.

Let

$$\nu =\frac{2\mathrm{log}\sqrt{2+(ab)^2}}{\mathrm{log}\gamma }.$$

Then the potential $`V_{\mathrm{sparse}}`$ gives rise to an operator satisfying

$$\beta _{\delta _1}^{}(p)\frac{p14\nu }{1+\nu }\text{ for every }p>0.$$

###### Proof.

Write $`C(a,b)=\sqrt{2+(ab)^2}`$. Then

$$T(a,E=a)^{2l+1}T(b,E=a)=\left(\begin{array}{cc}1& \hfill 0\\ ab& \hfill 1\end{array}\right)C(a,b).$$

For $`d_{n,m}=\mathrm{\#}\{mkn:V(k)=b\}`$, we have $`d_{n,m}\mathrm{log}|nm|/\mathrm{log}\gamma `$ and hence

$$T(n,m;E=a)C(a,b)^{d_{n,m}}C(a,b)^{\mathrm{log}|nm|/\mathrm{log}\gamma }=|nm|^{\mathrm{log}C(a,b)/\mathrm{log}\gamma }.$$

This yields the assertion. ∎

## 5. A Hierarchical Model

The hierarchical model is defined through the potential

(29)
$$V(n)=\lambda f(\mathrm{ord}n),$$

where $`f`$ is some real function and $`\mathrm{ord}n`$ is the number of factors 2 in the prime decomposition of $`n`$. The sequence (29) has some nice symmetries. Because $`\mathrm{ord}(n)=\mathrm{ord}n`$ for all $`n`$ and $`\mathrm{ord}(l2^m+k)=\mathrm{ord}k`$ for $`m1`$, all $`l`$ and $`|k|<2^m`$, analogous identities hold for $`V`$. In particular,

(30)
$$V(l2^m+k)=V(k)=V(k)=V(l^{}2^mk)$$

for any $`l`$ and $`l^{}`$, $`m1`$ and $`|k|<2^m`$. The Schrödinger operator with such a potential appeared first in the works  and  with the special choice

$$f(m)=\underset{k=0}{\overset{m1}{}}R^k,$$

where $`R`$ is a positive constant. The advantage of this choice is that in this case,

$$x_m=\mathrm{tr}M_m(E)\mathrm{tr}T(2^m,0;E)$$

satisfies an autonomous difference equation ,

(31)
$$x_{m+1}=x_m^22+Rx_m(x_mx_{m1}^2+2),m1.$$

The above recurrence and the symmetries (30) made it possible to obtain many rigorous results about the spectrum of the corresponding Schrödinger operator. A detailed mathematical study of this model was carried out by Kunz et al. . Among other things, it was shown that for every $`R>0`$, the spectrum is a Cantor set, and for $`R1`$, it is purely singular continuous. From the point of view of the present article, it is interesting that a countable infinite set of exceptional energies in the spectrum could be identified explicitly. The $`2^m`$ zeros $`E_{mk}`$, $`1k2^m`$, of $`x_m(E)`$ are simple and $`x_m=0`$ implies $`x_{m+1}=2`$ and $`x_{m+l}=2`$ for $`l>1`$; compare (31). From this it was possible to show that $`E_{mk}`$, for $`m0`$ and $`1k2^m`$, are lower (resp., upper) gap-edges in the spectrum of $`H`$ if $`\lambda >0`$ (resp., $`\lambda <0`$) and they are dense in the spectrum. For the corresponding gap-edge states, the following result was obtained (Proposition 15 in ).

###### Proposition 5.1.

Let $`x_m(E)=0`$ and let $`\psi `$ be a solution of $`H\psi =E\psi `$.

* If $`\psi (0)=0`$, then $`\psi (k+2^{m+1})=\psi (k)`$ for every integer $`k`$.
* If $`\psi (0)0`$, then $`\psi (2l2^m)=(1)^l\psi (0)`$ and asymptotically, as $`l\mathrm{}`$,

(32)
$$\psi ((2l+1)2^m)\psi (2^m)(1)^{l+1}\lambda _m\psi (0)f_R(l)$$

where

(33)
$$f_R(l)=\{\begin{array}{cc}\frac{2}{2R}l,\hfill & R<2\hfill \\ l\mathrm{log}_2l,\hfill & R=2\hfill \\ \left(\frac{2}{R}\right)^{\epsilon _l}\frac{R^2}{2(R1)(R2)}l^{\mathrm{log}_2R},\hfill & R>2.\hfill \end{array}$$

Here $`\lambda _m=\lambda R^mx_{m1}(E)\mathrm{}x_0(E)`$, $`\epsilon _l[0,1)`$ is the fractional part of $`\mathrm{log}_2l`$ and $``$ means equality in the leading order of $`l`$.

We use this proposition to prove the following theorem.

###### Theorem 6.

For every $`\lambda 0`$ and $`R>0`$,

$$\beta _{\delta _1}^{}(p)\frac{p14\alpha }{1+\alpha },$$

where

$$\alpha =\alpha (R)=\mathrm{max}\{1,\mathrm{log}_2R\}.$$

###### Proof.

We apply Proposition 5.1 with $`m=0`$ for which it provides the precise asymptotic form of the solutions. Because $`x_0(E)=E`$, these belong to $`E=0`$. Let $`\psi _\mathrm{D}`$ and $`\psi _\mathrm{N}`$ be the two solutions defined by the initial values

(34)
$$\psi _\mathrm{D}(0)=\psi _\mathrm{N}(1)=0,\psi _\mathrm{D}(1)=\psi _\mathrm{N}(0)=1.$$

According to part (i) of Proposition 5.1, $`\psi _\mathrm{D}`$ is a periodic solution with period 4, namely

(35)
$$\psi _\mathrm{D}(2l)=0,\psi _\mathrm{D}(2l+1)=(1)^l.$$

On the other hand,

(36)
$$\psi _\mathrm{N}(2l)=(1)^l,\psi _\mathrm{N}(2l+1)(1)^{l+1}\lambda f_R(l).$$

Equations (35) and (36) permit us to compute the asymptotic form of $`T(n,m;0)`$. Because of $`V(n)=V(n)`$, it suffices to consider $`nm0`$. In what follows, we use the simplified notation $`T(n,m)`$. Let $`\mathrm{\Psi }^i(n)=(\psi ^i(n+1)\psi ^i(n))^T`$ for $`i=0,1`$. Then $`T(n,0)=(\mathrm{\Psi }_\mathrm{D}(n)\mathrm{\Psi }_\mathrm{N}(n))`$. The determinant of any transfer matrix being unity, the inverse is easy to compute. We find

(37) $`T(n,m)`$ $`=T(n,0)T(m,0)^1`$
(42) $`=\left(\begin{array}{cc}\psi _\mathrm{D}(n+1)\hfill & \psi _\mathrm{N}(n+1)\hfill \\ \psi _\mathrm{D}(n)\hfill & \psi _\mathrm{N}(n)\hfill \end{array}\right)\left(\begin{array}{cc}\hfill \psi _\mathrm{N}(m)& \hfill \psi _\mathrm{N}(m+1)\\ \hfill \psi _\mathrm{D}(m)& \hfill \psi _\mathrm{D}(m+1)\end{array}\right).`$
With the short-hand notation

$$F(l)=(1)^l\psi _\mathrm{N}(2l+1),$$

equations (35), (36), and (42) then yield

(45) $`T(2l,2k)`$ $`=(1)^{k+l}\left(\begin{array}{cc}1& F(l)F(k)\\ 0& 1\end{array}\right)`$
(48) $`T(2l+1,2k)`$ $`=(1)^{k+l}\left(\begin{array}{cc}0& 1\\ 1& F(l)F(k)\end{array}\right)`$
(51) $`T(2l,2k+1)`$ $`=(1)^{k+l+1}\left(\begin{array}{cc}F(l)F(k)& \hfill 1\\ 1& \hfill 0\end{array}\right)`$
(54) $`T(2l+1,2k+1)`$ $`=(1)^{k+l+1}\left(\begin{array}{cc}1& \hfill 0\\ F(l)F(k)& \hfill 1\end{array}\right).`$
All these matrices have the same norm. Denoting the Hilbert-Schmidt norm by $`_2`$, for $`n=2l,2l+1`$ and $`m=2k,2k+1`$, we have

$$T(n,m)T(n,m)_2=\sqrt{2+[F(l)F(k)]^2}\sqrt{2+\lambda ^2[f_R(l)f_R(k)]^2}.$$

Therefore,

$$T(n,m;0)2\lambda f_R(n/2)$$

for any $`n`$ large enough and $`mn`$. If $`R2`$, the assertion of the theorem obviously follows from the definition (33) of $`f_R`$ and Theorem 1. If $`R=2`$, we note that for any $`ϵ>0`$,

$$T(n,m;0)\lambda n^{1+ϵ}$$

if $`n`$ is large enough. Therefore, by Theorem 1,

$$\beta _{\delta _1}^{}(p)\frac{p54ϵ}{2+ϵ}$$

for any $`ϵ>0`$ and, thus, for $`ϵ=0`$ as well. ∎

Remark. The proof shows that we can apply Corollary 3.1 and obtain that the dynamical bound in Theorem 6 is stable with respect to perturbations $`W`$ obeying $`|W(n)|C_2|n|^{12\alpha \epsilon }`$ for some fixed $`\epsilon >0`$ and every $`n`$.

We note that instead of $`m=0`$, we could have used Proposition 5.1 with any $`m>0`$ and any zero of $`x_m(E)`$. This holds because of the following:

###### Theorem 7.

For any $`\lambda 0`$, $`R>0`$, $`m0`$, and $`k\{1,2,\mathrm{},2^m\}`$, there exists a positive number $`C_{\lambda ,R}(m,E_{mk})`$ such that for any $`nn^{}0`$,

$$T(n,n^{};E_{mk})C_{\lambda ,R}(m,E_{mk})f_R(2^{m1}n).$$

###### Proof.

We fix $`m>0`$ and a zero $`E_{mk}`$ of $`x_m`$. From equation (42) it is clear that we have to bound the two particular solutions (34) of $`H\psi =E_{mk}\psi `$. According to Proposition 5.1, $`\psi _\mathrm{D}`$ is $`2^{m+1}`$-antiperiodic and, thus, bounded. On the other hand,

(55)
$$\psi _\mathrm{N}(2l2^m)=(1)^l,\psi _\mathrm{N}((2l+1)2^m)\psi _\mathrm{N}(2^m)(1)^{l+1}\lambda _mf_R(l).$$

Thus, the task is to bound $`\psi _\mathrm{N}(n)`$ in the intervals

(56)
$$2l2^m<n<(2l+1)2^m\mathrm{and}(2l+1)2^m<n<2(l+1)2^m.$$

To proceed with the proof, let us recall equation (3.29) of , according to which

$$\psi _\mathrm{D}(2^m)=x_{m1}\mathrm{}x_0$$

for any energy. Thus, $`\psi _\mathrm{D}(2^m)0`$ in the present case ($`E=E_{mk}`$), for otherwise $`x_i=0`$ for some $`i<m`$ would imply $`|x_j|=2`$ for every $`j>i`$, contradicting $`x_m=0`$. Then $`u_0:=\psi _\mathrm{D}/\psi _\mathrm{D}(2^m)`$ is a solution of the Schrödinger equation satisfying the boundary conditions $`u_0(0)=0`$, $`u_0(2^m)=1`$ and, according to Proposition 5.1, $`u_0(k+2^{m+1})=u_0(k)`$ for any $`k`$. From the general theory of second-order difference (differential) equations, it follows that there exists a linearly independent solution $`u_1`$ with boundary values $`u_1(1)=1`$, $`u_1(2^m)=0`$ and that we can write $`\psi _\mathrm{N}`$ for $`0n2^m`$ in the form

$$\psi _\mathrm{N}(n)=\psi _\mathrm{N}(2^m)u_0(n)+\psi _\mathrm{N}(0)u_1(n).$$

Next, we observe that $`u_1`$ can be expressed in terms of $`u_0`$. Indeed, from equation (30) we can see that the sequence $`V(1),\mathrm{},V(2^m1)`$ is a palindrome,

$$V(2^{m1}k)=V(2^{m1}+k)k=1,\mathrm{},2^{m1}1$$

and, hence,

$$u_1(n)=u_0(2^mn),n=1,\mathrm{},2^m1.$$

Furthermore, the translational symmetry of the potential,

$$(V(l2^m+1),\mathrm{},V((l+1)2^m1))=(V(1),\mathrm{},V(2^m1)),$$

valid for any $`l`$, implies that the translates of $`u_0`$ and $`u_1`$ can be used to give $`\psi _\mathrm{N}`$ in each of the intervals (56). Altogether we find

$$\psi _\mathrm{N}(n)=\psi _\mathrm{N}((2l+1)2^m)u_0(n2l2^m)+\psi _\mathrm{N}(2l2^m)u_0((2l+1)2^mn)$$

if $`2l2^mn(2l+1)2^m`$ and

$$\psi _\mathrm{N}(n)=\psi _\mathrm{N}(2(l+1)2^m)u_0(n(2l+1)2^m)+\psi _\mathrm{N}((2l+1)2^m)u_0(2(l+1)2^mn)$$

if $`(2l+1)2^mn2(l+1)2^m`$. Together with (55), in both intervals,

$$|\psi _\mathrm{N}(n)|\frac{\mathrm{max}|\psi _\mathrm{D}|}{|\psi _\mathrm{D}(2^m)|}(|\psi _\mathrm{N}((2l+1)2^m)|+1).$$

Since $`ln/2^{m+1}`$, we obtain that for $`n`$ large enough

$$|\psi _\mathrm{N}(n)|\frac{\mathrm{max}|\psi _\mathrm{D}|}{|\psi _\mathrm{D}(2^m)|}(|\lambda _m|f_R(2^{m1}n)+|\psi _\mathrm{N}(2^m)|+1).$$

Due to (42), the assertion of the theorem follows from this bound. ∎

## 6. Sturmian Potentials

In this section, we discuss dynamical bounds for the standard one-dimensional quasicrystal model which is given by a Schrödinger operator on the whole line whose potential is given by

(57)
$$V(n)=\lambda v_{\omega ,\theta }(n),\text{ where }v_\theta (n)=\chi _{[1\omega ,1)}(n\omega +\theta mod1),$$

where $`\lambda 0`$ is the coupling constant, $`\omega (0,1)`$ irrational is the rotation number, and $`\theta [0,1)`$ arbitrary is the phase. For more information on this family of operators, we refer the reader to the survey articles .

It is well known, and easy to see, that the spectrum of the operator $`H_{\lambda ,\omega ,\theta }`$ with potential $`V`$ from (57) is independent of $`\theta `$, that is, for every $`\lambda ,\omega `$, there is a set $`\mathrm{\Sigma }_{\lambda ,\omega }`$ with $`\sigma (H_{\lambda ,\omega ,\theta })=\mathrm{\Sigma }_{\lambda ,\omega }`$ for every $`\theta `$.

Consider the continued fraction expansion of $`\omega `$,

$$\omega =\frac{1}{a_1+{\displaystyle \frac{1}{a_2+{\displaystyle \frac{1}{a_3+\mathrm{}}}}}}$$

with uniquely determined $`a_n`$ (cf. ). The associated rational approximants $`p_k/q_k`$ are defined by

$`p_0`$ $`=0,`$ $`p_1`$ $`=1,`$ $`p_k`$ $`=a_kp_{k1}+p_{k2},`$
$`q_0`$ $`=1,`$ $`q_1`$ $`=a_1,`$ $`q_k`$ $`=a_kq_{k1}+q_{k2}.`$
The number $`\omega `$ is said to have bounded density if

(58)
$$d(\omega )=\underset{n\mathrm{}}{lim\; sup}\frac{1}{n}\underset{k=1}{\overset{n}{}}a_k<\mathrm{}.$$

The set of bounded density numbers is uncountable but has Lebesgue measure zero. The following was shown in  (see also  for the case of zero phase):

###### Theorem 8.

Suppose $`\omega `$ is a bounded density number. For every $`\lambda `$, there is a constant $`C`$ such that for every $`\theta `$, every $`E\mathrm{\Sigma }_{\lambda ,\omega }`$, and every $`n,m`$, we have

(59)
$$T_{\lambda ,\omega ,\theta }(n,m;E)C|nm|^{\alpha (\lambda ,\omega )},$$

with

(60)
$$\alpha (\lambda ,\omega )=Dd(\omega )\mathrm{log}C_\lambda ,$$

where $`D`$ is some universal constant, $`C_\lambda `$ is given by

(61)
$$C_\lambda =2+\sqrt{8+\lambda ^2},$$

and $`d(\omega )`$ is as in (58).

This yields the following:

###### Corollary 6.1.

Let $`\omega `$ be a bounded density number. Then, for every $`\lambda ,\theta `$, the operator $`H_{\lambda ,\omega ,\theta }`$ satisfies

$$\beta _{\delta _1}^{}(p)\frac{p3\alpha (\lambda ,\omega )}{1+\alpha (\lambda ,\omega )}\text{ for every }p>0,$$

with $`\alpha (\lambda ,\omega )`$ given by (60).

Since $`\mu (\mathrm{\Sigma }_{\lambda ,\omega })=1`$, this is an immediate consequence of (10). This bound is better than the corresponding result in  (which follows from (9), bounding from below $`|B_1(T)|`$). One should stress that as opposed to all the other examples discussed earlier, the dynamical bound in Corollary 6.1 is not stable with respect to perturbations of the potential. This is due to the fact that $`\mu (\mathrm{\Sigma }_{\lambda ,\omega })`$ may vanish for the perturbed measure. However, by Corollary 3.1, we have the following result:

###### Corollary 6.2.

Let $`\omega `$ be a bounded density number and let $`\lambda `$ be arbitrary. If $`\alpha (\lambda ,\omega )`$ is given by (60) and $`W`$ satisfies

$$|W(n)|C_2(1+|n|)^{12\alpha (\lambda ,\omega )\epsilon }\text{ for every }n$$

for some $`\epsilon >0`$, then, for every $`\theta `$, the operator $`H_{\lambda ,\omega ,\theta }+W`$ satisfies

$$\beta _{\delta _1}^{}(p)\frac{p14\alpha (\lambda ,\omega )}{1+\alpha (\lambda ,\omega )}\text{ for every }p>0.$$

As in the case $`\omega =(\sqrt{5}1)/2`$ and $`\theta =0`$, studied in , it is possible to improve this lower bound somewhat by exhibiting a suitable set $`A(N)`$ (stable under perturbation), studying its Lebesgue measure, and applying (9). The set $`A(N)`$ will again be given by the spectra of suitable periodic approximants, and the Lebesgue measure can again be bounded through a fine analysis of the trace map, akin to what is done in ; compare also . We leave the details to the interested reader.