Lyapunov spectra of Coulombic and gravitational periodic systems

We compute Lyapunov spectra for Coulombic and gravitational versions of the one-dimensional systems of parallel sheets with periodic boundary conditions. Exact time evolution of tangent-space vectors are derived and are utilized toward computing Lypaunov characteristic exponents using an event-driven algorithm. The results indicate that the energy dependence of the largest Lyapunov exponent emulates that of Kolmogorov-entropy density for each system at different degrees of freedom. Our approach forms an effective and approximation-free tool toward studying the dynamical properties exhibited by the Coulombic and gravitational systems and finds applications in investigating indications of thermodynamic transitions in large versions of the spatially periodic systems.


I. INTRODUCTION
One dimensional systems are of great interest to physicists in terms of their intrinsic properties and as a starting point in the analysis of their more-complicated higher-dimensional counterparts (see [1][2][3][4][5] and references therein). In the analysis of large systems considered in plasma and gravitational physics, periodic boundary conditions are preferred [6][7][8][9] and have been utilized in the study of one-dimensional plasma and gravitational systems [5,[10][11][12]. Such studies often rely on numerical simulations to validate the predictions made by the theory, and in cases where theoretical relations have not been mathematically formulated, numerical simulations serve as a powerful approach in characterizing the dynamical behaviors and thermodynamic properties [5].
Numerical studies usually employ the moleculardynamics (MD) approach in the study of the dynamical systems that undergo phase-space mixing to exhibit ergodic-like behavior. Phase-space mixing is a necessary condition for equilibrium statistical mechanics to apply and is often characterized by the existence of positive Lyapunov characteristic exponents (LCEs) [13]. LCEs represent the average rates of exponential divergence of nearby trajectories from a reference trajectory in difference directions of the phase space and quantify the degree of chaos in a dynamical system [14][15][16][17][18]. In addition, LCEs have also been reported to serve as indicators of phase transitions [19][20][21][22][23].
Numerical calculation of LCEs may be realized by studying the geometry of the phase-space trajectories [20,21] for smooth systems. In general, however, if the time evolution of each particles' position and velocity can be followed for a system, the largest LCE may be calculated by finding the rate of divergence between a refer- * b.miller@tcu.edu ence trajectory and a nearby test trajectory obtained by pertubing the former [17]. This numerical approach was extended to the case of systems with periodic boundary conditions by Kumar and Miller and was applied to find the largest LCE for a spatially-periodic one-dimensional Coulombic system [5].
While the largest LCE is a good indicator of the degree of maximum chaotic instability in a system, the mixing speed, that is, the rate at which a given phase-space volume element of the phase space diffuses across the allowed regions of the phase space is indicated by the Kolmogorov-Sinai (KS) entropy. For ergodic-like Hamiltonian systems, KS entropy is obtained as the sum of the positive LCEs [14]. In simulation, full spectrum of LCEs may be obtained by finding the time-averaged exponential rates of growth of perturbation vectors applied to the phase-space flows in the tangent space [16,[24][25][26]. An exact numerical method of calculating the full Lyapunov spectrum was proposed for the case of one-dimensional gravitation gas [16].
Even though a full spectrum is highly desirable, the calculation becomes computationally challenging for systems with large degrees of freedom and it is usually impractical to aim for a full spectrum through N -body simulations [16]. In this paper, we extend the numerical approach presented in Ref. [16] to compute the complete Lyapunov spectra of spatially-periodic one-dimensional Coulombic and gravitational systems and show that the energy dependence of the largest LCE emulates that of the sum of all the positive LCEs for the two versions of the system. The paper is organized as follows: In section II, we describe the Coulombic and gravitational versions of the model and discuss the potential interactions as well as their implications on the phase-space characteristics. In Sec. III, we recall some theoretical results for LCEs and the general numerical approach for their numerical computation. Section IV presents derivations for the time evolution of the tangent-vectors and the results of the N -body simulations. Finally, in Sec. III, we discuss the results and provide concluding remarks.

II. MODEL
We consider two versions of a spatially periodic lineal system on an x-axis with the primitive cell extending in [−L, L) and which contains N infinite sheets, each with a surface mass density m. The positions of the sheets are given by x 1 , . . . , x N with respect to the center of the cell and their corresponding velocities by v 1 , . . . , v N . In one version, the sheets are uncharged and are only interacting gravitationally. The other version is essentially a onedimensional Coulombic system in which the sheets are charged with a surface charge density q and are immersed in a uniformly distributed negative background such that the net charge is zero. For the case of the charged version of the system, we neglect the gravitational effects and take into account only the Coulomb interactions. If we denote momenta as p i ≡ mv i , the Hamiltonian of the system may be expressed as where κ = 2πkq 2 for the case of the Coulombic system [5] (with each sheet, henceforth referred to as a particle or a body, having a surface charge density q in addition to the surface mass density m) and κ = −2πGm 2 for the gravitational system [12].

A. Theoretical Overview
Here we provide a brief overview of dynamical system theory that will be helpful in developing the formulations in subsequent sections. While more general and comprehensive discussions are provided in Refs. [16,27,28], we restrict this overview to a smooth Hamiltonian system. We will see later how the concept may be extended to flows that take place on non-differentiable manifolds.
Let the phase-space flow, φ t (z) be a one-parameter group of measure-preserving diffeomorphisms M → M , where M is an n-dimensional compact differentiable manifold and z ∈ M . For a Hamiltonian system with a phasespace dimensionality of 2N , n = 2N − 1. If T z M is the tangent space to M at z, then we can define Dφ t z (w) as a linearized flow in the tangent space (T z M → T φ t z M ), where w is a vector in the tangent space. For a nonzero w, there are n independent eigenvectors e 1 , . . . , e n with χ 1 , . . . , χ n as the corresponding eigenvalues such that |χ 1 | ≥ |χ 2 | ≥ · · · ≥ |χ n |. For a periodic orbit with period t 0 , if we define λ i ≡ t −1 0 ln |χ i |, then where represents the Euclidean norm on T z M and k is a positive integer [16]. For a tangent-space vector w with a non-zero component along e 1 , the divergence for large t will be dominated by e λ1t , and therefore, λ 1 is usually called the largest Lyapunov characteristic exponent (LCE) of the orbit represented by the flow φ t (z), and is a measure of the overall stability of the orbit; if λ 1 ≥ 0, then the nearby trajectories diverge exponentially. Note that even though we have used a periodic orbit to define the largest LCE, it may be shown that the limit in the left hand side of Eq. (3) exists and is finite for any given dynamical system and the result applies rather generally under very weak smoothness conditions [29].
LCE defined in Eq. (3) may be thought of as the mean exponential growth rate of a one-dimensional "volume" (length of a vector w) in the tangent space. Therefore, λ 1 is often referred to as LCE of order 1. Similarly, λ p , that is, LCE of order p (where 1 ≤ p ≤ n, p ∈ Z + ), may be related to the mean exponential rate of growth of a p-dimensional hyperparallelepiped formed by the evolution of p linearly independent tangent-space vectors w 1 , . . . , w p . We first find the rate of volume divergence as where Vol P represents the volume spanned by a set of p tangent-space vectors. Finally, following Ref. [29], λ p is found as

B. Numerical Approach
We start with a randomly chosen set of n orthonormal tangent vectors {ŵ 0 1 , . . . ,ŵ 0 n }. Clearly, for each p ≤ n, Vol P ŵ 0 1 , . . . ,ŵ 0 p = 1. After a fixed time interval τ , the evolved tangent vectors-which we denote by {w 1 1 , . . . , w 1 n }, where w 1 i = Dφ t=τ z (ŵ 0 i )-are, in general, no longer mutually orthogonal. This is because the component of eachŵ 0 i along the direction of maximum divergence e 1 (that is,ŵ 0 i · e 1 ) will witness a disproportionately larger growth in its value as compared to the remaining components. In order to avoid numerical errors arising from one component getting increasingly large in comparison to the others, a new orthonormal set of tangent vectors {ŵ 1 1 , . . . ,ŵ 1 n } is defined after time τ through Gram-Schmidt reduction on the set of evolved {w 0 1 , . . . , w 0 n }. This new set of ornothormal tangent vectors are then used for the following iteration, and the process is recursively repeated until λ P has converged [16,28].
Numerical calculation of λ P involves finding the corresponding p-volume for each iteration. If at the end of the j-th iteration, ) represent the evolved versions of the orthornormal tangent vectorsŵ j−1 i , the p-volume may be found as the norm of the the exterior product involving the corresponding p vectors, that is, Finally, the average exponential growth rate of the pvolume is found as where l is the total number of iterations. A complete set of LCEs {λ 1 , . . . , λ n }, also known as Lyapunov spectrum, may then be obtained for the trajectory by utilizing Eq. (5) for all permissible values of p. Finally, an upper limit on the KS entropy h KS for Hamiltonian systems may be obtained as the sum of positive LCE's [14,16], where p + max is the largest value of p for which λ p is positive and where the equality h KS = λ S holds for ergodiclike systems. The sum of the positive LCEs λ S is often termed as the density of KS entropy.
In the following section, we discuss how we employ the theory and the numerical approach presented thus far to obtain Lyapunov spectra for the Coulombic and gravitational systems discussed in Sec. II. While most of the presented theory applies to the two systems in its original form, non-smoothness arising from the absolute valued linear terms in the potential demand additional consideration. It turns out that, as we shall see, following the time evolution of the tangent-space vectors involves treating the motion as a flow in between two consecutive events of interparticle crossings and as a mapping at each event of such crossings. Consequently, it becomes indispensable to have the ability to find the exact time corresponding to each crossing.

A. Equations of motion
Positions and velocities of the particles are obtained using event-driven algorithms based on the approaches proposed in Ref. [5] for the Coulombic system and in Ref. [12] for the gravitational system. The algorithms employ analytic expressions for the time dependencies of the relative separations Z j (t) and relative velocities W j (t) between two consecutive particles in the primitive , with x j and v j representing, respectively, the position and velocity of the j-th particle whereas x j+1 and v j+1 representing those of the (j + 1)-th particle. Combining the results of Refs. [5,12], we find that Crossing times may be found by solving Z j (t) = 0 for t. The corresponding positions x j and velocities v j are obtained using a matrix-inversion subroutine as described in Ref. [5].

B. Time evolution of tangent-space vectors
In order to follow the time evolution of the tangent vectors, we adopt an approach based on the "exact" numerical method proposed in Ref. [16]. The method invokes that, for a one-dimensional Hamiltonian system with N particles, one does not have to restrict to the (2N − 1)-dimensional manifold Γ E .
One may alternatively choose to represent the flow φ t in the entire 2N -dimensional phase space (say, Ω) whereby the tangent space T z Γ E becomes a subspace of T z Ω.
Let z(x, v) be a point in the phase space Ω, where x = (x 1 , . . . , x N ) and v = (v 1 , . . . , v N ). The equations of motion representing the system are given bẏ andv with Similarly, if we have a vector w(ξ, η) in the tangent space T z Ω, then the variational equations [16] governing the evolution of w are given by where I N is the N × N identity matrix and A(x) is an N × N matrix whose elements are given by For the potential expressed in Eq. (12), one finds that (15) Using Eqs. (13) and (15), we may deduce thaṫ anḋ Equations (IV B) and (17) imply that which implies that Solution to Eq. (21) yields where Hence, in between the events of crossings, the evolution of tangent vectors takes the form With the two values of κ, one for the Coulombic system and the other for the gravitational system, one may solve Eq. (24) to find the exact dependencies of ξ j and η j on time.

Coulombic system
Utilizing the initial conditions for the Coulombic system with κ = 2πkq 2 , we obtain the solutions to Eq. 24 for ξ j and η j in between crossings as and, where ω ≡ κN mL = 2πkq 2 N mL . If the r-th and s-th particles undergo a crossing at time t = t c , and t − and t + respectively denote the instants just before and after t = t c , then

Gravitational system
For the gravitational system, κ = −2πGm 2 . Similar to the Coulombic case, we utilize the initial conditions to find solutions to Eq. 24 for ξ i and η j as functions of time for the gravitational case in between the events of crossings: where, and Λ ≡ kN mL = 2πGmN L . For the r-th and s-th particle involved in a crossing, we get where we have used the same definitions of t c , t − , and t + as we did in the expressions' Coulombic counterparts.

C. p-volume and Lyapunov spectrum
We perform numerical computations using algorithms that are driven by tracking the events of interparticle crossings. Since the time derivatives of the velocities of the particles involved in a crossing undergo abrupt changes, tracking of crossing becomes indispensable. In other words, even if we were to sample the positions and velocities at fixed intervals of time, we would still need to track every single event of crossing [16]. Therefore, instead of choosing fixed time intervals for Gram-Schmidt orthonormalization of the tangent-space vectors, we choose to do it after each crossing. It should be noted that the duration of each iteration, fixed or variable, is irrelevant as long it remains short enough for the various components of the tangent vectors to remain comparable with a given precision offered by the computing platform utilized.
After an iteration, say, the (j − 1)-th iteration, ending at time t = t j−1 c , each of the 2N orthonormal tangent vectorsŵ j−1 1 , . . . ,ŵ j−1 N is allowed to evolve for the duration δt j c of the j-th iteration. δt j c may be thought of as the time elapsed between the instants right after the (j − 1)th crossing and right before the j-th crossing. Then the evolved vectors may be given by w j i = Dφ For the j-th iteration, p-volume is found as follows: We define a p × p symmetric matrix G j p whose elements G j p µν are inner products between w j µ and w j ν , where 1 ≤ µ, ν ≤ p (with µ, ν ∈ Z + ). That is, The matrix G j p , also referred to as Gram matrix, encapsulates the necessary geometric information about the subspace spanned by the set of vectors {w j 1 , . . . , w j p } such as the lengths of the vectors and the angles between them. The absolute value of the determinant of G j p , known as the Gramian, is essentially the square of the norm of the exterior product [30]. Using Eq. (6), the p-volume is found as With the ability to find each p-volume for a given iteration, the final value of the corresponding λ P are found as where t l c , as defined in Eq. (35), is the total time elapsed for l crossings to occur. Finally, LCEs λ p and Kolmogorov-entropy density λ S are obtained using Eqs. (5) and (8) respectively.

D. Results
In our simulation, the initial conditions are chosen as follows: For a given number of particles N and perparticle energy H, if H is lower than the maximum allowed value of the potential energy V max , then the particle positions are chosen randomly such that the potential energy is slightly smaller than the target value of H. For H greater than V max , the particle positions are randomly selected such that the potential energy is close to V max . Velocities are chosen randomly from a Gaussian distribution and are scaled such that the sum of the potential and kinetic energies exactly equals the target value of H [5].
Simulations are performed by rescaling the system parameters in a system of dimensionless units such that the number density N/2L = 1 and the characteristic frequencies ω, Λ equal unity [5,12]. Energies H c and H g (respectively for Coulombic and gravitational systems) are measured with respect to the minimum values of the corresponding potential energies allowed for each N .
Before running the Lyapunov algorithm, we allowed the system to evolve for a relaxation period of 500 time units (in terms of 1/ω or 1/Λ). In the calculation of LCEs, the total number of iterations l for each λ P in Eq. 38 are decided by an adaptive algorithm which runs a minimum preliminary number of iterations assigned beforehand and then continues running until the value of λ P has converged to within a pre-specified tolerance for the standard deviation. The tolerance value that we specified was 0.1 percent of the mean from the newest 500, 000 iterations, with a minimum of 1 million iterations. In each run, we found that the values converged to within our specified tolerance after the preliminary run of 1 million iterations.
We computed Lyapunov spectra for the Coulombic and gravitational systems with varying number of particles N with 5 ≤ N ≤ 20. Figure 1 shows examples of the dependence of the various LCEs on the per-particle energy for the Coulombic and gravitational systems with N = 11. The figure also shows the sum of all LCEs for each case. As expected, the sums were found to be close to zero.
The energy dependencies of the largest LCE and the Kolmogorov-entropy density with N = 5, 8, 11, 15, and 20 have been shown for the Coulombic and gravitational  systems in Figs. 2 and 3 respectively. It can be seen in each case that the behavior of λ S versus H/N resembles, in general, a scaled version of λ 1 versus H/N . To elucidate this, we have presented comparative plots of the normalized versions,λ 1 andλ S , of λ 1 and λ S against per-particle energy for N = 11 in Fig. 4, where we have divided λ 1 and λ S by their respective maximum values to get the normalized values.

V. DISCUSSION AND CONCLUSIONS
Our study provides interesting insight into the chaotic dynamics of the two versions of the periodic system under the conditions of varying energy and degrees of freedom. The results of the Coulombic system are consistent those provided in Ref. [5]. λ 1 for the Coulombic system stays zero as long as the energy is low enough for the particles to not undergo crossings. As the energy is increased  from a low value, λ 1 sees an initial increase, reaches a maximum and then decreases as the energy is progressively raised. However, with an increase in the number of particles, the right edge of "hill" near the maximum opens up to form a "plateau" and flattens out asymptotically as shown in Fig. 2(a). Interestingly, our results also suggest that all the other LCEs are, more or less, scaled (and inverted, for the case of negative LCEs) versions of the largest LCE, as we can see in Fig. 1(a) for N = 11.
Unlike the Coulombic version, the gravitational system shows a maximum degree of chaos at low energies. Our results suggest that the largest LCE starts off with a high value at low energies and decreases to a minimum as the energy is increased for a given N . With further rise in energy, λ 1 increases, reaches and maximum and then decreases asymptotically. As the number of particles is increased, the local trough near the local minimum and the hill near the local maximum start opening up to the right with the asymptotic edge rising up to a form a plateau as shown in Fig. 3(b). Moreover, it can also be seen in Fig. 3(b) that the local minimum and maximum themselves shift to the right as well with an increase in the number of particles. It should be noted that in the thermodynamic limit, that is, for versions of the system with sufficiently large N , the LCEs may exhibit discontinuities in their values or their slopes near the troughs and crests when plotted against temperature. Such an observation may indicate toward existence of phase transitions [22,23].
While the energy dependencies of λ 1 are drastically different for the two versions of the system, the other LCEs for the gravitational system are also, in general, scaled versions of λ 1 , similar to the Coulombic system, as exemplified in Fig. 1(b) for N = 11. Moreover, in both versions of the spatially-periodic system, the energy dependence of LCEs tends to approach a common limiting behavior as N is increased. This is in contrast to the freeboundary gravitating case in which the values of LCEs were shown to increase linearly with an increase in the number of particles [16].
From a dynamical perspective, the results are consistent with the theoretical predictions for Hamiltonian systems [16] as outlined below: (a). The sum of LCEs was found to converge to zero for all energies and N .
To summarize, we have provided an exact method to compute full spectra of LCEs for the spatially periodic versions of the one-dimensional Coulombic and gravitational systems. Analytic expressions for time-evolutions of the tangent vectors were derived and used in numerically computing LCEs using an efficient event-driven algorithm. While the resulting values of the largest LCE for the Coulombic system agree with reported previously [5] , our exact approach offers striking advantages over the method used in Ref. [5] in that it allows one to calculate a full spectrum of LCEs rather than just the largest LCE. Second, the results of the exact method do not depend on the size of the perturbation. In finding the largest LCE using finite perturbations to a reference orbit as discussed in Ref. [5], one has to first make sure that the value chosen for initial perturbation is small enough. This becomes challenging for systems in which particles tend stay clumped together. An example of such behavior is seen in the gravitational system at low energies. Finally, the exact approach circumvents the need for defining a test trajectory altogether, which as discussed in Ref. [5] poses difficulties in expressing phase-space separations for systems with periodic boundary conditions. Nevertheless, the method discussed in Ref. [5] still remains powerful, and perhaps the only resort, in dealing with spatially-periodic systems for which analytic evolution of tangent-space vectors may not be obtained.
The results of our study also indicate that the energy dependence of the largest LCE captures the general behavior of the dependence of the Kolmogorov-entropy density on energy for both Coulombic and gravitational systems. This result is particularly significant because of the numerical difficulties encountered while calculating the full Lyapunov spectra of large systems. For the two versions of the spatially-periodic system, our study suggests that one may gain insights into a full spectrum of LCEs by simply looking at the largest LCE, thereby allowing one to evade the computational complications faced when calculating a full spectrum.
It should be noted that, for a given number of particles, the energy dependence of the LCEs roughly followed the same behavior for any randomly selected initial conditions with only slight deviations. As the number of particles was increased, the deviations became smaller leading the behaviors to converge to a single universal one, indicating toward the approach of ergodic-like nature. Moreover, the convergence times of the LCNs for different randomly-chosen initial conditions also showed uniformity for larger number of particles, thereby pointing toward a consistent relaxation to equilibrium with increasing degrees of freedom. However, the exact dependence of relaxation time on the number of particles requires further investigation and we plan to pursue it in our future work.
Finally, it is worth emphasizing that if a phase transition occurs in either of two spatially-periodic systems, the temperature dependence of the largest LCE is expected to show a transitioning behavior in the thermodynamic limit (large N -limit) [22,23,31]. In our future work, we plan to utilize the various previosuly reported numerical techniques [5] as well as the approach presented in this paper to examine the chaotic and thermodynamic properties of the periodic gravitational and Coulombic systems for indications of phase transitions.