Relativistic quantum mechanics with trapped ions

We consider the quantum simulation of relativistic quantum mechanics, as described by the Dirac equation and classical potentials, in trapped-ion systems. We concentrate on three problems of growing complexity. First, we study the bidimensional relativistic scattering of single Dirac particles by a linear potential. Furthermore, we explore the case of a Dirac particle in a magnetic field and its topological properties. Finally, we analyze the problem of two Dirac particles that are coupled by a controllable and confining potential. The latter interaction may be useful to study important phenomena as the confinement and asymptotic freedom of quarks.


Introduction
The field of quantum simulators is one of the most rapidly growing in quantum information science. It was Feynman [1] who stated thirty years ago that a controllable quantum device could emulate the dynamics of another quantum system exponentially faster than a classical computer. Since then, this hypothesis has been confirmed [2] and, subsequently, intensive theoretical and experimental research has followed [3,4].
Quantum simulators of quantum relativistic systems are among the most elegant and, among a few others, the most promising candidates for going beyond classical computational capabilities. In the last years, the simulation of black holes in Bose-Einstein condensates (BEC) [5] and properties of the expanding universe [6,7] have been proposed, and an experiment on BEC sonic black hole has been realized [8]. A relevant milestone was given by the proposal for the simulation of Dirac equation and associated quantum relativistic effects in a single trapped ion [9], that was subsequently experimentally realized [10]. This was the first experimental observation of the physics associated with the Zitterbewegung phenomenon, the fast quivering motion stemming from the Dirac equation and predicted by Schrödinger in the early days of quantum mechanics [11]. Further developments produced the theoretical proposal [12] and experimental realization [13] of the simulation of Klein paradox [14], and a proposal for simulating the Majorana equation and unphysical operations like charge conjugation or time reversal with trapped ions [15].
Many other proposals and some experiments, in a wide variety of systems, have recently appeared, as the simulation of Zitterbewegung in semiconductor quantum wells [16] or in graphene [17,18], Klein paradox in graphene [19], Dirac oscillator in a trapped ion [20], Zitterbewegung and Dirac physics with ultracold atoms [21,22,23,24,25,26], Klein paradox with atomic ensembles [27], optical Zitterbewegung in metamaterials [28,29], delocalization of relativistic Dirac particles in cold atoms [30], photon wave function and Zitterbewegung [31], similarity of electron's Zitterbewegung to the Adler-Bell-Jackiw anomaly in QED and its manifestation in graphene [32], photonic analog of Zitterbewegung in binary waveguide arrays [33], Zitterbewegung theory in multiband Hamiltonians [34], classical Zitterbewegung in reduced plasma dynamics [35], Zitterbewegung analogs in nonlinear frequency conversion [36], experimental realization of an optical analog for relativistic quantum mechanics in an optical superlattice [37], relation between parity and Zitterbewegung and proposed simulation in trapped ions [38], Zitterbewegung in a magnetic field and proposal for trapped ion simulation [39], Wilson fermions and axion electrodynamics in optical lattices [40], the Schwinger effect for a possible implementation with atoms in optical lattices [41], Dirac equation for cold atoms in artificial curved spacetimes [42], a theoretical analysis of cold atom simulation of interacting relativistic quantum field theories [43], or an analysis of the photonic simulation of the quark model [44]. For a review of Zitterbewegung of electrons in semiconductors, see Ref. [45].
In this article, we make a step forward and analyze novel features of relativistic quantum mechanics simulations with trapped ions. We first analyze Klein paradox in 2+1 dimensions for different kinds of potentials. We address the problem of 2+1 Klein tunneling which contains novel and interesting features, like entanglement between the transmitted/reflected wave packets in the Klein dynamics and the transverse momentum wave function. We also point to the fact that more general nonlinear V (x, y) potentials will already demand computational classical resources that are approaching what is currently feasible. Thus, a general 2+1 Klein dynamics will already be an interesting problem to be addressed by a quantum simulator. Then, we analyze a Dirac particle in a magnetic field, making here an emphasis on its topological properties, and show how it can be simulated with trapped ions. Finally, we study the quantum simulation of two Dirac particles coupled by a confining potential implemented in a system of three trapped ions. This model represents a simplified version of the MIT bag model of nuclear physics [46,47]. This semianalytic approach is still frequently used for analyzing quantum chromodynamics (QCD) in the nonperturbative regime. In this manner, we show that one could implement interesting QCD features like asymptotic freedom and confinement in a table-top experiment, controlled just by turning on and off a laser.
The paper is organized as follows. In Section 2, we briefly review the quantum simulation of the Dirac equation in trapped ions. In Section 3, we analyze the Klein paradox in 2+1 dimensions. In Section 4, we study the simulation of a Dirac particle in a magnetic potential and its topological properties. In Section 5, we address the problem of two Dirac particles coupled via a confining potential simulated with three trapped ions. Finally, we present our concluding remarks in Section 6.

