Structural Transitions and Melting of Two-Dimensional Ion Crystals in RF Traps

Abstract

We investigate the structural properties and melting behaviors of two-dimensional ion crystals in an RF trap, focusing on the effects of ion temperature and trap potential symmetry. We identify distinct crystal structures that form under varying trapping conditions and temperatures through experimental observations and theoretical analyses. As the temperature increases or the trap potential becomes more symmetric, we observe a transition from a lattice arrangement to elongated ring-like formations aligned along the trap axes. Our experimental and theoretical efforts enhance our understanding of phase transitions in low-dimensional, confined systems, offering insights into the controlled formation of quantum crystals for applications in quantum simulations and many-body physics.

1. Introduction

Trapped-ion crystals provide a well-controlled platform for studying structural transitions and melting in low-dimensional systems. These systems, where ions self-organize under the influence of long-range Coulomb interactions and external trapping potentials, have been extensively used to explore fundamental aspects of phase transitions, many-body physics, and quantum simulations [1,2]. In particular, two-dimensional (2D) ion crystals confined in radiofrequency (RF) traps offer a unique setting where the interplay between trapping anisotropy and thermal effects leads to rich structural dynamics [3,4].

At low temperatures, ions in a 2D trap typically arrange into stable lattice-like structures resembling the Wigner crystal [5]. However, as the temperature increases, thermal fluctuations drive the system through a melting transition, where ions transition from localized vibrational motion around equilibrium positions to diffusive motion [6], as illustrated in Figure 1. Understanding the precise nature of this transition is critical for applications in precision metrology [7], quantum computing [8], and studies of non-equilibrium statistical mechanics [9].

In addition to temperature, the symmetry of the trap potential plays a crucial role in determining crystal structure and crystal melting. When the trapping potential is highly anisotropic, the ions form elongated, elliptical configurations, while more symmetric traps favor compact, circular arrangements. The dependence of melting behavior on trap anisotropy remains an open question, particularly regarding how structural deformations affect the critical temperature $T_c$ and whether standard melting criteria, such as Lindemann’s criterion [10], remain valid in this confined geometry.

In this work, we experimentally observe the structural properties and melting behaviors of 2D ion crystals in an RF trap and investigate them using a combination of theoretical and numerical approaches. We analyze the transition from lattice-like configurations to ring-like structures as a function of temperature and trap anisotropy. The critical temperature for melting is determined using both the Lindemann criterion and free energy calculations, allowing us to assess their consistency and limitations. Furthermore, we compute the diffusion coefficient of ions as a function of temperature and compare its behavior with the melting threshold obtained from structural fluctuations. Our results provide new insights into phase transitions in low-dimensional confined systems and have implications for optimizing ion-crystal-based platforms for quantum simulations and many-body experiments.

2. Experimental Techniques

The experiments focus on the generation, control, and characterization of two-dimensional (2D) Coulomb crystals of trapped barium ions (${\mathrm{Ba}}^{+}$) in a radiofrequency (RF) trap. These ion crystals are of interest for quantum computing, quantum simulations, and fundamental studies of phase transitions in low-dimensional systems [11,12].

The ion trap used in these experiments is a modification of the standard Paul trap, where the ring electrode is flattened and divided into eight sectors [13]. A diagram of the trap is shown in Figure 2. This configuration allows for the independent control of the trapping potential in the plane of the crystal. Two endcap electrodes, shaped as truncated cones (not shown in the figure for clarity), ensure strong transverse confinement and provide optical access for imaging and laser addressing. The trap operates at an RF frequency of approximately 10 MHz, with the amplitude in the range of 800–1500 V applied to the endcap electrodes via a double-coil helical resonator. The resulting pseudopotential supports the formation of 2D Coulomb crystals by balancing the in-plane repulsion from ion–ion interactions and the transverse confinement from the RF field [14]. For more details about this ion trap, please see [13].

Barium ions are loaded into the trap using a two-step photoionization process with a 791 nm external cavity diode laser (ECDL) and a 337 nm nitrogen laser. Once ionized, the ions are Doppler-cooled using 493 nm and 650 nm laser beams, which address the $6S_{1/2}\to 6P_{1/2}$ and $5D_{3/2}\to 6P_{1/2}$ transitions, respectively. The fluorescence of the ions is imaged using an electron-multiplying charge-coupled device (EMCCD) camera. The magnification of the imaging system is approximately 80. The laser-cooling process leads to crystallization, where the ions form a radial 2D structure due to their mutual Coulomb repulsion and the symmetry of the trapping potential [15].

