Long-range interacting pendula: A simple model for understanding complex dynamics of charged particles in an electronic curtain device

In this paper, we investigate the equilibrium and non-equilibrium properties of a model that shares several important characteristics with charged particles interacting in an Electric Curtain (EC) device. An EC comprises a periodic array of parallel electrodes, applied to each is an alternating electric potential. Depending on the applied potentials and the geometry of the electrodes, a wide variety of field structures above the plane of the electrodes are possible. The EC has multiple applications in the control and manipulation of small particles, but is under utilized in industry and science because of difficulties in predicting and understanding the particle dynamics. One particular challenge in understanding the dynamics is the many-body coulomb interactions. To better understand the role of the interactions, we study a one-dimensional analytically tractable model that encapsulates their long-range nature. Specifically, we study a Hamiltonian similar to that of the Hamiltonian mean field model but wi...


I. INTRODUCTION
In 1974, Masuda patented a device called the Electric curtain (EC) that would mitigate electrostatic paint particles from the surfaces in an electrostatic powder painting booth. 1 An EC is a device consisting of a series of parallel electrodes, often embedded in a dielectric surface, to which alternating electric potentials are applied. This device can create a large variety of electric field structures, resembling traveling or standing waves, above the plane of electrodes. Since Masuda's original proposal, the EC has been proposed for the manipulation and control of granular particles in many different applications including the separation of blood cells in aqueous solution, 2 separation of by-products from agricultural processes, 3 transport of toner particles in printing processes, 4 mitigation of charged dust from solar panels in extra-terrestrial environments, 5 and the separation of charged particles with different charge-to-mass ratios. 6 In general, the control and manipulation of charged granular materials can be challenging. Two prevalent obstacles are the complex particle dynamics of many interacting bodies, and the electrostatic adhesion of charged particles to surfaces. Electrostatic adhesion can make mechanical manipulations difficult because the particle(s) will stick to mechanical control surfaces resulting in a "sticky finger problem". In some applications, such as the mitigation of dust from a solar panel mentioned above, it is desirable to remove many particles from surfaces without any contact, unlike wiping which may abrade the surface. Since the electrostatic adhesion force is large for particles with large charge to mass ratios, simple devices that take advantage of particle's large charge-to-mass ratios are particularly useful in mitigation and control applications. Although the EC meets this criterion, it suffers in practical applications partly because of complex particle dynamics. For this reason, few of the proposed applications can be found in industry or science.
The lack of understanding alongside the potential promise of the EC has prompted a wide range of experimental and numerical studies 7-14 seeking to understand its particle dynamics and determine its average behaviors. Although a large body of work surrounds the EC, a comparatively small amount of literature addresses its study through simple models that allow for detailed investigations of its most basic dynamical features. In Refs. 15 and 16, we have presented investigations aimed at understanding some of the basic properties of an EC from a dynamical systems perspective. In a similar regard, in this paper we present a study of the statistical mechanics of a simple model with a focus on the interactions. We believe that the model presented here is not only important to understanding the behaviors of particles in an EC, but is also generally pertinent to spatially periodic systems with long-range interactions, a class of interactions to which the coulomb interaction belongs. Although the model presented here certainly oversimplifies the complexities of real EC systems, it is analytically tractable and its study serves two purposes: (i) to fundamentally improve understanding of long-range interactions in spatially periodic Published by AIP Publishing. 121, 154501-1 systems like an EC and (ii) to serve as a starting point for more detailed analytically tractable models of EC and related systems.
A. Motivating the model The simplest EC configurations are ones in which the alternating current (AC) potential applied to adjacent electrodes has a phase difference of p/2, called a 2-phase EC. The 2-phase EC is of particular interest because it can be easily fabricated with circuit board printing methods and requires only a single 2-phase AC power source. Above the surface of a 2-phase EC, the field resembles a standing wave with a wavelength equal to twice the electrode spacing. 1, 5 The dynamics of particles in a 2-phase EC can be classified into a few specific modes, the most interesting arguably being the "curtain mode." 12 In the curtain mode, particles can hover above the surface by oscillating between two electrodes. Impinging particles can become trapped in the curtain mode, effectively protecting the surface. When many particles are oscillating in curtain mode, it appears that the surface is behind an invisible waving curtain preventing them from settling on the "Electric Curtain".
Particles with small enough charge to mass ratios may not be levitated above the surface by the electric field. These particles can either enter the "hopping mode" 12 where they intermittently stick to the surface, sporadically freed by the electric field or collisions with other particles, or they may role/slide along the surface. In Ref. 15, we showed explicitly that the dynamics of a spherical particle constrained to the surface of a 2-phase EC are well described by a pendulum with a driven suspension point, a.k.a a Kapitza pendulum. For many particles constrained to the surface, one could think of multiple pendula with charged bobs interacting through the coulomb interaction. The pendulum analogy is particularly appropriate for many particles when they are constrained to a line, the pendulum angle becomes the position of the particle on a line. 15 Such a system with multiple particles can be experimentally realized using two ECs and a quadrupole trap as shown in Fig. 1. Particles are constrained to the axis of the quadrupole and the two ECs, when properly aligned, will have non-zero electric fields only along the direction parallel to the axis of the quadrupole. Such a configuration is a very simple way of creating an effective longitudinal pseudostatic (due to "small" frequency so the effect of the magnetic field on the charged particle can be ignored) standing wave. We call this a 1D-EC because of the restricted motions of particles in a quadrupole trap. However, in the right setup, particles could also pass each other during imperfect collisions. In a 1D-EC, particles will experience a time dependent potential, very similar to the potential experienced by the bob of a horizontal Kapitza pendulum. A 1D-EC is both a way to study coulomb interactions in time dependent potentials and might also be useful in certain transport and control applications. In the supplementary material section, we show a video of an experimental setup with six particles oscillating in the electric field.
An important property of the coulomb interaction is that it is long-ranged. Long-ranged interactions (as described in detail in Section I B) are known to exhibit interesting statistical and dynamical features. The goal of this paper is to study the equilibrium and non-equilibrium properties of a relevant long-range interacting model based on coulomb interactions. In this context, we restate the first objective, to understand these properties in a coordinate frame relevant to the EC, one where the electrodes serve as the reference points. This study, and model that is the focus of this study are also relevant to other systems with periodic potentials in which the important coordinates are the distance of particles from a set of potential minima or maxima such as particles on the surface of Faraday waves, or bosons in optical lattices. In this paper, we do not consider driving from an electric field so that we can isolate the role of the interactions. The interactions themselves are responsible for complex dynamics that must be understood in the appropriate coordinate frame.

