Why do some linear cellular automata over $Z_{2}$ on the torus have small order? At https://dmishin.github.io/js-revca/index.html, you can play around with reversible cellular automata. I noticed that on that site, that for the reversible linear cellular automata (which I have tested), there is always a small natural number $p$ ($p\leq 12$ in the cases which I have tested) where if the grid is a torus whose dimensions are $2^{n}\times 2^{n}$, then when one runs the cellular automaton $p\cdot 2^{n}$ times, the resulting pattern is the identity function. In other words, if $A$ is the matrix for the linear transformation computing the cellular automaton, then $A^{p\cdot 2^{n}}=I_{2^{2n}}$ but $A^{r}\neq I_{2^{2n}}$ whenever $0<r<p\cdot 2^{n}$.
What is the mathematical explanation of this phenomenon? How pervasive is this phenomena of low periods?
Many of these rules have some sort of Sierpinski triangle pattern, but some have more complicated patterns.
For example, for the rule 0,12,10,6,3,15,9,5,11,7,1,13,8,4,2,14, we apparently have p=7.
 A: Since nobody's posted a full answer yet, let me, for the record, provide my partial answer below.  I originally started writing this draft confident that I had the answer, since I was quite familiar with similar phenomena in ordinary linear two-state CA, but then I hit a snag right at the end.  Maybe you or someone else can pick up the argument from where I got stuck, though, or at least find something of use in it.

