Finite temperature spectrum at the symmetry-breaking linear-zigzag transition

We investigate the normal mode spectrum of a trapped ion chain at the symmetry-breaking linear to zigzag transition and at finite temperatures. For this purpose we modulate the amplitude of the Doppler cooling laser in order to excite and measure mode oscillations. The expected mode softening at the critical point, a signature of the second-order transition, is not observed. Numerical simulations show that this is mainly due to the finite temperature of the chain. Inspection of the trajectories suggest that the thermal shifts of the normal-mode spectrum can be understood by the ions collectively jumping between the two ground state configurations of the symmetry broken phase. We develop an effective analytical model, which allows us to reproduce the low-frequency spectrum as a function of the temperature and close to the transition point. In this model the frequency shift of the soft mode is due to the anharmonic coupling with the high frequency modes of the spectrum, acting as an averaged effective thermal environment. Our study could prove important for implementing ground-state laser cooling close to the critical point.

Amongst other things, the Kibble-Zurek mechanism [24][25][26] and creation of topological defects have been demonstrated [27][28][29][30]. Two widely discussed transitions are the linear to zigzag [31][32][33] and the pinning to sliding (Aubry) transition [34][35][36]. These were shown to be second-order phase transitions [33,37], that exhibit a soft mode with vanishing frequency at a critical point. The system however is critical solely at zero temperature. Therefore, the observation of critical behavior requires one to characterize and understand finite temperature effects at the transition.
The frequency spectrum at zero temperature is well described by the harmonic crystal approximation. Deviations to this analytical solution and in particular finite frequencies close to the critical point have been observed for the soft mode of the Aubry-type transition in trapped ion chains, at a temperature of around 1mK [12]. Here, we focus on the experimentally more accessible linear to zigzag transition and investigate the coupling of the soft mode to the thermal phonon environment. We * Corresponding author:Tanja.Mehlstaeubler@ptb.de develop a theoretical model that allows one to reproduce the presented spectroscopic measurements by means of a harmonic chain, whose normal-mode spectrum at low frequencies results from the temperature-mediated coupling with vibrational modes at high frequencies. In this sense, the high-frequency modes can be considered a thermal phonon environment. We discuss this result in connection to earlier works [38,39], that described finite temperature effects in terms of an effective shift of the transition point. Our findings deepen the understanding of the complex dynamics of normal modes. They could prove important, for instance, for laser cooling the linear ion chain to the ground state in the vicinity of the transition.
This paper is organized as follows: In Section II we briefly review the linear to zigzag transition. In Section III we present our experimental methods and results of vibrational mode measurements, using resonant light force modulation. Subsequently, in Section IV we compare our findings to molecular dynamics simulations. In Section V, we discuss a simplified analytical model which allows one to gain insight into the temperature dependence of the spectroscopic measurements. In Section VI the conclusions are drawn. The appendices provide supplementary material to the studies presented in Sec. IV and Sec. V.