A key challenge in 2D Coulomb crystals is excess micromotion, which results from the RF field and leads to time- and position-dependent Doppler shifts. Since the micromotion amplitude increases with radial distance from the trap center, ions in different regions of the crystal experience different effective cooling rates. To address this issue, augmented Doppler-cooling schemes [12,15,16] can be used to stabilize crystals of 50 ions, or alternate trap geometries that minimize or eliminate the micromotion in the direction of Doppler cooling. In the latter configuration, the stable manipulation of large 2D crystals of over 100 ions for quantum simulation has been achieved [11,17].

The trapped ions form concentric shells as pictured in Figure 2, with each shell corresponding to a different radial equilibrium position. The number of shells and their spacing are determined by the interplay between the Coulomb repulsion and the trapping potential. Experimental images confirm the expected shell structure, with ion positions matching theoretical predictions. The shape of the crystal can be dynamically tuned by adjusting the DC biases on the sector electrodes or modifying the endcap voltages, allowing for transitions between nearly circular and elliptical crystal geometries.

We use this experimental data to extract information about the ellipticity, as well as the melting behavior of trapped-ion crystals. During melting, the ions begin moving along the elliptical paths with increasing speed, eventually forming continuous rings in experimental images, as illustrated in Figure 1. Furthermore, the stability of the crystalline structure at different temperatures allows us to determine the threshold conditions for the transition between solid and melted states.

3. Theoretical Approach

We consider a system of N identical charged particles of mass m and charge e confined in an anisotropic quadratic potential trap. The effective confinement is created by the combination of the pseudopotential due to the high-frequency RF oscillating electric field and a static electric field [1,18]. In such traps, the rapid oscillations of the RF field average out, leading to a time-averaged potential that effectively confines the ions to a plane perpendicular to one of the principal axes of the trap (taken to be the x-axis in the following discussion). The interplay between the Coulomb repulsion and the external confinement results in a self-organized ion crystal, where the equilibrium positions of the ions are determined by a balance of these competing forces.

The Hamiltonian governing the system is given by

(1)$\begin{array}{c}\hfill H=\sum _{i=1}^N\left({\displaystyle \frac{p_i^2}{2m}}+{\displaystyle \frac{m}{2}}\left({\omega }_y^2{y}_i^2+{\omega }_z^2{z}_i^2\right)+\sum _{j>i}^N{\displaystyle \frac{\alpha }{\left|{\mathit{r}}_i-{\mathit{r}}_j\right|}}\right),\end{array}$

where the summation is over all the ions whose positions in the y-z plane are given by $y_i$ and $z_i$, and whose momenta are $p_i$, respectively; ${\omega }_y$ and ${\omega }_z$ are the secular frequencies in the plane of the crystal, and $\alpha =e^2/4\pi {ϵ}_0$ is the Coulomb interaction coefficient. The secular frequencies define the degree of anisotropy $\gamma ={\omega }_z/{\omega }_y$ in the potential well, which directly affects the shape and stability of the ion crystal.

3.1. Dimensionless Representation of Positions

To simplify the mathematical description, we introduce the characteristic length scale of the system using the Bohr radius of the ion:

(2)$\begin{array}{c}\hfill a_0={\displaystyle \frac{4\pi {ϵ}_0{\hslash }^2}{m{e}^2}},\end{array}$

where ${ϵ}_0$ is the vacuum permittivity and ℏ is the reduced Planck’s constant. This defines a natural unit of length associated with the quantum mechanical Coulomb problem. Then, using $a_0$, we introduce the dimensionless coordinates and time:

(3)$\begin{array}{c}\hfill \mathit{r}\to {\left({\displaystyle \frac{{\hslash }^2}{m^2{\omega }_y^2{a}_0}}\right)}^{1/3}\mathit{r},t\to t/{\omega }_y.\end{array}$

This rescaling allows us to express all quantities in terms of dimensionless parameters that highlight the fundamental interplay between the kinetic, potential, and interaction energies.

Then, the Hamiltonian of the system, in units of $\hslash {\left({\displaystyle \frac{\hslash {\omega }_y^2}{m{a}_0^2}}\right)}^{1/3}$, can be written as follows:

(4)$\begin{array}{c}\hfill H_{\mathrm{d}}=\sum _{i=1}^N\left({\displaystyle \frac{v_i^2}{2}}+{\displaystyle \frac{1}{2}}\left(y_i^2+{\gamma }^2{z}_i^2\right)+\sum _{j>i}^N{\displaystyle \frac{1}{\left|{\mathit{r}}_i-{\mathit{r}}_j\right|}}\right),\end{array}$

