Dynamics of an ion chain in a harmonic potential

Cold ions in anisotropic harmonic potentials can form ion chains of various sizes. Here, the density of ions is not uniform, thus the eigenmodes are not phononic-like waves. We study chains of N>>1 ions and evaluate analytically the long wavelength modes and the density of states in the short wavelength limit. These results reproduce with good approximation the dynamics of chains consisting of dozens of ions. Moreover, they allow to determine the critical transverse frequency required for the stability of the linear structure, which is found in agreement with results obtained by different theoretical methods [D.H.E. Dubin, Phys. Rev. Lett. 71, 2753 (1993)] and by numerical simulations [J.P. Schiffer, Phys. Rev. Lett. 70, 818 (1993)]. We introduce and explore the thermodynamic limit for the ion chain. The thermodynamic functions are found to exhibit deviations from extensivity.


I. INTRODUCTION
Coulomb crystals are organized structures of charged particles, which interact through the Coulomb repulsion and organize in regular patterns at sufficiently low temperatures in presence of a confining potential [1]. These potentials are realized by means of Paul or of Penning traps [2], and their geometry determines the crystal's structure. Several remarkable experiments have reported crystallization of ion gases in Paul traps [3,4,5,6,7], Penning traps [8] and ion-storage rings [9]. The Braggscattering in three-dimensional structures was studied providing information about the internal structure of the crystal [10,11]. These crystals represent a kind of rarefied condensed matter, the inter-particle distance being of the order of several micrometers, allowing to study the structure by means of optical radiation. Variation of the potential permits to control the crystal shape as well as the number of ions, thus offering the unique opportunity of studying the transition from few particles to mesoscopic systems. Besides, these structures have been object of growing interests as they provide promising applications for spectroscopy [12,13], frequency standards [14], study and control of chemical reactions [15], and quantum information processors [16,17,18,19].
In this work we investigate the dynamics of Coulomb chains. These are one-dimensional structures, obtained by means of strong transverse confinement and that usually consist of dozens of ions localized along the trap axis [5,6]. They represent a peculiar crystallized structure: In fact, due to the axial potential the equilibrium charge distribution is not uniform [20,21]. This is in contrast to the three-dimensional case, where the density of charges in a harmonic potential is uniform and, therefore, where the eigenmodes are phononic-like waves. In the Coulomb chain the non-uniformity of the density of ions combined with the long-range interaction result in excitations that are fundamentally different from the phonons in solids and lead to interesting thermodynamic properties. The exploration of these excitations and of the chain thermodynamics is the subject of the present paper.
Our starting point is the ions equilibrium configuration evaluated in [20]. We investigate the dynamics for small oscillations, when the harmonic approximation is valid, in the limit of large number of ions. We derive analytically the eigenfrequencies and the corresponding eigenmodes for the long wavelength excitations. These are compared with numerical results and good agreement is found. From the resulting dispersion relation the density of states of the long wavelength eigenmodes is determined. An analytic form of the density of states is also found for the short wavelength excitations. This result allows to evaluate the critical aspect ratio between the frequencies of the transverse and the axial confining potential, that determines the stability of the chain. The value we find agrees with numerical results [22], which have been experimentally verified for chains of few ions [23]. In particular, it is in agreement with the analytical estimate in [24], which was obtained under different requirements.
Using these results we discuss the statistical mechanics of the chain, and derive some thermodynamic quantities in a specific thermodynamic limit, that is defined here by keeping constant the density of ions in the chain center, as the number of ions tends to infinity and the axial frequency to zero, analogously to the definition for cold neutral gases in traps [25]. Non-extensive thermodynamic properties are found. We compare the thermodynamic functions with the ones of a chain of finite number of ions that are obtained numerically, and find reasonable agreement.
This work is organized as follows. In section II we introduce the basic equations and discuss the fundamental properties. In section III the spectrum of excitations is studied. In section IV we investigate the statistical mechanics of the system. Section V presents the conclusions and the outlook. In the appendices, several details of the calculations of section III are reported.