II. THE LINEAR-ZIGZAG TRANSITION
We consider N ions with charge e and mass m, which are confined by a linear Paul trap. The trap potential is described in ponderomotive approximation by three trapping frequencies ω z , ω x and ω y . The total potential energy V is the sum of the trap confinement and of the unscreened Coulomb interaction between the ions: where r i = (x i , y i , z i ) T denotes the position of the ion i (i = 1, . . . , N ) and ε 0 is the vacuum permittivity. For later convenience, we introduce the vector u = (x 1 , x 2 , . . . , x N , y 1 , y 2 , . . . , y N , z 1 , z 2 , . . . , z N ) T , which gives the configuration of the crystal.
At sufficiently low temperature, the ions localize at the equilibrium positions u(0) of the potential V , for which the equations ∂V /∂u j = 0 holds. In this configuration the dynamics of the chain is characterized by the matrix K ′ , with elements For stable equilibrium K ′ has finite and positive eigenvalues. In the rest of this paper we choose ω z < ω x < ω y , focusing particularly on the aspect ratio α = ω x /ω z for which the ions can either form an array along the z axis, the linear chain, or form a double array in the form of a zigzag configuration on the x − z plane [40]. Figure  1(a) displays an experimental photo of a linear and of a zigzag chain of 30 ions. The two structures are separated by the critical value of the aspect ratio α c ≈ 12. A numerical estimated of the scaling of the transition point with the number of ions gives 0.556N 0.915 [32,41,42]. The corrections due to quantum fluctuations have been determined in Refs. [43,44]. The linear to zigzag instability is a continuous phase transition in the thermodynamic limit, corresponding to letting N → ∞ while appropriately rescaling α [45]. It is associated with breaking of reflection symmetry about the z axis (for ω x = ω y the broken symmetry is rotational and the transition is characterized by a Goldstone mode) [33]. These properties are strictly valid in the limit of T = 0, while at finite temperature the transition becomes a crossover. In order to study and understand the effect of temperatures on the vibrational spectrum across the linear to zigzag transition and in a finite chain, we first discuss the normal mode spectrum at T = 0. The normal mode spectrum is determined by assuming that the thermal energy is much smaller than the potential energy, so that the ion array can be described by a harmonic crystal. The normal mode frequencies are related to the eigenvalues λ j of the matrix K ′ by the relation ω j = λ j /m (j = 1, . . . , 3N ). The corresponding mode vectors are given by the columns of the dynamical matrix λ ij , that diagonalizes K ′ and the mode amplitudes are denoted as Θ j . We use the notation (n, m) to identify the mode vectors by the number of nodal points (phase flips between ions) along the axial (n) and transverse direction (m) 1 . The lowest normal mode frequencies for N = 30 are displayed in Fig. 1(b) as a function of the aspect ratio α across the linear to zigzag transition.
At the transition point the frequency of one normal mode vanishes. In the linear phase, this mode is the zigzag mode and has a purely transverse oscillation with (0,N -1) nodal points, see Fig. 1(c). In the thermodynamic limit the zigzag mode of the linear chain is the soft mode of the phase transition [33]. When crossing α c from the linear (α > α c ) to the zigzag (α < α c ) phase, the soft mode gains an axial nodal point and becomes the new breathing mode (1,N -1) of the zigzag configuration. The axial breathing mode of the linear phase (1,0), increases in frequency at α c and gains 2 axial nodes, as well as a transverse zigzag pattern to become the (3,N -3) mode.
In the rest of this paper we analyze how the normal mode spectroscopy at the structural transition is modified at finite temperatures.

III. MEASUREMENT OF VIBRATIONAL MODES
In this Section we describe our experimental method for measuring vibrational mode frequencies, that makes use of a single laser beam with frequency near resonant to the Doppler cooling transition. We then present and discuss our measurements of the lowest axial modes near the linear to zigzag transition. This method was originally introduced in [12].

A. Setup
In order to have a well defined ordered structure, we trap N = 30 172 Yb + ions in a linear Paul trap with high control of electrical fields [46,47]. All ions are illuminated by a linearly polarized laser with central wavelength of 369.5 nm addressing the 2 S 1 2 ↔ 2 P 1 2 transition in Yb + and cooling the ions close to the Doppler cooling temperature of T D = 0.5 mK. As shown in Fig. 2, the laser beam forms an angle of about θ = 25°with the axial direction of the ion crystal and an angle of about ϕ = 45°with the transverse direction. We denote its wave vector as k 1 , which we use in the subsequent text to identify the laser beam itself. The beam has elliptic shape with waists of approximately 2.6 mm in the horizontal and 80 µm in the vertical direction, resulting in an almost uniform illumination of a 400 µm times 20 µm ion crystal in the z-x plane. Typical laser powers in the subsequent measurements are P 1 = 1 mW, corresponding to a saturation of s 1 ≈ 1.75 at the beam center. We also employ a second laser beam at the same wavelength and with the same angles to the crystal, that is focused to a beam waist of about 80 µm in both vertical and horizontal direction, addressing a smaller region of the crystal. Its wave vector is denoted by k 2 . The beam is amplitude modulated in order to excite the crystal's normal modes. The amplitude modulation is added by applying a sine wave with frequency ω e to the RF amplitude of an acousto-optic modulator used as a fast shutter. The modulation of the power is given by where P m is the maximum power in the beam. The saturation of k 2 at the beam center is then where s m = Pm/2 P2,s with saturation power P 2,s of k 2 . The total saturation of an ion at the center of the beams is than s = s 1 + s m + s m cos(ω e t). The ions fluorescence is imaged via a lens system of N/A = 0.2 and recorded by an electron-multiplying (EM)CCD camera, which can resolve individual ions.