where $v_i$ is the speed of the ith ion. In this form, the dimensionless Hamiltonian explicitly shows how the system’s behavior is controlled by the ratio of trapping frequencies, or the trap anisotropy $\gamma ={\omega }_z/{\omega }_y$, and the number of ions N. By varying these parameters, we can explore different structural regimes, from strongly confined Wigner crystals to fluid-like configurations at higher temperatures.

3.2. Melting Condition

For the finite systems, the concept of critical temperature is inherently ambiguous, as phase transitions are not sharply defined in such cases [19,20]. Nevertheless, structural changes occur, making temperature a key parameter in characterizing the system’s behavior. To describe this, we consider the partition function in dimensionless units of position r and momentum p:

(5)$Z\propto \int d^{2N}p{d}^{2N}r{e}^{-(T_0/T)H_{\mathrm{d}}},$

where $k_B$ is the Botzmann constant, T is the temperature, and the parameter $T_0={\displaystyle \frac{\hslash }{k_B}}{\left({\displaystyle \frac{\hslash {\omega }_y^2}{m{a}_0^2}}\right)}^{1/3}$ sets the temperature scale. Since the Hamiltonian in Equation (4) depends only on the ion positions and the anisotropy parameter $\gamma$, we expect the critical temperature to take the general form

(6)$\begin{array}{c}\hfill T_c=T_{0}g(\gamma ,N),\end{array}$

where function $g(\gamma ,N)$ encapsulates the effects of anisotropy and the finite number of ions. This scaling implies that, when expressed in units of $T_0$, the melting temperature should primarily depend on $\gamma$ and N, enabling direct comparisons between theory and simulations. Additionally, due to the system’s rotational symmetry, we have the property $g\left(\gamma ,N\right)=g\left({\gamma }^{-1},N\right)$, which is particularly useful, as it allows us to analyze the system’s behavior by focusing only on cases where $\gamma \le 1$.

To achieve a simple estimate of the critical temperature $T_c$, we use Lindemann’s criterion, which states that the amplitude of ionic vibrations must not exceed a certain fraction c of the interparticle spacing for the system to remain crystalline [10,21,22]. Assuming that the system is harmonic, the amplitudes of the ionic vibrations scale linearly with c, and the energy related to these vibrations scales as $c^2$. The average kinetic energy of ions scales as the system’s temperature—$k_{B}T$—while the average vibrational potential energy scales as $k_B{T}_0$, proportional to the square of the amplitude. For a harmonic system, the average kinetic energy is equal to the average potential energy, and, thus, we obtain for the critical temperature,

(7)$\begin{array}{c}\hfill T_c=c^2{T}_0=c^2{\displaystyle \frac{\hslash }{k_B}}{\left({\displaystyle \frac{\hslash {\omega }_y^2}{m{a}_0^2}}\right)}^{1/3}.\end{array}$

This equation provides a simple yet physically meaningful estimate. For typical ion trap frequencies (${\omega }_{y,z}\sim 2\pi \times 100$ kHz), and taking $c=1$, we obtain $T_c\sim 0.01K$. This corresponds to an approximation where the function $g(\gamma ,N)$ is treated as constant equal to $c^2$, and, therefore, it does not account for the dependence on the number of ions or the anisotropy parameter. In the following subsections, we develop a more precise approach based on Lindemann’s criterion, incorporating a more detailed analysis of structural fluctuations.

3.3. Toy-Model Theory of Melting

At low temperatures, ions in a 2D trap arrange into a series of concentric elliptical shells due to the anisotropic confinement. The equilibrium positions of ions $(y_i,z_i)$ within a given shell satisfy the constraint

(8)$\begin{array}{c}\hfill y_i^2+{\gamma }^2{z}_i^2=R_n^2,\end{array}$

where $R_n$ is the characteristic radius of the n-th shell.

A key idea of this toy model is that an ion in the outermost shell is effectively confined within a potential well created by the Coulomb interaction with the inner shell, which consists of $N_n$ ions. Experimental observations indicate that melting predominantly occurs via angular motion rather than radial motion, suggesting that the angular coordinate $\phi$ is the “soft” degree of freedom in the system. In contrast, the radial coordinate is effectively constrained by the ellipticity of the trap, meaning that an ion radius is dictated by the equation of the ellipse. As a result, the system’s dynamics can be effectively described using only the angular coordinate, reducing the problem to a one-dimensional potential landscape along $\phi$.

