{'fraction_non_alphanumeric': 0.032138000143153676, 'fraction_numerical': 0.009161835230119533, 'mean_word_length': 4.686609686609686, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Exact evolution of time-reversible symplectic integrators and their phase error for the harmonic oscillator

The evolution of any factorized time-reversible symplectic integrators, when applied to the harmonic oscillator, can be exactly solved in a closed form. The resulting modified Hamiltonians demonstrate the convergence of the Lie series expansions. They are also less distorted than modified Hamiltonian of non-reversible algorithms. The analytical form for the modified angular frequency can be used to assess the phase error of any time-reversible algorithm.


I. INTRODUCTION
Symplectic integrators 1 are methods of choice for solving diverse physical problems ranging from celestial mechanics 2 , molecular dynamics 3 , to accelerator physics 4 and lattice gauge theory 5 . In contrast to other numerical methods for solving Hamiltonian dynamics, symplectic integrators evolves the system according to a modified Hamiltonian 2,6 that can be made arbitrary close to the original Hamiltonian. For this reason, despite the fact that symplectic integrators can preserve all Poincaré invariants, they can never exactly conserve energy 7 . Otherwise, an integrator's modified Hamiltonian would have coincided with the original Hamiltonian, and the integrator's evolution would have been exact.
Symplectic integrators can be derived most easily by the method of factorization/composition 2,8 , where the algorithm is composed of many elementary ones with varying coefficients. In this case, one can systematically compute the modified Hamiltonian order-by-order via the Lie series 6 approach by using the Baker-Campbell-Hausdorff (BCH) expansion 2,6 . While this approach can demonstrate the existence of the modified Hamiltonian perturbatively, and show that it can be made arbitrary close to the original Hamiltonian, one has no sense of the modified Hamiltonian's global structure nor of the series' convergence.
In this work, we show that the evolution of any factorized, time-reversible integrators, can be exactly solved for the harmonic oscillator. The algorithm's evolution remained that of a harmonic oscillator but with modified mass and spring constant. Thus the model illustrates very simply, how reversible symplectic algorithms preserve invariant tori of a dynamical system. Moreover, the closed form expression for the modified Hamiltonian, when expanded, matches the Lie series order by order, thus demonstrating the latter's convergence. Finally, we show that the analytical modified oscillator frequency can be used to benchmark any time-reversible algorithm.
In the following, we will briefly summarize key features of the Lie-Poisson construction of symplectic integrators, the Lie series expansion for the modified Hamiltonian, the phase-space matrix approach, and comparative results for a number of fourth order reversible algorithms.

II. SYMPLECTIC INTEGRATORS
The evolution of any dynamical variable W (q i , p i ), is given by the Poisson bracket, and therefore by the corresponding Lie operatorĤ associated with the Hamiltonian function H(q i , p i ), i.e.
For any dynamical function Q, we can similarly define its associated Lie operatorQ via the Poisson bracket The operator equation (1) can be formally solved via Symplectic algorithms are derived by approximating the evolution operator e tĤ for a short time in a product form. For Hamiltonian function of the usual separable form, the Hamiltonian operator is also separable,Ĥ =T +V , with differential operatorsT andV given bŷ The Lie transforms 6 e εT and e εV , are then displacement operators which shift q i and p i forward in time via q → q + εp and p → p + εF.
Thus, if e εĤ can be factorized into products of Lie transforms e εT and e εV , then each factorization gives rise to an integrator for evolving the system forward in time. Existing literature on symplectic algorithms are concerned with decomposing e εĤ to high orders in the product form with a set of well chosen factorization coefficients {t i , v i }.
In most cases, we will consider only (left-right) symmetric factorization schemes such that either t 1 = 0 and In either cases, the algorithm is exactly timereversible, and the energy error terms can only be an even function of ε. Such a symmetric factorizations is then at least second order.

III. LIE SERIES EXPANSION
A well known symmetric algorithm is the second order Störmer/Verlet (SV) algorithm, whereĤ A is the approximate Hamiltonian operator of the form where e T T V , e V T T T V etc., are coefficients specific to a particular algorithm. We have used the condensed com- These are computed by repeated applications of the Baker-Campbell-Hausdorff (BCH) formula for combining exponentials of operators. Thus the algorithm evolves the system according toĤ A rather than the true Hamil-tonianĤ. KnowingĤ A then allows us to determine the actual Hamiltonian function H A which governs the algorithm's evolution. From the fundamental definition of Lie operator (2), one can convert operators back to functions For symmetric factorization with only even order commutators, the Hamiltonian function can be obtained simply by replacing condensed commutators with condensed Poisson brackets: the non-vanishing Poisson brackets are: There is thus a clear separation between the error coefficients, which are characteristics of the algorithm and the Poisson brackets, which are properties of the underlying Hamiltonian. If the SV algorithm is applied to the harmonic oscillator, then the resulting modified Hamiltonian function is with effective mass and spring constant, and approximate angular frequency ω A = k * m * . The algorithm therefore evolves the system according to withĤ A directly given by the BCH formula, produces the following expansion, Note that the overall factor is the same as (15). The odd order terms distort the harmonic oscillator fundamentally, producing a distinctive, 45 • -tilted phase-space ellipse, commonly seen in the standard map 9 , or when solving the pendulum 2 . This is an artifact of not being time-reversible, a property foreign to the original dynamics. Time-reversible algorithms with modified Hamiltonian such as (14) showed no such distortion 10 .

IV. MATRIX METHOD
The above expansion method for computing the modified Hamiltonian, while general, is severely limited to low orders. At higher order, the number of commutator and Poisson brackets proliferates greatly due to the complexicity of the BCH expansion, it is then difficult to assess the convergence of series (15), (16) and (20).
To compute H A exactly, we completely abandon the BCH-based approach. For the 1-D harmonic oscillator, let's denote the phase-space vector at the ith iteration of the algorithm as For the SV algorithm, let T and V denote the effect of Lie operators e εT and e 1 2 εV acting on these vectors. From (7), it is easy to see that they are upper and lower triangular matrices given by Hence, the algorithm corresponds to the product matrix, (23) Matrix representations have been widely used to study the stability 11,12 and convergence 3,13 of low order symmetric integrators. Our interest here is not to solve the harmonic oscillator per se, but to use it as a mean of obtaining the converged value of the Lie series. In contrast to previous use of the matrix method to study the harmonic oscillator 3,11,12,13 , we are not interested how each algorithm's approximate solution approaches the exact solution; we are interested only in the exact form of the approximate solution itself. Below, we also emphasize the distinction between time-reversible and time-irreversible algorithms. As illustrated in (22), for quadratic Hamiltonians, each elemental operator e tiεT and e viεVi can in general be represented as upper or lower triangular matrices, where the constant matrix elements σ i (ε) and µ i (ε) are both positive as ε → 0 but are odd functions of ε: Both T i and V i can be further decomposed as where t and v are nilpotent matrices, with property tv + vt = 1l. The nilpotence of t and v makes it obvious that both T i and V i are time-reversible, The time-reversibility of T i (ε) and V i (ε) can also be traced to the fact that they have equal diagonal elements even in ε and off-diagonal elements odd in ε. In fact such a general matrix with τ (−ε) = −τ (ε) and ν(−ε) = −ν(ε), is indeed timereversible, provided that (det G) = 1. Since this is true of any symplectic integrator, we conclude that any symplectic integrator of the form (30), is time-reversible. The generalization of (23), corresponding to the symmetric form of (8), can be now written as The advantage of the matrix representation is clear. Each matrix product can be computed exactly, whereas the corresponding product of exponential operators would have required an infinite BCH expansion. For symmetric factorization, M is time-reversible. We now prove by induction that M is indeed of the form (30). Let G be of the form (30), then which is again of the form (30). Similarly for V i GV i , Since the initial G is either T i or V i , which is of the form (30), the proof is complete. If the algorithm is not time-reversible, then the diagonal elements will not be equal and M will have the more general form We will consider this more general matrix and recover the time-reversible case at the end by setting h = g. For the harmonic oscillator, because the matrix M is constant for all time steps, we have To evaluate M N analytically, we diagonalize M via For ντ > (g − h) 2 /4, the eigenvalues are unitary, where The diagonal matrix and the transform are therefore (41) The N th power of M can now be computed via where we have substituted N = t/ε, yielding finally The exponentiation of Λ then gives, where R is a pure rotation and Σ, a translation: For time-reversible algorithms, h = g, Σ = 0 and ξ = √ ντ . By comparing (46) to (17) one can quickly identify and the modified Hamiltonian conserved by the algorithm, The analytical form of the algorithm's solution (46) is interesting in that the discrete character of the algorithm has disappeared! The solution is a continuous function of t, not necessarily a multiple of ε.
For non-reversible algorithm such as (18), Σ = 0. This operator translates q → q + s and p → p − s periodically, thus causing an elongation of the phase-space ellipse along the minor diagonal at -45 • . This is the phase-space distortion alluded to earlier and observed in the standard map 9 . While non-reversible algorithms may be of interest mathematically, from the perspective of physics, there is nothing to be gain in solving timereversible dynamics with non-reversible algorithms.

V. MODIFIED HAMILTONIANS AND THE PHASE ERROR
To illustrate how this formalism works, for the SV algoritm, we have from (23), From these we can determine in closed forms We can now expand these physical quantities to any order, but for the lack of space, we will stop at the 6th: The first three terms in each case are in agreement with (15) and (16). By comparing (54) to (51) one can immediately conclude that for 0 < εω < 2, the Lie series converges. Thus the Lie series converges up to the instability point. For another example, in the case of the fourth order Forest-Ruth integrator 14 , (We use factorization coefficients as derived by Creutz and Gocksch 5 and Yoshida 15 .) From these we can determine ω A , 1/m * and k * via (48). They can then be expanded to any order, and for lack of space, we will stop at the 6th: The modified Hamiltonian now deviates from the original Hamiltonian beginning at the fourth order.
In solving for periodic motion, the most serious error is not the energy, which is well conserved by symplectic integrators after each period. The most serious error is the phase (or angle) error, which accumulates after every period. The phase error thus grow linearly with integration time, even for symplectic integrators. The phase error is simply related to the fractional deviation of the the modified angular frequency from the original frequency: For an nth order algorithm, its phase error is essentially given by the coefficient c n , For the SV algorithm, we have from (54), c 2 = 1/24. Below we list the c 4 coefficient for some recent fourth order algorithms.
Since the dynamics of the harmonic oscillator is very simple, its phase error is only a necessary but not a sufficient criterion for gauging an algorithm efficiency. While it is well known that FR has relatively large error among fourth order algorithms, it is interesting that this simple analytical calculation can indeed distinguish the latter three algorithms as significantly better than FR. (In the case of the forward algorithm, we have deliberately chosen algorithm C to do the comparison. Within the parametrized family of fourth order forward algorithms 19,20 , one can actually make c 4 vanish. Thus a general fourth order forward algorithm can solve the harmonic oscillator to six order, but such a comparison would not have been fair without allowing other algorithms to be fine-tuned.)

VI. CONCLUDING SUMMARY
In this work, we have shown that the evolution of any time-reversible symplectic algorithm can be solved in closed form. We proved that any factorized, time reversible symplectic algorithm must have a phase-space matrix of the form (30), with equal diagonal elements. Based on the exact solution, we have further shown that 1) the Lie series expansion for the modified Hamiltonian converges in the case of the harmonic oscillator. The global structure of the modified Hamiltonian is known.
2) The phase-space structure of the harmonic oscillator is closely preserved by time-reversible algorithms but is distorted by non-reversible integrators.
3) The analytical form of the phase error can be used to broadly assess the effectiveness of any reversible algorithm without doing any numerical calculation. We verified this claim by computing the phase error coefficients for four fourthorder algorithms.