Buckling transitions and clock order of two-dimensional Coulomb crystals

Crystals of repulsively interacting ions in planar traps form hexagonal lattices, which undergo a buckling instability towards a multi-layer structure as the transverse trap frequency is reduced. Numerical and experimental results indicate that the new structure is composed of three planes, whose separation increases continuously from zero. We study the effects of thermal and quantum fluctuations by mapping this structural instability to the six-state clock model. A prominent implication of this mapping is that at finite temperature, fluctuations split the buckling instability into two thermal transitions, accompanied by the appearance of an intermediate critical phase. This phase is characterized by quasi-long-range order in the spatial tripartite pattern. It is manifested by broadened Bragg peaks at new wave vectors, whose line-shape provides a direct measurement of the temperature dependent exponent $\eta(T)$ characteristic of the power-law correlations in the critical phase. A quantum phase transition is found at the largest value of the critical transverse frequency: here the critical intermediate phase shrinks to zero. Moreover, within the ordered phase, we predict a crossover from classical to quantum behavior, signifying the emergence of an additional characteristic scale for clock order. We discuss experimental realizations with trapped ions and polarized dipolar gases, and propose that within accessible technology, such experiments can provide a direct probe of the rich phase diagram of the quantum clock model, not easily observable in condensed matter analogues. Therefore, this works highlights the potential for ionic and dipolar systems to serve as simulators for complex models in statistical mechanics and condensed matter physics.


I. INTRODUCTION
The statistical physics of two dimensional systems is fascinating, as it exhibits phases that are critical for a non-vanishing range of parameters. In these phases there is no true long range order and correlations decay as a power law. The pioneering work of Kosterlitz and Thouless 1 and of Berezinskii 2 laid the foundation for the theoretical analysis of such systems. This critical behavior requires the symmetry of the Hamiltonian with respect to rotation of the order parameter by some angle variable θ. Its physics is therefore captured by an easy-plane spin Hamiltonian known as the XY -model.
An interesting situation may be encountered if the planar rotation symmetry is broken by an additional term in the Hamiltonian, which favors specific evenly-spaced values of the angle θ i = 2πi/q, for i = 1, 2...q. The resulting model is dubbed the q-state (Z q ) clock model. At low enough temperatures, the phases get pinned at one of the favored discrete values θ i , signaling the onset of true long range order. Interestingly, for q > 4 the thermal transition between ordered and disordered phases is separated by an intermediate phase with critical correlations 3 . At zero temperature, quantum fluctuations take over. Then, the ordered and disordered phases are directly connected by a quantum critical point. [4][5][6][7] The unique features of the XY and clock models in 2D are manifested in a broad and expanding set of distinct physical systems. The most familiar realizations of the 2D XY-model are associated with superfluidity and superconductivity in thin films, where several predictions were verified experimentally [8][9][10] , and in closely analo-gous systems like melting of two dimensional crystals [11][12][13][14] and dusty plasmas 15 . More recently, features of the Kosterlitz-Thouless transitions have been revealed in cold atoms experiments with Bose-Einstein condensates in two dimensions 16,17 , including the observation of vortex proliferation in two-dimensional Josephson arrays 18 . Realizations of the clock model were also discussed in a variety of contexts, for example, supersolid order of hard core bosons in triangular lattices was predicted in Refs. [19][20][21]. This type of model provides a good description for some frustrated quantum magnet systems, where a twostep melting of the six-state clock order via an intermediate, critical phase was recently discussed 22 . We further note that the physics of the clock model and the existence of an intermediate phase has an analogue in lattice gauge theories in higher dimensions 23,24 , where the critical phase corresponds to quantum electrodynamics with massless photons.
In the present paper this rich critical behavior is studied in the framework of a concrete model of experimental relevance, namely, planar crystals of strongly interacting atoms, with a particular focus on realizations with trapped ions. Ion Coulomb crystals are an example of ordered states of matter which emerge from the competition between kinetic energy and repulsive forces in confined volumes 25,26 . Their realization in atomic physics laboratories opens the possibility to gain insight into the dynamics of diverse physical models, such as Wigner crystals 27 of electrons on the surface of liquid Helium 28 , white dwarf cores, and neutron star outer crusts 29,30 . Thanks to the advances in trapping and cooling 31  Left panel: Height pattern in the ordered phase. The hexagonal lattice is split into 3 sublattices, one of which stays in the z = 0 plane (gray circles), one is raised above it (red), and one is lowered below it (blue). There are 6 inequivalent configurations, corresponding to the 3! possible height assignments to the sublattices. Right panel: The first Brillouin zone of the undistorted hexagonal lattice is shown, together with the wave vectors K and K * = −K corresponding to the height pattern.
of the external confining potential. The different structural order results from the interplay between the repulsive interactions between the ions and the trap geometry at low thermal energies, and has been extensively analyzed in the classical regime by means of numerical simulations 26,32,33 and of impressive experimental realizations [34][35][36][37][38] .
The manifestation of statistical mechanics and fieldtheoretical models in ion crystal systems has already been discussed in the context of one-dimensional ion chains. For instance, the mechanical instability of a linear chain of ions to a zigzag structure can be classically mapped to a Landau model 39 , and the field theoretical model predicts a quantum phase transition of the Ising universality class [40][41][42] . Careful numerical analysis based on DMRGtype of programs evaluated the quantum critical exponents with high precision and showed that the quantum critical region has a very small size because of the heavy ion mass [42][43][44] .
The transverse instability of planar lattices of ions, by contrast, is largely unexplored, and offers the possibility of observing much richer physics. These two-dimensional structures are presently analyzed as platforms for quantum information processing 45,46 . Early theory work by Dubin predicted a buckling instability manifested as a continuous phase transition to a three-plane structure 32 . Similar results were obtained for a crystal of atoms interacting with repulsive dipolar interaction 47 . Features hinting to the formation of three planes were measured in trapped ions experiments, even though these could not be revealed by Bragg scattering 36,37 . The predicted structure is illustrated in Fig. 1.
Intriguing properties are to be expected when fluctu-ations are taken into account. The buckling instability in Fig. 1 can be described by a complex-valued order parameter (ψ = |ψ|e iθ ), namely, the Fourier transform of the height configuration, evaluated at the wave vector K of the tripartite order. Fluctuations in the phase θ of the order parameter are captured by the XY-model. At yet lower temperatures, the phase θ is expected to be pinned to one of six favored discrete values θ i , leading to true long-range order and spontaneous breaking of Z 6 symmetry. In this work we argue that the planar-to-buckled instability is captured by the six-state clock model 3,23,48 . As already noted above, the six-state clock model predicts the existence of two transition lines at finite temperatures. This behavior follows from a duality relation between the fully ordered (low-T ) and fully disordered (high-T ) phases 23,48 , dictating two distinct critical temperatures (T l , T h ) at opposite sides of the self-dual point. It implies the formation of an intermediate, critical phase with power-law correlations in the finite range of temperatures T l < T < T h . In our case, this physics is manifested as follows: at fixed transverse trap frequency, there exists a temperature T KT below which the planar crystal becomes unstable, as shown in the phase diagram in Fig. 2. Here, the ions' transverse displacements fluctuate between three layers, and exhibit a buckled pattern with quasi-long-range order. At a second temperature T 6 , below T KT , true long-range order is established, leading to the formation of three phase-locked layers, see Fig. 1. The values of T KT and T 6 , and hence the size of the critical intermediate phase, depend on a parameter r tuned by the trap frequency, and is shown to vanish at a single quantum critical point at T = 0. By the general theory of quantum phase transitions 49 the model is here in the universality class of the XY model in D=2+1 dimensions, up to dangerously irrelevant perturbations 4 (for the formal definition see 50). Finally, for T finite but sufficiently small, we predict a crossover within the ordered phase from classical to quantum behavior at a temperature T * 6 < T 6 . This signifies the emergence of an additional characteristic scale for clock order, ξ 6 ∼ 1/T * 6 , which is intimately related to the clock term being dangerously irrelevant 5 .
As we show in this work, the realization of the quantum clock model in cold ion systems provides unique access to features that are typically inaccessible in standard condensed matter systems (e.g. magnetic compounds). In standard condensed matter systems one usually measures global properties, such as susceptibilities and conductivities. By contrast, ionic traps and cold atomic systems allow for the measurement of local properties through direct imaging, and also offer continuous control of quantum parameters. In particular, a remarkable feature of the realization we propose in trapped ions is that the different phases separated by the buckling transitions are signaled by specific features of the Bragg diffraction patterns, measurements which are typically accessible in these systems 35 . The line-shape of the peaks, for example, provides a direct measurement of the T -dependent exponent η(T ) characteristic of the quasiordered Kosterlitz-Thouless phase. In addition, we show that by analyzing the shape of the transition lines, the quantum critical properties can be probed significantly far from the quantum critical point. We further argue that these features can be observed in polarized dipolar gases 51,52 , where the quantum critical behavior may be manifested more prominently.
The outline of the paper is as follows: in Sec. II we introduce the order parameter field and derive the fieldtheoretical model; technical details of the derivation are included in Appendix A. In Sec III, we analyze the phase diagram. In Sec. IV (supplemented by Appendix B) we derive the dependence of observable physical quantities in the ordered phase on T and the parameters of the model, particularly analyzing the signature of the classical-quantum crossover mentioned above. Our main conclusions are summarized in Sec. V.