Furthermore, we assume that the outermost shell melts before the inner ones, which is supported by experimental observations. This implies that the angular displacements become more pronounced closer to the boundary of the ion crystal. Consequently, in our model, we consider the inner ions to be nearly immobile compared to those in the outermost shell, as their motion is significantly restricted by the stronger Coulomb forces and the trap potential. This approximation simplifies the description of the melting process while retaining the essential features observed in the experiment.

Another important aspect of this model is the role of the trapping potential. While the external confinement plays a crucial role in shaping the overall crystal structure, its effect on individual ions within the same shell is approximately uniform. This is due to the fact that, for a given shell, all ions share the same elliptical constraint. Thus, their trapping potential energy remains nearly constant, and the variations in the total energy are dominated by the Coulomb interaction with the inner shell ions.

To estimate the effective potential barrier that an outer-shell ion must overcome to delocalize, we model the interaction energy of an ion at angular position $\phi$ and radial position r due to all $N_n$ ions in the lower shell as follows:

(9)$\begin{array}{c}\hfill V\left(\phi \right)=\sum _{i=1}^{N_n}{\displaystyle \frac{1}{\sqrt{r_i^2+r^2-2r_{i}rcos\left({\theta }_i-\phi \right)}}}.\end{array}$

Here, $r_i$ and ${\theta }_i$ are the radial and angular coordinates of the i-th ion in the inner shell, which satisfy the condition

(10)$\begin{array}{c}\hfill r_i^2={\displaystyle \frac{2R_n^2}{1+{\gamma }^2+\left(1-{\gamma }^2\right)cos\left(2{\theta }_i\right)}},\end{array}$

while the radial position r of the outer-shell ion satisfies

(11)$\begin{array}{c}\hfill r^2={\displaystyle \frac{2R_{n+1}^2}{1+{\gamma }^2+\left(1-{\gamma }^2\right)cos\left(2\phi \right)}}.\end{array}$

These expressions account for the elliptical distortion of the shells due to the anisotropic trapping potential.

To further simplify this toy model while retaining the essential physics, we assume that the ions within a given shell are distributed uniformly, such that

(12)$\begin{array}{c}\hfill {\theta }_i={\displaystyle \frac{2\pi i}{N_n}},i=1,\dots ,N_n.\end{array}$

This approximation is precise only for zero anisotropy, but it remains a good approximation even at modest non-zero anisotropy as well.

With these assumptions, the calculated total potential energy of an outer-shell ion as a function of its angular coordinate $\phi$ is plotted in Figure 3 with 10 ions in the shell, for the anisotropy of $\gamma =0.9$. This function exhibits two notable features: global minima at ${\phi }_{min}={\displaystyle \frac{\pi }{N_n}}\left(2⌊{\displaystyle \frac{N_n-2}{4}}⌋+1\right)$ and global maxima at ${\phi }_{max}=0,\pi$, representing the highest energy barrier an ion must overcome to become delocalized within the shell. Here, the square brackets $\lfloor ·\rfloor$ denote the floor function, which rounds down to the nearest integer.

Since melting occurs when the thermal energy $k_{B}T$ becomes comparable to the depth of this potential well, the critical temperature can be estimated by evaluating the potential energy difference between these two points:

(13)$\begin{array}{c}\hfill {\displaystyle \frac{T_c\left(\gamma \right)}{T_0}}=V\left({\phi }_{max}\right)-V\left({\phi }_{min}\right).\end{array}$

Using this equation, we can construct the theoretical phase diagram shown in Figure 4, which illustrates how the estimated critical temperature depends on the anisotropy parameter $\gamma$. As the anisotropy parameter deviates from 1 (fully isotropic case), the energy difference between the global minima and the global maxima increases (see Figure 3), and the critical temperature increases as well.

4. Numerical Simulations

In our study, we focus on systems with 14 and 28 ions because these configurations allow for a well-defined analysis of melting, as all ions form distinct shells. For other numbers of ions, equilibrium structures can exhibit significant variations, with some ions not fitting into well-formed shells. This structural variability complicates the definition of melting since different equilibrium states may correspond to different probability distributions for melting as a function of temperature and anisotropy. Moreover, the formation of a specific ion configuration in a 2D Coulomb crystal is not deterministic in these conditions, even when trapping parameters and ion numbers remain the same. Due to the presence of multiple energetically similar configurations, it is challenging to consistently reproduce an identical crystal structure across different experimental realizations. This intrinsic uncertainty in equilibrium states contributes to variations in melting behavior and critical temperature estimations. Similar challenges in achieving reproducible crystal configurations have been observed experimentally, as discussed in Ref. [17], where structural differences arising from small perturbations in initial conditions were analyzed. Understanding this variability is crucial for improving control over trapped-ion systems and refining theoretical models of their phase transitions.

