An orbit of symmetric polynomials Consider the ring of polynomials $R:=\mathbb{Z}[x_1,x_2,x_3]$. Define the operators $E, I:R\rightarrow R$ by $Ef(x_1,x_2,x_3)=f(x_1-1,x_2,x_3)$ and the identity $If=f$.
Let $\mathcal{L}:R\rightarrow R$ be the operator given by
$$\mathcal{L}f=[(x_1+x_2)(x_1+x_3)E-x_1^2I]f.$$
Let $1$ stand for the constant function $f(x_1,x_2,x_3)=1$ and $\mathcal{L}^2f=\mathcal{L}(\mathcal{L}f)$, etc.

CLAIM. Experiments suggest that $\mathcal{L}^n1$ is always a symmetric polynomial in $R$. Any proof?

EDIT. This has found a resolution (see Pietro Majer's answer).
For example, $\mathcal{L} 1=e_2$ and $\mathcal{L}^21=e_2^2-e_1e_2+e_3$ where $e_1=x_1+x_2+x_3, e_2=x_1x_2+x_1x_3+x_2x_3, e_3=x_1x_2x_3$ are the standard elementary symmetric polynomials.

QUESTIONS. (EDIT) These did not find a definitive answer (apart from Brendan McKay's evidence and argument).
(1) Are there other orbits of symmetric polynomials under $\mathcal{L}$?
(2) Are there other non-trivial operators with similar property over $R$?
(3) What about over rings of many more variables?

 A: Noting Pietro's finite check, I can report that all symmetric polynomials $p(x_1,x_2,x_3)$ up to degree 18 inclusive such that $\mathcal{L}p$ is symmetric are linear combinations of $1,\mathcal{L}1,\mathcal{L}^21,\ldots\,$.
This strongly suggests that there are no others.
An observation that might lead to an elementary proof is that, up to degree 18, all polynomials $p$ such that both $p$ and $\mathcal{L}p$ are symmetric are uniquely determined by the coefficients of the powers of $e_2$ (i.e. the terms in the representation in the base $\{e_1,e_2,e_3\}$ which have the form $c e_2^k$).
A: Suppose $f(x,y,z)$ is symmetric (in the following, symmetric tout court always means "symmetric w.r.to the three variables $(x,y,z)$") . Then  $\mathcal{L}f(x,y,z):=(x+y)(x+z)f(x-1,y,z)-x^2f(x,y,z)$ is already symmetric w.r.to $(y,z)$, so it is symmetric if and only if it is symmetric w.r.to $(x,y)$, that is, after simplifications, if and only if $f$ satisfies the functional equation
$$(x-y)f(x,y,z)-(x+z)f(x-1,y,z)+(y+z)f(x,y-1,z)=0.$$
Claim: if both $f$ and $\mathcal{L}f$ are symmetric, so is  $\mathcal{L}^2f$. In other words, the space of all symmetric solutions to the above functional equation is $\mathcal{L}$-invariant. 
Proof: According to the observation above, in order to prove the claim, we need to check that the following is identically zero:
$$(x-y)\mathcal{L}f(x,y,z)-(x+z)\mathcal{L}f(x-1,y,z)+(y+z)\mathcal{L}f(x,y-1,z) $$
and we exploit the symmetry of $\mathcal{L}f$ writing it
$$(x-y)\mathcal{L}f(x,y,z)-(x+z)\mathcal{L}f(y,x-1,z)+(y+z)\mathcal{L}f(x,y-1,z) $$
namely
$$(x-y)\big[(x+y)(x+z)f(x-1,y,z)-x^2f(x,y,z)\big]$$$$
-(x+z)\big[(y-1+x)(y+z)f(y-1,x-1,z)-y^2f(y,x-1,z)\big]$$$$
+(y+z)\big[(x+y-1)(x+z)f(x-1,y-1,z)-x^2f(x,y-1,z)\big]$$
$$=\big[(x-y) (x+y)(x+z) +(x+z) y^2 \big]f(x-1,y,z)$$$$-(x-y)x^2f(x,y,z) 
-(y+z) x^2f(x,y-1,z) =$$
$$=-x^2\big[(x-y)f(x,y,z)-(x+z)f(x-1,y,z)+(y+z)f(x,y-1,z)\big]$$
which is indeed zero, according to the initial remark, since $\mathcal{L}f$ was assumed symmetric.  