Dirac equation simulation in trapped ions
The Dirac equation is arguably the most fundamental wave equation [48]. It is historically considered as a significant step forward towards unifying quantum mechanics and special relativity, accurately describing the hydrogen atom spectrum, while predicting ab initio spin and antimatter. The Dirac equation acquires full significance in the context of quantum field theory and second quantization, where the number of quanta is not fixed. On the other hand, at the level of relativistic quantum mechanics, the Dirac single-particle solutions predict already intriguing phenomena. Among them, the best known are the Zitterbewegung [11] and Klein paradox [14].
The 3+1 Dirac equation reads where α and β are the Dirac matrices that obey the Clifford algebra structure, {α i , α j } = 2δ ij , {α i , β} = 0, β 2 = I 4 , that appear when linearizing the expression for the relativistic energy, E = p 2 c 2 + m 2 c 4 . The Hamiltonian operator in Eq. (1) expressed in its "supersymmetric" representation [48], takes the form where α := (α x , α y , α z ) = off-diag( σ, σ) is the velocity operator, σ := (σ x , σ y , σ z ) are the Pauli matrices, β := off-diag(−iI 2 , iI 2 ), c is the speed of light and mc 2 the electron rest energy. A quantum simulation of the Dirac equation is a useful playground to analyze different regimes and to verify in a table-top experiment predicted phenomena like Zitterbewegung [9,10] or Klein paradox [12,13]. To perform just the Dirac equation and Zitterbewegung, one considers a single trapped ion of mass M inside a Paul trap, with motional-mode frequencies ν x , ν y , and ν z [9]. In the 3+1 case, the fourcomponent Dirac bispinor will be codified in four internal levels of the trapped ion, |a , |b , |c , and |d .
For any spacetime dimensions, one should realize that the Dirac equation contains basically couplings σ i p i (spin-orbit), and mc 2 σ j (mass term). In general, one may implement the first kind of couplings, σ i p i , by a combination of Jaynes-Cummings (JC, also named red-sideband) and anti-Jaynes-Cummings (AJC, also named bluesideband) interactions. The JC interaction excites one quantum of vibration while deexciting the internal state of the ion. The AJC interaction, in turn, excites one quantum of vibration while exciting the internal state of the ion. These interactions may be obtained by suitably chosen lasers, either in Raman or quadrupole-transition configurations, depending on the chosen ion.
The resonant JC Hamiltonian can be written as whereΩ is the coupling strength, σ + and σ − are the raising and lowering spin-1/2 operators, and a and a † are the annihilation and creation operators associated with the corresponding motional degree of freedom, either x, y, or z. η = k /2M ν is the Lamb-Dicke parameter [49], where k is the wave number of the driving field. The AJC Hamiltonian reads By properly adjusting laser phases and frequencies, one may combine a JC and an AJC interactions to obtain the kind of couplings aimed for, namely where ∆ x := /2M ν x is the spread in position along the x-axis of the ground state harmonic-oscillator wavefunction and p x the corresponding dimensioned momentum operator. These interactions may be addressed to different internal levels and motional modes, in order to get the terms in Eq. (2) that are linear in the momenta.
For implementing the second kind of couplings, mc 2 σ j , one may consider the carrier interaction between two internal levels of the ion. This consists of a coherent excitation of the internal level, while leaving the motion unaffected. The carrier Hamiltonian is that, for appropriately chosen phases, gives rise to the term Ωσ y , needed for the mass terms for 3+1 dimensions. For the 3+1 dimension simulation [9], notice that the Hamiltonian is composed of basic JC, AJC and carrier units as exposed above. When this is expressed in matrix form, in the basis |a , |b , |c , and |d , An important property of quantum simulations of relativistic quantum mechanics with trapped ions is that the values of the speed of light and rest energy in Eqs. (9) may be controlled at will just by changing the laser intensities,Ω and Ω. Thus, one may explore the transition from massless to massive Dirac fermions, or, equivalently, the transition from ultrarelativistic to nonrelativistic physics.
One of the most astonishing predictions of the single free-particle solutions of the Dirac equation is the fast quivering motion called Zitterbewegung. It is unexpected because it predicts an oscillatory motion of a freely propagating electron. Thus, Galileo's inertia law is not fully verified for a free relativistic electron, contrary to the Schrödinger equation case, and the free Dirac particle is expected to quiver around in the absence of potentials. The reason for this is the non-commutativity of the components of the velocity operator, which is given by c α. Thus, the marriage between quantum mechanics and special relativity seems to contradict the inertia law, at least to some extent and within the subtle realm of relativistic quantum physics. The Zitterbewegung phenomenon has not been observed so far for a real relativistic electron, given that the predicted frequency, ∼ 10 21 s −1 , and amplitude, ∼ 10 −11 cm, are difficult to access experimentally. Moreover, its correct physical prediction and its validity has been constantly questioned in the last decades, whether we remain in the domain of relativistic quantum mechanics or quantum field theory.
The expression for the Dirac electron's position operator r = (x, y, z) in the Heisenberg picture, derived from d r/dt = [ r, H D ]/i , for 3+1 dimensions, reads [48] where the first term on the r.h.s. is the initial position, the second is just the inertia law term, and the third one is the term associated to the Zitterbewegung.
In the trapped ion simulation, the position of the ion mimics the position of the simulated Dirac particle. The ionic position operator evolution under the Dirac dynamics in Eq. (8), in 3+1 dimensions, is and the Zitterbewegung frequency can be estimated in the simulation as whereĒ D ≡ H D is the average energy, p 0 is the average momentum for a peaked distribution, and N ≡ a † a is the phonon number, respectively. Additionally, one can estimate the Zitterbewegung amplitude to be and R ZB ≈ ∆, if ηΩ ∼ Ω, which can be measurable in an experiment, as was shown in the 1+1 case in Ref. [10].