II. DERIVATION OF A FIELD-THEORETICAL MODEL
We consider a system consisting of repulsively interacting ions confined by a potential to move in the x−y plane and its vicinity. The ions form a triangular lattice of spacing a on the plane, as shown in Fig. 1. We will assume that their only freedom is to move a small distance in the direction perpendicular to the lattice. This relies on the assumption, supported by numerical data 32,47 , that the out-of-plane zigzag instability preempts a structural phase transition in the plane, e.g., from a triangular to a square lattice. Throughout this paper we focus our discussion on ions, but these considerations can be extended to ultracold dipolar atoms or molecules, as we argue in Sec. II E. Note that throughout the paper (and unless otherwise stated) we work with units for which the Boltzmann and Planck constants are k B = = 1. All length scales are in units of the lattice constant a.

A. Two dimensional ion crystals
Ions of charge Q and mass M are confined in the x − y plane by an anisotropic potential, whose strength is sufficiently large to overcome their mutual repulsion. We assume that the ion density n on the plane is homogeneous, and denote by U c the potential confining the ions along the z direction. The potential is harmonic and reads where z j is the dimensionless displacement of ion j from the plane (in units of the lattice spacing a), C is a parameter with the dimensions of an energy: and ω z is the harmonic-trap frequency. When the thermal energy is well below the Coulomb interaction energy the ions crystallize in a triangular lattice 25,26,33 . In finite systems, this is observed at a sufficient distance from the edges 36,37 . We denote by r j = (x j , y j ) the ions' positions on the plane, given in units of a. These are the solutions which minimize the in-plane potential energy U xy Q given by where K is the energy scale of the Coulomb interaction, with 0 the vacuum permittivity.
In order to set up some notation, we consider a two dimensional (2D) hexagonal lattice with lattice constant a = 1 and lattice vectors a 1 and a 2 , given by A general lattice vector can be written as with n {1,2} integers. The first Brillouin zone is a hexagon, with two inequivalent corners at the points ±K, where The other corner points are related to ±K by reciprocal lattice vectors. At fixed planar densities, the planar lattice is unstable against fluctuations in the z directions, which tend to minimize the potential energy. This instability occurs at a critical value of the transverse potential, which depends on K and has been determined by means of a stability analysis based on minimization of the interaction energy in Ref. 32 . For the purpose of our study, we consider the interparticle potential U Q resulting from the Coulomb repulsion after subtracting U xy Q : and derive an effective field theoretical model assuming small displacements |z j | 1 from the plane at z = 0.

