text:
# Spin Effects in the Local Density of States of GaAs

## Abstract

We present spin-resolved measurements of the local density of states in Si doped GaAs. Both spin components exhibit strong mesoscopic fluctuations. In the magnetic quantum limit, the main features of the spin-up and spin-down components of the local density of states are found to be identical apart from Zeeman splitting. Based on this observation, we introduce a mesoscopic method to measure the $`g`$-factor in a material where macroscopic methods are severely restricted by disorder. Differences between the spin-up and spin-down components are discussed in terms of spin relaxation due to spin-orbit coupling.

Spin-polarized electronic transport is currently attracting a lot of interest from both a fundamental and an applied point of view . Spin-polarized transport is essential for the operation of spin transistors or spin valves, and recently the spin injection from a ferromagnet into a semiconductor could be observed . Spin effects have also been proposed to be utilized for quantum computation. In this respect it is important to study spin relaxation and coherence in semiconductors.

Particularly interesting spin phenomena were observed in mesoscopic semiconductor structures. Prominent examples are the Kondo effect in quantum dots and spin-polarized tunneling through impurity levels. Resonant tunneling through impurities is also an established technique to image the local density of states (LDOS) of doped GaAs . The LDOS exhibits mesoscopic fluctuations which can be understood in terms of interference of elastically scattered electron waves. However, a spin-resolved measurement was impossible to date.

In this paper, we investigate the spin dependence of the LDOS via resonant tunneling through a spin-split impurity level. The strongest fluctuations of the spin-up and spin-down components of the LDOS are found to be identical apart from Zeeman splitting. This observation allows us to mesoscopically determine the $`g`$-factor in a material where macroscopic methods – such as photoluminescence or magnetotransport – are restricted by disorder. Differences between the spin-up and spin-down components are discussed in terms of spin relaxation due to spin-orbit coupling.

Our experiment is based on a strongly asymmetric double-barrier heterostructure, which was grown on $`n^+`$-type GaAs substrate. It consists of a 10 nm wide GaAs quantum well and two Al0.3Ga0.7As barriers of 5 and 8 nm width. The nominally undoped active region is sandwiched between 300 nm thick GaAs contact layers doped with Si. The donor concentration has been experimentally determined to $`3.3\times 10^{17}`$ cm-3 . From this material we fabricated a 2 $`\mu `$m diameter mesa with Ohmic contacts, which contains a small number of impurities in the quantum well.

We use the energetically-lowest impurity state $`S`$ as spectrometer for the LDOS $`\nu `$ in the emitter contact adjacent to the thick barrier, see Fig. 1. In a magnetic field $`B`$, the spectrometer exhibits a spin splitting $`\mathrm{\Delta }E_S=g_S\mu _BB`$, with $`g_S`$ the $`g`$-factor of the impurity and $`\mu _B`$ the Bohr magneton. The spin splitting in the emitter is $`\mathrm{\Delta }E=g\mu _BB`$. By applying a bias voltage $`V`$, the spin-split levels $`S_{}`$ and $`S_{}`$ are shifted in energy with respect to the emitter. They cross the Fermi level at $`V_{S_,}=(E_{S_,}\mu _E^{ch})/e\alpha `$ and two current steps result from resonant tunneling. Here, $`E_S_{}`$ and $`E_S_{}`$ are the energies of $`S_{}`$ and $`S_{}`$, $`\mu _E^{ch}`$ is the chemical potential in the emitter (measured from the band edge in the well and emitter, respectively), and $`\alpha `$ is the voltage-to-energy conversion coefficient. These two current steps can be clearly observed in Fig. 1(a), which shows the current-voltage characteristic $`I(V)`$ of our device at the base temperature of $`T=20`$ mK of a dilution refrigerator. The spin splitting between the steps increases as $`BI`$ is increased from 8.0 to 14.2 T. The overall shift of both steps to higher bias voltage is a consequence of the diamagnetic shift of $`E_S`$ to higher energy.

As the tunneling rate of the thick emitter barrier is much lower than that of the collector barrier, the tunneling current images the LDOS in the emitter below the Fermi level $`\mu _E`$ at the position of $`S`$. If spin is conserved during tunneling, the current on the short plateau of the first step is proportional to the spin-up component of the LDOS, $`I_1\nu _{}`$, while the current on the plateau of the second step reflects a superposition of both spin components, $`I_2\nu _{}+\nu _{}`$. At a bias voltage of $`V_,=(E_{S_,}E)/e\alpha `$, the spin-up and spin-down components of the LDOS at energy $`E`$ are probed by our experiment (measured from the band edge in the emitter). The oscillatory fine structure on the current plateaus in Fig. 1(a) reflects LDOS fluctuations in the heavily doped emitter . These fluctuations are more pronounced in the numerically obtained differential conductance $`G=dI/dV`$ shown in Fig. 1(b), where sharp peaks result from the current steps in the $`I(V)`$ data.

To determine the $`g`$-factor of the spectrometer, we plot $`\mathrm{\Delta }V_S=V_S_{}V_S_{}`$ as a function of magnetic field. Figure 2 shows that the data follows a straight line, which confirms that the two steps in the $`I(V)`$ curves reflect one spin-split impurity level. From the temperature dependence of the current step at $`B=0`$ T (not shown), we deduced a voltage-to-energy conversion coefficient of $`\alpha =0.50`$ for our device. Using this value, we obtain $`|g_S|=e\alpha (d\mathrm{\Delta }V_S/dB)/\mu _B=0.14`$. The $`g`$-factor of the impurity is different from the bulk value due to a change of the electron energy as a result of the double-barrier confinement and the penetration of the electron wavefunction into the barrier material. Note that the spin splitting of the spectrometer does not tell whether $`g_S`$ is positive or negative.

Now we turn our attention to the LDOS in the magnetic quantum limit, in which only one spin-split Landau band remains populated in the emitter. The transition to this limit occurs at $`B11.4`$ T. Figure 3(a) shows a color-scale image of the current as a function of bias voltage and magnetic field. It is dominated by two sharp current changes denoted by $`V_S_{}`$ and $`V_S_{}`$, which reflect the spin-split steps in the $`I(V)`$ curves in Fig. 1(a). The magnetic-field dependence of $`V_S_{}`$ has been eliminated by shifting the $`I(V)`$ curves in voltage. Between $`V_S_{}`$ and $`V_S_{}`$ as well as above $`V_S_{}`$, we observe pronounced fan-type resonances which reflect the local structure of Landau bands in the heavily doped emitter . These resonances in the LDOS are formed by quasi-one-dimensional interference of multiply-scattered electron waves subject to Landau quantization. The slope of the LDOS resonances is $`dV_E/dB=3.0`$ mV/T.

To study the spin dependence of the LDOS, we define two cuts at $`V_1`$ and $`V_2`$ \[dashed white lines in Fig. 3(a)\]. Figure 3(b) shows the corresponding currents $`I_1`$ and $`I_2`$ as a function of magnetic field. Assuming in a first attempt that the spin-up and spin-down components of the LDOS in the emitter are identical apart from a Zeeman shift $`\mathrm{\Delta }E`$, we can recalculate $`I_2(B)`$ from $`I_1(B)`$ as follows. From the sketch in Fig. 4, we deduce

$$I_1(B)\nu _{}(E_1,B)$$
(1)
and

$$I_2(B)\nu _{}(E_2,B)+\nu _{}(E_2+\mathrm{\Delta }E_S,B)$$
(2)
where $`E_1`$ and $`E_2`$ correspond to $`V_1`$ and $`V_2`$ and $`\mathrm{\Delta }E_S`$ denotes the spin splitting of the spectrometer. Using $`\nu _{}(E,B)=\nu _{}(E\mathrm{\Delta }E,B)`$, we replace the spin-down component of the LDOS in Eq. (2) by its spin-up counterpart. In addition, we exploit the linear, fan-type character of the LDOS resonances and use $`\nu _{}(E_1+E,B)=\nu _{}(E_1,B+[dB/dE_E]E)`$ to rescale all energy values to $`E_1`$. Defining $`E_{21}=E_2E_1`$, we obtain

$$\begin{array}{cc}\hfill I_2(B)& \nu _{}(E_1,B+\frac{dB}{dE_E}E_{21})\hfill \\ & +\nu _{}(E_1,B+\frac{dB}{dE_E}[E_{21}+\mathrm{\Delta }E_S\mathrm{\Delta }E])\hfill \end{array}$$
(3)
from Eq. (2). Using Eq. (1), we now rewrite Eq. (3) into

$$\begin{array}{cc}\hfill I_2(B)& =\gamma I_1\left(B+\frac{dB}{dE_E}E_{21}\right)\hfill \\ & +\gamma I_1\left(B+\frac{dB}{dE_E}[E_{21}+\mathrm{\Delta }E_S\mathrm{\Delta }E]\right).\hfill \end{array}$$
(4)
Here, $`\gamma `$ was introduced to compensate a difference in the coefficients of proportionality in Eqs.(1) and (2). It takes into account that the impurity can only be occupied by one electron at a time due to Coulomb interaction. For our analysis, we estimate $`\gamma `$ as $`I_2_B/2I_1_B=0.98`$ where $`\mathrm{}_B`$ symbolizes averaging over magnetic field. The slope $`dB/dE_E=1/[e\alpha (dV_E/dB)]=0.67`$ T/meV and $`|g_S|=0.14`$ have been measured and, in a first step, we assume $`g=0.44`$ .

