Defect generation and dynamics during quenching in finite size homogeneous ion chains

An equally spaced linear chain of ions provides a test-bed for studying the defect formation in a finite size 1D system. In particular, defect formation related to topological phase transition from a linear configuration to a zig-zag one is of interest here. A semi-empirical expression provides an excellent agreement to the numerical results. The non-adiabatic transition between the chain and zig-zag topologies for a finite size system of 30 ions shows clear distinction from non-uniformly distributed ion chain. Thus the underlying Homogeneous Kibble-Zurek model can be tested in presently accessible ion trap experiments. Furthermore, our study indicates collective defect behaviour appearing through the correlation length measurements.


Introduction
Phase transitions occur in a multitude of physical scenarios, from the Bose-Einstein condensation to boiling of water. Moreover, symmetry breaking during second-order phase transitions have been invoked to explain the observed inhomogeneity mass distribution of the universe through the Kibble-Zurek mechanism (KZM) [1,2]. Such a mechanism is universal, making some of the aspects of the dynamics independent of the underlying system. At the heart the KZM resides the concept of causality. The original KZM assumes that the system is initially homogeneous (also known as the Homogeneous KZM), leading to a phase transition occurring simultaneously across the whole system. If the system is crossing to a groundstate with multifold degeneracy, the system has to choose one among of the different possible final states. In a non-zero temperature scenario, some region of the system will cross first making necessary a local choice for the new state. If the speed at which this choice can be communicated,ŝ, is larger than the speed of propagation of the front velocity, v F , then the final state will be simply one of the different possible new ground-states in this adiabatic regime. However, ifŝ < v F , then, there is a probability that different causally disconnected regions will chose different states, leading to defect formation at the separation of the regions. Notice that for perfectly homogeneous systems v F diverges and therefore the probability to generate defects is always greater than zero.
The main measurable prediction of the KZM is the density of topological defects as a function of the rate at which the phase transition is crossed. Such speed is controlled experimentally through the time variation of a control parameter. Three main difficulties typically appears when studying the KZM in the laboratory. First the weak power law scaling of the defect densities requires to experimentally explore a large range of the control parameter. Second, the preparation of an homogenous initial system and finally, the preparation of a large enough system where finite size effects do not play a significant role.
An overview of the theoretical background of the KZM and of the experimental efforts regarding the KZM can be found in [3]. To our knowledge, the best agreement between the KZM and experimental work to date have been achieved using linear ion chains formed in a linear RF traps [4,5]. However, such experiments could not explore the homogeneous case as the ion density in a standard ion chain is not constant [6]. The initial theory was adapted to include this, leading to the Inhomogenous KZM [7] (IKZM), where the transition is crossed always at a particular place first, at the center of the ion chain, leading to a modification of the scaling law.
The HKZM can be studied using laser cooled ions if an homogeneous system could be generated. One possibility is the generation of ion rings in multipole traps [8]. The HKZM using an ion ring has been studied before numerically in [9]. However, the experimental realisation of ion rings in multipole traps has proven more difficult than initially thought [10,11], with only one reported ion ring generated in a multipole surface trap [12].
Another possibility is to use a linear ion chain with uniform inter ion spacing. Such ion crystal configuration has been proposed in the context of Quantum Computing with ion traps [13].
Here, we explore the possibility to use such topology, the homogeneous ion chain, for the study of the non-linear dynamics and HKZM. In particular, we have used molecular dynamics (MD) simulations of the 1D to 2D transitions of an homogeneous ion chain.
The article is organised as follows : first, the potential needed to generate an homogenous ion chain is discussed, followed by the derivation of the critical parameter at which the linear chain to zig-zag transition occurs. A discussion of the MD results in terms of defect dynamics in the context of the HKZM is elaborated.

Homogenous Linear Ion Chain
The first attempt to obtain a potential to generate a truly homogeneous ion chain is due to Johanning et. al. [14], where an analytical solution was found : where the ψ (n) (x) represents the polygamma function, d is the inter-ion spacing, N is the number of ions, Q is the ion's charge and k C is the Coulomb constant.
A more recent work [15] uses a different approach, where the ion chain is approximated as a uniform charge density ρ 0 = Q/d, leading to an analytical expression for the effective electric field acting on each ion due to the Coulomb interaction with the other ions on the chain. By integrating such electric field, a different expression than equation 1 was derived : where L is the half length of the ion chain. However, the above expression leads to infinite potential walls as noted by the authors. Moreover, such potential lead to a 5% variation on their studied case (number of ions, N = 50, and d = 3µm) [15]. In the following, it is shown that using equation 1 leads to a higher degree of homogeneity than equation 2. For such reason, equation 1 has been used as axial potential in the rest of this work.