B. Planar instability and order parameter
As the confining potential is reduced (or the planar density is increased), the system undergoes a structural transition, in which some of the atoms leave the z = 0 plane in order to reduce their mutual repulsion. In the energetically favored configuration, the hexagonal lattice divides into a tripartite lattice. This is indicated by direct minimization of the energy in numerical calculations 32,47 and by analytical studies based on the Taylor expansion of the potential in Eq. (2.7) to second order in the height z, which we report in App. A. Particles on one sublattice rise, z > 0, on another they are submerged, z < 0, and on the third they remain level at z = 0, as shown in Fig. 1. Such a configuration is described by the height z i of the particle at the location r i , and can be succinctly written as where K is the vector defined in Eq. (2.6), and ψ = |ψ|e iθ is a complex number that acts as the order parameter. Then, since the requirement that one of the three sublattices remains level at z = 0 implies that the phase of ψ can only take one of six values at the minimal energy configuration, where n ∈ {1, . . . , 6}. This suggests that the structural transition may be described in terms of a six-state clock model 3 and is corroborated by a symmetry analysis, as follows.
The hexagonal lattice is symmetric under a number of transformations, including translations, reflections, and rotations. These transformations act on both the lattice positions r i and on the heights z i . Equivalently, using Eq. (2.8), these transformations can be thought of as acting on the order parameter ψ. For example, a reflection about the z = 0 plane, R z , changes the sign of the heights z i → −z i without affecting the positions r i . This is equivalent to flipping the sign of ψ, (2.11) On the other hand, a translation T a1 by lattice vector a 1 can be absorbed into a phase shift of ψ, since Hence, A translation by a 2 acts the same way on ψ. A third example is R x , the reflection x → −x. Since K in Eq. (2.6) is parallel tox, R x changes the sign of K · r i , that is Therefore, i.e., R x acts as complex conjugation.
One can consider other symmetries of the hexagonal lattice, including rotations in the plane, the reflection y → −y, and rotations by 180 degrees along an axis lying on the z = 0 plane. However, in all cases the action of these transformations will reduce to one of the cases studied above, namely: 20) or to combinations of these actions performed in sequence. Hence, it is sufficient to consider these three basic actions.

C. Mapping to a six-state clock model
After coarse graining the lattice, one can write down a continuum Landau free energy for the order parameter ψ, which must be invariant under all the underlying symmetries of the lattice. To sixth order in ψ and ψ * , and without derivatives, this restricts the allowed terms in the free energy to: as well as 1 2 Note that the term 1 2i ψ 6 − (ψ * ) 6 is forbidden since it is not invariant under R x , Eq. (2.20). In addition, the lowest order derivative term is |∇ψ| 2 . (2.23) Thus, the Ginzburg-Landau (GL) free-energy density f GL reads where I n is a constant whose explicit form is given in Eq. (A19). For the coefficient of the quartic term, Eq. (A23) yields With the help of (A26) we also find (2.28) Here, I n and γ are constants explicitly given in Appendix A. The sign of w determines the values of the phase θ that minimize the free energy. In particular, for w > 0 this yields Eq. (2.10). We find that v is negative, a situation that may result in a first-order phase transition. In App. A it is checked that this is not the case. This behavior can be taken into account by adding a positive 8 th order term. Nevertheless, the correct physics can still be captured by Eq. (2.24) by setting v > 0, as assumed for simplicity in the discussion of next Section. When discussing quantum effects, we will need to add dynamics. The lowest order time-derivative term allowed by symmetry is is the speed of sound for transverse (z-polarized) phonons.

D. Critical point of the mean-field model
In order to connect with previous results, we remark that our theory naturally delivers the critical value of the ratio C/K of the mean field model, at which the planar crystal becomes unstable. This is found by setting r = 0 in Eq. (2.25) and yields the mean-field critical value that can be cast as a condition connecting the transverse trap frequency and the lattice constant: This value is consistent with that found in Ref. 32 (where it is reported in terms of the planar density σ). We further point out that C MF c is an upper bound to the critical value. As we will see below, at finite temperature the planar instability is at a value C c (T ) < C MF c , and decreases monotonically with T . At T = 0, moreover, quantum fluctuations lower the value of the critical point from C MF c to C c (0) [so that r c < 0, see Eq. (2.25)] by a quantity involving the ratio of kinetic and potential energies, as shown in Ref. 53.

E. Extension of the model to dipolar gases
The mapping leading to the free-energy density in Eq. (2.24) can be performed by using, instead of ions, dipolar atoms or molecules which are polarized by an external field perpendicular to the plane 51,52,54 . In Appendix A 3 we derive the coefficients of the Landau freeenergy density, corresponding to Eq. (2.24), as a function of the dipolar interaction. It must be kept in mind, however, that the dipolar interaction is not sufficiently long-range to warrant crystalline order in the plane, and the coupling with the planar modes can become relevant close to the mean-field critical point, changing the nature of the transition 55,56 . In principle, stability of the triangular lattice in the x − y plane can be enforced by a pinning potential, e.g., by an optical lattice. For this realization, the free energy in Eq. (2.24) describes the planar instability at small transverse displacements.

III. PHASE DIAGRAM
The mapping of the model to the Landau-Ginzburg free-energy in Eq. (2.24) allows us to determine the phase diagram of the system as a function of the temperature T and of the ratio K/C, as depicted in Fig. 2. In ion Coulomb crystals, the first parameter is typically tuned by means of the laser which cools the ions, while the second ratio is controlled by changing the aspect ratio between the planar and the transverse trap frequencies.
We discuss below the phase diagram for the transition at finite temperature and for the quantum phase transition at T = 0.