Figure 4 compares the experimental data for $`I_2`$ with $`I_2^R`$ recalculated according to Eq. (4). Assuming in addition $`g_S=0.14`$, the strongest features in the LDOS are nicely reproduced by our reconstruction. Thus we conclude that the spin dependence of large fluctuations in the LDOS can be understood in terms of Zeeman splitting. These features vary slowly in magnetic field (on a scale of about 0.5 T or larger). The agreement also corroborates that both the spin-up and spin-down components of the LDOS are measured at the same location. The reconstruction does not work for $`g_S=+0.14`$, which strongly suggests that $`g_S<0`$ as predicted by calculations of the $`g`$-factor in a 10 nm wide quantum well . Strong deviations between $`I_2`$ and $`I_2^R`$ are observed for features that vary rapidly in magnetic field (on a scale shorter than 0.1 T). This is more evident at the bottom of Fig. 4, where we plot $`\mathrm{\Delta }I_2=I_2I_2^R`$ vs magnetic field.

In order to determine the actual $`g`$-factor in the heavily doped GaAs emitter, we minimize the root-mean-square deviation $`\sigma _{\mathrm{\Delta }I_2}=(\mathrm{\Delta }I_2(B)\mathrm{\Delta }I_2_B)^2_B^{1/2}`$ with respect to $`g`$. Figure 5 shows a clear minimum at $`g=0.35`$. The accuracy of $`g`$ is limited by the uncertainty of $`g_S`$. Assuming an error of $`\pm 0.01`$ in $`g_S`$, we calculate the error of $`g`$ to be smaller than $`\pm 0.02`$. In our experiment, the $`g`$-factor is measured at the Fermi energy of about $`25`$ meV above the band edge and at magnetic fields of $`10.515`$ T, while the standard value $`g=0.44`$ was obtained at the band edge and at zero field . The observed reduction of $`|g|`$ is in accordance with theory predicting a reduction of the $`g`$-factor in high magnetic fields and at finite electron densities .

Now we address the origin of the difference between $`I_2`$ and $`I_2^R`$ in Fig. 4. One issue to be considered is energy relaxation. The “holes” generated by tunneling in the emitter Fermi sea have an energy deficit $`e\alpha (VV_{S_,})`$ with respect to $`\mu _E`$ (see sketch in Fig. 1). The higher the bias voltage the larger is this energy deficit, resulting in a faster relaxation via electron-electron scattering. Inelastic scattering reduces the phase-coherence length in the emitter and thus suppresses the LDOS fluctuations. This effect is clearly visible in Fig. 1, where the fluctuations become less pronounced with increasing bias voltage. Similarly one might expect that $`I_2`$ is different from $`I_2^R`$, since $`I_2`$ is measured deeper below the Fermi level than $`I_1`$ which is used to calculate $`I_2^R`$. To check this idea, we calculated $`\mathrm{\Delta }I_2=I_2I_2^R`$ for 0.5 mV $`V_21`$ mV corresponding to an energy range of 0.25 meV. Interestingly, $`\mathrm{\Delta }I_2`$ is fairly independent from $`V_2`$, which demonstrates that energy relaxation is irrelevant in our experiment.

Thus the difference between $`I_2`$ and $`I_2^R`$ either indicates an inaccuracy in our reconstruction – possibly resulting from a magnetic-field dependence of the parameters $`\gamma `$, $`dB/dE_E`$, $`g_S`$, $`g`$ – or it gives evidence of additional spin effects beyond Zeeman splitting. While we will not be able to give a definite answer, it is worth while to examine the influence of extra spin effects on the LDOS.

One possibility is spin relaxation due to spin-orbit coupling in the emitter. In this case, we can estimate the spin-relaxation time as follows. Semiclassically speaking, the LDOS fluctuations are a result of quantum interference of diffusive electron waves at length scales between the mean free path and the phase-coherence length . Long interference paths correspond to weak features that vary rapidly in magnetic field. By calculating the differential conductance $`G=dI/dV`$, we favor small and rapid variations resulting from long interference paths \[see Fig. 1(b)\]. Due to the limited energy resolution $`\mathrm{}\mathrm{\Gamma }`$ of our spectrometer, the longest and dominant length scale in the conductance is $`l_\mathrm{\Gamma }=\sqrt{D/\mathrm{\Gamma }}`$, where $`D`$ is the diffusion constant in the emitter . The correlation magnetic field of the conductance fluctuations is

$$B_{c,G}\varphi _0/l_\mathrm{\Gamma }^2=\varphi _0/(D/\mathrm{\Gamma }).$$
(5)
By calculating the difference $`\mathrm{\Delta }I_2=I_2I_2^R`$, we eliminate the spin-independent features of the LDOS. Hence the correlation magnetic field of the fluctuations in $`\mathrm{\Delta }I_2`$ contains information about spin relaxation,

$$B_{c,\mathrm{\Delta }I_2}\varphi _0/l_s^2=\varphi _0/D\tau _s,$$
(6)
with $`l_s=\sqrt{D\tau _s}`$ being the spin-relaxation length. From Eqs. (5) and (6), we now calculate the spin-relaxation time

$$\tau _s(B_{c,G}/B_{c,\mathrm{\Delta }I_2})/\mathrm{\Gamma }.$$
(7)
All parameters in this expression have been measured in our experiment. In particular, it does not depend on the diffusion constant. From the measured $`G(B)`$ and $`\mathrm{\Delta }I_2(B)`$ data we find $`B_{c,G}0.01`$ T and $`B_{c,\mathrm{\Delta }I_2}0.03`$ T (in analogy to Ref. ), and from the width of the conductance peaks in Fig. 1(b) we deduce an energy resolution of $`\mathrm{}\mathrm{\Gamma }e\alpha \times 80`$ $`\mu `$V $`40`$ $`\mu `$eV. Using these values, we obtain a spin-relaxation time of $`\tau _s5`$ ps. It is much longer than the elastic scattering time $`\tau =0.1`$–0.2 ps, which we determined according to Ref. .

Another interesting spin effect may result from the presence of magnetic impurities in the emitter. The exchange field of fully spin-polarized impurities leads to different disorder potentials for spin-up and spin-down electrons and thus influences the LDOS. As there may only be residual magnetic impurities, their density will be by orders of magnitude lower than the Si concentration. Therefore, the spin-dependent features of the disorder potential vary on longer length scales than the spin-independent features. Slow variations in space correspond to small features on short field scales in the LDOS. Thus a small amount of magnetic impurities could explain why small features in the LDOS exhibit an irregular spin dependence, while large features simply show a regular Zeeman splitting. To clarify the influence of magnetic impurities on the LDOS, we suggest to deliberately dope the GaAs contacts of a double-barrier heterostructure with magnetic impurities such as Mn.

In summary, we used a spin-split impurity level as spectrometer to study the interplay of Zeeman splitting and disorder in Si doped GaAs via resonant tunneling. Both the spin-up and spin-down components of the LDOS show fan-type resonances as a consequence of Landau quantization in the bulk semiconductor. While the main resonances in the LDOS exhibit a Zeeman shift that can be best described using a $`g`$-factor of $`g=0.35`$, the smaller resonances show a complicated irregular behavior. Spin relaxation due to spin-orbit coupling is one possible origin of this deviation from Zeeman splitting.

We gratefully acknowledge A. Förster and H. Lüth for growing the double-barrier heterostructure and V.I. Fal’ko for useful discussions.


format: text

text:
# Searching for Spaceships

## 1 Introduction

John Conway’s Game of Life has fascinated and inspired many enthusiasts, due to the emergence of complex behavior from a very simple system. One of the many interesting phenomena in Life is the existence of gliders and spaceships: small patterns that move across space. When describing gliders, spaceships, and other early discoveries in Life, Martin Gardner wrote (in 1970) that spaceships “are extremely hard to find” . Very small spaceships can be found by human experimentation, but finding larger ones requires more sophisticated methods. Can computer software aid in this search? The answer is yes – we describe here a program, gfind, that can quickly find large low-period spaceships in Life and many related cellular automata.

Among the interesting new patterns found by gfind are the “weekender” $`\mathrm{2}c/\mathrm{7}`$ spaceship in Conway’s Life (Figure 1, right), the “dragon” $`c/\mathrm{6}`$ Life spaceship found by Paul Tooke (Figure 1, left), and a $`c/\mathrm{7}`$ spaceship in the Diamoeba rule (Figure 2, top). The middle section of the Diamoeba spaceship simulates a simple one-dimensional parity automaton and can be extended to arbitrary lengths. David Bell discovered that two back-to-back copies of these spaceships form a pattern that fills space with live cells (Figure 2, bottom). The existence of infinite-growth patterns in Diamoeba had previously been posed as an open problem (with a $50 bounty) by Dean Hickerson in August 1993, and was later included in a list of related open problems by Gravner and Griffeath . Our program has also found new spaceships in well known rules such as HighLife and Day&Night as well as in thousands of unnamed rules.