II. A STRING OF CHARGES IN A HARMONIC POTENTIAL
The Hamiltonian describing the dynamics of a chain of N ions of mass m and charge Q, which are confined by a harmonic potential, is given by where r j = (x j , y j , z j ), p j are the positions and conjugate momenta (j = 1, . . . , N ). The term V accounts for the oscillator's potential and the Coulomb repulsion, where the harmonic oscillator has rotational symmetry around the x-axis with axial and transverse frequencies ν, ν t , respectively.
For sufficiently low kinetic energy crystallization occurs. The temperature, at which the gas is crystallized, corresponds to large plasma parameters Γ = Q 2 /a WS k B T ≫ 1. Here, a WS is the Wigner-Seitz radius, that is a function of the mean density n and is defined as a WS = (3/4πn) 1/3 [1]. In this regime the ions are localized at the classical equilibrium positions r (0) j , that satisfy the equations ∂V /∂r j | r (0) j = 0, and such that the potential energy is minimal. When the harmonic potential is sufficiently asymmetric, i.e. for ν ≪ ν t , the ions equilibrium positions are confined to the trap axis [22], namely r j , 0, 0), and satisfy the equation describing the equilibrium of the forces, where the numbering convention is x i > x j for i > j. The stability of these points with respect to the transverse vibrations depends on the number of ions N and on the ratio ν t /ν, and it is discussed in Section III B. In this section, we assume that the configuration is stable and approximate the potential by its second order Taylor expansion around the points r (0) j . We denote by the displacements in thex-direction, and approximate the Hamiltonian (1) as N ) is the classical minimum energy, while H har describes the (classical) harmonic oscillations around the equilibrium points [26,27,28], and the coefficients K i,j = ∂ 2 V /∂x 2 j | {x 0 l } are positive and take the form Equation (4) shows that the axial motion is decoupled from the transverse motion in the harmonic expansion. The corresponding equations of motion arë and describe a system of coupled oscillators, with long range interaction and position-dependent coupling strength. In this paper eigenmodes will be calculated. For this we assume q i (t) = e iωtq i (ω)dω/2π. To simplify notations we replaceq i (ω) by q i . This results in equations for the eigenmodes of frequency ω that are similar to (6), but withq i replaced by −ω 2 q i . The same replacement will be performed for y i (t) and z i (t).
It can be easily verified that the center-of-mass motion is an eigenmode of the secular equations (6)(7)(8). The axial center-of-mass mode is q 1 = . . . = q N at the characteristic frequency is ν, while the transverse center-of-mass modes are y 1 = . . . = y N and z 1 = . . . = z N at frequency ν t . We remark that the axial and the transverse coupling terms appearing in Eqs. (6)(7)(8) have opposite signs. Due to this property, in the axial direction the collective excitations are of higher frequencies than the center-of-mass frequency ν, while in the transverse plane the collective excitations are of lower frequencies than ν t .

A. Properties and symmetries
Hamiltonian (4) is not translationally invariant, and this is a consequence of the non-uniformity of the ions equilibrium distribution, due to the harmonic force appearing in Eq. (3). The Hamiltonian (4) is however invariant under reflection with respect to the center of the trap, which coincides with the origin of the axes. In particular, where i = 1, . . . , N ′ (here, N ′ = N/2 for even N , while N ′ = (N − 1)/2 for odd N ). Hence, the normal modes of the chain are symmetric (even) or antisymmetric (odd) under reflection with respect to the center [27,28], such that i , and n labels the mode. Some general properties can be inferred from this simple consideration. For instance, the even modes of the axial motion are characterized by constant length L of the chain, since q (n) For the odd modes, on the other hand, the center of mass of the chain, which coincides with the chain center, does not move. Clearly, the center-of-mass mode, which we denote by w (1) i , is an even mode characterized by equal displacements at the positions x It is remarkable that also the lowest axial odd mode (stretch mode) and its frequency can be exactly determined. In fact, taking q where we have used relation (3). Therefore, the frequency of the axial stretch mode q (2) i is √ 3ν and its value is independent of the number of ions N of the chain. This property has been first demonstrated in [29]. It has also been observed by numerical evaluation of the normal modes of chains up to 10 ions [17,26]. Analogously, the transverse stretch mode, which is the highest odd transverse excitation, satisfies y We remark that the invariance under reflection imposes different boundary conditions than the ones that are usually chosen for a crystal with uniform ion distribution. In a crystal that is translationally invariant even and odd modes are degenerate and one may choose periodic boundary conditions [30]. In presence of an external potential with central symmetry this invariance is broken, apart for the mirror symmetry with respect to the center. Hence, at the edges the eigenmodes fulfill the relation w −N ′ = ±w N ′ , where the sign is determined by the parity.