A. Thermal transitions
Let us first consider the phase diagram at finite temperature. At mean-field level, assuming u, v, w > 0 and v > w, the free energy in Eq. (2.24) undergoes a second order phase transition when r changes sign. For r > 0 the system is disordered, ψ = 0, whereas for r < 0 the system orders in one of six minima of the free energy, ψ = |ψ|e iπ(2n+1)/6 .
Beyond mean-field, in two spatial dimensions thermal fluctuations in the phase of the order parameter affect the phase diagram significantly. These fluctuations are captured by writing ψ = |ψ 0 |e iθ , where |ψ 0 | is the mean-field value of the order parameter and θ is its phase. Then, Eq. (2.24) yields where ρ 0 s is the bare superfluid stiffness, and is a term that tends to pin θ in one of the six values given by Eq. (2.10). On the level of mean-field theory for the model in Eq. (2.24), for r < 0 the order parameter is different from zero and given by: (a similar expression is found for dipolar interactionssee App. A). The resulting values of the parameters in Eq. (3.1) are and (3.5) Equation (3.1) defines the six-state clock model, which has been shown to exhibit two separate phase transitions 3,23 . Let us consider r fixed, r < 0. Starting from high temperatures (but still below the crystallization temperature), the system is in a disordered phase in the transverse direction, with exponentially decaying correlations. Then, as the temperature is reduced below a critical temperature, which we denote by T KT , vortices bind and the system undergoes a Kosterlitz-Thouless (KT) transition 1,2,8 . The transition temperature is given by where the coefficient ρ s is the renormalized superfluid stiffness, that depends on r and on temperature. It is smaller than the bare superfluid stiffness ρ 0 s , but is typically of the same order of magnitude, ρ s ρ 0 s . At the transition, the clock coupling h 6 is irrelevant and an intermediate phase is reached in which the phase correlations are power-law in nature, Here, η is a continuous parameter which is monotonically decreasing with temperature, and can be expressed as It attains the universal value η = 1/4 at T KT . As the temperature is further reduced, when η reaches the value η = 1/9, i.e., at the temperature T 6 < T KT , with the clock term becomes relevant and the system undergoes an additional transition. Below T 6 , long-range Z 6 order is established, in which the phase of the order parameter is locked at one of the values in Eq. (2.10).
Let us now comment on the meaning of superfluid stiffness in the present context. While the crystal of ions is clearly not superfluid, this terminology is borrowed from superconductivity and is broadly applied to XY models. In this case, it is a helicity modulus, which measures the free-energy increment associated with "twisting" the direction of the order parameter, see Eq. (3.1). Hence, it expresses the rigidity to deformations of the ordered pattern 58 .

B. Quantum critical point and its vicinity
We next focus on the quantum phase transition at T = 0. The temperature is a relevant variable at the quantum phase transition and, therefore, the transition is in a different universality class from the thermal transitions at T > 0. Furthermore, in the region surrounding the quantum critical point, the finite temperature phase diagram can be deduced from the scaling properties of the quantum critical point 59,60 , as described below.
The quantum phase transition is described by the Lagrangian in Eq. (2.30). Using the quantum to classical correspondence 49 this can be mapped to a classical sixstate clock model, Eq. (3.1), but now living in three dimensions. The three-dimensional clock model has been found to have a direct continuous transition between the long-range ordered Z 6 phase and the disordered phase. The transition lies in the XY universality class [4][5][6][7] , with the addition of the clock term h 6 which is found to be dangerously irrelevant at the critical point 4,50 . This implies that on both sides of the T = 0 transition, the quantum problem has a critical energy scale that vanishes at the transition as where ν ≈ 0.671 is the correlation length critical exponent of the D = 3 XY transition 61 . In addition, when the transition is approached from the ordered phase, the dangerously irrelevant term h 6 introduces an additional critical energy scale, whose effects will be discussed at the end of this section.
In the region of the phase diagram where h 6 is irrelevant (i.e., in the disordered and the critical phases above T 6 in Fig. 2), ∆ is the only critical energy scale. Hence, in these regions, scaling near the quantum critical point implies that the two dimensional superfluid density behaves as 59 where Φ is a universal function of ∆/T . In particular, this implies that the contours of constant exponent η in Eq. (3.8) are also contours of constant ∆/T . Hence, where A KT and A 6 are two non-universal constants. As the stiffness ρ s is monotonic in T at fixed ∆ (it is reduced by thermal fluctuations), and since η is larger at T KT than at T 6 , it follows that A KT > A 6 . Thus, we see that both transition temperatures go to zero at the quantum critical point, with the same exponent but different coefficients, and that the intermediate phase exists arbitrarily close to the quantum critical point, although it becomes narrow. We finally note that, according to Ref. 5, within the ordered phase there is a length scale ξ 6 larger than the XY correlation length ξ. In particular, near the quantum critical point ξ 6 diverges faster, ξ 6 ∼ ξ a6 where a 6 ≈ 9/4. For distances smaller than ξ 6 , the order parameter has effective XY symmetry, and only for distances beyond ξ 6 does the true long-range Z 6 order become apparent 5 .
By the quantum to classical correspondence, this implies that in the Z 6 phase at T = 0 there is a critical scale ∆ 6 which is softer than ∆, ∆ 6 ∼ ∆ a6 .
At first glance, one may expect the scaling of T 6 near the quantum critical point to be determined by ∆ 6 instead of ∆. However, this is not the case since, when T 6 is approached from above, the relevance of the clock term is determined by the value of the superfluid stiffness alone. This implies Eq. (3.13), as outlined above. One may instead define a crossover temperature T * 6 , such that (3.14) Since a 6 > 1, the temperature T * 6 is parametrically smaller than T 6 near the quantum critical point. For finite temperatures in the range T * 6 < T < T 6 the Z 6ordering is characteristic of the classical d = 2 model, while for T < T * 6 it is dominated by the quantum XYmodel in D = d + 1 = 3.

