Effect of Nonreciprocal Forces on the Stability of Dust Clusters

Results are presented from studies of the stability of the plane dust clusters in the form of a regular polygon with the number of particles from two to five. It is assumed that the particles are placed in the plasma consisting of Maxwellian electrons and a directed flow of cold ions. It is shown that, in such clusters, the oscillatory instabilities can develop along with the aperiodic instabilities. The ranges of plasma parameters are determined, within which the oscillatory instability of the five-particle cluster becomes saturated at the weakly nonlinear stage. As a result, the cluster forms a time crystal, which can be a chiral crystal.


INTRODUCTION
In the physics of dusty plasma [1][2][3], the ensembles of a small number of macroparticles levitating in the near-electrode layer of the gas-discharge plasma are called the clusters [4]. The equilibrium levitation height is determined by the balance between the gravity and electrostatic repulsion of charged particles from the lower electrode. In the horizontal direction, the particle confinement is ensured by the special electrode profiling. Thus, each individual particle is localized in the asymmetric potential well characterized by the oscillation frequencies (in the vertical direction) and (in the horizontal direction). The distinctive feature of particle interactions is the fact that due to the presence of the directed ion flow in the near-electrode region, the electrostatic potential in the vicinity of the charged particle turns out to be an asymmetric function of coordinates (the z axis is directed vertically). As a result, the interparticle interaction turns out to be nonreciprocal, i.e. the third Newton law is violated. This, in particular, results in the occurrence of correlation between the displacements of individual particles in the vertical and horizontal directions.
Theoretical studies of the dynamics and stability of the plane plasma clusters (for example, [5][6][7][8][9]) were performed under the assumption that the potential of particle interaction is the Coulomb or some other model potential that does not take into account the effect of nonreciprocal forces. This article is the continuation of a series of studies [10][11][12][13][14], in which the dynamics of the two-dimensional ensembles of particles was investigated using the interaction potential calculated numerically. The approach being developed makes it possible to fully take into account the nonreciprocity of interparticle forces.
The dynamics of clusters in the form of a regular polygon, consisting of a small number ( ) of particles, is considered. As is known, in the quasi-twodimensional plasma crystals, the nonreciprocity of interparticle forces leads to the development of the instability of coupled waves [3]. The nonlinear theory of this instability [14] shows that a time crystal is formed in a certain range of plasma parameters, i.e., in the ground state, the dust particles oscillate near the equilibrium positions. Such instability, which, in this article, is called the oscillatory instability, can also develop in the case of finite clusters.
The article is organized as follows. In Section 2, the model used is briefly described. Next (in Section 3), the linearized equations of motion for the N-particle cluster are reduced to the block-diagonal form. In Section 4, the eigenfrequencies of different clusters are investigated, and the ranges of parameters are determined, within which one or another type of instability develops. Finally, in Section 5, the nonlinear theory is developed for the oscillatory instability of a five-particle-cluster which is of the most interest.

MODEL
An ensemble is considered consisting of N particles located in the external potential well . It is believed that the particles have the same constant charges Q and are placed in the plasma consisting of the Maxwellian electrons where is the flow velocity. In this case, the dielectric constant of plasma is equal to , where is the ion plasma frequency, is the electron Debye radius, and is the infinitely small collision frequency. The interaction is isotropic in the plane , i.e., potential (1) has the following form: ( ). The presence of the directed ion flow results in the fact that , i.e., the interparticle forces are nonreciprocal.
We use the dimensionless variables with the length scale ; the interparticle forces are normalized to , and the time scale for particles with the mass is . When we use these variables, potential (1) is a function of the distance and also of the only one parameter characterizing the plasma, which is proportional to the Mach number of the ion flow, , where are the equilibrium densities of electrons and ions, and is the ion mass. The derivatives of the potential are calculated using the numerical methods described in [13].
In equilibrium, the particles are located at the vertices of a regular polygon The general equations of motion are as follows: where . Expanding the equations of motion in powers of the small deviations , in the zero-order approximation, we obtain the following equilibrium conditions: (4) where we used the following notations for the derivatives of the interaction potential: (5) and is the distance from the particle with number j to the particle with number 0.
In further consideration, the trap frequencies and the Mach number M are considered as the external control parameters. The equilibrium cluster radius is calculated using the first of Eqs. (4); it depends on и M. Several examples of how the parameter R depends on M at the fixed frequency are shown in Fig. 1. We note that at potential (1) turns out to be the attractive potential at considerably large distances [11]. The shape of the potential well bottom in the coordinates (R, M) is shown in Fig. 1 (curve 1). In this case, there may be no external confinement in the transverse direction, and the cluster radius reaches its maximum value at . At , the cluster radius indefinitely increases with decreasing frequency .