As well as providing a useful tool for discovering cellular automaton behavior, our work may be of interest for its use of state space search techniques. Recently, Buckingham and Callahan  wrote “So far, computers have primarily been used to speed up the design process and fill gaps left by a manual search. Much potential remains for increasing the level of automation, suggesting that Life may merit more attention from the computer search community.” Spaceship searching provides a search problem with characteristics intriguingly different from standard test cases such as the 15-puzzle or computer chess, including a large state space that fluctuates in width instead of growing exponentially at each level, a tendency for many branches of the search to lead to dead ends, and the lack of any kind of admissable estimate for the distance to a goal state. Therefore, the search community may benefit from more attention to Life.

The software described here, a database of spaceships in Life-like automata, and several programs for related computations can be found online at http://www.ics.uci.edu/~eppstein/ca/.

## 2 A Brief History of Spaceship Searching

According to Berlekamp et al. , the $`c/\mathrm{4}`$ diagonal glider in Life was first discovered by simulating the evolution of the R-pentomino, one of only 21 connected patterns of at most five live cells. This number is small enough that the selection of patterns was likely performed by hand. Gliders can also be seen in the evolution of random initial conditions in Life as well as other automata such as B3/S13 , but this technique often fails to work in other automata due to the lack of large enough regions of dead cells for the spaceships to fly through. Soon after the discovery of the glider, Life’s three small $`c/\mathrm{2}`$ orthogonal spaceships were also discovered.

Probably the first automatic search method developed to look for interesting patterns in Life and other cellular automata, and the method most commonly programmed, is a brute force search that tests patterns of bounded size, patterns with a bounded number of live cells, or patterns formed out of a small number of known building blocks. These tests might be exhaustive (trying all possible patterns) or they might perform a sequence of trials on randomly chosen small patterns. Such methods have found many interesting oscillators and other patterns in Life, and Bob Wainwright collected a large list of small spaceships found in this way for many other cellular automaton rules . Currently, it is possible to try all patterns that fit within rectangles of up to $`\mathrm{7}\times \mathrm{8}`$ cells (assuming symmetric initial conditions), and this sort of exhaustive search can sometimes find spaceships as large as $`\mathrm{1}\mathrm{2}\times \mathrm{1}\mathrm{5}`$ (Figure 3). However, brute force methods have not been able to find spaceships in Life with speeds other than $`c/\mathrm{2}`$ and $`c/\mathrm{4}`$.

The first use of more sophisticated search techniques came in 1989, when Dean Hickerson wrote a backtracking search program which he called LS. For each generation of each cell inside a fixed rectangle, LS stored one of three states: unknown, live, or dead. LS then attempted to set the state of each unknown cell by examining neighboring cells in the next and previous generations. If no unknown cell’s state could be determined, the program performed a depth first branching step in which it tried both possible states for one of the cells. Using this program, Hickerson discovered many patterns including Life’s $`c/\mathrm{3}`$, $`c/\mathrm{4}`$, and $`\mathrm{2}c/\mathrm{5}`$ orthogonal spaceships. Hartmut Holzwart used a similar program to find many variant $`c/\mathrm{2}`$ and $`c/\mathrm{3}`$ spaceships in Life, and related patterns including the “spacefiller” in which four $`c/\mathrm{2}`$ spaceships stretch the corners of a growing diamond shaped still life. David Bell reimplemented this method in portable C, and added a fourth “don’t care” state; his program, lifesrc, is available at http://www.canb.auug.org.au/~dbell/programs/lifesrc-3.7.tar.gz.

In 1996, Tim Coe discovered another Life spaceship, moving orthogonally at speed $`c/\mathrm{5}`$, using a program he called knight in the hope that it could also find “knightships” such as those in Figure 5. Knight used breadth first search (with a fixed amount of depth-first lookahead per node to reduce the space needs of BFS) on a representation of the problem based on de Bruijn graphs (described in Section 4). The search took 38 cpu-weeks of time on a combination of Pentium Pro 133 and Hypersparc processors. The new search program we describe here can be viewed as using a similar state space with improved search algorithms and fast implementation techniques. Another recent program by Keith Amling also uses a state space very similar to Coe’s, with a depth first search algorithm.