B. Long range interactions
Systems with long-range interactions, such as the Coulomb interactions studied here, are unique in the field of statistical mechanics and thermodynamics, most notably due the absence of additivity. A system is additive when the total energy can be estimated by summing over smaller pieces of the system. 17 Directly related to the lack of additivity in longrange interacting systems is the fact that they are not extensive, i.e., the energy of a long-range interacting system increases much faster than the number of particles as one approaches the thermodynamic limit. 18 Therefore, caution is needed when applying the usual tools of statistical mechanics.
Long-range interactions are known to be a source of many interesting dynamical and statistical features. For example, an important property of some long-range interacting systems is that the canonical and microcanonical ensembles may be inequivalent. In such cases, the microcanonical ensemble exhibits a negative specific heat. 19,20 Additionally, some long-range interacting systems can contain interesting states that are "quasistationary." 21,22 Such states are different from metastable ones, which lie on local extrema of equilibrium potentials, in that they are stationary only in the thermodynamic limit. The lifetime of a quasistationary state increases with the number of particles in the system. Finally, non-equilibrium states in long-range interacting systems can create interesting dynamical features such as the spontaneous creation of macroscopic structures. 23 FIG. 1. Two electric curtains can create an effective pseudostatic longitudinal electric field standing wave along the axis of the quadrupole trap.
For some time, the primary motivation for the study of long-range interactions was to understand galaxies, galactic clusters, the general thermodynamic properties of selfgravitating systems, and mean field models. Interests have been reinvigorated since the observation of modified scattering lengths in Bose-Einstein Condensates (BEC) through the use of the Feshbach resonances. 24 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. 25 have 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 tabletop 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 inherent in understanding systems with long-range interactions have driven the development of solvable models that could provide new insights into the aforementioned phenomena. 17 One particularly significant model that is directly relevant to this study is the Hamiltonian Mean Field (HMF) XY spin model, 23 often written in the form where h i is the angular position of the i th particle (spin), c sets the scale of the interactions, p i is its conjugate angular momentum, and we set the mass to unity. The HMF model is generally used to describe two different classes of systems: (i) a mean field XY classical spin model, and (ii) a onedimensional periodic system of itinerant particles with longrange interactions. Although the connection between Eq. (1) and the latter system could be thought of as contrived given the simplifications under which the model is realized, it has provided useful insights into how non-neutral plasmas and self-gravitating systems behave. 23 Some of the essential features of the EC can be captured by a simple model of N pendula with charged bobs, meant to represent the positions of charged particles. Within this approximate description, the HMF model of Eq. (1) provides a framework to study a system with explicit index dependence, i.e., where the distance of each pendula bob from individual pivot points is physically meaningful, measuring the distance from the potential minima of a field created by a series of electrodes.
In this paper, we study such a model of long pendula with flat charged bobs undergoing small oscillations and having parallel planes of rotation as illustrated in Fig. 2. It is shown that this model can be related to the HMF model through a coordinate transformation that introduces indices of the particle labels. The index dependence in the Hamiltonian, which will be described in detail in Section II, 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 angular 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 angular coordinate frame clustering is only found for the initially ordered angles and, unlike the biclusters, these are clearly quasistationary states that can persist for a length of time which goes to infinity as the thermodynamic limit is approached. 17 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 Ref. 26 on a model sometimes called the HMF a-model where a 1=r a ij dependence between the classical spins is introduced, 17,27-30 with r ij being the distance between the i th and j th spins on a lattice. Although the physical motivations behind studying these various models may be different, it is interesting that their equilibrium behavior is nearly or almost equivalent. 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. 2, 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 where d is the distance between the pivots of neighboring pendula and ' is the length of each pendulum. The small h regime makes the problem one-dimensional in x. We choose periodic boundary conditions and rescale the system by 2p/Nd, so that making the total system length a dimensionless 2p, where N is the number of particles in one period. We refer to a periodic space with length 2p as a unit circle. The position of the i th particle (bob) can then be written as For reasonable choices of ' and d where '/d ( N, the second term is suitably small such that the Hamiltonian can be written with terms that are quadratic in h. However, we are primarily interested in a regime where '/d ! 1 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 right hand side of Eq. Given the choice of large '=d; / i can take any value in the range [0, 2p). This is only true because '/d is large, the h i s remains small. In terms of / i , the positions can thus 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 simply 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 one-dimensional ring. 23,31 The distribution of the bobs is such that the density, q(x), is given by where subtraction of the constant 1/2p is necessary to produce a meaningful expression for the potential U 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 onedimensional 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 c contains the particle (bob) charge or mass and becomes the interaction strength in the Hamiltonian, and the zeroth-order term in the Fourier series canceled the constant neutralizing background. The most important simplification in this paper is truncation of the Fourier series used to represent the delta function at n ¼ 1. Antoni and Ruffo 23 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. The truncation of the sum is physically 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.

