The velocity Verlet method and variable time steps

Does the velocity Verlet handle variable time steps? I found controversial statements about it.

In the paper Skeel, R. D., "Variable Step Size Destabilizes the Stömer/Leapfrog/Verlet Method", BIT Numerical Mathematics, Vol. 33, 1993, p. 172–175. the author proves that the Verlet methods, including the leapfrog (velocity) Verlet

$$x_{i + 1} = x_i + v_i \Delta t + \frac{1}{2} a_i \Delta t^2$$

$$v_{i + 1} = v_i + \frac{1}{2} (a_i + a_{i + 1}) \Delta t$$

with variable time steps have stability issues.

The Wikipedia article of the Leapfrog method and this stackoverflow answer presents the alternative formulation

$$v_{i + 1/2} = v_i + \frac{1}{2} a_i \Delta t$$

$$x_{i + 1} = x_i + v_{i + 1/2} \Delta t$$

$$v_{i + 1} = v_{i + 1/2} + \frac{1}{2} a_{i + 1} \Delta t$$

of velocity Verlet for variable time steps.

Well, this is weird. The two formulas are mathematically equivalent, the second one is only rearranged for an explicit half-step velocity term $$v_{i + 1/2}$$ - if we substitute this term back, then we get back the original velocity Verlet formula.

Why would the first formula fail while the second one remains correct with variable time steps?

The second formula is just the velocity verlet, and it's correct but if you adapt time steps then it's not symplectic. In a separate answer I describe in quite detail that symplecticness is a global property of the integration, it's not a stepwise property. Because of this, local error estimates and local changes of $$\Delta t$$ do not necessarily preserve the property. In fact, most choices will break this property.
It's not necessarily easy to fix, but I can mention the gist of it. What you need to do is extend your Hamiltonian system so that the independent variable s(t) is a part of the Hamiltonain, and then you can change s in a way that stays on the manifold. This choice of function for s, known as the "step size function", requires appending to your Hamiltonian some condition on the time. The resulting methods are then usually implicit in the step size function, so it's not necessarily simple to do or easy to code, which is why it's still an active area of research and not typically used. Chapter 8 Structure-Preserving Implementation of Hairer's Geometric Numerical Integration is a great resource on the topic, so I'll defer the details to that. (The preview for the chapter might get you hooked).