ROTOR-ROUTER AUTOMATA (Spring, 2002) I've done some work recently that combines my interests in probability, cellular automata, number theory, and strongly convergent games. The following is a status report on this work. In the one-dimensional setting, the newest results are findings of my student Lionel Levine, whose boundedness theorem for the one-dimensional version of the "rotor-router" model implies that these automata are able to generate good approximations to quadratic irrationals. (See section 4.) In the two-dimensional setting, I give both experimental and theoretical evidence for the proposition that the asymptotic growth shape for the growth of the occupied region is a perfect circular disk. (See section 6.) Non-probabilists may wish to skip sections 1 and 2 and the first two paragraphs of section 3; this material is only included for motivation. Number theorists will probably find section 4 the most interesting section. CONTENTS: 1. Diffusion-limited aggregation: external and internal 2. The Diaconis-Fulton model 3. The rotor router model 4. The one-dimensional case 5. The directed case 6. A more general circularity conjecture 7. The cellular automaton version 1. Diffusion-limited aggregation: external and internal Diffusion-limited aggregation (DLA for short) is a (probabilistic) growth process that creates a sequence of finite subsets of the lattice Z^d, which I'll denote by S_1, S_2, ... . Here S_1 is the singleton consisting of the origin, and for all n>1, S_n is equal to S_{n-1} with a single lattice-point x_n adjoined, where x_n is adjacent to some element of S_{n-1}. To choose x_n, one imagines a particle that wanders in "from infinity", and lets x_n be the first site that the particle visits that is adjacent to S_{n-1}. See "Diffusion-Limited Aggregation: A Model for Pattern Formation", http://www.aip.org/pt/vol-53/iss-11/p36.html (Physics Today On-Line). Internal DLA is a variant of DLA in which particles are added at the origin rather than at infinity. Unlike external DLA, which gives rise to feathery dendritic shapes of a conjecturally fractal nature that have so far mostly defied rigorous analysis, internal DLA is mathematically tractable: Lawler, Bramson, and Griffeath ("Internal diffusion limited aggregation", Ann. Prob. 14 (1992), 2117-2140) proved that the set S_n, rescaled by sqrt(n), converges in probability to a disk. 2. The Diaconis-Fulton model Diaconis and Fulton, in an unjustly neglected article ("A growth model, a game, an algebra, Lagrange inversion, and characteristic classes", Commutative algebra and algebraic geometry, II (Italian) (Turin, 1990). Rend. Sem. Mat. Univ. Politec. Torino 49 (1991), no. 1, 95--119 (1993)), define an amusing way of "adding" two finite subsets of Z^d. If the sets S and T are disjoint, their sum is just their union. If the intersection of S and T consists of a single point x, the sum of S and T is S union T together with a random point y obtained by letting a particle execute a random walk in Z^d starting from x until it first lands on a point y not in S union T. More generally, if S and T overlap in a set {x_1,...,x_k}, the sum of S and T is a random set of cardinality |S|+|T|, constructed by letting k particles do random walk starting from the respective points x_1,x_2,...,x_k, terminating in distinct points y_1 not in S union T, y_2 not in S union T union {y_1}, y_3 not in S union T union {y_1,y_2}, ... y_k not in S union T union {y_1,y_2,...,y_{k-1}} and adjoining these k new points. This might not appear to be a well-defined definition of addition, since one might suppose that the probability law governing the resulting set S union T union {y_1,...,y_k} would depend on the choice of ordering of the points x_i in the intersection of S and T. However, it turns out that the ordering of the x_i's does not matter. A good way to see this is to use a "stacks" picture, in the manner of David Wilson's work on loop-erased random walk and random spanning trees (see J. Propp and D. Wilson, "How to get a perfectly random sample from a generic Markov chain and generate a random spanning tree of a directed graph", Journal of Algorithms 27 (1998), 170-217). In the case of the Diaconis-Fulton operation, one associates a stack with each vertex v, and the stack at v consists of an infinite sequences of vertices, each of which is a random neighbor of v; all the elements of all the stacks are independent of one another. Then, to do a random walk from a vertex v, one pops the top element of the stack of v (call it v') and walks to v', then one pops the top element of _that_ stack (call it v") and walks to v", etc., until the stopping-condition is satisfied (that is, the particle arrives at a vertex that is not already occupied). Relative to the randomness embodied in the entries of the stacks, the Diaconis-Fulton operation is deterministic, and one need only prove that the outcome does not depend on the order in which one resolves conflicts (i.e. situations in which two particles start out attempting to occupy the same site). Proving this claim can be reduced to the simpler claim that (S + {x}) + {y} is the same set as (S + {y}) + {x}. In the Diaconis-Fulton framework, internal DLA can be described very succinctly as the result of adding the singleton {0} to itself repeatedly. 3. The rotor-router model Let us focus on the case of the graph Z^2 (with the usual way of viewing the four cardinal directions as North, South, East, and West), and replace the random stacks by deterministic stacks, so that the stack at v consists of the Eastward neighbor of v, the Northward neighbor of v, the Westward neighbor of v, the Southward neighbor of v, the Eastward neighbor of v, the Northward neighbor of v, the Westward neighbor of v, the Southward neighbor of v, ... etc. This turns internal DLA into a deterministic game satisfying the "strong convergence property": if you add a finite number of particles to the system (at different locations), letting each particle find an empty site before the next particle is added, the resulting state does not depend on the order in which the particles were added. (Note that this is an example of a strongly convergent game, of the sort studied by Kimmo Erickson; see e.g. "Strong convergence and the polygon property of 1-player games", Discrete Math. 153 (1996), 105-122.) Here's a description of this evolution rule that doesn't depend on your having read anything that appears above: Particles are added to sites in the grid. Each site can be either empty or occupied. When a particle is added to an empty site, the site becomes occupied. When a particle is added to a site that is already occupied, the new particle is ejected to one of the four neighboring sites. The first time a specific site ejects a particle, ejection occurs to the EAST; the second time, ejection occurs to the NORTH; then to the WEST; then to the SOUTH; then to the EAST again, and so on, rotating ad infinitum. As long as there are only finitely many particles on the grid, the new particle will eventually be routed to an empty site it can call its own (this requires proof but is not hard). Note that each site has its own mod-4 counter, telling how many particles (mod 4) it has already ejected, thereby determining the direction into which the next visiting particle will be routed. We could call these mod-4 numbers "states" rather than "counters", and loosen the restriction that their initial values must be 0. I call this the "rotor-router" model (though these are not the same sort of rotors as have been studied by Serge Troubetzkoy and others). I conjecture that: If one adds n particles at the origin (allowing the system to equilibrate after the addition of each new particle), the set of occupied sites, when rescaled by the square root of n, converges to a disk as n goes to infinity. This guess is suggested by three things: 1) The analogy with internal DLA (defined above). Recall that in the original (stochastic) version, a particle that arrives at an occupied site is ejected in a random direction, with each of the four possible directions being equally likely. Here we rotate among the four possible directions in cyclic fashion. Recall that Lawler et al. proved that the set of occupied sites, rescaled by sqrt(n), converges (in probability) to a disk. 2) The proven result that in one dimension, the deterministic version of internal DLA has the same limiting behavior as the stochastic version (see section 4). 3) The simulation I've done with n going from 1 to 100,000. The blob one sees at the 100,000th step appears to be extremely round! What do I mean by "very round"? Well, one way to measure whether a set of lattice points is round is to measure its radius in a bunch of directions from a putative "center" and see how closely the radii agree. I did this with the blob in sixteen different directions, and the radii were all 179 plus or minus 1. More specifically: Direction of radius Length of Radius ------------------- ----------------- (+1, 0) 179 (+4,+3) 179 (+1,+1) 180 (+3,+4) 179 ( 0,+1) 178 (-3,+4) 178 (-1,+1) 179 (-4,+3) 178 (-1, 0) 178 (-4,-3) 178 (-1,-1) 179 (-3,-4) 178 ( 0,-1) 179 (+3,-4) 179 (+1,-1) 180 (+4,-3) 179 4. The one-dimensional case It is natural to consider an analogous system in one dimension, where the routers flip between pointing east and pointing west. Say we start with all rotors pointing east, with one particle at each site in the interval [-x(0),y(0)] with -x(0) <= 0 <= y(0), and we add one particle at a time, letting the system equilibrate before a new particle is added; let [-x(n),y(n)] denote the set of occupied sites after n new particles have been added to the system. Then, as Harvard undergraduate Lionel Levine proved (as part of a much more general result), y(n)^2 - x(n)^2 is bounded below by a constant, while (y(n)-1)^2 - (x(n)+1)^2 is bounded above by a constant. Why is this good news? Because it corresponds to the continuum limit of IDLA! I'll explain... Imagine we add the particles one at a time, at a rate that is much slower than the rate at which the particles find unoccupied sites. If we let -x denote the location of the leftmost particle, and +y denote the location of the rightmost particle (with both x and y functions of time), then for IDLA we have dx:dy::y:x. This is an easy consequence of the fact that if you do Brownian motion in the interval [a,c] with absorption at a and c, and you start from b in (a,c), then P(absorption at a):P(absorption at c)::(c-b):(b-a). Solving x dx = y dy we get y^2 - x^2 = constant. So the near-constancy of (y(n))^2 - (x(n))^2 that we see in the discrete case mirror the exact constancy we see in the continuous case. Here's an amusing aside on the discrete case, also proved by Levine: let A,B,C be any Pythagorean triple with A^2+B^2=C^2. If we begin with x(0)=0 and y(0)=A, then after n = B+C-A additional particles have been added, we get x(n)=B and y(n)=C, _exactly_. Levine's work has other interesting consequences. Suppose we change the rules so that when a particle that we've added to the system at the origin ends up at a hitherto unoccupied site to the left (resp. right) of the other occupied sites. we add r (resp. s) occupied sites at the left (resp. right) of the growing interval. (Thus the rules described above correspond to the case r=s=1.) If we let -x and y denote the positions of the leftmost and rightmost occupied sites, the heuristics for the continuum limit predict that x/y will converge to sqrt(r/s). My experiments in the case r=2, s=1 strongly suggested that in fact x - y sqrt(r/s) is bounded, and Levine has proved this. I suspect that there is a link between what the system does and the continued fraction for sqrt(r/s), though neither Levine nor I has explored this. (Any of you who wish to look into these matters are very welcome to do so, provided you let us know what you find!) 5. The directed case I also ran a 2-dimensional directed IDLA in which each site can eject particles only in the (+1,0) and (0,+1) directions, and in which you put down lots of particles at (0,0) and see where they go. At large times, the set of occupied sites seems to be the intersection of the first quadrant with an increasingly eccentric ellipse that passes through the origin and whose major axis has slope +1. It is not surprising that the occupied sites stay so close to the line y=x, since directed random walk stays close to this line. However, I see no reason why the set of occupied sites should be elliptical. (Perhaps a look at the corresponding stochastic IDLA model would be enlightening.) 6. A more general circularity conjecture One way to prove the circularity conjecture stated above (for the rotor- router model in two dimensions) would be to find some regularities for the behavior of the model, analogous to (but more complex than) the ones found by Levine for one dimension. However, simulations have not been encouraging. A different approach would be to try to prove something much more general and robust, in which what matters is not the fact that the ejections from each vertex v occur in the order East, North, West, South, etc., but the fact that the ejections from each vertex v occur in each of the four directions in nearly equal numbers. Here is a broadened version of the circularity conjecture: Suppose we replace the rotating routers by some more general routing mechanism, with the property that at any stage, and for any vertex v, the number of particles that have been ejected from v at that stage in the process is the same in each of the four directions, to within an error of plus or minus one (or some other fixed bound). Then I claim that the shape of the occupied region (after the system has equilibrated), rescaled by the square root of the number of particles, converges to a disk. A (related) suggestion is that one could prove the original circularity conjecture by focussing on the function N(v,w) = # of particles ejected from vertex v to vertex w and using the fact that N(v) = sum_w N(v,w) is approximately harmonic away from the origin. One might develop a theory of "approximately harmonic" Z-valued functions and use it to show that circularity holds, at least to within an error of o(sqrt(n)). I suspect (based on empirical evidence) that the actual deviations from circularity are much smaller, but even a crude circularity theorem would be worth having as a first step. 7. The cellular automaton version One can convert the rotor-router into a more standard CA with 24 states. With 24=6x4 states we can keep track of the occupancy level of a site (which can be anything between 0 and 5), as well as the direction of the rotor (EAST, NORTH, WEST, or SOUTH). At each tick of the clock, all the sites that are occupied by more than 1 particle eject all the excess particles (and at the same time may receive some particles from neighboring sites). If at time t no site has more than 5 particles, the same will be true at time t+1, since ejection of particles brings a site's occupancy down to 1 and at most 4 new particles can be received during that same time-step. The direction of the rotor determines which neighbors of a site receive the ejected particles. The 24-state model does not allow for new particles to be added to the system, but we can modify the model to allow for particle creation. My approach is to make the origin a special site, with the property that whenever its occupancy is 1, four new particles get created there (which will be sent to the site's neighbors at the next time-step). The set of occupied sites becomes more and more disk-like as time passes. This growth can be modeled by an ordinary nearest-neighbor-interaction cellular automaton (using the 24 states mentioned above along with extra states for the origin). I've run the CA version of the rotor-router model for just under five million steps, and the blob formed by the occupied sites is extremely well approximated by a disk. Specifically, after the 4,706,517th global update, when the millionth particle has been created at the origin, the blob of occupied sites has radius 548-plus-or-minus-1 in each of sixteen different directions. (I remind you that the blob grows asymmetrically, so the radii in these sixteen different directions are all independent of one another.) Are there other CA rules that do as good a job of making circles? The directions in which I measured the radius were N,S,E,W,NE,NW,SE,SW, and the eight directions (+-.6,+-.8), (+-.8,+-.6). Starting with the E-direction, and proceeding in counterclockwise order, the radii were (+1.0, 0.0) direction: 549 (+0.8,+0.6) direction: 549 (+0.7,+0.7) direction: 549 (+0.6,+0.8) direction: 549 ( 0.0,+1.0) direction: 548 (-0.6,+0.8) direction: 549 (-0.7,+0.7) direction: 549 (-0.8,+0.6) direction: 549 (-1.0, 0.0) direction: 548 (-0.8,-0.6) direction: 547 (-0.7,-0.7) direction: 548 (-0.6,-0.8) direction: 547 ( 0.0,-1.0) direction: 547 (+0.6,-0.8) direction: 547 (+0.7,-0.7) direction: 548 (+0.8,-0.6) direction: 547 The deviations from circularity, with the biggest radii in the (+,+) quadrant and the smallest radii in the (+,-) quadrant (more or less), reproduce on a macroscopic scale the rules of the CA (where particles get ejected from occupied sites in the E, N, W, and S directions, cyclically, with the first ejection occurring in the E direction). I'd be interested to know whether there are other CA models with simple update rules and a simple initial state that are either rigorously proved or plausibly conjectured to have a circular asymptotic shape. (Norm Margolus has sent me some references to articles in the physics literature that describe models of wave propagation and models of diffusion that have this property, such as http://xxx.lanl.gov/abs/nlin.CG/0106038, http://xxx.lanl.gov/abs/comp-gas/9811002, and http://xxx.lanl.gov/abs/comp-gas/9403004, but I haven't had a chance to look at them yet.) Jim Propp There's been a new outgrowth of my work on the rotor-router model, involving the use of rotors to simulate not random aggregation but random walk. Define an "over-random walk" on a graph G as an infinite path on G with the property that there exists a universal bound B such that if v is any vertex of G and w_1 and w_2 are any two of its neighbors, then for every N, the number of steps the walk takes from v to w_1 between time 0 and time N is within B of the number of steps the walk takes from v to w_2 between time 0 and time N. That is to say, a walk (viewed as a sequence of directed edges) is over-random if the directed edges emanating from a given vertex are traversed nearly equally often (where "nearly" means O(1)). Note that true random walk on a finite graph is (with probability 1) NOT over-random in my sense, since for random walk one expects the quantities in question to differ by an unbounded amount as N gets large. (So "over-random" seems like a more apt term than "quasi-random", which I used when I posted an earlier version of this message to the domino forum.) In the case where G is finite, it's easy to show that over-random walk visits each site its fair share of the time asymptotically (where "fair" means "in accordance with the steady-state distribution for random walk). Are there any theorems that imply that over-random walk (in the presence of additional assumptions) must satisfy suitable laws of large numbers when G is _infinite_? Note that additional assumptions are clearly needed, since a walk in an infinite graph that never visits any vertex twice counts as over-random under my definition. As a special case of an over-random walk, one could consider a particle travelling in Z2 (starting from (0,0)) where the particle, upon its first visit to a vertex v, chooses a random neighbor of v as the vertex to visit next, and upon each subsequent visit to v, always moves to whichever neighbor of v is "next" in cylic (say counterclockwise) order. Has anyone studied this kind of walk? I have some startling experimental evidence involving a closely related random walk, much like the one described in the preceding paragraph, but completely deterministic: on its first visit to a vertex v, the particle moves to the neighbor of v that is closest to the origin (with a simple rule for breaking ties). The other difference is that when the particle arrives at (1,1), it is placed back at (0,0) before continuing its walk. Let A(n) be the number of times that the particle visits (1,1) before its nth visit to (0,0). Then I find that A(n)/n tends quite rapidly towards p = pi / 8 (which is the probability that a random walker starting at (0,0) visits (1,1) before returning to (0,0)). Specifically, the discrepancy A(n) minus pn stays stays in the interval [-.22, +3.39] for all n up to 10,000 ! This is much smaller than the discrepancy one would see for true random walk, which would be on the order of the square root of n. (E.g., for n=10,000, the root-mean-square deviation would be about 50: an order of magnitude larger than 3.39.) The interval [-.22, +3.39] seems even more surprisingly narrow when one considers that for n between 1 and 5, the discrepancy A(n) minus pn ranges between 0 and 2.963495. So the discrepancy interval doesn't grow much during the particle's next 9,995 visits to (0,0). It seems that A(n) is trying really hard to stay close to pn, with error on the order of log(n) or maybe even less. Has anyone looked at things like this before? I could try to prove some theorems, but first I'd like to make sure I'm not just rediscovering things people already know. (E.g., do numerical analysts know things about approximately harmonic functions on the grid?) Related questions could be formulated about weaker forms of over-randomness, such as requiring discrepancies to be O(sqrt(n)) rather than O(1), or maybe just o(n). In fact, I suspect that I would find results analogous to the ones I'm looking for in the o(n) case if I looked in the early literature on the foundations of probability (specifically, theorems that say that if a sequence is "von-Mises random" then it satisfies various sorts of laws of large numbers). But the O(1) case seems so silly that maybe nobody has studied it, since (e.g.) an "over-random" sequence of coin flips would have the property that the cumulative numbers of heads and tails stay within a constant bound of each other. This is the sort of axiomatization of randomness that might be devised by a frequentialist at a casino when his dwindling supply of cash and his fifth free drink lead him to revert to primitive belief in a "law of averages" that surely must rescue him at any moment. Any thoughts you may have will be appreciated! Thanks, Jim