III. THE SECULAR PROBLEM
The systematic derivation of an analytic solution of Eqs. (6-8) is a challenging problem, since it requires to take systematically into account the position dependent coupling constant and the long range interaction. Nevertheless, in the limit of large number of ions N ≫ 1 we can make some simplifying assumptions. In this limit, in fact, the interparticle spacing a L ( is a smooth function of the position, and it is inversely proportional to the density of ions per unit length [20], The density of charges for unit length can be evaluated by applying the Gauss theorem to a continuous distribution of charges, that is assumed to be uniformly distributed in an elongated ellypsoid. The resulting one dimensional density is [24] n L (x) = 3 4 which is defined for |x| ≤ L, where 2L is the length of the crystal at equilibrium. The density (12) gives a good estimate of the charge distribution in the center of the chain for N sufficiently large [20]. The length L is evaluated by minimizing the energy of the crystal, and at leading order in N fulfills the relation [20] In the following, we use these quantities to derive an approximate solution for the long-and short-wavelength modes in the limit of N ≫ 1 ions. Furthermore, we compare the results of the derivation with the numerical calculations, obtained by solving the eigenvalue equations (6-8) for a finite number of ions.
A. The eigenmodes in the long wavelength limit We use the ansatz q i (x, t) = e iωtq i (x) in Eqs. (6) and define the rescaled positions ξ i /L and the rescaled interparticle distances a(ξ i ) = a L (x i )/L. With these definitions, denoting for simplicity of notationq i → q i , Eqs. (6) take the form where we have introduced the dimensionless constant If the number of ions is large (N ≫ 1), for the longwavelength modes one can approximate the chain by a continuous distribution of charges. In this limit, ξ is a continuous variable varying in the interval (−1, 1), while the displacement q i = q(ξ i ) is a continuous function, here denoted by q(ξ). Then, Eqs. (14) take the form where while n(ξ) = 1 − ξ 2 is the density of charges normalized to 4/3. Equations (16) and (17) are valid away from the edges of the chain and for long-wavelength excitations, where the continuum approximation is reasonable. The continuum limit for Eqs. (7) gives A similar type of equations is obtained for z i . It is remarkable that the axial trap frequency enters only as a prefactor on the right hand side of Eqs. (16). Consequently the axial eigenfrequencies are proportional to ν. The transverse eigenfrequencies, instead, do not show this behaviour, as one can see from Eq. (18): Here, it is the quantity ω 2 − ν 2 t which is proportional to ν 2 . The proportionality constants depend on the modes and will be calculated in what follows within some approximations. The results are independent of the charge Q and of the mass m.
According to Eqs. (16)(17)(18), the secular problem consists of solving the eigenvalue equation where w (n) (ξ) can be the axial or the transverse modes, that satisfy the orthogonality relation By partial integration Eq. (17) can be written as the sum of two terms where I 0 contains the contributions of the ions around the point ξ, while the term ∆I is determined by the value of the density and of the eigenmode function q(ξ) and their derivatives at the end points of the chain. In appendix A we derive the explicit form of the two terms and discuss their order of magnitude. For the long wavelength modes and at the points ξ sufficiently far away from the chain end-points ξ = ±1, we find I ≈ I 0 where and a(ξ) is an infinitesimal quantity of the order 1/N . Using (11) and (12) in (21) we obtain where In the limits of validity of Eq. (21) and for N sufficiently large, such that log N ≫ 1, Eq. (22) can be approximated by the leading order in log N and the eigenvalue equation (19) reduces to where λ n = − log Nλ n . Equation (25) is the differential equation fulfilled by the Jacobi polynomials P 1,1 ℓ (x) at the eigenvalues [31]λ with ℓ = 0, 1, . . . and n = ℓ + 1. After substitution of (23-26) into (19), the eigenfrequencies of the axial excitations are found from (15)(16) and take the form Analogously, the eigenfrequencies of the transverse modes are obtained from Eq. (18) with (19), resulting in It should be noted that the solutions of (23) for ℓ = 0 and ℓ = 1 are exact solutions of the original problem (14): The corresponding eigenmodes w (1) (x) =constant (center of mass) and w (2) (x) ∝ x (stretch mode), which are the continuum limits of the eigenmodes we have found for the discrete case, are in fact Jacobi polynomials. The result for the center of mass is obvious. The exact result for the stretch mode can be understood, noting that   (27) and (28) derived for the axial and the transverse modes. Figure 2(a) exhibits the part of the spectrum with the long-wavelength axial modes: Here, one sees that Eq. (27) approximates well the lowest part of the axial spectrum, where the limit of continuous charges distribution is reasonable. The short-wavelength modes are not reproduced correctly by (27). These are better evaluated by using a more proper approximation for this regime, and will be discussed in the next section.
We remark that, apart for the first two eigenmodes, the Jacobi polynomials describe the eigenmode excitation at leading order in log N and near the center of the chain, where the interparticle separation is of order 1/N and the distribution of charges can be treated as a continuum for sufficiently long wavelengths. The continuum approximation fails at the edges, where the interparticle spacing is significantly larger and Eq. (12) is not meaningful. In particular, Eq. (13) gives the upper bound for the chain length, which would be obtained in the limit of N ≫ 1 particles. Hence, a reasonable boundary condition is to assume that the eigemodes and their derivatives vanish at x = ±L, where there are no charges and hence the energy density is zero. The solution (23) taken at the center of the chain neglects the charges at the edges, on the basis of the observation that there the number of ions is much smaller than at the center, and their contribution to the integral (17) can therefore be neglected.
The evaluation of the correction to the results (27) and (28) should be done in perturbation theory in the parameter 1/ log N , following an analogous procedure to the one applied in [20] for evaluating the correction to the density of ions (12) and to the equilibrium length (13). In practice, this expansion has a very slow convergence, and does not allow for a simple analytical expression. Nevertheless, the comparison with the spectrum evaluated numerically, by solving equations (3) and (6)(7)(8), shows that Eq. (27) and (28) give already a good estimate of the eigenfrequencies for a chain of 10 ions, as it can be seen in Fig. 3(a), where the relative deviation of the frequency given by Eq. (27) from the numerical result is plotted for chains with different number of ions. This agreement is in general valid, respectively, for the axial low modes and the transverse high modes, and exhibits a very slow improvement as the number of ions in the chain increases, due to the slow convergence of the 1/ log N expansion.
B. The density of states in the short-wavelength limit Simple physical considerations show that the shortwavelength eigenmodes are characterized by relatively large displacements of the ions around the center of the chain, while the ions at the edges nearly do not move. In fact, the interparticle distance is minimal in the middle of the chain, and for N ≫ 1 it is of order 1/N , while it is consistently larger at the edges. Hence, a wave cannot propagate in a region where the interparticle spacing is larger than the wavelength. Here, we can make some simplifying assumptions in solving the eigenvalue equations (6)(7)(8) for the short wavelength modes. In fact, in the center of the chain we expect that the relevant contributions to the force on an ion originates from the neighbouring charges, as nearby groups of charges move in opposite directions, resulting in forces that mutually cancel. By this hypothesis, in (6-8) we keep only the nearest-neighbour interaction, so that the equations to solve arë where the tridiagonal matrix Λ is defined by its action on the vector (w 1 , . . . , and Λ i = K i,i+1 /m where K i,j is given in (5). The matrix Λ is symmetric, and we denote its characteristic frequencies by −ω 2 . Following the derivation of Dyson [34], that is summarized in appendix B, the density of states D(ω 2 ) for N ≫ 1 is found from the characteristic function Ω(θ) according to while Ω(θ) is explicitly evaluated by using the properties of antisymmetric matrices and takes the form Here, ζ(θ, j) is the infinite continued fraction, andΛ 2i−1 =Λ 2i = Λ i . For N ≫ 1 and around the center of the chain we may assume Λ(ξ) to be a slowly-varying function of the position ξ, such that Λ i+1 = Λ i + δΛ i and δΛ i /Λ i ≪ 1. This allows to evaluate explicitly ζ(θ, j) at first order in δΛ i . For this purpose, we defineΛ j−1 =Λ j = Λ, andΛ j+1 =Λ j+2 = Λ + δΛ, where δΛ is the first order variation. Furthermore, we denote ζ(θ, j) = ζ and assume ζ(θ, j + 2) = ζ + δζ, where δζ is a first order variation. We substitute these quantities into (35) keeping only the terms up to first-order and look for a consistent solution. The resulting equation is which is integrated to Here, we have taken the integration constant to be zero since at the boundaries of the chain Λ → 0. The resulting solution has the form leading to the characteristic function where we have used the rescaled variable ξ and the fact that the integrand is even in the interval (−1, 1). Equation (39) corresponds to the continuum limit of the discrete summation in Eq. (34), and it is valid away from the edges for N ≫ 1. Substituting (39) into (33) we obtain the equation for the density of states as a function of the physical parameters, where f (z) with at leading order in log N .
Equation (40) withω 2 = 1/z gives the density of states of the short wavelength modes. The corresponding spectrum 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. (40) provides a good approximation in this regime. In particular, this result allows to evaluate explicitly the value of the maximal axial frequency ω max , the minimal transverse frequency ω ⊥ min , and the spectrum of the eigenfrequencies in their neighbourhood. The maximal axial frequency and the minimal transverse frequency are found from the maximal value of z for which the integrand in (40) is real, corresponding to f (z) = 0. For larger values of z the density of states vanishes. The corresponding eigenvalue of the matrix Λ isω 2 = 4Λ 0 . From (29-31) one finds ω max = √ ν 2 + 4Λ 0 and ω ⊥ min = ν 2 t − 2Λ 0 . Therefore, the largest value of the axial frequency is determined by the largest value of the spring constant, that is the value of the spring constant at the center of the chain, and at leading order in log N is Analogously, the smallest value of the transverse modes frequency is From Eq. (44) one sees that ω ⊥ min can vanish for certain values of ν, ν t , and N . We denote by ν t,cr = 9 16 the value of the transverse frequency, such that for ν t < ν t,cr the linear chain is unstable with respect to excitations of the transverse vibrations. Using the notation introduced in [22] we define the critical value of the aspect ratio between the trap frequencies α cr = ν 2 /ν 2 t,cr . It takes the value α cr = 16 9 log N N 2 which fixes the condition on ν t for the chain stability according to the inequality ν 2 t > α −1 cr ν 2 , and it is in agreement with the analytical estimate in [24], which was obtained under different requirements. In Fig. 4 we compare the result (46) with the relation cN β1 , with β 1 = −1.73 and c = 2.53, obtained in [22] by fitting points calculated with molecular dynamics simulations, and verified experimentally in [23]. Our result reproduces approximately this curve. At α cr the crystal undergoes a structural transition from a linear string to a zig-zag configuration, as it has been observed in [6]. It is interesting to ask whether this structural transition can be considered as a phase transition, as discussed in [22,23,24]. A systematic investigation in this direction requires a proper definition of the thermodynamic limit for this kind of system. A natural thermodynamic limit, that will be discussed in the next section, is one where as N → ∞, the axial frequency ν → 0 so that the density in the center is fixed. From (12) and (13) this requires that the ratio ν 2 N 2 / log N is kept constant. From (45) it is clear that in this limit the critical transverse frequency tends to a well defined value. The exploration of the properties of this transition in the thermodynamic limit will be  [22] (see also [23]).
object of further studies. Particularly interesting is the comparison with standard phase transitions [35].

C. The phonon-like approximation
It is natural to introduce a phonon-like approximation for the eigenmodes of Eqs. (6)(7)(8). In this approximation, a phonon-like solution is superimposed by a slowly varying amplitude, that takes into account the slow variation of the coupling strength K i,j of Eq. (5) as a function of both i and j. This approximation, that is outlined in appendix C, is reasonable for a relatively wide range of long wavelength excitations compared to the Jacobi polynomials solution, discussed in Sec. III A (see Fig. 7). It is inferior compared to the Jacobi polynomials in the very long wavelenght limit, and it is a bad approximation for the short wavelength regime. Therefore, the phononlike approximation, that is natural in condensed matter physics, is not an asymptotic approximation in the thermodynamic limit N → ∞, neither for the long wavelength nor for the short wavelength part of the spectrum.

IV. STATISTICAL MECHANICS
In this section we use the density of states, that was evaluated in the preceding section, in order to derive the thermodynamic quantities of a linear crystal of N ions. The linear chain is assumed to be in the regime of stability and to be in equilibrium with a thermal bath at temperature T . The oscillations around the classical equilibrium points are quantized using standard procedures [26]. It should be noted that in this limit the quantum statistics of the atoms are irrelevant: In fact, the single particle wavepacket is much smaller than the interparticle distance [36]. The dynamics of the crystal is thus described by 3N harmonic oscillators of frequencies ω n , ω ⊥ n , where the frequencies ω ⊥ n are doubly degenerate. It is modeled by the HamiltonianĤ, obtained from H after quantizing the eigenmodes of H har in (4). Here,Ĥ = H 0 +Ĥ n , where H 0 is the ground state energy, ω n + 2ω ⊥ n whileĤ n describes the contribution of the collective excitationŝ whereN n ,N ⊥ n,y ,N ⊥ n,z are the operators counting the number of excitations. The term V 0 corresponds to the classical minimum energy of the Coulomb crystal. For an infinite chain, it is obtained by minimizing the classical energy with respect to the length L using the density of charges in Eq. (12), and it is evaluated to be [20] where L(N ) is given in (13). We remark that we investigate the thermodynamic quantities for crystals characterized by a finite number of particles N and finite axial frequencies ν, i.e. crystals of finite size, that may be close to experimental situations. It is however instructive to consider the definition of the thermodynamic limit for this kind of system characterized by strong correlations, where the effect of the charges at the edges cannot be neglected a priori in evaluating the statistical properties. Here, the thermodynamic limit can be appropriately defined by assuming constant interparticle spacing (thus constant linear density) at the center of the chain x = 0, namely requiring a L (0) to be constant. Denoting by a 0 = a L (0), it scales as where g = (3Q 2 /m) 1/3 . This requirement corresponds to a vanishing axial frequency, according to ν ∼ √ log N /N , as N → ∞. With this definition, in the thermodynamic limit the maximum axial frequency (43) and the critical transverse frequency (45) are independent of N and of ν, taking the values ω max = 3(g/2a 0 ) 3/2 and ν t,cr = 3/4(g/a 0 ) 3/2 . In the following, we derive the thermodynamic quantities for ion chains of finite size characterized by a fixed and finite value of the number of particles N and of the axial frequency ν, and discuss how these quantities behave when we take the specific thermodynamic limit that was defined above.
Assuming thermal equilibrium with the bath, the state of the system at constant number of ions N is described by the density matrix of the canonical ensemble where β = 1/k B T and Z is the partition function, which determines the Helmoltz free energy We identify the thermodynamic variables with T , the temperature, N , the number of atoms, and ν, the axial trap frequency, whose variation corresponds to a variation in the length of the chain [1]. These are not a complete set of thermodynamic variables, but they fully determine the state of the crystal for the thermodynamic quantities we investigate in the following. In particular, we take ν t as constant and assume ω max ≪ ω ⊥ min , i.e. that there is a large gap between axial and transverse excitations. In this limit, we consider temperatures T such that k B T ≪hω ⊥ min . In this regime the transverse modes can be considered frozen, hence the contribution due to their excitations to the crystal's thermodynamic properties can be neglected, and the dynamics of the system is one-dimensional. The thermal energy of the crystal is given by The heat capacity C a = ∂U th /∂T | ν,N is The behaviour at high temperatures, such that k B T ≫ hω max (but k B T ≪hω ⊥ min ), is given by the Dulong-Petit law C a = N k B as is clear from (50). On the other hand, for k B T ≪hν all modes are frozen and the energy of the chain is the zero point energy H 0 . For large number of particles N ≫ 1 we can approximate the sum in (50) by the integral where g(ω) = ∂n/∂ω is the density of states. In particular, at temperatures such that the contribution of the long-wavelength excitations to the sum is predominant and is given by (27), the resulting density of states is and the heat capacity is given by the integral where we have defined x 0 = βhν. Hence, the integrand and the integration limits depend on the temperature through x 0 . In this regime and for k B T ≫hν we can set x 0 ∼ 0 in (51), and recover the result C a =cT , with c = √ 8ζ(2)k 2 B /hν, where ζ(2) is the Riemann's zeta function. Therefore, for the considered regime the heat capacity is proportional to the temperature, which is a characteristic behaviour encountered in a one-dimensional Debye crystal [30]. In Fig. 5 the specific heat c a = C a /N for N = 1000 ions is plotted as a function of the temperature T . In the inset, the low temperatures behaviour is shown, and compared with the curvecT /N estimated above. The figure shows that the evaluated behaviour, valid in the asymptotic limit of an infinite number of ions, provides a reasonable description of the specific heat at low temperatures.
It is remarkable that the heat capacity in (51) scales like C a ∼ 1/ν. The specific heat per particle c a = C a /N behaves thus like c a ∼ 1/N ν at low temperatures, and in the thermodynamic limit it vanishes as c a ∼ 1/ √ log N . It thus depends on the number of ions, and this is a manifestation of the deviation from extensivity of the system's behaviour. Note that the relative energy fluctuations are of the order ( √ log N /N ) 1/2 , therefore the usual equivalence of ensembles holds.
The pressure P in the axial direction is defined as and it is the variation of the free energy with the length of the string at constant N and T . Under these conditions, the length of the crystal is varied by changing the axial trap frequency ν, according to as obtained from (13). Substituting the explicit form of the free energy into (52), we find is the pressure at zero-temperature, and P T gives the contribution of the excitations, In deriving (54) and (55) we have used ∂V 0 /∂L = −V 0 /L that results from (47), and the relations that are obtained from the functional dependence of ω n and ω ⊥ n on ν, as can be extracted from (16) and (18). Note that a variation of the axial trap frequency implies also a variation of the transverse eigenfrequencies, which give a contribution of opposite sign to the pressure, as it is obvious from (54). However, the contribution due to the quantum mechanical zero-point energy is very small compared to the classical term V 0 /L. Therefore, P 0 is dominated by V 0 /L and, using (47), in the thermodynamic limit it scales as P 0 ∼ log N . The term P T depends on the temperature. For low temperatures, such that k B T ≫hν, it scales as P T ∼ 1/Lν ∼ 1/ √ log N . At high temperatures, in the Dulong-Petit regime, P T does neither depend on N nor on ν. In the regime where the chain is thermally stable, which we consider here, V 0 ≫ U th giving P ≈ P 0 . Thus, the pressure is dominated by the zero-temperature contribution and in the thermodynamic limit P ≈ P 0 ∼ log N . A useful relation for the following discussion is The isothermal compressibility κ T is evaluated from the pressure according to where B is the bulk modulus. Using (54), (55) in (57) we find In the thermodynamic limit the bulk modulus B is dominated by the zero-temperature contribution −L∂P 0 /∂L, that in turn is dominated by the term −L∂V 0 /∂L ∼ V 0 /L ∼ log N . Therefore, B ∼ log N and the compressibility κ T vanishes as 1/ log N . The coefficient of thermal expansion α T can be evaluated from the knowledge of κ T and C a according to [30] where we have used (56). Since the compressibility is dominated by the zero temperature term, the behaviour of the coefficient of thermal expansion α T as a function of T is determined by the heat capacity: Linear dependence at low temperatures, and saturation at high temperatures. 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 . In Fig. 6 the coefficient of thermal expansion for a chain of 1000 ions is presented. Calculations made with different numbers of ions, taking the trap frequency ν such that the linear density at the center of the chain remains constant, show that α T decreases with N . The numerical results are consistent with the behaviour we expect in different temperature regimes, according to the above considerations. In particular, for finite N it is significantly different from zero. This is in contrast to the behaviour found for uniform harmonic solids, where the coefficient of thermal expansion vanishes [30]. The thermodynamic quantities of a linear string of charges confined by a harmonic trap are affected by the way the thermodynamic limit is taken. Nevertheless, the behaviour of the system is intrinsically non-extensive. The non-extensitivity is due to the strong correlation and the dimensionality of the crystal, which determine a regime where the correlation energy, associated with the discreteness of the individual charges, cannot be neglected [20]. It manifests in particular in the log N terms appearing in the thermodynamic quantities. A representative example is the dependence of the specific heat per particle on N .
We finally remark on the thermal stability of the chain. The derivation presented in this section, in fact, relies on the assumption that the thermal excitations do not affect the stability of the system and thus the validity of the physical model we are considering, namely of ions oscillating around their equilibrium positions. This condition is equivalent to the statement that the thermal energy at the considered temperatures is considerably lower than the equilibrium energy, and the displacements are much smaller than the respective interparticle distance. Hence, this condition is valid when This relation is amply satisfied for the parameters of the numerical evaluation of the thermodynamic quantities considered in this section.

V. SUMMARY AND CONCLUSIONS
We have investigated the dynamics of a linear crystal excitations with Coulomb interaction and in presence of an external potential. We have derived an analytical formula for the density of states at long-and shortwavelengths. In the long-wavelength part of the spectrum, we have calculated analytically the eigenmodes and eigenfrequencies. The eigenmodes and eigenfrequencies for the center-of-mass and the first excitation of the axial and of the transverse motion are exact and independent of the number of ions. Apart for these modes, the results we derive are valid in the limit of infinite number of ions. Nevertheless, they already give a good description of the spectrum of excitations of chains of dozens of ions. Using our results we study the statistical mechanics and thermodynamics of the linear chain.
Our derivation allows to find an analytical formula for the critical transverse frequency required for determining the stability of the linear chain. It agrees with the analytical estimate by [24], that was obtained under different assumptions, and it is consistent with the formula fitted from numerical data [22] and verified experimentally [23]. It was suggested that this instability of the ion chain, resulting in a transition from a linear to a zig-zag equilibrium configuration, can be treated as a phase transition [22,23]. In future works we will explore, using the formulation developed in this work, whether the thermodynamic quantities exhibit singularities characteristic of phase transitions [35]. This system differs from systems that are traditionally studied in the framework of statistical physics, since it is not extensive. In particular the specific heat per particle depends on the number of ions.
The results presented in this work show a statistical mechanics approach applied to a strongly-correlated mesoscopic system. They contribute to the on-going research on low-dimensional cold gases [37,38] and may be relevant to studies of the quantum dynamics of manyparticle Coulomb systems like ion crystals in storage rings [9] and cold neutral plasmas [39,40], which at sufficiently low temperatures are predicted to crystallize [41]. The connection with Wigner crystals, where the quantum statistics may play a relevant role [42], will be explored.
Moreover, the spectra of excitations here evaluated are relevant for the implementations of quantum logic with ion traps, where the knowledge of the long-wavelength modes is important for the realization of logic gates [16,17,18,19], as well as realization of solid-state models [43].

Acknowledgments
It is our great pleasure to thank Andreas Buchleitner for hospitality at the Max-Planck-Institut für Komplexe Systeme in Dresden, where most research was done. S.F. is grateful to Wolfgang Schleich and the department of quantum physics in Ulm for hospitality. G.M. acknowledges the kind hospitality of the Institute for Theoretical Physics at Technion in Haifa. The authors acknowledge stimulating discussions and helpful comments of T. Antonsen, M. Drewsen, J. Eschner, D. Habs, P. Kienle, E. Ott, R. Prange, T. Schätz, J. Schiffer, E. Shimshoni, and H. Walther. We thank the referee for bringing to our attention the article in [29].
This work was partly supported by the US-Israel Binational Science Foundation (BSF), by the Minerva Center of Nonlinear Physics of Complex Systems, by the Institute of Theoretical Physics at the Technion, and by the EU-networks QUEST and QGATES.
Using the matrixΛ, in the limit N ≫ 1 the spectrum of characteristic frequencies M (ω 2 ), that is defined as the proportion of roots withω 2 j <ω 2 , and the corresponding density D(ω 2 ) = dM (µ)/dµ| µ=ω 2 , can be derived from the characteristic function Ω(θ), which is defined as [34] where θ is a complex variable. The density of characteristic frequencies D(µ) and the spectrum M (µ) are found by using the properties of the analytic continuation of Ω(θ), In this appendix we discuss an ansatz based on the phononic solution for a translationally invariant crystal, where we assume that the envelope, superimposed on the amplitude and the phase of the phononic solution, varies slowly as a function of the position. This ansatz is supposed to be valid for long-wavelength excitations and chains with N ≫ 1 atoms. The exact solution of (6-8) in the limit of uniform charge distribution with constant spacing a, K i,j /m = κ/|i − j| 3 with κ = 2Q 2 /a 3 , is [30] q j ∼ e i(kja−ωt) (C1) Here, ω is the eigenfrequency, k the wavevector, and the boundary conditions for the uniform chain are fixed so that q N = q 0 , giving e ikN a = 1. Hence, ka = 2πl/N with l = 0, 1, . . . , N − 1. The dispersion relation is obtained by substituting (C1) into (6) with K i,j /m = κ/|i − j| 3 . If the interaction is determined by the nearest-neighbours, while the interaction with the other ions is neglected, then the eigenfrequencies are given by Note that ω depends on ka, which takes the values ka = 2πl/N .
We now construct from this solution an ansatz for the inhomogeneous chain, in the limit of slowly-varying spring constant K j,j+1 over the wavelength of the propagating perturbation. We assume thus the local solution: where ka is a constant, i.e. we have assumed k j x (0) j = jka. The non-linear variation of the phase with the site is included in A j . The ansatz (C3) assumes that it is possible to write the solution as the product of a slowly varying amplitude and phase, represented by the factor A j , and a fast oscillating part, that is the exponential term in (C1). The validity of this assumption is checked later on.
By substituting (C3) into (6), keeping only the nearest neighbour interaction, we obtain: where κ j = K j,j+1 /m. We make use of the fact that κ j varies slowly with the position, which allows to expand A j and κ j in δj, of order unity, according to A j±1 ≈ A ± δA (C5) κ j = κ + ∆κ + δκ j /2 (C6) κ j−1 = κ + ∆κ − δκ j /2 (C7) where δA = ∂A ∂j δj with δj = ±1. The spring constant has been divided into three contributions: κ is constant and fulfills (C2), ∆κ is also constant and δκ j = κ ′ δj, with κ ′ = ∂κj ∂j , such that κ + ∆κ ≫ δκ j . The chosen variation of κ j reflects the symmetry under reflection of the crystal, such that the ion at the center is characterized by equal couplings on both sides (which implies the condition δκ = 0). The term ∆κ does not depend on j, and it is the correction to κ, the spring constant fulfilling the dispersion relation for the case of a uniform chain: ∆κ ≈ κ j − κ. Substituting (C5-C7) into (C4) yields 4∆κA sin 2 ka 2 − 2iκδA sin ka − 2iAδκ sin ka = 0 where we have applied the dispersion relation (C2). Since there is only one zero-order term in j, it is evident that ∆κ = 0: Thus, there is no zero-order correction to the spring constant κ of the uniform chain, and its value in this order is accounted for by (C2). The first order expansion gives thus a differential equation (valid for sin ka = 0) which admits the solution The eigenmode has then the form: and the wave vector takes the same values as for a uniform crystal with half-periodic boundary conditions, taking into account the symmetry under reflection, where k = k(ω) is given by Eq (C2), and l = 0, 1, . . . , N − 1. Therefore, this ansatz gives simply a variation of the amplitude of the displacement of the ion as a function of the local spring constant, but no change in the phase, which is the same as the one of a uniform crystal. The amplitude is smaller at the center of the chain, where the ions are closer and the coupling constant is larger. Taking κ 0 = Λ 0 as given by (42), i.e. the maximum value of the spring constant in the chain, and using (C2), we obtain the maximal frequency as given in (43). A similar argument can be developed for the transverse frequency. Figure 7 compares the prediction of the slowly-varying ansatz with the numerical results and the Jacobi polynomials results given by Eq. (27). Here, it is obvious that the Jacobi polynomials provide a better approximation for the spectrum at the lowest eigenfrequencies. The curve evaluated from Eq. (C2) using Eqs. (C3) and (C11) lies close to the spectrum evaluated numerically for a larger range of modes, but it does not reproduce its asymptotic behaviour for N → ∞.