C. Practical considerations for experimental realizations on ion traps
We now make a number of quantitative estimates of parameters that may be useful in comparisons with experiments on ion traps. First, it is useful to connect the temperatures appearing in Sec. III A with experimental values for trapped ions. The features we predict will be observed below the crystallization temperature T cryst , where k B T cryst K/Γ and Γ is the plasma parameter at which crystallization occurs, Γ ≈ 140 in two dimensions 57 . Equations (3.6) and (3.4) give k B T KT 0.011 K(I 2 −C/2K). The temperature T KT can thus be a fraction of the crystallization temperature depending on the distance of the ratio C/K from the mean field critical value. For interparticle distances of the order of 15µm the crystallization temperature is in the milliKelvin range 25 . For these parameters, sufficiently far away from the mean field critical value, T KT can be of the order of hundreds of µKelvin. Comparing Eqs. (3.6) and (3.9), the transition to long-range order is at a temperature T 6 of the same order as T KT .
These considerations suggest that the intermediate critical phase could be observed using sub-Doppler cooling techniques 62,63 . We remark that in general laser cooling can lead to different effective temperatures for the planar and transverse modes of the Coulomb crystal 37 . During the experiment, thermalization between the vibrations is typically not observed, which indicates that there is a sufficiently small coupling between axial and transverse vibrations, so that it can be neglected. This leads us to conjecture that corrections due to in-plane thermal fluctuations can be neglected also sufficiently close to the transition, since in two dimensions the Coulomb interaction is sufficiently long-ranged to make the crystal incompressible along the plane 55 . This also corroborates the validity of our model, in which we as-sume that planar and transverse vibrations are consistently decoupled.
Next, we can estimate the size of quantum fluctuations in Sec. III B by computing the shift in quantum critical point r c that they induce. For this, we refer to Eq. (2.12) in Ref. 53 which, in the notation of the present paper becomes where u is given by Eq. This corresponds to a fractional shift in the trap frequency For a trap of frequency 1 MHz, the shift is 17 Hz, and it could be measured when δω z T h 1, with T h the heating time scale of the ion trap 37,64 .
The shift in r c in Eq. (3.18) can also be compared with the width of the buckled phase. Based on the stability analysis of Ref. 32, the three-layer phase is expected to exist over a range of planar densities σ, given by σa 2 0 ∈ [1.11, 1.15], where a 0 is a length scale extracted from the trapping potential. This somewhat narrow range is limited at high densities by a first order transition into a staggered square lattice structure, which occurs at σ = 1.15/a 2 0 . In terms of the parameter r, the lower density corresponds to r = 0, where the mean-field onset of the three-layer occurs, whereas the higher density corresponds to where δσ/σ = (1.15 − 1.11)/1.11 is the relative width of the phase. Hence, the fluctuation-induced shift in Eq. (3.18) is very small relative to the full width of the phase.
Another practical issue to address is the effect of the finite system size. For harmonic confinement, the triangular crystalline structure can be observed with 1000 ions at the center of the plane 26 . Present experiments can realize planes with about 1000 ions, corresponding to a linear dimension of L ≈ 30 ions. Deep in the ordered and disordered phases, when the correlation length ξ is well below L, the finite size of the system is not important. By contrast, as the critical phase is approached within either of these phases, ξ diverges and the finite system size must be taken into account. In particular, when ξ is of order L, the system will appear to be critical. This places a practical limitation in mapping out the phase transition lines. However, by studying the dependence on system size, it is possible to locate these transitions with high accuracy using relatively modest system sizes. For example, in Ref. 65, the Kosterlitz-Thouless temperature was determined with accuracy of three significant digits by simulating systems of size 12 × 12 and smaller.
On the other hand, within the critical phase itself, the finite system size does limit the ability to accurately measure critical exponents, since L = 30 gives access to only a relatively short segment of a power law. A rough rule of thumb is that one can extract one significant digit for each decade of scaling, so that here we may expect to obtain the first digit of the critical exponents. Finally, the observation of the scaling form for T * 6 , given in Eq. (3.14), will require larger system sizes than 1000 atoms. Close to the quantum critical point, in fact, this crossover occurs at a long correlation length which is associated with the dangerously irrelevant clock term. Far from the critical point, T * 6 might be extracted from the temperature dependence of the clock order parameter, as discussed below.

IV. MEASURABLE SIGNATURES OF THE PHASE DIAGRAM
We now focus on physical quantities that are in principle accessible to experimental probes. We first summarize our findings. The mechanical instability of the twodimensional Coulomb crystal gives rise to a transition to order in the transverse direction, which occurs in two stages: for temperatures T between T 6 and T KT there is a critical phase with quasi-long-range order 3 . For T < T 6 , in the ordered phase, the system spontaneously selects one of 6 equivalent patterns, i.e. breaks a Z 6 -symmetry. This structural phase transition is therefore a realization of the ordering in a six-state clock model, manifested here by splitting of the single layer of ions into three layers, at z = 0, ±h. Hence the layer separation h serves as a measurable probe of the order parameter. In addition, the resulting enlargement of the unit cell in the ordered pattern Fig. 1 is expected to induce new Bragg peaks in a diffraction experiment. Their height and line-shape as functions of T and the system parameters can provide a further test of the predictions of the theoretical model. In this Section we also propose a possible observable for the classical-quantum crossover within the ordered phase at a characteristic temperature T * 6 . For temperatures below T * 6 , the layer separation h(T ) is expected to saturate to a constant value. This criterion will serve as an operational definition of T * 6 . The discussion of this section is centered on determin-ing the signatures of the phase diagram by direct measurement of the crystal planes and by a Bragg diffraction experiment. Bragg scattering could be performed at the asymptotics of laser cooling, when the ions have reached the stationary state and the emitted photons carry the information about the crystal structure and excitations 35,66 . In addition, it is worth mentioning that the transverse shift of the crystal planes could also be measured by mapping the crystal structure into the electronic degrees of freedom of the ions 37,67-69 . These methods can in principle provide precise information of the transverse structure and fluctuations. Moreover, thanks to the advances in high precision spectroscopy, such measurements can provide the precision required to access the quantum regime.