B. Method
We excite the crystal's collective motion with the help of the amplitude modulated cooling laser. Both cooling lasers k 1 and k 2 exert a constant light force on the crystal, that shifts the minimum of the trap potential. The amplitude modulated laser adds an oscillating force F m with excitation frequency ω e . This oscillating force is roughly linear, if s m is smaller than s 1 . For multiple ions, the saturation s 2 and s 1 will depend on the ion positions with respect to the laser beam center. Specifically, the saturation power P s (r i ) will depend on the position of the ith ion.
In principle, all normal modes can be excited by means of this technique. In the measurements we present below the waist of laser k 2 was focused to only 80 µm. It illuminates several ions at the same time, as illustrated in Fig. 2(a). This prevented the excitation of modes with a higher number of nodal points, due to the small overlap of their mode vector with the laser intensity profile. We note that the excitation of an arbitrarily chosen modes can be realized by implementing single ion addressing.
On resonance, the amplitude of the driven mode increases linearly with F m [48]. In combination with a constant linear damping γ due to laser cooling of k 1 , a steady state with a constant, frequency-dependent mode amplitude Θ j (ω e ) can be reached after several oscillations. In order to detect an excitation, we record the ions fluorescence with an EMCCD over an exposure time of typically 100 ms. This is long compared to the normal mode oscillation periods, which are on the order of the center of mass oscillation period of about 40 µs. Therefore, light from all possible ion positions during the oscillations is recorded, leading to an apparent increase of the ions' size at the resonance ω e ≈ ω j . The imaged spatial extent of each ion i is proportional to the amplitude of the driven normal mode and the ions vector element of the corresponding mode vector λ ij . The resonance frequency is found by identifying the frequency at the maximum amplitude of the ions oscillation. An experimental photo of the excited (1,0) mode and (2,0) mode in the linear phase is shown in Fig. 2(b), for which we used P 0 ≈ 20 µW (s m ≈ 0.53 for an ion at the beam center). Similar to the ion amplitude, the velocity of the ions increases on resonance, leading to a drop in fluorescence due to the Doppler shift. For a single ion this decrease in fluorescence can be measured with a PMT and enables one to identify the motional resonance.