If I understand you correctly, you're looking at a subset of two-state Margolus-type block cellular automata that are both reversible (in the sense that the block transition map is bijective) and linear (in the sense that, with the block state space viewed as a vector space over F2, the transition maps are linear maps).
Let's first look at the behavior of such CA on a 2×2 toroidal lattice (that being the smallest lattice where the Margolus neighborhood makes sense).  In this case, the entire lattice consist of a single block.  Here it's convenient to assume that, instead of the block boundaries alternating between time steps, the entire lattice configuration is shifted one step diagonally after each time step — which, on a 2×2 cell toroidal lattice, is equivalent to rotating it by 180°.
In this way, we can directly read the time evolution of any configuration of cells on the 2×2 lattice from the block transition table, if we just rotate all the outputs configurations by 180°.  For example, for your example automaton on a 2×2 lattice, the single-cell initial configurations evolve as follows:
$$
\substack{\blacksquare \square \\ \square \square} \mapsto
\substack{\blacksquare \blacksquare \\ \square \square} \mapsto
\substack{\square \blacksquare \\ \blacksquare \square} \mapsto
\substack{\blacksquare \square \\ \square \blacksquare} \mapsto
\substack{\square \blacksquare \\ \blacksquare \blacksquare} \mapsto
\substack{\square \square \\ \blacksquare \square} \mapsto
\substack{\square \square \\ \blacksquare \blacksquare} \mapsto
\substack{\blacksquare \square \\ \square \square}
\\
\substack{\square \blacksquare \\ \square \square} \mapsto
\substack{\blacksquare \square \\ \blacksquare \square} \mapsto
\substack{\blacksquare \blacksquare \\ \blacksquare \blacksquare} \mapsto
\substack{\blacksquare \blacksquare \\ \blacksquare \square} \mapsto
\substack{\square \blacksquare \\ \square \blacksquare} \mapsto
\substack{\square \square \\ \square \blacksquare} \mapsto
\substack{\blacksquare \square \\ \blacksquare \blacksquare} \mapsto
\substack{\square \blacksquare \\ \square \square}
\\
\substack{\square \square \\ \blacksquare \square} \mapsto
\substack{\square \square \\ \blacksquare \blacksquare} \mapsto
\substack{\blacksquare \square \\ \square \square} \mapsto
\substack{\blacksquare \blacksquare \\ \square \square} \mapsto
\substack{\square \blacksquare \\ \blacksquare \square} \mapsto
\substack{\blacksquare \square \\ \square \blacksquare} \mapsto
\substack{\square \blacksquare \\ \blacksquare \blacksquare} \mapsto
\substack{\square \square \\ \blacksquare \square}
\\
\substack{\square \square \\ \square \blacksquare} \mapsto
\substack{\blacksquare \square \\ \blacksquare \blacksquare} \mapsto
\substack{\square \blacksquare \\ \square \square} \mapsto
\substack{\blacksquare \square \\ \blacksquare \square} \mapsto
\substack{\blacksquare \blacksquare \\ \blacksquare \blacksquare} \mapsto
\substack{\blacksquare \blacksquare \\ \blacksquare \square} \mapsto
\substack{\square \blacksquare \\ \square \blacksquare} \mapsto
\substack{\square \square \\ \square \blacksquare}
$$
Since your automata are reversible, all patterns must eventually cycle back to their original state, as seen above.  Since your automata are also linear, adding together any two patterns (modulo 2) and letting the resulting pattern evolve simply yields the sum (modulo 2) of the corresponding time evolutions of the original patterns.
As a result, the time evolution of any pattern must return to its original state after $p_1$ steps, where $p_1$ is the least common multiple of the periods of the single-cell patterns.  In particular, for your example automaton this period is 7 steps, as shown by the diagram above.
(Note that, if the 2×2 lattice period $p_1$ obtained by the method above is odd, then using the standard Margolus block division scheme with alternating block grids will cause the pattern after $p_1$ steps to be flipped 180° compared to the original pattern, so that $2p_1$ steps are actually needed for the original pattern to reappear exactly.  This should not affect the further analysis of 4×4 and larger grids in any substantial way, so I'll ignore that particular detail below.)

Now, given that the time evolution of the automaton on a $2^n \times 2^n$ lattice recurs after $p_n$ steps (which it of course, due to reversibility, must do for some period $p_n$), let's consider the same automaton on a $2^{n+1} \times 2^{n+1}$ lattice.  Here, it's natural and useful to divide the $2^{n+1} \times 2^{n+1}$ lattice into four $2^n \times 2^n$ quadrants, like this:
\begin{array}{|c|c|}
\hline
A & B \\
\hline
C & D \\
\hline
\end{array}
and to look at the evolution of those quadrants over $p_n$ steps.
Since the automata we're considering are all linear, we can limit ourselves to considering the evolution of patterns that are initially confined to, say, the upper left quadrant $A$; the evolution of more complicated initial patterns follows from this by linearity and shift invariance.  (Indeed, we could even restrict ourselves to only considering initial patterns consisting of a single live cell in one of the four upper left corner cells, but for the argument given below we don't really need to be that restrictive.)
Let the initial pattern in the upper left $2^n \times 2^n$ quadrant, represented as a $2^{n^2}$-element vector over $\mathbf F_2$, be $\vec x$.  By linearity, the patterns in the four quadrants after $p_n$ time steps will then corespond to the vectors $\vec y_A = \mathbf A \vec x$, $\vec y_B = \mathbf B \vec x$, $\vec y_C = \mathbf C \vec x$ and $\vec y_D = \mathbf D \vec x$, where $\mathbf A$, $\mathbf B$, $\mathbf C$ and $\mathbf D$ are some $2^{n^2} \times 2^{n^2}$ matrices that together describe the time evolution of the automaton over $p_n$ steps on the $2^{n+1} \times 2^{n+1}$ cell lattice.  Since we know that the same automaton on a $2^n \times 2$ lattice is periodic with period $p_n$, it follows from linearity that $\vec y_A + \vec y_B + \vec y_C + \vec y_D = \vec x$.  And since this holds for any initial configuration $\vec x$, we in fact have:
$$\mathbf A + \mathbf B + \mathbf C + \mathbf D = \mathbf I,$$
where $\mathbf I$ is the $2^{n^2} \times 2^{n^2}$ identity matrix.
Given this, let's run the automaton another $p_n$ steps further to obtain the quadrant configurations:
\begin{aligned}
\vec z_A &= \mathbf A \vec y_A + \mathbf B \vec y_B + \mathbf C \vec y_C + \mathbf D \vec y_D, \\
\vec z_B &= \mathbf B \vec y_A + \mathbf A \vec y_B + \mathbf D \vec y_C + \mathbf C \vec y_D, \\
\vec z_C &= \mathbf C \vec y_A + \mathbf D \vec y_B + \mathbf A \vec y_C + \mathbf B \vec y_D, \\
\vec z_D &= \mathbf D \vec y_A + \mathbf C \vec y_B + \mathbf B \vec y_C + \mathbf A \vec y_D.
\end{aligned}
By substituting in the expressions for $\vec y_A$, $\vec y_B$, $\vec y_C$ and $\vec y_D$ given above, we get:
\begin{aligned}
\vec z_A &= (\mathbf A^2 + \mathbf B^2 + \mathbf C^2 + \mathbf D^2) \vec x, \\
\vec z_B &= (\mathbf B \mathbf A + \mathbf A \mathbf B + \mathbf D \mathbf C + \mathbf C \mathbf D) \vec x, \\
\vec z_C &= (\mathbf C \mathbf A + \mathbf D \mathbf B + \mathbf A \mathbf C + \mathbf B \mathbf D) \vec x, \\
\vec z_D &= (\mathbf D \mathbf A + \mathbf C \mathbf B + \mathbf B \mathbf C + \mathbf A \mathbf D) \vec x.
\end{aligned}
Now, here's the missing part: if the matrices $\mathbf A$, $\mathbf B$, $\mathbf C$ and $\mathbf D$ all commute with each other, then (given that $1+1=0$ in $\mathbf F_2$) the last three expressions all simplify to $\vec z_B = \vec z_C = \vec z_D = 0$.  Meanwhile, the first expression can be rewritten using the same properties as:
\begin{aligned}
\vec z_A &= (\mathbf A^2 + \mathbf B^2 + \mathbf C^2 + \mathbf D^2) \vec x \\
 &= (\mathbf A + \mathbf B + \mathbf C + \mathbf D)^2 \vec x & (!) \\
 &= I^2 \vec x = \vec x.
