A Simple Model for Long-Range Interacting Pendula

We show that the Hamiltonian mean field (HMF) model describes the equilibrium behavior of a system of long pendula with flat bobs that are coupled through long-range interactions (charged or self gravitating). We solve for the canonical partition function in the coordinate frame of the pendula angles. The Hamiltonian in the angles coordinate frame looks similar to the form of the HMF model but with the inclusion of an index dependent phase in the interaction term. We also show interesting non-equilibrium behavior of the pendula angles, namely that a quasistationary clustered state can exist when pendula angles are initially ordered by their index.


I. INTRODUCTION
Systems with long-range interactions are a source of unique problems in the field of statistical mechanics and thermodynamics. This is due to several properties of long-range systems which fall outside of the conditions normally needing to be satisfied when applying the methodologies of thermodynamics. Simply from the words "long-range" the first infringement can be deduced, that long-range systems are not additive. If two systems with short-range interactions are brought together to form a larger system then the energy difference between the conglomerate system and the sum of its constituents is the new potential energy from the boundary between them. In the thermodynamic limit, the potential energy of the boundary is small compared to the bulk and can be neglected, making short-range systems additive. In the case of long-range interactions, one particle will feel a significant potential created by every other particle, so the additional potential energy of two systems added together does not scale as the boundary but in a more complicated way that depends on the specific nature of the interactions [1]. Directly related to the lack of additivity is the fact that systems with long-range interactions are not extensive because their energy diverges in the thermodynamic limit [2]. Although * oweenm@gmail.com these characteristics compel cautious use of the usual tools of statistical mechanics, they are also the source of many interesting dynamical and statistical features. Depending on the system of interest, such features include canonical and microcanonical ensemble inequivalence and related negative specific heat [3], quasistationary states (different than metastable states which lie on local extrema of equilibrium potentials) whose lifetimes increase with the number of particles [4], an interesting dependence of the largest Lyapunov exponent on particle number in a longrange Fermi-Pasta-Ulam model [5], and spontaneous creation of macroscopic structures in nonequilibrium states [6]. In some cases, long-range interactions can greatly simplify problems. For instance, mean field models depend on one of two premises: (i) interactions are short-range but the system is embedded in a space of infinite dimension so that all bodies in the system are nearest neighbors, or (ii) interactions are infinitely long.
For some time, the primary motivation for the study of long-range interactions was to understand galaxies, galaxy clusters and the general thermodynamic properties of selfgravitating systems. Aside from mean field models, interest has further built since the observation of modified scattering lengths in Bose-Einstein Condensates (BEC) through the use of Feshbach resonances [7]. Using this technique, a BEC can be made to be almost non-interacting by tuning the scattering length to zero. One could even tune the scattering length to a negative value, making the BEC collapse. More recently, O'dell et al. [8] has shown that it may be possible to produce an attractive 1/r potential between atoms in a BEC by applying an "extremely off resonant" electromagnetic field. This has opened the possibility of creating table-top methods which physically model aspects of cosmological behavior on a laboratory scale, as well as the possible development of entirely new dynamics in BEC.
The challenges in understanding long-range systems drive the development of solvable models that could help better explain some of the aforementioned phenomena. Campa et al. [9] have recently published a collection of important solvable models. One particularly significant model, which is important to this work, is the Hamiltonian Mean Field (HMF) XY spin model [6], often written in the form where θ i is the angular position of the i th particle (spin), as shown in Fig. 1, and p i is its conjugate angular momentum. The HMF model is generally used to describe two different classes of systems: 1) a mean field XY classical spin model, and 2) a one dimensional periodic system of itinerant particles with long-range interactions. Though the connection between the HMF model and the second class of systems mentioned could be thought of as contrived given the simplifications under which the model is realized, it has been shown that the model produces useful insights into how non-neutral plasmas and self gravitating systems behave [6].
In this paper, we study the dynamics of an array of N pendula with long-range interacting bobs. By considering long pendula with flat bobs undergoing small oscillations and having parallel planes of rotation, we produce a model related to the HMF model through a coordinate transformation. The transformation introduces a dependence on the indices of the particle labels. A cartoon of the physical picture is shown in Fig. 1.
The index dependence in the Hamiltonian, that Figure 1. N pendulum system with parallel planes of rotation. The i th pendulum angle at some time t is θ i (t).
will be described in detail in the next section, is a consequence of the pendula pivots being slightly offset from one another and appears as a phase in the cosine term of the HMF model. It inspires the investigation of non-equilibrium "repulsive" behavior in the angle coordinate frame where we find an interesting quasistationary state when the angles of the pendula are initially ordered according to their indices. We find the clustered positions in the usual HMF coordinate frame (biclusters), but in the angle coordinate frame clustering is only found for the initially ordered angles and, unlike the biclusters, these are clearly quasistationary states. A quasistationary state is defined as a dynamical state that can persist for a length of time which goes to infinity as the thermodynamic limit is approached [9]. In addition to discussing the clustered angle states exhibited by the system, we also solve for the canonical partition function in the pendulum angle coordinate frame, finding that in equilibrium with a heat bath, the probability distributions of the angles can be described by the original HMF model. This finding is similar to the work done by [10] on a model sometimes called the HMF αmodel. In the HMF α model, a 1/r α ij dependence between the classical spins is introduced [9,[11][12][13][14], where r ij is the distance between the i th and j th spins on a lattice. Though the physical motivations behind studying these various models can be very different, it is interesting that their equilibrium behavior is the same or nearly the same. We believe that the work in this paper further suggests that the HMF model universally describes an entire class of long-range interacting systems in equilibrium.

A. Coordinates
In Fig. 1, we show an array of pendula rotating in the same plane with bobs that interact through a long-range potential. If we consider the case where all the pendula only undergo small oscillations, we may write the horizontal location of the i th particle, x i , as x i = id + θ i , where d is the distance between the pivots of neighboring pendula and is the length of each pendulum. The small θ regime makes the problem one dimensional in x. We choose periodic boundary conditions and rescale the system by 2π/N d so that making the total system length a dimensionless 2π where N is the number of particles in one period. We will refer to a periodic space with length 2π as a unit circle. The position of the i th particle (bob) is now For reasonable choices of and d ( /d << N ), the second term on the RHS is suitably small such that the Hamiltonian can be written with terms that are quadratic in θ. However, we are primarily interested in a regime where /d → ∞ as the thermodynamic limit is approached. Physically this corresponds to the small oscillations of very long pendula with suspension points that are close together compared to their lengths. In order to simplify the calculations that follow, we define φ i to be the last term on the RHS of Eq. (3), namely φ i ≡ 2π θ i /N d.
Given the choice of large /d, φ i can take any value in the range [0, 2π). This is only true because /d is large, not because the θ i s are. In terms of φ i , the positions can be rewritten as

B. Density Approximation
We have not yet explicitly stated the physical mechanism through which the bobs interact. Connecting the interactions with specific physical motivations should be discussed with some discretion because the development of the model leaves these motivations up to some freedom of interpretation. Imagine that the bobs all carry some charge. We will not distinguish between particles in any other way than their indices, so in the case where all particles carry the same charge, repulsive behavior is expected. On the other hand one could make the bobs attract one another, which could be thought of as the self-gravitating case. To be solvable, the model requires some simplifications. For the sake of brevity we will speak of the particle charge or mass density as the "density".
The approximation that we invoke is similar to that used when justifying the HMF model (Eq. (1)) to describe free particles in a onedimensional ring [6,15]. The distribution of the bobs is such the mass density, ρ(x), is given by The constant 1/2π subtracted from the delta function is necessary to produce a meaningful expression for the potential Φ and corresponds to the inclusion of a neutralizing (of opposite sign) homogeneous background density. Restricting the problem further to that of solving Poisson's equation for a one-dimensional potential physically amounts to choosing large and flat bob geometries oriented with their smallest axis parallel to the x axis. Writing the delta function as a cosine Fourier series, Poisson's equation becomes The parameter γ contains the particle (bob) charge or mass and becomes the interaction strength in the Hamiltonian. We can see that the zeroth-order term in the Fourier series canceled the constant neutralizing background that was superficially added. The most important simplification in this paper is truncating the sum of the Fourier coefficients used to represent the delta function after the n = 1 coefficient. Antoni et al. defend the truncation by asserting that the "large scale collective properties" do not greatly change when higher order terms of the sum (including interactions at the smaller length scales) are included, and discuss the consequences of the approximation in some detail [6]. The simplification also warrants a brief discussion of the way that it could be physically interpreted. The truncation of the sum is equivalent to smearing out the density of each particle over the system so that it is peaked at its given location, x i , but also having a negative density peak on the opposite side of the unit circle. This could be thought of as doubling the number of particles and enforcing that each particle has a negative partner that always remains on the opposing side of the unit circle. After this doubling, the now nebulous masses are dispersed such that a pair's density is described by a cosine function with the positive peak centered at x i .

C. Solving Poisson's Equation
Integrating Poisson's equation once, we obtain: In order to determine the constant c 1 from the integration, the physical picture should be examined. A sensible requirement is that when all of the bobs are hanging at their equilibrium positions, directly below their pivot (all φ i = θ i = 0), the net force experienced by any bob is zero. This is a valid requirement if the bobs are attractive or repulsive, the only difference being that the configuration would be unstable or sta-ble, respectively. The force that the j th particle experiences when φ j and all φ i are zero is given by The sum i sin [2π(j − i)/N ] equals zero for any j, so c 1 must be zero. Integrating once more to obtain the potential yields To determine c 2 we stipulate that if all φ i = 0, then Φ(0) = 0. Inserting Eq. (3) (or Eq. (4)) for The sum over the cosine is zero, therefore c 2 = 0 and we can now write the potential energy of the j th particle as

D. The Hamiltonian
The Hamiltonian can be written as where H 0 is the kinetic energy piece and is the interaction piece, so The mass of the bobs has been set to unity, γ is the interaction strength, a factor of 1/2 accounts for the double counting, and the 1/π coefficient in the potential energy has been absorbed into γ. The factor of 1/N is a rescaling of the potential energy that ensures that as the thermodynamic limit is approached, the potential energy of the system does not diverge. The 1/N scaling is known as the Kac prescription [16]. The Kac serves to keep both the energy and entropy of a system proportional to the number of particles in the system, an important prerequisite for phase transitions [9].

E. Relationship to the HMF model and the Spin Interpretation
Due to the simple bijective relationship between x i and φ i one can simply solve the equations of motion for the HMF model and find the dynamics for φ i via the coordinate transform x i → φ i . Previously it was mentioned that the HMF model is used to describe free particles on a ring with long-range repulsion or attraction, as well as describing a classical XY spin model. The θ i played the role of either the position of the i th particle on the ring or the orientation of the i th spin. Therefore, it is interesting to speculate about the type of spin system the model describes in the φ i picture. Thus far, the rescaled angle φ i = 2π θ i /N d describes the distance of a pendulum bob from the point directly below its pivot, but it could also be interpreted as the orientation of spin. In the spin interpretation of Eq. (15), the potential energy of the i th and j th spin pair depend on both their relative orientation as well as the difference between their indices. In the following discussion it will sometimes be convenient to speak about φ i in the spin language.
We will prove that in the φ picture, the system in equilibrium with a heat bath is equivalent to the HMF model (the x picture) in equi-librium with a heat bath by solving the partition function in the φ i coordinate frame. In the process of simplifying the Hamiltonian to solve for the partition function, we will find expressions of the form cos φ i and sin φ i which we talk about as the horizontal and vertical components of a magnetization m i = (cos φ i , sin φ i ). It could easily be stated that in the spin analogy, the φ i are orientations of the spins, but we should make a more concrete connection between this idea and the original presentation of the model. We would like to remind the reader that even though the angles θ i of the pendula are small, the long suspensions of the bobs ( ) allow φ i to cover the entire system which, rescaled, has dimensionless length 2π. The system is also periodic, so the bobs can be thought of as moving on a unit circle where the position of the i th bob is x i = 2πi/N + φ i . In order to think of φ i as the spin orientations, we start by considering each bob as living on its own individual unit circle. An example of these unit circles is shown in Fig. 2, a visual aid to the following. Imagine stacking horizontal circles in the vertical direction and rotating each by an angle 2π/N from the one below. The projection of these circles onto the horizontal plane would be the system viewed in x, i.e. the HMF model. If we twist the stack so there is no rotation between adjacent circles and then project onto the horizontal plane it creates the picture viewed in φ, where the pivot points are all aligned. The reason for this artificial construction of stacked circles is partly to pictorially depict the transformation between x and φ and partly to show how m i (as defined) is just the orientation of the i th spin in the φ picture. Said differently, each circle in the φ picture represents a spin with an orientation in the horizontal plane determined by φ i ; an infinite-range classical mean field spin model described by Eq. (15).

III. EQUILIBRIUM
In this section, we solve for the canonical partition function, in the φ coordinate frame using the Hamiltonian in Eq. (15) and show that, in equilibrium, the HMF model describes the an- gles of long pendula with long-range interacting bobs. In order to solve the configurational piece of the partition function the Hamiltonian must be modified. Using the cosine and sine sum and difference identities twice, the potential interaction piece of the Hamiltonian H I can be written as The coefficients in the Hamiltonian cos [2π(i − j)/N ] and sin [2π(i − j)/N ] should be thought of as matrices with components C ij and S ij respectively. The Hamiltonian can now be written in the form (17) which is suggestive because it can be regarded as the matrix equation It is helpful to consider the particles positions on the unit circle with respect to their pivot (φ) as magnetizations. Defining and with m T µ = (m 0,µ , m 1,µ , ..., m N −1,µ ) where µ holds the place of x or y, the Hamiltonian becomes A closer examination of the structure of the coefficient matrices C and S indicates that they take the special form of circulant matrices, and thus can be simultaneously diagonalized by a unitary matrix U . A circulant matrix has the form a special kind of Toeplitz matrix, where each subsequent row is a cyclic permutation of the row above or below it. Any matrix A with elements a ij that can be written in terms of some function f (i − j) is a circulant matrix. Because a circulant matrix is a normal matrix it can be diagonalized by a unitary matrix. We show that The matrices C and S can be rewritten as C = U † D C U and S = U † D S U , where D C and D S are diagonal matrices with diagonal elements that are the eigenvalues of C and S, respectively, which we denote as λ C i and λ S i . From here on we label the indices i from 0 to N − 1. It is worth pointing out that due to S being antisymmetric, U must be complex. Equation (22) becomes We will move back to the index notation using the following relations: and where X is not to be confused with x. Using the Kronecker delta, we set all i = j since these are the only nonzero terms. The Hamiltonian is now given by with X = XX * and Y = Y Y * . The inclusion of the eigenvalues λ C and λ S simplifies the Hamiltonian further. We will now solve for λ C and λ S . Looking at the form of a circulant matrix shown in Eq. (21) reminds us that the elements of a circulant matrix can be defined with a single label. We write the single labeled elements of the cosine and sine matrices respectively as and where l = 0, 1, 2, ..., N − 1. The eigenvalues, λ A , of a N × N circulant matrix A can be written in terms of the single label elements a l . The j th eigenvalue of A is known to be λ A j = N −1 l=0 a l exp (2πilj/N ), where i is √ −1 (not an index) and l = 0, 1, 2, ..., N − 1. Therefore, and Writing cosine and sine in their exponential forms gives The above representations of the eigenvalues show that C and S each have only two non-zero eigenvalues corresponding to j = 1, N − 1 given by λ C 1 = λ C N −1 = N/2 and λ S 1 = (λ S N −1 ) * = iN/2. The Hamiltonian simplifies greatly to The representation of H I in Eq. (33) must be further modified before the partition function can be found. We do this by splitting the Fourier matrix U into its real and imaginary components, a ik and b ik , given by This was done to write the absolute squares in Eq. (33) in terms of the squares of a ik and b ik . By noticing that a 1k = a (N −1)k and b 1k = −b (N −1)k we write the configurational partition function as where β = 1/k B T .
The Hubbard-Stratonovich transformation is now applied twice, once to each quadratic quantity in the partition function. The integration variables introduced through this transformation are z 1 and z 2 with subscripts for first and second quadratic quantities, respectively. After after switching the order of integration, we find The integration can be performed using the identity π −π dφe ξ cos φ+η sin φ = 2πI 0 where which simplifies when a and b are included to It is convenient to make a change to polar coordinates by introducing z = z 2 1 + z 2 2 , following which the partition function can be written as Equation (41) is recognized to be an intermediate step of the solution to the canonical partition function for the HMF model. From here we jump to the main results, the details of which are included in the HMF literature [6,9,15] . The integration over z in Eq. (41) can be preformed using the saddle point approximation. The rescaled free energy per particle follows as The expression above permits a convenient path to finding the phase transition. Solving for the minimum values of z in order to satisfy the last term in Eq. (42) results in the equation which can be solved self consistently for z and represented graphically for different values of β as in Fig. 3. The reader will see that after β is increased passed the critical value (β = 2) there are two well-defined solutions. The Hubbard-Stratonovich transformation decouples spin-spin (squared terms in the Hamiltonian) contributions to the partition function at the price of needing to create a linear interaction between each spin with an auxiliary field z [17]. Again, a more detailed procedure can be found in [6,9,15] where discussion of the internal energy in the equilibrium state is followed by non-equilibrium behavior of the system prepared in microcanonical ensembles. Here we will simply touch on the most important point of the equilibrium behavior, being that for β < 2 the system is paramagnetic but for β ≥ 2 a pitchfork bifurcation occurs resulting in two stable solutions. At this point there is a discontinuity in the free energy, a second order phase transition occurs and the system can maintain finite magnetization. In this case, the order parame- Showing that the canonical partition function in the φ coordinate frame model and the HMF model are equivalent necessitates a more detailed discussion of the equilibrum behaviors in the φ frame. Campa et al. [9], in their review of the HMF model, rigorously show ensemble equivalence between the canonical and microcanoical ensemble of the HMF model. In light of this fact, a large N microcanonical simulation should be able to produce equilibrium behavior like the phase transition mentioned above. The temperature in a numerical simulation of a system with many particles would be "set" through a choice of the initial momenta distribution. In this type of simulation, it is common practice to compute the order parameter and free energy [1,6,18], begging the question: does a large microcaonical simulation of Eq. (15) approximate the expected equilibrium behavior? Also, since the index-dependent model in equilibrium with a heat bath can be described by the HMF model, would the dynamics of such a simulation qualitatively resemble those in the HMF model? The answer to both of these questions is no if one were to find the equations of motions in φ i for some large N and then compare them to an HMF model or the x coordinate frame. As stated, this discrepancy may appear to detract from our result. Indeed, it uncovers a conceptual omission in the model, but it is one whose rectification gives insight into the models ensemble equivalence property of the model, or lack thereof. The omission was in the arbitrary scaling of x which we will now rectify.
We introduce the parameter L which generalizes the scaling in Eq. (2) to making the position of the i th particle and changing the definition of φ i to φ i ≡ 2πL θ i /N d. It can be shown that the introduction of L only changes the final result of the partition function by a constant factor of L due to the enlarged limits of integration. Numerically, we find is that if L >> 1, then the simulations in φ closely reproduce the dynamics of HMF model simulations (dynamics in x). Therefore, for large L the mirocanonical simulations can approximate equilibrium and the answers to the previous questions -does a large microcaonical simulation of Eq. (15) approximate the expected equilibrium behavior, and since the index-dependent model in equilibrium with a heat bath can be described by the HMF model, would the dynamics of such a simulation qualitatively resemble those in the HMF model? -becomes yes. Alternatively, the coordinate frame inequivelence is most extreme for small L. These numerical results were found using initial conditions that are randomly distributed φ i about the domain [−Lπ, Lπ). It should be stated that for the rest of this paper we work with L = 1 becuase we are inetersed in cases where the φ coordinate frame is markedly differnet than the x coordinate frame.

IV. NON-EQUILIBRIUM RESULTS
For a system of pendula, it is interesting to study an initial configuration where all pendula are set to random small displacements from φ i = 0. Specifically we initialize the i th pendulum angle, φ i , randomly in the range [−π/N, π/N ). In x the indices are ordered in x such that x 1 < x 2 < x 3 < ... < x N and the i th bob is randomly distributed in the range [2πi/N − π/N, 2πi/N + π/N ). It is possible to make some general statements about the dynamics of this configuration in x using the equations of motion. Expressing the Hamiltonian with terms that are quadratic in φ yields With this expression, the equations of motion for the i th particle can be written as from which we obtain In the above equation, the last term and the cos [2π(i − j)/N ]φ i term sum to zero, leading tö Using the difference formula, we write the acceleration as where µ 1 = φ j cos (2πi/N ) and µ 2 = φ j sin (2πi/N ). The mass (moment of inertia) has been set to unity so the above expression is the force as a function of index,φ i = F (i).
If µ 1 and µ 2 are known, then the initial dynamics of the system are elucidated by Eq. (50), but in the case of randomly initialized φ i the µ 1 and µ 2 are also random and can be different from one another in both magnitude and sign. However, a general description of the results can be given without exactly knowing these coefficients. Equation (50) shows that the initial force on a given particle depends on its position because the indices are ordered in x. In the continuum (thermodynamic limit), the force takes the form Therefore µ 1 and µ 2 partly play the role of the amplitude of this force as a function of x, but also can shift the cos x + sin x spatial dependence, which is periodic over the system length.
In Fig. 4, we show a fit of the force as a function of x using Eq. (51) as well as the actual force calculated for an example set of initial conditions. The domain in Fig. 4 can be split into two pieces (independent of µ)-one where the particles experience a positive force, the other in which the particles experience a negative force. As time is increased, the movements of the particles evolve the coefficients µ 1 and µ 2 in such a way that the magnitude of the force decreases to zero for all particles and then switches sign complementary to the original force. This results in a standing compression wave of the particles with a wavelength 2π. The compression wave is not stable and eventually two clusters form about each node. These two clusters are often referred to as a "bicluster", or the antiferromagnetic state in the HMF model, and have been explained by Barré et al. by analysing the Vlasov equation. They find that the initial com- pression wave (referred to by a different name) creates an effective double-well potential giving rise to the bicluster [19]. The question of the bicluster stability has not yet been definitely answered, but for a detailed discussion we refer the reader to Leyvraz et al. [20]. Given the simple mapping between the φ and x coordinate frames, we should also be able to show the initial form of the force in x as well. As presented in Eq. (15), the Hamiltonian in the x coordinate frame only differs from the HMF model by a constant γ/2. In x, H I is Using the difference identity, we find the equations of motion for the i th particle to bë (53) The sums over cosine and sine of x j play the same role as µ 1 and µ 2 , and the force at a given position x i is clearly of the same form as that shown in Eq. (50).
Depending on µ 1 and µ 2 , all φ i oscillate about zero with amplitudes and phases that depend on their location x i as discussed above. As the clustering in x begins, the φ i begin to spread out over the full domain [0, 2π) and continue to do so until it is covered. The more interesting case in φ is when all φ i are initially randomly distributed in ranges that depend on their index, specifically when φ i are chosen in the ranges. [2πi/N − π/N, 2πi/N + π/N ) so that φ 1 < φ 2 < φ 3 < ... < φ N . It should be noted that in this new configuration the dynamics in x are nearly identical to the configuration previously discussed for ordered x i . The dynamics in φ differ drastically between the two cases though. In this ordered angle case, we find some interesting grouping of the scaled angles.
Initially the bobs oscillate with an amplitude that depends sinusoidally on their position in φ, similar to the previous discussion in the x picture but with four nodes where the φ i remain relatively stationary. Once again this behavior could be thought of as a standing compression wave, but in φ i it has a wave length of π whereas in the x picture it had a wavelength of 2π. As the system evolves, all φ i slowly begin to shift towards the nodes of this standing wave until there are four clusters of the angles. After some time, the angles begin to re-distribute themselves randomly about the domain. The distribution of φ i in these three regimes is summarized in three histograms shown in Fig. 5. Aside from the number of clusters, there are two primary differences between the clustering in φ and the clustering in x: (i) The clustering in φ only occurs when the angles are ordered in the method described above, whereas the dynamics in x look identical regardless of the distribution in x, presuming it is somewhat homogeneous about the domain. (ii) The clustering in φ is a quasistationary state whereas the clustering in x exists for much longer times regardless of the system size. Since the clustering in φ is quasistationary, a properly prepared system could exist in the clustered angle state for an arbitrarily long time but only for large N . We can view the effect of increasing N and therefore the lifetimes of the clustered states by observing the order of the particle index as a function of time. In Fig. 6(a)-(c), we show that as N is increased, the time it takes for particles to fully mix increases. This is shown by plotting the indices on a color scale from 0 (blue) to N − 1 (red) along the horizontal axis as time is increased along the vertical axis. In Fig. 6(d), we show how the ordering of the particles changes at the very beginning of clustering for N = 100. Figure. 6 also shows that the compression wave is not quasistationary since it quickly reduces to the clustered state regardless of N .

V. CONCLUSION
Though the application of statistical mechanics and thermodynamics to systems with longrange interactions may not always be appropriate, we find that the canonical partition function improves our understanding of a system of pendula with long-range interacting bobs. Solving for the canonical partition function of the Hamiltonian in Eq. (15), we show that the equilibrium behavior in the φ coordinate frame is equivalent to the x coordinate frame, i.e. the HMF model. As we have argued that the Hamiltonian in Eq. (15) describes the behavior of the angles of repulsive or attracting pendulum bobs (see Fig. 1), then the proven equivalence of the canonical partition function of Eq. (15) and the Hamiltonian mean field model suggests that the Hamiltonian mean field model sufficiently describes the angles of a system of pendula in equilibrium. Ensemble equivalence between the microcanonical ensemble and the canonical ensemble is known for the Hamiltonian mean field model model [9] and because of this, the microcanonical simulations could be used to approximate equilibrium behavior. We find numerically that in the case of large system lengths, L, the dynamics of the system in φ resemble the dynamics of the Hamiltonian mean field model, equivalently the behavior of the system in x. Therefore for large system sizes of long pendula in equilibrium, the HMF model describes their dynamics and statistics.
In this paper we also briefly discuss two particular sets of non-equilibrium results. In one case, the system is initialized with small φ i so that x i are distributed relatively evenly throughout the x domain. This initial configuration essentially gives rise to the "repulsive" low temperature HMF model which exhibits interesting non-equilibrium behavior and is described in great detail by [19,21]. In the second case, in which φ i are ordered by their index i, we show there is a compression wave in φ, followed by clustering, and finally a mixed index state displaying no apparent order or structure. This is in contrast to the dynamics produced by a randomly distributed set of initial φ i which begins and then remains in a random disordered state. The clustering that can occur in φ is different from the clustering in x because it only occurs when the angles are initially ordered and because it is quasistationary; the lifetime increases with the number of particles in the system.