Quantum simulations of relativistic quantum physics in circuit QED

Feynman [1] is credited for introducing the idea of a quantum simulator. 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 10 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 reappeared later in 1996, when Lloyd [3] presented a general proof for Feynman's conjecture. 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 100-digit number with Shor's algorithm would need thousands of qubits to become non-trivial. 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 the most varied type of systems one can think of, including condensed-matter physics, high-energy physics, cosmology and quantum chemistry [4, 5]. When compared with 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

$$\begin{eqnarray} &&\mathrm{i}\,\hbar \frac{\mathrm{d} \psi}{\mathrm{d}t} = ( mc^2 \sigma_z + c \hat{P} \sigma_y ) \psi , \end{eqnarray} \tag{ 1 }$$

where σ_y and σ_z are Pauli matrices, $\hat {P}$ represents the momentum of a Dirac particle of mass m and c is the speed of light. Soon after its introduction, Schrödinger realized that there was something odd in the solution of this equation [8]. The position for a free Dirac particle can be integrated to give [9]

$$\begin{eqnarray} &&\hat{x}(t) = \hat{x}_0 + c^2 \hat{P} {\cal H}^{-1} t + \hat{Z}(t) \end{eqnarray} \tag{ 2 }$$

with ${\cal H}= mc^2 \sigma _z + c \hat {P} \sigma _y$ and

$$\begin{eqnarray} &&\hat{Z}(t) = \frac{\mathrm{i}}{2} \hbar c ( \sigma_y - c \hat{P} {\cal H}^{-1} ) {\cal H}^{-1} \left(\exp\left(- \mathrm{i} \,2 {\cal H} t/\hbar \right) - 1\right) . \end{eqnarray} \tag{ 3 }$$

This third term appears because the time derivative of the position operator is not a constant of motion, i.e. $[ \dot { \hat {x} } , {\cal H} ] \neq 0$, even for a free particle. Given a general initial wavepacket, the expectation value $\langle \hat {Z}(t) \rangle$ shows an oscillating pattern with amplitude ∼ℏ/mc and frequency ∼2mc²/ℏ that is known as Zitterbewegung or jittering motion [9]. This behavior 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⁻⁴ nm and a frequency of 10²¹ Hz, making its eventual detection very demanding.

There is another interesting property of the Dirac equation related to the behavior of relativistic Dirac particles under the effect of an external scalar potential,

$$\begin{eqnarray} &&\mathrm{i}\,\hbar \frac{\mathrm{d} \psi}{\mathrm{d}t} = \{ mc^2 \sigma_z + c \hat{P} \sigma_y + \Phi (\hat{X}) \} \psi . \end{eqnarray} \tag{ 4 }$$