C. Experimental Results
We measured several mode frequencies near the linear to zigzag transition in two experiment series, with two different modulation powers: P For measurement run (A), we determined the center frequency ω j of any resonance by scanning the excitation frequency ω e manually and searching for the maximum amplitude of the ions for the excited mode. The uncertainties were estimated by finding a region, in which the amplitude of the excitation was still maximal. This resulted in a typical error of about 100 Hz to 300 Hz for each resonance. The power of the excitation laser was set to P (A) m = 20 µW, chosen such that a resonance of the three to four lowest modes could be observed.
At first, the trapping ratio α is determined by measuring the axial and transverse center-of-mass mode frequencies, i.e. the trapping frequencies. This is followed by searching for the low lying modes with 1,2, and 3 axial nodal points. In Fig. 3 we show the measured vibrational mode frequencies, in comparison to the normal mode frequencies expected from the second-order approximation. Away from the phase transition the experimental results agree with the theoretical predictions. However, close to the phase transition, the frequency of the zigzag mode does not vanish. The measured frequency of mode (1,N -1) (blue circles) increases when α approaches α c until it reaches the expected frequency of the breathing mode of the linear phase. While the (2,0) mode frequency (purple diamonds) was observed over the complete phase transition region, close to transition the breathing mode in the zigzag phase (red diamonds) was not detected.
In the second measurement run, a single ions fluorescence was recorded with a region of interest (ROI) on the EMCCD, while sweeping the excitation frequency. Near resonance, a decrease in fluorescence in the ROI is observed, because the excited ion moves partially out of the ROI during exposure and it gains a Doppler shift due to its increased velocity. We fit the fluorescence drop to a Lorentzian line shape in order to determine the resonance frequencies of the axial center of mass and the normal mode with one axial nodal line.
To obtain a finer resolution in series (B) the maximum power of the amplitude modulated laser was about P (B) m = 6 µW, sufficient to excite the center of mass and breathing mode. In Fig. 3, we show the results of these measurements as cyan diamonds. The results agree qualitatively with the measurements from series (A). We verify an increased vibrational mode frequency of the zigzag mode in the zigzag phase close to α c .
The quantitative difference between (A) and (B) is due to the smaller power of the amplitude modulated laser used for measurement run (B). A larger driving force increases the mode amplitudes in the steady state and therefore enhances non-linear frequency shifts due to the Coulomb interaction. Increased power leads to a higher observed frequency at maximum excitation. Additionally, the lineshape becomes increasingly asymmetric with increased power. We refer the interested reader to Appendix C, where we discuss he influence of P 0 on the frequency of the (1,N -1) mode, see

IV. MOLECULAR DYNAMICS SIMULATIONS
As seen in the last section, close to the phase transition the measured excitation frequencies significantly deviate from the vibrational spectrum in the harmonic approximation. This can only be explained by anharmonicities. In the experiment there are two excitation sources: the thermal noise from laser cooling and the sinusoidal driving force. In order to gain deeper insight into the impact of the mode populations on the measurable frequencies, we carry out molecular dynamics simulations of the crystal under a stochastic force.
We simulate the dynamics of the ion crystal by numerically solving the classical equations of motion of the ions in the presence of damping and of the Langevin force describing thermal noise [28]. This approach is complementary to the Fokker-Planck equation for Doppler cooling of an ion crystal [49]. The equation of motion for the ith degree of freedom takes the form where γ is a damping term from laser cooling and V is the total potential energy, see Eq. (1), and ξ i (t) is the stochastic force, with moments: Here, . . . indicates ensemble averaging. The second equation links the amplitude of the stochastic force with temperature and damping coefficients and is commonly known as the fluctuation-dissipation theorem. In the simulations, laser cooling is treated as isotropic for all degrees of freedom. We first detail the simulation procedure and the spectral analysis. The results of the molecular dynamics simulation are then reported and discussed in Section IV C.

A. Simulation procedure
The ground state configuration for each trapping ratio α is found by simulating a crystal with N = 30 ions choosing strong damping mγ = 2.5 × 10 −19 kg s −1 and with T = 0. The resulting equilibrium positions so obtained are the initial configuration for the simulation at finite temperature T . The simulation is run with lower damping mγ = 2.5 × 10 −21 kg s −1 , comparable to experimental conditions [28], and over a time of 100 µs to thermalize the system. This result is used as a starting point for the final simulations, which run for 10 ms in total. All simulations have an integration time step of 19 ns, which is much smaller than the expected period of the vibrational mode with the largest frequency, which is here the transverse vibration of the center of mass mode at about 3 µs. Every 100th value is saved, resulting in time resolution of the ions evolution of 1.9 µs.