To analyze the melting transition of trapped-ion crystals, we employ molecular dynamics simulations. This approach allows us to investigate the system’s behavior at different temperatures and extract key thermodynamic properties, such as the Lindemann criterion for melting and the probability of phase transition. The melting probability is determined by the fraction of simulations in which the crystalline structure is lost at a given temperature and anisotropy parameter.

4.1. Molecular Dynamics Simulations

The molecular dynamics simulations are based on integrating the classical equations of motion:

(14)$\begin{array}{c}\hfill m_i{\ddot{\mathit{r}}}_i=-{\nabla }_{i}V,\end{array}$

where V is the total potential energy, consisting of the harmonic trapping potential and the Coulomb interaction between ions.

To numerically integrate the equations of motion, we use the Verlet algorithm [23], which ensures energy conservation and stability over long simulation times. The time step is chosen to be sufficiently small to resolve the characteristic vibrational frequencies ${\omega }_y$ and ${\omega }_z$.

Unlike traditional thermostats such as Langevin dynamics [24], we set the temperature through the initial kinetic energy of the ions:

(15)$\begin{array}{c}\hfill K=N{k}_{B}T.\end{array}$

The initial positions of the ions are set at the calculated equilibria, corresponding to the minimum of potential energy. This configuration ensures that the system starts in the crystalline state at low temperatures. This process is much faster than the Langevin dynamics, in which one has to generate random numbers on every time step, but provides identical results.

To determine whether the system undergoes melting, we apply Lindemann’s criterion for each ion. If the displacement exceeds the threshold for a given ion, that ion is considered to have “melted”.

The melting probability is then estimated as the fraction of simulations in which at least one ion meets this criterion at a given temperature and anisotropy parameter. This probability is statistically averaged over multiple simulation runs to obtain reliable phase transition data, which are plotted in Figure 5. For 14 ions, we performed 400 simulation runs for each pair of T and $\gamma$, whereas for 28 ions, we only performed 25 runs due to the increased computational complexity and limited computational power.

4.2. Diffusion Coefficient Estimation

The transition to the liquid phase is signaled by a sudden increase in diffusion coefficient D, marking the melting temperature. This approach has been successfully applied to analyze phase transitions in trapped-ion systems [1,18]. Similar methodologies have been employed in molecular dynamics studies of Coulomb crystals [25]. To quantify the melting transition, we compute the mean squared displacement (MSD) for all ions:

(16)$\begin{array}{c}\hfill \mathrm{MSD}\left(t\right)={\displaystyle \frac{1}{N}}\sum _{i=1}^N⟨|{\mathit{r}}_i\left(t\right)-{\mathit{r}}_i\left(0\right){|}^2⟩.\end{array}$

For the crystalline phase, the ion motion remains bounded, and the MSD saturates at a finite value, while in a liquid-like state, ions exhibit diffusive motion, leading to a linear increase in MSD at long times. The diffusion coefficient is extracted from this long-time behavior as follows:

(17)$\begin{array}{c}\hfill D=\underset{t\to \infty }{lim}{\displaystyle \frac{\mathrm{MSD}\left(t\right)}{4t}}.\end{array}$

By evaluating D as a function of temperature T and anisotropy parameter $\gamma$, we obtain the diffusion behavior of the system, as shown in Figure 6. While the diffusion coefficient provides an independent measure of melting, it is important to note that extracting $T_c$ from $D\left(T\right)$ requires computing numerical derivatives, which introduces additional sources of errors due to finite-difference approximations.

Additionally, the diffusion-based approach is computationally more expensive compared to the direct probability-based method. Calculating the diffusion coefficient requires long-time averaging for accurate estimation; thus, the resulting phase boundary has a larger error. At the current level of statistical uncertainty, with a maximum error of 0.00334, the melting temperatures derived from the diffusion analysis and those obtained from the melting probability method may either coincide for different anisotropies or exhibit significant discrepancies.