LINEARIZED EQUATIONS OF MOTION
In the first order of expansion of the general equations of motion (3) in powers of the small deviations from equilibrium, we obtain the following equations: In all, there can be oscillations with frequencies determined by the eigenvalues of the force matrix . The analysis of the linearized equations of motion can be greatly simplified, if we perform the decomposition into irreducible representations of the symmetry group of the cluster [5]. For the hori-   zontal and vertical displacements, this decomposition looks differently: (6) where . Using the formula of finite geometric progression, it is easy to obtain the following inverse transformations: At the perturbation with corresponds to the homogeneous cluster expansion or compression in the xy plane, and the perturbation with results in the cluster rotation around the z-axis. Since the displacements are real, the transformed parameters meet the obvious correlations and . After using transformations (6) and (7), the linearized equations of motion (3) can be rewritten in the block-diagonal form: where the matrix has the following form: The matrix elements included here can be writted as the finite real sums: In the general case, matrix (8) does not have any symmetry properties. However, in the absence of nonreciprocal forces, the derivatives are equal to zero, , , and the matrix becomes the Hermitian one. In this case, the equations for the displacements along the z-axis and in the xy plane can be uncoupled and studied independently.
The use of transformations (6) and (7) makes it possible to considerably simplify the problem of searching for the oscillation frequencies and analyzing the stability of dust clusters. All oscillations are divided into the groups of three oscillations; each group is characterized by the number l, which corresponds to the wave vector in the infinite systems. We denote the eigenvalues of the matrix (8) as ( ). These frequencies are the roots of the cubic polynomial (I is the unity matrix), and, as can be seen from the structure of matrix (9), the polynomial coefficients are real. It follows from explicit expressions (9) that , so the relation is satisfied, i.e., all oscillations with are doubly degenerate. This degeneracy occurs due to the fact that some waves propagating in different directions of the ring structure should have the same frequencies. The oscillations with and (for even numbers N) are not degenerate. For this reason, for each cluster consisting of particles, it is necessary to study the eigenvalues of the matrix for . For each oscillation mode, the displacements of particles in the cluster are determined by the eigenvectors of the matrix and can be calculated using relations (6).
Two different scenarios are possible for the loss of stability when the external parameters change. First, one of the roots of the characteristic polynomial can change its sign, which corresponds to the development of the aperiodic instability. Second, two real roots can merge and form a pair of the complex conjugate roots. In this case, we deal with the oscillatory instability similar to the instability of coupled waves developing in the infinite plasma crystal.
First of all, we notice some general properties of matrix (8). At , the eigenfrequencies can be obtained explicitly for an arbitrary number N. First, there is the neutrally stable mode with , which corresponds to the rigid rotation of the cluster without changing its shape ( , , ). Another oscillation with ( , , ) corresponds to the vertical oscillations of the center of mass of the cluster without changing its shape. Finally, the third mode has the following frequency: (10) For this oscillation mode, the particles are displaced along the straight lines (Fig. 2). In this case, the horizontal components of the displacements are parallel to vectors (2), and all their vertical components are equal. Although the frequency (10) does not depend on , the slope of the displacements is determined just by . In all the cases that have been studied, the square of frequency (10) turns out to be positive, and no instabilities arise associated with the excitation of the symmetric oscillations. In what follows, the oscillations with are not specially discussed.
In the case of oscillations with , one of the frequencies of a cluster consisting of particles coincides with the frequency of the trap , i.e., the characteristic polynomial can be written as a product . In this case, the shape of the cluster does not change; all particles are in the xy plane, and the center of mass oscillates in an arbitrary direction with the frequency . The frequencies and polarizations of the other two oscillation modes with depend on the number of particles in the cluster.

Two Particles
In the case of a cluster consisting of two particles, the oscillation with the frequency turns out to be doubly degenerate. The frequency of the remaining oscillation is equal to , and the particle displacements are proportional, and . The condition that the quantity should be positive determines the frequency range of the cluster stability: . An example of the (M) dependence constructed using potential (1) is shown in Fig. 3 (solid line). As the frequency decreases to the frequency, there occurs a tendency to the formation of the vertical chain of particles; in this case, due to nonreciprocity of the interparticle forces, both particles are displaced in the same direction along the -axis.

Three Particles
In the case of three particles, the only nontrivial mode is characterized by the number . The roots of the quadratic equation are as follows: where . Since the parameter D is positive, , both solutions (11) are real, but one of them can become negative with decreasing frequency . The corresponding critical frequency is shown in Fig. 3 (dashed curve). When the aperiodic instability develops, the vertical displacements of two arbitrary particles are directed opposite to the displacement of the third one, for example, , , and . In this case, the horizontal displacements have the same sign. Thus, during the development of the aperiodic instability of the three- ( 0 ) x z r = , ,−