B. Spectral analysis
In order to extract the normal mode spectrum of the ion crystal, we carry out the a Fourier transform (FT) of the trajectories of the ions' axial and transverse degrees of freedom. Due to the simulation length and time resolution the FT has a frequency resolution of 100 Hz and a maximum observable frequency of about 263 kHz, which covers the frequency range of interest. For our analysis, the simulation procedure described above is repeated five times, due to the stochastic nature of the thermal noise and the FTs are averaged over all simulations with identical parameters. Then the absolute value A(ω) χ,i = F (x χ,i ) of the averaged FTsF is calculated, where i is the ion index and χ is either x or z. We are interested in collective motions of the crystal, i.e. the normal modes, but we do not make any assumptions on possible mode vectors. Therefore, the A(ω χ,i ) for all ions are added together to get a signal S(ω) χ : where the degrees of freedom along x and z are treated separately. An example of such a signal is shown in Fig. 4. The width of the resonances depends on the damping γ, which is here fixed to the value mγ = 2.5 × 10 −21 kg s −1 in order to compare the data to our experiment.
We extract the resonances and the peaks' widths from S(ω). Without prior knowledge of the complete model of the peak functions, we estimate the positions based on a peak search algorithm. This method cannot treat noisy signals well. Therefore, S is smoothed before starting a peak search, using a running mean over n values S rm (k) = n i=0 S k+i /n. For T = 0.1 mK an average over 20 points (≈ 2 kHz) and for the other shown temperatures an average over 30 points (≈ 3 kHz) is used. The running mean of the signal is plotted as a function of the respective running mean frequency f rm (k) = n i=0 f k+i /n. In Fig. 4 the smoothed signal S rm is shown as a red dotted line. The orange line shows the peak positions and estimated half-maximum widths which we identified. Not all peaks are captured by the algorithm, especially for ω > 2π · 130 kHz, but this is outside the frequency range that we are interested in. Due to the noisy data, additional small peaks might be found close to strong resonances, e.g. the axial center of mass mode. These false positives have to be removed manually.