There are regions where the diffusion coefficient decreases with increasing temperature (for example for $\gamma =1$), which may seem counterintuitive since higher temperatures generally lead to greater disorder and increased ion mobility. However, because melting in this system occurs ring by ring, ions are not free to move arbitrarily but are instead constrained to motion along their respective rings. At certain temperatures, this constraint may effectively reduce diffusion, as ions experience increased interactions with their neighbors while still being partially confined. Additionally, structural rearrangements within the rings could temporarily stabilize ion positions before full melting occurs, further suppressing diffusion. This behavior suggests that melting is not a simple monotonic process but rather involves intermediate structural transitions where diffusion can momentarily decrease.

5. Discussion

Despite the relatively small size of the simulated system, our results successfully capture the phase transition behaviors of two-dimensional ion crystals in an RF trap. By systematically analyzing the melting process, we have obtained phase diagrams for systems with 14 and 28 ions. A key observation is that the critical temperature for melting is consistently lower for the 28-ion system across all values of anisotropy, as shown in Figure 5. As seen in Equations (13) and (9), the critical temperature decreases with increasing $R_n$ and $R_{n+1}$. Since the radii of the outermost and second-to-last rings grow as the number of ions in the system increases, the overall effect leads to a reduction in the critical temperature for larger ion crystals. This aligns with the observed trend that increasing system size results in a lower melting threshold, as the spatial expansion of the outermost rings weakens the effective confinement and reduces the characteristic energy scale of structural fluctuations.

Our theoretical model also predicts an increase in the critical temperature with the increasing anisotropy of the system. However, the numerical results exhibit non-monotonic behavior in certain regions. While some discrepancies exist between the simple theoretical model and the numerical simulation results, the overall trends are in good agreement. The dependence of the critical temperature on anisotropy was fitted in Figure 7 using Equation (13). Although the toy model does not fully capture the non-monotonic intervals of this dependence, it still provides a close approximation to the overall trend. One possible explanation for this discrepancy is that the toy-model approximation is only valid for low anisotropy, where the agreement is indeed quite good.

To further validate our findings, we computed the diffusion coefficients for the 14-ion system across a range of temperatures and anisotropies. In principle, the diffusion coefficients provide an independent means of estimating the melting transition, complementing the structural analysis. However, due to significant uncertainties in the computed values, extracting precise critical temperatures from the diffusion data remains challenging. This highlights the need for the further refinement of numerical techniques or alternative methodologies to improve the accuracy of diffusion-based melting estimations.

Overall, our study provides valuable insight into the melting dynamics of confined 2D ion crystals. The agreement between theoretical and numerical approaches supports the robustness of our model, while the observed deviations emphasize the complexity of melting in finite systems.

Future work should explore larger ion ensembles and refine theoretical approximations to improve quantitative accuracy. One key objective is to refine the computation of the diffusion coefficient by reducing statistical uncertainties, thus enabling a more precise determination of the critical temperature. This will involve increasing the number of independent trajectories and optimizing averaging techniques.

Another promising approach is the calculation of the free energy, which would provide a more fundamental thermodynamic perspective on the melting transition. Implementing free energy methods, such as thermodynamic integration, will allow for a direct comparison with structural and dynamic criteria for melting.

Yet another alternative method for determining the $T_c$ involves analyzing the velocity autocorrelation function (VACF) [26]. This approach provides insight into the dynamical properties of the system, particularly the transition from localized oscillatory motion in the crystalline phase to the diffusive motion in the liquid phase. By examining how the VACF evolves with temperature and anisotropy, it will be possible to establish an independent criterion for the melting transition, thereby complementing both structural and diffusion-based analyses.

These future directions aim to enhance the accuracy and robustness of critical temperature estimation in 2D trapped-ion crystals, providing deeper insights into phase transitions in confined Coulomb systems.

6. Conclusions

We have explored the structural transitions and melting behaviors of two-dimensional ion crystals in an RF trap, combining theoretical analysis, molecular dynamics simulations, and experimental insights. Our results show that the melting transition strongly depends on ion number and trap anisotropy, with larger systems exhibiting lower critical temperatures. We identified that melting occurs through a ring-by-ring mechanism, with outer ions losing structural order first.

Despite some deviations, our theoretical model provides a reasonable approximation of the melting behavior, capturing key trends in the phase diagram. Additionally, we compared different approaches for determining the melting transition, showing that diffusion-based methods, while useful, introduce larger uncertainties than probability-based calculations.

These findings contribute to the broader understanding of phase transitions in confined Coulomb systems and offer insights for optimizing ion-crystal-based platforms for quantum simulations and many-body physics. Future work will refine the theoretical model and extend numerical analyses to improve accuracy in determining critical temperatures.

