Quantum Simulations of Relativistic Quantum Physics in Circuit QED

We present a scheme for simulating relativistic quantum physics in circuit quantum electrodynamics. By using three classical microwave drives, we show that a superconducting qubit strongly-coupled to a resonator field mode can be used to simulate the dynamics of the Dirac equation and Klein paradox in all regimes. Using the same setup we also propose the implementation of the Foldy-Wouthuysen canonical transformation, after which the time derivative of the position operator becomes a constant of the motion.


Introduction
R. P. Feynman is credited for introducing the idea of a quantum simulator [1]. At his keynote speech during the 1st Conference on Physics and Computers in 1981, he claimed that because the laws of Nature are not those of classical mechanics, it would be better to build a quantum simulation of a physical problem under the same quantum mechanical laws. For the following ten years, his contribution was somewhat hindered by the work of Deutsch [2], who showed the existence of a universal quantum computer, a quantum mechanical generalization of the classical Turing machine using quantum logic gates and algorithms.
The concept of quantum simulators will reappear later in 1996, when Lloyd presented a general proof for Feynman's conjecture [3]. Quantum simulators can surpass the capabilities of classical ones, where the complexity of a problem scales exponentially. Lloyd argues [3] that whereas the simulation of a general 40 spin−1/2 particle system would be enough to outperform existing classical computers, the factorization of a 100digit number with Shor's algorithm would need of thousands of qubits to become nontrivial. Therefore, it is likely that the first quantum device to beat classical computers would be a quantum simulator. But there is another reason why quantum simulators are attracting so much attention lately. They provide powerful analogies to understand the physics of most varied type of systems one can think of, including condensed-matter physics, high-energy physics, cosmology and quantum chemistry [4,5]. When compared to the original systems, quantum simulators show an unprecedented degree of controllability over all physical parameter regimes.
In this work, we show how the physics of relativistic quantum mechanics can be simulated in a very different setting. That of circuit quantum electrodynamics (cQED) and superconducting qubits. We will focus our attention on the fundamental equation in relativistic quantum physics, the Dirac equation. Proposed by Dirac himself in an attempt to combine quantum mechanics and special relativity in 1928 [6], this equation is the first one to describe elementary spin−1/2 particles-such as the electron-and the existence of antimatter. Here, we restrict our discussion to the simplest case of 1+1 dimensions for which the algebra of the Dirac spinors can be reduced to that of standard Pauli matrices [7]. In the standard representation [7], the 1+1 Dirac equation can be written as where σ y and σ z are Pauli matrices,P represents the momentum of a Dirac particle of mass m and c the speed of light. Soon after its introduction, Schrödinger realised that there is something odd in the solution of this equation [8]. The position for a free Dirac particle can be integrated to give [9] x(t) =x 0 + c 2P H −1 t +Ẑ(t), with H = mc 2 σ z + cP σ y and This third term appears because the time derivative of the position operator is not a constant of motion, i.e. [ẋ, H] = 0, even for a free particle. Given a general initial wavepacket, the expectation value Ẑ (t) shows an oscillating pattern with amplitude ∼ /mc and frequency ∼ 2mc 2 / that is known as Zitterbewegung or jittering motion [9]. This behaviour stems from the interference between positive and negative energy components. Indeed, it is well established that initial states with only positive or negative components will not show any Zitterbewegung [7,9]. In spite of the appealing physics behind the Dirac equation, there has never been an experiment to test the existence of such an effect. In fact, for relativistic electrons, one can estimate a Zitterbewegung with an amplitude of about 10 −4 nm and a frequency of 10 21 Hz, making its eventual detection very demanding.
There is another interesting property of the Dirac equation related to the behaviour of relativistic Dirac particles under the effect of an external scalar potential, In case of a nonrelativistic particle impinging against a semi-infinite step potential Φ(X), we know that if the kinetic energy of the particle is smaller than the potential height, no propagating solutions exist at the other side of the potential barrier. In other words, the particle penetrates slightly, but it will be eventually reflected. Klein found out that Eq. (4) for a relativistic Dirac particle allows for propagating solutions beyond the potential barrier [10]. And in the case of an infinitely high potential, the tunneling occurs with unit probability. More precisely, while for a small potential Φ < E−mc 2 the particle with kinetic energy E will be mostly transmitted, for a larger one E−mc 2 < Φ < E+mc 2 it will be reflected. However, making even larger the potential Φ > E + mc 2 the particle will propagate through the barrier. To explain this apparent paradox, one can make use of the idea of particle-antiparticle creation by the potential [7]. A particle can go through the potential by turning into its antiparticle. Although the Dirac equation is considered a milestone of relativistic quantum physics, playing a crucial role towards the later development of quantum field theories, in the last decades it has driven a renewed interest in different playgrounds. Examples include: superconductivity, where Bogoliubov quasi-particles can exhibit Zitterbewegung [11]; semiconductor physics, with a close analogy between the k · p theory of energy bands in narrow-gap semiconductors and the Dirac equation [12]; the physics of graphene, characterised by quasi-electrons behaving like relativistic Dirac particles in two dimensions [13]; a proposal for generating relativistic scattering in optical latices [14]; photonic arrays with implementation of the quantum Rabi model and the Dirac equation [15,16]; and the physics of ion traps, where the vibrational and internal degrees of freedom of the ions can be made to follow relativistic quantum dynamics [17,18,19,20,21], or even behave as interacting bosons and fermions [22,23]. In the spirit of quantum simulations, the latter represents an important advancement in terms of full quantum controllability of the dynamics and the retrieval of information through measurement.
Here, we present a scheme for performing quantum simulations of relativistic quantum mechanics using state-of-the-art cQED technology. This is motivated by two reasons. First, since 2004 [24,25], cQED has become the quantum platform with the most promising perspectives in terms of scalability and coherence. Second, there are crucial physical differences with respect to any previous implementation of the Dirac equation and Klein paradox. In our case, the physics of a relativistic spin−1/2 particle is simulated by the interaction of two degrees of freedom of different physical systems, i.e., a standing wave in a superconducting resonator and the two lowest levels of a superconducting qubit, where none of them move at all. The position and momentum of the Dirac particle are then encoded in its phase-space or field quadrature representation. This opens up new possibilities for combining intracavity fields with propagating quantum microwaves [26,27,28] in scalabale quantum network architectures [29] with delocalised and/or sequential interactions.