C. Numerical Results
We evaluate the simulations for temperatures T = [0.1, 0.5, 2.0, 3.5] mK and for several trapping ratios α. Figure 5 displays the estimated peak positions. The expected normal mode frequencies at T = 0 is also shown for comparison. We observe good agreement between the numerical results and the harmonic spectrum at low temperatures. For higher temperatures, deviations appear near the phase transition for the (1,0) and (0,N -1) modes, that qualitatively agree with the experimental measurement. The simulations for T = 3.5 mK match best to the measurement series (B), which are shown in Fig. 5(d) for comparison. Here, we point out, that the frequency of the zigzag mode of the linear chain remains finite at α c and increases with the temperature.
In the simulations we do not observe a significant deviation of the (2,0) mode from the harmonic approximation, as we did in the experiments. This is most likely due to the high excitation power P (A) m used in series (A), as we described in Section III C. The additional increase in mode amplitude leads to non-linear mode coupling on top of temperature effects.
The simulations reveal that in the zigzag phase thermal effects give rise to collective jumps of the ions between the two degenerate zigzag configurations. We illustrate the mechanism as thermal switching between the two minima of the Landau free-energy in the symmetry broken phase [33], see Fig. 6 (a). We can estimate the corresponding switching rate k est by counting the number of sign changes of the transverse coordinate of the central ion P x N/2 over the simulation length ∆t: k est = P x N/2 /∆t .
The inverse of this rate is the average dwelling time in one crystal configuration τ e = k −1 est . We identify two regimes with the help of τ e . In the first regime τ e is the largest time scale of the system, in particular it exceeds the characteristic period of oscillations of the (1,N -1) mode, denoted as T zz . Here the two crystalline configurations are well defined and thermal effects give rise to approximately instantaneous jumps, whose net effect is to broaden the linewidth of the resonance lines. In the second regime, where τ e T zz , the non-linearities of the system are dominant and expected to modify the normal mode spectrum. Figure 6 (b) displays τ e as a function of the aspect ratio and for different temperatures, the horizontal line indicates τ e = T zz . We have verified that frequency deviations from the harmonic solution are observed when τ e T zz . This is visible, for instance, in In order to gain insight into these dynamics, we analytically estimate the switching rate using the Landau free energy of Ref. [33]. In the symmetry broken phase the free-energy is a double well potential, the two minima correspond to to the two zigzag configurations. We then interpret the switching rate as the rate of thermal activation overcoming the barrier that separates the two minima [50], as shown in Fig. 6(a). For this purpose, we first calculate the potential energy along the adiabatic path connecting the two equilibrium configurations [51]. Since the two ground states stem from the breaking of the mirror symmetry, we parameterize the path by the transverse crystal size: where x N/2 is the transverse position of the ion left of the crystal center for even N . For the calculation we minimize the crystal energy using a Lagrange multiplier with a constraint for the crystal transverse size g(u) = X. Afterwards the total potential energy for this configuration is taken as the energy U (X) of the potential at size X. In the 2D phase it has the shape of a double well with two minima, symmetric about X = 0. The energy barrier E B is then given as the difference between the potential energy at X = 0 and the minimum potential energy Sufficiently close to the transition the energy barrier increases with |α − α c | 2 , see inset in Fig. 6 (b), in agreement with the predictions of Ref. [33] and with the numerical simulations of clusters of metallic beads [52].
The trajectory of the collective coordinate of the crystal that jumps between the minima of the bistable potential results from the interplay of driving, damping, and noise. Quantitatively accounting for the prefactors in the Kramer's escape formula [53] is beyond the scope of the current work. Here, we perform an estimate using transition-state theory [50]: with ω a = U ′′ (X min )/m, where U ′′ is the second derivative with respect to X and X min is the transverse crystal size in equilibrium. For α = 11.8 and T = 2.0 mK the transition-state theory predicts a rate about 16 000 s −1 . Taking into account, that in the simulations the particles can also return to each minima the escape rate in the simulation from one minimum is about 14 000 s −1 . In Appendix B we compare the rates extracted from the simulations with the predictions of transition-state theory rates over a range of parameters. We find agreement within a factor of 2.
The molecular dynamics simulation validate the thermal fluctuations as the source of the observed frequency deviations. Moreover, they supply a deeper insight into the exact dynamics behind the non-linear mechanism at hand, i.e. the frequent crossing of the potential barrier between the two degenerate ground states of the zigzag phase.

V. EFFECTIVE MODEL FOR THE MODES DYNAMIC
In this Section we use a simplified model in order to determine the temperature dependence of the mode spectrum close to the linear to zigzag instability. For this purpose we use a complementary approach to the one based on dwelling times and consider the normal mode expansion around the linear chain for aspect ratios at which the linear chain is mechanically unstable. We show that the effect of anharmonicities at finite temperatures is to effectively stabilize the linear chain. The resulting normal mode spectrum agrees with the numerical results close to the transition point, as we discuss below and summarize in Fig. 7.

Normal modes at the instability
We first review the normal modes of the linear chain and the equations for the structural instability in the absence of damping and noise. Close to the linear to zigzag instability we expand the total potential energy V of Eq. (1) to fourth order around the equilibrium positions of the linear chain, where q i are the displacements around the equilibrium positions u i (0), q i = u i − u i (0), and the tensors K ′ , L ′ and M ′ are given by the expressions: Note that V 4 approximates the total potential V , Eq. (1), in the limit in which the displacements around the equilibrium positions are much smaller than the interparticle distance at equilibrium. Let Θ j denote the normal modes of the linear chain, which diagonalize the matrix K ′ and have eigenvalues mω 2 j . The linear chain is stable provided that all eigenvalues are positive. In this regime the ω j are real and correspond to the normal mode frequencies. The condition min j ω j = 0 identifies the classical transition point of the linear to zigzag instability. Potential (8) is cast in terms of the normal mode by means of the dynamical matrix λ ij such that q i = j λ ij Θ j and takes the form where now the tensors L and M are related to the tensors L ′ and M ′ by the relations: The total Lagrangian for the normal modes takes the . In an appropriately defined thermodynamic limit the potential V 4 can be cast in terms of a Landau potential at the linear to zigzag instability [33]. The one dimensional model strictly exhibits a phase transition at zero temperature, where a quantum description becomes appropriate [43,44,54,55]. In what follows, instead, we consider a finite system and do not scale the physical parameters with N . In our analysis, therefore, the threshold of the instability as well as the gap between excitations depend on the system size.