A. Temperature dependence of the layer separation
Assuming that direct imaging of the particles is possible, the ordering established in the dark shaded area of Fig. 2 would be manifested as splitting of the plane to three layers separated by a distance where z i is given by Eq. (2.8) with i on the + or − sublattice of the broken-symmetry pattern (Fig. 1), and ... denotes the thermal expectation value. Expressing the phase of the order parameter field ψ as θ = θ 0 + δθ, where θ 0 is one of the clock states π 6 (2n + 1), one obtains h = |ψ 0 | cos(π/6) cos(δθ) . where ψ 0 is given in Eq. (3.3). While in the intermediate critical phase the ions exhibit thermal fluctuations between the three layers, for r < r c and T < T 6 the expectation value cos(δθ) = 0 and consequently h is finite in the thermodynamic limit. It can be evaluated using an effective Gaussian theory for the phase fluctuations δθ (see Appendix B for details). This yields here m(T ) is a temperature dependent coefficient, whose explicit form is given in Appendix B, and is a non-universal amplitude which depends on the high momentum cutoff Λ ≈ π/a. In the quantum regime, i.e., for T T * 6 (below the dashed line in Fig. 2), Eq.  Note that although m 0 of this magnitude is too small to be accessed directly in experiments, its value affects the temperature dependence of h(T ) at T m 0 . Hence, it may be possible to extract an estimate for m 0 (and therefore of T * 6 ) from a fit to h(T ), using Eq. (4.9).

B. Diffraction pattern
Another measurable quantity is the diffraction pattern I k 3d at wave-vector k 3d = (k, k z ), where k is a vector in the x − y plane. The diffraction pattern probes the structure factor S k 3d , that is the Fourier transform of the density operator, via the quantity (4.12) For k z h 1, this expression can be cast in the form with the leading k z -dependent contribution  9). The zero temperature value, h0, is of order one micron, as shown in Eq. (4.4). Note that the functional form of h(T ) is sensitive to the value of m0, even at temperatures well above m0 itself. Throughout, the onset temperature for clock order is taken to be T6 = 10 mK.
(details are given in Appendix C). This contribution exhibits peaks at new wave vectors k 0 = ±K + G, where K is given by Eq. (2.6) and G is a reciprocal lattice vector, reflecting the enlarged unit cell. Defining the small deviation q = k − k 0 , we get Within the Gaussian effective theory (Eq. (C6) in Appendix C), the asymptotic form of C(r) for large r is C(r) ≈ cos(δθ) 2 (1 + δθ(r)δθ(0) ) = Const. + δC(r) , (4.16) where in the entire ordered phase (T < T 6 , r < r c ) the correlation function δC(r) is exponentially decaying over a length c/m [see Eq. (C8)]. Consequently, the line-shape of the new diffraction peaks [F (q, k z ) defined in Eq. (4.15)] will be given by δ(q) + F fluc (q) where the fluctuations contribution F fluc (q) has the form of a Lorentzian with a characteristic width ∼ m(T )/c. We finally focus on the critical phase established in the higher T regime T 6 < T < T KT (the light-shaded region in Fig. 2), where quasi-long-range order is anticipated. The hallmark of such phase is a critical behavior of the correlation functions C(r), which should manifest the power-law decay of Eq. (3.7) with the Tdependent exponent η. As T is reduced from the upper critical temperature T KT , this exponent varies continuously from η = 1/4 at T KT to η = 1/9 at T 6 . We expect that this behavior will be reflected by the appearance of smeared Bragg peaks at the new wave-vectors k 0 . Their line-shape F (q, k z ) should exhibit the power-law dependence 11 F (q, k z ) ∼ |q| η−2 . (4.17) Hence, in principle a diffraction measurement can serve as direct probe of the exponent η(T ).