Other techniques for finding cellular automaton patterns include complementation of nondeterministic finite automata, by Jean Hardouin-Duparc ; strong connectivity analysis of de Bruijn graphs, by Harold McIntosh ; randomized hill-climbing methods for minimizing the number of cells with incorrect evolution, by Paul Callahan (http://www.cs.jhu.edu/~callahan/stilledit.html); a backtracking search for still life backgrounds such that an initial perturbation remains bounded in size as it evolves, by Dean Hickerson; Gröbner basis methods, by John Aspinall; and a formulation of the search problem as an integer program, attempted by Richard Schroeppel and later applied with more success by Robert Bosch . However to our knowledge none of these techniques has been used to find new spaceships.

## 3 Notation and Classification of Patterns

We consider here only outer totalistic rules, like Life, in which any cell is either “live” or “dead” and in which the state of any cell depends only on its previous state and on the total number of live neighbors among the eight adjacent cells of the Moore neighborhood (Figure 4). For some results on spaceships in cellular automata with larger neighborhoods, see Evans’ thesis . Outer totalistic rules are described with a string of the form B$`x_\mathrm{1}x_\mathrm{2}\mathrm{}`$/S$`y_\mathrm{1}y_\mathrm{2}\mathrm{}`$ where the $`x_i`$ are digits listing the number of neighbors required for a cell to be born (change from dead to live) while the $`y_i`$ list the number of neighbors required for a cell to survive (stay in the live state after already being live). For instance, Conway’s Life is described in this notation as B3/S23.

A spaceship is a pattern which repeats itself after some number $`p`$ of generations, in a different position from where it started. We call $`p`$ the period of the spaceship. If the pattern moves $`x`$ units horizontally and $`y`$ units vertically every $`p`$ steps, we say that it has slope $`y/x`$ and speed $`max(|x|,|y|)c/p`$, where $`c`$ denotes the maximum speed at which information can propagate in the automaton (the so-called speed of light). Most known spaceships move orthogonally or diagonally, and many have an axis of symmetry parallel to the direction of motion. Others, such as the glider and small $`c/\mathrm{2}`$ spaceships in Life, have glide-reflect symmetry: a mirror image of the original pattern appears in generations $`p/\mathrm{2}`$, $`\mathrm{3}p/\mathrm{2}`$, etc.; spaceships with this type of symmetry must also move orthogonally or diagonally. However, a few asymmetric spaceships move along lines of slope 2 or even 3/2 (Figure 5). According to Berlekamp et al. , there exist spaceships in Life that move with any given rational slope, but the argument for the existence of such spaceships is not very explicit and would lead to extremely large patterns.

Related types of patterns include oscillators (patterns that repeat in the same position), still lifes (oscillators with period 1), puffers (patterns which repeat some distance away after a fixed period, leaving behind a trail of discrete patterns such as oscillators, still lifes, and spaceships), rakes (spaceship puffers), guns (oscillators which send out a moving trail of discrete patterns such as spaceships or rakes), wickstretchers (patterns which leave behind one or more connected stable or oscillating regions), replicators (patterns which produce multiple copies of themselves), breeders (patterns which fill a quadratically-growing area of space with discrete patterns, for instance replicator puffers or rake guns), and spacefillers (patterns which fill a quadratically-growing area with one or more connected patterns). See Paul Callahan’s Life Pattern Catalog (http://www.cs.jhu.edu/~callahan/patterns/contents.html) for examples of many of these types of patterns in Life, or http://www.ics.uci.edu/~eppstein/ca/replicators/ for a number of replicator-based patterns in other rules.

We can distinguish among several classes of spaceships, according to the methods that work best for finding them.

* Spaceships with small size but possibly high period can be found by brute force search. The patterns depicted in Figures 3 and 5 fall into this class, as do Life’s glider and c/2 spaceships.
* Spaceships in which the period and one dimension are small, while the other dimension may be large, can be found by search algorithms similar to the ones described in this paper. The new Life spaceships in Figure 1 fall into this class. For a more extreme example of a low period ship which is long but narrow, see Figure 6. “Small” is a relative term—our search program has found $`c/\mathrm{2}`$ spaceships with minimum dimension as high as 42 (Figure 10) as well as narrower spaceships with period as high as nine.
* Sometimes small non-spaceship objects, such as puffers, wickstretchers, or replicators, can be combined by human engineering into a spaceship. For instance, in the near-Life rule B37/S23, the R-pentomino pattern acts as an (unstable) puffer. Figure 7 depicts a $`\mathrm{9}c/\mathrm{2}\mathrm{8}`$ spaceship in which a row of six pentominoes stabilize each other while leaving behind a trail of still lifes, which are cleaned up by a second row of four pentominoes. Occasionally, spaceships found by our search software will appear to have this sort of structure (Figure 10). The argument for the existence of Life spaceships with any rational slope would also lead to patterns of this type. Several such spaceships have been constructed by Dean Hickerson, notably the $`c/\mathrm{1}\mathrm{2}`$ diagonal “Cordership” in Conway’s Life, discovered by him in 1991. Hickerson’s web page http://math.ucdavis.edu/~dean/RLE/slowships.html has more examples of structured ships.
* The remaining spaceships have large size, large period, and little internal structure. We believe such spaceships should exist, but none are known and we know of no effective method for finding them.

## 4 State Space

Due to the way our search is structured, we need to arrange the rows from all phases of the pattern we are searching for into a single sequence. We now show how to do this in a way that falls out naturally from the motion of the spaceship.

Suppose we are searching for a spaceship that moves $`k`$ units down every $`p`$ generations. For simplicity of exposition, we will assume that $`gcd(k,p)=\mathrm{1}`$. We can then think of the spaceship we are searching for as living in a cellular automaton modified by shifting the grid upward $`k/p`$ units per generation. In this modified automaton, the shifting of the grid exactly offsets the motion of the spaceship, so the ship acts like an oscillator instead of like a moving pattern. We illustrate this shifted grid with the turtle, a $`c/\mathrm{3}`$ spaceship in Conway’s Life (Figure 8).

Because of the shifted grid, and the assumption that $`gcd(k,p)=\mathrm{1}`$, each row of each phase of the pattern exists at a distinct vertical position. We form a doubly-infinite sequence

$$\mathrm{}r[\mathrm{2}],r[\mathrm{1}],r[\mathrm{0}],r[\mathrm{1}],r[\mathrm{2}]\mathrm{}$$

of rows, by taking each row from each phase in order by the rows’ vertical positions (Figure 9). For any $`i`$, the three rows $`r[i\mathrm{2}p],r[ip],r[i]`$ form a contiguous height-three strip in a single phase of the pattern, and we can apply the cellular automaton rule within this strip to calculate

$$r[ip+k]=evolve(r[i\mathrm{2}p],r[ip],r[i]).$$
($``$)
Conversely, any doubly-infinite sequence of rows, in which equation ($``$) is satisfied for all $`i`$, and in which there are only finitely many live cells, corresponds to a spaceship or sequence of spaceships. Further, any finite sequence of rows can be extended to such a doubly infinite sequence if the first and last $`\mathrm{2}p`$ rows of the finite sequence contain only dead cells.

Diagonal spaceships, such as Life’s glider, can be handled in this framework by modifying equation ($``$) to shift rows $`i`$, $`ip`$, and $`i\mathrm{2}p`$ with respect to each other before performing the evolution rule. We can also handle glide-reflect spaceships such as Life’s small $`c/\mathrm{2}`$ spaceships, by modifying equation ($``$) to reverse the order of the cells in row $`ip+k`$ (when $`k`$ is odd) or in the two rows $`ip`$ and $`ip+k`$ (when $`p`$ is odd). Note that, by the assumption that $`gcd(k,p)=\mathrm{1}`$, at least one of $`k`$ and $`p`$ will be odd. In these cases, $`p`$ should be considered as the half-period of the spaceship, the generation at which a flipped copy of the original pattern appears. Searches for which $`gcd(k,p)>\mathrm{1}`$ can be handled by adjusting the indices in equation ($``$) depending on the phase to which row $`r[i]`$ belongs.

Our state space, then, consists of finite sequences of rows, such that equation ($``$) is satisfied whenever all four rows in the equation belong to the finite sequence. The initial state for our search will be a sequence of $`\mathrm{2}p`$ rows, all of which contain only dead cells. If our search discovers another state in which the last $`\mathrm{2}p`$ rows also contain only dead cells, it outputs the pattern formed by every $`p`$th row of the state as a spaceship.

As in many game playing programs, we use a transposition table to detect equivalent states, avoid repeatedly searching the same states, and stop searching in a finite amount of time even when the state space may be infinite. We would like to define two states as being equivalent if, whenever one of them can be completed to form a spaceship, the same completion works for the other state as well. However, this notion of equivalence seems too difficult to compute (it would require us to be able to detect states that can be completed to spaceships, but if we could do that then much of our search could be avoided). So, we use a simpler sufficient condition: two states are equivalent if their last $`\mathrm{2}p`$ rows are identical. If our transposition table detects two equivalent states, the longer of the two is eliminated, since it can not possibly lead to the shortest spaceship for that rule and period.

We can form a finite directed graph, the de Bruijn graph , by forming a vertex for each equivalence class of states in our state space, and an edge between two vertices whenever a state in one equivalence class can be extended to form a state in the other class. The size of the de Bruijn graph provides a rough guide to the complexity of a spaceship search problem. If we are searching for patterns with width $`w`$, the number of vertices in this de Bruijn graph is $`\mathrm{2}^{\mathrm{2}pw}`$. For instance, in the search for the weekender spaceship, the effective width was nine (due to an assumption of bilateral symmetry), so the de Bruijn graph contained $`\mathrm{2}^{\mathrm{1}\mathrm{2}\mathrm{6}}`$ vertices. Fortunately, most of the vertices in this graph were unreachable from the start state.

Coe’s search program knight uses a similar state space formed by sequences of pattern rows, but only uses the rows from a single phase of the spaceship. In place of our equation ($``$), he evolves subsequences of $`\mathrm{2}p+\mathrm{1}`$ rows for $`p`$ generations and tests that the middle row of the result matches the corresponding row of the subsequence. As with our state space, one can form a (different) de Bruijn graph by forming equivalence classes according to the last $`\mathrm{2}p`$ rows of a state. Coe’s approach has some advantages; for instance, it can find patterns in which some phases exceed the basic search width (as occurred in Coe’s $`c/\mathrm{5}`$ spaceship). However, it does not seem to allow the fast neighbor-finding techniques we describe in Section 7.

## 5 Search Strategies

There are many standard algorithms for searching state spaces , however each has some drawbacks in our application:

* Depth first search requires setting an arbitrary depth limit to avoid infinite recursion. The patterns it finds may be much longer than necessary, and the search may spend a long time exploring deep regions of the state space before reaching a spaceship. Further, DFS does not make effective use of the large amounts of memory available on modern computers.
* Breadth first search is very effective for small searches, but quickly runs out of memory for larger searches, even when large amounts of memory are available.
* Depth first iterative deepening  has been proposed as a method of achieving the fast search times of breadth first search within limited space. However, our state space often does not have the exponential growth required for iterative deepening to be efficient; rather, as the search progresses from level to level the number of states in the search frontier can fluctuate up and down, and typically has a particularly large bulge in the earlier levels of the search. The overall depth of the search (and hence the number of deepening iterations) can often be as large as several hundred. For these reasons, iterative deepening can be much slower than breadth first search. Further, the transposition table used to detect equivalent states does not work as well with depth first as with breadth first search: to save space, we represent this table as a collection of pointers to states, rather than explicitly listing the $`\mathrm{2}p`$ rows needed to determine equivalence, so when searching depth first we can only detect repetitions within the current search path. Finally, depth first search does not give us much information about the speed at which the search is progressing, which we can use to narrow the row width when the search becomes too slow.
* Other techniques such as the A algorithm, recursive best first search, and space-bounded best first search, depend on information unavailable in our problem, such as varying edge weights or admissable estimates of the distance to a solution.

Therefore, we developed a new search algorithm that combines the best features of breadth first and iterative deepening search, and that takes advantage of the fact that, in our search problem, many branches of the search eventually lead only to dead ends. Our method resembles the MREC algorithm of Sen and Bagchi , in that we perform deepening searches from the breadth first search frontier, however unlike MREC we use the deepening stages to prune the search tree, allowing additional breadth first searching.

Our search algorithm begins by performing a standard breadth first search. We represent each state as a single row together with a pointer to its predecessor, so the search must maintain the entire breadth first search tree. By default, we allocate storage for $`\mathrm{2}^{\mathrm{2}\mathrm{2}}`$ nodes in the tree, which is adequate for small searches yet well within the memory limitations of most computers.

On larger searches (longer than a minute or so), the breadth first search will eventually run out of space. When this happens, we suspend the breadth first search and perform a round of depth first search, starting at each node of the current breadth first search queue. This depth first search has a depth limit which is normally $`\delta `$ levels beyond the current search frontier, for a small value $`\delta `$ that we set in our implementation to equal the period of the pattern we are searching for. (Setting $`\delta `$ to a fixed small constant would likely work as well.) However, if a previous round of depth first searching reached a level past the current search frontier, we instead limit the new depth first search round to $`\delta `$ levels beyond the previous round’s limit.

When we perform a depth first search from a BFS queue node, one of three things can happen. First, we may discover a spaceship; in that case we terminate the entire search. Second, we may reach the depth limit; in that case we terminate the depth first search and move on to the next BFS queue node. Third, the depth first search may finish normally, without finding any spaceships or deep nodes. In this case, we know that the root of the search leads only to dead ends, and we remove it from the breadth first search queue.

After we have performed this depth first search for all nodes of the queue, we compact the remaining nodes and continue with the previously suspended breadth first search. Generally, only a small fraction of the previous breadth first search tree remains after the compaction, leaving plenty of space for the next round of breadth first search.

There are two common modes of behavior for this searching algorithm, depending on how many levels of breadth first searching occur between successive depth first rounds. If more than $`\delta `$ levels occur between each depth first search round, then each depth first search is limited to only $`\delta `$ levels beyond the breadth first frontier, and the sets of nodes searched by successive depth first rounds are completely disjoint from each other. In this case, each node is searched at most twice (once by the breadth first and once by the depth first parts of our search), so we only incur a constant factor slowdown over the more memory-intensive pure breadth first search. In the second mode of behavior, successive depth first search rounds begin from frontiers that are fewer than $`\delta `$ levels apart. If this happens, the $`i`$th round of depth first search will be limited to depth $`i\delta `$, so the search resembles a form of iterative deepening. Unlike iterative deepening, however, the breadth first frontier always makes some progress, permanently removing nodes from the actively searched part of the state space. Further, the early termination of depth-first searches when they reach deep nodes allows our algorithm to avoid searching large portions of the state space that pure iterative deepening would have to examine.

On typical large searches, the amount of deepening can be comparable to the level of the breadth-first search frontier. For instance, in the search for the weekender, the final depth first search round occurred when the frontier had reached level 90, and this round searched from each frontier node to an additional 97 levels. We allow the user to supply a maximum value for the deepening amount; if this maximum is reached, the state space is pruned by reducing the row width by one cell and the depth first search limit reverts back to $`\delta `$.

There is some possibility that a spaceship found in one of the depth first searches may be longer than the optimum, but this has not been a problem in practice. Even this small amount of suboptimality could be averted by moving on to the next node instead of terminating the search when the depth first phase of the algorithm discovers a spaceship.

## 6 Lookahead

Suppose that our search reaches state

$$S=r[\mathrm{0}],r[\mathrm{1}],\mathrm{}r[i\mathrm{1}].$$

The natural set of neighboring states to consider would be all sequences of rows

$$r[\mathrm{0}],r[\mathrm{1}],\mathrm{}r[i\mathrm{1}],r[i]$$

where the first $`i`$ rows match state $`S`$ and we try all choices of row $`r[i]`$ that satisfy equation ($``$).

However, it is likely that some of these choices will result in inconsistent states for which equation ($``$) can not be satisfied the next time the new row $`r[i]`$ is involved in the equation. Since that next time will not occur until we choose row $`r[i+pk]`$, any work performed in the intermediate levels of the search between these two rows could be wasted. To avoid this problem, when making choices for row $`r[i]`$, we simultaneously search for pairs of rows $`r[i]`$ and $`r[i+pk]`$ satisfying both equation ($``$) and its shifted form

$$r[i]=evolve(r[ipk],r[ik],r[i+pk]).$$
(L)
Note that $`r[ipk]`$ and $`[ik]`$ are both already present in state $`S`$ and so do not need to be searched for. We use as the set of successors to row $`S`$ the sequences of rows $`r[\mathrm{0}],r[\mathrm{1}],\mathrm{}r[i\mathrm{1}],r[i]`$ such that equation ($``$) is true and equation (L) has a solution.

One could extend this idea further, and (as well as searching for rows $`r[i]`$ and $`r[i+pk]`$) search for two additional rows $`r[i+p\mathrm{2}k]`$ and $`r[i+\mathrm{2}p\mathrm{2}k]`$ such that the double lookahead equations

$$\begin{array}{ccc}\hfill r[ik]& =& evolve(r[ip\mathrm{2}k],r[i\mathrm{2}k],r[i+p\mathrm{2}k])\text{ and}\hfill \\ \hfill r[i+pk]& =& evolve(r[i\mathrm{2}k],r[i+p\mathrm{2}k],r[i+\mathrm{2}p\mathrm{2}k])\hfill \end{array}$$
(LL)
have a solution. However this double lookahead technique would greatly increase the cost of searching for the successor states of $`S`$ and provide only diminished returns. In our implementation, we use a much cheaper approximation to this technique: for every three consecutive cells of the row $`r[i+pk]`$ being searched for as part of our single lookahead technique, we test that there could exist five consecutive cells in rows $`r[i+p\mathrm{2}k]`$ and $`r[i+\mathrm{2}p\mathrm{2}k]`$ satisfying equations (LL) for those three cells. The advantage of this test over the full search for rows $`r[i+p\mathrm{2}k]`$ and $`r[i+\mathrm{2}p\mathrm{2}k]`$ is that it can be performed with simple table lookup techniques, described in the next section.

For the special case $`p=\mathrm{2}`$, the double lookahead technique described above does not work as well, because $`r[i+p\mathrm{2}k]=r[i]`$. Instead, we perform a reachability computation on a de Bruijn graph similar to the one used for our state space, five cells wide with don’t-care boundary conditions, and test whether triples of consecutive cells from the new rows $`r[i]`$ and $`r[i+pk]`$ correspond to a de Bruijn graph vertex that can reach a terminal state. This technique varies considerably in effectiveness, depending on the rule: for Life, 18.5% of the 65536 possible patterns are pruned as being unable to reach a terminal state, and for B27/S0 (Figure 10) the number is 68.3%, but for many rules it can be 0%.

## 7 Fast Neighbor-Finding Algorithm

To complete our search algorithm, we need to describe how to find the successors of each state. That is, given a state

$$S=r[\mathrm{0}],r[\mathrm{1}],\mathrm{}r[i\mathrm{1}]$$

we wish to find all possible rows $`r[i]`$ such that (1) $`r[i]`$ satisfies equation ($``$), (2) there exists another row $`r[i+pk]`$ such that rows $`r[i]`$ and $`r[i+pk]`$ satisfy equation (L), and (3) every three consecutive cells of $`r[i+pk]`$ can be part of a solution to equations (LL).

Note that (because of our approximation to equations (LL)) all these constraints involve only triples of adjacent cells in unknown rows. That is, equation ($``$) can be phrased as saying that every three consecutive cells of $`r[i]`$, $`r[ip]`$, and $`r[i\mathrm{2}p]`$ form a $`\mathrm{3}\times \mathrm{3}`$ square such that the result of the evolution rule in the center of the square is correct; and equation (L) can be phrased similarly. For this reason, we need a representation of the pairs of rows $`r[i],r[i+pk]`$ in which we can access these triples of adjacent cells.

Such a representation is provided by the graph depicted in Figure 11. This graph (which like our state space can be viewed as a kind of de Bruijn graph) is formed by a sequence of columns of vertices. Each column contains 16 vertices, representing the 16 ways of forming a $`\mathrm{2}\times \mathrm{2}`$ block of cells. We connect a vertex in one column with a vertex in the next column whenever the two blocks overlap to form a single $`\mathrm{2}\times \mathrm{3}`$ block of cells; that is, an edge exists whenever the right half of the left $`\mathrm{2}\times \mathrm{2}`$ block matches the left half of the right $`\mathrm{2}\times \mathrm{2}`$ block.

If we have any two rows of cells, one placed above the other, we can form a path in this graph in which the vertices correspond to $`\mathrm{2}\times \mathrm{2}`$ blocks drawn from the pair of rows. Conversely, the fact that adjacent blocks are required to match implies that any path in this graph corresponds to a pair of rows.

Since each triple of adjacent cells from the two rows corresponds to an edge in this graph, the constraints of equations ($``$), (L), and (LL) can be handled simply by removing the edges from the graph that correspond to triples not satisfying those constraints. In this constrained subgraph, any path corresponds to a pair of rows satisfying all the constraints. Thus, our problem has been reduced to finding the constrained subgraph, and searching for paths through it between appropriately chosen start and end terminals. The choice of which vertices to use as terminals depends on the symmetry type of the pattern we are searching for.

This search can be understood as being separated into three stages, although our actual implementation interleaves the first two of these.

In the first stage, we find the edges of the graph corresponding to blocks of cells satisfying the given constraints. We represent the set of 64 possible edges between each pair of columns in the graph (as shown in Figure 11) as a 64-bit quantity, where a bit is nonzero if an edge exists and zero otherwise. The set of edges corresponding to blocks satisfying equation ($``$) can be found by a table lookup with an index that encodes the values of three consecutive cells of $`r[ip]`$ and $`r[i\mathrm{2}p]`$ together with one cell of $`r[ip+k]`$. Similarly, the set of edges corresponding to blocks satisfying equation (L) can be found by a table lookup with an index that encodes the values of three consecutive cells of $`r[ik]`$ and $`r[ipk]`$. In our implementation, we combine these two constraints into a single table lookup. Finally, the set of edges corresponding to blocks satisfying equations (LL) can be found by a table lookup with an index that encodes the values of five consecutive cells of $`r[i\mathrm{2}k]`$ and $`r[ip\mathrm{2}k]`$ together with three consecutive cells of $`r[ik]`$. The sets coming from equations ($``$), (L), and (LL) are combined with a simple bitwise Boolean and operation. The various tables used by this stage depend on the cellular automaton rule, and are precomputed based on that rule before we do any searching.

In the second stage, we compute a 16-bit quantity for each column of vertices, representing the set of vertices in that column that can be reached from the start terminal. This set can be computed from the set of reachable vertices in the previous column by a small number of shift and mask operations involving the sets of edges computed in the previous stage.

In the third stage, we wish to find the actual rows $`r[i]`$ that correspond to the paths in the constrained graph. We perform a backtracking search for these rows, starting with the end terminal of the graph. At each step, we maintain a 16-bit quantity, representing the set of vertices in the current column of the graph that could reach the end terminal by a path matching the current partial row. To find the possible extensions of a row, we find the predecessors of this set of vertices in the graph by a small number of shift and mask operations (resembling those of the previous stage) and separate these predecessors into two subsets, those for which the appropriate cell of $`r[i]`$ is alive or dead. We then continue recursively in one or both of these subsets, depending on whichever of the two has a nonempty intersection with the set of vertices in the same column that can reach the start vertex.

Because of the reachability computation in the second stage, the third stage never reaches a dead end in which a partial row can not be extended. Therefore, this algorithm spends a total amount of time at most proportional to the width of the rows (in the first two stages) plus the width times the number of successor states (in the third stage). The time for the overall search algorithm is bounded by the width times the number of states reached.

A simplified example of this graph representation of our problem is depicted in Figure 12, which shows the graph formed by the rows from the Life turtle pattern depicted in Figure 9. To reduce the complexity of the figure, we have only incorporated equation ($``$), and not the two lookahead equations. Therefore, the vertices in each column represent the states of two adjacent cells in row $`r[i]`$ only, instead of also involving row $`r[i+pk]`$, and there are four vertices per column instead of sixteen. Each edge represents the state of three adjacent cells of $`r[i]`$, and connects vertices corresponding to the left two and right two of these cells; we have marked each edge with the corresponding cell states.

Due to the symmetry of the turtle, the effective search width is six, so we need to enforce equation ($``$) for six cells of row $`r[ip+k]`$. Six of the seven columns of edges in the graph correspond to these cells; the seventh represents a cell outside the search width which must remain blank to prevent the pattern from growing beyond its assigned bounds. Below each of these seven columns, we have shown the cells from previous rows $`r[ip]`$, $`r[i\mathrm{2}p]`$, and $`r[ip+k]`$ which determine the set of edges in the column and which are concatenated together to form the table index used to look up this set of edges.

The starting vertices in this example are the top and bottom left ones in the graph, which represent the possible states for the center two cells of $`r[i]`$ that preserve the pattern’s symmetry. The destination vertex is the upper right one; it represents the state of two cells in row $`r[i]`$ beyond the given search width, so both must be blank. Each path from a start vertex to the destination vertex represents a possible choice for the cells in row $`r[i]`$ that would lead to the correct evolution of row $`r[ip+k]`$ according to the rules of Conway’s life. There are 13 such paths in the graph shown. The reachability information computed in the second stage of the algorithm is depicted by marking unreachable vertices and edges in gray; in this example, as well as the asymmetric states in the first column, there is one more unreachable vertex.

For this simplified example, a list of all 13 paths in this graph could be found in the third stage by a recursive depth first search from the destination vertex backwards, searching only the black edges and vertices. Thus even this simplified representation reduces the number of row states that need to be considered from 64 to 13, and automatically selects only those states for which the evolution rule leads to the correct outcome. The presence of lookahead complicates the third stage in our actual program since multiple paths can correspond to the same value of row $`r[i]`$; the recursive search procedure described above finds each such value exactly once.

## 8 Conclusions

To summarize, we have described an algorithm that finds spaceships in outer totalistic Moore neighborhood cellular automata, by a hybrid breadth first iterative deepening search algorithm in a state space formed by partial sequences of pattern rows. The algorithm represents the successors of each state by paths in a regularly structured graph with roughly $`\mathrm{1}\mathrm{6}w`$ vertices; this graph is constructed by performing table lookups to quickly find the sets of edges representing the constraints of the cellular automaton evolution rule and of our lookahead formulations. We use this graph to find a state’s successors by performing a 16-way bit-parallel reachability algorithm in this graph, followed by a recursive backtracking stage that uses the reachability information to avoid dead ends. Empirically, the algorithm works well, and is able to find large new patterns in many cellular automaton rules.

This work raises a number of interesting research questions, beyond the obvious one of what further improvements are possible to our search program:

* The algorithm described here, and the two search algorithms previously used by Hickerson and Coe, use three different state spaces. Do these spaces really lead to different asymptotic search performance, or are they merely three different ways of looking at the same thing? Can one make any theoretical arguments for why one space might work better than another?
* Is it possible to explain the observed fluctuations in size of the levels of the state space? A hint of an explanation for the fact that fluctuations exist comes from the idea that we typically run as narrow as possible a search as we can to find a given period spaceship. Increasing the width seems to increase the branching factor, and perhaps we should expect to find spaceships as soon as the state space develops infinite branches, at which point the branching factor will likely be quite close to one. However this rough idea depends on the unexplained assumption that the start of a spaceship is harder to find than the tail, and it does not explain other features of the state space size such as a large bulge near the early levels of the search.
* We have mentioned in Section 4 that one can use the size of the de Bruijn graph as a rough guide to the complexity of the search. However, our searches typically examine far fewer nodes than are present in the de Bruijn graph. Further, there seem to be other factors that can influence the search complexity; for instance, increasing $`k`$ seems to decrease the running time, so that e.g. a width-nine search for a $`c/\mathrm{7}`$ ship in Life would likely take much more time than the search for the weekender. Can we find a better formula for predicting search run time?
* What if anything can one say about the computational complexity of spaceship searching? Is the problem of determining whether a given outer totalistic rule has a spaceship of a given speed or period even decidable?
* David Bell and others have had success finding large spaceships with “arms” (Figure 13) by examining partial results from an automatic search, placing “don’t care” cells at appropriate connection points, and then doing secondary searches for arms that can complete the pattern from each connection point. To what extent can this human-guided search procedure be automated?
* What other types of patterns can be found by our search techniques? For instance, one possibility would be a search for predecessors of a given pattern. The rows of each predecessor satisfy a consistency condition similar to equation ($``$), and it would not be difficult to incorporate a lookahead technique similar to the one we use for spaceship searching. However due to the fixed depth of the search it seems that depth first search would be a more appropriate strategy than the breadth first techniques we are using in our spaceship searches.
* Carter Bays has had some success using brute force methods to find small spaceships in various three-dimensional rules , and rules on the triangular planar lattice . How well can other search methods such as the ones described here work for these types of automata?
* To what areas other than cellular automata can our search techniques be applied? A possible candidate is in document improvement, where Bern, Goldberg, and others  have been developing algorithms that look for a single high-resolution image that best matches a given set of low-resolution samples of the same character or image. Our graph based techniques could be used to replace a local optimization technique that changes a single pixel at a time, with a technique that finds an optimal assignment to an entire row of pixels, or to any other linear sequence such as the set of pixels around the boundary of an object.

## Acknowledgements

Thanks go to Matthew Cook, Nick Gotts, Dean Hickerson, Harold McIntosh, Gabriel Nivasch, and Bob Wainwright for helpful comments on drafts of this paper; to Noam Elkies, Rich Korf, and Jason Summers for useful suggestions regarding search algorithms; to Keith Amling, David Bell, and Tim Coe, for making available the source codes of their spaceship searching programs; and to Richard Schroeppel and the members of the Life mailing list for their support and encouragement of cellular automaton research.


text:
# Molecular Dynamics Studies of Sequence-directed Curvature in Bending Locus of Trypanosome Kinetoplast DNA

## Introduction

It is now generally accepted that the double helical DNA can somehow translate its base pair sequence in tertiary structural forms. The simplest such form is a bend. Large bends in natural DNA were discovered nearly twenty years ago for sequences containing regular repeats of $`\mathrm{A}_\mathrm{n}\mathrm{T}_\mathrm{m},\mathrm{with}\mathrm{n}+\mathrm{m}>3`$, called A-tracts . Since then this intriguing phenomenon has been thoroughly studied, with several profound reviews of the results published in different years . Every A-tract slightly deviates the helical axis towards its minor groove, which results in significant macroscopic curvature when they are repeated in phase with the helical screw. However, in spite of considerable efforts spent in attempts to clarify the the physical origin of this phenomenon it still remains a matter of debate because all theoretical models proposed until now contradict some of the experimental results. This problem is of general importance because the accumulated large volume of apparently paradoxical observations on the subject points to our lack of understanding of the fundamental principles that govern the DNA structure.

A variety of theoretical methods based upon computer molecular modeling have been used in order to get insight in the mechanism of DNA bending . The most valuable are unbiased simulations where sequence dependent effects result only from lower level interactions and are not affected by a priori empirical knowledge of relevant structural features . Such calculations can reveal the essential physical factors and eventually shed light upon the true mechanism of DNA bending. We have recently reported about the first example of stable static curvature emerging spontaneously in free molecular dynamics simulations a B-DNA fragment with A-tracts phased with the helical screw, and proposed a new mechanism of bending that could explain our results as well as experimental data . However, the sequence used in these computational experiments was artificial in the sense that it was designed empirically so as to accelerate the development and maximize the amplitude of bending. It was never studied experimentally, therefore, even though it was similar to canonical DNA benders, one could not exclude that the static bending observed in calculations was of a different origin than that found in experiments. Here we report the results of a similar study of a complementary 25-mer DNA duplex $`\mathrm{AAAATGTCAAAAAATAGGCAAATTT}`$. This fragment is found in the center of the first bent DNA locus investigated experimentally in vitro . It belongs to a minicircle DNA from the kinetoplast body of Leishmania tarentolae and, together with several additional A-tracts, provides planar curvature that apparently serves in vivo to facilitate the loop closure. We have only replaced the 3’-terminal $`\mathrm{A}_6`$ tract in the original fragment  by $`\mathrm{A}_3\mathrm{T}_3`$ because our preliminary empirical observations suggested that 3’-terminal $`\mathrm{A}_\mathrm{n}`$ tracts usually need larger time to converge to characteristic conformations with a narrow minor groove.

We show that two independent long MD trajectories starting from straight conformations corresponding to canonical A and B-DNA forms both converge to statically bent structures with similar bending direction and magnitude, thus giving the first example of a natural DNA fragment where this phenomenon is reproduced in simulations. The results are qualitatively similar to our first report  as regards the kinetics of convergence and comparison with different theoretical models of bending. At the same time, along with convergence of the overall macroscopic DNA from, we find here a remarkably larger than earlier degree of convergence in local backbone and base pair conformations. These results confirm that A-tract induced DNA bending found in calculations corresponds to the experimental phenomenon. They provide additional information concerning its mechanism and the intriguing relationship between the base pair sequence and the DNA structure.

## Results

Two long MD trajectories have been computed for the complementary duplex $`\mathrm{AAAATGTCAAAAAATAGGCAAATTT}`$. The first trajectory referred to below as TJB started from the fiber canonical B-DNA structure and continued to 10 ns. The second trajectory (TJA) started form the fiber canonical A-DNA conformation and continued to 20 ns. The longer duration of TJA was necessary to ensure structural convergence. The computational protocols are detailed in Methods.

Figure 1 shows two series of representative structures from the two trajectories. Each structure is an average over a nanosecond interval, with these intervals equally spaced in time. In TJB the molecule was curved always in the same direction, but both the shape of the helical axis and the magnitude of bending varied. The first two structures shown in Fig. 1a are locally bent between the upper two A-tracts while their lower parts are nearly straight. In contrast, in the last three structures, the planar curvature is smoothly distributed over the entire molecule. In TJA, distinguishable stable bending emerged only after a few nanoseconds, but after the fourth nanosecond all average conformations were smoothly bent. In contrast to TJB, however, the bending direction was not stable, and by comparing the two time series of perpendicular projections one may see that during the first 15 nanoseconds the bending plane slowly rotated. During the final five nanoseconds the overall bending direction was stable. In the last conformation, an S-shaped profile of the helical axis is found in the perpendicular projection, which indicates that there are two slightly misaligned bends located between the three A-tracts. The orientations of the helices in this figure was chosen separately for the two trajectories, and one can notice that the left projection in plate (a) is close to the right one in plate (b), in other words, the final bend directions differed by approximately 90. In TJA, the intersection of the minor groove with the bending plane occurs close to the center of the middle A-tract while in TJB this point is shifted towards its 3’ end. In both cases, however, the narrowed minor grooves of the A-tracts appear at the inside edge of the curved axis.

Figure 2 displays fluctuations of various parameters that characterize the overall bending of the double helix. In both trajectories the rmsd from the canonical B-DNA usually fluctuates between 2 and 4 Å and correlates with the three parameters shown in Fig. 2 that measure the bending magnitude. Note that in TJA there was a short period of strong initial bending which have not been detected in Fig. 1. The most significant difference between the two plates in Fig. 2 is in the dynamics of the bending direction. In TJB, the final orientation of the bend was found early and remained quite stable, although the molecule sometimes straightened giving, simultaneously, a low bending magnitude and large temporary fluctuations of direction. In contrast, during the first 15 nanoseconds of TJA, the bending plane made more than a half turn with respect to the coordinate system bound to the molecule, that is a transition occurred between oppositely bent conformations by passing through a continuous series of bent states. Temporary transitions to the straight state were short-living and always reversed, with the bending resumed in the previous direction. Owing to this turn the bend orientations in TJA and TJB converged considerably although not exactly. The zero direction in Fig. 2 corresponds to bending towards the minor groove at the fifth AT base pair of the middle A-tract. We see that, at the end of TJA, it is shifted from zero by an angle corresponding to rotation in one base pair step, which gives a bend towards the minor groove close to the center of the middle A-tract. In TJB, the direction deviates from zero in the opposite sense by an angle corresponding to roughly two base pair steps, resulting in a residual divergence of approximately 90 between the two trajectories. The slow kinetics of convergence exhibited in Fig. 2 indicates, however, that still better accuracy, if ever possible, would require much longer time.

Figure 3 displays the time evolution of the profile of the minor grooves in TJB and TJA. At the end of both trajectories the minor groove width exhibits modulations phased with the helical screw. It is significantly widened between the A-tracts and narrowed within them by approximately 1 Å with respect to the canonical level. This magnitude of narrowing corresponds well to the values observed in experimental structures of A-tract containing B-DNA oligomers . In TJB, the overall waving profile established early and remained more or less constant. Interestingly, during two rather long periods, a secondary minimum occurred close to the 5’ end of the middle A-tract, and at the same time the main central minimum sifted towards the 3’ end of this A-tract. These motions involve almost an entire helical turn and, apparently, are concerted, which demonstrates the possibility of medium range structural correlations along the double helix. Comparison of Figs. 1a and 3a suggests that there is no simple one-to-one relationship between bending and minor groove modulations. Notably, the right smaller and narrower widening corresponds to a stable and strong bending point of the helical axis, while the left one, which is evidently larger, gives less or no bending. In TJA, the final configuration of the minor groove established only during the last few nanoseconds, but the final profile has the same number and similar positions of local maxima and minima as that in TJB. The overall minor groove dynamics in TJA looks rather complicated and its relationship with the quasi-regular rotation of the bending plane demonstrated in Figs. 1b and 2b is not readily seen.

Figure 4 displays dynamics of $`\mathrm{B}_\mathrm{I}\mathrm{B}_{\mathrm{II}}`$ backbone transitions in the two trajectories. A few common features can be noticed that have been encountered previously . For instance, in A-tracts, the BII conformers are commonly found in ApA steps and almost never in TpT steps. They tend to alternate with BI in consecutive ApA steps. $`\mathrm{B}_\mathrm{I}\mathrm{B}_{\mathrm{II}}`$ transitions often occur concertedly along the same strand as well as in opposite strands. The $`\mathrm{B}_\mathrm{I}\mathrm{B}_{\mathrm{II}}`$ dynamics comprises all time scales presented in these calculations and clearly involves slower motions as well. Note, for instance, a particularly stable BII conformer in the only GpA step available in the two strands. On the other hand, there is some similarity in the distributions of the BII conformers in the two trajectories, which is a new feature compared to our previous report . It is seen in Fig. 4 that in total ten BII conformers were found at the end of TJB and eight in TJA. Among them six and five, respectively, occurred in non A-tract sequences. In three base pair steps the BII conformers are found in both final conformations, with all of them in non A-tract sequences. A careful examination the two plates in Fig. 4 shows that, although in A-tracts the preferred sites of BII conformations differ, in the intervening sequences their dynamics is rather similar in TJB and TJA. This trend is demonstrated in Fig. 5 where inter-trajectory correlations are examined for the specific base pair step propensities to BI and BII conformers. We have not included here the TpT steps in A-tracts because they strongly prefer the BI conformation and, therefore, are trivially correlated. It is evident that, except for ApA steps in A-tracts, there was certain correlation in the average populations of BII conformers for each specific base pair step in the two trajectories. The ApA steps apparently can adopt both conformations with little effect of the sequence context and the overall structure.

Figure 6 shows variation of several helicoidal parameters along the duplex in three representative one nanosecond averaged structures. Two of them were taken from TJA, namely, the 16th and 18th nanosecond averages which we refer to as TJA16 and TJA18. They illustrate the scale and the character of fluctuations of these parameters in the course of the same dynamics. The third conformation is the last average of TJB (TJB10) and it illustrates convergence of helical parameters in independent trajectories. We have chosen TJA16 and TJA18 because, as seen in Fig. 1, the corresponding two structures are particularly similar. They are both smoothly bent in a virtually identical direction and their rmsd is only 0.95 Å. All parameters shown in the figure, except the inclination, exhibit jumping alternations between consecutive base pair steps. Although they look chaotic, there is a considerable similarity between TJA16 and TJA18 and less significant, but still evident similarity of the two with TJB10. Notably, a remarkable correspondence of alterations of the twist is observed in the right-hand half of the sequence. At the same time, even the TJA16 and TJA18 plots sometimes diverge. Note, for instance, that the alteration in their roll traces are phased in the central A-tract, but dephased in the other two, with a particularly significant difference around the TpA step. These results show that, in a statically curved double helix, the base pair helicoidal parameters fluctuate around certain specific values that are stable in the nanosecond time scale. There is, however, more than one combination of these parameters for the same overall bend. At the same time, the evident convergence of the corresponding distributions in TJA and TJB suggests that, at least for this particular base pair sequence, the number of such redundant combinations should not be very large.

### Discussion

The computational experiments described here give the first example of a natural DNA fragment with phased A-tracts which in free unbiased MD simulations adopts spontaneously a statically bent shape. In the analogous earlier calculations the static curvature has been demonstrated for a different A-tract containing sequence constructed artifically and never tested in experiments . The qualitative similarity between these two simulations is evident. Trajectories starting from canonical A and B-DNA forms both travel within the B-DNA family and, in straight, states yield rmsd from B-DNA of around 2 Å. TJA enters the B-DNA family with a significant temporary bending during the first 500 ps. Later it becomes bent in an arbitrary direction and next changes the bend direction by slowly rotating the bending plane. This rotation slows down after 10 ns, with the final orientation much closer to that in TJB than the initial one. In both cases the residual discrepancy was in the range of 60 – 90 . The final minor groove profiles are not identical, although similar for TJA and TJB, as well as the distributions of the $`\mathrm{B}_\mathrm{I}`$ and $`\mathrm{B}_{\mathrm{II}}`$ backbone conformers and base pair helicoidal parameters. The present results, therefore, suggest that the A-tract induced DNA bending observed in calculations here and before corresponds to the experimental phenomenon.

At the same time, there are several clear differences. Notably, the preferred bending direction here is closer the centers of the minor grooves of the A-tracts, whereas the magnitude of bending is somewhat less than in the previous calculations. The bending angle in the average structures shown in Fig. 1 fluctuates between 12 and 25 in TJB and between 7 and 28 in TJA, with the maximal values reached at the end in both cases. The previous artifical sequence was constructed to maximize the bending and it showed the corresponding values beyond 35 . According to the experimental estimates made for “good benders” in an optimal sequence context, the magnitude of bending is around 18 per A-tract , which in our case gives 36 for the overall bend because the principal bending elements are the two intervening zones between the A-tracts. The bends observed here are somewhat below this estimate. However, in experiments, the bending magnitude differs significantly between different A-tract sequences and depends upon many environmental parameters that are not controlled in simulations. One can expect to observe in calculations sequence variations of the bending magnitude that may not exactly follow those in experiments. Therefore, whatever the possible reasons of the apparent discrepancy, the overall correspondence of the computed bending magnitudes to the experimental estimates should be considered as surprisingly good. Yet another difference is a larger than earlier degree of similarity in the profiles of the minor groove, in the distributions of $`\mathrm{B}_{\mathrm{II}}`$ conformers, and helicoidal parameters in trajectories starting from A and B-DNA forms. It was the most surprising observation of our previous report that reasonably good convergence in terms of the macroscopic bent shape of the double helix was not accompanied by the parallel convergence of microscopic conformational parameters. Here the two trajectories manifest clear signs of convergence for base pair step parameters as well as for the backbone conformations. Additional studies are necessary in order to tell if this difference is a sequence specific property or just an occasional effect.

In spite of these differences the results presented here support our previous conclusion concerning the qualitative disagreement of the computed structural dynamics of the double helix and the most popular earlier theories of bending . In Figs. 3, 4, and 6, multiple examples are found of strong variability of local helical parameters in bent substates, which argues against the local interactions and the preferred base pair orientations as the cause of bending. All three A-tracts are characterized by a narrowed minor groove and local minima in the traces of the propeller twist and the inclination of base pairs. Nevertheless, their internal structures are not homogeneous and vary from one base pair step to another. Moreover, the structures of the three are not the same and present another example of an ensemble of conformational substates with a common overall shape. This pattern is qualitatively different from that implied by the junction models of bending . At the same time, the present results are well interpreted in the framework of our model that sees the principal cause of the bending in the backbone stiffness and the postulated excess of its specific length over that in the idealized regular B-DNA double helix . Several non trivial observations support this view.

The first observation is the microheterogeneity of the ensemble of conformations that provide the same bent form of the double helix during the last nanoseconds of both trajectories. Once the backbone have found its preferred waving shape on the surface of the B-DNA cylinder it fixes the bending direction. The thermal motion of bases is allowed, but they are forced to respect this mechanical constraint, giving rise to an ensemble of conformations with different base orientations, but the same bent form of the double helix

The second observation is an always waving minor groove profile which does not change during temporary short-living straightening. The waving profile is the direct consequence of the postulated excess of the specific backbone length over that in the regular B-DNA with even grooves. The main immediate cause of bending is the necessity to compress the backbone in the minor groove widenings if the parallel base pair stacking is to be preserved . The backbone stiffness tends to cause destacking from the minor groove side, which results in bending towards the major groove. Symmetrical destacking is also possible, however, and transitions between various types of stacking perturbations makes possible time variations of the magnitude of bending with a constant backbone profile.

Finally, our model explains well the persistent bends in incorrect directions and the rotation of the bending plane observed in TJA. According to this view the excessive backbone length and its stiffness force the backbone to wander on the surface of the B-DNA cylinder whatever the base pair sequence is. In the dynamics starting from the A-DNA structure the duplex enters the B-DNA family in strongly non-equilibrium conditions, with rapidly changing different energy components. The backbone quickly extends to its preferred specific length taking some waving profile and causing bending in an arbitrary direction. During the subsequent slow evolution it remains always waving, and that is why there is always a preferred bend direction which is not lost during occasional straightening.

### Methods

As in the previous report , the molecular dynamics simulations were carried out with the internal coordinate method (ICMD) . The minimal model of B-DNA was used, with fixed all bond length, rigid bases, fixed valence angles except those centered at sugar atoms and increased effective inertia of hydrogen-only rigid bodies as well as planar sugar angles . The time step was 10 fsec. AMBER94  force field and atom parameters were used with TIP3P water  and no cut off schemes. The heating and equilibration protocols were same as before . The Berendsen algorithm  was applied during the production run, with a relaxation time of 10 ps, to keep the temperature close to 300 K. The coordinates were saved once in 2.5 ps.

The initial conformation for TJB was prepared by vacuum energy minimization starting from the fiber B-DNA model constructed from the published atom coordinates . 375 water molecules were next added by the hydration protocol designed to fill up the minor groove . The initial conformation for TJA was prepared by hydrating the minor groove of the corresponding A-DNA model  without preliminary energy minimization. The necessary number of water molecules was added after equilibration to make it equal in TJA and TJB.

During the runs, after every 200 ps, water positions were checked in order to identify those penetrating into the major groove and those completely separated. These molecules, if found, were removed and next re-introduced in simulations by putting them with zero velocities at random positions around the hydrated duplex, so that they could readily re-join the core system. This time interval was chosen so as to ensure a small enough average number of repositioned molecules which was ca 1.5.

## Appendix

This section contains comments from anonymous referees of a peer-review journal where this paper was been considered for publication, but rejected (see also ).

### A Journal of Molecular Biology

#### 1 First referee

These companion manuscripts describe a series of molecular dynamics trajectories obtained for DNA sequences containing arrangements of oligo dA - oligo dT motifs implicated in intrinsic DNA bending. Unlike previous MD studies of intrinsically bent DNA sequences, these calculations omit explicit consideration of the role of counterions. Because recent crystallographic studies of A-tract-like DNA sequences have attributed intrinsic bending to the localization of counterions in the minor groove, a detailed understanding of the underlying basis of A-tract-dependent bending and its relationship to DNA-counterion interactions would be an important contribution.

Although the MD calculations seem to have been carried out with close attention to detail, both manuscripts suffer from some troubling problems, specifically:

The sequence investigated here is a 25-bp segment of the well-characterized L. tarentolae kinetoplast-DNA bending locus. Two trajectories, TJA and TJB, were computed starting from canonical A-form and B-form structures, respectively. Although the author argues that greater structural convergence between TJA and TJB has taken place in these simulations, there is still a significant disparity concerning the observed bending directions in these two structures. Moreover, the extent of bending in this simulated helix is significantly less than that observed in the previous study, which is unexpected because of out-of-phase placement of the third A tract in the previous sequence. This behavior is not explained and seems difficult to rationalize.