General derivation
For our protocol, we require a two-level superconducting qubit, as a flux qubit [30], working at its symmetry point, strongly coupled to a single electromagnetic field mode of the resonator. This interaction will be described by the Jaynes-Cummings model (JCM) [31,24,25]. In addition, we introduce three external classical microwave drivings, two transversal to the resonator [32], coupling only to the qubit, and one longitudinal that only sees the resonator mode. The time dependent Hamiltonian of the complete system is given by with σ y = iσ − iσ † = i |g e| − i |e g| and σ z = |e e| − |g g|, |g and |e being the ground and excited states of the qubit. Here, ω and ω q stand for the photon energy and qubit energy splitting, whereas g is the coupling constant. Likewise, the two orthogonal drivings have real amplitude Ω, λ, phase ϕ, and frequencies ω, ν. Finally, the longitudinal driving is characterised by an amplitude ξ and a frequency ω. Note that while we have chosen two of the drivings to be resonant with the resonator field mode, the other parameters will be set later on. For simplicity in our presentation, we will also assume that ω q = ω, i.e. qubit and resonator field interact resonantly as well.
Our derivation consists of two straightforward transformations. First, Hamiltonian in Eq. (5) can be simplified using the reference frame rotating with the resonator frequency ω Second, we will write this Hamiltonian in another interaction picture with respect to the term H L 1 0 = − Ω e iϕ σ + e −iϕ σ † . The resulting expression reads where we have introduced the rotated spin basis |± = (|g ± e −iϕ |e ) / √ 2. To simplify this Hamiltonian expression further, we can now set ω − ν = 2Ω, and assume the amplitude of the first driving Ω to be large as compared to the rest of frequencies in (7). This implies we can apply the rotating-wave approximation (RWA), yielding the effective Hamiltonian where ϕ = π/2 and we make use of standard quadratures of the electromagnetic field, at variance withX andP in Eq. (4), with the commutation relation [x,p] = i. Note that Ω does not appear in the effective Hamiltonian. This is because the Hamiltonian was obtained in an interaction picture with Ω playing the role of a dominant frequency in the strong driving regime.
The Schrödinger equation of our system resembles that of the 1+1 Dirac equation as in Eqs. (1) and (4), where the terms g/ √ 2 and λ/2 are related to the speed of light and the mass, respectively. In addition, we have an external linear potential Φ = ξ √ 2x, that depends linearly on the position of the particle. Clearly, in the simulated dynamics, one can cover a wide range of physical parameters. While, for a fixed coupling strength g, the simulated mass is proportional to the amplitude of the weak orthogonal driving λ, the strength of the linear potential can be tuned through the amplitude ξ of the longitudinal driving. This is an interesting advantage over the implementation in ion traps, where a second ion is needed to simulate the external potential [19,20]. Note that in the case of a massless particle, λ = 0 and ν = 0, so that ω = 2Ω in Eq. (7).
As already mentioned, the study of relativistic quantum effects, such as Zitterbewegung or Klein tunneling, should be done in the phase-space representation of the electromagnetic field in the superconducting resonator. In all numerical cases studied below, the initial state of the bosonic degree of freedom of the Dirac particle is assumed to be a wavepacket with average position x 0 and average momentum p 0 , This coincides with the x−quadrature representation of a coherent state of the electromagnetic field x 0 +i p 0 |0 , |0 being the vacuum state of the bosonic field, and D(α) = exp αâ † − α * â the coherent displacement operator.
For the numerical simulations discussed below, we have chosen a value of the superconducting qubit and resonator frequencies of ω q = ω = 2π × 9 GHz, a qubit-field coupling of g = 2π × 10 MHz and a strong transverse driving amplitude of Ω = 2π × 200 MHz. In addition to being experimentally realistic, this value takes into account the limitation imposed by the qubit-resonator coupling constant g on the strong driving fields. The interaction time needed is of about 60 nsec, which is well below standard decoherence times in cQED experiments.