\end{aligned}
(Yes, the square of the sum indeed does equal the sum of the squares here, as can be confirmed by writing it out and observing that all the cross terms cancel out due to the assumed commutativity and cancellation properties noted above.)

Unfortunately, I have not been able to show that the matrices $\mathbf A$, $\mathbf B$, $\mathbf C$ and $\mathbf D$ indeed do commute.  If we were dealing with an ordinary, fully shift-invariant CA over $\mathbf F_2$, then I believe I could show this by appealing to the shift-invariance and the basic fact that addition of cell coordinates modulo $2^{n+1}$ commutes.  But the Margolus lattice is only invariant with respect to shifting by an even number of cells, and that messes things up.
In particular, for a fully shift-invariant CA, knowing that a cell at position $x$ gives rise $p_n$ steps later to a cell at $y = y+\delta_1$, which in turn gives rise $p_n$ steps later to a cell at $z = y+\delta_2 = x+\delta_1+\delta_2$, would imply that the same cell at $x$ must also give rise to a cell at $y' = x+\delta_2$, which must then give rise to a cell at $z = y'+\delta_1 = x+\delta_2+\delta_1$ (and therefore, by linearity, cancel it out, except in the special case where $\delta_1 = \delta_2$ and therefore $y = y'$).  But this argument doesn't work as stated on the Margolus grid, since it's quite possible for $x$ and $y$ to have different parity, and therefore obey different transition laws.
Given that, experimentally, the property you've observed does hold, it seems that I must be missing something essential, but I can't figure out exactly what.  Perhaps an explicit appeal to reversibility, in some form, might be needed here.  Or something.
A: So I have a proof of the low period phenomenon that works for all cellular automata prime characteristics $p$ and all finite dimensions (i.e. for linear cellular automata  whose alphabet is a vector space of characteristic $p$ over a finite $p$-group). The general result can quite easily be verified empirically. These cellular automata could be applied to the reversible symmetric cryptosystems of the future since functions $f$ where $f^{m}=\textrm{Id}$ for some $m$ which are easily computable by a reversible computer will be important components of cryptographic hash functions of the future.
Let $G$ be a group and $A$ be a set. Then define $\phi_{g}:A^{G}\rightarrow A^{G}$ by letting $\phi_{g}(x_{h})_{h\in G}=(x_{g+h})_{h\in G}$.
A cellular automaton over the group $G$ and the alphabet $A$ is a map $\tau:A^{G}\rightarrow A^{G}$ such that
$\tau\phi_{g}=\phi_{g}\tau$ for all $g\in G$.
If $G$ is a group and $V$ is a vector space over a finite field $K$, then let $\textrm{LCA}(G;V)$ denote the collection of all linear cellular automata over the group $G$ with alphabet 
$V$. Let $c_{a}=(a)_{g\in G}\in V^{G}$ for each $a\in A$. If $V$ is a vector space over a field $K$ and $G$ is a group, the define a homomorphism
$\Gamma:\textrm{LCA}(G;V)\rightarrow\textrm{End}_{K}(V)$ by letting $c_{\Gamma(\tau)(a)}=\tau(c_{a})$.
Theorem: Let $K$ be a finite field of prime characteristic $p$ and let $V$ be a vector space over $K$. Let $G=Z_{p^{a_{1}}}\times Z_{p^{a_{r}}}$. Let $m=p^{a_{1}}+...+p^{a_{r}}$. Let $n'$ be a natural number such that
$p^{n'}\geq m$ and let $n=p^{n'}$.


*

*Suppose that $r\in \textrm{LCA}(G;V)$ and $\Gamma(r)=0$. Then $r^{m}=0$.

*Suppose that $r\in \textrm{LCA}(G;V)$ and $\Gamma(r)=1$. Then $r^{n}=1$.

*Suppose that $r\in\textrm{LCA}(G;V)$ and $\Gamma(r)^{k}=0$. Then $r^{km}=0$.

*Suppose that $r\in\textrm{LCA}(G;V)$ and $\Gamma(r)^{k}=1$. Then $r^{kn}=1$.
Proof: Statement 2 follows from statement 1; if $\Gamma(r)=1$, then $\Gamma(r-1)=0$, so $r^{n}-1=(r-1)^{n}=0$, hence $r^{n}=1$.
Statements 2-4 therefore all follow from statement $1$, so it suffices to prove statement 1.
Suppose that $\Gamma(r)=0$. Now, suppose that
 $$r(x_{\alpha})_{\alpha\in G}=(\sum_{\beta\in G}R_{\beta}x_{\alpha+\beta})_{\alpha\in G}.$$
Then since $\Gamma(r)=0$, we know that $\sum_{\beta\in G}R_{\beta}=0$, so $$R_{0}=-\sum_{\beta\in G^{+}}R_{\beta}.$$ Therefore,
$$r^{m}(x_{\alpha})_{\alpha\in G}=(\sum_{\beta_{1},...,\beta_{m}\in G^{+}}\sum_{u_{1},...,u_{m}\in\{0,1\}^{m}}(-1)^{m+u_{1}+...u_{m}}R_{\beta}x_{\alpha+u_{1}\beta_{1}+...u_{m}\beta_{m}})_{\alpha\in G}$$
$$=(\sum_{\beta_{1},...,\beta_{m}\in G^{+}}R_{\beta}\sum_{u_{1},...,u_{m}\in\{0,1\}^{m}}(-1)^{m+u_{1}+...u_{m}}x_{\alpha+u_{1}\beta_{1}+...u_{m}\beta_{m}})_{\alpha\in G}.$$
It therefore suffices to show that 
$$\sum_{u_{1},...,u_{m}\in\{0,1\}^{m}}(-1)^{m+u_{1}+...u_{m}}x_{\alpha+u_{1}\beta_{1}+...u_{m}\beta_{m}}=0$$
for all $\alpha,\beta_{1},...,\beta_{m}\in G$.
To verify this claim, we shall verify that $$\sum_{u_{1},...,u_{m}\in\{0,1\}^{m}}(-1)^{m+u_{1}+...u_{m}}\cdot(\alpha+u_{1}\beta_{1}+...u_{m}\beta_{m})=0$$ where the element $$(-1)^{m+u_{1}+...u_{m}}\cdot(\alpha+u_{1}\beta_{1}+...u_{m}\beta_{m})$$
is considered as an object in the group ring $K[G]$, Take note that
$$\sum_{u_{1},...,u_{m}\in\{0,1\}^{m}}(-1)^{m+u_{1}+...u_{m}}\cdot(\alpha+u_{1}\beta_{1}+...u_{m}\beta_{m})=\alpha\cdot(\beta_{1}-1)\cdot...\cdot(\beta_{m}-1).$$
Now, take note that the group-ring $K[G]$ is isomorphic to the ring $$K[x_{1},...,x_{r}]/\langle x_{1}^{p^{a_{1}}}-1,...,x_{r}^{p^{a_{r}}}-1\rangle,$$ and by this isomorphism, $\alpha\cdot(\beta_{1}-1)\cdot...\cdot(\beta_{m}-1)$ is sent to
$$\mathbf{x}^{\alpha}\cdot(\mathbf{x}^{\beta_{1}}-1)\dots(\mathbf{x}^{\beta_{m}}-1)+\langle x_{1}^{p^{a_{1}}}-1,...,x_{r}^{p^{a_{r}}}-1\rangle.$$
 Now, 
 $$K[x_{1},...,x_{r}]/\langle x_{1}^{p^{a_{1}}}-1,...,x_{r}^{p^{a_{r}}}-1\rangle=
K[x_{1},...,x_{r}]/\langle (x_{1}-1)^{p^{a_{1}}},...,(x_{r}-1)^{p^{a_{r}}}\rangle$$ and the automorphism $\Phi:K[x_{1},...,x_{r}]\rightarrow K[x_{1},...,x_{r}]$ defined by $\Phi(x_{i})=x_{i}+1$ induces an isomorphism on the quotient rings $$\Phi^{*}:K[x_{1},...,x_{r}]/\langle x_{1}^{p^{a_{1}}}-1,...,x_{r}^{p^{a_{r}}}-1\rangle
\rightarrow K[x_{1},...,x_{r}]/\langle x_{1}^{p^{a_{1}}},...,x_{r}^{p^{a_{r}}}\rangle$$
and
$$\Phi^{*}(\mathbf{x}^{\alpha}\cdot(\mathbf{x}^{\beta_{1}}-1)\dots(\mathbf{x}^{\beta_{m}}-1)+\langle x_{1}^{p^{a_{1}}}-1,...,x_{r}^{p^{a_{r}}}-1\rangle)$$
$$=p(\mathbf{x})f_{1}(\mathbf{x})\dots f_{m}(\mathbf{x})+\langle x_{1}^{p_{a_{1}}},\dots,x_{r}^{p_{a_{r}}}\rangle$$
for some polynomials
$f_{1},...,f_{m}$ with $f_{1}(\mathbf{0})=\dots=f_{m}(\mathbf{0})=0$. Since $$f_{1}(\mathbf{0})=\dots=f_{m}(\mathbf{0})=0,$$
we conclude that
$f_{i}=\sum_{\beta\neq 0}a_{i,\beta}\mathbf{x}^{\beta}$. However,
$$f_{1}(\mathbf{x})\dots f_{m}(\mathbf{x})=\sum_{\beta_{1}\neq 0,...,\beta_{m}\neq 0}a_{1,\beta_{1}}...a_{m,\beta_{m}}\mathbf{x}^{\beta_{1}+...+\beta_{m}}\in
\langle x_{1}^{p_{a_{1}}},\dots,x_{r}^{p_{a_{r}}}\rangle.$$
Therefore,
$$p(\mathbf{x})f_{1}(\mathbf{x})\dots f_{m}(\mathbf{x})+\langle x_{1}^{p_{a_{1}}},\dots,x_{r}^{p_{a_{r}}}\rangle=0.$$
We therefore conclude that
$r^{m}(x_{\alpha})_{\alpha\in G}=0$. $\mathbf{QED}$
