The following is a summary of work on the cube recurrence through 5/12/2002.
These (edited) emails describe groves and show that the cube recurrence with
standard initial conditions:
f(i,j,k) = x(i,j,k) (i+j+k = 1, 0, 1) ;
f(i,j,k) = [f(i1,j,k)f(i,j1,k1) + f(i,j1,k)f(i1,j,k1) +
f(i,j,k1)f(i1,j1,k)] / f(i1,j1,k1) (i+j+k > 1)
produces Laurent polynomials which count groves. We intend to wait on a
completely formal writeup until we have at least a picture of how the same
recurrence works with general initial conditions, but in the meantime, if anyone
finds mistakes or unclear statements in the following proofs, please let me
know.
 Gabriel
Groves of order n are defined as follows: Begin with a triangular lattice, n
vertices on a side, and label the boundary vertices so that two vertices receive
the same label iff they are exchanged by a reflection through a median of the
triangle passing through the vertex of the triangle closest to either of them
(see example below).
a
/ \
bb
/ \ / \
c*c
/ \ / \ / \
d**e
/ \ / \ / \ / \
fdceg
Groves are forests contained in this graph such that each tree contains
precisely the lettered vertices of one letter, and every tree contains at
least one lettered vertex. For example,
a
bb
c* c
\ /
d* * e
/ / /
f d c e g
Groves can be identified with Laurent monomials in the initial conditions
x(i,j,k) as follows. To get the exponents of the face variables x(i,j,k) with
i+j+k=0, draw them in the trangular grid in the natural way.
(0, 0, 0)
(1,1, 0)(1, 0,1)
(2,2, 0)(2,1,1) (2, 0,2)
\ /
(3,3, 0)(3,2,1) (3,1,2) (3, 0,3)
/ / /
(4,4, 0) (4,3,1) (4,2,2) (4,1,3) (4, 0,4)
For the nonlettered vertices, the exponent is (# of incident edges)2. For the
letter vertices other than the three corners (that is, for b, c, d and e), the
exponent is (# of incident edges)1 and for the corners a, i and j, it can be
called (# of incident edges) but is simply always 0.
Thus, we get the exponents
0
00
00 0
\ /
00 1 0
/ / /
0 0 0 0 0
So the contribution from the x(i,j,k) with i+j+k=0 is x(3,1,2)^1.
Now, we have to get the exponents of the f[i,j,k] with i+j+k=+/ 1. These we
write in the faces of the triangular grid. There is no way I can do ASCII art
large enough to show this, but I'll draw the two types of triangles that occur
and the labels they receive
(i,j,k+1)(i,j+1,k)
\ /
\ (i,j,k) /
\ /
\ /
\ /
(i+1,j,k)
(i1,j,k)
/ \
/ \
/ \
/ (i,j,k) \
/ \
(i,j1,k)(i,j,k1)
The exponent of a triangle is 1(# of adjacent edges used) (even if it is on the
boundary.) Thus, we get
*
0
**
0
0 1
** *
0 \1 /
0 0 \ / 0
** * *
1 / 0 / 0 /
1 / 0 / 0 / 0
* * * * *
Thus, we get a contribution of x(2,0,1)^1 x(2,1,2)^(1) x(3,3,1)^(1)
x(4,3,0) from these terms.
Putting it all together, the monomial assosciated to this grove is x(2,0,1)^1
x(2,1,2)^(1) x(3,3,1)^(1) x(4,3,0) x(3,1,2) which, sure enough, occurs
in f(4,0,0). What we want to show is that each cuberecurrence polynomial is
the sum of all monomials corresponding to groves (of the appropriate order).
To convert back from a monomial to a grove, write the exponents of the monomial
inside the "Star of David" graph. We do our previous example, using  to signify
1 due to save space and using no symbol for a 0.
**
/ \
* *
\ /
****
/ \ / \
* * *
\ / \ /
******
/ \ / \1/ \
* * * *
\ / \ /\ /
********
/ \ / \ / \ / \
* * * 1 * *
\ /\ / \ / \ /
**********
/ \1/ \ / \ / \ / \
* * * * * *
\ / \ / \ / \ / \ /
** ****** **
Thus, the numbers that were written on the vertices of the triangular graph are
written on the hexagons of the star of David graph and the numbers that were
written on the faces of the triangular graph are written on the triangles of the
star of David graph. An edge of the triangular graph corresponds to a vertex of
the star of David graph.
For any vertex of the star of David graph, draw the four rays extending from it.
* *
* *
* * * *
* * *
\
* * * * * *
\ 1
* * * *
\ 1
********
\
* * * 1 * *
1 \
* * * * * * * * * *
1 \
* * * * * *
\
* * * * * * * * * *
The sums of the numbers in the four wedges will either be 0, 0, 0 and 1 (as in
this case) or 1, 1, 1 and 0. In the former case, we do not draw the edge
corresponding to this vertex; in the latter case we do.
(David, 5/5)
Next, here's a description of the groveshuffling algorithm, which turns ordern
groves into order(n+1) groves, and a proof that it works.
The basic idea is that we turn upwardpointing triangles in an ordern grove
into downwardpointing triangles in an order(n+1) grove. To see how this
happens, let's first superimpose the vertex grids of the two groves:
+
*
+ +
* *
+ + +
* * *
+ + + +
* * * *
+ + + + +
* * * * *
+ + + + + +
Here the *'s are vertices of the smaller grid, and the +'s are vertices of the
larger grid. We can see that the upwardpointing triangles of *'s are
concentric with the downwardpointing triangles of +'s. Now suppose we have a
grove on the * grid. For ASCIIgraphic purposes, I will mark edges according to
the grid they come from; hopefully this won't be confusing.
+
*
+ +
*******
+ + +
* * *
* * *
+ * + * + * +
* * *
* ******* *
* * *
+ * + * + + * +
* * *
* * * * *
+ + + + + +
The algorithm is as follows: First we "switch" the uppointing * triangles  if
an up * triangle has two edges, remove both of them; if it has no edges, choose
two edges arbitrarily.
+
*
+ +
*******
* *
+ * + * +
* *
*************
* *
+ * + + * +
* 0 *
* * * *
* * * *
+ * + * + * + * +
* * * *
* * * * *
+ + + + + +
Then shift every horizontal edge up, shift every \shaped edge down and left,
and every /shaped edge down and right. Thus, each edge of a * triangle (in
which it must now be the unique edge) is turned into a parallel edge of the
corresponding + triangle. (My drawings probably lose lucidity at this point; if
so, the interested reader is encouraged to work through the example on paper.)
+
*
+++++++
*******
+++++++++++++
+ +
* + * + *
* + + *
+ * + + * +
+ * 0 * +
* + * * + *
* + * + *
+ * + * +++++++ * +
+ * + * + * +
* + * + * + * + *
+ + + +
+ + + + + +
The claim is that we now have a grove on the + grid.
If an initial * (upward) triangle has one edge, the corresponding + triangle
will have one edge, in the same orientation, produced by shifting the given *
edge. If our * triangle has two edges, the corresponding + triangle will have
no edges. If our + triangle has no edges, the corresponding * triangle will not
have any edges after the shift step, so it will be assigned (randomly) two
edges. Translating this into faceexponent terms, this tells us that the up
triangle exponents in the * graph are the negatives of the corresponding down
triangle exponents in the + graph, since a triangle's face exponent is 1 
(# of edges used).
This actually has two important consequences:
(a) The * edges and + edges _never intersect_. In fact, from looking at the
graphs, we can see that any possible intersection would have to involve a * edge
and a + edge on corresponding triangles. The nature of the shuffling algorithm
guarantees that this won't happen: for any two corresponding triangles, either
the * triangle has two edges and the + triangle has none, or vice versa, or each
triangle has one edge, but in the last case, the two edges used are parallel and
so don't intersect.
(b) The total number of + edges used is forced. This follows from the
inversion of the faceexponents. More specifically: The ordern grove (the *
graph) has 3n/2 trees if n is even, (3n1)/2 trees if n is odd, on C(n,2)
vertices. Each facevariable (corresponding to a * triangle) is of the form (1
 # of edges), and summing over all C(n1,2) triangles counts each prospective
edge exactly once, so the sum S of the facevariables is
S = C(n1,2)  (C(n,2)3n/2) = 3n/2  n = n/2 if n is even,
S = C(n1,2)  (C(n,2)(3n1)/2) = (3n1)/2  n = (n1)/2 if n is odd.
The sum of the facevariables in the corresponding order(n+1) graph is then S.
But this is again sum (1  # edges) over all the + triangles, which is (# of +
triangles)  (# of edges), since we can easily check that the only potential
edges not counted in the sum are those along the border of the + grid, which our
algorithm never produces. The number of triangles is still C(n1,2), so the
number of + edges is
C(n1,2) + S = C(n1,2) + n/2 = C(n+1,2)  (3n+2)/2 if n is even,
C(n1,2) + S = C(n1,2) + (n1)/2 = C(n+1,2) + (3n+3)/2 if n is odd.
In particular, if we can show that the + graph is acyclic, then (since it has
C(n+1,2) vertices) it has (3n+2)/2 or (3n+3)/2 trees, respectively, which is the
correct number for a grove of order n+1.
However, observation (a) shows that the + graph is acyclic. Indeed, we can
check that any possible cycle (except for the trivial cycle made of a single
downwardpointing triangle, which can never occur) would have to contain at
least one * vertex in its interior. Because + edges never cross * edges, the +
cycle would then contain an entire * component. But each * component touches
two of the boundaries of the * grid, so the + cycle cannot contain it without
using boundary edges of the + grid. Since our algorithm never uses these
boundary edges, acylicity follows. The number of + components then follows as
well.
Now for the coup de grace: Let's draw rays coming out of the boundary
vertices of the * grid.
*** + ***
*** ***
***
*** ***
*** + + ***
*** ***
** **
*** ***
*** + + + ***
*** ***
** * **
*** + + + + ***
*** ***
** * * **
+ + + + +
*** ***
** * * * **
* * * * *
+ * + * + * + * + * +
* * * * *
* * * * *
* * * * *
When we insert the * grove graph, the plane is divided into "sectors" by virtue
of the way that the boundary * vertices are connected by the trees.
*** + ***
*** ***
***
*** ***
*** + + ***
*** ***
*********
*** ***
*** + + + ***
*** ***
** * **
* * *
*** + * + * + * + ***
*** * * * ***
** ******* **
* * *
+ * + * + + * +
*** * * * ***
** * * * **
* * * * *
+ * + * + * + * + * +
* * * * *
* * * * *
* * * * *
By lemma (a) and the fact that the + boundary edges are never used, we conclude
that each + component is contained inside one sector. However, it is not hard
to check that we have (3n+2)/2 sectors when n is even, (3n+3)/2 when n is odd 
exactly the number of components of the + graph. Therefore, each sector
contains precisely the vertices of one component of the + graph. It is now a
straightforward induction to see that the boundary vertices are matched up in
the desired manner  if the ordern grove has the property that two boundary
vertices are in the same tree iff they have the same label, then it follows that
two boundary vertices of the order(n+1) grid lie in the same sector iff they
have the same label, and we just saw that being in the same sector is equivalent
to being in the same + tree. (The induction step actually has two cases,
depending on the parity of n, but both are straightforward).
At this point, we've shown that the + graph is acyclic, that it has the
correct number of components, and that the boundary vertices are properly
matched  in other words, it is a grove.
(Gabriel, 5/11)
Shuffling backwards works much the same way; the only added detail to the proof
is that, if an edge lies along the boundary, as below, there is no place to
shuffle it to.
*
* *
\
* * *
* * * *
However, such an edge would join two vertices that are supposed to be in
different trees, so it cannot be present.
So we can conclude that every ordern grove shuffles into order(n+1) groves
and, conversely, every order(n+1) grove is obtained by a shuffle from ordern
groves. In other words, the set of objects produced by the shuffling algorithm
is the same as the set of groves, for any given order.
(David, 5/11)
Given any grove of order n, we can extend it to a forest on the entire plane by
adding some edges. Here we have the edges to be added:
* * * * * ****
\ \ \ \ \
* * * * *****
\ \ \ \
* * * * *****
\ \ \ \
* * * * *****
\ \ \
* * * * * ****
\ \ \
* * * * * ****
\ \
* * * * * * ***
\ \
* * * * * * ***
/ / / / / / /
* * * * * * * **
/ / / / / / /
* * * * * * * **
/ / / / / / / /
* * * * * * * * *
Of course, this pattern continues outward to infinity. In case it's not clear,
the procedure is simply to extend each side of the triangle to a ray in one
direction, thus dividing the exterior of the triangle into three regions, and
then to fill each region with edges that are all of the same orientation.
The central triangle is where you fill in the edges of your grove (in the
example above, of order 6). It's straightforward to check that we can now set
the exponent of any triangle to be 1  (# of surrounding edges) and the exponent
of any * to be (# of edges)  2. We get the same values as before for all the
exponents (and those that lie outside the triangle are all zero). It's also
straightforward to check that if we shuffle an extended ordern grove, the
result is an extended order(n+1) grove.
This extension will make subsequent computations a bit easier, since we
don't have to worry about special stuff happening on the boundary. It's
analogous to a procedure used in EKLP: to shuffle a domino tiling of the Aztec
diamond, we first extend the tiling to the entire plane by using horizontals.
(Gabriel, 5/12)
So now let's show how shuffling groves is equivalent to a change of variables in
the grovecounting polynomials. From here on we will assume that groves are
extended to the whole plane as described in my previous email, so we don't need
to worry about boundaries, and we're still assured that only finitely many
exponents are nonzero for any given grove.
I'm going to use the a notational convention David recently established:
the initial conditions are
x(i,j,k) for i + j + k = 1 (the "uppointing triangles")
y(i,j,k) for i + j + k = 0 (the vertices of the grove)
z(i,j,k) for i + j + k = 1 (the "downpointing triangles")
So we denote the grovecounting polynomial at (I,J,K), which is a Laurent
polynomial in infinitely many variables (although only finitely many of them
actually appear), as
f(I,J,K) = P_I,J,K ({x(i,j,k)},{y(i,j,k)},{z(i,j,k)})
Notice that lowercase letters index all the formal variables, and uppercase
letters index the particular polynomial in question.
We want to show that
f(I+1,J,K) = P_I,J,K (
{ ( x(i,j,k)y(i+1,j1,k1) +
x(i+1,j1,k)y(i,j,k1) +
x(i+1,j,k1)y(i,j1,k) ) / z(i,j1,k1) },
{x(i+1,j,k)},
{y(i+1,j,k)} ).
Because the cuberecurrence polynomials satisfy this same recurrence (I'm
not sure how to explain this if it isn't clear  basically, you get this by
treating the planes i+j+k = 0, 1, 2 as your initial conditions instead of 1, 0,
1), and it's straightforward to check that they agree for the intial cases, this
will prove that the cuberecurrence polynomials and the grovecounting
polynomials coincide.
So let's look at the process of shuffling groves. Recall that it operates
as follows:
 First "switch" the up triangles: if a up triangle has two edges, remove
them; if it has no edges, add two arbitrarily.
 Then, "shift" every edge of every up triangle to the parallel edge of the
corresponding down triangle, e.g.
* *
+ \ + + +
\ > \
* * * \ *
+ +
(\*'s are the vertices of the ordern grove, and +'s are the vertices of the
order(n+1) grove, as before.)
Let's see what each of these steps does to the grove polynomial. In the
process, we'll work through an example. Pretend that the following groves have
been extended to the entire plane as described in my previous message.
* 0
0
** 0 0
0
1 0
* ** 0 0 1
\ / 0 
\ / 0 1 0
* * * * 0 0 0 0
Switching can produce the following (or any of eight other possibilities):
* 0
0
** 1 0
/  0
/  0
*** 2 1 1
\ \ /  2
\ \ / 0  0
* ** * 0 1 2 0

Notice that we are using the same rules as before to turn a graph into a
monomial, even though the graph is not a grove at this stage. Also notice that
we may temporarily have a nonzero exponent outside of our triangular grid, but
since we're actually working with infinite grids, this isn't a problem.
When we do the shift, we get:
+ 0
0
++ 0 0
0
0 0
+++ 0 1 0
/  0
/ 1  1
+ ++ + 0 1  0
\ \ / 0  0
\ \ / 0 0 1 0
+ + + + + 0 0 0 0 0
It should be noted that there is a somewhat arbitrary choice in the
correspondence between the vertices of the ordern grid and those of the order
(n+1) grid; we have effectively made this choice by choosing to compare P_I,J,K
with P_I+1,J,K (as opposed to, say, P_I,J+1,K).
How do we express the new exponents in terms of the old ones?
First, let's deal with the switching step. If an up triangle has exponent 1
(no edges), we can replace it by any of the following three configurations:
* * *
/ / \ \
/ / \ \
** * * **
These correspond to changing the 1 to a 1 and performing the following net
changes on the six neighboring exponents:
1 2 1
/ 0 / \ 0 \
/ / \ \
21 1 1 12
 0 
It can be checked in the example above that the switched order4 grove monomial
is obtained from the original by two operations of this sort.
So if a monomial in P_I,J,K has a variable x(i,j,k) with exponent 1, it gets
effectively replaced by the trinomial
B(i,j,k) =
y(i1,j,k)^2 y(i,j1,k) y(i,j,k1) / z(i1,j1,k) z(i1,j,k1) x(i,j,k) +
y(i1,j,k) y(i,j1,k)^2 y(i,j,k1) / z(i1,j1,k) z(i,j1,k1) x(i,j,k) +
y(i1,j,k) y(i,j1,k) y(i,j,k1)^2 / z(i,j1,k1) z(i1,j,k1) x(i,j,k).
We claim that the same substitution describes the effect of a switch on up
triangles with exponent 1. Indeed, if a up triangle has exponent 1 in
P_I,J,K, then we have one of the three twoedge configurations above. We
claim that the monomials obtained by replacing this configuration with either
of the other two also occur in P_I,J,K. To see this, notice that the claim is
equivalent to the assertion if an up triangle has two edges used, we may replace
them by either other pair of edges and still obtain a valid grove. But this is
easy to see.
Therefore, the polynomial consisting of those monomials in P_I,J,K where
x(i,j,k) has exponent 1 is actually divisible by the same trinomial B(i,j,k) as
above. Moreover, we can write B(i,j,k) = C(i,j,k)/x(i,j,k), where C(i,j,k) is a
trinomial independent of the variable x(i,j,k). In particular, the substitution
x(i,j,k) < B(i,j,k) is an involution, so applying the switch move to a triangle
of exponent 1 (which should entail dividing by B(i,j,k) and multiplying by
x(i,j,k)) is represented by this substitution.
Finally, if an uptriangle has exponent 0 (one edge used), it will not be
affected by the substitution, which is exactly what we want from the switching
step. Thus, if we let Q_I,J,K({x(i,j,k)},{y(i,j,k)},{z(i,j,k)}) denote the
polynomial which counts switched (but not yet shifted) groves, we have shown
that Q_I,J,K is equal to P_I,J,K({B(i,j,k)},{y(i,j,k)},{z(i,j,k)}).
Now let's deal with the shifting step. The new down triangles have the same
exponents as the corresponding up triangles. For the + vertices we have the
following picture:
*
a f + +
a f f a
* * f a
b e > +ccc+ddd+
b e b e
*ccc*ddd* b e
+ +
where each distinct lowercase letter denotes an edge that may be present or
absent. If we denote the exponents in the picture on the left (order n) by
uppercase letters as follows:
*
A
* *
X
B C
* * *
then the number of edges used from {a,b,c,d,e,f} is (1A) + (1B) + (1C) 
(1X) = X  (A+B+C) + 2. Therefore, the exponent of the central * (in order
n+1) is X(A+B+C).
For the up triangles, we have the following picture:
* *
a c +
a c a c
* * * > a c
+bbb+
*bbb*
If we again label the three up triangles' and the central vertex's exponents by
uppercase letters:
* *
A C
* X *
B
* *
then the number of edges used from among {a,b,c} is (1A)+(1B)+(1C)(X+2) =
1(A+B+C+X). Therefore, the exponent of the resulting up triangle (in order
n+1) is A+B+C+X.
Now, for each exponent in the ordern picture, let's see which order(n+1)
exponents it contributes to. If an old up triangle had exponent X, then it
contributes X to the corresponding new down triangle exponent and the three
neighboring up triangle exponents, and X to the three neighboring vertices (in
the order(n+1) picture).
The other ordern exponents contribute only to the corresponding order(n+1)
exponents. That is, a down triangle of exponent X in order n contributes X to
the exponent of the corresponding vertex in order n+1 (and doesn't contribute to
any other exponents), and a vertex of exponent X in order n contributes X to the
exponent of the corresponding up triangle in order n+1 (and doesn't contribute
to anything else). All of this can be verified in the example shown above  if
we take each nonzero exponent in the postswitch order4 grove, record its
contributions to the nearby exponents in the order5 graph according to the
rules just stated, and add, we do get the same values that were observed above.
Let's write these rules down algebraically. After performing the shift move
on the alreadyswitched ordern groves, we get the order(n+1) groves, which are
counted by P_I+1,J,K. Therefore,
P_I+1,J,K ({x(i,j,k)},{y(i,j,k)},{z(i,j,k)})
= Q_I,J,K ({ z(i,j1,k1)x(i,j,k)x(i+1,j1,k)x(i+1,j,k1) /
y(i+1,j1,k1)y(i,j1,k)y(i,j,k1) },
{ x(i+1,j,k) },
{ y(i+1,j,k) } )
= P_I,J,K ({ ( x(i,j,k)y(i+1,j1,k1) +
x(i+1,j1,k)y(i,j,k1) +
x(i+1,j,k1)y(i,j1,k) ) / z(i,j1,k1) },
{x(i+1,j,k)},
{y(i+1,j,k)} )
where the second equality comes from our earlier expression for Q_I,J,K. And
we're done.
(Gabriel and David, 5/12)