Free Dirac particle
In the absence of external potential, Φ = 0, the Dirac equation for a free particle (1) does not couple the different spinor components. The Dirac Hamiltonian has two important limits depending on the value of the mass of the particle. First, in the massless case the Hamiltonian reduces to H D = g/ √ 2σ yp allowing for straightforward interpretation of its dynamics in terms of phase-space representation. Fig. 1(a) shows the evolution of a massless particle with initial state |+, 0 . Being its spin in an eigenstate of σ y , the effect that the interaction has is just to generate a coherent field displacement from its original position, the vacuum state |0 , to the final one in which the state of the field is |gt/2 . This simulates the evolution of a massless free Dirac particle moving at a speed of light stemming from g/ √ 2. As we expect, the dynamics of this massless particle would be essentially unchanged if one assumes the initial state to have a nonzero value of the expectation value of momentum. This is depicted in Fig. 1(b), where the particle starts its trajectory with p 0 = 2. The particle moves exactly the same amount of space in the x−quadrature, ending in the coherent state gt/2 + √ 2i . Second, when the mass of the particle becomes large enough, the dynamics should tend towards that of a non-relativistic one evolving under a free Schrödinger equation. For λ/2 ≫ g/ √ 2 p , i.e. mc 2 ≫ c P , one can derive a second-order effective Hamiltonian [17] for Eq. (8) to yield H NRel = σ z (g 2 /λ)p 2 . This is indeed the nonrelativistic Schrödinger Hamiltonian we expected, but with the addition of the spin operator σ z , which is reminiscent of the spinor-momentum coupling of the original Dirac equation. This spin operator ensures the existence of two decoupled equations, valid for positive and negative kinetic energies. Note that, for suitable initial states, entanglement between the spinorial and kinematic degrees of freedom is always possible, irrespective of the system being in the nonrelativistic or relativistic regimes. In Fig. 1(e), we depict the evolution of the initial state |e, 0 under the time dependent Hamiltonian (6) for the case in which the mass of the particle is such that λ/2 = 4 × g/ √ 2. This generates a squeezed vacuum state with degree of squeezing increasing linearly in time. This behaviour is in agreement with the aforementioned nonrelativistic approximation for the Hamiltonian, proportional top 2 . Such a simulation realizes a standard example in all quantum mechanics textbooks. Namely, the free evolution under the nonrelativistic  Schrödinger equation of an initial wavepacket with p 0 = 0 produces a particle that remains centered at the same initial position, with a wavepacket spreading over time in the x−quadrature. However, the effect of the operatorp 2 is not simply squeezing. This can be clearly seen in Fig. 1(f) where we have taken a particle with initial state e, √ 2i , that is a wavepacket with non-zero kinetic energy such that p 0 = 2. Now the wavepacket not only spreads over time, but its centre moves linearly in time to the right as the expectation value of the position of a nonrelativistic particle would do.
While for the limits of zero and large mass, the field state remains Gaussian during the evolution, this is no longer true when the particle mass is such that mc 2 becomes comparable to c P . This is the case of Figs. 1(c) and (d), where λ/2 = g/ √ 2. Starting from initial states |+, 0 and +, √ 2i , respectively, the Gaussian states of the field evolve into a nonclassical ones with negative Wigner function through a highly nonlinear interaction. The effect is more evident for the case of an initial state |+, 0 in Fig. 1(c). Initial states with non-zero kinetic energy, such as +, √ 2i in Fig. 1(d) makes the qubit-field coupling term in Eq. (8) dominant. This interesting finding shows that a Dirac evolution generates nonclassical features mainly in intermediate relativistic regimes, assuming that one traces out over the qubit degree of freedom.
With our scheme we can also study the appearance of Zitterbewegung for a free spin−1/2 Dirac particle. In Fig. 2, we plot the time evolution of x for a particle prepared in an initial state |+, 0 and with three different values of the mass. For better comparison, plots (a) and (b) correspond to calculations done using the time dependent Hamiltonian (6) and the effective one (8). The good agreement found indicates that choosing a value of the strong driving amplitude Ω = 2π × 200 MHz is enough for our purpose. The dashed line in Fig. 2 shows the evolution of a massless particle, λ = 0, whose expectation value of the position increases monotonically at a speed of light related to g/ √ 2. If we assume a mass such that λ/2 = g/ √ 2, then we observe that the evolution of the x followed by dotted line departs considerably from the constant rectilinear behaviour one would expect for a free particle. Likewise, for a mass even larger, choosing for instance λ/2 = 4 × g/ √ 2, the effect seen in the dash-dotted line is even more significant. Now the particle hardly moves forward. Instead, it develops an almost stationary oscillation around its initial position.
To gain in physical insight, we can derive an analytical expansion for the evolution operator of a free Dirac particle, where sinc(x) = sin(x)/x. Here, a larger initial average momentum gt p(0) /2 leads to a faster collapse of the trembling motion for a particle prepared in a wavepacket [11,12]. As already discussed, Zitterbewegung stems from the fact that the time derivative of the position operator is not a constant of the motion for the Dirac equation. This can be seen if one writes down the evolution of this operator (2). The Zitterbewegung termẐ(t) will be exactly zero if the initial state of the particle contains purely positive or negative components of the spinor. However, for a completely arbitrary initial wavepacketẐ(t) does not necessarily vanish. In addition, the preparation of a purely positive spinor for a massive Dirac particle is not straightforward. In Ref. [18], a heuristic method for preparing a state with only negative components was used. This shows that even for a massive Dirac particle, evolutions without Zitterbewegung are possible. Another alternative to study the existence/absence of Zitterbewegung would be to implement a measurement of the Foldy-Wouthuysen position operator [9,33], where∆ is the time independent operator σ x = iσ z σ y ,Ê = c 2P 2 + m 2 c 4 . The Foldy-Wouthuysen position operator in (12) has the property of being also canonically conjugate to the momentum operator. The main drawback ofX FW (t) from an experiential point of view is that is not diagonal in the coordinate space. This means that measuring such an operator might require a full quantum tomography of the evolved state of the system. Here, we follow a different route to the problem. Foldy and Wouthuysen are also credited for having introduced a unitary transformation [33] that allows for the direct diagonalization of the spin degree of freedom in the Dirac Hamiltonian. The physical relevance of this transformation is that it allows for a simple correspondence with classical mechanics [34]. Using the parameters of our simulation, this new Hamiltonian is given by whereŜ If one were able to follow the evolution of the system in the Foldy-Wouthuysen representation, but without transforming the measurement apparatus, then no Zitterbewegung would be seen. That is precisely, what we have done in the solid line of Fig. 2(b). This shows the evolution of x for a massive particle with λ/2 = 4× g/ √ 2, just like the dash-dotted one, but using the evolution operator instead of U D (t) = exp (−iH D t/ ). Now the particle remains centered at x = 0 during the evolution (solid line). More importantly, for large masses we can propose a physical implementation of the Foldy-Wouthuysen transformation using the same scheme we have developed. Clearly, when g p /(λ √ 2) is small enough, we can linearize the interaction in the unitary operatorŜ FW . This means that where the first and last step involve the realization of a Dirac equation for a massless particle with an evolution time t = λ −1 . Such an implementation would require some additional local rotations and Ramsey-like driving pulses [32]. It would be particularly straightforward if the setup allows for a variable coupling from −g to +g, being g in strong-coupling regime. This is the case for instance, of a gradiometer type of coupling seen in flux qubits [30]. The solid line in Fig. 2(a) depicts the evolution of x , computed under these assumptions, for the same particle of large mass, with λ/2 = 4 × g/ √ 2. Even for such moderate value of the mass, it is evident how oscillations seen for the Dirac evolution become significantly reduced.

