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, unlessexcept 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.
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, unless $\delta_1 = \delta_2$!). 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.
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.
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 $n$$2^{n+1}$ commutes. But the Margolus lattice is not fully shift-invariantonly 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, unless $\delta_1 = \delta_2$!). 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.
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 $n$ commutes. But the Margolus lattice is not fully shift-invariant, 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, unless $\delta_1 = \delta_2$!). 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.
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, unless $\delta_1 = \delta_2$!). 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.
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 $n$ commutes. But the Margolus lattice is not fully shift-invariant, 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, unless $\delta_1 = \delta_2$!). 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.