In the case of a non-relativistic particle impinging against a semi-infinite step potential $\Phi (\hat {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 that equation (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² the particle with kinetic energy E will be mostly transmitted, for a larger one E − mc² < Φ < E + mc² it will be reflected. However, on making the potential Φ > E + mc² even larger 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 toward the later development of quantum field theories, in the last few decades it has driven 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, characterized by quasi-electrons behaving like relativistic Dirac particles in two dimensions [13]; a proposal for generating relativistic scattering in optical lattices [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–21], or even behave like 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. Firstly, since 2004 [24, 25], cQED has become the quantum platform with the most promising perspectives in terms of scalability and coherence. Secondly, 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–28] in scalable quantum network architectures [29] with delocalized and/or sequential interactions.

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 [24, 25, 31]. In addition, we introduce three external classical microwave drivings, two transverse 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

$$\begin{eqnarray} &&\fl {\cal H} = \frac{\hbar \omega_q}{2} \sigma_z +\hbar \omega \hat{a}^\dagger \hat{a} -\hbar g (\sigma^\dagger \hat{a} + \sigma \hat{a}^\dagger ) - \hbar \Omega ( \mathrm{e}^{ \mathrm{i} ( \omega t +\varphi) } \sigma + \mathrm{e}^{-\mathrm{i} ( \omega t +\varphi ) } \sigma^\dagger ) \nonumber \\ &&- \hbar \lambda ( \mathrm{e}^{\mathrm{i} \left( \nu t +\varphi \right)} \sigma + \mathrm{e}^{-\mathrm{i} \left( \nu t +\varphi \right) } \sigma^\dagger ) + \hbar \xi ( \mathrm{e}^{\mathrm{i} \,\omega t} \hat{a} + \mathrm{e}^{-\mathrm{i}\, \omega t} \hat{a}^\dagger ) \end{eqnarray} \tag{ 5 }$$

with $\sigma _y=\mathrm {i}\, \sigma - \mathrm {i}\,\sigma ^\dagger =i \left \vert g\rangle \langle e\right \vert - i \left \vert e\rangle \langle g\right \vert$ and $\sigma _z= \left \vert e\rangle \langle e\right \vert - \left \vert g\rangle \langle g\right \vert$, $\left \vert g\right \rangle$ and $\left \vert e\right \rangle$ 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 characterized 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. Firstly, the Hamiltonian in equation (5) can be simplified using the reference frame rotating with the resonator frequency ω:

$$\begin{eqnarray} &&\fl {\cal H}^{L_1} = - \hbar g (\sigma^\dagger \hat{a} +\sigma \hat{a}^\dagger ) - \hbar \Omega ( \mathrm{e}^{ \mathrm{i} \,\varphi } \sigma + \mathrm{e}^{-\mathrm{i}\, \varphi } \sigma^\dagger ) + \hbar \xi ( \hat{a} + \hat{a}^\dagger ) \nonumber \\ &&- \hbar \lambda ( \mathrm{e}^{\mathrm{i} \left[ (\nu-\omega) t + \varphi \right] } \sigma + \mathrm{e}^{-\mathrm{i} \left[ (\nu-\omega) t + \varphi \right] } \sigma^\dagger ). \end{eqnarray} \tag{ 6 }$$

Secondly, we will write this Hamiltonian in another interaction picture with respect to the term ${\cal H}_0^{L_1} = - \hbar \Omega ( \mathrm {e}^{ \mathrm {i}\, \varphi } \sigma + \mathrm {e}^{-\mathrm {i}\, \varphi } \sigma ^\dagger )$. The resulting expression reads

$$\begin{eqnarray} &&\fl {\cal H}^{\rm I} = - \frac{\hbar g }{2} ( \{ \left\vert+\rangle\langle+\right\vert - \left\vert-\rangle\langle-\right\vert + \mathrm{e}^{-\mathrm{i} 2 \Omega t} \left\vert+\rangle\langle-\right\vert - \mathrm{e}^{\mathrm{i} 2 \Omega t}\left\vert-\rangle\langle+\right\vert \} \mathrm{e}^{\mathrm{i} \varphi} \hat{a} + {\rm h.c.} ) \nonumber \\ &&- \frac{\hbar \lambda}{2} ( \{ \vert+\rangle\langle+\vert - \vert-\rangle\langle-\vert - \mathrm{e}^{-\mathrm{i}\, 2 \Omega t} \vert+\rangle\langle-\vert \nonumber\\ &&+ \mathrm{e}^{\mathrm{i} 2 \Omega t}\vert-\rangle\langle+\vert \} \mathrm{e}^{\mathrm{i} (\nu-\omega) t} + {\rm h.c.} )+ \hbar \xi ( \hat{a} + \hat{a}^\dagger ) , \end{eqnarray} \tag{ 7 }$$

where we have introduced the rotated spin basis $\left \vert \pm \right \rangle = \left (\left \vert g\right \rangle \pm \mathrm {e}^{-\mathrm {i}\, \varphi } \left \vert e\right \rangle \right )/\sqrt 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 the frequencies in (7). This implies that we can apply the rotating-wave approximation, yielding the effective Hamiltonian

$$\begin{eqnarray} &&{\cal H}_{\rm eff} = \frac{\hbar\lambda}{2} \sigma_z + \frac{\hbar g}{\sqrt2} \sigma_y \hat{p} + \hbar \xi \sqrt2 \, \hat{x} , \end{eqnarray} \tag{ 8 }$$

where φ = π/2 and we make use of standard quadratures of the electromagnetic field, i.e. dimensionless $\hat {x} =(\hat {a}+\hat {a}^\dagger )/\sqrt 2$, $\hat {p} =-\mathrm {i}(\hat {a}-\hat {a}^\dagger )/\sqrt 2$, at variance with $\hat {X}$ and $\hat {P}$ in equation (4), with the commutation relation $[ \hat {x} , \hat {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 equations (1) and (4), where the terms $\hbar g/\sqrt 2$ and ℏλ/2 are related to the speed of light and the mass, respectively. In addition, we have an external linear potential $\Phi = \hbar \xi \sqrt 2 \, \hat {x}$ 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 equation (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 the 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 $\langle \hat {x}_0 \rangle$ and average momentum $\langle \hat {p}_0 \rangle$,

$$\begin{eqnarray} &&\psi (x) = \pi^{-1/4} \exp\left\{ \mathrm{i} \langle \hat{p}_0 \rangle x \right\} \,\, \exp \left\{ - \frac{ ( x - \langle \hat{x}_0 \rangle )^2 }{2} \right\} . \end{eqnarray} \tag{ 9 }$$

This coincides with the x-quadrature representation of a coherent state of the electromagnetic field $\vert \frac { \langle \hat {x}_0 \rangle + \mathrm {i}\, \langle \hat {p}_0 \rangle }{\sqrt 2} \rangle = {\cal D} ( \frac { \langle \hat {x}_0 \rangle + \mathrm {i}\, \langle \hat {p}_0 \rangle }{\sqrt 2} ) \vert 0\rangle$, $\vert 0 \rangle$ being the vacuum state of the bosonic field and ${\cal D} ( \alpha ) = \exp \{ \alpha \hat {a}^\dagger - \alpha ^* \hat {a} \}$ 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 ns, which is well below standard decoherence times in cQED experiments.

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

$$\begin{eqnarray} &&{\cal H}_{\rm D}= \frac{\hbar\lambda}{2} \sigma_z + \frac{\hbar g}{\sqrt2} \sigma_y \hat{p} \end{eqnarray} \tag{ 10 }$$

has two important limits depending on the value of the mass of the particle. Firstly, in the massless case the Hamiltonian reduces to ${\cal H}_{\mathrm { D}} = \hbar g/\sqrt 2 \sigma _y \hat {p}$ allowing for straightforward interpretation of its dynamics in terms of phase-space representation. Figure 1(a) shows the evolution of a massless particle with initial state $\left \vert +,0\right \rangle$. With 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 $\left \vert 0\right \rangle$, to the final one in which the state of the field is $\left \vert g t /2\right \rangle$. This simulates the evolution of a massless free Dirac particle moving at the speed of light stemming from $\hbar g/\sqrt 2$. As we expect, the dynamics of this massless particle would be essentially unchanged if one assumes the initial state to have a non-zero value of the expectation value of momentum. This is depicted in figure 1(b), where the particle starts its trajectory with $\langle \hat {p}_0 \rangle = 2$. The particle moves exactly the same amount of space in the x-quadrature, ending in the coherent state $\vert g t /2 + \sqrt 2 \,{i} \rangle$.

Zoom In Zoom Out Reset image size

Secondly, when the mass of the particle becomes large enough, the dynamics should tend toward that of a non-relativistic one evolving under a free Schrödinger equation. For $\hbar \lambda /2 \gg \hbar g/\sqrt 2 \langle \hat {p} \rangle$, i.e. $mc^2 \gg c \langle \hat {P} \rangle$, one can derive a second-order effective Hamiltonian [17] for equation (8) to yield ${\cal H}_{\mathrm { NRel}} = \sigma _z \hbar (g^2 / \lambda ) \hat {p}^2$. This is indeed the non-relativistic 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 non-relativistic or relativistic regimes. In figure 1(e), we depict the evolution of the initial state $\vert e,0\rangle$ under the time-dependent Hamiltonian (6) for the case in which the mass of the particle is such that $\hbar \lambda /2= 4\times \hbar g/\sqrt 2$. This generates a squeezed vacuum state with the degree of squeezing increasing linearly in time. This behavior is in agreement with the aforementioned non-relativistic approximation for the Hamiltonian, proportional to $\hat {p}^2$. Such a simulation realizes a standard example in all quantum mechanics textbooks. Namely, the free evolution under the non-relativistic Schrödinger equation of an initial wavepacket with $\langle \hat {p}_0 \rangle =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 operator $\hat {p}^2$ is not simply squeezing. This can be clearly seen in figure 1(f), where we have taken a particle with initial state $\vert e,\sqrt 2\, {i} \rangle$, that is a wavepacket with non-zero kinetic energy such that $\langle \hat {p}_0 \rangle =2$. Now the wavepacket not only spreads over time, but its center moves linearly in time to the right as the expectation value of the position of a non-relativistic 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² becomes comparable to $c\langle \hat {P} \rangle$. This is the case of figures 1(c) and (d), where $\hbar \lambda /2=\hbar g/\sqrt 2$. Starting from initial states $\left \vert +,0\right \rangle$ and $\vert +,\sqrt 2\,{i} \rangle$, respectively, the Gaussian states of the field evolve into non-classical ones with a negative Wigner function through a highly nonlinear interaction. The effect is more evident for the case of an initial state $\left \vert +,0\right \rangle$ in figure 1(c). Initial states with non-zero kinetic energy, such as $\vert +,\sqrt 2\,{i} \rangle$ in figure 1(d), make the qubit–field coupling term in equation (8) dominant. This interesting finding shows that a Dirac evolution generates non-classical 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 figure 2, we plot the time evolution of $\langle \hat {x} \rangle$ for a particle prepared in an initial state $\left \vert +,0\right \rangle$ 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 figure 2 shows the evolution of a massless particle, λ = 0, whose expectation value of the position increases monotonically at the speed of light related to $\hbar g/\sqrt 2$. If we assume a mass such that $\hbar \lambda /2= \hbar g/ \sqrt 2$, then we observe that the evolution of $\langle \hat {x} \rangle$ followed by the dotted line departs considerably from the constant rectilinear behavior one would expect for a free particle. Likewise, for a mass even larger, choosing for instance $\hbar \lambda /2=4 \times \hbar g / \sqrt 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.

Zoom In Zoom Out Reset image size

To gain physical insight, we can derive an analytical expansion for the evolution operator of a free Dirac particle,

$$\begin{eqnarray} &&\fl \exp\left(-\mathrm{i}\, {\cal H}_{\rm D}t/\hbar\right) = {\rm exp}\left(-\mathrm{i}\,\mathrm{ }\,t \sigma_y \sqrt{ \frac{\lambda^2 }{4} + \frac{ g^2 }{2} \hat{p}^2 } \right) - \mathrm{i}\,\sigma_z \frac{\lambda t}{ 2} {\rm sinc}\left( \sqrt{ \frac{\lambda^2 t^2}{4} + \frac{g^2 t^2}{2} \hat{p}^2 } \right) \nonumber \\ &&- \mathrm{i}\,\sigma_y \left( \frac{ g t \hat{p} /\sqrt2 }{ \sqrt{ \frac{\lambda^2 t^2}{4} + \frac{ g^2 t^2}{2} \hat{p}^2 }} - 1 \right) {\rm sin}\left( \sqrt{ \frac{\lambda^2 t^2}{4} + \frac{g^2 t^2}{2} \hat{p}^2 } \right) \nonumber \\ &&= \exp\left(-\mathrm{i} \frac{gt}{\sqrt2}\sigma_y \hat{p}\right) - \mathrm{i} \frac{\lambda t}{2} \sigma_z {\rm sinc}\,(gt\hat{p}/\sqrt2) \nonumber \\ &&+ \frac{\lambda^2 t^2}{4} \left( - \frac{ {\rm sinc}\,(gt\hat{p}/\sqrt2)}{2} - \mathrm{i}\,\sigma_y \frac{ {\rm cos}\,(gt\hat{p}/\sqrt2) - {\rm sinc}\,(gt\hat{p}/\sqrt2)}{ gt\hat{p} \sqrt2 } \right) \nonumber \\ &&- \mathrm{i} \frac{\lambda^3 t^3}{8} \sigma_z \frac{ {\rm cos}\,(gt\hat{p}/\sqrt2) - {\rm sinc}\,(gt\hat{p}/\sqrt2)}{ g^2t^2\hat{p}^2 } + \cdots , \end{eqnarray} \tag{ 11 }$$

where sinc (x) = sin (x)/x. Here, a larger initial average momentum $g t\langle \hat {p}(0) \rangle /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 $\hat {Z}(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, $\hat {Z}(t)$ does not necessarily vanish. In addition, the preparation of a purely positive spinor for a massive Dirac particle is not straightforward. In [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],

$$\begin{eqnarray} &&\hat{X}_{\rm FW}(t) = \hat{X}_0 + c^2 \hat{P} {\cal H}_{\rm D}^{-1} t + \hat{\Delta} , \end{eqnarray} \tag{ 12 }$$

where $\hat {\Delta }$ is the time-independent operator

$$\begin{eqnarray} &&\hat{\Delta} = \frac{ \hbar c \sigma_x }{2 \hat{E}} - \frac{ \hbar c^3 \sigma_x \hat{P}^2 }{2 \hat{E} ^2 ( \hat{E} + m c^2)} , \end{eqnarray} \tag{ 13 }$$

σ_x = i σ_zσ_y, $\hat {E}=\sqrt {c^2 \hat {P}^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 of $\hat {X}_{\mathrm { FW}}(t)$ from an experiential point of view is that it 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 [33] are also credited for having introduced a unitary transformation 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

$$\begin{eqnarray} &&{\cal H}_{\rm FW} = \hat{S}_{\rm FW} \, {\cal H}_{\rm D} \, \hat{S}_{\rm FW}^\dagger = \sigma_z \sqrt{ \frac{\hbar^2\lambda^2}{4} + \frac{\hbar^2\,g^2}{2} \hat{p}^2 }, \end{eqnarray} \tag{ 14 }$$

where

$$\begin{eqnarray} &&\hat{S}_{\rm FW} = \exp \left\{ - \mathrm{i} \quad {\rm atan} \left( \frac{g}{\lambda\sqrt2} \sigma_x \hat{p} \right) \right\} . \end{eqnarray} \tag{ 15 }$$

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 figure 2(b). This shows the evolution of $\langle \hat {x} \rangle$ for a massive particle with $\hbar \lambda /2=4\times \hbar g/\sqrt 2$, just like the dash-dotted one, but using the evolution operator

$$\begin{eqnarray} &&{\cal U}_{\rm FW}(t)=\exp\left(-\mathrm{i}\, {\cal H}_{\rm FW}t/\hbar\right)= \hat{S}_{\rm FW} \,\, \exp\left(-\mathrm{i} \,{\cal H}_{\rm D}t/\hbar\right) \,\, \hat{S}_{\rm FW}^\dagger \end{eqnarray} \tag{ 16 }$$

instead of ${\cal U}_{\mathrm { D}}(t)=\exp \left (-\mathrm {i}\, {\cal H}_{\mathrm { D}}t/\hbar \right )$. Now the particle remains centered at $\langle \hat {x} \rangle =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 \langle \hat {p} \rangle /(\lambda \sqrt 2)$ is small enough, we can linearize the interaction in the unitary operator $\hat {S}_{\mathrm { FW}}$. This means that

$$\begin{eqnarray} &&{\cal U}_{\rm FW}(t) \approx \exp\left(-\mathrm{i} \frac{g }{\lambda \sqrt2} \sigma_x \hat{p }\right) \,\, \exp\left(-\mathrm{i}\, {\cal H}_{\rm D}t/\hbar\right) \,\, \exp\left(\mathrm{i}\, \frac{g}{\lambda \sqrt2} \sigma_x \hat{p }\right) , \end{eqnarray} \tag{ 17 }$$

where the first and last steps involve the realization of a Dirac equation for a massless particle with an evolution time t = λ⁻¹. 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, g being in the strong-coupling regime. This is the case, for instance, of a gradiometer type of coupling seen in flux qubits [30]. The solid line in figure 2(a) depicts the evolution of $\langle \hat {x} \rangle$, computed under these assumptions, for the same particle of large mass, with $\hbar \lambda /2=4\times \hbar g/\sqrt 2$. Even for such a moderate value of the mass, it is evident how oscillations seen for the Dirac evolution become significantly reduced.

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 ${\cal H}_{\mathrm { K}} = \hbar g/\sqrt 2 \sigma _y \hat {p} + \hbar \xi \sqrt 2 \, \hat {x}$. In figures 3(a) and (b) the evolution of the initial state $\vert +,0 \rangle$ and $\vert +,\sqrt 2\,{i} \rangle$, respectively, is shown. We use the same set of parameters as in figure 1, with only the addition of a linear potential of strength ξ = g/2. Again, starting with the qubit state, $\left \vert +\right \rangle$, allows for a simple interpretation. The initial state of the field remains coherent while being subjected to two independent displacements: namely, as for the free-particle case, there is one in the x-quadrature proportional to $g/\sqrt 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. The rate at which this happens is given by $\xi \sqrt 2$. This explains why the only noticeable difference between the plots in figures 1 and 3 is in the value of $\langle \hat {p} (t) \rangle$.

Zoom In Zoom Out Reset image size

If the mass of the particle is large enough, then the non-relativistic Schrödinger equation should be a good approximation. Its Hamiltonian can be written as ${\cal H}_{\mathrm { NRel}} = \hbar \sigma _z \hat {p}^2/\lambda + \hbar \xi \sqrt 2 \, \hat {x}$, which guarantees that any initial Gaussian state should remain Gaussian during the evolution. This is indeed what we see in figures 3(e) and (f) where we have prepared initial states $\left \vert e,0\right \rangle$ and $\vert e,\sqrt 2\,{i}\rangle$, respectively, with a mass such that $\hbar \lambda /2=4\times \hbar g/\sqrt 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 ns 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 figure 3(c), which corresponds to the initial state $\left \vert +,0\right \rangle$ of a massive particle with $\hbar \lambda /2= \hbar g/\sqrt 2$. The wavepacket 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 figure 3(d) where we prepare $\vert +,\sqrt 2\,{i} \rangle$, then the particle tunnels slightly to stop and break up eventually.

As we did in equation (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

$$\begin{eqnarray} &&\fl \exp\left\{-\mathrm{i} \left( \frac{\lambda t}{2} \sigma_z + \frac{ g t }{2} \sigma_y \hat{p} + \xi \sqrt2 \,t\, \hat{x} \right) \right\} \nonumber \\ &&= \exp\left(-\mathrm{i} \frac{{\cal H}_{\rm D}}{\hbar} t \right) \exp\left(-\mathrm{i} \,\xi \sqrt2 \,t\, \hat{x} \right) \exp\left(\mathrm{i}\, \frac{ g \xi t^2}{2} \sigma_y\right) + O\left(t^3 \right) , \end{eqnarray} \tag{ 18 }$$

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 one to characterize the type of dynamics 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–28], already available in cQED.

We have shown that cQED provides a powerful platform for simulating 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 the presence of a linear potential. The degree of controllability offered by this scheme would allow one 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.

We thank Leo DiCarlo, Jorge Casanova and Lucas Lamata for insightful discussions. The authors acknowledge funding from EPSRC EP/H050434/1, Spanish MINECO FIS2012-36673-C03-02; UPV/EHU UFI 11/55; UPV/EHU PhD fellowship; Basque Government IT472-10; SOLID, CCQED, PROMISCE and SCALEQIT European projects.