Structure, dynamics and bifurcations of discrete solitons in trapped ion crystals

The suggestion to investigate discrete solitons (kinks) in trapped ion Wigner crystals has been put forward in [1], where structural and spectral properties of different types of topological defects were studied, and an analysis of the quantization of modes of discrete solitons was presented. Such discrete kinks were first realized in [2], and recently they were further studied both experimentally and theoretically in [3], where novel types of kinks with intriguing characteristics were discovered. The Kibble–Zurek mechanism [4–7], which allows us to predict the density of structural defects formed during a second-order phase transition, was analyzed for the finite and inhomogeneous ion trap setting [8, 9], and was recently studied in experiments observing the formation of kinks with trapped ions [10–12].

A basic property of solitons is their particle-like mobility. Understanding transport properties of discrete solitons is of importance for studying nonequilibrium processes in many nonlinear models widespread in physics [13], where discrete solitons (and more generally, discrete breathers [14]) form basic excitations of the system. Even in other fields such as chemistry and biology, discrete solitons have been proposed for describing dynamical processes in, e.g., DNA and biological macromolecules, or in chains of ions or molecules interacting through hydrogen bonds, where transfer of protons is described by propagating solitons [15].

One of the most important concepts with respect to propagating discrete solitons is the effective potential landscape that the kinks moving in a lattice are subject to, known as the Peierls–Nabarro (PN) potential [13]. The PN potential is typically periodic and can often lead to pinning of solitons in lattice sites, in particular because traveling solitons can radiate phonons and slow down. However, recently there has been rising interest in the study of novel phenomena associated with the PN potential in different models and schemes. Enhanced mobility of localized optical excitations in a coupled array of waveguides was discovered in [16], in connection with bifurcations of kink configurations. Exact vanishing of the effective PN potential for specific parameter values and velocities of the propagating solitons was found in some models corresponding to nonlinear optical media, in one dimension (1D) and two dimensions (2D) [17–22]. With discrete solitons of cold condensed atoms in optical lattices, mobility within the PN potential has also been studied [23], and possibilities of controlled 'routing' of matter-wave solitons have been suggested [24, 25].

Cold trapped ions offer a high degree of experimental control [2] allowing the study of the nonlinear physics of discrete solitons, possibly even in the quantum regime [1], which is usually experimentally inaccessible. Discrete solitons have also been recently proposed for high-fidelity entanglement generation in ion traps [26, 27], in the presence of noise and with prospects for studying open quantum system dynamics. Such applications and further fundamental studies of solitons in ion traps could benefit substantially from accurate control over the solitonic motional degrees of freedom, and in particular, from a deeper understanding of soliton mobility dynamics in the lattice.

In this paper, we study the structure and dynamics of discrete solitons in inhomogeneous, trapped Wigner crystals. In particular, we focus on the formation and properties of the PN potential as a function of trap parameters for various kink types. We analyze the bifurcations in the space of crystal configurations that lead to formation of the PN potential. We investigate further the unique phenomenon of the PN potential becoming a global trapping potential for certain kinks, with strong anharmonic modification of the kink motion. We study, both experimentally and theoretically, single-kink configurations, two-kink interactions and a novel type of kink, modified by a large mass defect.

In this section we study theoretically different kink solutions with 31 ions, restricted to a planar structure, where a single parameter controls the model, namely the ratio of the radial and axial confinement strength. We investigate the structure of the solution space of crystal configurations, and study how kinks of different types are related with other solutions. We analyze bifurcations which lead to kink solutions as the model parameter is varied, and explain their structure from considerations of symmetry and topology. All of these types of kinks were observed experimentally.

Crystals with N ions of a single species trapped in a quadrupole trap [28], can be described within the pseudopotential approximation by using non-dimensional potential energy in the form

$$\begin{equation} {V} = \sum _{i}^{N}\frac{1}{2} (x_{i} ^{2} +\omega_{y}^2 y_{i} ^{2} +\omega_{z}^2 z_{i} ^{2} ) + \sum _{i\ne j}\frac{1}{2} \parallel\!{ \skew3\vec{R}_{i} -\skew3\vec{R}_{j}}\! \parallel^{-1}, \end{equation} \tag{ 1 }$$

