Optical simulation of the free Dirac equation

TL;DR

This paper demonstrates a classical optics setup that simulates the one-dimensional Dirac equation, allowing observation of relativistic quantum phenomena like zitterbewegung using tunable optical components.

Contribution

It introduces a novel optical simulation method for the Dirac equation that does not require specially engineered media, enabling study of relativistic effects in a controllable classical system.

Findings

Observed zitterbewegung oscillations in the optical setup

Measured mean-position operator without oscillations

Demonstrated tunable simulation of relativistic quantum phenomena

Abstract

We present a classical optics simulation of the one-dimensional Dirac equation for a free particle. Positive and negative energy components are represented by orthogonal polarizations of a free propagating beam, while the spatial profile represents the spatial wave function of the particle. Using a completely tunable time parameter, we observed the oscillation of the average value of the Dirac position operator - known as zitterbewegung (ZB). We are also able to measure the so called mean-position operator which is a single-particle observable and presents no oscillations. Our work opens the way for optical simulation of interesting phenomenon of relativistic systems, as well as condensed-matter physics, without any requirement for specially engineered medium.

Tables

Table 1. Table 1: Summary of the optical analogy

Optical System	Simulated System
Vertical transverse position	$x$
Transverse profile of horizontal polarization	${\psi }_{D_1}\left(x\right)$
Transverse profile of vertical polarization	${\psi }_{D_2}\left(x\right)$
QWP@45°	$\hat{Q}=e^{i\frac{\pi }{4}}/\sqrt{2}\left(𝟙-i{\sigma }_1\right)$
HWP@45°	$\hat{H}={\sigma }_1$
SLM printing phase $f(x,y)$	$\hat{P}\left[f(x,y)\right]=$ $\text{exp}\left[i\text{f}(x,y)\right]{\sigma }_{+}{\sigma }_{-}+{\sigma }_{-}{\sigma }_{+}$
Normalized horizontal polarization intensity at $x$	${\left|{\psi }_{D_1}\left(x\right)\right|}^2$

Equations

$$i\hslash \frac{\partial }{\partial t}{\psi }_D={\hat{\mathcal{H}}}_D{\psi }_D=(c\hat{p}{\sigma }_1+m{c}^2{\sigma }_3){\psi }_D,$$

$${\hat{x}}_D(t)=$$

$$-\frac{c\hslash {\hat{\mathcal{H}}}_D^{-1}}{2i}\left(e^{2i{\hat{\mathcal{H}}}_{D}t/\hslash }-1\right)\left(c\hat{p}{\hat{\mathcal{H}}}_D^{-1}-{\sigma }_1\right).$$

$${\hat{x}}_{FW}(t)={\hat{x}}_{FW}(0)+c^2\hat{p}{\hat{\mathcal{H}}}_D^{\prime -1}t$$

$$\hat{U}(p)=\hat{Q}\hat{P}\left[-\theta (p)\right]\hat{H}\hat{P}\left[\theta (p)\right]\hat{Q},$$

$${\psi }_D(x,t=0)=\frac{e^{-i\pi x^2/(\lambda R)}{e}^{-x^2/(4{\Delta }^2)}}{(\sqrt{2\pi }2\Delta {)}^{1/2}}\left(\begin{array}{c}a\\ b\end{array}\right),$$

$$H=c({\sigma }_1{p}_x+{\sigma }_2{\mathrm{p}}_y)+m{c}^2{\sigma }_3.$$