C. Solving Poisson's equation
where the determination of the constant c 1 requires the examination of the physical picture. A sensible requirement is that when all of the bobs are hanging at their equilibrium positions, directly below their pivot (all / i ¼ h i ¼ 0), the net force experienced by any bob is zero. This is a valid requirement if the interaction is attractive or repulsive and is satisfied by c 1 ¼ 0. Integrating once more to obtain the potential yields To determine c 2 , we stipulate that if all / i ¼ 0, then U(0) ¼ 0. This condition is satisfied by c 2 ¼ 0 and we can now write the potential energy of the j th particle as D. The Hamiltonian Using Eq. (9), the Hamiltonian can be written as where the mass of the bobs has been set to unity, and the 1/p coefficient in the original potential energy has been absorbed into interaction strength c. 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. This 1/N scaling is known as the Kac prescription 32 and 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. 17 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 particles on a ring with longrange repulsion or attraction, as well as describing a classical XY spin model. The h i played the role of either the position of the i th particle on the ring or the orientation of the i th spin. Thus far, the rescaled angle / i ¼ 2p'h i =Nd 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 language of Eq. (10), the potential energy of the i th and j th spin pair depends 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 and define expressions in terms of cos / i and sin / i , the horizontal and vertical components of a fictitious magnetization m i ¼ ðcos / i ; sin / i Þ.