Dirac particle in an external potential: Klein tunneling
The addition of an external potential has the effect of coupling the different energy components of the spinor. As before, in the evolution under an external potential we can identify two limiting cases. A massless Dirac particle is described by the Hamiltonian Figs. 3 (a) and (b) the evolution of the initial state |+, 0 and +, √ 2i , respectively, is shown. We use the same set of parameters as in Fig. 1, with the only addition of a linear potential of strength ξ = g/2. Again, starting with the qubit state, |+ , allows for a simple interpretation. The initial state of the field remains coherent while being subject to two independent displacements. Namely, as for the free-particle case, there is one in the x−quadrature proportional to g/ √ 2, and the second one induced by the potential that drags it down in the p−quadrature. Hence, the existence of such an external potential cannot alter the rectilinear movement in position representation at the speed of light. This is precisely the effect of the Klein paradox. A massless Dirac particle may propagate through the potential barrier with finite probability, regardless of the potential strength. The tunneling is accompanied by the 'rotation' of the positive energy components of the initial spinor. In our simple case, this means that the components with positive momentum components in the Gaussian  wavepacket will be pushed downwards by the potential becoming negative. And the rate at which this happens is given by ξ √ 2. This explains why the only noticeable difference between the plots in Figs. 1 and 3 is in the value of p(t) .
If the mass of the particle is large enough, then the nonrelativistic Schrödinger equation should be a good approximation. Its Hamiltonian can be written as H NRel = σ zp 2 /λ + ξ √ 2x, which guarantees that any initial Gaussian state should remain Gaussian during the evolution. This is indeed what we see in Figs. 3(e) and (f) where we have prepared initial states |e, 0 and e, √ 2i , respectively, with a mass such that λ/2 = 4 × g/ √ 2. In the first of them, the particle is being pushed backwards by the potential. The same can be said for the second. Although having an initial positive kinetic energy has allowed it to enter the barrier further, after 60 nsec it is already moving backwards.
While for these two limiting cases we see complete transmission or reflection, a particle with an intermediate mass will present only partial transmission/reflection. This is shown in Fig. 3(c) which corresponds to the initial state |+, 0 to a massive particle with λ/2 = g/ √ 2. The wave packet has now broken up into spinor components of different sign which move away from the center of the potential. If the wavepacket has some initial kinetic energy, as in Fig. 3(d) where we prepare +, √ 2i , then the particle tunnels slightly to stop and break up eventually.
As we did in Eq. (11), one can obtain an analytical expansion for the evolution operator of a Dirac particle under the effect of a linear potential. Using the formula we see that for short time intervals the evolution corresponds to that of the free Dirac particle (11), but with an initial state with rotated spin and displaced initial momentum with respect to the prepared wavepacket. The reconstruction of the Wigner function will allow to characterise what type of dynamics is being simulated. This can be done either using intracavity techniques [35] or allowing the quantum field to leak the superconducting resonator for its subsequent measurement. This will require the use of methods for reconstructing the moments of a propagating quantum microwave signal [26,27,28], already available in cQED.

Conclusions
We have shown that circuit QED provides a powerful platform for simulating the relativistic quantum physics, establishing a useful dialogue between unconnected fields. Tuning the parameters of three classical microwave drives, the system is made to evolve according to the Dirac equation for a particle in presence of a linear potential. The degree of controllability offered by this scheme would allow to study interesting physical phenomena far beyond what any direct experiment has achieved. A relevant example is the possibility of implementing the Foldy-Wouthuysen transformation, which for decades has remained a purely abstract construction.