V. CONCLUSIONS
In conclusion, we have argued that the mechanical instability of a planar crystal of trapped ions can be mapped to the six-state clock model, and thus manifests a two-step buckling transition. On the basis of symmetry consideration we derived the field theoretical model and determined the coefficients starting from the microscopic Hamiltonian. The phase diagram is characterized by several phases, depicted in Fig. 2 as a function of the temperature T and of the external trap frequency, here parametrized by the coefficient r. It exhibits an intermediate, critical phase between disorder and long-range order in the buckled pattern of transverse displacements, within a finite range of parameters. This intermediate phase shrinks to zero as T → 0, where a unique critical point is expected. Finally, we argued that within the ordered phase, a crossover from a classical to a quantum behavior occurs at a yet lower temperature T * 6 (r) < T 6 . This signifies the emergence of an additional characteristic scale for clock order, ξ 6 ∼ 1/T * 6 , which is intimately related to the clock term being dangerously irrelevant.
These behaviors are manifested in the functional dependence of the interlayer distance h and of the Bragg signal, accessible in experiments with trapped ions, on the trap frequency and on the temperature of the crystal. In particular, in a pancake-shaped trap of the type used in Ref. 36, a gradual change of particle density towards the center of the x − y plane can yield a "wedding cake" structure 70 , and would enable the probing of h vs. the tuning parameter r as it gradually varies with the distance from the center.
We finally note that according to our estimate of the parameters for ion crystal systems, the thermal transition lines T KT (r) and T 6 (r) are expected to be observable with present trapping and cooling techniques. Direct probing of the quantum regime (by cooling to temperatures of order T * 6 ) is more challenging. However, information on quantum effects can be drawn from measurements done at accessible, higher temperatures (T T * 6 ). For example, by fitting the curve h vs. T to the theory, the energy scale T * 6 can be extracted. In addition, by analyzing the critical temperatures T KT , T 6 as a function of r, it may be possible to verify the critical exponent for the quantum transition [see Eqs. Long-range order results in the increase of the unit cell, therefore we turn first to find the wave vectors of the long-range order. Since the interaction is the strongest for nearest neighbors we assume the nature of the longrange order is determined by the nearest neighbors, while further neighbors determine the numerical values of parameters. This assumption is tested numerically.

The order determined by nearest neighbor interactions
The potential energy for nearest-neighbor interactions can be written after truncating the Coulomb sum, and takes the form where i, j denotes the sum over nearest neighbors.
Since the long-range order is characterized by a wave vector, it is convenient to introduce the Fourier Expansion of the transverse displacement, where r i is a point on the two dimensional lattice end z i is the distance of ion at this point from the plane. In Fourier space the first term in Eq. (A1) reads where we have used that z −k = z * k , being z i real. The second term on the right-hand side of Eq. (A1) is given by with Here, δr i are the vectors that join some site of the planar lattice to the corresponding six nearest neighbors, they are combinations of a 1,2 with integer coefficients ±1, 0. Using their explicit form, Finally, in Fourier space potential (A1) takes the form For sufficiently large confinement C the ions are confined to the x, y pane. This is the high-symmetry disordered phase. According to Landau theory, long-range order is expected for a vector k = (k x , k y ) for which the sign of the corresponding term in the sum (A7) first changes sign as C is increased. This happens for the value of k that maximizes F (k). A straightforward calculation shows that the appropriate wave vectors are the corners of the Brillouin zone ±K where K is given by Eq. (2.6), namely as well as for other vectors related to them by the reciprocal lattice vectors. In particular, for K in Eq. (A8), F (K) = 9 2 , and in the vicinity of the critical point the energy related to the long-range order reads The contribution of the nearest neighbors to the gradient term in the LG free-energy density is where δk is the deviation from K, or and which thus delivers the value of the constant γ in Eq. (2.24) for nearest-neighbour interactions. In the following, the contribution of farther neighbors will be calculated.

Parameters of the GL free energy for the Coulomb potential
We now derive the parameters of the Ginzburg-Landau free energy, Eq. (2.24), taking the full Coulomb interaction of Eq. (2.7). Using the order parameter (2.8) the distance from the x, y plane takes the form where r i ≡ r(n 1,i , n 2,i ) = 1 2 (n 1,i + n 2,i ), with n 1,i , n 2,i = 0, ±1, ±2, . . ., and K is given in Eq. (2.6). Substituting Eq. (A12) into Eq. (2.1), the contribution of the harmonic potential to the GL free energy is written as where N is the number of ions and is the area of the sample (in units of a 2 ). In order to obtain the GL free energy of Eq. (2.24) we expand the Coulomb energy U Q , Eq. (2.7), in powers of z i . Being the plane an equilibrium configuration, the first order term vanishes. The second-order term reads: and in terms of the order parameter ψ of Eq, (2.8) can be recast in the form: Performing the sum in Eq. (A16) over i for fixed (r j − r i ) (that can be considered a lattice vector of the form (A13)), we find where I n = n1,n2 sin n π 3 (n 1 + n 2 ) |r n1,n2 | n+1 , with |r n1,n2 | = 1 4 (n 1 + n 2 ) 2 + 3 4 (n 1 − n 2 ) 2 .
The sum is over all integers n 1 and n 2 , excluding the point n 1 = n 2 = 0. This sum can be calculated numerically, and one finds I 2 = 6.6830, where the contribution of nearest neighbors is 9 2 . Combining the contributions (A14) and (A18) one finds where we rescaled the term by K, as in Eq. (2.24). Note that U 2,Q does not depend on θ, as expected from the general symmetry considerations. This is also the property of the coefficient for the quartic term, as we show now. The fourth-order contribution to the Coulomb energy is Using similar manipulations as for the second order, one finds with I 4 = 3.5596. After rescaling by K, this leads to the coefficient The sixth-order contribution to the Coulomb energy is Using similar manipulations again one finds This is the nonvanishing term at lowest order that depends on θ, as expected from symmetry. From Eq. (A26) we identify the coefficients and where numerical evaluation yields I 6 = 2.5577. Note that v is negative. This may indicate a first order phase transition. However, we have found that this is not the case by numerical evaluation of the energy at fixed wave vector K, as a function of the order parameter ψ. For large confining potential, above the critical confining potential, C c = 2KI 2 , the minimum energy lies at ψ = 0. At C c , the minimum is still at ψ = 0, and energy is quartic in ψ at small ψ. Below C c , the minimum energy shifts away from ψ = 0 and evolves smoothly as a function of C. Hence, there is no first-order phase transition, at least at mean field level, despite the negative value of v. This results from positive higher order terms, such as |ψ| 8 , which stabilize the second-order phase transition.
In order to calculate the gradient term of Eq. (2.24) taking the contribution of all neighbors into account, we consider the second-order term of the expansion in the Coulomb interaction for deviation of the wave number δk from K: (A29) where r n1,n2 is given by (A13). Expanding in δk, the zeroth order term yields Eq. (A18), whereas the first order term vanishes in the sum. The second-order term is (A30) which has the form δU 2,Q = M ab δk a δk b , where the coefficient matrix M ab is symmetric and real. At the K point, the system is invariant under 120 degree rotations. This implies that the eigenvalues of M ab must be degenerate, i.e., M ab must be proportional to the 2 × 2 identity matrix. Otherwise, the eigenvectors of M ab would pick special directions in the xy plane, thus breaking the rotational symmetry by 120 degrees. Hence, M xy = 0 and M xx = M yy , allowing us to write M ab = Mxx+Myy 2 δ ab . This gives, The resulting gradient term is which determines the coefficient γ in Eq. (2.24).

Parameters of the GL free energy for the dipolar potential
We now calculate the parameters of the Ginzburg Landau free energy for dipolar interactions. Let us assume dipoles of magnitude P are oriented by an external field to point in the z direction. The interaction energy for a pair of dipoles connected by a vector r is in units of the lattice spacing where whileẑ andẑ are unit vectors in the z and r directions respectively. For dipoles in the x−y plane then (ẑ·r) = 0. The energy resulting of displacement of dipoles from the x − y plane is The vector K of long-range order is determined by the nearest neighbors as in the Coulomb case and takes the values of Eqs. (2.6), (A8). In what follows we expand U D in powers of z i and use similar manipulations as in the Coulomb case. The second-order term is Writing the expression in the terms of the order parameter ψ of Eq. (2.8) performing a calculation similar to the one resulting in Eq. (A18) we find where I D n = n1,n2 sin n π 3 (n 1 + n 2 ) |r n1,n2 | n+3 , with |r n1,n2 | given by (A20). The sum is over all integers n 1 and n 2 , excluding the point n 1 = n 2 = 0. This sum can be calculated numerically, and one finds I D 2 = −4.746. Combining the contributions (A14) and (A39) one finds The fourth-order contribution to the dipolar energy is that gets the form with I D 4 = 3.410. This leads to The sixth-order contribution to the dipolar energy is that gives This is the lowest order term that depends on θ, as expected from symmetry, and defines the coefficients and where numerical evaluation yields I D 6 = −2.537. The gradient term is of a similar form as in Eq. (A34), and its explicit form reads where Appendix B: Effective theory for phase fluctuations In the Z 6 -ordered phase, the effective theory describing the phase fluctuations θ(r) can be approximated by a massive quadratic form characterized by a mass m. In this Appendix we derive this effective theory, allowing us to express the physical observables h [Eq. (4.2)] and C(r) [Eq. (4.14)] in terms of θ-correlation functions. To find concrete expressions for these quantities as functions of T , we derive a self-consistent expression for m using a variational principle. In particular, this provides its T dependence, generated by fluctuation corrections.
Our starting point is the Euclidean action, which reads where τ denotes imaginary time and Ω the spatial area.
Here, c is the speed of sound, given by Eq. (2.31), and ρ s is the renormalized stiffness, which in the ordered phase is approximately equal to the bare stiffness, Eq. (3.4).
Defining the Fourier components where q ≡ (ω n , ck) with ω n = 2πnT the Matsubara frequencies, the first (free) term in Eq. (B1) is recast as S free = T ρ s 2c 2 ωn d 2 k (2π) 2 q 2 |θ q | 2 (B3) while the second (clock) term cannot be written in a simple form. However, within the ordered phase (dark shaded region in Fig. 2) this term is relevant and generates a gap, which can be approximated by a mass term correction to Eq. (B3). We therefore introduce a variational Ansatz for this massive theory in the form where the propagator G q,q assumes the form G q,q = T Ω δθ q δθ q = δ q,−q c 2 ρ s (q 2 + m 2 ) , (B5) Here δθ denotes a deviation from the configuration (θ 0 ) minimizing the clock term, and m is a variational parameter, subsequently adjusted to minimize the free energy. Employing the exact expression we find that F ≤ F var where The minimum of F var with respect to the variational ansatz S 0 is therefore an optimal approximation of the free energy F . The condition for a minimum yields ∂F var ∂G(q) = ∂F 0 ∂G(q) + T ∂ S 0 ∂G(q) = 0 (B10) (where G(q) = G q,q ) for all q. Employing Eq. (B6) we obtain S 0 = T Ω 1 2c 2 q ρ s q 2 G(q) + h 6 cos(6θ) 0 , where q ≡ kx,ky ωn>0 , thus avoiding double counting. When substituted in Eq. (B10), one gets a selfconsistent equation for the variational Green's function G(q): G −1 (q) = ρ s q 2 c 2 + 2h 6 ∂ cos(6θ) 0 ∂G(q) . (B12) Recalling that θ = θ 0 + δθ where θ 0 is one of the minima of the clock term [i.e. cos(6θ 0 ) 0 = −Sign(h 6 )], and using the quadratic form of S 0 , the expectation value cos(βθ) 0 for arbitrary β can be expressed as cos(βδθ) = e − β 2 2 (δθ) 2 , is T -independent. The self-consistent solution of Eq. (B17) yields the T -dependent mass m(T ). In the low T regime T m (below the dashed line in Fig. 2), one obtains a constant and m ≈ m 0 . However in the classically ordered regime m T T 6 we get where we have used T 6 = (2πρ s /9) (see Sec. IIIA).
Appendix C: Evaluation of the structure form factor To obtain the diffraction intensity I k 3d from Eq. (4.12), we recall the definition of the structure factor S k 3d = i e i(k·r i +kzzi) . (C1) Inserting into Eq. (4.12), one obtains For k z h 1, the last exponent can be expanded, yielding with the leading k z -dependent contribution where we have omitted a term ∝ δ(k). Employing Eq. (2.8), we express the heights z i as where K is given by Eq. (2.6). When substituted in Eq. (C4), one obtains a summation over terms with phase-factors e ±iK·(r i +r j ) or e ±iK·(r i −r j ) . Performing the sum over the center-of-mass coordinates first, the former type yields zero; from the latter type we get F (k, k z ) of the form Eq. (4.14), which relates it to the correlation function C(r).
We now evaluate C(r) within the Gaussian theory for the fluctuating field δθ derived in Appendix B, which relates it to the correlation function δθ(r)δθ(0) : C(r) = cos(δθ) 2 e δθ(r)δθ(0) (C6) where we have used Eq. (B13) for cos(δθ) . Inside the ordered phase, the Gaussian theory for δθ is massive and the propagator in momentum space is given by Eq. (B5) with m the mass associated with the clock-ordering term, evaluated as a function of T and the parameters of the model using a variational principle in Appendix B. We thus obtain δθ(r)δθ(0) = T Ω q e ik·r G q,q (C7) From this general expression, one can obtain asymptotic expressions for δθ(r)δθ(0) and consequently C(r) for large |r| in two limits -for T above or below the crossover line: δθ(r)δθ(0) ≈ c 4πρ s |r| e −m|r|/c (T m) (C8) δθ(r)δθ(0) ≈ η(T )K 0 m|r| c (m T < T 6 ) ,