Bidimensional relativistic scattering simulated in trapped ions
The Klein paradox is another counterintuitive prediction of the relativistic quantum mechanical solutions of Dirac equation. A Dirac particle may tunnel through a steep barrier and propagate indefinitely regardless of the fact that the potential may be extended until arbitrarily long distances. The reason for this is that the positive energy electron may become, upon collision with the barrier, a negative energy electron and as such propagate freely throughout the barrier, wherever the potential barrier is V ≥ 2mc 2 . The standard explanation in quantum field theory is that there is enough energy in the system to produce a particle/antiparticle pair. The theoretical proposal [12] of the simulation of the Klein paradox in 1+1 dimensions was based on the Hamiltonian where α is the potential gradient constant. The experimental simulation of the Klein paradox relied on the manipulation of two trapped ions [13]. The center of mass mode together with the internal degrees of freedom of one of them simulates the positive and negative energy spinors, while the second one implements the external potential that produces the Klein paradox [12,13]. In order to measure the motional state of the ions in the previous simulations, the standard approach is to couple the motion to the internal states of one ion and, then, apply fluorescence detection to measure the internal state with novel techniques [50,51,52,9,10].
Here, we propose an extension of the Klein paradox simulation to 2+1 dimensions. We will show that, for some potentials, this problem may be addressed with the techniques already developed for 1+1 dimensions. On the other hand, for arbitrary potentials in 2+1 dimensions, one is already approaching the limit of classical computational power and, in this case, a quantum simulator could already predict physics beyond current capabilities. At the end of this section, we give a heuristic analysis of this issue.
We begin by analyzing the Dirac equation in 2+1 dimensions with a linear potential along the x coordinate and constant along the y coordinate. The expression reads where α is the potential gradient constant. In order to analyze this dynamics, one can point out that p y is a constant of motion: [H Kl,2+1 D , p y ] = 0, such that one may solve the problem for each value of p y . One can realize that cσ y p y + mc 2 σ z , with p y fixed, can be expressed as a new Pauli matrix with an effective mass coefficient, withmc 2 = p 2 y c 2 + m 2 c 4 . Hereσỹ = n y σ y + n z σ z is an effective Pauli matrix in the y − z plane (with a basis change in that plane), with (n y , n z ) = (p y /mc, m/m). Accordingly, the problem, for each p y , is reduced to a Klein Hamiltonian in 1+1 dimensions, In Fig. 1, we illustrate the dependence of the Klein tunneling on the incoming momentum orthogonally to the potential gradient, p y . To do so, we simulate a wavepacket which initially has a Gaussian profile on the x direction and a well defined momentum p y . We plot the resulting wavefunction (normalized just in x and constant in p y for the sake of comparison between different p y components), |ψ(x, p y , t)| 2 , for different time snapshots t = 0, 20, 40, 60, 80, 100, in units of = c = α = 1, m = 0.5. These plots evidence the rapid suppression of Klein tunneling for increasing p y . Given that the transmitted amplitude decreases exponentially with the effective massm as exp(−πm 2 c 4 / cα) [12], we find that already at p y = 1 almost all the wavefunction gets reflected.
The plots in Fig. 1 admit another interpretation, based on the entanglement between the transmitted and reflected wave packets in the x coordinate and the transverse momentum along the y direction. In the language of quantum field theory, the electron-positron pair production along x direction will be conditional to the specific value of each momentum wave packet component along y direction. Thus, one could control pair production just by changing incidence angle, p y /p x , and even perform quantum logic on pair creation conditional to the value of the transverse momentum.
Based on the behaviour at fixed values of p y we can reconstruct the dynamics of arbitrary wavepackets, fully in position space. More precisely, we can write where ψ 0 (x, p y , s) is the initial state of the wavepacket, expressed in position space for the X coordinates, |x , momentum for the Y coordinate, |p y , and spinor state, |s = 0, 1 . This reconstruction formula evidences the conditional evolution mentioned before, and in particular the entanglement between the x coordinate transmitted and reflected wave packets, and the y component wave packet.
In Figs. 2(a-c), we illustrate solutions obtained using the reconstruction procedure (18), plotting the wavefunction on position space, |ψ(x, y, t)| 2 , at three instants of time t = 0, 60, 100. In addition, to ease comparison, Fig. 2d shows a plot that combines all three density distributions. The most obvious effect is the squeezing of the transmitted and the broadening of the reflected wave functions after the collision with the potential barrier. To understand this phenomenon we must realize that the energy band curvature induced bym(p y ) is associated to a momentum-dependent group velocity which in turn leads to a spreading of wavepackets in position space. This spreading is most relevant for small values of p x , that is when the particle crashes against the potential barrier. At this point, which is when the particle splits into a transmitted and a reflected component, a narrow band of p x components will be filtered and become part of the antiparticle branch, moving on with very little spreading and almost uniform group velocity. The reflected component, on the other hand, will get to see the whole curvature of the energy band and become very broad very quickly, until it reaches relativistic velocities (v g c) and stops spreading. For performing the protocol, we envision to follow the lines developed for the 1+1 case [12,13] and use two trapped ions. One of them, together with two motional modes, will codify the spinor, including the x and y dependence of the wave packet. The other one will produce the linear potential in x. To be more specific, we consider two trapped ions trapped in a linear Paul trap, with two motional modes, with frequencies ν cm = ν and ν st = √ 3ν. To implement the Dirac-Klein dynamics given by Eq. (15), we proceed in an analogous way as in Refs. [9,10,12,13] and use red and blue sidebands for each of the terms cσ i p i , i = x, y, as well as an AC-Stark shift for term mc 2 σ z . The αx potential term will be implemented with a red and blue sidebands, with appropriate phases, applied to the second ion [12,13]. The corresponding trapped-ion Hamiltonian will read which coincides with Eq. (15) when one performs the analogies 2η cm ∆ cmΩcm = 2η st ∆ stΩst ↔ c, Ω ↔ mc 2 , and η cmΩ0 /∆ cm ↔ α, when ion 2 is initialized in the eigenstate + of σ x,2 .
For measuring the wave packet after performing the quantum simulation, one would proceed by mapping the motional state to the internal degrees of freedom of ion 2, and subsequently detecting the internal state by fluorescence detection, as exposed in Refs. [50,51,52,9,10]. In these references, it was shown how to perform tomography of the ion motional state in an efficient way. We point out that, when considering a nontrivial potential V (x, y) in x and y coordinates, say V (x, y) = α x x 2 + α y y 2 , which could be implemented with three ions in the dispersive limit, its simplification into a set of 1+1 problems cannot be done anymore, and the full bidimensional problem has to be considered. This means that, for reasonable sizes of motional Hilbert spaces of around 200 phonons per mode, one would have a Hilbert space dimension of 2 × 200 for the 1+1 dimensional case. On the other hand, for the 2+1 case considered here, one would have a Hilbert space dimension of 8 × 10 4 . The size of the states and the size of the associated Hamiltonian matrix force us to adopt particularly efficient schemes in the simulation, such as using finite differences or Fourier representations for the Hamiltonian and Trotter methods for the unitary evolution. However, none of these techniques are without errors. They are severely affected by discretizations of time, and also by the boundary conditions, and even the simulations presented here are well within the limits of what can be accurately simulated, both in time and in space. These are the reasons why an actual experimental simulation of the 2+1 Klein scattering may have such an important value: without these constraints -memory, computational power, narrow boundary conditions-, it allows for an independent and perhaps more reliable verification of our predictions.