Four Particles
In contrast to the examples considered above, the dynamics of a four-particle cluster is determined by the potential derivatives at two distances: and . For this reason, the explicit formulas for the frequencies turn out to be more cumbersome, and in further consideration, we will restrict ourselves to a qualitative discussion of the cluster stability regions.
In the case , the eigenfrequencies of oscillations that remain after we have considered the solution with are determined by the roots of the quadratic equation . These roots will be the complex ones, if the following two conditions are satisfied: (12) and (13) where (14) At the boundary of the domain defined by inequalities (12) and (13), two real roots merge, and the frequencies of coupled oscillations are determined by the following relations: For the coupled oscillations at the boundary of the instability domain, the particle trajectories are shown in Fig. 4. The particles move circularly around the equilibrium position in the plane tilted at some angle   to the xy plane. As already noted, the oscillation frequencies with and coincide, but in the latter case, the direction of particle motion is opposite.
The similar dynamics is observed for oscillations with . There is the stable oscillation with the frequency , for which the particle displacements shown in Fig. 5a are in the horizontal plane and do not depend on the external control parameters. The two remaining coupled oscillations can become unstable, if inequalities (12) and (13) become satisfied. In this case, the critical values of the parameter will change as follows: (16) The frequencies of coupled oscillations at the boundary parameters of the oscillatory instability development are as follows: (17) and the displacements of particles occur along the straight lines inclined at the same angles to the xy plane (Fig. 5b).
If one of inequalities (12) or (13) is violated, the roots of the characteristic equation become real. However, for both and , the square of one of the frequencies can become negative, which results in the development of the aperiodic instability. For the fourparticle cluster, the resultant stability diagram is shown in Fig. 6. At the boundary of the region , the aperiodic instability of the mode with develops; an example of the particle displacements is shown in Fig. 7. The region of the aperiodic instability development for the mode with (the region in Fig. 6) is entirely inside the region .
The regions of the oscillatory instability development are small. They are shown on the enlarged scale in Fig. 6 on the right. These regions correspond   to the same domain of variation of the parameter M, determined by the first inequality (12), and partially overlap.

Five Particles
In the case of five particles, the cluster dynamics is determined by the derivatives of the interaction potential at two distances: and . The explicit expressions for the frequencies turn out to be very cumbersome, but can be easily calculated using the numerical methods.
Just as in the case of the four-particle cluster, it is necessary to study the modes with and . The cluster stability diagram is shown in Fig. 8. At the boundary of the region , the aperiodic instability of the mode with develops; in this case, two arbitrary neighboring particles leave the xy plane in the course of their displacement, and the other three displace in the opposite direction. The region of the aperiodic instability development for the mode with is entirely inside the region . In contrast to the case of the four-particle cluster (Fig. 6), the regions corresponding to the development of the oscillatory instability occupy a considerable part of the stability diagram. For both types of the coupled oscillations with (the regions in Fig. 8), the particles move in circular orbits in the vicinity of the equilibrium positions (Figs. 9 and 10).
Their motions differ in the amplitudes and phases of particle displacements. It turns out that both regions consist of two parts. The two parts of the region are separated by a narrow "slit" corresponding to the range ; therefore, it is rather difficult to reach the regions under conditions of the poorly controlled changes in the control parameters.

Open-Boundary Clusters
As already noted, in the subsonic or pure ion flow with , in the horizontal direction, interaction potential (1) turns out to be either the attractive potential (at the considerably large interparticle distances), or the repulsive one (at small distances). In this case, the clusters can exist even in the absence of the radial confinement, . The structure and dynamics of such clusters have some specific features.
The bottom of the potential well is located at the distance determined as the solution to the equation . The shape of the potential well in the (R, M) coordinates is shown in Fig. 1 (curve 1). : 1 0 1 Fig. 7. Particle displacements at the boundary of the region in Fig. 6; , , and .   In accordance with conditions (4), at , in the case of the two-particle cluster, the maximum radius is , and for the three-particle cluster, it is . The formulas for the oscillation frequen- cies given in Sections 4.1, 4.2 are also applicable in the case of . When the parameter decreases to a certain critical value (Fig. 11), the aperiodic instability develops. An interesting feature of this instability is that, in contrast to the case of , at , the vertical displacements of the particles vanish. In this case, at the boundary of the region of the aperiodic instability development, the two-and three-particle clusters begin accelerating in the horizontal direction without changing their shape, that is, they behave like molecular motors.
Because of equilibrium condition (4), in the case of four particles located at the vertices of a square, at , the inequalities are satisfied, i.e., the nearest particles repel each other, and the opposite ones attract each other. Obviously, such a configuration is unstable even in the case of twodimensional motion ( ), and the square tends to transform into a rhomb. Similar situation arises when the cluster has the form of a regular pentagon.
Thus, in the absence of the radial confinement, the clusters in the form of a regular polygon may exist only when the number of particles is two or three.
In conclusion of this section, we note that the clusters with a large number of particles ( ) always turn out to be unstable. Apparently, this is a general law: all known examples indicate that the two-dimensional clusters in the form of a regular polygon exist only for . For example, in the case of the Coulomb interaction and the two-dimensional motion ( ) of six particles, the only stable configuration has the form of a regular pentagon with an additional particle in its center [5].

