What does the time-reversibility of Verlet or any other integration method mean? The wikipedia article about it is very complex, unclear and confusing. And how can I determine whether a method is time reversible or not?
For example the classical Störmer-Verlet method
$ x_{n+1} = 2x_n - x_{n-1} + a(x_n) · dt^2 $
How is this time-reversible? If I change the sign of the timestep dt, because of the square nothing changes.
Yes, it is exactly that.
The time reverse of the explicit Euler method is the implicit Euler method. $y_{n+1}=y_n+f(y_n)\,dt$ gets reversed to $y_{n-1}=y_n+f(y_n)\,(-dt)$ and after index shift $y_{n+1}=y_n+f(y_{n+1})\,dt$.
The same for the symplectic Euler methods. \begin{array}{lll} forward:&x_{n+1}=x_n+v_n\,dt,& v_{n+1}=v_n+a(x_{n+1})\\ reverse:&x_{n-1}=x_n-v_n\,dt,& v_{n-1}=v_n+a(x_{n-1})\\ shifted:&v_{n+1}=v_n+a(x_n)\,dt,& x_{n+1}=x_n+v_{n+1}\,dt \end{array} (Velocity) Verlet is a combination of an explicit and implicit symplectic Euler step, thus invariant under time reversal. \begin{align} v_{n+1/2}&=v_n+a(x_{n})\,dt/2\\ x_{n+1/2}&=x_n+v_{n+1/2}\,dt/2\\ x_n&=x_{n+1/2}+v_{n+1/2}\,dt/2\\ v_n&=v_{n+1/2}+a(x_{n+1})\,dt/2 \end{align} Elimination of $x_{n+1/2}$ gives velocity Verlet, elimination of the integer-indexed velocities gives the Leapfrog method, elimination of all velocities gives the basic Stoermer-Verlet method.
The first reversion-invariant Runge-Kutta methods are the (necessarily implicit) trapezoidal and midpoint methods. \begin{align} &trapezoid:&y_{n+1}&=y_n+\tfrac12(f(y_n)+f(y_{n+1}))\,dt\\ &midpoint:& y_{n+1}&=y_n+f(\tfrac12(y_n)+y_{n+1}))\,dt \end{align}
The advantage of reversion-invariant methods is that the error in scalar invariants is also reversion invariant and thus an even function. For the second-order Verlet method one would expect from the order alone that the first term of the local error in energy and momentum is $O(dt^3)$, but since the odd powers are missing the actual error is $O(dt^4)$.
Of course, the same or better can be achieved with the classical RK4, but symplectic Euler and Verlet are better suited for real-time simulations and provide, because of "symplectic", global (very slowly eroding) error bounds on conserved quantities of the dynamical system.