The comments give one way of understanding, but I figured I could maybe explain the original proof better.
Suppose our plane has normal vector $n$ (which is normal to the plane at every point on it) and an example point on the plane: $p$. These two vectors, one a direction and the other a point, determine any plane.
Then, our plane equation is $$ n\cdot(x-p)=0 $$
that is, every $x$ that satisfies this equation is a member of the plane. This is equivalent to $n\cdot x = n\cdot p$.
Notice that $x=p$ is a member of the plane. But, $x=p+n$ is not a member, nor is $x=\vec{0}$, assuming $p\ne \vec{0}$.
Another way to write this is as a point set $\Pi$:
$$
\Pi = \{ \;x\in\mathbb{R}^3\;|\;n\cdot x= n\cdot p\;\}
$$
i.e. these are the set of points making up the plane. Again, notice that if $p=\vec{0}$, then the origin is a point in $\Pi$. But if $p$ does not vanish, then the origin is not in the plane.
The reason is that $p$ is a translation or shifting parameter. That is, a plane has an orientation parameter, $n$, which "rotates" it, and a position parameter $p$, which slides the plane around. When $p=\vec{0}$, we have slid the plane so that it intersects that origin. The plane equation in this case is $n\cdot x = 0$.
Here's a different approach. Every plane is determined by giving 3 unique points. Let's take $p,a,b$. Define $T_1=a-p$ and $T_2=b-p$. We can suppose $T_1$ and $T_2$ are orthogonal; if they are not, we can use Gram-Schmidt orthonormalization.
Now, suppose we walk around between $p$ to $a$ or $b$. This is the same as adding some multiple of $T_1$ or $T_2$ to $p$. So a parametric equation for the plane is
$$
x(s,t) = p + sT_1 + tT_2
$$
so that if you input any $s$ and $t$, your output is a point on the plane.
See also here.
Notice that a normal vector is simply $n=(T_1\times T_2)/||T_1\times T_2||_2$.
So we get equivalence to $n\cdot(x-p)=0$ as before.
One possible proof could be;
Since reflections and rotations are all orthogonal, reflections have determinant $-1$, and rotations have determinant $1$, then the product of two reflections is the product of two orthogonal matrices, hence it is orthogonal, and since
$$\det(AB) = \det(A)\det(B)$$
then the determinant will be $(-1)(-1) = 1$, so the product of two reflections is an orthogonal matrix with determinant $1$, therefore it is a rotation.
However I would say that your proof of multiplying the matrices together and using trig identites certainly counts as a "non-geometric" proof. It just requires algebraic manipulation, not any geometric arguments.
Best Answer
Let's assume the particle passes the origin at some iteration. It follows a straight line and hits the sphere. Notice that a straight line passing the origin and the sphere is a radius, and therefore normal to the sphere. Because of the normality to the surface, the particle returns in the same straight line (this can be proven by symmetry arguments)- and it must pass the origin again.
We conclude that if at some iteration we pass through the origin, we are forever confined to the same line, passing the origin at every iteration after the initial one.
If we use the same logic in reverse-time we conclude that it passed the origin at every iteration before the initial one.
We are told that the particle did not pass the origin at the first iteration, so we can conclude that it never will.