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

J. Pedregosa-Gutierrez, M Mukherjee

► To cite this version:

J. Pedregosa-Gutierrez, M Mukherjee. Defect generation and dynamics during quenching in finite size homogeneous ion chains. New Journal of Physics, 2020, 22 (7), pp.073044. 10.1088/1367-2630/ab9a86. hal-02930370

HAL Id: hal-02930370

https://amu.hal.science/hal-02930370

Submitted on 4 Sep 2020

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
PAPER • OPEN ACCESS

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

## Recent citations

- Coulomb crystal splitting effect in a linear electrodynamic trap
Olga Kokorina et al

To cite this article: J Pedregosa-Gutierrez and M Mukherjee 2020 *New J. Phys.* **22** 073044

View the [article online](#) for updates and enhancements.
<figure><img src="image_2.png" /></figure>

CrossMark

OPEN ACCESS

RECEIVED

3 February 2020

REVISED

14 May 2020

ACCEPTED FOR PUBLICATION

8 June 2020

PUBLISHED

24 July 2020

Original content from
this work may be used
under the terms of the
Creative Commons
Attribution 4.0 licence.

Any further distribution
of this work must
maintain attribution to
the author(s) and the
title of the work, journal
citation and DOI.

<figure><img src="image_3.png" /></figure>

PAPER

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

J Pedregosa-Gutierrez¹,²,³ and M Mukherjee²,³,⁴

¹ Aix Marseille University, CNRS, PIIM, Marseille, France

² MajuLab, International Joint Research Unit UMI 3654, CNRS, Université Côte d'Azur, Sorbonne Université, National University of Singapore, Nanyang Technological University, Singapore

³ Centre for Quantum Technologies, National University Singapore, Singapore, 117543, Singapore

⁴ Department of Physics, National University Singapore, Singapore, 117551, Singapore

E-mail: jofre.pedregosa@univ-amu.fr

**Keywords:** RF ion trap, Kibble-Zurek, symmetry-breaking, phase transition

## Abstract

An equally spaced linear chain of ions provides a test-bed for studying the defect formation during a topological phase transition from a linear to a zig-zag configuration. By using a particular axial potential leading to an homogeneous ion chain, the boundary conditions are not needed, allowing new rich defect dynamics to appear on an homogeneous system. A semi-empirical expression for the critical transition frequency provides an excellent agreement to the numerical results for low ion numbers. The non-adiabatic crossing of the phase transition shows different power-laws for the defect probability density for different quench rates regions. Information regarding defect dynamics is obtained through the measurement of the defect density at different times during the quench. By comparing the defect density and the correlation length dynamics among the different number of trapped ions, the role of the different defect loss mechanism can be deduced. An excellent agreement with the predictions given by the homogeneous Kibble-Zurek model is found on a finite size system of 30 ion system which can be tested in present ion trap experimental set-ups.

## 1. Introduction

Phase transition is ubiquitous in nature as well as in the laboratories, from boiling of water to the formation of Bose-Einstein condensate. Symmetry breaking during the second-order phase transitions have been invoked to explain the observed inhomogeneity of mass distribution of the Universe through the Kibble-Zurek mechanism (KZM) [1, 2]. However, such a mechanism is universal to any phase transition crossing a critical point, making some of the aspects of the dynamics independent of the underlying system.

The original KZM (also known as the homogeneous KZM) assumes that the system is initially homogeneous (uniform density) and infinite (no edge boundaries). If the system is driven over a second-order phase transition to a ground-state with multifold degeneracy, the system has to choose one among the different possible final states. In an infinite and homogeneous system the whole system crosses to the new phase simultaneously. The speed at which the new phase choice can be communicated is finite, from which follows that there is a finite probability that different causally disconnected regions will choose different final states, leading to defect formation at the boundaries of different choices. Since these defects are structural in nature, they are known as topological defects.

One of the quantifiable predictions of the KZM is the number of topological defects formed as a function of the rate at which the critical point of a phase transition is crossed, also called the defect probability density. Such speed is controlled experimentally through the time variation of a control parameter like temperature in case of boiling water. Three main difficulties typically appear when
investigating 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 around the critical point. Second, the preparation of an homogeneous 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 predicted power law and experimentally measured one is achieved using a linear ion chains formed in linear RF traps [4, 5]. In this experiment, a finite size chain of ions trapped in a harmonic potential and laser cooled to form a linear chain of 1D Coulomb crystal was used. The chain was driven to a doubly-degenerate configuration of Coulomb crystal, namely zig-zag and zag-zig, based on whether the ions are push up or down with respect to the linear chain equilibrium point. However, such experiments could not explore the homogeneous case as the ion density in a harmonic ion trap is not constant [6]. This experiment was faithfully explained by modifying the original homogeneous case to the inhomogeneous KZM [7] (IKZM) in which the transition is crossed always at the place of highest density first i.e., at the center of the ion chain, leading to a modification of the scaling law.

As the ion trap technology is advancing rapidly, the HKZM could be studied using an uniform distribution of laser cooled ions. 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 expected [10, 11], with only one reported ion ring generated in a multipole surface trap [12]. A ring satisfies both the criteria namely, homogeneity and infinity.

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 ion trap quantum computing [13]. Here, we explore the possibility to use such a topology, the homogeneous ion chain of finite size, for the study of the defect dynamics under the ambit of HKZM. The KZM in general and the homogenous KZM in particular are becoming more and more relevant in the context of quantum technology. As an example, mesoscopic systems of uniform lattice structure are generated via classical phase transitions in well controlled experiments in order to develop super-lattice by molecular self-assembly for single photon generation [14, 15]. However, such systems are often bugged by unwanted defects leading to inefficient photon generations. In order to understand similar complex systems, it is important to design controlled experiments where defect formation during a phase transition can be better understood and controlled. One of the simplest but highly controllable system is a linear chain of uniformly spaced trapped and laser cooled ions. A structural or topological phase transition can be driven by varying only the ratio of the confinement strength in two dimensions in a well controlled manner, thus providing a toolbox to verify and control defect formation in any second order phase transition due to the universality of the KZM.

Unlike condensed matter systems, where the size often can be considered infinite but provides limited degree of controllability, an ion chain is a well-controlled Coulomb crystal that can be made homogeneous but finite in size. Finite size systems are going to be relevant in the near term quantum technology, thus knowing finite size effects in HKZM with edge boundaries, is likely to play a significant role. The universality of defect formation as a function of the rate at which the critical point of the phase transition is crossed is rather well known from the theoretical side. Such universality has also been tested for finite size systems with periodic boundary conditions [9, 16] and re-scaling has been suggested [9]. However, such studies do not include the possibility of defect loss which can be significant, particularly for the type of small systems which are experimentally accessible for near term quantum technology implementation. Moreover, such losses may even alter the universality of the scaling when experiments are performed to test these universality classes. In particular, we have used molecular dynamics (MD) simulations of the 1D to 2D (linear to zig-zag) transitions of a homogeneous ion chains to study the behaviour of defect formation as a function of the velocity by which the transition is crossed for a finite size system with edges. To the best of our knowledge, this investigation for the first time provides an insight on the losses, both at the boundaries and through annihilations.

The article is organised as follows: first, the electric potential needed to generate an homogeneous ion chain is discussed and verified through MD simulations. A new derivation of the critical parameter at which the linear chain to zig-zag transition occurs and MD simulations to verify the value of the critical parameter in function of the ion-ion distance and the number of ions are discussed next. Followed by a discussion in terms of defect generations of MD simulations of the linear to zig-zag phase transition over a large range of
the rate at which the critical parameter is changed. The results are then analysed in the context of the HKZM, followed by a discussion of the correlation length evolution.

## 2. Homogeneous linear ion chain

In a segmented linear ion trap the radial potential is generated by a radio-frequency field while the axial confinement is achieved by appropriate electrostatic voltages applied to the segmented electrodes along the axis. Without segmentation the axial potential is approximately harmonic in nature. However, the electrostatic potential required to generate an homogeneous ion chain can be obtained by finding an ion position dependent potential that compensates all Coulomb forces among the ions at each ion site [17]. The resulting analytical solution is:

$$ \phi_N(z) = \frac{Qk_C}{d} [2\psi^{(0)}(N_+) - \psi^{(0)}(\tilde{z}_+) - \psi^{(0)}(\tilde{z}_-)] \quad N_+ = (N+1)/2; \quad \tilde{z}_\pm = N_+ \pm z/d, \quad (1) $$

where the $\psi^{(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 [18] uses a different approach, where the ion chain is approximated as a uniform charge density $\rho_0 = Q/d$, leading to an expression for the effective electric field acting on each ion due to the Coulomb interaction with all the other ions on the chain. By integrating such an electric field, a different expression than the equation (1) was derived:

$$ \phi_L(z) = k_C \rho_0 \ln \left( \frac{L^2}{L^2 - z^2} \right) \quad (2) $$

where $L$ is the half length of the ion chain. The above expression leads to infinite potential walls as noted by the authors. Moreover, such potential leads to a 5% variation in homogeneity on their studied case (number of ions, $N=50$, and $d=3 \ \mu m$) [18].

In the following, through the use of MD simulations, it is shown that equation (1) leads to a degree of inhomogeneity much lower than 5%. For such reason, equation (1) has been used as axial potential in the rest of this work.

The experimental realisation of an homogeneous ion chain has been attempted in Xie et al [19]. A chain of 21 ions with an ion-ion distance of $8.5 \pm 0.39 \ \mu m$ was achieved using a surface ion trap with multiple electrodes. While the degree of homogeneity achieved at Xie et al [19], 4.6%, it is probably not enough for the KZM experiments, it is indeed reasonable to expect an improvement in the near future.

## 3. Molecular dynamic simulations

The problem at hand involves solving the equations of motion of $N$ interacting ions in the presence of the trap potential. The laser cooling is represented through a friction term while the different possible heating mechanisms in a real trap are taken into account through a single thermal bath. Such a problem can be described by the Langevin equation for each ion $i$ evolving as:

$$ \begin{aligned} m \partial_{tt} x_i &= Q^2 k_C \sum_{j=1, j \neq i}^{N} \frac{x_i - x_j}{|\vec{r}_i - \vec{r}_j|^3} - m \omega_x^2 x_j - \Gamma \partial_t x_i + \sqrt{2 \Gamma k_B T} \theta_{xj} \\ m \partial_{tt} y_i &= Q^2 k_C \sum_{j=1, j \neq i}^{N} \frac{y_i - y_j}{|\vec{r}_i - \vec{r}_j|^3} - m \omega_y^2 y_j - \Gamma \partial_t y_i + \sqrt{2 \Gamma k_B T} \theta_{yj} \\ m \partial_{tt} z_i &= Q^2 k_C \sum_{j=1, j \neq i}^{N} \frac{z_i - z_j}{|\vec{r}_i - \vec{r}_j|^3} - Q \left| \frac{dV(z)}{dz} \right|_{z_i} - \Gamma \partial_t z_i + \sqrt{2 \Gamma k_B T} \theta_{zj}, \end{aligned} \quad (3) $$

where $\vec{r} = (x, y, z)$, $\omega_x$ and $\omega_y$ are the secular frequencies along the transverse directions, $x$ and $y$ respectively, $V(z) = \phi_N(z)$ is the axial potential given by equation (1), $\Gamma$ 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 $\theta_{xj}$, $\theta_{yj}$ and $\theta_{zj}$ are a collection of independent standard Wiener processes [20]. The equations of motion are numerically solved using the vGB82 algorithm as described in [20].

MD simulations are used to test the degree of homogeneity of the ion-ion distance on the final ion distribution when the potential given by equation (1) is used. The parameters used are: $\omega_y/2\pi = 1$ MHz,
<figure><img src="image_4.png" /><figcaption>Figure 1. Molecular dynamic simulation of a quench for a system of 32 $^{138}$Ba$^{+}$. (a) Radial position of each ion vs the $\omega_x(t)$ evolution. (b) Ion configuration at $\omega_x(t = 0) = 352$ kHz. (c) Ion configuration at end of the simulation: $\omega_x(t = \tau_Q) = 302$ kHz.</figcaption></figure>

$\omega_x/2\pi = 500$ kHz, $\Gamma = 1.5 \times 10^{-20}$ kg s$^{-1}$ and $N = 128$ ions of $^{138}$Ba$^{+}$. The temperature is set to the typical Doppler limit value of $T = 0.5$ mK. The ions are initialised as a perfectly homogeneous chain with inter-ion distance $d = 10$ µm, followed by uniformly distributed random displacement between $\pm 1$ µm in both transverse and axial directions. The ions are evolved during 1 ms with a time step of 2 ns. If the average position of each ion for the last 0.5 ms of evolution is taken as their mean position, the final average ion-ion distance is $\langle \Delta z \rangle = 10.000 \pm 0.002$ µm, (one standard deviation is used as error measurement) which represent an excellent agreement with the designed inter-ion distance of 10 µm. If the temperature is increased to 5 mK, the homogeneity remains high with an average ion distance of $\langle \Delta z \rangle = 9.999 \pm 0.003$ µm.

## 4. Linear to zig-zag transition

The transition from linear to zig-zag occurs when the control parameter (in this case the radial trap frequency) is modified such that the system crosses the 1D to the 2D topological phase. In the following, we obtain both analytically and numerically, the value of the control parameter at which such a phase transition occurs.

In the present system, it is the value of $\omega_x$ that is reduced over time while keeping $\omega_y$ fix, 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 $\omega_y/2\pi = 1$ MHz is kept constant through the rest of the present work. Such value has been chosen so that $\omega_y > \omega_x$ in the range of $\omega_x$ used. Therefore, the transition always occurs along the $x$ direction. An example of such transition for 32 $^{138}$Ba$^{+}$ ions is shown on figure 1, where the transition has been crossed adiabatically. The phase transition is clearly visible in the figure 1(a), where the radial position as a function of $\omega_x(t)$ for all the 32 ions is plotted. The initial, figure 1(b), and final positions, figure 1(c), are also shown. In this example, the simulation parameters are: $d = 10$ µm, $\Gamma = 1.5 \times 10^{-20}$ kg s$^{-1}$ and $T = 1$ nK.

If periodic boundary conditions are used, the critical frequency, $\omega_c$, at which the transition occurs can be obtained analytically [21]:

$$ \omega_c^2 = \frac{Q^2 k_c}{md^3} 4 \sum_{j=1}^{N} \frac{1}{j^3} \sin^2 \frac{j\pi}{2} \quad (4) $$

$$ \omega_c^2(N \to \infty) = \frac{Q^2 k_c}{md^3} \frac{7\gamma(3)}{2} \quad (5) $$
<figure><img src="image_5.png" /><figcaption>Figure 2. Full symbols: ion configuration assumed for the derivation of the equation (8). Orange squares: ions with a net contribution on the total force on the central ion. Green empty circles: final ion configuration after an MD simulation where the transition chain to zig-zag is crossed adiabatically.</figcaption></figure>

where $\gamma(p)$ is the Riemann-zeta function.

The periodic boundary condition do not correctly represent a true finite size system. In the following, a more realistic situation has been considered, which leads to a better agreement with numerical results as shown in the following. Let us assume that we have a homogeneous zig-zag ion configuration, as given by the full symbols in figure 2.

We assume an odd number of N ions aligned in perfect zig-zag configuration along the x plane ($y = 0$), as shown by full circles in figure 2 with an 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 all the other ions cancels out. Therefore, the total force on the central ion can be written as:

$$|\vec{F}_c| = 2Q^2 k_c \sum_{j=1}^{N'} \frac{2h}{(4h^2 + (2j-1)^2 d^2)^{3/2}}, \qquad (6)$$

where $N' = \lfloor \frac{N+1}{4} \rfloor$.

The force due to the confining transverse potential is given by $|\vec{F}_{\text{trap}}| = -m\omega^2h$. By imposing equilibrium of forces we thus obtain:

$$\omega^2 = 4 \frac{Q^2 k_c}{m} \sum_{j=1}^{N'} \frac{1}{(4h^2 + (2j-1)^2 d^2)^{3/2}}. \qquad (7)$$

Taking the limit $h \to 0$, we arrive at:

$$\omega_c^2(N, d) = \frac{Q^2 k_c}{md^3} \sum_{j=1}^{N'} \frac{1}{(2j-1)^3}. \qquad (8)$$

Notice that equation (8) leads to a non continuous increase of $\omega_c$ with the number of ions N. For example, on figure 2, we see that, in the idealised final configuration, the net force on the central ion is exactly the same for the case N = 7 and 9 ions.

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, $\omega_a = \omega_x(t=0) > \omega_c$, with zero velocity. They are first thermalised at the initial radial frequency during 2.0 ms. Then the $\omega_x$ frequency is lowered linearly to a final frequency in the zig-zag phase, $\omega_b = \omega_x(t = \tau_Q) < \omega_c$ at a speed given by the quench rate:
<figure><img src="image_6.png" /><figcaption>Figure 3. Evolution of the maximum radial extension of the ions configuration ($N = 30, d = 10 \mu m$ and $T = 1$ nK) as a function of the $\omega_x$ for several quench rates, $\frac{\Delta\omega}{\tau_Q}$. (a) Dashed orange line: $10^{11}$ Hz s$^{-1}$, dash-dot green line: $10^{10}$ Hz s$^{-1}$, dotted red line: $10^9$ Hz s$^{-1}$ (b) Dotted red line: $10^9$ Hz s$^{-1}$, dash-dot green line: $10^8$ Hz s$^{-1}$, dashed orange line: $10^7$ Hz s$^{-1}$, solid blue line: $10^6$ Hz s$^{-1}$</figcaption></figure>

$\frac{\Delta\omega}{\tau_Q}$, where $\Delta\omega = \omega_b - \omega_a$ and $\tau_Q$ is the quench duration:

$$ \omega_x(t) = \begin{cases} \omega_a & \text{if } t \le 0 \\ \omega_a + \frac{\Delta\omega}{\tau_Q}t & \text{if } 0 < t < \tau_Q \\ \omega_b & \text{if } t \ge \tau_Q. \end{cases} \qquad (9) $$

The following MD simulations have been performed for $N = 30$ ions and $d = 10 \mu m$, using a thermal bath temperature of 1 nK. While this temperature represents a non-realistic experimental value, it seemed appropriate to use a very low temperature as equations (4) and (8) are derived assuming zero temperature. The critical frequency, obtained from equation (8) is $\omega_c = 327.38$ kHz. The values used for the initial and the final transverse frequencies are $\omega_b/2\pi = 140$ kHz and $\omega_a/2\pi = 500$ kHz.

The maximum transverse distance between the two ions $h$, as a function of the instant transverse secular frequency, $\omega_x(t)$, for different values of $\tau_Q$ is shown in figure 3. Note that the time axis runs from higher values of $\omega_x$ to lower ones. Topological phase transitions are clearly observed: for $\omega_x(t) > \omega_c$, the value of $h$ is very close to zero (for the slowest quench, the mean value of the $h$ before the transition is $\langle h \rangle = 0.34 \pm 0.17$ nm), which corresponds to a chain configuration, while for $\omega_x(t) < \omega_c$ there is a jump on the $h$ value, corresponding to the jump to the zig-zag phase.

Figure 3 shows a strong dependency between the exact value of the frequency at which the phase transition occurs, $\omega_c$ and the quench rate, $\frac{\Delta\omega}{\tau_Q}$. For the fastest quench studied, $\frac{\Delta\omega}{\tau_Q} = 10^{12}$ Hz s$^{-1}$, or $\tau_Q = 2.26 \mu 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 corresponds to the adiabatic case, where $\frac{\Delta\omega}{\tau_Q} \to 0$. The oscillations observed on $h$ just after the transition had been already observed in different system by Shimizu et al [22, 23] when studying the dynamics of a Mott insulator to a superfluid crossing.

A fourth order polynomial fit was found to correctly describe the evolution of $h$ after the transition has happened. For example, a fit on the slowest quench studied $\frac{\Delta\omega}{\tau_Q} = 10^6$ [Hz s$^{-1}$], see figure 3(c), leads to critical frequency of $\omega_c^{MD}/2\pi = 326.96$ kHz, obtained by finding the roots of the polynomial. The same numerical value is obtained if the data from the quench $\frac{\Delta\omega}{\tau_Q} = 10^7$ [Hz s$^{-1}$] is used instead. In order to save computational time, a quench of $\frac{\Delta\omega}{\tau_Q} = 10^7$ [Hz s$^{-1}$] has been used to perform MD simulations for different initial values of inter-ion distance, using $N = 32$ and $T = 1$ nK, leading to figure 4. The solid line
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Comparison between the critical frequency obtained from molecular dynamics (circles) and the theoretical value obtained from equation (8), for *N* = 32 ions (solid line) as a function of the initial inter-ion distance.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 5.** Comparison between the critical frequency obtained from MD simulations (full circles) and the analytical values given by equation (5) (red solid line), equation (4) (green stars) and equation (8) (orange dashed line). The purple dotted line corresponds to a fit to the MD simulations results.</figcaption></figure>

of figure 4 corresponds to equation (8), while the circles correspond to the critical frequencies deduced from the numerical simulations. The agreement between them is excellent.

If now the inter-ion distance is fixed to *d* = 10 µm and the number of ions is modified from 5 to 64, we obtain figure 5. In this figure the *ω*<sub>c</sub>(*d*) from equation (4), the computed values for *ω*<sub>c</sub>(*N*, *d*) from equation (8) and the results of the MD simulations are compared.

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 2. The configuration used for the derivation of
<figure><img src="image_9.png" /><figcaption>Figure 6. Random examples of final ion distribution after a quench. Several defects are indicated through boxes. Solid line box: intermediate defect; dotted line box: odd defect; dashed line box: extended defect.</figcaption></figure>

equation (8) is an oversimplification and therefore the total transverse force experienced by the central ion has extra terms which were not accounted for, leading to a lower value of the critical frequency. Also, the numerical results show a continuous increase of $\omega_c$ with the number of ions, in disagreement with equations (8) and (4) that have a non continuous increase. Nevertheless, equation (8) gives a hint on the nature of the N dependence. The following function has been use to fit the MD results:

$$ \omega_{\text{fit}}^2 = \frac{Q^2 k_c}{md^3} a_1 \sum_{j=1}^{N} \frac{1}{(j + a_2)^{a_3}} \quad (10) $$

leading to $a_1 = (6.8 \pm 0.3) \times 10^{-5}$, $a_2 = 1.90 \pm 0.03$ and $a_3 = 3.22 \pm 0.02$. The fitted function shows an excellent agreement with MD simulations results.

## 5. Defect generation during fast quenching in homogeneous chains

The focus is shifted now to the generation of defects during a non-adiabatic but linear crossing of the transition. As discussed in the introduction, different causally disconnected regions of the chain can make different choices of the two possible configurations at the new phase. Figure 1(c) is an example of zig-zag, but the symmetric (respect $x = 0$) configuration could equally have happened. At the boundary of two different topologies, a defect necessarily arises.

Two examples of the final configurations for $N = 64$ are shown in figure 6. The defects can be classified in three categories [24]: intermediate defect, with one ion at the trap axis ($|x| \approx 0$); odd defect, where two consecutive ions positioned 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$^5$. Such an approach detects all three types of defects simultaneously. Its implementation is efficient and universal as it does not need manual input or threshold.

The results shown in figure 3 indicates that the exact moment at which the phase transition takes place, and therefore the creation of eventual defects, depends on the quench rate used. Although the role of the temperature was not systematically explored, it is highly probable that it also plays a role to some extent on the exact frequency at which the transition takes place. The phase transition occurs after the adiabatic critical frequency has been crossed. It can even take place after the quench has finished. It implies that computing the number of defects at the end of the quench will lead to incorrect results. Therefore, in order to properly compare between different quench rates, we monitor the mean absolute transverse displacement, defined as $\langle x \rangle = \sum |x_i|/N$, and record the ion positions for post analysis when $\langle x \rangle$ reaches a particular fraction, $\epsilon$, of the adiabatic one, $\langle x \rangle = \epsilon \langle x \rangle_{\text{adiabatic}}$.

Another important aspect to take into account is defect loss which, unlike other previous works, it is clearly observable in the present study. Indeed, there are two mechanisms by which the number of defects can decrease during and after a quench: annihilation of the defect pairs and losses through the edges of the

<footnote>$^5$ The algorithm can be written in a single line when using Python/Scipy: $n = size$ (where (abs(diff(sign(x))) < 1)).</footnote>
<figure><img src="image_10.png" /><figcaption>**Figure 7.** Number of defects versus the dimensionless quench rate (see main text). (a): $\langle r \rangle = 0.25 \langle r_a \rangle$; (b): $\langle r \rangle = 0.85 \langle r_a \rangle$. Small full circles: $N = 30$; large full circles: $N = 64$; small empty circles: $N = 128$; large empty circles: $N = 256$. The standard deviation of each point is smaller than the symbols used.</figcaption></figure>

system. In order to gain insight on the evolution of the number of defects during a single quench, the average defect number was measured for $\epsilon = 0.25$ to $\epsilon = 0.85$ with steps of $d\epsilon = 0.05$. For the smaller values of $\tau_Q$ studied, the quenches are too fast for the ions to reach the required values of $\epsilon$ during the quench duration itself. In those situations, the MD simulation is continued with $\omega = \omega_b$ until $\epsilon = 0.85$ is reached.

For the simulation results discussed henceforth, we have used $^{138}$Ba$^+$ ions, using a thermal bath of $T = 0.5$ mK and a friction coefficient of $\Gamma = 1.5 \times 10^{-20}$ kg s$^{-1}$. The quench starts at $\omega_a/2\pi = 500$ kHz and finishes at $\omega_b/2\pi = 140$ kHz. The values of the quench duration, $\tau_Q$, have been explored over a large enough range. Four different ion numbers have been studied: $N = 30, N = 64, N = 128, N = 256$. In all the cases, the ions are initialised with an ion-ion distance of $d = (10 + \text{rand})\mu\text{m}$, where rand is a random number in the range $\pm 1$ and evolved for 2 ms before the quench starts.

The evolution of the logarithm of the averaged normalised number of defects, $\ln(n/N)$, versus the logarithm of the dimensionless quench rate, $\Lambda = -\ln(\tau_Q\omega_0)$, with $\omega_0^2 = \frac{Q^2k_B}{md^2}$ [9], is shown in figure 7 for $\epsilon = 0.25$ and $\epsilon = 0.85$. The results for the different ion numbers, $N$, are shown. Each point on figure 7 corresponds to the average of 2000 independent simulations. The standard deviations on the normalised number of defects, $n/N$, ranges from $10^{-4}$ and $10^{-3}$. Those values are much smaller than the symbols used on figure 7.

Three regions can be identified in figure 7. Region III where fast quench occurs, does not show a significant difference on the number of defects between $\epsilon = 0.25$ and $\epsilon = 0.85$. The values for $\ln(n/N)$, and therefore the slope, are similar among the different studied ion numbers. This is in contrast with the slow quench regime, region I, where although there is no significant defect loss between $\epsilon = 0.25$ and $\epsilon = 0.85$, there is a clear dependence on the ion number.

The absence of losses in region I can be explained by a very slow of dynamics of the defects generated, slower than the time difference between $\epsilon = 0.25$ and $\epsilon = 0.85$ or that any defect dynamics, and therefore any loses, have already taken place at $\epsilon = 0.25$ due that the slow quench regime of this region. For fast quenches, region III, the absence of losses could be explained by the fact that the times involved are too short for any dynamics of the defects to take place.

Region II presents a clear evolution of the number of defects with the value of $\epsilon$: a significant portion of the defects present at $\epsilon = 0.25$ have been lost when $\epsilon = 0.85$ has been reached. Moreover, while the power-law is similar for region II at $\epsilon = 0.25$ across the different values of the ion number $N$ studied, see figure 7(a), the slope increases at $\epsilon = 0.85$ in region II as the ion chain size decreases, see figure 7(b).

This behaviour can be interpreted using the work of Partner et al [24] on the dynamical behaviour of the defects created in the context of single defects in laser cooled ions on ‘standard’ ion traps, where the harmonic potential along the axis leads to the creation of defects in the middle of the chain. Through the calculation of the Peierls–Nabarro (PN) potential, Partner et al [24] was able to explain qualitatively the dynamics of such defects. The PN potential corresponds to the potential seen by the defect along the chain
<figure><img src="image_11.png" /><figcaption>**Figure 8.** Plot of $\frac{nN^2}{(\tau_Q\omega_0)^{2/3}}$ vs $\frac{N^3}{\tau_Q\omega_0}$. (a) $\epsilon = 0.25$, (b) $\epsilon = 0.85$. Blue circles: $N = 30$; orange squares: $N = 64$; green circles: $N = 128$; empty red circles: $N = 256$.</figcaption></figure>

as a function of the defect center (see [24] for a definition of defect center and how to compute the PN potential). The PN potential rises towards the edges for the *extended* defects while decreases for the other two types, *odd* and *intermediate* defects. Moreover, Partner *et al* [24] also showed how defect dynamics imply the transformation among the different types of defects and, as a consequence, the defect sees alternatively a confining and de-confining PN potential. While the above results were obtained in an ion chain trapped using an axial potential and with a non-uniform ion distribution, the PN potential on an homogenous chain, although flatter at the center, should be qualitatively similar.

Full ion trajectories were recorded for some single quench rates. Their analysis showed axial displacements of the defects, sometimes leading to losses through the edges while other times the defect oscillated along the chain without being lost. Such observations have no statistical meaning, but give us some hint of the type of dynamics that can be expected and they are in accordance to what was reported in [24]. It also provides insight on how edge plays a role in defect dynamics.

Therefore, the difference in the power laws observed in region II when $\epsilon = 0.85$, can be explained by the fact that, as we are working in an homogeneous system, the probability of defect creation is uniform along the chain and therefore the probability that a defect is created at the edge is higher for smaller systems than for larger systems. Once the defects are created, the ones closer to the edge are more easily lost if their corresponding PN potential is a de-confining one.

The same reason can be invoked for the *knee* observed between region II and region I observed for $\epsilon = 0.85$ and not for $\epsilon = 0.25$. The strong *knee* observed for $N = 256$ at $\epsilon = 0.85$, coupled with its absence at $\epsilon = 0.25$ can be explained due to a decrease of the absolute number of defects between both values of $\epsilon$, implying a decrease of defect density in the chain, which in turns means that the probability of two defects annihilating each other during a quench should decrease. Such interpretation also explains the decrease of the *knee* with the decrease of the number of ions: the relative annihilation contribution to the defect losses during a quench is necessarily lower for smaller chains.

Finally, it should be noted that by normalizing the number of defects with respect to the total number of ions, $n/N$, clearly converges as the ion number increases. The numerical values for $N = 128$ and $N = 256$ are indeed very close, specially for region II and III when $\epsilon = 0.25$, indicating the approximation of an infinite system should be valid at this regions and for such value of $\epsilon$. The finite size universality for an homogeneous system can be recovered by using a different scaling if $\frac{nN^2}{(\tau_Q\omega_0)^{2/3}}$ vs $\frac{N^3}{\tau_Q\omega_0}$ is plotted, as it was shown by Nigmatullin *et al* [9]. Our results, plotted using such scaling, are shown in figure 8. The different ion number results collapse, as an overall trend, into a single line thus confirming the results from Nigmatullin *et al* However, a closer look clearly indicates some deviations, coming from slow quenches set of results. The differences from our results and the ones reported in [9] is that in Nigmatullin *et al* periodic boundary conditions are imposed. These boundary conditions forbids any type of defect loss as the PN potential will not decrease at the edges, since there exists no edges, and therefore, no defect dynamics are expected.
<figure><img src="image_12.png" /><figcaption>**Figure 9.** Number of defects versus $\Lambda$ and $\epsilon$ for several values of the friction coefficient, $\Gamma$, for $N = 64$ and two values of $\epsilon$: (a) $\epsilon = 0.25$ and (b) $\epsilon = 0.85$. Blue small dots: $\Gamma = 1.5 \times 10^{-21}$ kg s; orange empty squares: $\Gamma = 1.5 \times 10^{-20}$ kg s; green empty circles: $\Gamma = 1.5 \times 10^{-19}$ kg s$^{-1}$.</figcaption></figure>

## 6. Homogeneous KZM

The present work differs from the standard HKZM in two main aspects. First, the scaling law derived from HKZM theory for the number of defects generated does not take into account possible defect losses. Secondly, the presence of laser cooling, modelised here as a damping term, was not present in the original KZM theory. The introduction of a damping term has been studied within the Ginzburg–Landau description by several authors, for example Chiara et al [25]. Two specific limits have been analysed on the context of infinite homogeneous ion chains [25]: the over-damped and the under-damped, leading to $n \propto (\tau_Q\omega_0)^{-1/4}$ and $n \propto (\tau_Q\omega_0)^{-1/3}$, respectively.

It is this second scenario, the under-damped limit, that leads to the same scaling as the standard HKZM. In Nigmatullin et al [9], it is argued that in the under-damped limit, a small change on the damping term should not significantly affect the results. For this reason, we have performed MD simulations for three different values of $\Gamma$: $1.5 \times 10^{-21}$ kg s$^{-1}$, $1.5 \times 10^{-20}$ kg s$^{-1}$ and $1.5 \times 10^{-19}$ kg s$^{-1}$. For this comparison only $N = 64$ have been used. The rest of the parameters and procedure remain the same as before. The results, figure 9, show how the results corresponding to $\Gamma = 1.5 \times 10^{-21}$ kg s$^{-1}$ and $1.5 \times 10^{-20}$ kg s$^{-1}$ are essentially the same. Therefore, our simulation parameters over the quench rate range explored can be safely assumed to belong to the under-damped limit.

In order to minimise the effects of defect loss, the curve corresponding to $\epsilon = 0.25$, figure 7(a), has been used to obtain the coefficient of the power law for the different regions and different number of ions, $n \propto (\tau_Q\omega_0)^{-\alpha}$, see figure 10.

Figure 10(c) shows how region II agrees with the expected value corresponding to the HKZM, $\alpha = 1/3$ [3]: $\alpha_{II}[N = 30] = 0.338 \pm 0.005$, $\alpha_{II}[N = 64] = 0.324 \pm 0.004$, $\alpha_{II}[N = 128] = 0.324 \pm 0.005$ and $\alpha_{II}[N = 256] = 0.322 \pm 0.004$. While the agreement with the expected value is relatively good, the values are consistently but only slightly lower. The reason is probably that the results are already affected by defect losses at $\epsilon = 0.25$. In addition, it should be noticed that these values are sensitive to the exact interval used and they can change significantly.

Interestingly, a new scaling law seems to appear at region III $\alpha_{III} \approx 1/6$. The values are $\alpha_{III}[N = 30] = 0.174 \pm 0.003$, $\alpha_{III}[N = 64] = 0.160 \pm 0.005$; $\alpha_{III}[N = 128] = 0.163 \pm 0.003$ and $\alpha_{III}[N = 256] = 0.162 \pm 0.003$. The origin of this new scaling, lower to than the standard HKZM, is unknown although some possible explanations are hinted in the next section.

Finally, three different regions were already observed in the case of the IKZM, see figure 2 of Pyka et al [4], corresponding to 1/3, 4/3 and 8/3 for (qualitatively speaking) fast, intermediate and slow quenches respectively. The change from 1/3 to 4/3 appears from the restriction of the region where defects can appear, due to the highest density of ions at the ion chain centre when a harmonic potential is used to confine the
<figure><img src="image_13.png" /><figcaption>**Figure 10.** Number of defects versus the adimensional quench rate for different values of N: (a) N = 30; (b) N = 64; (c) N = 128; (d) N = 256. The power low fit to each region is shown. (e) The values of the fit (multiplied by 3) are shown in function of the number of ions. Circles: region I; squares: region II; diamonds: region III. Notice that the error bars associated with the fits are much smaller than the symbols.</figcaption></figure>

ions [7]. In the present case, the homogeneity of the initial ion chain assures that the defect formation is not restricted spatially, independently of the quench rate used.

For slow quenches, the correlation length, ξ, defined in the present context as the distance between two consecutive defects, becomes comparable with the size of the system and the density of defects equals the probability of a single defect, leading to a doubling of the previous power law to 8/3 [7]. The same argument could be invoked here as a factor 2.5 appears for the smaller chain studied, N = 30. As the length of the ion chain increases, the power law of this slow quench region approaches 1/3 of the HKZM. This is consistent with the evolution of the correlation length respect the system size, $\xi_L = \frac{\xi}{d(N-1)}$ (see next section), where the correlation length decreases as the ion number increases for the slowest quenches studied.

## 7. 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, as the low-energy Ginzburg-Landau theory is no longer valid. At the other extrema, for very slow quenches, the correlation length could be of similar order as our finite size system, the ion chain length.

Therefore, the evolution of the correlation length has also been extracted from the MD simulations. The mean values of the correlation length for each quench time (corresponding to $\epsilon = 0.25$ for full symbols and to $\epsilon = 0.85$ for empty symbols), normalised by the ion-ion distance, $\xi = \xi/d$ are shown in figure 11(a) where the standard deviation is smaller than the size of the symbols used. The distances from the first and last defect to the edges have been excluded. This new figure shows how the average correlation length increases with the parameter $\epsilon$, and therefore, with time.

Once more, several regions with different behaviours are observed. The regions do not share the same boundaries for the different ion numbers and therefore, the regions marked with different tones of blue
<figure><img src="image_14.png" /><figcaption>**Figure 11.** (a) Normalised correlation length vs normalized quench rate for the different ion numbers studied at two different values of $\epsilon$: full symbols, $\epsilon = 0.25$; empty symbols $\epsilon = 0.85$. Blue small circles: $N = 30$; orange squares: $N = 64$; green large circles: $N = 128$; red diamonds: $N = 256$. The two vertical lines indicate the boundaries between the regions observed in figure 7. The (b)–(d) figures correspond to the coefficient of the power law fits for region I, region II and region III respectively.</figcaption></figure>

indicate the data sets used for the power law fits. The two vertical lines indicate the boundaries between the regions observed in figure 7.

When comparing with figure 7, we remark that the regions where the $\epsilon = 0.25$ and $\epsilon = 0.85$ differs, are not the same for the $\Lambda$ and $\xi$. This is specially marked for slow quenches. For example, for $N = 30$ and $\Lambda = -5$, where the averaged correlation length does not change as $\epsilon$ increases but the number of defects does decrease as $\epsilon$ increases. If there was annihilation of pair of defects, the correlation length should increase, which is not the case. With the current information at hand, a constant correlation length coupled with a loss of defects can only be understood as a collective behaviour happening at the edges. Further theory as well as simulation work is needed to elucidate the observed behaviour which is beyond the scope of this work.

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 $L/3$. If the $L/3$ should also be applied to the other studied cases, a plateau should be observed at: $\approx 21$ for $N = 64$, $\approx 43$ for $N = 128$ and $\approx 85$ for $N = 256$. We see that $N = 64$ was just reaching $\hat{\xi} \approx 21$ for the slowest of the simulated quench rates, while the higher values of $N$, were still far from it. The theoretical value obtained using using Ginzburg–Landau theory [16]: $\frac{1}{2\sqrt{6}} \approx 0.204$, which is lower than our measured value of $1/3$.

The scaling law of $1/6$ observed in region III of figure 7 corresponds approximately to $\hat{\xi} < 5.5$. Such low value of the $\hat{\xi}$ could explain the change on the scaling as the model used for the derivation of the HKZM loses its validity [9].

A power law fit to the correlation length for the different regions defined by the different blue regions in figure 11(a) has been done. The HKZM predicts an scaling of $\hat{\xi} \sim (\tau_Q\omega_0)^{1/3}$. Region II, represented as squares in figure 11(c), seem to approach such a value as the ion chain increases. The correlation length appears to be more sensitive to finite size effects than the number of defects as the value has not yet converged to the infinite case, $1/3$ when $N = 256$. The power law in region III, figure 11(c), seem to converge again to the same scaling of $1/6$, already observed in the defect density in this region.
## 8. Conclusion

We investigated the defect dynamics during a 1D to 2D phase transition of an uniform ion chain confined in a linear ion trap using MD simulations. This allowed us for the first time to study a homogeneous but finite size system without the need of periodic boundary conditions, and therefore introducing the possibility of defect dynamics in presence of defect losses.

A new analytic expression for the critical frequency at which a uniform linear chain undergoes a structural phase transition to a zig-zag has been derived. Using MD simulations, the critical frequency has been obtained for a range of different ion numbers and ion-ion distances. The numerical results showed a better agreement with the new expression but the deviation for small ion numbers is still significant. However, a fit to a function inspired by this new expression, leads to an excellent agreement. The MD simulations showed that the rate at which the control parameter is changed, modifies significantly the critical frequency at which the transition takes place, motivating the introduction of the $\epsilon$ parameter.

The role of the $\epsilon$ parameter has proven critical in order to understand defect losses dynamics and to identify the quench rate parameter range over which the HKZM scaling can be obtained. Furthermore, it allowed us to identify three different regimes regarding the defect generation scaling and the defect losses for a given quench rate. In the region that better agrees with the expected HKZM scaling, defect loss have been shown to depend on the system size, given here by the ion number $N$.

We further investigated the dependence of the correlation length evolution with the quench rate. This combined with the loss dynamics for the same quench rate indicates the appearance of new collective defect behaviour at edges. Finally, we have shown that the HKZM regime emerges with a relative low ion number, $N = 30$. Such ion number value is experimentally accessible and it could open the door to an HKZM experimental verification on a system with an extreme controllability and with a high degree of repeatability, essential for the building up of the needed statistics of such types of experiments.

## Acknowledgment

This work is supported by Singapore Ministry of Education Academic Research Fund Tier 2 (Grant No. MOE2016-T2-2-120) and NRF-QEP (Grant No. QEP-P6).

## ORCID iDs

J Pedregosa-Gutierrez https://orcid.org/0000-0001-7527-4522

## References

[1] Zurek W H 1985 *Nature* **317** 505–8

[2] Kibble T W B 1976 *J. Phys. A: Math. Gen.* **9** 1387–98

[3] del Campo A and Zurek W H 2014 *Int. J. Mod. Phys. A* **29** 1430018

[4] Pyka K et al 2013 *Nat. Commun.* **4** 2291

[5] Ulm S 2013 *Nat. Commun.* **4** 2290

[6] Kamsap M R, Champenois C, Pedregosa-Gutierrez J, Mahler S, Houssin M and Knoop M 2017 *Phys. Rev. A* **95** 013413

[7] del Campo A, De Chiara G, Morigi G, Plenio M B and Retzker A 2010 *Phys. Rev. Lett.* **105** 075701

[8] Champenois C, Marciante M, Pedregosa-Gutierrez J, Houssin M, Knoop M and Kajita M 2010 *Phys. Rev. A* **81** 043410

[9] Nigmatullin R, del Campo A, De Chiara G, Morigi G, Plenio M B and Retzker A 2016 *Phys. Rev. B* **93** 014106

[10] Pedregosa-Gutierrez J, Champenois C, Kamsap M R, Hagel G, Houssin M and Knoop M 2018 *J. Mod. Opt.* **65** 529–37

[11] Pedregosa-Gutierrez J, Champenois C, Houssin M, Kamsap M R and Knoop M 2018 *Rev. Sci. Instrum.* **89** 123101

[12] Li H K et al 2017 *Phys. Rev. Lett.* **118** 053001

[13] Lin G D, Zhu S L, Islam R, Kim K, Chang M S, Korenblit S, Monroe C and Duan L M 2009 *Europhys. Lett.* **86** 60004

[14] Rainò G, Becker M A, Bodnarchuk M I, Mahrt R F, Kovalenko M V and Stöferle T 2018 *Nature* **563** 671–5

[15] Li S, Duan S, Zha Z, Pan J, Sun L, Liu M, Deng K, Xu X and Qiu X 2020 *ACS Nano* **14** 6331–8

[16] Puebla R, Nigmatullin R, Mehlstäubler T E and Plenio M B 2017 *Phys. Rev. B* **95** 134104

[17] Johanning M 2016 *Appl. Phys. B* **122** 71

[18] Leung P H and Brown K R 2018 *Phys. Rev. A* **98** 032318

[19] Xie Y, Zhang X, Ou B, Chen T, Zhang J, Wu C, Wu W and Chen P 2017 *Phys. Rev. A* **95** 032341

[20] Skeel R D and Izaguirre J A 2002 *Mol. Phys.* **100** 3885–91

[21] Fishman S, De Chiara G, Calarco T and Morigi G 2008 *Phys. Rev. B* **77** 064111

[22] Shimizu K, Kuno Y, Hirano T and Ichinose I 2018 *Phys. Rev. A* **97** 033626

[23] Dutta T, Mukherjee M and Sengupta K 2012 *Phys. Rev. A* **85** 063401

[24] Partner H L, Nigmatullin R, Burgermeister T, Pyka K, Keller J, Retzker A, Plenio M B and Mehlstäubler T E 2013 *New J. Phys.* **15** 103013

[25] Chiara G D, Campo A d, Morigi G, Plenio M B and Retzker A 2010 *New J. Phys.* **12** 115003