Simulation of Dirac particles in magnetic fields
So far, we have seen examples in which the trapped-ion simulation may be used to learn more about the Dirac equation, but the converse is also true: working in the quantum simulations may offer a new perspective on quantum optical problems. In particular, in this section, we will relate the simulation of Dirac particles in magnetic fields with the Jaynes-Cummings model, and use this to shed light on the topological features of the eigenstates of both models -the simulated one and the one used to implement the simulation.
The topological properties of free Dirac fermions has become a very active research area in the last years [53]. Roughly, in this context the energy band of a 2D free Dirac fermion can be seen as a mapping from momenta, (k x , k y ) to a particular spinor state on the sphere, S ∈ S, obtained from the Dirac wavefunction. The topologically nontrivial models are those for which this mapping winds up one or more times on the same sphere. The intriguing result is that such models, when embedded on surfaces with boundaries, present robust topologically protected "edge states", that is states that may transport charge or spin even when the bulk is placed in an insulating state (for instance with a large mass).
We will see that some of these ideas may be cast in the context of trapped-ion simulations. For that, let us recall the model of a Dirac particle moving on a uniform two-dimensional magnetic field [20,54], We may choose an axial gauge in which A(x, y) = Bx. In this setup, the momentum along the y direction becomes a constant of motion, [p y , H] = 0. Changing units so that eB = 1, and displacing the "x" coordinate by an appropriate amount, x + p y → x, we may rewrite the previous model as a simple Jaynes-Cummings model with a detuning, This model has a ladder of discrete energy eigenstates, the so called Landau levels. The wavefunctions of these levels are confined along "x" and form plane waves along the "y" direction where φ n (x) is the n-th level of a quantum harmonic oscillator and the coefficients (α ±n , β ±n ) are the corresponding eigenvectors of the matrix At this point, instead of focusing as usual on the topology of energy bands, we will try to understand the topological properties of these individual states. For that we introduce the Wigner function of the spinor, given by This quasiprobability distribution may be interpreted as a mapping from phase space ( x, p) to points on the Bloch sphere, S ∈ S, where N = Tr{W σ} is a normalization factor and we have dropped the y coordinate. Figure 3 shows the Wigner function of the two simplest Landau levels, n = ±1 and n = ±2, decomposed into two quantities, the polar angle φ = tan −1 (s y /s x ), and the separation from the equator, s z . For the n = 0 state we see that the pseudospin points North in the center of the phase space and rotates to the South at r = x 2 + p 2 → ∞. This is equivalent to a single covering of the sphere, or skyrmion, with a topological charge ‡ ν = 1. To achieve greater complexity we have to increase the number of quanta, moving on to state n = 2, which performs a double covering of the sphere. This is done via an extended singularity of the angle φ at radius r = 1, which is when the spin points North. In general we have verified that the number of coverings increases linearly with the number of quanta, ν = ±n, leading to a family of topological defects with unprecedented richness. It may seem that these topological defects are mere artifacts of the Dirac or Jaynes-Cummings equations and that under normal circumstances, dephasing, spontaneous decay or other perturbations, will rapidly disappear. However, this is not the case. Pure dephasing, for instance along the "z" direction, is equivalent to a shrinking of the transverse components s(t) = (s x e −γt , s y e −γt , s z ). (28) This map contracts the total surface in Bloch space, which now looks closer to a rugby ball instead of a sphere. However, after suitable rescaling via a new normalization factor in Eq. (27), one finds that the total solid angle covered is the same and the charge, ν, is not modified.
It is even more surprising to see that spontaneous emission has also a mild effect on the winding number. Applying the positive map that corresponds to this process onto the ion state, we see that the winding on the {S x , S y } plane remains unaffected by the same reasons as before, and the only change is a shift of the weights of the different harmonic oscillator eigenstates, or equivalently s z . In Fig. 3d, we summarize that process for the n = 2 Landau level, showing how dissipation tries to "unwrap" the sphere, removing one layer and falling down to a lower winding number. Note how, despite the effort in unwrapping the sphere, the defects remained pinned and the total number constant even for long times. One may consider a final form of losses: total dissipation in the form of friction in "x" or "p". In either case, we expect population being transferred down the ladder of harmonic oscillator states, which this time would lead to a continuous unwrapping of the sphere ν = n, n − 1, n − 2 . . . 0.
Regarding the experimental implementation, the Hamiltonian for the Dirac particle in the field has, as we have shown, a direct reinterpretation in terms of a detuned Jaynes-Cummings model (22), which makes its simulation rather straightforward. Measurement of the Wigner function is also possible using standard trapped ion techniques [49], but now the reconstruction has to be performed without tracing out the internal state of the ion. To do this will require using an ancillary ion for measurement and a different ion to implement the spinor. The ancillary ion will remain in the ground state until the end of the experiment, when it will be used to gather statistics about the vibrational mode that implements the Jaynes-Cummings Hamiltonian [49]. The novelty is that this statistics will have to be performed conditionally on the spinor-ion, which should be simultaneously measured in one of the σ x , σ y or σ z basis, as it is usual for wavefunction reconstruction. The resulting statistics and accuracy should be more than enough to allow for accurate reconstruction of the ν = 1 and ν = 2 skyrmions.