NONLINEAR STAGE OF OSCILLATORY INSTABILITY
As can be seen in Figs. 6 and 8, when varying the external control parameters , the boundary of the region of the oscillatory instability development can be most easily reached for the mode of the five-particle cluster (the region in Fig. 8). For this reason, we restrict ourselves to discussing the nonlinear dynamics of the cluster with . The standard Krylov-Bogolyubov method is used to derive nonlinear equations. We set , where is a small parameter. When approaching the boundary of the region (Fig. 8) from above, the parameter δ vanishes, and at , the instability develops. At the lower boundary of the region the reversed situation occurs: the instability develops at δ > 0.
After substituting into the initial equations of motion (3), they are expanded in powers of up to the third order, and then transformations (6) Fig. 9. Particle trajectories at the boundary of the region of the oscillatory instability development; , , , , and .    (7) are used. As a result, we obtain equations in the following form: (18) where the matrix is given by formula (8). The matrices are the quadratic functions of , and the matrices , in addition to the cubic expansion terms, contain a linear term proportional to δ.
We try the solution to Eqs. (18) in the form of the following expansion: where is the slow time. In expansion (19), we set and is the real slowly varying amplitude, i.e., the first term of the expansion describes the slow rotation of the cluster around the z-axis. The complex amplitudes correspond to the perturbations propagating in different directions of the ring structure. The frequency is a doubly degenerate root of the dispersion equation at , and is the corresponding eigenvector, i.e., . The corrections and to be determined should be real, which follows from the correlations . After substituting expansion (19) into Eqs. (18), the zero term of the expansion in powers of vanishes. The next term of the expansion makes it possible to express the corrections in terms of the derivatives of the amplitudes and their quadratic combinations. The existence condition for the solution for the corrections leads to the following set of equations for the amplitudes: The explicit expressions for the coefficients α, β, and are too cumbersome, but they can be calculated using the numerical methods. We note that the β coefficient is always positive at the upper boundary of the region in Fig. 8 and always negative at its lower boundary. The rest of the coefficients can change sign when changing the parameters M and .
The physical meaning of Eqs. (20) and (21) is quite clear. In the linear stage, the presence of the terms proportional to β leads to (depending on the sign of δ) either an exponential increase in the amplitudes or the = ⋅ + ε + ε , Thus, during the development of oscillatory instability of the cluster, different stationary states arise. The parametric interaction can result in the fact that the amplitudes of the perturbations are equal and the rigid rotation of the cluster does not occur. In this case, the cluster particles oscillate along the straight lines, as it is shown in Fig. 2; however, the directions of their trajectories depend on the external parameters.
It is also possible that only one of the perturbation amplitudes is nonzero; then the particles circularly move around the equilibrium positions (Fig. 10), and the rigid rotation of the cluster around the z-axis occurs at the constant angular velocity. Which of the amplitudes is non-zero and the direction of rotation are determined by random factors. In this case, the cluster becomes the chiral one.
The coefficients and at the upper boundary of the O 2 region (Fig. 8) as functions of the Mach number M (calculated numerically using potential (1)) are shown in Fig. 12. At , the coefficient is negative, , i.e., in accordance with inequalities (25), there are no stationary states and the instability development is explosive. In particular, in the entire left part of the region (Fig. 8) this coefficient is negative.
In the ranges and M > , inequality (29) is true, i.e., the instability becomes saturated in the weakly nonlinear stage. In the stationary state, the clusters become chiral, and the amplitude of the oscillations is determined by expressions (28). In the range , the amplitudes of both oscillations are determined by expression (26), and the chirality disappears.

CONCLUSIONS
In this work, the theory is developed of the stability of finite dust particle clusters in the near-electrode plasma layer, which takes into account the nonreciprocal nature of the forces. The clusters in the form of a