$${\theta }^{\prime }(\mathbf{p})=\{\begin{array}{ll}{\mathrm{tg}}^{-1}(p_y/p_x),& p_y>0\text{ or }{p}_y=0,p_x>0\\ \pi +{\mathrm{tg}}^{-1}(p_y/p_x),& p_y<0\text{ or }{p}_y=0,p_x<0,\end{array}$$

$$e^{iS_{2}}=\left(\begin{array}[]{cc}\cos\theta(\mathbf{p})&e^{-i\theta^{\prime}}\sin\theta(\mathbf{p})\\ -e^{i\theta^{\prime}}\sin\theta(\mathbf{p})&\cos\theta(\mathbf{p})\end{array}\right).$$

$$i\hslash \frac{\partial \psi (\mathbf{x},t)}{\partial t}=\left[c\mathbf{\alpha }\cdot \mathbf{p}+m{c}^2\beta \right]\psi (\mathbf{x},t)$$

$$\psi(\mathbf{x},t=0)=\begin{cases}\frac{1}{\sqrt{L}}\left(\begin{array}[]{c}\tilde{\phi}_{1}(x,y)\\ \vdots\\ \tilde{\phi}_{4}(x,y)\end{array}\right)&-\frac{L}{2}<z<\frac{L}{2}\\ 0&\mathrm{elsewhere}\end{cases}.$$

$$\tilde{\psi}(t=0,\boldsymbol{p})=\left(\begin{array}[]{c}\tilde{\phi}_{1}(p_{x},p_{y})\\ \vdots\\ \tilde{\phi}_{4}(p_{x},p_{y})\end{array}\right)\sqrt{\frac{L}{2\pi\hbar}}\left(\frac{\sin{\frac{p_{z}L}{2\hbar}}}{\frac{p_{z}L}{2\hbar}}\right),$$

$${\tilde{\psi }}^{\prime }(\mathbf{p},t=0)\approx \left(\cos \left(\theta (p_x,p_y,0)\right)+\beta \frac{{\alpha }_x{p}_x+{\alpha }_y{p}_y}{|(p_x,p_y,0)|}\sin \theta (p_x,p_y,0)\right)\tilde{\psi }(\mathbf{p},t=0),$$

$$\psi^{\prime}(\boldsymbol{p},t=0)=\left(\begin{array}[]{cccc}\cos{\theta}&0&0&\frac{p_{x}-ip_{y}}{|\mathbf{p}|}\sin{\theta}\\ 0&\cos{\theta}&\frac{p_{x}+ip_{y}}{|\mathbf{p}|}\sin{\theta}&0\\ 0&\frac{-p_{x}+ip_{y}}{|\mathbf{p}|}\sin{\theta}&\cos{\theta}&0\\ \frac{-p_{x}-ip_{y}}{|\mathbf{p}|}\sin{\theta}&0&0&\cos{\theta}\end{array}\right)\left(\begin{array}[]{c}\phi_{1}\\ \phi_{2}\\ \phi_{3}\\ \phi_{4}\end{array}\right)\sqrt{\frac{L}{2\pi\hbar}}\left(\frac{\sin{\frac{p_{z}L}{2\hbar}}}{\frac{p_{z}L}{2\hbar}}\right)$$

$$i\hslash \frac{\partial }{\partial t}{\psi }_D=(c\hat{p}{\sigma }_1+m{c}^2{\sigma }_3+V(x)){\psi }_D,$$

$$V(x)=V_1(x){\sigma }_1.$$

$${\psi }_D=e^{-i\frac{1}{\hslash c}\int V_1(x^{\prime })1d{x}^{\prime }}{\varphi }_D.$$

Full text

I Introduction

Although the Dirac equation Dirac1928 represents a historical landmark in the quantum description of relativistic systems – satisfactorily explaining the origin of spin and predicting the existence of antimatter Anderson1933 – it still provokes a lot of discussion about its interpretation, even when applied to describe the simplest physical system, that is, a free particle. In this case, the Dirac equation predicts intriguing phenomena, for instance the ZB Schrodinger1930 and Klein’s paradox Klein1929 , which impede the single-particle (SP) interpretation of the Dirac equation. As fundamental effects in the understanding of relativistic influence over quantum theory, they have contributed to the transition to the many-body quantum field approach Krekora04 .

The important technical difficulties involved in the direct observation of several relativistic quantum predictions have led to an increased interest in their simulations in trapped ions Gerritsma10 ; Gerritsma11 , photonic crystals Zhang08 , confined light Otterbach09 , graphene Katsnelson06 , optical superlattices Dreisow10 , Bose-Einstein condensates Salger11 ; LeBlanc13 and ultracold atoms Vaishnav08 . Among the unexpected effects of Dirac equation, the ZB – the flickering motion of a free relativistic quantum particle described by a Dirac wavefunction with positive and negative energy components – is one of most investigated due to its interesting counterintuitive nature.

Despite the previous studies Gerritsma10 ; Dreisow10 ; LeBlanc13 , an important feature was not explored: the Dirac position operator related to the flickering motion is not a SP observable, i.e. it cannot be written as the direct sum of its positive and negative energy projections. From a SP perspective, Dirac’s equation must be interpreted as the simultaneous solving of two independent problems, for instance, the single free evolution of both electron and positron. Therefore, there is no meaning in assigning a physical interpretation to an operator that is not SP. This fact can be well understood by making explicit the connection between the theory of irreducible unitary representations of the Poincaré group and the explicitly covariant reducible representation associated to Dirac’s equation Caban03 . In the SP approach, physical results must be obtained by projecting SP observables over the subspace corresponding to the problem of interest.

For Dirac’s theory, a SP position observable exists and is obtained using the so-called Foldy–Wouthuysen transformation (FWT) Foldy50 – a momentum dependent unitary transformation that diagonalizes Dirac’s Hamiltonian and is at the kernel of important algorithms used to obtain quantum relativistic corrections Obukhov01 ; Quach15 . This observable is often called mean-position operator and does not exhibit the oscillatory behavior characteristic of the ZB, a result that generates doubts concerning its actual existence.

Here, the simulation of the Dirac one-dimensional free evolution and the ZB is performed using the transverse degrees of freedom of a paraxial light beam, where different components of the spinor are represented by different polarization components of the beam. This physical setup is well suited for the purpose of quantum simulation, as it allows for implementation of the dynamical phases with easy tuning of the important physical parameters Lemos12 . Besides being a proof-of-concept for the optical simulation of Dirac particles, the present approach differs from others in the theoretical procedure adopted to perform the simulation. Previous works perform a direct simulation of the 1+1D Dirac Hamiltonian, while the present approach performs the evolution in the diagonalized Foldy–Wouthuysen representation (FWR) and permits one to switch back and forth between this and Dirac’s representation, allowing us to investigate the behavior of both the Dirac position and the mean-position operators.

II Dirac Equation and position operator

Consider the 1D Dirac equation

[TABLE]

where $c$ is the speed of light, $\hat{p}$ is the momentum operator, $m$ is the mass of the particle, and $\sigma_{i}$ are the usual Pauli matrices. The information of this system is encoded in the spinor $\psi_{D}$ that has only two components which are related to positive and negative energy states in the particle’s rest frame, i.e. spin degrees of freedom are eliminated by the dimensional constrain Schwabl . In Dirac’s coordinate representation, the momentum operator $\hat{p}$ assumes the usual form $-i\hbar\frac{\partial}{\partial x}$, where $x$ is the so-called Dirac coordinate associated to the multiplication operator $\hat{x}_{D}\psi_{D}(x)=x\psi_{D}(x)$. Since the Hamiltonian operator is not diagonal in this representation, the positive and negative energy eigenstates are non-trivial and assume, respectively, the forms $\psi_{p}^{+}(x,t)=u(p)e^{-i\varepsilon(p)t/\hbar}e^{ipx/\hbar}$ and $\psi_{p}^{-}(x,t)=v(p)e^{i\varepsilon(p)t/\hbar}e^{-ipx/\hbar}$, with $\varepsilon(p)\equiv\sqrt{(mc^{2})^{2}+(pc)^{2}}$, $u(p)=[2mc^{2}(\varepsilon(p)+mc^{2})]^{-1/2}\begin{pmatrix}\varepsilon(p)+mc^{2}&cp\end{pmatrix}^{T}$ and $v(p)=[2mc^{2}(\varepsilon(p)+mc^{2})]^{-1/2}\begin{pmatrix}cp&\varepsilon(p)+mc^{2}\end{pmatrix}^{T}$.

The non-diagonal form of $\hat{\mathcal{H}}_{D}$ in Dirac’s representation is evinced by the commutator $[\hat{\mathcal{H}}_{D},\hat{x}_{D}]=-ic\hbar\sigma_{1}$ and leads to the Heisenberg picture evolution given by Thaller

[TABLE]

The first two terms on the right represent the expected linear time evolution of a free particle, the last term being associated to the ZB. This flickering motion is accompanied by other particularities of the $\hat{x}_{D}$ operator. Indeed, the evolution given in (2) is derived from the equation of motion $\dot{\hat{x}}_{D}=\frac{i}{\hbar}[\hat{\mathcal{H}}_{D},\hat{x}_{D}]=c\sigma_{1},$ which implies that, although $\langle\dot{\hat{x}}_{D}\rangle=\langle c^{2}\hat{p}\hat{\mathcal{H}}_{D}^{-1}\rangle$, the eigenvalues associated to the velocity ${\dot{\hat{x}}_{D}}$ are restricted to $\pm c$, a remarkable result which contributes to raise doubts as to the correct interpretation of $\hat{x}_{D}$ as definition of position. These peculiarities of the operator $\hat{x}_{D}$ arise from the fact that this is not a SP observable, i.e. $\hat{x}_{D}\neq\hat{P}_{+}\hat{x}_{D}\hat{P}_{+}^{\dagger}+\hat{P}_{-}\hat{x}_{D}\hat{P}_{-}^{\dagger}$, where $\hat{P}_{\epsilon}=\frac{1}{2mc^{2}}\left(\begin{smallmatrix}mc^{2}+\epsilon\varepsilon(p)&-\epsilon cp\\ \epsilon cp&mc^{2}-\epsilon\varepsilon(p)\end{smallmatrix}\right)$ is the projection operator over the subspace of states with energy sign $\epsilon$.

To obtain a SP position, the FWT must be applied to diagonalize the Dirac Hamiltonian. For the 1D Dirac free particle, this canonical transformation is given by the momentum dependent unitary operator $\hat{U}(\hat{p})=e^{i\hat{S}(\hat{p})}$ with $\hat{S}(\hat{p})\equiv\frac{\sigma_{2}}{2}\mathrm{tg}^{-1}\left(\frac{\hat{p}}{mc}\right)$. In the resulting FWR, the original Dirac Hamiltonian is given by $\hat{\mathcal{H}}_{D}^{\prime}=\sigma_{3}\varepsilon(\hat{p})$, and the former $\hat{x}_{D}$ operator by $\hat{x}_{D}^{\prime}=\hat{x}_{FW}+\frac{\hbar mc^{3}}{2\varepsilon(p)^{2}}\sigma_{2}$, where $\hat{x}_{FW}$ is the new multiplication operator in the FWR. The operator $\hat{x}_{FW}$ is the so called mean-position operator and, unlike the operator $\hat{x}_{D}$, it is a SP observable since $\hat{x}_{FW}=\hat{P}_{+}^{\prime}\hat{x}_{FW}\hat{P}_{+}^{\prime}+\hat{P}_{+}^{\prime}\hat{x}_{FW}\hat{P}_{+}^{\prime}$, where $\hat{P}_{\epsilon}^{\prime}=\left(\begin{smallmatrix}\delta_{\epsilon+}&0\\ 0&\delta_{\epsilon-}\end{smallmatrix}\right)$ are the energy projectors in the new representation.

Aside from being SP, the operator $\hat{x}_{FW}$ also satisfies the equation $\dot{\hat{x}}_{FW}=c^{2}\hat{p}\hat{\mathcal{H}}_{D}^{\prime}$, resulting in the Heisenberg picture evolution Thaller

[TABLE]

that is linear in time, as expected for a free particle. Thus, as stated before, the ZB does not occur for this operator.

Here, as a proof-of-concept for the simulation of relativistic systems using free propagating light beams, the simulation of both the Dirac and FWRs will be performed in a single setup. This difference with other simulation procedures open the possibility for future investigations on more complex FWTs associated to relativistic scenarios involving interactions.

III Simulation Protocol and Experiment

One way to simulate the dynamics associated to Eq. (1) is to directly implement the evolution operator $\exp\left(-\frac{i\hat{\mathcal{H}}_{D}t}{\hbar}\right)$, which is usually a tough task due to the non-diagonal character of $\hat{\mathcal{H}}_{D}$. This difficulty can be overcome by using the FWT, since this transformation allows the time evolution operator to be written as $\exp\left(-\frac{i\hat{\mathcal{H}}_{D}t}{\hbar}\right)=\hat{U}^{-1}(\hat{p})\exp\left(-\frac{i\hat{\mathcal{H}}_{D}^{\prime}t}{\hbar}\right)\hat{U}(\hat{p}).$ This operator can be implemented in an optical beam by considering the vertical coordinate on the transverse plane as the particle’s position and the horizontal (vertical) polarization as the superior (inferior) component of the spinor. The horizontal spatial degrees of freedom on the transverse plane play no relevant role in the experiment. Although a spinor is a mathematical object which transforms very specifically under a reference frame change, it is not a concern for this simulation since the reference frame is assumed to be fixed.

The optical transformations required for the simulation are polarization transformations (acting as nondiagonal operators) and phase shifts (used to introduce momentum dependent phases). The former are obtained with the suitable application of wave plates and the last are realized by spatial light modulators (SLMs), which are able to imprint programmable position dependent phases $\textrm{exp}[i\>\textrm{f}(x,y)]$ in the horizontal polarization. To transform these position-dependent phases into momentum-dependent ones, optical Fourier transforms are used. The action of a quarter wave plate (QWP) set to $45^{\circ}$ is given by the operator $\hat{Q}=\nicefrac{{e^{i\tfrac{\pi}{4}}}}{{\sqrt{2}}}\left(\mathbb{1}-i\sigma_{1}\right)$, while $\hat{H}=\sigma_{1}$ describes a half wave plate at $45^{\circ}$. The action of a SLM is equivalent to applying $\hat{P}[f(x,y)]=\textrm{exp}[i\>\textrm{f}(x,y)]\sigma_{+}\sigma_{-}+\sigma_{-}\sigma_{+}$ over the transverse profile spinor. Using this operator representation for the optical devices, it follows that

[TABLE]

with $2\theta(p)=\mathrm{tan}^{-1}\left(\frac{p}{mc}\right)$. We express the inverse FWT in an analogous fashion. As the Hamiltonian is diagonal in the FWR, the transformed time evolution operator is achieved via the application of the dynamical phase $\exp\left[\pm it\varepsilon(p)/\hbar\right]$ in each spinor component using waveplates and the SLM, which concludes the simulation. A summary of the analogy between the optical simulator and the simulated system is given in Table I.

The experimental scheme is shown in Fig.1. A He-Ne laser with wavelength 632.8 nm and two Holoeye reflective SLMs, each of which divided into halves to operate twice, are used. The Fourier transforms are made by plano-convex cylindrical lenses with 150 mm focal distance such that the position space (mirrors and camera) is in one focal plane and the momentum space is in the opposite focal plane where the SLM is located. The momentum $p$ and the position on the SLM, $X$, are connected by $X=\frac{\lambda f}{h}p$, where $f$ is the focal distance and $\lambda$ is the laser wavelength Saleh . In terms of $X$, the applied phases become $2\theta(X)=\mathrm{tan}^{-1}\left(\frac{h}{mc}\frac{X}{\lambda f}\right)$ and $t\varepsilon(p)/\hbar=2\pi t\sqrt{\left(\frac{X}{\lambda f}\right)^{2}+\left(\frac{mc}{h}\right)^{2}}$, so the parameters we need to set are the speed of light $c$ and the Compton wavelength $\lambda_{C}=h/mc$, which are easily tunable since they enter as programmable parameters in the imprinted phases. Since the time coordinate also comes up as a programmable parameter, we could in principle take measurements for as many time values as we wish inside a time interval. This also implies that the unit of measurement for time is an arbitrary $\tau$. In this realization we chose $\Delta t/\tau=1$.

The laser produces a Gaussian spatial profile separable in the $x$ and $y$ coordinates, so the initial spinor is

[TABLE]

where $a$ and $b$ are the proportions of horizontal and vertical polarizations ($|a|^{2}+|b|^{2}=1$), $\Delta$ is the beam width in the vertical direction, and $R$ is the vertical radius of curvature of the beam in the initial position plane. The propagation and Gouy terms of the Gaussian beam only introduce global phases which do not affect the dynamical evolution Saleh . We start with $a=-b=1$, but changing $a$ and $b$ would enable us to prepare different positive and negative energy superpositions. Two cylindrical lenses are placed before the first position space in order to manipulate the initial momentum distribution which depends on $R$ and therefore on $\Delta$. Using a beam profiler, we determined the initial state parameters to be $\Delta=48,6\,\mu m$ and $\lambda R/\pi=(2.2\Delta)^{2}$.

The average position of the simulated particle is calculated as $\langle\hat{x}_{D}\rangle(t)=\sum_{i=1,2}\int dx\,x|\psi_{Di}(x,t)|^{2},$ where $\left|\psi_{Di}(x,t)\right|^{2}$ is proportional to the light intensity of polarization component $i$ at position $x$ on the transverse plane measured by a CCD camera placed at the output position space. Each instant of time corresponds to one programmable-phases configuration and one intensity-profile measurement. It is worth noting that the evolved state is accessible for any time value.

IV Results

A summary of our experimental procedure for the particular case $\lambda_{C}=5\Delta$ is depicted in Fig.2. In Fig.2-(a) we present the mean position $\langle\hat{x}_{D}\rangle(t)$ as a function of $t$, the ZB being evidenced by the oscillatory behavior. The solid red line is the theoretical prediction, while points are experimental results obtained from the images shown in Fig.2-(b). The shaded region represents uncertainty of one $\sigma$. Fig.2-(b) shows samples of the data collected by the CCD camera for some instants of time, the $x$ distribution being obtained by considering only a fixed $y$ coordinate at the center of the beam. The $x$ distributions used to calculate $\langle\hat{x}_{D}\rangle(t)$ are shown in Fig.2-(c). This procedure assumes that $x$ and $y$ intensity distributions remain separable troughout all the apparatus. This is true in the ideal case, however the cylindrical lenses can introduce some non-separability as one can see in the slightly tilted eliptical intensity pattern shown in Fig. 2-(b). The non-separability causes the initial state to be not entirely pure. Since our experimental results agree well with theory, we conclude that these effects are negligible for the present experiment.

For a fixed initial state and speed of light $c=0.1\,\nicefrac{{\Delta}}{{\tau}}$, we measured the average position in Dirac’s representation for different values of the Compton wavelength, as shown in Fig.3. We fitted the average position with the function $vt+A\sin(\omega t+\delta)$ for each $\lambda_{C}$ to estimate the mean velocity, amplitude and frequency of the oscillation. These quantities are shown in Fig.4. Our experimental results are in agreement with the expected linear dependence of amplitude and inverse dependence of frequency on $\lambda_{C}$ for small $\lambda_{C}$ nota , as can be seen in Fig.4-(b). This is consistent with the fact that the ZB visibility in Fig.3 increases for smaller values of $\lambda_{C}$.

The different inclinations exhibited in Fig.3 are due to the fact that each mass, i.e. Compton wavelength, is associated to a different velocity distribution, even the momentum distribution being the same for all values of $\lambda_{C}$. Although the initial state (5) has zero average momentum, this is not true for the mean velocity in Fig.4-(a). As is expected, the mean velocity falls quadratically with $\lambda_{C}$ for large masses (small $\lambda_{C}$), while it is close to the speed of light for very small masses ($\lambda_{C}=100\Delta$).

The agreement between the experimental ZB data and the theoretical predictions confirm that our optical setup is well suited for the study of 1+1D relativistic dynamical systems, the theoretical extension to larger dimensions being discussed in the Appendix. Beside serving as a proof-of-concept, the proposed setup permits to investigate the system in the FWR, an interesting possibility since it allows to describe the dynamics of the system according to the single-particle perspective, i.e. assigning physical sense only to the projections of single-particle operators over the subspaces of definite sign of energy.

From the SP perspective, operators that are not block-diagonal in the FWR, as is the case for $\hat{x}_{D}$, have no physical meaning, since they mix components of positive and negative energy that are associated to two distinct problems. On the other hand, operators that are block-diagonal in the FWR, as $\hat{x}_{FW}$ or the Dirac Hamiltonian, may have a physical sense assigned to their positive and negative projections. In this sense, the correct description of the dynamics of a single electron (positron), for example, should be given by $\hat{x}_{FW,+}=P^{\prime}_{+}\hat{x}_{FW}P^{\prime}_{+}$ ($\hat{x}_{FW,-}=P^{\prime}_{-}\hat{x}_{FW}P^{\prime}_{-}$). In our setup, the average $\langle\hat{x}_{FW}\rangle$ measured using both positive and negative components of the spinor can be obtained by measuring the transverse profile of the beam before the inverse FWT. However, since the positive and negative components of the spinor are encoded in the horizontal and vertical polarizations of the beam in the FWR, the single-particle position dynamics described by $\hat{x}_{FW,+}$ (particle) and $\hat{x}_{FW,-}$ (anti-particle) is also accessible by simply selecting one of the polarizations prior to the CCD measurement in the FWR.

Experimental results for the mean-position operator are shown in Fig. 5 for $\lambda_{C}=5\Delta$. The experimental data concerning the ZB effect for the $\hat{x}_{D}$ operator is also plotted for comparison (blue points). Measurements of $\langle\hat{x}_{FW,+}(t)\rangle$ and $\langle\hat{x}_{FW,-}(t)\rangle$ are shown as black dots. As is expected from the independence of the two problems in the SP description, we have two independent mean trajectories corresponding to the free evolution of the particle and the corresponding anti-particle. The ZB is not present for these mean trajectories and a linear behavior in time is observed, as it was expected. As mentioned earlier, we were also able to measure the mean value of the FW mean-position operator $\hat{x}_{FW}$. The results are plotted as the red circles and represent an average of the positive and negative projections cases. The small deviation from a perfect linear behavior can be explained assuming that the SLMs do not modulate all the incident light but a fraction of it, as is shown in the inset picture which shows the same mean values as the experimental plot but obtained from a numerical simulation of the experiment for modulation efficiency of $95\%$ in each SLM.

V Discussion and Conclusions

Our experiment demonstrates how relativistic dynamics can be studied using classical optics, and opens the way to more sophisticated investigations. For this purpose it would also be desirable to produce more general initial states. This can be accomplished using intensity and phase masks in the initial state preparation. The state produced in this experiment had zero avarage momentum, but simply shifting the momentum in all SLMs phases by the same $\Delta p$ can be interpreted as if the state has non-vanishing average momentum.

In principle the method implemented in this simulation using the FWT could be applied for other simulation schemes of the Dirac equation, however this transformation requires applying a phase shift that is proportional to the inverse tangent of momentum. In our approach the application of this phase is fairly easy, thanks to the spatial light modulator (SLM). However, in other systems, this is quite challenging. Typically in continuous variable quantum simulators one can implement Gaussian Hamiltonians, but non-Gaussian operations (third order and above) are quite difficult Tasca11 . Thus, we believe that our approach is quite interesting in this regard, as it allows one to employ the FWT and investigate relevant aspects of it.

Albeit here we present the theory and results for the 1+1D case, the extension for 2+1D and for some initial states in 3+1D is straightforward as shown in the Appendix. The first is a direct extension considering the second transverse coordinate of the beam as the second spatial degree of freedom of the simulated particle. A 2+1D simulation also allows for investigation of electronic behavior in bidimensional condensed matter systems such as graphene Katsnelson06 , but still do not present any spin effect. For a general 3+1D simulation it would be necessary a third beam coordinate what is not available in this scheme. In spite of this limitation, in the Appendix we show a class of initial states which dependence on the third coordinate does not alter the time evolution. In this case, the two extra spinor components are provided by the polarization components of a second beam.

It is well known that there is no exact FWT for the non-free Dirac equation, i.e., if we add a potential to the free Dirac equation (1) the Dirac Hamiltonian becomes no longer diagonalizable with one single unitary transformation Schwabl . This seems to be a very limiting factor of our simulation technique and indeed it is if we try to implement the actual FWT for a potential problem. Instead of doing so, we can try to find other kinds of unitary transformations which reproduce the time evolution and are experimentaly feasible with the available optical elements. Up to now we know that at least for a particular class of potentials it is possible to break the time evolution operator into a position dependent phase which carries all the information about the potential followed by the free evolution presented in this work. This particular case is discussed on Appendix B.

In conclusion, we have presented an all-optical simulation of the dynamics of a one-dimensional relativistic free point particle, where the beam’s spatial profile plays the role of the particle’s wavefunction, and its orthogonal polarization components are associated to spinor components. Our experiment is based on the diagonalization of the Dirac Hamiltonian using the FWT, which allowed for the decomposition of the unitary evolution into operations that are realizable with off-the-shelf optical components. Adjusting the tunable time parameter we observed the oscillatory ZB phenomenon for Dirac’s position operator. Using our experimental FWT, we were also able to address this phenomenon from a single-particle perspective, where the position description is given by the positive and negative energy projections of the single-particle mean-position operator. This approach allowed us to observe the absence of ZB oscillations for the particle and anti-particle single-particle dynamical evolutions.

Acknowledgements.

We thank the Brazilian agencies CAPES (PROCAD2013 project), CNPq (#459339/2014-1), FAPEG (PRONEX # 201710267000503) and the Instituto Nacional de Ciência e Tecnologia-Informação Quântica (INCT-IQ) for partial support.

Appendix A Simulation for more spatial degrees of freedom

The goal of this Appendix is to show the simulation protocols for particles in two and three spatial dimensions.

A.1 Simulation of 2+1 dimensional Dirac Equation

Consider the Dirac Hamiltonian for a free particle existing in a 2D space

[TABLE]

The unitary transformation $e^{iS_{2}}=e^{-i\frac{\sigma_{1}p_{x}-\sigma_{2}p_{y}}{|\mathbf{p}|}\theta(\mathbf{p})}$ is the FWT which diagonalizes the Hamiltonian in this case, with $\theta(\mathbf{p})=\frac{1}{2}\textrm{tg}^{-1}\frac{|\mathbf{p}|}{mc}$. For the same reasons as in the 1D situation, there is no spin if the space is restricted to two dimensions and the particle state is a spinor with two components. Once we manage to construct FWT from optical device operators, the simulation protocol is made possible identifying again the transverse profile of a laser beam in the two orthogonal polarizations with the spinor components and identifying the two transverse coordinates with particle position. To show that in fact there is such a decomposition let us define a new momentum dependent phase

[TABLE]

motivated by the polar expression $p_{x}-ip_{y}=|\mathbf{p}|e^{-i\theta^{\prime}}$. In terms of the two phase functions the FWT reads

[TABLE]

Using the definitions given in the main text it is easy to see that the operator sequence $\hat{P}\left[-\theta^{{}^{\prime}}(\mathbf{p})\right]\hat{Q}\hat{P}\left[-\theta(\mathbf{p})\right]\hat{H}\hat{P}\left[\theta(\mathbf{p})\right]\hat{Q}\hat{P}\left[\theta^{{}^{\prime}}(\mathbf{p})\right]$ is equal to the FW unitary. The transformed diagonalized time evolution is then given by the product $\hat{H}\hat{P}\left[i\varepsilon(\mathbf{p})t/\hbar\right]\hat{H}\hat{P}\left[-i\varepsilon(\mathbf{p})t/\hbar\right]$. Just to conclude the protocol, the inverse FWT is given analogously by the same set of devices as the FWT with different imprinted phases and different angles for the wave plates.

A.2 Simulation of 3+1 dimensional Dirac Equation for a particular class of initial states

Let us consider the 3+1 Dirac equation

[TABLE]

with the standard choice of Dirac matrices $\beta=\begin{pmatrix}\mathbb{1}_{2\times 2}&0\\ 0&-\mathbb{1}_{2\times 2}\end{pmatrix}\quad,\qquad\alpha_{k}=\begin{pmatrix}0&\sigma_{k}\\ \sigma_{k}&0\end{pmatrix}$. The four components of the spinor $\psi(\mathbf{x},t)$ accounts for the two signs of the energy and for the two spin projections along a fixed direction. The FWT reads $e^{iS_{3}}=\cos{\theta(\mathbf{p})}+\beta\frac{\boldsymbol{\alpha\cdot}\mathbf{p}}{|\mathbf{p}|}\mathrm{sen}\,\theta(\mathbf{p})$ with the same definition for $\theta(\mathbf{p})$ as before.

A general simulation of the above equation, besides of requiring a four dimensional object to emulate the four spinor components, it would also require three spatial degrees of freedom, while the presented setup allows for just two. Instead of proposing a complete simulation, let us consider only the particular family of initial states given by

[TABLE]

In momentum space this state becomes

[TABLE]

which depedence on $p_{z}$ behaves like a $\delta(p_{z})$ for large values of $L$.

Thus we can approximate the FW transformed state by

[TABLE]

and all the momentum dependent phases only depend on two coordinates and can be aplied with SLMs.

The transformed state is explicitly written as

[TABLE]

we notice that the FWR only mixes the components two by two what makes possible to simulate it using two beams without any joint transformation between them. Then we can address the transverse profiles of horizontal and vertical polarizations of the first beam to $\phi_{1}$ and $\phi_{4}$ and the FWR as well as the subsequent diagonal time evolution do not mix this components with the two remaining. Moreover, each pair of mixed components transforms like the two spatial dimensions case (Eq. (8)) with the suitable phase signs.

The interesting thing about three spatial dimensions simulation is that it would enable us to investigate also spin effects like the spin analogous to Zitterbewegung.

Appendix B Simulation of 1+1 dimensional Dirac Equation for a particular class of potentials

Employing the strategy of Sabin12 , we show in this Appendix that, if the initial state is conveniently prepared, our approach is able to simulate the Dirac equation for a particular class of potentials. To this end, consider the 1D Dirac equation

[TABLE]

where $V(x)$ is a spinorial potential of the form

[TABLE]

For the above particular case, we define $\phi_{D}$ such that

[TABLE]

The substitution of this state on Eq.(14) shows that the spinor $\phi_{D}$ evolves according to the free Dirac equation (Eq.(1)). Thus, if we want to simulate the time evolution of the initial state $\psi_{D}(x,t=0)$, we need to prepare the state $\phi_{D}(x,t=0)=\exp\left[i\frac{1}{\hbar c}\int V_{1}(x^{\prime})\mathbb{1}dx^{\prime}\right]\psi_{D}(x,t=0)$ and the dynamics of the free evolution. From the experimental point of view, this corresponds to applying a local phase in position space to all the components of the spinor before performing the free evolution in the way as it is shown in the main text.