Thermal effects
We now discuss the low frequency spectrum of the linear chain across the linear to zigzag instability and in the presence of laser cooling. We consider a finite chain and, starting from the Fokker-Planck equation [49], we model the dynamics of laser Doppler cooling in terms of Langevin equations. We denote the damping (cooling) rates of the normal modes by γ i and write the corresponding Langevin equations as [56] where Ξ i (t) is the Langevin force for the normal mode Θ i , and we neglect here mode-mode correlations due to the dissipative dynamics.
For finite chains and in the zigzag phase the lowest frequency mode is a superposition of the zigzag mode and of the axial breathing mode of the linear chain as seen in Fig. 1c). The gap between the mode driving the instability and all other normal modes remains finite. Thus, whenever thermal excitations and the line broadening are smaller than the gap, normal-mode spectroscopy of the chain shall provide in first approximation the mode spectrum obtained by diagonalizing the quadratic term of potential (1) about the stable equilibrium positions.
We now determine the effects of thermal excitation on the lowest energy spectrum by considering the equations of the lowest energy mode, here (0,N -1), which we label by Θ 1 , and the mode which is closest in frequency and to which it couples. This mode is labeled by Θ 2 and is according to our notation (1,0). We then make the simplifying assumption that ω 1 , ω 2 ≪ ω ℓ , where ω ℓ are here the frequencies of the modes Θ ℓ to which Θ 1 and Θ 2 appreciably couple through the anharmonicities. In this regime we can identify the time scale δt for which ω 1 δt, ω 2 δt ≪ 1 and ω ℓ δt ≫ 1. Moreover, we assume that the modes Θ ℓ are at thermal equilibrium. We now perform the time average of Eq. (15) over the grid with step δt. Under this assumption anharmonic terms in Eq. (15) become, e.g., Θ ℓ Θ ℓ ′ Θ 1 → δ ℓ,ℓ ′ Θ 2 ℓ δt Θ 1 , where · δt denotes the time average. This procedure leads to the two coupled equations: whereω i , ν 12 and η i are explicitly dependent on the temperature. In particular, the frequency squaredω 2 i now readsω and it contains a shift proportional to the temperature with proportionality constant 11.6 11.7 11.8 11. The second and third terms on the right-hand side of Eqs. (16)- (17) describe an effective coupling between the two modes and a mean displacement force, respectively, with The effective constants ν 2 eff,i , ν 2 eff,12 and η eff,i are determined by carrying out the summation over the other modes.
Equations (16) and (17) describe mode mixing and frequency shifts induced by the thermal excitation of the chain. Within this classical model these terms are directly proportional to the temperature. We can now determine the resulting normal mode frequencies. For this purpose we note that the term ν 12 = 0 for the expansion about the linear chain equilibrium positions, see Table  I of Appendix A. This is a consequence of the fact the breathing mode is an exact eigenmode of the linear chain [45]. In the underdamped limit the characteristic frequencies are now given by Eq. (18) and Eq. (19). Figure  7 displays the frequenciesω i as a function of the aspect ratio α and for four increasing values of the temperature, ranging between 0.1 mK and 3.5 mK. For comparison, the results of the numerical simulation of Eq. (3) are reported. The prediction of the analytical model and the result of the numerical simulation agree for aspect ratios close to the transition point α c : this is the regime where our model is plausible since the truncation of the Taylor expansion is justified. Loosely speaking, the thermal effects stabilize the linear chain also for aspect ratios beyond the critical point. This behavior might be interpreted as a shift of the transition point [38]. However, it is actually the manifestation that at finite temperatures the linear to zigzag transition is a crossover region [52,55,57].