where $\skew3\vec{R}_{i} = \{x_{i},y_{i},z_{i} \}$ is the vector coordinate of ion i, time and distances are measured in units of 1/ω_x (with ω_x the axial trapping frequency) and $d= (e^{2} /m{\omega _x}^{2} )^{1/3}$, respectively, (with m and e the ion's mass and charge) and the energy is in units of mω²_xd².

For a small number of ions and if the axial confinement is much smaller than the radial, sufficiently cooled ions in a quadrupole trap form a crystalline chain along the trap symmetry axis. For a given number of ions and with lower radial confinement, the linear chain configuration becomes unstable, and the global minimum configuration takes the shape of a zigzag [29–37], which extends from the center of the trap, with possibly a part of the ions on the sides still in an axial chain. For parameters close enough to the transition from a linear chain, the zigzag is a planar structure. For radial confinement which is further decreased, a transition from the zigzag to a helix occurs for (nearly) degenerate radial frequencies [29, 30, 38–40], while planar and three-dimensional (3D) structures can also form for suitable trapping parameters (see e.g. [41–45]).

In this section we focus on planar crystals in the xy plane, i.e. such that ω_z ≫ ω_y, omitting the dependence on the z coordinate from the discussion. Within this description the structure of the planar crystals is therefore determined by the single parameter γ_y ≡ ω²_y/ω²_x. Let us define a crystal configuration as any critical point of the potential of equation (1), i.e. a solution of the equation ∇V =0. In other words, the static configurations are zeros of the vector field defined by ∇V . If we define the Hessian of a general configuration by

$$\begin{equation} K_{i,\alpha;j,\beta} = \frac{\partial^2V}{\partial \skew3\vec{R}_{i,\alpha}\partial\skew3\vec{R}_{j,\beta}}, \end{equation} \tag{ 2 }$$

where i,j = 1,...,N and $\alpha , \beta \in \left \{x,y\right \}$, then the configuration will be stable if all eigenvalues of K are non-negative, while if there are negative ones (which correspond to imaginary frequencies of the normal modes), the configuration is unstable.

We emphasize that, in this discussion, we are considering on equal footing all of the possible (stable and unstable) solutions, and not only the actual configuration that an ion crystal might take in an experiment. We are interested in exploring the structure of the space of crystal configurations, and for this we consider what solutions exist and how they are related for different values of the model parameter γ_y. Thus, we consider the solution space as a function of γ_y, but here we do not refer to the dynamics involved in an actual experiment which realizes a variation of the trapping parameters. The crystal configurations and also the eigenvalues of K vary continuously with the model parameter γ_y. In general, when varying one model parameter, only one eigenvalue of the Hessian reaches zero at certain isolated values of the parameter. This is the case for all configurations considered here. Under these conditions, and for γ_y > 0, when an eigenvalue of the Hessian reaches zero, a bifurcation of the solutions occurs, with solutions either changing stability, getting 'created' or 'disappearing'. These bifurcations can be classified and the qualitative properties of the solutions on both sides of the bifurcation can be understood from considerations of symmetry and conservation of topological numbers, as we detail below.

For any solution with no zero eigenvalue, the local index of a solution is defined by $i_0\equiv {\mathrm { sign}}\det K=\left (-1\right )^n$, where n is the number of negative eigenvalues. The sum of the local indices of all solutions emanating from the bifurcation point, is conserved separately on both sides of the bifurcations we consider, a result of topological degree theory [46]. At each bifurcation, the soft mode whose frequency is equal to zero is the mode that generates the spatial configurations after crossing the bifurcation point. The spatial symmetries of this mode determine the breaking of symmetry by the bifurcating solutions, as we detail below.

Figure 1 depicts schematically parts of the solution space of ion configurations for 31 ions and a given range of γ_y. There are many possible solutions that are not indicated, however most of them are unstable. For γ_y ≳ 152 the only stable configuration is the linear chain along the x-axis. This configuration remains a solution for any value of γ_y > 0, and its spatial structure is determined only by the number of trapped ions. Its axial and radial normal modes are exactly decoupled and the eigenvectors of the Hessian are independent of γ_y. In fact, the axial eigenvalues are also independent of γ_y while the radial eigenvalues depend linearly [47] on γ_y. Therefore, as γ_y is decreased below the critical value for stability of the linear chain solution, the radial normal modes of the chain lose stability one after the other at isolated bifurcation points, and these bifurcations generate new solutions in the solution space, as can be seen in figures 1 and 3.

Zoom In Zoom Out Reset image size

As the lowest-frequency radial mode of the linear chain crystal crosses zero (bifurcation a in figure 1), the linear chain solution bifurcates in a 'pitch-fork' bifurcation, well studied in the thermodynamic limit [33]. The linear chain becomes an unstable solution with one negative eigenvalue, and two new stable configurations appear in the solution space, the zigzag and its radial mirror image ('zagzig'), which become the new degenerate global minimum energy configurations.

The potential energy V of equation (1) is separately invariant under two symmetry operations, inversion about the x-axis and inversion about the y-axis, symmetries which we denote as $\mathbb {Z}_2^x$ and $\mathbb {Z}_2^y$. The linear configuration has both symmetries. At the zigzag bifurcation, when the number of ions is odd (bifurcation a in figure 1), the radial symmetry $\mathbb {Z}_2^y$ is broken, but not the axial symmetry, and this directly results from the structure of the soft mode, depicted in figure 1 as well. Since the chain solution is stable with index $i_0=\left (-1\right )^0=1$, and becomes unstable with i₀ = −1 at the bifurcation point where a single eigenvalue crosses zero, the symmetry 'protects' the spontaneously broken symmetric solutions, which must have index 1 (since the sum of the local indices is conserved), and therefore are stable in all modes.

Following the unstable axial chain configuration to lower values of γ_y, the odd kink solution is created as the next eigenvalue of the axial chain crosses zero (bifurcation b in figure 1). Applied to the kink configuration, the axial and radial inversions $\mathbb {Z}_2^x$ and $\mathbb {Z}_2^y$ coincide and are individually broken, but not their combination (the kink is invariant only under total inversion). This is a result of the structure of the mode generating this bifurcation, which shows the same symmetry (for an odd number of ions). Since the chain solution is unstable with index −1 on the left of the bifurcation and a further mode turns negative, necessarily all bifurcating solutions are unstable. At bifurcation c, which occurs at lower γ_y, on each branch of the unstable odd kink solution, the remaining total inversion symmetry is broken, and this time the kink solution becomes stable, with the two bifurcating (necessarily unstable) solutions now breaking the symmetry.

Therefore we find that only for parameter values γ_y ≲ 106, following bifurcation c, can kinks exist as metastable excitations in the crystal. The centered kink solutions arise from a bifurcation of the axial chain solution, when the second-lowest mode crosses zero, and another bifurcation is required to endow the kink with stability. The unstable solutions of bifurcation c, which are created in parallel with the kink solution gaining stability, are also of fundamental importance for understanding kink dynamics in the crystal. These unstable solutions are in fact 'saddle-point' solutions, which separate a centered kink from a displacement by one lattice site [48].

The kink solution that is stable at the center of the crystal cannot be found at any other site, neither stable or unstable—there is no kink configuration except at the center for γ_y ≳ 96.5. The zigzag is not wide enough to support any other kink in the chain. Only at γ_y ≈ 96.5 we find a bifurcation e which creates a solution of an odd kink displaced one lattice site from the center. Configurations of kinks displaced further away from the center will be created at lower values of γ_y via similar bifurcations.

Bifurcation e is different from the pitch-fork bifurcations described previously, due to symmetry considerations. The solutions that 'appear' for γ_y ≲ 96.5 and do not exist for a larger γ_y, have none of the inversion symmetries of the potential and do not emanate from any solutions having such a symmetry. The local index to the left of the bifurcation is 0 (since no solution exists), and therefore two solutions are created and split apart in solution space, one is stable and the other unstable⁵. Here as well, together with the (stable) displaced kink configuration, an unstable saddle-point configuration is created, which sets the barrier for the kink translation by one lattice site.

Following the saddle solutions of the centered odd kink configuration, bifurcation d is reached at γ_y ≈ 45.6. At this value of the parameter the crystal is wide, and the configuration and mode generating this bifurcation can be seen to have no inversion symmetries. Therefore, as in bifurcation e, this is a bifurcation where the saddle-point meets a stable solution, and both disappear to the right of the bifurcation. The stable solution is a displaced 'almost-extended' kink. For γ_y ≲ 45.6, the extended kink can only be found at the center of the chain, and here we see that this occurs as the displaced extended kinks disappear by colliding with their saddle point solutions. We will elaborate on this phenomenon in the following section, where we look into the properties of the PN potential in the setting of an inhomogeneous crystal.

The saddle solution that bifurcates together with every stable kink configuration (see figure 1) corresponds in fact to the PN potential barrier [13] which separates the kink from neighboring lattice sites. Thus, we see how the PN potential landscape forms in this model of an inhomogeneous discrete nonlinear system with a global trapping potential, through the same bifurcations in which the kink solutions are created. Moreover, examining the saddle solutions allows obtaining accurate quantitative properties of the PN potential barriers and local minima at each lattice site, as will be shown in the following.

Figure 2 shows the potential energy of three static kink configurations of 31 ions in a planar crystal as in figure 1. The figure shows the dependence of the energy on γ_y, in the regime where the configurations are stable. The energy is measured in nondimensional units, corresponding to V of equation (1), and the energy of the global minimum zigzag configuration has been subtracted at each point. From the data in the figure, the landscape of the PN potential can be accurately defined and extracted (for the depicted centered and single-site displaced odd kink), as presented schematically in inset (b). A measure of the energy in units of the Doppler-cooling limit temperature for the trapping parameters of the experiment is also given.

Zoom In Zoom Out Reset image size

For the odd kinks, which are defects contained within the zigzag boundaries, we see that in a large range of γ_y, the kink can be translated in the crystal, and successive bifurcations allow for stable kinks within local PN potential minima at lattice sites along the crystal. However, in some parameter regimes there exist stable solutions of kinks of different form, e.g. the extended planar kinks, that can reside only at the center of the crystal. This 'trapping' of kinks at the crystal center is the result of the PN potential being modified from a (modulated) periodic potential to a global trapping potential [3]. The PN potential is studied using a constrained energy minimization approach in [49].

The centered extended kink is stable in the regime 63 ≳ γ_y ≳ 25. For 63 ≳ γ_y ≳ 45 it appears as depicted near bifurcation d of figure 1 (experimentally imaged in figure 1(f) of [12]), while for lower γ_y values it becomes more and more extended as the crystal becomes wider and wider (as in [3, 10, 12]). The centered odd kink solution (of bifurcation c) becomes unstable through a bifurcation which occurs at γ_y ≈ 65.2. In this bifurcation the center ion of the chain leaves the origin, breaking the last remaining composite symmetry $\mathbb {Z}_2^x \mathbb {Z}_2^y$, and the resulting kink configurations (stable in a small parameter regime 64 ≲ γ_y ≲ 65.2) are similar to the (displaced) kink imaged in figure 13 of [2]. In fact, although not directly connected by a bifurcation, these solutions are very close in shape to the 'almost-extended' kink which is found for nearby parameters, as discussed above.

The transformations of different kink types as a function of the trapping parameters is an important property of kinks in Coulomb crystals. It has been used in [10] to evolve odd kinks into extended kinks for stabilization, and investigated in detail also in [12, 49]. It has been suggested in [27] to facilitate mobility of kinks in the lattice, since the mobility of kinks is enhanced near the bifurcation points, with the translational mode becoming soft or the configuration becoming unstable, or even the PN potential landscape changing qualitatively.

In figure 3 some configurations in the full solution space of 50 ions are considered, now without the restriction to planar crystals. Analyzing bifurcations which depend on two parameters is significantly more complicated than with one parameter, as simple bifurcations (with one mode crossing zero) occur generically along curves in parameter space. Bifurcations involving two modes simultaneously are also typical, occurring at points in two-dimensional parameter space. We restrict ourselves to just a few bifurcations leading to the formation of kinks, and vary only one parameter at a time.

Zoom In Zoom Out Reset image size

At bifurcation a the linear chain becomes unstable and the zigzag configurations are formed in the xy plane, for ω_z > ω_y. Here, for an even number of ions, the zigzags are invariant only under total planar inversion, i.e. under the composition of $\mathbb {Z}_2^x$ and $\mathbb {Z}_2^y$. At bifurcation b a centered kink is created. Now, only $\mathbb {Z}_2^y$ is broken. This kink is not the odd kink configuration discussed above, because the odd kink would require a breaking of $\mathbb {Z}_2^x$ symmetry. We follow this unstable configuration down to γ_y ≈ 121, while keeping ω_z ≫ ω_y, thus restricting the crystal to the xy plane. Now, lowering ω_z while keeping ω_y constant, at ω_z/ω_y ≈ 1.131 bifurcation d is hit, where the soft mode is transverse to the crystal plane (along the z direction), breaking $\mathbb {Z}_2^z$.

The mode that generates bifurcation d is a localized mode, with the eigenvector having significant components only on the center ions of the kink. The bifurcating solution is no longer 2D (for γ_y ≈ 121 held fixed, and ω_z/ω_y < 1.131), but rather the kink-core ions extend out of the crystal plane. This configuration becomes stable following a nearby bifurcation at e, giving rise to the 'blurred' kink solution found in [3]. The blurred kink at the center of the crystal is the only kink solution stable in this regime of parameters, for reasons which will be discussed below. Bifurcations c and f , which are related to the formation of two-kink configurations, will be discussed in section 5.

In this section we extend the analysis of the 'blurred' kinks, which were first discovered and studied in [3]. The center ions of these kinks appear blurred in the experimental images and this apparent blurring is reproduced in the numerical simulations using the trapping parameters of the experiment (see figure 4). In order to derive the trapping parameters, we reproduce the crystalline structures observed via a CCD camera, by using a molecular dynamics simulation and an optimization routine as detailed in appendix A. The blurring results from large amplitude and highly nonlinear oscillations of the center kink ions, which correspond to an excitation of the kink's lowest-frequency localized mode. This large and anharmonic oscillation is a direct probe of the shape of the local PN potential which becomes globally trapping, with the kink-core ions extending out of the crystal plane (see figure 5).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

As detailed in appendix A (see also [3]), at the considered experimental parameters, ω_z/ω_y ≈ 1.047. For this small radial indegeneracy, we find that the zigzag configurations starting with 53 ions are 'quasi-3D' structures, with a few of the center ions extending out of the zigzag plane in both directions, while the rest of the zigzag structure remains purely 2D. Since this transition starts with the ions at the center, and the ions' change in position is small and nearly transverse to the crystal plane, the zigzags of 53–57 ions still seem to be planar zigzags when projected on the camera plane (which is almost parallel to the crystal plane). In addition, all kink configurations, starting with 46 ions, are also quasi-3D structures, with at least the kink-core ions extending out of the crystal plane, in opposite directions. This happens because a kink is a defect of increased ion density, which can be reduced by pushing the two kink ions, which are closest, oppositely out of the crystal plane, if the trap confinement in this direction is not strong enough.

4.1. Linearized normal modes of oscillations

Figure 5 shows the configuration of a blurred kink in crystals of 50 ions and 51 ions, projected onto the crystal plane and transversally to it. It can be seen how the ions at the kink center are pushed out of the crystal plane in opposite directions. In order to understand the dynamics of the kink configuration, we expand the potential of equation (1) to second order about the ion positions in a minimum-energy configuration, $\{\skew3\vec{R}_{i} ^{0}\}$, and obtain the normal modes Θ_j in a well-known procedure (see e.g. [1]), by setting

$$\begin{equation} \skew3\vec{R}_{i,\alpha} \left(t\right)=\skew3\vec{R}_{i,\alpha} ^{0} +\sum _{j}^{3N}D_{i,\alpha}^{j} \Theta _{j} \left(t\right), \end{equation} \tag{ 3 }$$

where D^j_i,α is the matrix of normal mode vectors, with rows indexed by the N ion numbers i and three directions α, and columns by the 3N normal mode numbers j.

The frequency spectrum of normal modes of the kink configuration of figure 5(a), with 50 ions, is shown in figure 6(a), where the spectrum of the 50 ion zigzag is also plotted for comparison. The spectra of the two configurations overlap for the most part (there are differences in the normal modes). However, the notable localized low-frequency kink mode, which has in this case frequency ω_low ≈ 0.4ω_x, lies well below the lowest frequency of the zigzag configuration. The components of ion oscillations in this low-frequency mode vector are depicted by the arrows in figure 5(a). The lowered frequency of this mode is a result of the near degeneracy of the two radial trapping frequencies, and the motion of ions corresponds to oscillations in directions for which the restoring forces are minimal—i.e. which are 'empty' of ions, pointing away from the locally high density of the defect.

Zoom In Zoom Out Reset image size

For the kink with 51 ions in figure 5(b), it can be seen how due to the small radial indegeneracy, the kink-core ions are pushed out of the crystal plane in slightly asymmetric directions, with the localized low-frequency normal mode components expressing the asymmetry in the configuration. This explains the different blurring patterns of the kinks with an even and odd number of ions as observed in the experiment (figure 4).

In figure 6(c) the dependence of the localized low-mode's frequency is plotted versus the radial indegeneracy, ω_z/ω_y. Structural transitions occur at ω_low = 0, where this kink configuration loses stability. For intermediate values of the indegeneracy, ω_low increases and can surpass the axial c.o.m. mode's frequency (constant at ω_x ≡ 1), and even the second axial mode.

We note that for the low-frequency mode of the kink, there may be some deviations in frequency, introduced by the effect of micromotion [36, 50, 51], as compared with the pseudopotential analysis presented here. We have verified that these deviations do not alter the general results of the analysis, except for some fine tuning of the parameters, of the order of a few per cent [36], depending on the exact configuration.

4.2. Nonlinear oscillations

The strong Coulomb nonlinearity, which creates the kink and enables the high ion density at the kink center, means that the description of ion oscillations in terms of linearized decoupled modes will break at a certain amplitude of motion, when the nonlinear terms become significant. It was shown in [1] that different localized modes are coupled much more strongly than is typical in ion crystals. The blurring of the kink ions, observed in the experiment and in numerical simulations that were run at the Doppler cooling limit temperature (T_D), corresponds to a sufficiently large amplitude motion of the kink ions to leave the description in terms of a linear oscillator (of the low-frequency normal mode) inadequate.

An exact numerical simulation of the motion of the ions (in the potential of equation (1)) was run, with initial conditions corresponding to an excitation only of localized low-frequency mode of the kink, within the pseudopotential approximation. Figure 7 shows phase-space orbits of the oscillator corresponding to the low-frequency kink mode (obtained by inverting equation (3)), for four different runs with increasingly larger motional energy. Still, the outermost orbit corresponds to an initial excitation energy (only for the low-frequency mode) of ∼ 2/3 of T_D. The low-energy orbits are nearly elliptical, corresponding to harmonic motion. However, the oscillator is clearly anharmonic for motion amplitudes which are relevant for the experiment, with a temperature in the range T_D to 5T_D. Figure 7(b) shows an integration image corresponding to the outermost orbit, demonstrating that the typical blurring observed in the experiment can be obtained solely as excitation of the low-frequency mode, in a dynamic simulation without micromotion or laser-related effects. It should be noted that this does not mean that at T_D, the anharmonic normal mode remains an uncoupled oscillator—it is coupled to other degrees of freedom of the chain. However, when ω_low is tuned to a relatively high value (as in the middle region in figure 6(c)), which is an indication of the local potential turning stiffer, the large amplitude anharmonic motion is suppressed and the blurring of the kink ions disappears in the simulations.

Zoom In Zoom Out Reset image size

In the numerical simulations, blurred kinks, which are created displaced from the center of the crystal, immediately start to oscillate about its center. Figure 8(a) shows a simulated integrated fluorescence image of a kink 'released' from a few lattice sites to the side of the center. We note that in this simulation the system is initialized at 'zero' temperature, and the Doppler cooling is neglected. The blurring across the lattice reveals the oscillation of the kink along the chain. In the experiment, events like this require a careful choice of the parameters, since the motion gets easily damped, both by the cooling laser and by the radiation of phonons from the kink [13] as it moves in the PN potential.

Zoom In Zoom Out Reset image size

Experimental evidence of a structural defect hopping back and forth along the crystal on a timescale of ∼100 ms step⁻¹ is shown in figure 8(b), where the image consists of an integration accumulating light scattered during ∼6 s. A single CCD image (exposure 200 ms) is shown in figure 8(c), revealing the contribution of a different, dark ion species to the kink. From the exact numerical simulations including the time-dependent potential we deduce a dark ion mass of at least 40 atomic mass units. This could be explained by collisions with background gas (∼ 0.1 s⁻¹), such as ${\mathrm { Mg}}^+ + {\mathrm { H}}_2{\mathrm { O}} \rightarrow {\mathrm { Mg}}^+\left ( {\mathrm { H}}_2{\mathrm { O}} \right )$, due to a vacuum around 10⁻⁹ mbar.⁶ We observe melting events of the crystal, most likely induced by a collision with the residual gas, followed by the occurrence of a dark ion. Since we find that the dark ion has a high probability to turn into a bright ion after a subsequent collision, we conjecture it to be a molecular ion incorporating an Mg⁺ ion. The lack of closed electronic transitions in such ions explains why they are not directly observable via flourescent light. However, the dark ion gets sympathetically cooled by its coupling to the cold bath of already laser cooled Mg⁺ ions via mutual Coulomb interaction. The different mass/charge ratio causes it also to experience a reduced confinement by the mass-charge-dependent pseudopotential as compared with the Mg⁺ ions. As a consequence, the Coulomb repulsion additionally displaces the molecular ion out of the zigzag crystal to roughly twice the radial extent of the other ions, so in our case to ∼10 μm away from the trap axis. Furthermore, this displacement occurs also out of the crystal plane, therefore deforming the structure of the Coulomb crystal, either combined with a kink or as a configuration that corresponds to a ground state distorted zigzag.

The hopping of the massive kink in the crystal of figure 8(b) occurs on a long timescale, and characterizing the exact dynamics requires further studies. The controlled vibrational excitation of the massive kink exploiting its again altered and gapped spectrum in the PN potential and a fast imaging setup will allow us to study the kink with a mass defect and the related dynamics [35, 49]. In addition, the kink with a mass defect occurs also in combination with a second kink in the crystal (figures 9(a) and (b)), which leads us to consider configurations of two kinks [11, 12, 49] and their interaction.

Zoom In Zoom Out Reset image size

Returning to the bifurcation diagram of 50 ions in figure 3, at bifurcation c the third-lowest mode of the axial chain crosses zero and generates a solution, which is not invariant under inversions about x or y, but only under the combined operation. This solution is planar, incorporating two kinks, however it is unstable. The solution of two kinks side by side, which are stable at the experimental trapping parameters, is predicted after bifurcation f, where the z-inversion symmetry is broken, and the two solutions are not connected at the shown 'cut' through parameter space. An isolated simple bifurcation (with one mode crossing zero) connecting these two solutions is ruled out by considerations of index conservation (see section 2) and the single z-inversion symmetry which must be broken.

Moreover, from calculations similar to those related to figure 2, it is possible to find the energy gap compared to the ground state of the two-kink configuration. For 50 ions in the experimental parameters γ_y = 11² and ω_z/ω_y = 1.047, it is E_two-kink = 0.27, which should be contrasted with 2E_kink = 2 × 0.1265 = 0.253, where the energy of the single kink is taken at the center. Kink configurations with a larger separation can also be stable, and the dependence of the interaction energy of the kinks on their distance can be computed. Also of importance is the energy barrier for collision of the two kinks (which results in their annihilation). Using the bifurcation method presented in this paper, the barrier can be exactly defined and is computed to be E_barrier = 0.078, which corresponds for the experiment to ∼35k_BT.

Therefore, we find numerically and experimentally that at a sufficiently low temperature, which renders improbable high-amplitude oscillations leading to annihilation, two kinks can exist side by side, repelling each other and stabilized at a defined position by their interaction. Figures 9(c) and (d) show two such experimentally realized configurations.

Finally, the interaction of a kink with a mass defect and different blurred kinks (see figures 9(a) and (b)) also results in a change of the local potential landscape. In the case shown here, the presence of the massive kink shapes the local PN potential. As a directly detectable result, the kink no longer appears blurred, indirect evidence for the interaction of the two kinks causing an additional modification of the local PN potential. Along the lines of sections 2 and 3, a study of the PN potential bifurcations for such configurations can be performed, for example, by taking the mass of the defect ion to vary continuously as another parameter. In addition, the spectrum of the normal modes of two interacting kinks can be studied as in section 4. The localized modes of the two kinks will interact and allow for exchange of energy between these two kinks as well as with the rest of the mesoscopic ensemble of ions and collective modes of the crystal. This might become exploitable to study open quantum system dynamics [27].

To conclude, the concepts introduced in this paper, revealing the bifurcations of configurations in inhomogeneous trapped ion crystals, have enabled us to give an unambiguous definition of the PN potential landscape for one or more kinks. By following the bifurcations of relevant solutions in dependence on the trapping parameters, the existence and stability of different kink solutions can be predicted. We experimentally observe and numerically confirm how for different trapping parameters and related configurations, the kinks can be either highly mobile, pinned locally at lattice sites along the crystal or even trapped globally by a modified PN potential. Variation of the trapping parameters in a dynamic way allows for different kinks to be transformed into each other, or set to motion along the crystal. Controlled manipulation of kinks in Coulomb crystals will enable further studies into the dynamics of discrete solitons with unmatched experimental accessibility.

Studying the localized modes of kinks is a further valuable tool for understanding kink dynamics, both linear and nonlinear. Depending on the trapping parameters, kinks may have localized modes that, for cold crystals, can be well described in terms of linearized small oscillations, while for some configurations, the nonlinearity may play a crucial role even for laser-cooled crystals. The isolation of such collective degrees of freedom, or, on the other hand, the coupling of these nonlinear modes to their environment, is of fundamental interest. As studied in detail in [1, 27], cooling one or a few modes close to the quantum ground state may be feasible (in particular with the odd kinks, which manifest a gap separated high-frequency internal mode). Therefore, the ion trap system may allow reaching the quantum regime for motional modes localized to unique topological solutions, together with the ability to accurately measure the dynamics of these modes. This could allow for a variety of novel experimental studies—starting from open quantum system dynamics, to using solitons for conveying quantum information. In addition, quasi-2D zigzag phases similar to those realized in the ion trap system, have recently also attracted attention in the context of dipolar gases [52, 53] and 1D quantum wires [54–56] and the study of discrete solitons in such systems may be an intriguing further step.

The configurations presented here of two interacting kinks, either of the same type or combined with a mass defect, are a first step toward studying the interaction of solitons. We have shown that two-kink configurations realizable in the ion trap are stable for times which are orders of magnitude longer than the natural timescale of the dynamics in the system, and the barrier for annihilation of two solitons can be computed accurately, as has been demonstrated for a specific case. We also demonstrated that the interaction between kinks may cause a modification of the PN potential in their vicinity, leading to observable effects in their motion. Various promising possibilities exist for pursuing studies of kink interactions, first classically and eventually perhaps quantum mechanically, in this rich and complex nonlinear system.

TS was funded by the Deutsche Forschungsgemeinschaft (SCHA973). MM, TS and BR acknowledge the support of the European Commission (STREP PICC). BR acknowledges the support of the Israel Science Foundation and the German–Israeli Foundation. HL, TS and BR thank A Retzker, T Mehlstäubler and M B Plenio for fruitful discussions, held in particular during the 'Ulm Workshop on Theoretical and Experimental Aspects of Nonlinear Physics in Ion Traps'. HL thanks R Geffen.

In this appendix, we describe the exact form of the time-dependent trapping potential modeling the trap used in the experiments, and formulate how the orientation of the crystal plane can be calculated, taking into account the time dependence of the solution. The fitting routines which allow the accurate analysis of crystal structures and trap parameters are described. The calculations presented here were also used to derive the structure of crystals in [3].

A.1. A tilted trap potential

We start with the trapping potential of a linear Paul trap [28], taking x as the axial direction, and (in this section) measuring time and distances in units of 2/Ω (where Ω is the rf frequency) and $d= (e^{2} /4m{\Omega }^{2} )^{1/3}$, respectively (with m and e the ion's mass and charge). The non-dimensional trapping potential for a single ion is

$$\begin{equation} {V}_{\rm 4-rod} =\textstyle\frac{1}{2} (a_{x} x ^{2} +a_{y} y^{2} + a_{z} z^{2} -2q\cos2t(y^2 -z^{2} ) ) , \end{equation} \tag{ A.1 }$$

with q and $a_{\alpha }, \alpha \in \left \{x,y,z\right \}$, being the non-dimensional Mathieu parameters of the respective coordinates [28]. Since the Laplace equation requires that $\sum _{\alpha }a_{\alpha } = 0$, we take a_z = −a_x − a_y.

To the potential of equation (A.1) we add further harmonic terms pertaining to a dc trapping voltage in the radial plane, however whose principal axes are different, i.e.

$$\begin{equation} V_{\rm rot} = \textstyle\frac{1}{2} (a'z'^2 - a'y'^2), \end{equation} \tag{ A.2 }$$

where the y'z' plane is rotated with respect to the y and z coordinates by some angle θ, and V_rot obeys the Laplace equation by construction. Writing y',z' in terms of the original coordinates, we get the trapping potential in the form

$$\begin{equation} {V}_{\rm trap} = \textstyle\frac{1}{2} (a_{x} x ^{2} +a_y y ^{2} +a_{z} z ^{2} +a_{yz} y z - 2q\cos2t(y^2 -z^{2} ) ), \end{equation} \tag{ A.3 }$$

where the parameters have changed by

$$\begin{eqnarray*} &&a_y \to a_y -(\cos^2\theta - \sin^2\theta)a', a_z \to a_z +(\cos^2\theta - \sin^2\theta)a'\end{eqnarray*}$$

and we have a_yz ≡ 2a'cos θsin θ. We note that for θ = 45°, the original parameters a_y and a_z are unchanged, and a_yz ≡ a'.

For a single cold ion trapped at the origin of the potential V_trap of equation (A.3), the parameters of the potential will affect the frequencies and eigenvectors of its modes of oscillation about the origin. It is instructive to consider a few simple cases. The x motion is of course decoupled and depends only on a_x. Setting first a_yz = 0, the motion in the two radial directions (along the y- and z-axis) decouples and the modes will be aligned with those axes. For a_y = a_z = −a_x/2, the two radial frequencies will be exactly degenerate. This degeneracy can be lifted by varying a_y (while a_z is not an independent parameter, but determined by a_z = −a_x − a_y). Setting a_yz ≠ 0 couples (linearly) the motion in the radial plane, causing the mode vectors to mix, such that their direction is no longer aligned with the y- and z-axis [57–59]. The frequencies of the radial modes will change as well. The magnitude of the radial modes' indegeneracy, and the rotation of the principal axes, depend also on the value of the other parameters (q, a_y and a_z) in a complicated manner.

With crystals of multiple identical ions, the effects described above hold for the three c.o.m. modes, which are decoupled from other modes of the crystal. The pseudopotential approximation is constructed by taking the trapping potential to be time-independent, with the principal axes defined by the orthogonal directions of the c.o.m. modes, and trap strength given by the respective frequencies . In many cases, the structure of the crystal predicted within the pseudopotential approximation is very close to the crystal which is a solution of the full time-dependent potential [36]. However, it is known that there exist regions in parameter space, in which crystals may be found in 'peculiar' configurations, which do not exist within the pseudopotential approximation [50, 60–62]. The modes of oscillation of the ions in the time-dependent potential are in many cases also very close to the modes calculated from a pseudopotential, but there may be corrections to the eigenvectors or eigenfrequencies, which can be calculated [36, 50, 51]. In the cases examined here, the time dependence of the potential results in some tuning of the parameters, which was accurately accounted for.

A.2. Fitting the trapping parameters

We fit the experimental trapping parameters using an optimization routine which works in the following way. The free parameters in V_trap of equation (A.3) are taken to be a_x,a_y,a_yz and q. The theoretical crystal configuration for the given trapping parameters is obtained by a molecular dynamics simulation of the equations of motion of the ions with some small damping, which is adiabatically turned off. The resulting solution is verified to be a stable crystal configuration. If the initial conditions for the simulation correspond to a chain, the zigzag configuration is easily obtained. Finding a kink configuration requires starting with a zigzag configuration in which the radial coordinates of half the crystal ions were flipped.

The theoretical crystal at given parameters is then compared to the experimental image, by the following steps. The coordinates of all the bright ions are extracted from the experimental image. Then the ion coordinates are measured from the c.o.m., which coincides with the center of the trap for crystals with no dark ions. For crystals with dark ions, the center is taken to be the average of several all-bright configurations. The orientation of the camera plane with respect to the trap radial axes is parameterized by two additional free parameters, the azimuth and elevation angles, which are known to be close to 0° and −45°. The theoretical crystal's 3D coordinates are projected onto the plane of the camera, and then the sum of the distances between each ion in the experimental image and the theoretical projection is computed. A least squares fit is performed to minimize this sum in order to find the trapping parameters in combination with the point of view of the camera which best reproduce the experimental images.

Running the above analysis for a single-crystal image does not determine the parameters uniquely. However, already for the crystal with 39 ions (figure A.1(a)) it is immediately found that the crystal plane cannot coincide with the trap's y = 0 or z = 0 planes, which are at ±45° to the camera, since in this case, the width of the crystal is $\sqrt {2}$ wider than appearing in the image, and a matching configuration does not exist. Quantitatively however, the image could be generated by different combinations of trapping parameters, which correspond to a crystal plane which is either facing the camera (i.e. lying in the y = −z plane), or rotated from it at up to ∼ 22°. A nonzero a_yz is necessary for reproducing the observed crystals, in a plane which is viewed top-on from the camera. A zigzag crystal which lies in the y = −z plane can be reproduced in the simulations by using nearly degenerate radial trapping frequencies, which result in a 'peculiar' planar zigzag in the y = −z plane for 39 ions (the slight indegeneracy which is still present distinguishes this plane from the y = z plane). However, most of the other crystals which were observed in the experiment could not be explained as 'peculiar' crystals, since the observed configurations do not exist or are unstable when the radial frequencies are near-degenerate—they all require that the out-of-plane trapping frequency be larger then the in-plane one, by at least a few per cent.)

Zoom In Zoom Out Reset image size

In the experimental data, images of a zigzag configuration are found with up to 57 ions. For 58 and 59 ions, there is a variety of images of configurations with defects, and starting with 60 ions the configurations in the 2D projection seem 'wavy' (figure 9(e)). These crystals can be reproduced in the simulations as various 3D structures. Trying to fit unique trapping parameters for the different crystals between 39 and 57 ions, it is found that the outermost ions of the longest crystals observed are a bit too expanded sideways, an effect that can be explained by a small (∼  0.5%) decrease in the axial trapping on the sides, due to nonlinearity effects. Extensive searches in parameter space using dynamical simulations, show that all of the observed images can be accounted for quantitatively and the trapping parameters are derived. The average deviation between the experimental data and theory is smaller than 0.5 μm per ion, which is less than the pixel size (0.8 μm) and well within 1σ for all of the ions. The following parameters were obtained from the theoretical fitting : a_x = 0.000328, a_y = −0.0002, a_yz = 0.0019 and q = 0.286. The rf frequency from the experiment is Ω = 2π·6.22 MHz, and the ions are ²⁴Mg⁺. The camera plane is found to be at azimuth angle of ∼ − 1.92° and the elevation ∼ − 44.5°. The crystals lie almost exactly in the y = −z plane (in rotated coordinates y' = 0), and the secular frequencies are ω_x ≈ 2π × 56.7 kHz, ω_y' ≈ 2π × 623.3 kHz and ω_z' ≈ 1.047ω_y'.