III. EQUILIBRIUM PARTITION FUNCTION
In this section, we solve for the canonical partition function, in the / coordinate frame using the Hamiltonian in Eq. (10) and show that, in equilibrium, the HMF model describes the angles 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 Considering the coefficients in the Hamiltonian cos½2p ði À jÞ=N and sin½2pði À jÞ=N to be matrices with components C ij and S ij , respectively, the interaction can be written in the simplified form where m > l ¼ ðm 0;l ; m 1;l ; …; m NÀ1;l Þ with l ¼ x, y. The coefficient matrices C and S are circulant, a special kind of Toeplitz matrix where each subsequent row is a cyclic permutation of the row above or below it. They have the general form a 1 a 2 a 3 … a NÀ1 a N a N a 1 a 2 … a NÀ2 a NÀ1 a NÀ1 a N a 1 … a NÀ3 a NÀ2 . .
and can be diagonalized by a unitary matrix. In particular, C and S are simultaneously diagonalizable as it is easily shown Starting with the second term, ÀSC ¼ S > C > ¼ ðCSÞ > which is found by arguing that C is symmetric since the cosine is an even function and does not change under the exchange of i and j, whereas S is odd because sine is an odd function and does change sign under exchange of i and j. The commutation becomes ½C; S ¼ CS þðCSÞ > . Also, an odd function multiplied by an even function results in an odd function so the entire matrix CS is odd. Therefore, ðCSÞ > ¼ ÀCS bringing us to the final expression [C, S] ¼ CS -CS ¼ 0.
We have shown that C and S can be simultaneously diagonalized by a complex unitary matrix U which is known for circulant matrices and is called the Fourier Matrix.
The matrices C and S can be rewritten as C ¼ U † K C U and S ¼ U † K S U, where K C and K S are diagonal matrices with elements that are the eigenvalues of C and S, respectively, which we denote as k C i and k S i . From here on we label the indices i from 0 to N -1. Inserting the spectral decomposition of C and S, Eq. (12) becomes We may return from the matrix to index notation using the following relations: where X is not to be confused with x and the interaction is now given by where kXk ¼ XX Ã . The inclusion of the eigenvalues k C and k S simplifies the Hamiltonian further as shown in Appendix A and leads to This representation must be further modified before the partition function can be found and by splitting the Fourier matrix U into its real and imaginary components, U ik ¼ a ik þ ib ik , given by We would like to remind the reader that i in the above equations refers to the index of each particle. Noticing that a 1k ¼ a ðNÀ1Þk and b 1k ¼ Àb ðNÀ1Þk , we write the configurational partition function as where b ¼ 1/k B T. The Hubbard-Stratonovich transformation 33,34 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 switching the order of integration, we find The integration can be performed using the identity where and we have used Eqs. (20) and (21) to obtain the second line. It is convenient to make a change to polar coordinates by introducing z ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffi ffi z 2 1 þ z 2 2 p , following which the partition function can be written as Equation (26)  The integration over z in Eq. (26) can be performed using the saddle point approximation. The rescaled free energy per particle is given by 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. (27) results in the equation which can be solved self consistently for z for a given value of b. For b < 2, the system is paramagnetic but for b ! 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 parameter is the total magnetization M ¼ ð1=NÞ Showing that the canonical partition function in the / coordinate frame model and the HMF model are equivalent necessitates a more detailed discussion of the equilibrium behaviors in the / frame. Campa et al., 17 in their review of the HMF model, rigorously show ensemble equivalence between the canonical and microcanonical ensemble of the HMF model. It has also been shown that a second order phase transition in the canonical ensemble implies ensemble equivalence. 35,36 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, 23,37,38 begging the question: does a large microcanonical simulation of Eq. (10) 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 ensemble equivalence property of the model, or lack thereof. The omission was in the arbitrary scaling of x which we now rectify.
We introduce the dimensionless parameter L which generalizes the scaling in Eq. (2) to making the position of the ith particle and changing the definition of / i to / i 2pL'h i =Nd. 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, what we find is that if L ) 1, then the simulations in / closely reproduce the dynamics of HMF model simulations (dynamics in x). Therefore, for L ) 1 a large microcanonical simulation of Eq. (10) approximates the expected equilibrium behavior, and the index-dependent model in equilibrium with a heat bath can be described by the HMF model. Alternatively, the coordinate frame inequivalence is most extreme for small L. These numerical results were found using initial conditions that are randomly distributed / i about the domain [-Lp, Lp). It should be stated that for the rest of this paper we work with L ¼ 1 because we are interested in cases where the / coordinate frame is markedly different than the x coordinate frame.

