UPOsHam: A Python package for computing unstable periodic orbits in two-degree-of-freedom Hamiltonian systems

Hamiltonian models are used in a diverse array of systems in natural and engineering sciences, for example, celestial mechanics, ship dynamics, chemical reactions, and structural mechanics. In two-degree-of-freedom Hamiltonian systems, the fundamental phase space structures that partition dynamically distinct trajectories and mediate transitions between multi-stable regions are stable and unstable invariant manifolds of an unstable periodic orbit (UPO). In two-degree-of-freedom systems, the phase space is four dimensional and the dynamics is constrained to the three dimensional energy surface which is partitioned by the two dimensional stable and unstable manifolds of the UPO around an index-1 saddle equilibrium point more details). the trajectories, computation and stability analysis form the starting point for dynamical systems analysis. meant to serve this by providing a module of numerical methods, along with Hamiltonian the periodic orbits at any total long guaranteed. Even though there are existing numerical methods for computing UPOs, we that they either lack in reproducibility, have a steep using the software, or have written using closed source software, and at times combination of these Child, Pechukas Farantos 1998). Our aim is to provide an open source package that implements some of the standard methods and shows the results in the context of example Hamiltonian systems. This is meant as a starting point to integrate other numerical methods in an open source package such that UPOs computed in dynamical systems papers can be reproduced with minimal tweaking while providing an exploratory environment to further develop the underlying methods.


Summary
The Python package, UPOsHam, is a collection of three methods in the form of submodules under uposham for computing UPOs around index-1 saddles in the bottleneck of Hamiltonian systems. When the form is kinetic (purely momenta-dependent terms) plus potential energy, the UPOs project as lines on the configuration space (x, y) and connect opposite points of an equipotential line V (x, y) = E. The three methods described below have been implemented as example Hamiltonian systems (also available as submodules under uposham) and are described in §:Examples. The demonstration scripts available in the package show how to import each of the methods and implement system-specific functions for computing the UPO.
The computed UPOs using the three methods for the coupled quartic Hamiltonian are compared in Figure 1.

Features: Available Methods
In this package, the user has the option to choose between the three methods described below. These are implemented in separate scripts with functions that can be modified to define the total energy (Hamiltonian), potential energy, vector field, Jacobian, and variational equations (Parker and Chua 1989).

Turning point
This method is based on finding the UPO by detecting trajectories initialized on the equipotential contour (V (x, y) = E where V (x, y) is the potetial energy function and E is the total energy) that turn in the opposite directions (Pollak, Child, and Pechukas 1980). This method relies on the fact that for Hamiltonians of the form kinetic plus potential energy, the UPO is the limiting trajectory that bounces back and forth between the equipotential contour corresponding to the given total energy. To converge on this limiting trajectory, the turning point method iteratively decreases the gap between the bounding trajectories that turn in the opposite directions. Detection of the turning point is done using a dot product condition which leads to stalling of the method beyond a certain tolerance (typically 10 −6 in the examples here.)

Turning point based on configuration difference
Based on the turning point approach, we have implemented a new method which shows stable convergence and does not rely on the dot product formula. Suppose we have found two initial conditions on a given equipotential contour and they turn in the opposite directions. If the difference in x-coordinates is small (≈ 10 −2 ), the generated trajectories will approach the UPO from either sides. If the difference in x-coordinates is large, we can integrate the Hamilton's equations for a guess time interval and find the turning point (event using ODE event detection) at which the trajectories bounce back from the far side of the equipotential contour in opposite directions. We choose these two points as our initial guess and the difference of x-coordinates becomes small. Without loss of generality, this method can be modified to either pick the difference of y-coordinates or a combination of x and y coordinates. This choice will depend on the orientation of the potential energy surface's bottleneck in the configuration space.

Differential correction
This method is based on small (≈ 10 −5 ) corrections to the initial conditions of an UPO and continuing to desired total energy. The procedure is started from the linear solutions of the Hamilton's equations and which generates a small amplitude (≈ 10 −5 ) UPO. This is fed into the procedure that calculates corrections to the initial condition based on errors in the terminal condition of the UPO. This leads to convergence within 3 steps in the sense of the trajectory returning to the initial condition. Once a small amplitude UPO is obtained, numerical continuation increases the amplitude and, correspondingly, the total energy, while a combination of bracketing and bisection method computes the UPO at the desired energy for a specified tolerance (Naik and Wiggins 2019;Koon et al. 2011).

Example systems
Consider the following two-degree-of-freedom Hamiltonian model where x, y are configuration space coordinates and p x , p y are corresponding momenta, V (x, y) is the potential energy, and T (x, y) is the kinetic energy.

Quartic Hamiltonian
This Hamiltonian can be considered as a low dimensional model of a reaction in a bath where the coupling is controlled using a parameter. The potential energy is a double-well surface and the bath is modeled using a harmonic oscillator.
where α, β, ω, are free parameters. When = 0, the system is referred to as the coupled quartic Hamiltonian, and uncoupled quartic Hamiltonian otherwise.

De Leon-Berne Hamiltonian
This Hamiltonian has been studied as a model of isomerization of a single molecule that undergoes conformational change (De Leon and Berne 1981;De Leon and Marston 1989) and exhibits regular and chaotic dynamics relevant for chemical reactions.
where the potential energy function V DB (x, y) is The parameters in the model are m A , m B which represent mass of the isomers, while s , D x denote the energy of the saddle, dissociation energy of the Morse oscillator, respectively, and will be kept fixed in this study, λ, ζ denote the range of the Morse oscillator and coupling parameter between the x and y configuration space coordinates, respectively.

Visualization of UPOs
In Fig. 1, we compare the results for the three methods for the coupled quartic Hamiltonian to show that they reproduce each other upto visual inspection. In Fig. 2, we compare the results for the turning point based on configuration difference method for the three example Hamiltonians and find they are consistent for different total energies.

Relation to ongoing research projects
We are developing geometric methods of phase space transport in the context of chemical reaction dynamics that rely on identifying and computing the UPOs. Manuscripts related to the Quartic Hamiltonian and De Leon-Berne Hamiltonian are under preparation.