Molecular Dynamic Simulations
The problem at hand involves solving the equations of motion of the N interacting ions, in the presence of the trap potential. The laser cooling is represented through a friction term while the heating parameters are taken into account through a single thermal bath. Such a problem is described by the following Langevin equation for each ion i evolving in a 3D space : where F c represents the Coulomb force due to the other ions, ω x and ω y are the secular frequencies along the transverse direction, x and y respectively, V z is the axial potential given by 1, Γ is the friction coefficient, k B is Boltzmann's constant, T is the temperature, which has been assumed to be the same on all three spatial dimensions and θ xj and θ zj are a collection of independent standard Wiener processes [16]. The equation of motion are numerical solved using the vGB82 as described in [16].
First, MD simulations is used to test the degree of homogeneity of the ion-ion distance of the final ion distribution when the potential given by eq 1 is used. The parameters used are : ω x /2π = 430kHz, Γ = 1.5 · 10 −20 kg/s and N = 128 ions of 138 Ba + . The temperature is set to the typical Doppler limit value of T = 0.5mK. The ions are initialised as a perfectly homogeneous chain with inter-ion distance d = 10νm, followed by uniformly distributed random displacement between ±5mum is applied in both transverse and axial directions. After 6ms of evolution time (with time step of 2ns), we obtain a final average ion-ion distance of ∆z = 10.00±0.13µm, (the standard deviation is used as error measurement) which represent an excellent agreement with the designed inter-ion distance of 10νm. If we repeat the simulation for a higher temperature, T = 5mK, we obtain ∆z = 10.0 ± 0.4µm.

Linear to Zig-Zag transition
The next step is to obtain the control parameter at which the phase transition occurs. In the present system, it is the value of ω x that is reduced over time while keeping ω y fix and ω y ω x , leading to a decrease in the transverse confinement which eventually generates the topological phase jump to a zig-zag structure along the x direction. The value of ω y /2π = 1MHz is kept constant through the rest of the present work.
If periodic boundary conditions are used, the critical frequency, ω c , at which the transition occurs can be obtained analytically [Fishman 2008] : where γ(p) is the Riemann-zeta function. The periodic boundary conditions do not represent correctly finite size system. In the following, a more realistic situation has been considered, which lead to a better agreement with numerical results as it is shown shortly. Let's assume that we have an homogeneous zig-zag ion configuration, as given by the full circles in figure 1.
We concentrate on the case where N is odd. We assume a perfect zig-zag along the x plane (y = 0), see full circles on figure 1, with a axial inter ion distance, z j+1 −z j = d and with |x j | = h. By symmetry considerations, only the ions marked with a full square contribute to a net transverse force on the central ion as the axial force contribution of the all the other ions cancels out. Therefore, the total force on the central ion can be written as : where N = | N +1 4 |. The force due to the confining transverse potential is given by F trap = −mω 2 h. By imposing equilibrium of forces we obtain : Taking the limit lim h → 0, we arrive at : The MD simulations have been performed in order to check the validity of the above expressions, equations 4 and 8. The ions are initialised in the chain phase, ω a > ω c , with zero velocity. They are first thermalised at constant radial frequency, ω r = ω a , during 2.0ms. Then the ω frequency is lowered linearly to a final frequency in the zig-zag phase, ω b < ω c at a speed, given by the quench rate : ∆ω τ Q , where ∆ω = ω b − ω a and τ Q is the quench duration : Simulations have been performed for N = 30 ions and d = 10µm, using a thermal bath temperature of 1nK. While this temperature is a non-realistic experimental value, it seemed appropriated to use a very low temperature as equation 4 and equation 8 are derived assuming zero temperature. The values used for the initial and final transverse frequencies are ω b /2π = 140kHz and ω a /2π = 500kHz. The critical frequency, obtained from equation 8 is ω c = 327.38kHz.
The maximum transverse distance between two ions, h, as a function of the instant transverse secular frequency, ω, for different values of τ Q is shown in figure 2. The reader is reminded that time runs from higher values of ω to lower ones. Topological phase transitions are clearly observed : for ω > ω c , the value of h is very close to zero (for the slowest quench, the mean value of the h before the transition is < h >= 0.34 ± 0.17nm), which corresponds to a chain configuration, while for ω r < ω c there is a jump on the h value, corresponding to the zig-zag phase. Figure 2 shows a strong dependency between the frequency at which the phase transition occurs and the quench rate. For the fastest quench, ∆ω τ Q = 10 12 [Hz/s], or τ Q = 2.26µs, the ions are still in the chain phase at the end of the quench. The observed oscillations following the phase transitions are quickly reduced when the quench time increases. The critical frequency obtained analytically correspond to the adiabatic case, where ∆ω τ Q → 0. Similar behaviour was observed by Shimizu et al when studying the dynamics of a Mott insulator to a superfluid crossing [17].
A fit on the slowest quench studied ∆ω τ Q = 10 6 [Hz/s] of the type h 2 (ω) = a 0 ω 4 + a 1 ω 3 + a 2 ω 2 + a 3 ω + a 4 on the values corresponding to h > 0,see fig 2c, leads to critical frequency of ω M D c /2π = 326.96kHz, obtained by finding the roots of the fit function. The same value is obtained if the data from the quench ∆ω τ Q = 10 7 [Hz/s] is used instead.
If we repeat the above analysis for different initial values of inter-ion distance (using N = 32, R = 10 7 [Hz/s] and T = 1nK), figure 3 is obtained. The solid line of figure 3 correspond to eq 8, while the circles correspond to the critical frequencies deduced from the numerical simulations. The agreement between them is excellent.
If now we fix the inter-ion distance to d = 10um and modify the number of ions from 5 to 64, we obtain figure 4. In this figure the ω c (d) from eq 4, the computed values for ω c (N, d) from eq 8 and the results of the MD simulations are compared.
The equation 4 fails to reproduce the numerical results as expected. Although closer,equation 8 also fails in predicting the right critical frequencies, particularly for smaller ion numbers. The reason becomes clear by comparing the structure assumed for  the force calculation (full circles) and the structure at the end of the quenches (empty circles), as shown in figure 1. The configuration used for the derivation of equation 8 is an over simplification and therefore the total transverse force experienced by the central ion has extra terms which were not taken into account, leading to a lower value of the critical frequency. Nevertheless, eq 8 gives a hint on the nature of the N dependence. The following function has been use to fit the MD results : leading to (a 1 = (6.8 ± 0.3) · 10 − 5, a 2 = 1.90 ± 0.03 and a 3 = 3.22 ± 0.02). The fitted function shows an excellent agreement with MD simulations results.