IV. NON-EQUILIBRIUM RESULTS
We now focus on the non-equilibrium dynamics of the system of pendula. For this discussion, we will only work with initial conditions where all initial _ x i ¼ 0 and _ / i are set to zero. We choose to investigate the lowest energy behaviors of a few specific sets of initial positions x i and / i and do not discuss the temperature dependence of our results in this paper.
We begin with an initial configuration where all interacting bobs are set to random small displacements from / i ¼ 0. Specifically, we initialize the i th pendulum angle, / i , randomly in the range [-p/N, p/N) (all _ / i ¼ 0). If 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 [2pi/N -p/N, 2pi/ N þ p/N) it is possible to make some general statements about the dynamics of this configuration in x using the equations of motion found in the / coordinate frame. Expressing the interaction with terms that are quadratic in / yields and the equations of motion for the ith particle are found to be where l 1 ¼ / j cosð2pi=NÞ and l 2 ¼ / j sinð2pi=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 hl 1 i and hl 2 i are known, then the initial dynamics of the system are elucidated by Eq. (32), but in the case of randomly initialized / i the hl 1 i and hl 2 i 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 (32) 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 i ! x, so the force takes the form Therefore, hl 1 i and hl 2 i 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. 3, we show a fit of the force as a function of x using Eq. (33) as well as the actual force calculated for an example set of initial conditions. This result holds even for initial positions in x which are not ordered by their indices as long as they are homogeneously distributed about the x domain. The equations of motion in the HMF do not depend on the indices but we arrive at Eq. (33) from the / coordinate using the assumption of ordered indices. The domain in Fig. 3 can be split into two pieces (independent of l): 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 hl 1 i and hl 2 i 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 2p. 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 e et al. by analysing the Vlasov equation. They find that the initial compression wave (referred to by a different name) creates an effective double-well potential giving rise to the bicluster. 39 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. 40 Depending on hl 1 i and hl 2 i, 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 begins to spread out over the full domain [0, 2p) and continues 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. ½2pi=N À p=N; 2pi=N þ p=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 p, whereas in the x picture it had a wavelength of 2p. 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 homogeneously about the domain. The distribution of / i in these three regimes is summarized in three histograms shown in Fig. 4. 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 / appears to be quasistationary state, as we will show, whereas the clustering in x exists for much longer times regardless of the system size. We find that the clustering in / is directly related to the order in the indices. We show in Fig. 5 that the lifetime of the order in the indices and therefore the clusters appear to depend on N. 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. Initially, some mixing of the indices happens quickly regardless of N as the particles first collect at the nodes of the standing compression wave. The initial mixing is shown in Fig. 5 for N ¼ 100, after which the system retains some order in the indices. The order observed after the initial mixing corresponds directly to the presence of the clusters. In fact, it may be possible that the clusters and their stability could be described by this rearrangement of the indices because of the index dependence in the potential energy of the system (Eq. (9)). In Figure 5, N ¼ 15, N ¼ 30, and N ¼ 60 are all shown for the same time axis. For N ¼ 15, no visible order exists whereas for N ¼ 30 some order can be seen up to t $ 6000. For N ¼ 60, the order is much more clear and can be seen past t $ 8000. As N is increased, the time it takes for particles to fully mix increases.
In astrophysics, there is often an abrupt change that occurs early in a systems dynamics, known as violent relaxation, before a slower evolution takes place towards equilibrium. 41 The time scale for the violent relaxation is independent of N and should be recognizable as a drastic change in a systems total potential and kinetic energy. In Fig. 6, we show the total potential energy, U tot , as a function of time for two different values of N, N ¼ 50 and N ¼ 100. Initially, the potential energy only fluctuates about a constant value with small amplitude oscillations. For N ¼ 50 around t ¼ 200, and for N ¼ 100 around t ¼ 500, the behavior changes drastically and U tot starts undergoing much larger fluctuations with a much different time scale than the small fluctuations that continue to occur. The small amplitude oscillations of U tot occur too rapidly to be seen clearly in Fig. 6 and are only perceivable by the apparent thickness of the line in the plot which would be much thinner in their absence. The dramatic shift in the behavior of U tot is similar to a violent relaxation in that it immediately precedes the more slowly evolving clustered state, but Fig. 6 also clearly shows its dependence on N. In this case, the N dependence is a consequence of the inhomogeneous initial conditions and the particle interaction's 1/N scaling. Because the dynamical evolution of the initial compression wave slows for increasingly large values of N, so too does the initial rearranging of the indices which is responsible for the abrupt change in potential energy.