Author Contributions

Conceptualization, B.V.P.; Methodology, B.V.P. and A.K.; Software, B.V.P.; Investigation, B.V.P. and A.K.; Resources, A.K.; Writing—original draft, B.V.P.; Writing—review & editing, B.V.P., A.K. and B.B.B.; Supervision, B.B.B.; Project administration, B.B.B.; Funding acquisition, B.B.B. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement

Data sets generated during the current study are available from the corresponding author upon reasonable request.

Conflicts of Interest

Alexander Kato was employed by the IonQ, Inc. The authors declare no conflicts of interest.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures

Figure 1 An example of melting transition in a crystal or 8 trapped Ba⁺ ions confined in a 2D trap as the ion temperature increases. Left: the ions form a crystalline structure with two distinct shells at a lower temperature. Right: the 6 ions of the outer shell have melted at a higher temperature. The temperature is controlled by tuning the Doppler-cooling laser frequencies.

Figure 2 The schematic of the top view of the ion trap used in the experiments. The eight sectors of the ring electrode are shown, with two electrodes highlighted in red that were used to control the anisotropy of the trapping potential. An image of a 14-ion 2D crystal is superimposed (not to scale), with the red lines outlining the shell structure of the crystal.

Figure 3 Dimensionless potential energy $V\left(\phi \right)$ of the outer shell ion as a function of its angular position for the trap anisotropy $\gamma =0.9$. The red dashed line indicates one of the global minima.

Figure 4 Toy-model phase diagram. The critical temperature $T_c$ in units of $T_0$ is plotted against the anisotropy parameter $\gamma$. The light blue-shaded region corresponds to the crystalline state, while the pink-shaded region corresponds to the melted state. The minimum of the critical temperature is at the point of $\gamma =1$, as expected.

Figure 5 Melting probability for different temperatures and anisotropies for crystals contained 14 ions (left) and 28 ions (right). Fewer numerical simulation runs were carried out for the latter, resulting in noisier results. Nevertheless, the lower overall melting temperature for the larger crystal is evident.

Figure 6 Diffusion coefficient for different temperatures and anisotropies for the 14-ion crystal.

Figure 7 Comparison of the critical temperature $T_c$ as a function of the crystal anisotropy $\gamma$ between toy-model theory predictions and the molecular dynamics simulations for the crystal containing 14 ions. While there is general agreement, the toy model does not predict the structure observed in the simulations.

References

1. Leibfried, D.; Blatt, R.; Monroe, C.; Wineland, D. Quantum dynamics of single trapped ions. Rev. Mod. Phys.; 2003; 75, 281.

2. Häffner, H.; Roos, C.F.; Blatt, R. Quantum computing with trapped ions. Phys. Rep.; 2008; 469, 155.

3. Diedrich, F.; Peik, E.; Chen, J.M.; Quint, W.; Walther, H. Observation of a Phase Transition of Stored Laser-Cooled Ions. Phys. Rev. Lett.; 1989; 62, 403.

4. Schiffer, J.P. Phase Transitions in Coulomb and Dipole Systems. Phys. Rev. Lett.; 1993; 70, 818. [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/10054211]

5. Wigner, E. On the Interaction of Electrons in Metals. Phys. Rev.; 1934; 46, 1002.