Defect generation during quenching in Homogeneous Chains
The focus is shifted now to the generation of defects during a non-adiabatic crossing of the transition. To our knowledge, only two numerical works have been done regarding KZM in homogeneous ion crystals [9,18]. However, both cases used the periodic boundary conditions (the facto, simulating an infinite chain) and therefore the simulated system differs with the present one.
We have performed MD simulations of the of quench from ω a /2π = 500kHz to ω b /2π = 140kHz for different values of the quench time, τ Q and different values of ion number, N . The ions are initialised with a d = 10µm and evolved for 2ms, before the quench starts. All the following results were obtained for the 138 Ba + ion, T = 0.5mK and Γ = 1.5 · 10 −20 kg/s. Two examples of the final configurations for N = 64 are shown in figure 5. The defects can be classified in three categories [19] : intermediate defect, with one ion at the trap axis (|x| ≈ 0 ; odd defect, where two consecutive ions sits beside each other and extended defects where two ions have nearly the same axial position. An algorithm that detects the number of defects by counting how many consecutive pairs have the same x sign is used. Such an approach detects all three types of defects simultaneously. Its implementation is efficient and universal as it does not need to impose any threshold. Unlike ref. [9,18], the present system has edges where defects can be lost. For this reason information of the defect evolution during the quench is needed. This is achieved by measuring the total number of defects on the crystal when the mean absolute transverse displacement, x = |x i |/N reaches a particular fraction, , of the adiabatic one, r = x adiabatic . In this way, the average defect number was measured for ranging from 0.25 to 0.85 in intervals of 0.5. For small values of τ Q , the quench is too fast for the ions to reach the required values of during the quench time itself. In those situation, the MD is continued with ω = ω b until = 0.85 is reached. Figure 6 shows the evolution of the logarithm of the normalized number of defects, ln(n/N ), versus the logarithm of the adimensional quench rate, ln(ν), where ν = 1 τ Q ω0 , with ω 2 0 = Q 2 k C md 3 [9]. Two cases are presented, figure 6a correspond to = 0.25 while figure 6b correspond to = 0.85. In both cases, the evolution for different values of N is shown. Each value on figure 6 correspond to the average resulting from 2000 independent simulations.
For the case = 0.25, figure 6a, three different regimes are observed, as indicated approximately in the figure. For slow quenches, region I, clearly shows finite size effects,   as the power law depends on the number of ions. Moreover, the power law converges as the N is increased. For regions II, and III, the slope is similar among different N indicating that finite size effects are less relevant for fast quench rates. Figure 6b correspond to = 0.85 and therefore to a later time of the quench. While regions I and III are essentially the same as for = 0.25, region II presents a clear difference : a significant portion of the defects have been lost. Moreover, a significant difference on the slope is now observed among the four ion number studied : the ion losses through the extremes are greater for smaller ion number chains.
In order to gain insight into the quench losses, the defects versus ln ν and have been plotted in a 2D pseudo-colour plot, figure 7.
Previous work done on harmonic chains have shown the presence of dynamics on the defects after their formation [20,19] and their manipulation has been recently reported [21]. In our analysis of single quenches where the trajectories were fully recorded shows axial displacements of the defects, sometimes leading to losses through the edges and defect oscillations. More importantly, if the density and the kinetic energy of the defects are very high, the probability that they annihilate each other becomes significant. The kinetic energy of the defects is not deterministic in nature but stochasticly driven.
Two possible phenomena can explain the lack of defect loses between figure 6a and 6b for regions I : that slow quenches do not generate energetic defects, or that when the number of defects is computed at the first value of , the energetic kinks have already left the system and only the defects with low kinetic energy remains. For fast quenches, region III, the absence of losses could be explained by the fact that the short times involved are not long enough for them to escape the system.
Such an interpretation agrees with the observations of figure 2, where the oscillations of the maximum perpendicular extension was clearly dependent on the quench rate. Another important remark is that the qualitative description provided above do seem to apply to the different ion numbers explored, as im figure 6.

Homogeneous KZM
The scaling law derived from the HKZM concerns the number of defects generated after crossing the critical point and therefore it does not taken into account possible annihilation of defects through defect recombination. It assumes, in addition to the homogeneity, an infinite system, and therefore does not into account the loss of defects through the edges either.
In order to minimise the effects of defect loss when comparing with HKZM, the curve corresponding to = 0.25, figure 6a, has been used to obtain the coefficient of the power law for different regions and different number of ions, n ∝ ν α , see figure 8.  numerical values do not match the predicted exponent of 1/3 [3]. : α N =30 = 0.338 ± 0.005, α N =64 = 0.324 ± 0.004, α N =128 = 0.324 ± 0.005 and α N =256 = 0.322 ± 0.004. These values are sensitive to the exact interval used and they can change significantly.
The result of the power law fits for the particular case of N=30, seem to hint to a (5/6, 2/6, 1/6) for the three different regions. It is also worth mentioning that, although the region, in the limit of large N, tends to 1/3, there is a clear distinction between the regions.

Correlation length
Deviations from the standard HKZM are expected to happen when the correlation length is comparable to the length scale of the system, in this case the ion-ion distance. This scenario should happen at very fast quenches. At the other extrema, for very slow quenches, the correlation length could be of similar order as our finite size system.
The correlation length, ξ, is defined in this context as the distance between two consecutive defects. The distances from the first and last defect to the edges have been excluded. The mean values for each quench time (corresponding to = 0.25 for full symbols and to = 0.85 for empty symbols), normalised by the ion-ion distance, ξ = ξ/d are shown in figure 9a.
Again, multiple regions appear. The most relevant information is the fact that the correlation length increases with time (or when losing defects), which would indicate that the defects are being annihilated inside the crystal rather than escaping through the edges. This hypothesis is reinforced by the plateaus observed for = 0.85. However, when comparing with figure 6, we observe that the regions where the = 0.25 and = 0.85 differs, are not the same for the − ln ν andξ. This is specially marked for slow quenches. For example, for N=30 and − ln ν = −5, where the correlation length does not change with but the number of defects changes. This can only be explained by a collective behaviour of the defects that are moving together towards an edge, hence being lost, while keeping the same correlation length.
Finally, it seems as the adiabatic correlation length is only achieved for N=30, where the correlation length shows a plateau at slow quenches corresponding roughly to 1/3 of the length of the chain, higher than the value obtained using using Ginzburg-Landau theory [18] : 1 2 √ 6 ≈ 0.204. A fit to the correlation length for the different regions (as defined in figure 8) has been done and the results are shown in figure 9b. The HKZM predicts an scaling of ξ ∼ ν −1/3 . Region II, represented as squares in figure 9b, approaches such a value.

Conclusion
We have presented a numerical study of an iso-spaced ion chain to study the HKZM. An new analyttic expression for the critical frequency, at which a uniform chain undergoes a structural phase transition to a zig-zag has been obtained and tested using MD and a Langeving integrator. It has been shown that the rate at which the control parameter is changed, modifies significantly the frequency at which the transition takes place.
By measuring the average defect formation for different quench rates and ion numbers, we could identify three different regimes with different power laws. In particular, we have shown that the HKZM regime could be measured experimentally with a relatively low number of ions (N=30). Moreover, by monitoring the number of defects and the correlation length as a function of the spatial extension of the global structure, we have found a significant defect losses.
The smallest chain studied N=30 is experimentally realizable. It has hinted to few interesting behaviours that will be numerically studied in detail in the near future : first, the collective defect behaviour which is deduced by comparing time evolution of the defect number and the correlation length and second, the hinted (5/6, 2/6, 1/6) scaling laws which, if confirmed, could lead to a better understanding of homogeneous but finite size systems.