V. DISCUSSION
Although the application of statistical mechanics and thermodynamics to systems with long-range 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. (10), we show that the equilibrium behavior in the / coordinate frame is equivalent to the x coordinate frame, i.e., the Hamiltonian mean field model. We have demonstrated that the Hamiltonian in Eq. (10) describes the behavior of the angles of repulsive or attracting pendulum bobs (see Fig. 2) and shown the equivalence of the canonical partition function of Eq. (10) and the canonical partition function of the Hamiltonian mean field model. This equivalence suggests that the partition function of the HMF model describes the statistics of the pendula angles in equilibrium. Furthermore, an important subtlety in the equilibrium behavior of the pendula angles (/) as compared to the equilibrium behavior of the positions (x) is found. For the HMF model, ensemble equivalence between the microcanonical ensemble and the canonical ensemble is known, 17 and because of this, the microcanonical simulations could be used to approximate equilibrium behavior. The subtlety in the equivalent equilibrium behaviors lies in the scaling of the system length in the coordinate transformation from x i to / i . In the case of large total system length, L, the dynamics of the system in / resemble the dynamics of the regular Hamiltonian mean field model. Therefore, for large system sizes of long pendula in equilibrium, the HMF model describes their dynamics and statistics.
In this paper, we have also briefly discussed two particular sets of non-equilibrium results for L ¼ 1. 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 nonequilibrium behavior and is described in great detail by Refs. 39 and 42. In the second case, in which / i are ordered by their index i, we show that 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 it appears to be quasistationary; the lifetime increases with the number of particles in the system.
The understanding of the HMF model has been very much improved through the study of the Vlasov equation which can approximate the dynamics in the N ! 1 limit. 23,39,[43][44][45] For this reason, we should comment on the difficulties in applying it to the model studied in this paper. A foundation for the use of the Vlasov equation to study the HMF model, and other pairwise long-range interactions, is a theorem due to Braun and Hepp. 46 This theorem addresses potentials of the form Uð/ j À / i Þ, whereas the model studied in this paper, with its explicit dependence on the indices, is of the form Uð/ j À / i ; j À iÞ. For the case when all particles are initially ordered by their indices, the index dependence in the potential can be approximated by a spatial dependence, especially in the N ! 1 limit, which is how Eq. (33) is derived. In this case, the Vlasov dynamics should be a good approximation for the model. However, because the indices mix during the violent relaxation, the dynamics would quickly depart from the Vlasov dynamics. As the clustered states and the mixing of the indices are our primary interests, we do not investigate the initial behavior of the model using the Vlasov equation.
We show that one of the interesting features of an EC is that, for particles in a simple curtain mode, the particles can be considered as have a pair-wise index dependent potential energy. Although we have not discussed the effects of a driving electric field, this view of the particle iterations maybe a useful tool in exploring curtain mode oscillations. It may be interesting to construct an EC system where these effects are clearly visible as a way to explore dynamics in a system with novel index dependent interactions. It is also interesting to speculate on other physical systems that might exhibit a similar potential energy. One possibility would be to consider a trapped ultra-cold system of quasi-one-dimensional boson tubes consisting of polarized electric dipole moments perpendicular to the tubes. 47 An index dependence in the dipole-dipole potential energy between tubes would be found if each tube was located at a step of a spiral staircase wound around a central axis. The labeling of the indices would start at one end of the boson tube "staircase" and increase as the staircase is climbed. This particular geometry would result in the index dependent phase in the potential energy because the tubes are offset from one another in two ways, the direction along the axis of the staircase and the rotation about the axis. In this system, the / coordinate frame would correspond to the angular position of the dipole moments looking down a given tube.
In conclusion, we have shown that the angles of a system of pendula, consisting of charged or self-gravitating bobs, can be described by a coordinate transformation of the Hamiltonian mean field model that introduces indices of the particle labels. We believe that this model adds to the understanding of the behavior of particles in an electric curtain and other systems with long range interactions and spatially periodic potentials. For this model, we have presented results for the equilibrium and non-equilibrium behaviors of the pendula angles and find that in equilibrium with a heat bath, the pendula angles are described by the Hamiltonian mean field model. In the non-equilibrium case, for a particular set of initial conditions, we find spontaneous ordering of the angles into four clusters. These angle clusters appear to be quasistationary and only appear to stabilize after some initial reordering of the particles that happens quickly compared to the lifetime of the clusters. More work is required to understand the clustering in the angular coordinate frame and it would be interesting to study explicitly how the stability of the four clusters depends on the ordering of the particle indices by varying the initial configurations of pendula bobs. A better theoretical understanding of the dynamics of particles with long range interactions may be essential to improve the viability of the proposed electronic curtain devices and applications.

SUPPLEMENTARY MATERIAL
See supplementary material for a video of six particles oscillating in an experimental setup similar to that shown in Fig. 1. The particles are glass beads with diameters of approximately 35lm. They are confined to nearly onedimensional movement by a quadrupole trap and two arrays of parallel electrodes (electric curtains) produce the field causing the oscillations. In this Appendix, we solve for k C and k S and show how they simplify the expression of the Hamiltonian in Eq. (18). Looking at the form of a circulant matrix shown in Eq. (13) reminds us that the elements can be defined with a single label which we write as where ' ¼ 0, 1, 2,…, N -1. The eigenvalues, k 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 k A j ¼ P NÀ1 '¼0 a ' expði2p'j=NÞ and ' ¼ 0, 1, 2,…, N -1. Therefore The above expressions of the eigenvalues show that C and S each have only two non-zero eigenvalues corresponding to j ¼ 1, N -1 given by k C 1 ¼ k C NÀ1 ¼ N=2 and k S 1 ¼ ðk S NÀ1 Þ Ã ¼ iN=2. Using these, we arrive at the simplified Hamiltonian in Eq. (19).