6. Duca, L.; Mizukami, N.; Perego, E.; Inguscio, M.; Sias, C. Orientational Melting in a Mesoscopic System of Charged Particles. Phys. Rev. Lett.; 2023; 131, 083602. [DOI: https://dx.doi.org/10.1103/PhysRevLett.131.083602] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/37683173]

7. Keller, J.; Partner, H.L.; Burgermeister, T.; Mehlstäubler, T. Precise Determination of Thermal Radiation Properties of an Ion in a Paul Trap. Phys. Rev. A; 2016; 93, 013409.

8. Monroe, C.; Meekhof, D.; King, B.; Wineland, D. A “Schrödinger Cat” Superposition State of an Atom. Science; 1996; 272, 1131.

9. Schaetz, T. Trapped-ion quantum simulation and experimental quantum information. C. R. Physique; 2011; 12, 5.

10. Lindemann, F.A. Über die Berechnung molekularer Eigenfrequenzen. Phys. Z.; 1910; 11, 609.

11. Guo, W.X.; Zhang, C.; Qi, B.X.; Zhou, Z.C.; He, L.; Duan, L.M. A site-resolved two-dimensional quantum simulator with hundreds of trapped ions. Nature; 2024; 630, pp. 613-618. [DOI: https://dx.doi.org/10.1038/s41586-024-07459-0] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/38811737]

12. D’Onofrio, M.; Xie, Y.; Rasmusson, A.J.; Wolanski, E.; Cui, J.; Richerme, P. Radial Two-Dimensional Ion Crystals in a Linear Paul Trap. Phys. Rev. Lett.; 2021; 127, 020503. [DOI: https://dx.doi.org/10.1103/PhysRevLett.127.020503]

13. Ivory, M.K.; Kato, A.; Hasanzadeh, A.; Blinov, B. A Paul Trap with Sectored Ring Electrodes for Experiments with Two-Dimensional Ion Crystals. Rev. Sci. Insruments; 2020; 91, 053201. [DOI: https://dx.doi.org/10.1063/1.5145102]

14. Yoshimura, B.; Stork, M.; Campbell, W.C.; Freericks, J.K. Creation of Two-Dimensional Coulomb Crystals of Ions in RF Paul Traps. EPJ Quantum Technol.; 2015; 2, 2. [DOI: https://dx.doi.org/10.1140/epjqt14]

15. Kato, A.; Nomerotski, A.; Blinov, B. Micromotion-synchronized pulsed Doppler cooling of trapped ions. Phys. Rev. A; 2023; 107, 023116. [DOI: https://dx.doi.org/10.1103/PhysRevA.107.023116]

16. Kato, A.; Goel, A.; Lee, R.; Ye, Z.; Karki, S.; Liu, J.J.; Nomerotski, A.; Blinov, B.B. Two-tone Doppler cooling of radial two-dimensional crystals in a radiofrequency ion trap. Phys. Rev. A; 2022; 105, 023101. [DOI: https://dx.doi.org/10.1103/PhysRevA.105.023101]

17. Kiesenhofer, D.; Hainzer, H.; Zhdanov, A.; Holz, P.C.; Bock, M.; Ollikainen, T.; Roos, C.F. Controlling Two-Dimensional Coulomb Crystals of More Than 100 Ions in a Monolithic Radio-Frequency Trap. PRX Quantum; 2023; 4, 020317. [DOI: https://dx.doi.org/10.1103/PRXQuantum.4.020317]

18. Wineland, D.J.; Monroe, C.; Itano, W.M.; Leibfried, D.; King, B.E.; Meekhof, D.M. Experimental Issues in Coherent Quantum-State Manipulation of Trapped Atomic Ions. J. Res. Natl. Inst. Stand. Technol.; 1998; 103, 259. [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/28009379]

19. Guardiola, R.; Navarro, J. On the Lindemann Criterion for quantum clusters at very low temperature. J. Phys. Chem. A; 2011; 115, pp. 2734-2740. [DOI: https://dx.doi.org/10.1021/jp1111313]

20. Robbins, M.O.; Grest, G.S.; Kremer, K. Effect of finite system size on thermal fluctuations: Implications for melting. Phys. Rev. B; 1990; 42, 5579. [DOI: https://dx.doi.org/10.1103/PhysRevB.42.5579]

21. Zheng, X.H.; Earnshaw, J.C. On the Lindemann criterion in 2D. Europhys. Lett.; 1998; 41, pp. 635-640. [DOI: https://dx.doi.org/10.1209/epl/i1998-00205-7]

22. Zhou, M.Y.; Sheng, P. Lattice softening and melting characteristics of granular particles. Phys. Rev. B; 1991; 43, 3460. [DOI: https://dx.doi.org/10.1103/PhysRevB.43.3460] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/9997658]

23. Verlet, L. Computer “Experiments” on Classical Fluids. I. Thermodynamical Properties of Lennard-Jones Molecules. Phys. Rev.; 1967; 159, 98.

24. Berendsen, H.J.C.; Postma, J.P.M.; van Gunsteren, W.F.; DiNola, A.; Haak, J.R. Molecular dynamics with coupling to an external bath. J. Chem. Phys.; 1984; 81, 3684.

25. Dubin, D.H.E.; O’Neil, T.M. Trapped nonneutral plasmas, liquids, and crystals (the thermal equilibrium states). Rev. Mod. Phys.; 1999; 71, 87.

26. Ott, T.; Stanley, M.; Bonitz, M. Non-invasive determination of the parameters of strongly coupled 2D Yukawa liquids. Phys. Plasmas; 2011; 18, 063701. [DOI: https://dx.doi.org/10.1063/1.3592659]