Buckling of two-dimensional plasma crystals with nonreciprocal interactions

Laboratory realizations of two-dimensional (2D) plasma crystals typically involve monodisperse microparticles conﬁned into horizontal monolayers in radio-frequency (rf) plasma sheaths. This gives rise to the so-called plasma wakes beneath the microparticles. The presence of wakes renders the interactions in such systems nonreciprocal, a fact that can lead to a quite different behavior from the one expected for their reciprocal counterparts. Here we examine the buckling of a hexagonal 2D plasma crystal, occurring as the conﬁnement strength is decreased, taking explicitly into account the nonreciprocity of the system via a well-established point-wake model. We observe that for a ﬁnite wake charge, the monolayer hexagonal crystal undergoes a transition ﬁrst to a bilayer hexagonal structure, unrealizable in harmonically conﬁned reciprocal Yukawa systems, and subsequently to a bilayer square structure. Our theoretical results are conﬁrmed by molecular dynamics simulations for experimentally relevant parameters, indicating the potential of their observation in state-of-the-art experiments with 2D complex plasmas.


I. INTRODUCTION
Monodisperse negatively charged microsized particles, levitating in a plasma sheath above a powered radio-frequency (rf) electrode, form, under sufficiently strong confinement, a monolayer hexagonal complex plasma crystal [1,2].Such crystals offer the possibility to study complex kinetic phenomena in solids [3][4][5], such as crystal melting [2,[6][7][8][9][10][11] and the dynamics of dislocations [12], on a particle-resolved level.Not only have these studies shed light on the generic behavior of typical solids, but also they have revealed alternative mechanisms such as crystal melting induced by a mode coupling instability (MCI) [9,13,14].
The MCI makes its appearance in monolayer hexagonal complex plasma crystals due to the existence of wakemediated interactions between the dust particles [15][16][17].The electric field, essential for the levitation of the particles in the rf discharge, perturbs the ionic cloud around them, rendering it highly asymmetric.These asymmetric ion clouds, known as "plasma wakes" [18][19][20][21], exert attractive interactions on the negatively charged dust particles, making their pair interactions nonreciprocal and consequently prohibiting a description of the system in terms of a Hamiltonian [22,23].
Published by the American Physical Society under the terms of the Creative Commons Attribution 4.0 International license.Further distribution of this work must maintain attribution to the author(s) and the published article's title, journal citation, and DOI.Open access publication funded by the Max Planck Society.
Nevertheless, the MCI is inhibited when the gas damping is strong enough [15].Then the observation of further unconventional effects introduced by the plasma wakes, such as distinctive structural changes, can be facilitated.
The non-Hamiltonian nature of the 2D complex plasma becomes more evident for sufficiently weak confinements, where the wake-charge interactions become more prominent.In such a case, the charged particles create vertical pairs [24,25].When the interactions' nonreciprocity is large enough, the particles constituting these pairs can jointly selfpropel as a doublet, providing a remarkable example of emerging activity in complex plasma systems [26].
Regarding the complex plasma crystal, the vertical pairs of charged particles trigger, for a weak confinement, the formation of vertically aligned hexagonal layers [27,28].In individual cases of neither very weak nor very strong confinement, the formation of multilayer structures has also been observed [4,29], but so far a systematic experimental study of the structural transitions in the system is missing.
Meanwhile, several theoretical and numerical studies have systematically examined the instability of the hexagonal monolayer Yukawa crystal in the context of complex plasma [30][31][32][33][34][35][36], classical Wigner crystals [37][38][39][40][41] and charged colloids [42][43][44][45][46][47][48][49][50][51][52].It has been shown that under harmonic confinement, for a decreasing confinement strength or increasing density, the hexagonal monolayer is expected to buckle first to a hexagonal triple layer crystal and subsequently to a bilayer square [4].For even weaker confinement or higher densities, multilayer structures are expected to prevail, found occasionally also in complex plasma experiments [4,29].Although the stability of the bilayer square crystal is prominent in all the aforementioned theoretical and numerical studies, it still remains elusive in 2D complex plasma experiments.A possible reason could be that these theoretical predictions are questionable in the context of 2D plasma crystals, since they do not take into account the inherent non-Hamiltonian nature of complex plasma caused by the plasma wakes [53].
In this paper, we systematically study the buckling of 2D plasma crystals following the structural instability of the hexagonal monolayer crystal.In our theoretical and numerical investigation, we explicitly take into account the nonreciprocal character of pair interactions in the system, employing a simplified but successful model of wakes as pointlike positive charges located below the dust particles [15,28,54,55].Within this model, we show that for decreasing confinement strength, the hexagonal monolayer gives its place to a hexagonal triple or bilayer crystalline structure, depending on the value of the effective wake charge.Moreover, we identify a large stability region of the bilayer square structure, along the lines of the theoretical predictions for reciprocal Yukawa interactions [4,30,38,45].A part of this region overlaps with the stability region of the bilayer (or triple) layer hexagonal crystal, marking a regime of bistability which proves to be a source of hysteresis.Our detailed phase diagram for experimentally relevant values of the parameters not only extends the existing theoretical results to the case of nonreciprocal interactions but also allows for estimating the conditions under which each of the investigated structures can be found, facilitating their observation in experiments of 2D complex plasma crystals.

A. The model
Ignoring the thermal agitation, the equations of motion (EOM) for the dust particles in a 2D monodisperse complex plasma read where r i is the position of the particle i, m is the particle mass, and ν is the damping rate originating from the gas friction.Each particle i is subjected to two different kinds of forces, namely the interaction force F int (r i − r j ), exerted by any other particle j, and the force of the external confinement F ext (r i ).
The external force results from the assumed parabolic confinement of the particles in the vertical direction z and therefore reads where it is implied that r = xn x + yn y + zn z and con stands for the eigenfrequency of the confining potential well.The interaction force is much more involved since, except from the direct reciprocal interaction between the dust particles i and j, we should take also into account the effect of nonreciprocal interactions, stemming from the asymmetry of their surrounding ionic clouds.For the description of the latter nonreciprocal interactions, we employ a simple model, used so far successfully in literature [15,28,54,55] in various situations.Within this model, the focusing of ions downstream of the dust particles with a negative charge Q, leads to the formation of plasma wakes which can be approximated by pointlike positive charges q, located directly below each particle at a fixed distance δ [Fig.1(a)].
In this picture, the total interaction force F int (R i j ) exerted on particle i by particle j can be written as a sum of the direct reciprocal interaction between the particles i and j, F i j (F i j = −F ji ), and the nonreciprocal interaction between the particle i and the wake of particle j, F q i j (F q i j = −F q ji ), i.e., (in Gaussian units).
Here we have assumed that both the reciprocal and the nonreciprocal forces are of the Debye-Hückel (Yukawa) form with λ denoting the effective screening length, and we have introduced the notation R i j = r i − r j and R q i j = r i − r j + δn z .

B. The emergence of the structural instability
The complex plasma system described by Eq. ( 1) in view of Eqs. ( 2)-( 4) is known to relax at a hexagonal monolayer crystal [2,14] for strong enough parabolic confinement.The reference frame we use for the description of this crystal is shown in Fig. 1(b), where the equilibrium lattice constant is denoted by and the angle between the wave vector k and the x axis by θ .The reciprocal lattice is shown in Fig. 1(c), in which also the first Brillouin zone is indicated.Under these assumptions the particles' positions in the lattice are given by r (0) j = s (H )

j
, where m j , n j are specific integers, and n x , n y are the unit vectors in the real space.Hereafter, we will use dimensionless units, measuring distance in units of (i.e., = 1) and frequency in units of the dust lattice frequency We also introduce the dimensionless quantities of the effective wake charge q = |q/Q|, the effective wake distance δ = δ/ , the screening parameter κ = /λ, and the normalized confinement frequency ˜ con = con / DL .The stability of the hexagonal monolayer crystal can be investigated on the linear level by assuming a plane-wave ansatz for each particle's displacement d j from its equilibrium position r (0) j as follows: The linearization of the EOM (1) in terms of d j and the use of the above plane-wave ansatz lead, as shown in Refs.[9,14,15], to the dynamical matrix D with eigenvalues 2 j = ω j (ω j + iν), where ω j denote the system's eigenfrequencies.Under the assumption ν ω j , the hexagonal equilibrium is stable, if all 2 j are positive.As discussed in Refs.[9,14,15], two of the 2 j , in addition to their positive real part, attain a finite imaginary part as the confining frequency ˜ con decreases below approximately 3.5.This is an imprint of the mode coupling instability (MCI) of the hexagonal monolayer, for weaker confinement strengths, causing the crystal to melt.Importantly, both this instability and the induced melting can be suppressed by increasing the damping rate ν, i.e., increasing the gas pressure [15].
For an even weaker confinement, in addition to the MCI, a structural instability of the monolayer is expected to set in, causing the formation of a multilayered structure.In order to gain a deeper understanding into the nature of this structural instability, we examine the behavior of the squared frequency of the out-of-plane mode 2 v in the reciprocal space.When 2 v < 0, the out-of-plane eigenvector grows exponentially in time and the system becomes structurally unstable in the vertical direction.Thus, at the critical value of the confinement frequency for the structural instability, ˜ (SI )  con , the minimum value of 2 v (for k = 0) becomes zero.For lower confinement frequencies, con < (SI ) con , 2 v becomes negative at a certain k.Such a case, with con slightly smaller than (SI )  con , is depicted in Fig. 1(d), where we observe that 2 v attains its lowest (negative) value at θ = 30 • and a wave-vector magnitude |k (cr) | = 4π/3.This points to the fact that the structural instability sets in at k 3 e y , where e x and e y are the unit vectors in the reciprocal space.
The value of k (cr) when complemented with the information of the particular crystalline configuration [Eq.( 5)] leads to the identification of the softening mode as The crystalline structure resulting from this mode is the one described by s z, j = A cos [ 2π 3 (2m j + n j )] [Fig.1(e)], which can be identified with the hexagonal bilayer (21) with a doubly occupied bottom layer (see Sec. III A).This will be one of our first candidates for the multilayer structures resulting from the buckling of the hexagonal monolayer.

III. MONOLAYER BUCKLING: THEORY
Having seen that with decreasing confinement frequency ˜ con the monolayer hexagonal crystal undergoes not only a MCI but also a structural instability at a value ˜ (SI )  con , we examine in this section the structure of the 2D complex plasma crystal as ˜ con lowers below ˜ (SI )  con .We remark here that due to the non-Hamiltonian nature of the system, caused by the presence of wakes, the energy-minimization techniques which are commonly used to tackle theoretically such problems [30,[42][43][44][45], cannot be applied in our case.Therefore, we take the cumbersome route of first solving the force equilibrium equations to identify some of the system's crystalline equilibria and subsequently using the dynamical matrix to determine the regimes of their stability.

A. Candidate structures
As discussed above, the first unstable mode of the hexagonal monolayer yields a bilayer hexagonal structure [Fig.1(e)].Having in mind also the results for reciprocal Yukawa systems [30,43,45] which predict a bifurcation of the hexagonal monolayer (1 ) to a symmetric triple layer hexagonal structure (3 ), we proceed to the investigation of a more general structure of hexagonal order, i.e., an asymmetric hexagonal triple layer, termed hereafter as (111) and presented schematically in Figs.2(a) and 2(b).In this configuration and within the frame of reference of Fig. 1(b), the position of the jth particle in the lattice is given by the vector with s (H ) j defined above in Eq. ( 5) and FIG.
2. Sketch of the different candidate equilibrium structures of hexagonal order, expected to emerge at the buckling point and investigated in this work: [(a), (b)] a general hexagonal triple layer structure (111) and [(c), (d)] a hexagonal bilayer structure with a doubly occupied bottom layer (21).The left part shows a top view of each structure, while the right part is a side view.For the triple layer structure, the separation between the middle layer and the top or the bottom layer is denoted by h 1 or h 2 respectively.For the bilayer structure, the separation between the two layers is denoted simply by h.Note that all the particles are identical and the color distinguishes only between particles in the different layers.Also, in all the subfigures, the numbers 0,1,2 stand for the corresponding values of l as discussed in the text [Eqs.(10)] and the shaded green regions mark the unit cells of the corresponding lattices.
where in order to simplify the notation we have introduced l = mod(2m j + n j , 3).(10) According to this notation, when l = 0 the jth particle belongs to the middle layer at a height h 0 = 0 [blue particles in Figs.2( .Note that due to the presence of wakes the layers are generally not equidistant, i.e., The general triple layer hexagonal configuration (111), depicted in Figs.2(a) and 2(b), can lead, for suitably chosen h 1 and h 2 , to different configurations of hexagonal order.The most important for this study are those of an enhanced symmetry.In particular, we obtain the bilayer (21) configuration with a doubly occupied bottom layer for h 1 = h, h 2 = 0 [Figs.2(c) and 2(d)], its reverse bilayer configuration (12) for h 1 = 0, h 2 = h, and the equidistant triple layer configuration (3 ) for h 1 = h 2 .The hexagonal monolayer configuration (1 ) is the trivial case h 1 = h 2 = 0. Note that all these multilayer configurations, similarly to the (111) structure, possess a unit cell consisting of three atoms [Figs.2(a) and 2(c)].
Motivated by the results regarding the buckling of the hexagonal monolayer crystal in colloids [43][44][45], we use here two additional candidate structures, namely the bilayer staggered honeycomb (2 ) and the bilayer square (2 ), depicted  10), (13)] and the shaded green regions mark the unit cells of the corresponding lattices.in Figs.3(a) and 3(b) and Figs.3(c) and 3(d), respectively.Both these structures consist of two layers, separated by a distance h and possess a unit cell consisting of two atoms.The staggered honeycomb structure can be described by Eqs. ( 5)- (10), if we ignore the particles with l = 0 and keep only those with l = 1, 2. The position of the jth particle in the square bilayer lattice is given in a similar way by r (2 )   j = m j n x + n j n y + s (2 )  z, j n z (11) with m j , n j specific integers and where we have introduced in order to distinguish between the particles belonging to the two different layers.

B. Equilibrium structures
Evidently all the candidate structures considered here (Figs. 2 and 3) are equilibrated in the x-y plane, since due to their symmetry the x and y components of the total interaction force for each particle vanish.Thus, the only free parameters to be determined are the interlayer separations h 1 , h 2 , or h.Their equilibrium values h (0) 1 , h (0) 2 , or h (0) can be found by demanding a force equilibrium in the z direction.Since we only care about the interlayer separations, our equilibrium condition is that the z component of the net force acting on each layer should be the same with the z component of the net force acting to each other layer.For the triple-layer configuration of Figs.2(a) and 2(b), the equilibrium condition thus reads where Ft p,u is the total interaction force in the z direction exerted on the layer with l = p from the layer with l = u, both from its charges and wakes, and F q p,u is the net interaction force in the z direction acting on the layer with l = p due to the wakes of the layer with l = u.The full expressions for Ft p,u and F q p,u are provided in Appendix A, but here we remark that both are basically functions of h 1 and h 2 , so that the solutions of the system (14) supply us with the equilibrium values h (0) 1 , h (0) 2 .Of course, for a bilayer structure, such as the ones in Figs. 3, the equilibrium condition ( 14) reduces to a single equation.
Using a Newton method with different initial guesses of the values of h 1 , h 2 in order to solve numerically Eq. ( 14), we obtain all the possible equilibrium values h (0) 1 , h (0) 2 for different confinement frequencies ˜ con and different values of the effective wake charge q.Our results are presented in Figs.4( Starting with the case of zero wake charge, q = 0, we observe that for high values of ˜ con the only possible equilibrium is that of the hexagonal monolayer with h (0) 1 = h (0) 2 = 0 [Fig.4(a)].Below a critical frequency ˜ (SI )  con ≈ 2.7, the situation changes, and the hexagonal monolayer bifurcates to three qualitatively different solutions of hexagonal order.Namely, we identify with the help of Fig. 4 1 > 0, h (0) 2 = 0 (brown bold dashed line), the bilayer (12) configuration with h (0) 1 = 0, h (0) 2 > 0 (orange bold dashed line), and the symmetric triple layer (3 ) configuration with h (0) 1 = h (0) 2 (light blue bold line).Importantly, all the three different configurations emerge at the same critical value ˜ (SI )  con where the monolayer becomes structurally unstable [Fig.4(a)] and the solution space displays a high degree of symmetry, since the ( 21) and ( 12) bilayer solutions are mirror symmetric [Fig.4(a)] and the triple-layer solution is equidistant [Fig.4(b)].
This symmetry breaks for a finite wake charge as shown in Figs.4(c) and 4(d), where q = 0.2.The bilayer solutions, although they still arise at the same value of the confining frequency ˜ (SI )  con , cease to be mirror symmetric, since the interlayer separation h (0) of ( 21) is larger than that of (12) [Fig.4(c)].This can be seen a consequence of the nonreciprocity of the interactions, considering that in the presence of wakes a double occupation of the top layer (12) exerts to the bottom layer a larger attraction than the one exerted by a doubly occupied bottom layer (21) to the top layer.This wake-induced symmetry breaking affects also the triple-layer solution.For a finite wake charge q = 0.2, this becomes asymmetric (111) with h (0) 1 > h (0) 2 [Fig.4(d)] and emerges from a bifurcation of the bilayer (21) solution at a lower value of the confinement frequency ˜ con < ˜ (SI )  con [Fig.4(c)]. 1 and (f) h (0) 2 (depicted by color) as a function of ˜ con and q.The dashed-dotted white line marks the critical confinement frequency (SI )  con for the structural instability of the monolayer hexagonal structure, while the dashed red line marks the critical confinement frequency for the bifurcation of the (21) to the (111) structure, along the lines of panels (c) and (d).In all the cases, we have used κ = 1 and δ = 0.3.

The values of h (0)
1 and h (0) 2 in the branch of the triple layer (111) solution are presented as a function of the confinement frequency ˜ con and the effective wake charge q in Figs.4(e) and 4(f).We can see that the monolayer configuration bifurcates first to the bilayer configuration (21) with h (0) 2 = 0 and subsequently to a triple layer (111) with h (0) 1 > h (0) 2 for all values of q except for q = 0, where it bifurcates directly to the (3 ) structure with h (0) 1 = h (0) 2 .Also the critical confinement frequencies for both bifurcations shift to lower values for increasing q.
Regarding the other two candidate structures (Fig. 3), i.e., the bilayer honeycomb (2 ) and the bilayer square structure (2 ), they only lead to a single equilibrium solution, since they possess only one free parameter, the interlayer distance h.In contrast to the multilayer hexagonal structures discussed above, these solutions are structurally incompatible with the hexagonal monolayer solution (1 ) realized in strongly confined 2D plasma crystals.Therefore, a transition to the bilayer honeycomb (2 ) or the bilayer square (2 ) configuration could only occur in a discontinuous manner, as discussed in Sec.IV.We note that the dependence of both the (2 ) and the (2 ) structures on the effective wake charge q and the reduced confinement frequency ˜ con is very similar to that of the hexagonal bilayer structure (21) [Fig.4(e)], with the interlayer separation h increasing with decreasing ˜ con or q.

C. Stability and phase diagram
Having explored the character of the different equilibrium configurations of Figs. 2 and 3, we next investigate their stability, which is essential for their physical realization.For this reason, we linearize our equations of motion (1) around each of the equilibria {h (0) 1 , h (0) 2 } or h (0) and assume a planewave ansatz similarly to Eq. ( 7) for the displacement d j of the jth particle.Following this procedure, we construct the corresponding dynamical matrix In this expression, the D uv , with u, v = x, y, z are square submatrices of dimension equal to the number of atoms in the unit cell of the examined structure.As an example, for the case of the hexagonal triple layer structure (111) of Figs.2(a for L, l = 0, 1, 2, δ i j being the Kronecker δ, and the formulas for a uv,L , b uv,Ll presented in Appendix B. Note that in Eq. ( 16) the 0,1,2 stand for the different values of l [Eq.(10)].A similar expression can be derived also for the other equilibrium structures explored here, e.g., the square bilayer (2 ), with a replacement of l with w [Eq.( 13)].The resulting dynamical matrix D has always a dimension d equal to three times the number of particles in the crystal's unit cell.
In order to investigate the stability of each configuration, we examine the eigenvalues 2 j of D for the respective equi-libria {h (0) 1 , h (0) 2 } or h (0) .In a strict sense, stability is provided if and only if the imaginary part of all eigenfrequencies j is equal to zero, i.e., Im j = 0, ∀ j.Therefore, a suitable measure of instability can be provided by the quantity S int = d j=1 |Im j |, with it being zero corresponding to the respective configuration being stable and otherwise unstable.
As we have discussed, however, our system can display apart from structural instability (SI) also mode coupling instability (MCI) connected to the appearance of an imaginary part to some of the eigenvalues of D, 2 j [14,15].This would yield as well Im j = 0 for some j, and therefore S int = 0.In order to capture only the structural instability, we need to filter out the cases with Re 2 j > 0 and Im 2 j = 0 from S int , i.e., use the measure In the following, we judge stability by the values of Sint .If this is zero for all wave vectors k, the corresponding configuration is regarded as structurally stable and otherwise as structurally unstable.
It turns out that from the structures discussed above only the (12) hexagonal bilayer, with a doubly occupied top layer is always unstable for the 2D complex plasma system.All the other structures explored here (Figs. 2 and 3) turn to be stable in different parameter regimes.
The stability regions of each of these structures in terms of the confinement frequency ˜ con and the effective wake charge q are shown in Fig. 5, which can be interpreted as the zero-temperature structural phase diagram of our system.Here we observe that for a finite wake charge q with decreasing confinement frequency ˜ con the complex plasma crystal is expected to pass sequentially through the following structures: The reverse sequence is expected for an increasing frequency ˜ con with a hysteresis in the transition (2 ) → (21), stemming from the evident bistability region between the ( 21) and ( 2) structures.We note here that for typical finite values of q the (21) structure becomes unstable before bifurcating to the hexagonal (111) structure [Figs.4(c)-4(f)].Thus, the asymmetric triple layer structure does not appear before the transition to the bilayer square (2 ) apart from very low values of q (Fig. 5).
Overall, the buckling behavior of the hexagonal monolayer is found to be quite similar to that of the recipocal system ( q = 0), where a decreasing frequency triggers the transitions as known also from Refs.[43][44][45].The major difference is that the hexagonal monolayer (1 ) gives its place to the bilayer (21) hexagonal structure for q > 0 instead of the symmetric triple-layer hexagonal structure (3 ) for q = 0.However, when these structures become unstable, for lower ˜ con values, a bilayer square (2 ) is realized for any value of q, suggesting that the nonreciprocity of complex-plasma crystals does FIG. 5.The structural phase diagram of a quasi-2D plasma crystal with nonreciprocal interactions, according to our theoretical findings.The stability regions of the different phases, discussed in the text (Figs. 2 and 3), are depicted here by different colors.In this sense, the overlap of two different colors (e.g., orange region) indicates a regime of bistability.The control parameters are the confinement frequency ˜ con and the effective wake charge q.The other parameters are kept constant at values ρ = 2/ √ 3, κ = 1, and δ = 0.3.The horizontal dashed line marks the value of q = 0.2 used in our simulations (Fig. 6).In the white region, none of the examined phases is stable.Note, however, that overall the existence of further stable phases is not precluded.not prevent the hexagonal-to-square transition known to take place in reciprocal Yukawa systems [4,30,38,45].
Nevertheless, we note here that neither the bilayer square (2 ) nor the bilayer honeycomb (2 ) structure are stable for high values of the wake charge q > 0.5 (Fig. 5).In this region, we expect that the wake-particle attraction will be very large, causing the formation of vertical pairs which at the end will destabilize the 2D complex plasma crystal [24,25].

IV. MONOLAYER BUCKLING: SIMULATIONS
In order to check the validity of our theoretical results discussed above, we have performed molecular dynamics (MD) simulations for our point-wake model [Fig.1(a)].During our simulation time, we have been changing the confinement strength con in small steps (for a more detailed description, see Appendix C).We performed two different simulations: one starting from the hexagonal monolayer (1 ) configuration and slowly decreasing the confinement strength and another starting from a square bilayer (2 ) configuration and slowly increasing con .
Our results for the two simulations are presented in Figs.6(a) and 6(b) and Figs.6(c) and 6(d), respectively.Here we focus on the transition from the hexagonal (1 ) to the bilayer square (2 ) structure (decreasing ˜ con ) and backwards (increasing ˜ con ).As a measure of the degree of hexagonal and square order in the system, we use the quantities 6 and 8 respectively, where the global bond angular order FIG. 6. [(a), (b)] MD simulation results for the buckling of the hexagonal monolayer (1 ) for experimentally relevant parameters (see Appendix C) as the confinement strength decreases slowly in time.In panel (a), the evolution of 6 and 8 [Eq.( 21), used as measures for the hexagonal and the square order respectively] is depicted as a function of ˜ con .Panel (b) shows the evolution of the average interlayer separation h in units of (light blue markers with error bars) with decreasing ˜ con (from right to left).For comparison, we also show the theoretically expected values for the interlayer separation of the bilayer hexagonal (21) structure h (21) (yellow solid line) and the bilayer square (2 ) structure h (2 ) (purple solid line).The square insets are enlarged snapshots of the system at the indicated by the arrows points.The color of the particles encodes their vertical positions z in units of .[(c), (d)] Same as panels (a) and (b) but for the reverse process in which we start by a square bilayer (2 ) and increase slowly ˜ con (left to right).The simulation results presented here are for an effective wake charge q = 0.2, a screening parameter κ = 1, and a coupling parameter = Q 2 /(k B T ) = 1.187 × 10 4 .The vertical lines represent the theoretically predicted (Fig. 5) transition points for (1 ) → ( 21) and ( 21  (21) In the above expression, N is the total number of particles, c j denotes the coordination number of the particle j, θ j,k is the angle of the bond between the adjacent particles j and k with respect to a fixed reference direction, and k(nn) denotes the summation over all particles k which are nearest neighbors of the particle j.Within this definition, the 6 ( 8) is expected to be one for a perfect hexagonal (square) lattice and zero if no hexagonal (square) order is present.
As we observe in Fig. 6(b), for a confinement frequency ˜ con decreasing slowly with time, the initial hexagonal monolayer (1 ) bifurcates at a certain value of ˜ con to a clear bilayer hexagonal structure (21) with a doubly occupied bottom layer.At even lower values of ˜ con , the hexagonal structure becomes unstable and several domains of a bilayer square structure (2 ) appear which fill eventually our simulation box.This transition from the bilayer hexagonal to the bilayer square configuration is accompanied by an abrupt increase of the interlayer separation h [Fig.6(b)], a steep increase of 8 , and a subsequent steep decrease of 6 [Fig.6(a)].This behavior is in line with our theoretical predictions, discussed above [Fig.5, sequence (19)].Even more, the numerical values for the interlayer separation h and the two transition points are in an excellent agreement with the ones predicted by our theory.
The results for the reverse scenario [Figs.6(c) and 6(d)], in which we start from a bilayer square configuration (2 ) and increase the confinement frequency ˜ con , reveal a clear hysteretic behavior, signified by the decrease of the square order ( 8 ) and increase of the hexagonal order ( 6 ) at a larger ˜ con value than the one found for the (21) → (2 ) transition [Figs.6(a) and 6(b)].The transition point for the loss of square order proves to be well estimated by our theoretical results, in view of which the cause of the hysteresis is the bistability between the (2 ) and ( 21) structure in the regarded region (Fig. 5).The transition (2 ) → (21) appears to be overall slower than the (21) → (2 ), especially in the vertical direction, where the separation of the particles remains close to the expected one for the bilayer square (h (2 ) ) for some interval past the transition [Fig.6(d)].In addition, the structures realized in our simulations after the destabilization of the (2 ) and before the stabilization of the hexagonal monolayer (1 ), i.e., between the two transition lines, are quite disordered, rendering the identification of the different layers difficult and resulting in the large deviations in the interlayer separation h [Fig.6(d)].
Upon a decrease of its confinement frequency ˜ con below 2, the system exhibits an even more complex behavior.Its detailed description, however, goes beyond the scope of the present paper.Here we only mention that after the destabilization of the square bilayer (2 ) at a certain ˜ con , we find in our simulations the occurrence of a bilayer staggered honeycomb (2 ) structure (Fig. 7), along the lines of our theoretically calculated phase diagram (Fig. 5).Beyond the destabilization of the (2 ) structure, the system presumably enters a transient regime and at even lower confinement frequencies we expect FIG. 7. Snapshot of the particles' positions, resulting from MD simulations with slowly decreasing confinement strength, for the case ˜ con = 1.5.The color of the particles encodes the values of their vertical positions z in units of .The parameters used read q = 0.2, κ = 1 and = 1.187 × 10 4 (for more information, see Appendix C).
the successive formation of more layers, similar to the results of Ref. [29].

V. CONCLUSIONS AND OUTLOOK
In this paper, we have investigated the buckling of 2D monodisperse complex plasma crystals as their vertical confinement weakens.Unlike the Yukawa systems whose buckling has been explored in literature to a large extent [4,30,38,45], the 2D complex plasma crystals feature nonreciprocal interactions, due to the presence of the plasma wakes.This fact prohibits their description in terms of a Hamiltonian [9] and renders the standard minimization techniques unsuitable for the investigation of the structural transitions in the system.
Employing a simple point-wake model, we have solved the force equilibrium equations for different lattice structures and determined their stability through the corresponding dynamical matrices, in order to construct the structural phase diagram of the system.Here we have focused on the exploration of the regime close to the instability of the hexagonal monolayer (1 ).
For a finite wake charge, we find that below a critical confinement frequency, the hexagonal monolayer (1 ) gives its place to a bilayer hexagonal structure (21) with a doubly occupied bottom layer.This is different from what has been found for reciprocal Yukawa systems [45], where instead of the (21) a symmetric triple layer hexagonal (3 ) structure is realized and can be viewed as an imprint of the nonreciprocity, which breaks the system's symmetry.Decreasing further the confinement frequency, we observe a discontinuous transition to a bilayer square (2 ) and subsequently to a bilayer staggered honeycomb (2 ) structure, similar to the reciprocal case.Our theoretical results are confirmed by molecular dynamics simulations which also show a clear hysteresis of the (21) → (2 ) transition, owing to the existence of a bistability region of the ( 21) and (2 ) structures.
The results presented here for experimentally relevant parameters confirm our intuition that a bilayer square structure should be realized for 2D monodisperse complex plasma crystals at weak enough confinement frequencies.We believe therefore that they can be a useful guide for future experiments aiming at observing the long-sought bilayer square structures at monodisperse complex plasma systems.The major challenge thereof will be to avoid the vertical pairing which destabilizes entirely the plasma crystal for high wake charges.Thus, it would be important to determine the experimental parameter values of the wakes, e.g., by employing existing self-consistent kinetic models for the wake [21].
Future studies could also explore the buckling of a quasione-dimensional chain of microparticles.With reciprocal interactions alone, this has lead to interesting zigzag structures with nontrivial associated dynamics [56][57][58][59].With nonreciprocal interactions added, an even more intricate behavior is expected.

APPENDIX A: EXPRESSIONS FOR THE INTERLAYER FORCES
Here we provide the full expression for the forces appearing in the equilibrium condition (14).As a reference example, we use the case of a triple layer hexagonal (111) configuration with interlayer separations h 1 and h 2 , depicted in Figs.2(a) and 2(b).The equilibrium conditions for the other configurations are obtained following a very similar procedure.
As already mentioned in the text, the particles for the (111) configuration can be distinguished into three classes, indicative of their position in the unit cell (Fig. 2), according to the value of l [Eq.( 10 Using these, along with the formula for the interaction forces in our system [Eq.( 3)], we can write the expression for the z component of the total force Ft p,u exerted on the layer l = p from the layer with l = u as follows: where we have also used the DL unit [Eq.( 6)].The first sum in this expression corresponds to the interaction between the charges of layer p with the charges of layer u.The second sum, denoted by F q p,u = − m, n mod [p − (2m + n), 3] = u q κs (δ)  p s (δ) z,p , (A4) refers to the interactions between the charges of layer p with the wakes of layer u.Note that through their dependence on s p , s z,p , s (δ) p , and s (δ) z,p , both Ft p,u and F q p,u are functions of h 1 and h 2 .
For the case of the bilayer square (2 ), the situation is very similar to the one described above.The only difference is in the expressions for s x = m, s y = n, and the use of w = 0 or 1 [Eq.(13)] in the place of l, which lead to The corresponding expressions for Ft p,u and F q p,u in this case read

FIG. 1 . 3 ,
FIG. 1.(a) Schematic illustration of the interaction forces F int (r i − r j ) between the particles i and j.(b) Sketch of the elementary hexagonal lattice cell of the monolayer configuration with the frame of reference.Here denotes the lattice constant and the angle θ is the angle between the wave-vector k and the x axis.(c) Sketch of the reciprocal lattice in k space with basis vectors b 1,2 = 2π −1 (1/ √ 3, ±1).The shaded region corresponds to the first Brillouin zone with boundaries |k| = 2π −1 / √ 3 for θ = 0 • and |k| = 4π −1 /3 for θ = 30 • .(d) Color plot of the squared eigenfrequency 2 v of the out-of-plane mode, pointing (see discussion in the text) to the onset of the structural instability for q = 0.2, ˜ con = 2.48, δ = 0.3, and κ = 1.The white spots indicate where 2 v < 0 and the structural instability sets in (near θ = 30 • and |k (cr) | = 4π/3).(e) Color plot of the z value of each particle for the real part of the first unstable eigenmode d soft as discussed in the text.The particles' positions in the x-y plane, s x, j and s y, j , are depicted, with the color indicating the position of each particle in the z direction.
a) and 2(b)], when l = 1 it belongs to the top layer at a height h 1 [yellow particles in Figs.2(a) and 2(b)], whereas when l = 2 it belongs to the bottom layer with a height −h 2 [purple particles in Figs.2(a) and 2(b)]

FIG. 3 .
FIG. 3. Sketch of two further candidate equilibrium structures expected to develop from the structures shown in Fig. 2: [(a), (b)] a bilayer staggered honeycomb structure (2 ) and [(c), (d)] a bilayer square structure (2 ).The left part shows a top view of each structure, while the right part is a side view.The separation of the two layers in each configuration is denoted by h.As in Fig. 2, all the particles here are identical and the color distinguishes only between particles in the different layers.The numbers 0,1,2 stand for the corresponding values of l [(a), (b)] or w [(c), (d)] as discussed in the text [Eqs.(10),(13)] and the shaded green regions mark the unit cells of the corresponding lattices.

FIG. 4 .
FIG. 4. [(a)-(d)] Bifurcation diagrams of the hexagonal monolayer configuration in terms of the interlayer distances h (0) 1 and h (0) 2 for two qualitatively different cases [(a), (b)] q = 0 and [(c), (d)] q = 0.2.In the first column [(a), (c)], all the possible equilibrium values of h (0) 1 (and consequently also h (0) 2 ) are depicted as a function of the confinement frequency ˜ con .In the second column [(b), (d)], all the equilibrium pairs (h (0) 1 , h (0) 2 ) obtained when varying ˜ con are shown, so that we can identify the character of the different equilibria.All the solid light blue (light grey) lines correspond to triple layer configurations [(111) or (3 )], all the dashed brown (dark grey) lines to bilayer (21) configurations, and all the dashed orange (grey) lines to bilayer (12) configurations.Note that the different lines of the same color correspond essentially to the same equilibrium structure and account for all the different permutations of particles 0,1,2 in Figs.2(a) and 2(b).The thicker lines (dashed or not) guide the eye to one of the equivalent structures.[(e), (f)] The values of the equilibrium separations (e) h (0)1 and (f) h (0) 2 (depicted by color) as a function of ˜ con and q.The dashed-dotted white line marks the critical confinement frequency(SI )  con for the structural instability of the monolayer hexagonal structure, while the dashed red line marks the critical confinement frequency for the bifurcation of the (21) to the (111) structure, along the lines of panels (c) and (d).In all the cases, we have used κ = 1 and δ = 0.3.

) and 2
(b), we have that D uv = D uv,00 D uv,01 D uv,02 D uv,10 D uv,11 D uv,12 D uv,20 D uv,21 D uv,22 (16) with D uv,Ll = a uv,L δ Ll + b uv,Ll for uv = zz, D zz,Ll = a zz,L + ˜ 2 con δ Ll + b zz,Ll FIG.6.[(a), (b)] MD simulation results for the buckling of the hexagonal monolayer(1 ) for experimentally relevant parameters (see Appendix C) as the confinement strength decreases slowly in time.In panel (a), the evolution of 6 and 8 [Eq.(21), used as measures for the hexagonal and the square order respectively] is depicted as a function of ˜ con .Panel (b) shows the evolution of the average interlayer separation h in units of (light blue markers with error bars) with decreasing ˜ con (from right to left).For comparison, we also show the theoretically expected values for the interlayer separation of the bilayer hexagonal (21) structure h(21) (yellow solid line) and the bilayer square (2 ) structure h(2 ) (purple solid line).The square insets are enlarged snapshots of the system at the indicated by the arrows points.The color of the particles encodes their vertical positions z in units of .[(c), (d)] Same as panels (a) and (b) but for the reverse process in which we start by a square bilayer (2 ) and increase slowly ˜ con (left to right).The simulation results presented here are for an effective wake charge q = 0.2, a screening parameter κ = 1, and a coupling parameter = Q 2 /(k B T ) = 1.187 × 10 4 .The vertical lines represent the theoretically predicted (Fig.5) transition points for (1 ) → (21) and (21) → (2 ) in panels (a) and (b) and for (2 ) → (21) and (21) → (1 ) in panels (c) and (d).

ACKNOWLEDGMENTS
This work was supported by the German Research Foundation (DFG) under Grants No. LO 418/23-1 and No. IV 20/3-1.H.H. and C.-R.D. acknowledge the support by the National Natural Science Foundation of China (NSFC), Grant No. 11975073.The authors would like to thank E. Oǧuz, V. Nosenko, and H. Thomas for useful discussions.