VI. CONCLUSION
In this work we investigated experimentally and theoretically the effect of thermal noise on the low-frequency spectrum of an ion chain near the symmetry-breaking linear to zigzag transition.
In the experiment we employed resonant light force modulation with an amplitude modulated laser beam to excite collective oscillations in a crystal. The method is simple to implement and can also be used to measure trapping frequencies, replacing established excitation methods, such as modulation of the trapping potentials [58]. This allows for stronger filters in the RF and DC electronics of the Paul trap [46,47], reducing the heating by electrical noise of the trapped ion crystals. While we used an excitation beam profile encompassing multiple ions, a more focused beam or a spatially engineered beam profile, e.g. generated by an spatial light modulator [59,60], would allow for arbitrary mode excitations.
The experimental measurements did not show the softening of the zigzag mode that is predicted at the structural phase transition. Also the frequency of the breathing mode was nearly constant and independent of α when sweeping into the zigzag regime, instead of increasing as expected in the absence of thermal noise.
With the help of molecular dynamics simulations we could reproduce the experimental observations within the uncertainties, thus confirming that this behavior is majorly due to thermal excitations. In particular, inspection of the trajectories show that finite temperature effects induce collective jumps of the ions between the two degenerate zigzag configurations. This microscopic picture is at the basis of the expected crossover behavior at finite temperatures.
We developed a simple analytical model that builds on these findings and predicts the experimentally observed frequency spectrum. This model shows that the shift of the zigzag mode at the transition point is due to anharmonic coupling with high frequency modes, which act as an effective phonon environment. Separation of timescales between the low frequency soft mode and the higher frequency modes allows taking the averaged higher frequency modes as an effective potential that influences the soft mode. Note, that the thermally excited phonon environment in our model can be replaced by non-thermal excitations.
Our results suggest that a similar model can be developed in order to describe the experimental measurements of the Aubry-type transition in ion Coulomb crystals.
Here the soft mode of the pinning to sliding transition exhibited a finite frequency at the critical point, when a finite temperature allowed the system to switch between different minima of the Peierls-Nabarro potential [12,61].
Our work is relevant for experiments operating close to phase transitions in ion chains such as studies of energy transport [62] or quantum information using the gapped topological defect mode [63,64]. According to the results presented the cooling of high-frequency modes is crucial to avoid the heating of the soft mode due to higher-order coupling, showing the experimental complexity of these plans. Similar challenges were recently discussed for laser cooling a 2D planar crystal confined in a Penning trap [65].

ACKNOWLEDGMENTS
We gratefully thank J. Keller for fruitful discussions. This project has been supported by the Deutsche Appendix B: Switching rates: Simulation and Transition-state Theory In Fig. 8 we show the comparison between the escape rate obtained by transition-state theory from the double well potentials calculated in Section IV C and the rates estimated from the simulation results. The latter have been corrected by the factor two for the comparison to account for the possibility to come back to a potential minimum.

Appendix C: Power dependency of experimental signal
Driving the intrinsically non-linear Coulomb system can lead to additional frequency shifts due to further increased amplitudes. The excitation method employed in this paper can lead to such shifts, depending on the laser power of the modulated laser P 0 .
In Fig. 9 we show the power dependence of recorded resonance features of measurement series (B), when sweeping over the axial center of mass and the breathing mode at α ≈ 11.72. It can be clearly seen, that the position of maximum excitation shifts, when increasing the driving force. Additionally, the resonances become increasingly asymmetric with higher forces. The center of mass mode resonance around 31 kHz increases in width with higher forces.
Measurements presented in Section III C, were carried out with P m = 20 µW for series (A) and with P m = 6 µW for series (B).