Two Dirac particles interacting via a classical potential: Simulation of an analogue of MIT bag model
The quantum field theory for the strong interaction, known as quantum chromodynamics (QCD), accurately predicts the behavior of nuclei, hadrons, and quarks and gluons. On the other hand, for low energies and long distances the theory becomes non-perturbative (because the coupling becomes too large) such that Feynman-diagram expansion cannot be applied.
Many decades ago, some effective theories were developed that qualitatively and quantitatively described the physics of the strong interaction even in the nonperturbative regime. Among those theories, arguably, the best known is the MIT bag model [46,47]. In its simplest version, it consists of two or more Dirac particles contained inside a volume, a bag, whose energy grows linearly with the volume. In this way, this model can simulate behavior like asymptotic freedom and confinement: when the fermions go far away from each other, the energy grows steadily and there is a tendency to come back to the original position, i.e., confinement takes place. On the other hand, when the fermions are nearby, the potential energy becomes negligible and the fermions behave as if they were free, i.e., asymptotic freedom holds. Although these models can be solved in a classical computer, and are phenomenological, they have quite a good agreement with experiments and predictive power. Even nowadays, this kind of methods, either with color charge or without it, string models, etc., are the basic tools to analyze QCD at low energies in the nonperturbative regime, with semianalytic tools (without resorting to lattice gauge theory, which requires huge computational power) [55]. We believe it would be interesting to implement analogues of the MIT bag model in a quantum simulator, even if it is just for reproducing the physics of quark confinement and asymptotic freedom in a table top tunable experiment. In addition, quantum simulations of these models with state-of-theart technology could already go beyond what can be computed classically (for large motional Hilbert space dimensions). In this Section, we propose the simulation of a simplified analogue of MIT bag model in 1+1 dimensions. Thus, we propose a simulation of two Dirac equations, one for each Dirac particle, that are coupled by a potential which grows with the distance, in our case quadratically for the sake of experimental feasibility. This model could describe a quark and antiquark coupled by gluons inside a meson, and cannot be solved analytically. By appropriately tuning the lasers in an experiment, one could perform phase transitions between the asymptotically free and the confined meson phases, like the one that supposedly took place in the early Universe.
The model we aim to simulate is In order to simulate two Dirac particles (e.g., quark-antiquark constituting a meson), each obeying a 1+1 Dirac equation and with a potential that grows with the separation between the particles, we envision to use three ions: two of them will codify the positive and negative energy spinors of the two Dirac particles, while the third one will be used to generate the attractive potential among them. The normal modes for a three ion system with equal mass and charge are We will consider the two first modes to codify the two free Dirac equations plus the potential. By considering the Hamiltonian we simulate an equivalent dynamics to two free Dirac equations with a potential in the relative coordinate, wherep r = (p 1 −p 3 )/2, andP cm = p 1 +p 3 are the relative and center of mass momenta of the system of two particles that we want to simulate, andx r = x 1 − x 3 is the relative coordinate. To make the analogy complete, one has to make the substitutions P cm ↔p r , p r ↔P cm /2, Q cm ↔x r , with 2η cm ∆ cmΩcm = 2η r ∆ rΩr = c, Ω = mc 2 , and the potential is obtained in the dispersive limit, for ∆ 3 η cm Ω 3 a † cm + a cm , considering ion 2 in the + eigenstate of σ z,2 , and with V 0 = ( η cm Ω 3 /∆ cm ) 2 / ∆ 3 . With this simulation, we may analyze confinement and asymptotic freedom by playing with V 0 potential. One would expect in principle, though predictions are complicated in systems with no analytical solutions, to have two independent dynamics. One of them will happen for the case in which the fermions are nearby (x r 0) and the momentum |p r | is small (asymptotic freedom). And another more complicated dynamics, coupling the two fermionic motions, for those cases in which the two fermions fly away (x r large) (confinement). Additionally, when V 0 x 2 r ≥ 2mc 2 , one would expect Klein tunneling to take place. In this model, we interpret this Klein tunneling as a pair creation of a new quark and a new antiquark: when the original quark and antiquark go very far away from each other, there is energy available for creating a new antiquark (that gets bound to the original quark) and a new quark (that gets bound to the original antiquark), such that the original pair can split up without violating the color neutrality standard in QCD for free composite particles. This is the natural interpretation, given that the usual explanation for Klein tunneling in quantum field theory is the creation of a particle/antiparticle pair. In QCD, the splitting of a hadron into new hadrons is well known and takes place, for example, in the jet emission in high-energy colliders.
We have realized numerical simulations of this model, and we show our results in Fig. 4. There, we plot the probability |ψ(x r , t)| 2 as a function ofx r and t, for the two Dirac fermion system as given by Eq. (34), forP cm = 2 (a constant of motion), mc 2 = 1, and V 0 = 0.5, for the initial values (a) p r = 0, Π = +1, (b) p r = 2, Π = +1, (c) p r = 0, Π = −1, and (d) p r = 2, Π = −1. In this notation, Π denotes the expected value of the operatorΠ = 1 2 (σ x,1 − σ x,3 ) on the initial state. We consider initial Gaussian states inx r , normalized in this coordinate, and constant P cm , that will remain unchanged during the evolution, for the sake of simplicity and ease of the numerical simulation.
Despite the fact that this problem is not analytically solvable, and it is difficult to get a clear intuitive prediction of the behavior that Eq. (34) will produce, we observe interesting features. For Π = +1, and p r = 0 the dynamics is quasifree, besides a small Klein tunneling appearance. This may be a signature of asymptotic freedom, because the Dirac fermions do not couple much, as seen in (a). For Π = +1 and p r = 2, the dynamics of the two-fermion system is more involved: the wave function evolves in a more complex way for larger relative momentum, splitting up, bouncing back and getting distorted. This may be related to confinement taking place, when the two fermions try to separate, as seen in (b). On the other hand, for Π = −1 (another spinor initial state), we observe that, for finite relative momentum, p r = 2 (d), the two fermions have a larger Klein tunneling than for zero relative momentum, p r = 0 (c). This may also indicate that, for large enough relative momentum, the Dirac fermion pair has enough energy to create another quark-antiquark pair and escape from each other while keeping color neutrality (d). For relative momentum equal to zero, there is little Klein tunneling (c) associated to asymptotic freedom. Summing up, we have shown than an experiment with three ions can simulate an analogue of the well-known MIT bag model [46,47] for two Dirac particles in 1+1 dimensions. Despite their simplicity, such experiments would allow us to study signatures of asymptotic freedom and confinement. In addition, those experiments would easily reach simulation regimes which are beyond current numerical simulation possibilities. For example, here, we were severely limited by the problem size and complexity, preventing us from studying things like adiabatic preparation of eigenstates, a full diagonalization of the spectrum, or longer time dynamics.

