Eigenmodes and thermodynamics of a Coulomb chain in a harmonic potential

The density of ions trapped in a harmonic potential in one dimension is not uniform. Consequently the eigenmodes are not phonons. We calculate the long wavelength modes in the continuum limit, and evaluate the density of states in the short wavelength limit for chains of $N\gg 1$ ions. Remarkably, the results that are found analytically in the thermodynamic limit provide a good estimate of the spectrum of excitations of small chains down to few tens of ions. The spectra are used to compute the thermodynamic functions of the chain. Deviations from extensivity of the thermodynamic quantities are found. An analytic expression for the critical transverse frequency determining the stability of a linear chain is derived.

The density of ions trapped in a harmonic potential in one dimension is not uniform. Consequently the eigenmodes are not phonons. We calculate the long wavelength modes in the continuum limit, and evaluate the density of states in the short wavelength limit for chains of N ≫ 1 ions. Remarkably, the results that are found analytically in the thermodynamic limit provide a good estimate of the spectrum of excitations of small chains down to few tens of ions. The spectra are used to compute the thermodynamic functions of the chain. Deviations from extensivity of the thermodynamic quantities are found. An analytic expression for the critical transverse frequency determining the stability of a linear chain is derived.
Cold atomic gases are systems allowing an extremely high degree of control, opening new directions and challenges both for experiments and theory. The realization of low dimensional ultracold gases enables one to design and investigate, among others, systems with peculiar excitations spectra and thermodynamic properties [1,2]. In this context, Coulomb crystals constitute a unique form of an ordered structure, formed by cold ions in a confining potential, that balances the Coulomb repulsion [3]. The ions vibrate about fixed positions in analogy to the situation in an ordinary solid, while the interparticle distance is usually of the order of several micrometers, constituting an extremely rarefied type of condensed matter [4,5,6,7]. Variation of the potential permits to control the crystal shape as well as the number of ions, allowing to explore crystals of very different sizes, thus offering the opportunity of studying the transition from few particles to mesoscopic systems dynamics [3]. Besides, these structures provide promising applications for spectroscopy [8], frequency standards [9], study and control of chemical reactions [10], and quantum information processors [11,12,13].
In this letter we investigate the dynamics of Coulomb chains, i.e. one-dimensional structures obtained by means of strong transverse confinement and that usually consist of dozens of ions localized along the trap axis [4,5]. The equilibrium charge distribution for such chains is not uniform [14,15]. This is in contrast to the three-dimensional ordering, where the density is uniform and the eigenmodes are phonons. In this letter, the chain excitations are calculated. These are fundamentally different from the phonons in solids because of the non-uniformity of the ion distribution and the long range of the Coulomb interaction. We derive some thermodynamic functions of the chain, and define a specific thermodynamic limit. Interestingly, the linear ion chain differs from systems that are traditionally studied in the framework of statistical physics, since the thermodynamic quantities are not extensive. This behaviour is due to the strong correlation and to the crystal dimensionality. It manifests, for instance, in the dependence of the specific heat on the number of ions, as we show below. The detailed analytical derivation of the results presented here will be published elsewhere [16].
We consider N ions of mass m and charge Q, which are confined by a harmonic trap of cylindrical symmetry and are crystallized along the trap axis at the positions (x (0) j , 0, 0), at which the harmonic force and the Coulomb force, exerted by the other ions, balance. At sufficiently low temperatures the vibrations around these points can be considered harmonic and described by [17,18] where and y i are the displacements in the axial and transverse directions (the equations for z i are similar to the ones for y i ), while the harmonic confinement is characterized by the axial and transverse frequencies ν, ν t , respectively (ν t ≫ ν). The coupling matrix is with w i = q i , y i , z i . Equations (1-2) describe a system of coupled oscillators, with long range interaction and position-dependent coupling strength. In what follows we calculate the eigenmodes of this system. For this we assume w i (t) = e iωtw i (ω)dω/2π. To simplify notations we replacew i (ω) by w i . This results in equations for the eigenmodes of frequency ω that are similar to (1-2), but withẅ i replaced by −ω 2 w i . The problem reduces to the eigenvalue equation where n = 1, . . . , N labels the mode. The axial and transverse eigenfrequencies are given by the relations ω 2 n = ν 2 + λ n and ω ⊥2 n = ν 2 t − λ n /2. Since λ n > 0, the axial collective excitations are of higher frequencies than the trap frequency ν, while the transverse excitations are of lower frequencies than ν t . In general, Eqs. (1-2) are invariant under reflection with respect to the center of the trap, which coincides with the origin of the axes. Hence, the normal modes of the chain are symmetric (even) or antisymmetric (odd) under reflection with respect to the center, w i = ±w N −i [18]. It can be easily verified that the center-of-mass motion is an even eigenmode, w (1) , of the Eqs. (1) and (2) at frequency ν, ν t , respectively. Remarkably, also the first odd eigenmodes can be exactly determined. They are w 3ν for the axial, ν 2 t − ν 2 for the transverse excitation, independent of the number of ions, as can be found by substitution of this ansatz in Eq. (3)(4), combined with the equilibrium condition [16].
In order to solve Eq. (4), we assume N ≫ 1. In this limit, we may consider the interparticle spacing a smooth function of the position, which is inversely proportional to the density of ions per unit length n L (x i ) = 1/a L (x i ). The density n L (x) is determined assuming a uniform distribution of charges inside an oblate ellipsoid of axial length 2L, and takes the form n L (x) = 3N/4L 1 − x 2 /L 2 , which is defined for |x| ≤ L, while the length L is found by minimizing the energy of the crystal, and fulfills the relation L(N ) 3 = 3 Q 2 /mν 2 N log N [14]. These approximations provide a good estimate of the ions distribution in the center of the chain for N sufficiently large. In fact, n L (x) and L(N ) used here are the leading terms of an expansion in powers of 1/ log N [14]. For N ≫ 1 one can approximate the chain by a continuous distribution of charges for the evaluation of the long-wavelength modes. In this limit, x (0) i is a continuous variable in the interval (−L, L), the displacement w i = w(x 0 i ) is a continuous function denoted by w(x), and the sum in (3) is approximated by an integral, where the density n L (x) appears as a weight function. In particular, the eigenmodes fulfill the orthogonality relation Using an integration by parts the sum (3) can be approximated by where ξ = x/L and The approximation is valid at leading order in log N and sufficiently far away from the edges of the chain. The eigenfunctions of Eq. (5) are the Jacobi polynomials P 1,1 ℓ (ξ) [19] with the eigenvalues λ Jac n = −ℓ(ℓ + 3)ν 2 /2 with ℓ = 0, 1, . . . and n = ℓ + 1. Substituting this result into the eigenvalue equations (3)(4), one finds ω Jac n = ν n(n + 1)/2.
Analogously, the eigenfrequencies of the transverse modes are ω ⊥ Jac n = ν 2 t − ν 2 (n − 1)(n + 2)/4. It should be noted that the Jacobi polynomials at ℓ = 0, 1 are exact solutions of Eqs. (1)(2). In general, they are the correct eigenmodes of Eq. (4) in the asymptotic limit N → ∞. For finite N a perturbative expansion in the parameter 1/ log N is required, resulting in a very slow convergence. Nevertheless, the behaviour λ n ∝ ν 2 is exact in the continuum limit [14,16]. Figure 1 presents the comparison between the spectrum of eigenfrequencies, calculated numerically for 1000 ions, and the solutions ω Jac n , ω ⊥ Jac n . Figure 2 shows the relative deviation of the frequencies (6) from the numerical results for chains of different number of ions. From Fig. 2(a) one sees that Eq. (6) gives already a good estimate of the eigenfrequencies for a chain of 10 ions. Its prediction is valid for the axial low-and the transverse high-frequencies, and it improves slowly as N increases, as expected from the slow convergence of the 1/ log N expansion.
The short wavelength eigenmodes are characterized by relatively large displacements of the ions at the chain center, while the ions at the edges nearly do not move. In fact, the interparticle distance is minimal around the center, and for N ≫ 1 it is of order 1/N , while it is significantly larger at the edges. Hence, a wave cannot propagate to the edges, where the interparticle spacing is larger than the wavelength. Moreover, in the center and for the short-wavelength excitations we expect that the relevant contributions to the force on an ion originate from the neighbouring charges, as nearby groups of ions move in opposite directions, resulting in forces that mutually cancel. By this hypothesis we keep in (3) only the nearest-neighbour interaction and apply the formalism developed by Dyson [20] to derive the density of states as a function of the squared frequencies λ. This is given by where Λ(x) = 2Q 2 n L (x) 3 /m is the nearest-neighbour coupling constant for slowly-varying interparticle distance, and f (λ) = L 1 − (λ/4Λ(0)) 1/3 [16]. Here, Λ(0) ≈ 9ν 2 N 2 /32 log N at leading order in log N . The spectrum obtained from Eq. (7) is shown in Fig. 1. The deviations from the numerical result are due to the assumption of nearest-neighbour coupling, and are small in the short-wavelength part of the spectrum, showing that Eq. (7) provides a good approximation in this regime. In particular, it provides a good estimate of the eigenvalue λ N , which determines the maximal axialand minimal transverse-frequencies through the relations ω max = √ ν 2 + λ N and ω ⊥ min = ν 2 t − λ N /2. This eigenvalue is found from f (λ) = 0, as for λ > λ N the function D(λ) vanishes, and λ N = 9ν 2 N 2 /(8 log N ) at leading order in log N . The minimal transverse frequency is It can vanish for certain values of ν, ν t and N . We identify the critical value ν t,cr = λ N /2, such that for ν t < ν t,cr the linear chain is unstable. The critical aspect ratio between the trap frequencies α cr = ν 2 /ν 2 t,cr takes the value α cr = 16 9 log N N 2 (9) which fixes the condition on ν t for the chain stability according to the inequality ν 2 t > α −1 cr ν 2 . The result (9) is in agreement with the analytical estimate in [21], which was obtained under different requirements. Figure 3 shows that it also approximately reproduces the curve obtained by fitting points calculated by molecular dynamics simulation [22], which has been experimentally verified for  3: αcr(N ) as a function of the number of ions. The solid curve gives the result according to (9). The dashed line is the curve 2.53 N −1.73 , that fits the values obtained by molecular dynamics simulations [22,23].
chains of tens of ions [23]. It was suggested that this instability, resulting in a structural transition from a linear to a zig-zag equilibrium configuration [4,22], can be treated as a phase transition [22,23].
Using these results we now discuss the thermodynamics of a linear crystal of N ions, whose dynamics are described by 3N harmonic oscillators of frequencies ω n , ω ⊥ n , in the regime where the axial frequencies are much smaller than the transverse excitations, i.e. ω max ≪ ω ⊥ min . The eigenmodes are quantized using standard procedures [17]. The chain is assumed to be in thermal equilibrium with a bath at temperature T , which is sufficiently low so that only the axial modes are excited (k B T ≪hω ⊥ max ). Hence, the dynamics of the system is one-dimensional. The thermodynamic limit for this kind of system can be defined by assuming constant interparticle spacing (thus constant linear density) at the center of the chain as N increases, namely requiring a L (0) to be constant, in analogy with the definition for cold gases in traps [24]. Since a L (0) ∝ √ log N /νN 2/3 , this requirement corresponds to a vanishing axial frequency as N → ∞, according to ν ∼ √ log N /N . Using this definition, ω mac and ω ⊥ min take finite values in the thermodynamic limit, since they depend on N and ν only through the combination νN/ √ log N . We identify the thermodynamic variables with T , N , and ν, whose variation results in a variation in the length of the chain [3], and take constant ν t . First we study the heat capacity C a = ∂U th /∂T | ν,N , where the thermal energy U th is the average energy of the excited state relative to the ground state energy. At high temperatures C a exhibits the Dulong-Petit law C a = N k B , while at low temperatures it takes the form where x 0 = βhν and we have used the density of states derived from (6). For k B T ≫hν we can set x 0 ∼ 0 and obtain C a ∝ T . Hence, in this temperature regime the heat capacity is linear in T , like for a one-dimensional Debye crystal [25]. Remarkably, for these temperatures C a is inversely proportional to ν. Hence, the specific heat per particle c a = C a /N behaves like c a ∼ 1/N ν and vanishes as c a ∼ 1/ √ log N in the thermodynamic limit defined above. Thus, at low temperatures it depends on the number of ions, thereby manifesting a deviation from extensivity of the system's behaviour.
The coefficient of thermal expansion α T also exhibits a remarkable behaviour as a function of the temperature. It is related to C a and to the isothermal compressibility κ T according to the relation α T = 3κ T C a /2L [25]. In the regime of thermal stability, the isothermal compressibility does not vanish and is practically independent of the temperature, the pressure being dominated by the zero-temperature component. In particular, it scales as κ T ∼ 1/ log N . Hence, the temperature dependence of α T is solely determined by the heat capacity C a . At low temperatures C a ∼ 1/ν, and as the thermodynamic limit is approached α T ∼ (log N ) −3/2 . At higher temperatures, when the heat capacity manifests the Dulong-Petit behaviour C a = N k B , the coefficient of thermal expansion vanishes like α T ∼ 1/ log N . For large but finite number of ions α T does not vanish, in contrast to the behaviour found in a uniform harmonic solid [25]. The coefficient of thermal expansion α T is plotted in Fig. 4 as a function of T for a chain of 1000 ions.
Finally, we remark on the regime of thermal stability, which holds when the thermal energy is considerably smaller than the equilibrium energy, or, equivalently, the displacements are much smaller than the respective interparticle distances. The condition for thermal stability is Q 2 log N/a L (0) ≫ k B T . Outside this regime, thermal excitations may cause structural transitions, accompanied by a critical behaviour of the thermodynamic functions.
In conclusion, we studied the low temperature dynamics of a strongly-correlated and one-dimensional crystal inside a harmonic potential. We investigated how the dynamics of the system is affected by the dimensionality and the long-range correlations and how this is manifest in the thermodynamics functions. In future works we will explore the behaviour of these functions at the instabil-ity points, in order to characterize the crystal structural transitions and their relations to the standard theory of phase transitions [26].
The excitations spectra we presented are of interest for spectroscopy and quantum information with ion traps, where quantum logic is implemented with the lowest collective modes [11,12,13]. This work may be relevant to studies of Coulomb systems like ion crystals in storage rings [2] and cold neutral plasmas [27], which at sufficiently low temperatures are expected to crystallize [28], and contributes to the on-going research on lowdimensional cold gases [1].
It is our great pleasure to thank Andreas Buchleitner for hospitality at MPIPKS in Dresden, where most part of the research was done, and J. Eschner for careful reading of the manuscript. Support by BSF, the Minerva Center of Nonlinear Physics of Complex Systems, the Center for Theoretical Physics at the Technion, QUEST, and QGATES, is acknowledged.