Conclusions
In this article, we have shown that trapped-ion physics provides a flexible platform to simulate a number of interesting cases and effects in relativistic quantum mechanics. In particular, we have studied bidimensional relativistic scattering, topological effects of Dirac particles in magnetic fields, and two Dirac particles coupled by a potential, an analogue of nuclear physics bag models. Summarizing, we have described the bidimensional relativistic scattering for an x-dependent potential, showing that this case can be solved with the help of the one-dimensional case. In addition, we showed it contains interesting novel features, like entanglement between transmitted/reflected wave packets and the transverse momentum. We have also given a heuristic analysis of the complexity of more complicated potentials that will already approach the classical computational limit of resources. We have introduced a relativistic quantum mechanics model -a Dirac particle in a magnetic field-that could shed light and establish interesting analogies with an equivalent quantum optical model, namely, the Jaynes-Cummings model with detuning. We showed a possible quantum simulation of this case, including its topological features in trapped ions. Finally, we have shown that three trapped ions together with two motional modes suffice to implement the dynamics of two Dirac particles that are coupled by a confining potential, like in the bag models of nuclear physics. This could allow the implementation in a table top experiment of quark confinement and asymptotic freedom associated to QCD, while moving from one regime to the other by switching on and off a laser. This would amount to a phenomenological simulation of the phase transition between the quark-gluon plasma and the hadronic confined phases that took place in the early universe.
In a physical implementation, strings of Ca + ions, for instance, could be used and manipulated with laser light. The spinor degrees of freedom can be encoded in long-lived electronic states of this ion, as was done before in quantum simulations of relativistic quantum mechanics [10,13]. Such states can have coherence times of several ms, while coupling strengths of a few 100 kHz between the spinor states can be obtained. For each motional degree of freedom, a normal motional mode of the ion string can be used and manipulated with lasers coupling the internal and motional degrees of freedom. The implementation of the proposed quantum simulations requires the validity of the Lamb-Dicke approximation, which for Ca + ions in a trap with secular frequency of ∼ 1 MHz and manipulated with 729 nm laser light is well justified (Lamb-Dicke parameter η ∼ 0.05 1) [13] . In recent experiments [13], it has been shown that the motional state of the ions can be coherently manipulated while reaching states that have >100 phonons on average, and the coherence can be sustained for more than 10 ms. These parameters could be further improved by using, e.g., ions that are less sensitive to decoherence sources and by further decreasing the Lamb-Dicke parameter with the choice of a different beam geometry.
We believe that quantum simulation is one of the most promising fields inside quantum information science. At the same time, among other research lines, quantum simulations of relativistic quantum mechanics may soon allow us to consider problems that are difficult or even impossible for classical computers. In this sense, quantum simulators will give us the opportunity to enter into unexplored regimes of the physical world. Trapped ions offer one of the most promising platforms for achieving this goal.