Simulation of Zitterbewegung by modelling the Dirac equation in Metamaterials

We develop a dynamic description of an effective Dirac theory in metamaterials, in which the wavefunction is modeled by the corresponding electric and magnetic field in the metamaterial. This electro-magnetic field can be probed in the experimental setup, which means that the wavefunction of the effective theory is directly accessible by measurement. Our model is based on a plane wave expansion, which ravels the identification of Dirac spinors with single-frequency excitations of the electro-magnetic field in the metamaterial. The characteristic Zitterbewegung is shown to emerge in simulations of the effective theory and we verify this signature with an analytic solution.


Introduction
Fundamental properties in quantum electrodynamics can be studied in strong field physics, in which for example intense laser beams can create electron-positron pairs in the unstable vacuum [1]. The study of these processes analytically [2] or numerically [3][4][5] demands extended efforts. Therefore modelling systems are discussed in the literature which allow to imitate physics in extreme conditions [6,7]. The capability of simulating intrinsic relativistic quantum processes like pair-creation or Zitterbewegung can therefore improve and test the understanding of quantum field theory.
The latter mentioned Zitterbewegung, first considered by Schrödinger [8], is a small quivering motion of particles in relativistic quantum mechanics. The theoretical existence of this quivering motion has been evidenced by numerical simulations of the Dirac equation and of quantum field theory [9][10][11]. However, the motion itself is very small and very fast. It can be estimated to be at the time scale of electron-positron pair creation (10 −20 s) and has an amplitude on the order of the Compton wavelength of an electron (10 −12 m). Therefore, the experimental observation of this counter-intuitive property of nature is challenging.
However, it is difficult to access the wave function itself in the listed experiments, where in particular higher order processes in quantum field theory are at least based on the knowledge of the wave function [37]. Recent experimental investigations demonstrated that the wave function of the Dirac equation can be imitated by the electromagnetic field of a waveguide structure with designable electrodynamic properties [38]. x y r z in the photonic crystal is adapted from an effective model for the electromagnetic field in a metamaterial waveguide which has been developed recently [38]. Here, the electric and magnetic field is propagating in the xdirection and ò 0 and μ 0 are the vacuum permittivity and vacuum permeability, respectively. According to the circuit model of the composite right/left handed transmission line, the electromagnetic parameters of our metamaterial obey the Drude dispersion relations [39,40]. In details, the frequency dependent relative permittivity and relative permeability are given by p C L a 1 1 d and 2 r 0 where we assume that damping of the electromagnetic field is negligible and the metamaterial is homogeneous along the x-direction in space. Damping terms may lead to a decay of electromagnetic excitations, which might be subject to future studies. We assume specific material parameters for our model which are listed in table 1. As a result, the electromagnetic metamaterial properties assume the specific relations

Effective Dirac equation
As discussed in [38], the Maxwell equations can be transformed to a one-dimensional, two-component Dirac equation by introducing a , 4 If one identifies the effective mass and the effective energy one can rewrite the Maxwell equations, such that they are of the same structural form as the one-dimensional Dirac equation Here, j is the two-component wave function is the vacuum speed of light and x s and σ z are the first and third Pauli matrices

Time evolution of Maxwell equations
In order to establish a dynamic Dirac theory, we first consider the time-evolution of the Maxwell equations (1).
Since the Maxwell equations are specified in frequency space, we evolve the electromagnetic field in the metamaterial by the time-evolution of an expansion of plane wave momentum eigenfunctions e i k x . Applying these plane waves yields the Maxwell equations in frequency-and momentum-space where the E z k , ( ) w and H y k , ( ) w are the electric and magnetic field amplitudes for a certain field frequency ω and a certain momentum k. From these Maxwell equations one can derive the dispersion relation The dispersion relation implies that there exists a positive and negative momentum k for each frequency ω. Furthermore, if one inserts the permittivity (2a) and permeability (2b), one obtains a polynomial of second order in ω 2 . This means that there exist two frequencies (upper band ω + and lower band ω − ) and the identical negative values −ω + and −ω − for each momentum k. This is plotted in figure 1 (a) for the two positive frequencies ω + and ω − .
Upper and lower band are divided by the band gap, in which the product of r r ( ) ( )  w m w turns negative. Outside of the band gap, the product r r ( ) ( )  w m w is always positive, and if also k 2 is positive (i.e. k  Î ), it follows that ω 2 is positive too and therefore .  w Î The upper and lower boundary of the band gap in frequency space is given by the frequencies, at which ò r and μ r are zero. This happens at the frequency 10.4 GHz, 12 at which ò r is zero and at the frequency Since the Maxwell equations are linear differential equations, the superposition of solutions is a solution again. Therefore, the time-evolution of the electromagnetic field can be written down by the expansion in which the sum runs over all discrete momenta k and all frequencies for each momentum. The factor k D is a measure for the summation in momentum space and is considered to be analogous to the measure dk of an integral. Since momentum space is spaced equidistantly in our approach, we write down a sum in (14) with measure Δ k and focus on the plane wave solution in this more fundamental theory part. The spacing of Δ k and its relation to position space is discussed in the section on numerical implementation 4.
In the expansion (14), the magnetic field can be expressed by the electric field by using one of the two Maxwell equations (10), yielding Therefore, the sum (14a) is already sufficient for a unique specification of the electromagnetic field in the metamaterial.

Time evolution of the effective Dirac equation
In order to write down the time-evolution in Dirac theory, the plane wave solutions of the Dirac equation (7) are to be considered. These are usually obtained by applying the plane waves e kx i at the Hamiltonian and solving the eigenvector and eigenvalue problem of the Hamiltonian matrix. The left-hand side of the Dirac equation (7) turns into the matrix when applying the plane wave e i k x . The definition of the mass (5) and energy (6) together with the dispersion relation (11) ensure the identity for the matrix (16). The solutions (18) are commonly known as bi-spinors of the Dirac equation in one dimension [41]. We have normalized the bi-spinors in our notion, such that they form an orthonormal base. Note that the spinors (18) have the identical analytic expressions Even though the eigensolutions of the matrix (16) are determined by (18,19), one has to keep in mind that this matrix is an integrated part of equation (7). Therefore, the eigenvalues and eigenvectors of the Hamiltonian matrix which is identical to (21a), except the sign of k. The k sign( ) factor in (21b) is necessary, because we want to keep the convention of the spinor solutions (18) as close to common relativistic quantum dynamics as possible (see also (A7)). Note that the fraction of upper and lower component of the spinors (21) corresponds to the factor (15a) between electric and magnetic field, if one accounts for the transformation law (4) of the electromagnetic field of the Dirac wave function. Therefore, the free eigensolutions of the Dirac equation are linearly dependent on the plane wave solutions of the Maxwell equations and it just remains to determine the proportionality constant to make a complete, unique link between the Maxwell equations of the electromagnetic field in the metamaterial and the corresponding dynamic Dirac equation. To do so, we write down the time evolution of the wave function by expanding it with respect to the obtained eigenfunctions u e is a normalization constant and Δ k is the measure for the momentum space summation, which is already mentioned in section 2.3 and further discussed in section 4. Since (22) is a unitary transformation of the initial state of the wave function, the normalization constant will keep constant in time and x t , ( ) j will always stay normalized at one. Note that according to the considerations above, the positive eigensolutions which are proportional to u k + will evolve with the lower band frequency ω − and the negative eigensolutions which are proportional to u k will evolve with the upper band frequency ω + .
If one compares the expansion of the electric field (14a) with the first component of the wave function (22) and requires equality for the prefactors of the exponential e kx i w if one requires equality for the prefactors of e .

The metamaterial's scaled Dirac equation
According to the considerations in the sections 2.3 and 2.4 the metamaterial's emulated quantum dynamics is expected to evolve as in Dirac theory, but natural constants like the electron mass and the effective speed of light depend on the metamaterial properties. For deducing these properties, we first consider the time-evolution equation in relativistic quantum mechanics. Here, m 0 is the electron rest mass and is the momentum operator. After a Fourier transformation of this equation into frequency space, the time- If equation (7) could be cast in a similar shape, such that the prefactor of the function x, ( ) j w was a linear function of ω at the right-hand side, its generalization to the time-domain would be straight forward and the mass and speed of light parameters of the emulated Dirac equation could be read from the corresponding Hamiltonian on the left-hand side. This can be achieved by a Taylor expansion of the effective energy ( )  w and effective mass m ( ) w of the metamaterial Dirac equation (7), at the frequency , 0 w at which ( )  w is zero. Thus, ω 0 is determined by setting the definition (6) of ( )  w to zero and solving for ω, resulting in For the parameters chosen in table 1 ω 0 has the value 11.0 GHz. The first non-vanishing order of the Taylor expansion is sufficient for our consideration. For the mass m , ( ) w the first non-vanishing order is the zeroth order, whereas for the effective energy , ( )  w the first non-vanishing order is the first order. The reason is that ω 0 is chosen such that the zeroth order of ( )  w vanishes. Accordingly, one obtains for equation (7). If one defines the shifted frequency This is consistent with the Hamiltonian in equation (28) and a negative time evolution, according to equation (25). Comparing the corresponding Hamiltonian with the Hamiltonian of Dirac theory (26), one deduces the substitution which scales the Dirac theory of electrons such that it is similar to the emulated Dirac dynamics of the metamaterial.
We point out that equation (34) implies a negative time-evolution, because its right-hand side contains a minus sign in front of the  wj D expression. In a Dirac equation with positive time-evolution, this minus sign would be a plus sign (see equation (28)). The minus sign originates from our chosen convention (4, 5, 6). A different convention might require an inversion of the x-axis. Therefore, if we want to compare the metamaterial dynamics with results from exact Dirac theory, we have to reverse time. In particular we have to add a minus sign in the time argument of the semi-analytic expression (A16) and (A17) of the Zitterbewegung, which is derived in appendix A.

Boundary conditions
In this section, we consider electrodynamic boundary conditions, which are used for the injection of electromagnetic excitations at the boundaries of the metamaterial. The boundaries are at the positions x a and x b , such that the physical space is divided into : region 1 , the left input wave guide, : region 2 , the metamaterial, : region 3 , the right input wave guide. 37 The boundary conditions are obtained by integration of equation (1) We make the reasonable assumption that the right-hand side of equation (1) is bound. In this case the integrals over the intervals with infinitely small ò imply that the electric and magnetic fields must be smooth at the boundaries x a and x b , which means that We assume that the relative permittivity and permeability are just 1 for the propagation along the input wave guides of the regions (1) and (3), like it is for electromagnetic waves in vacuum. Then, the dispersion relation (11) simplifies to c k .
the field expansion (14) can be simplified to (1) and (3). Similarly, for right propagating waves, in which E H 0 (14) can be simplified to for the regions (1) and (3). The simplified expansions imply a dispersionless translation of the signal in the regions (1) and (3). In other words: once the left and right propagating input and output signals are known in the regions (1) and (3) at any position x, they can be trivially deduced from the above equations. It is most convenient to know the input signal at the metamaterial boundaries x a and x b . Therefore, we introduce the new coordinates x x x a for region 1 and 43 in terms of these new coordinates. If one inserts the field expansion (14) in these continuity conditions and keeps in mind that the exponentials e t iw are linearly independent on the time interval , ] [ -¥ ¥ for each frequency ω, one obtains the boundary conditions in frequency space The relations (15) allow for the substitution of the magnetic field expansion coefficients in (45c) and (45d) with the electric field expansion coefficients, resulting in is an abbreviation for convenience.

Numerical implementation
The numerical implementation makes use of an equidistant grid of states in momentum space, which allows for a discrete Fourier transformation of the grid points. The corresponding position space grid points are therefore equidistant as well. Assume, there are n grid points, each of them spaced by the extension x 8 mm D = of the metamaterial's unit element size (see table 1). This implies a length l of the metamaterial which is l n x.
= D Since the plane waves e kx i should be 2π periodic after this extension l, the spacing in momentum space must be k l 2 . [ ] -Furthermore, the simulation has periodic boundary conditions at the interval limits, because all exponentials are periodically continuing. Therefore, the simulations in the sections 5.1 and 5.2 are valid, as long as the wave function's probability density stays within the boundaries of the simulation area.
In section 5.3 the initial condition of the simulation is not specified at an initial time, but imposed by additional boundary conditions instead. These boundary conditions, which are discussed in section 3, are modelling the injection of input pulses at the physical boundaries of the metamaterials. However, the positions x a and x b , at which the injection boundary conditions are imposed are within the simulation area (see illustrative figure 2), such that they don't interfere with each other.
The time-evolution is computed according to equation (22), which is implemented numerically by first multiplying the product of each expansion coefficient f k and spinor u k with its time-evolving phase factor e t iw at the certain time t and then summing over k. Since each summand of the summation contains the factor e kx i and the discrete set of momenta k is equidistantly spaced, the sum is executed by performing a Fourier transformation of the upper and the lower component of the u e k k k t i ( ) f w array.

Results from simulation
In this results section, we apply the theory considerations from above in order to demonstrate a Zitterbewegung of the emulated Dirac dynamics in the metamaterial. In a series of three different simulations, we first consider a Gaussian wave packet excitation as the simplest possible setup in section 5.1. For easier experimental implementation we also consider a moving Gaussian wave packet in section 5.2. In the last section 5.3, we implement boundary conditions and account for the influence of these boundaries in the simulation.

Simulation
The simulation is carried with a resolution of 401 grid points in momentum space. Therefore the simulation area has an extension of 3.2 m from the 1st till the 401th index, according to the numerical considerations section.
We have plotted the probability density of the metamaterial simulation with initial condition (48) in figure 3, in which we also compute the position expectation value Note that we have inserted a minus sign in the argument of the sine function by hand, according to the negative time-evolution of the metamaterial simulation which we have discussed at the end of section 2.5.

Initial condition
The Gaussian wave packet excitation in the above subsection is entering and exiting the simulation region very slowly, which means that the required electromagnetic input pulse at the metamaterial interfaces will have an infinitely long head and tail. This is no useful property for an experimental realization and therefore, we demand a wave packet which performs a Zitterbewegung but evolves more suitably in a prospective experiment. This property can be achieved by shifting the Gaussian wave packet (48) to the right in momentum space by the value k 20 m , Here, the width of the wave packet in frequency space is chosen to be m c .  ñ of (51), which is also shown as solid line in figure 6  (a). (b) The time-evolution of the same initial condition, but simulated with exact Dirac theory. The corresponding position expectation (red line) is also plotted as black dotted line in figure 6 (b). Like in figure 4, the initial condition for the simulation is specified at t 0 ns, = from which the backward and forward time-evolution is computed. might expect that positive and negative states move in the same direction, because we have shifted the excitation by momentum k 0 to the right in momentum space. But since the upper and lower band of the dispersion relation have the opposite slope at k 0 the shift in momentum space results in two counterpropagating wave packets. The metamaterial simulation of the effective Dirac equation and the simulation of the exact Dirac equation differ significantly from each other. Therefore, we plot them in two different plots (a) and (b) in figure 5, respectively. In the case of the effective Dirac theory of the metamaterial, there is an asymmetry which we explain with the asymmetry of the dispersion relation (see figure 1 (a)). The positive energy eigenstates are moving with group velocity at k 0 of the lower band, while the negative energy eigenstates are moving with group velocity at k 0 of the upper band. Since the upper and lower band are differently shaped, the positive and negative energy-eigenstates are propagating at different speed through the medium. On the other hand, in the case of the exact Dirac equation in figure 5 (b), the dynamics appears symmetric which we explain with the symmetric upper and lower band of the relativistic energy-momentum relation in figure 1 (b). As a result of the differences of the dispersion relations, there is an additional effective motion superimposed to the Zitterbewegung (see red line in figure 5 (a) and solid black line in figure 6 (a)), as compared to the wave packet in the exact Dirac theory (see red line in figure 5 (b) and dotted black line in figure 6 (b)).
Similarly as in the above section 5.1, we want to verify the exact Dirac theory by comparison with the semianalytic solution in appendix A. The new initial condition (53), substituted in equation (A17) yields Note that we inserted a minus sign by hand in the argument of the sine function, as it has been done for equation (52) already. We plot the value of the solution (55) as plus-marked line in figure 6 (b).

Discussion
The lines 'exact Dirac' of the exact Dirac simulation and the 'analytic' solution in figure 6 (b) are coinciding. We are expecting this property, because the analytic solution is based on exact Dirac theory. Therefore, the coincidence is a mutual validation of both descriptions. For a comparison of the metamaterial simulation with the exact Dirac theory, we subtract the value t 49.7 ns ( )from the position expectation value in figure 6 (a), and plot it as the solid black line 'modified simulation' in figure 6 (b). We obtain very good agreement of both descriptions and take this match as further confirmation for the reliability of the simulation of the effective Dirac dynamics in the metamaterial.

Initial condition
In the case of a metamaterial with physical boundaries, the initial wave function (54) has to be replaced by an electric and magnetic field which is propagating from the regions (1) and (3) into region (2) of the metamaterial. This can be done by specifying the incoming electric field of the regions (1) and (3) at the boundaries by which according to Parseval's theorem, is equivalent to normalization in position space (23). In contrast to the initial condition at time t = 0 along the whole x-axis in the above two subsections, this most realistic simulation specifies the 'initial condition' at the boundary positions x 1.0 m a =for (56a) and x 1.0 m b = for (56b) along the time-axis. However, through the relations (40,41,42), the alignment of the 'initial condition' can be recast from the time-axis to the x-axis regions (1) and (3) at a certain time t 0 ns.  . Therefore, even though the time-evolution has periodic boundary conditions, the wave packet does not immediately reenter the simulation region (2) from the other side, after it went out before.

Discussion
The simulation of a Zitterbewegung with real metamaterial boundaries is shown in figure 7 and the position expectation value of this simulation is overlayed as red line and plotted a second time in figure 8. A similar simulation for the exact Dirac equation with the same initial condition differs much from the effective Dirac dynamics of the metamaterial (see also figure 5) and avoids the feasibility of a direct comparison. Therefore, there is no comparison with the corresponding dynamics of the exact Dirac equation in this subsection. In figure 7 one can see that the electromagnetic excitation in the setup is propagating from the physical boundaries towards the central region of the simulation area. In the overlap region of both excitations an oscillatory pattern emerges (see figure 8) which we attribute as Zitterbewegung of the effective Dirac equation in the metamaterial. This demonstrates that in the case of an ideal metamaterial, given by the equations (1) and (2) a Zitterbewegung of a simulated Dirac equation may occur for a real experiment.

Conclusion and outlook
The article discusses the feasibility on an exactly mappable, time-dependent Dirac simulation which is actually emulated by electrodynamics of Maxwell equations with designed electromagnetic properties of a metamaterial structure. We present a time-evolution equation in a Dirac-like fashion (22) and give arguments of why it evolves very similarly to exact Dirac theory. Furthermore, we consider scaling relations in section 2.5, with which the slower dynamics of the metamaterial can be mapped to the faster dynamics of elementary particles in Dirac theory. Based on the scaling rules and the exact Dirac equation, we have derived an explicit, semi-analytic expression for the position expectation of a general wave function in appendix A.
We demonstrate the dynamic description with three simulations: in sections 5.1 and 5.2, we compare our metamaterial simulations with exact Dirac theory and with the analytic solution which is derived within this exact Dirac theory. We get very good agreement between the emulated Dirac dynamics and the exact Dirac dynamics, which we interpret as validation of our description. It is also a proof for the emergence of Zitterbewegung, which is a characteristic property of the free Dirac equation. Subsequently, we present a more  realistic scenario, in which boundary conditions are included, such that the simulation can be realized in the experiment in subsection 5.3. In this most realistic simulation we also observe an oscillatory pattern which matches the Zitterbewegung.
The differences between the effective metamaterial description of the Dirac equation and the exact Dirac equation can be listed as follows: • The dispersion relation of the metamaterial description is only an approximation of the hyperbolas of the relativistic energy-momentum relation in exact Dirac theory (see figure 1). The differences emerge in quantum dynamics, which can seen distinctly in the comparing figures 5 and 6 of section 5.2. The differences between effective Dirac theory and exact Dirac theory get large for high momenta.
• The mass of the effective Dirac theory is frequency-and wave-length dependent. For the right choice of material parameters this frequency dependence can be small. On the other hand, in contrast to exact Dirac theory the frequency dependent mass may even change its sign in certain metamaterial configurations.
• Since the mass depends on frequency and wave-length, also the plane wave spinors (18) are slightly deviating from the spinors in exact Dirac theory.
• The wave-function, which is emulated by electrodynamics in the metamaterial has no imaginary part, in contrast to exact Dirac theory. However, the relation between a real-valued and a complex-valued wavefunction in the metamaterial is unique, which is explained in appendix B.
• In this work we treat the wave-function as classical field. In contrast, the exact Dirac theory may be quantized. The quantization properties of the metamaterial system are subject to future studies.
We conclude that a Dirac wave function can be emulated by the electromagnetic field in a metamaterial, such that the Zitterbewegung occurs. Next steps are an experimental implementation of this theoretical study of the effect and an investigation of the quantization of the system.  (24) in momentum and frequency space. Therefore, the question why the electromagnetic field does not contain an imaginary part in reality, but that in contrast the Dirac wave function does contain an imaginary part is still to be answered! The difference of the description of the Maxwell equations and the Dirac equation can be found in the expansions (14) and (22). While the sum in (22) only runs over the positive frequencies ω + and ω − , it also runs over the negative frequencies −ω + and −ω − in (14). Therefore, on one hand, the simulated Dirac equation only makes use of the upper two bands of the dispersion relation (11) and is consistent with the two bands of the positive and negative energy eigenstates of the effective Dirac equation (34) (see also figure 1). On the other hand, the dispersion relation (11) has four bands, of which the two negative ones just have the negative value of the positive ones. So the expansion of the electromagnetic field runs additionally over these two negative energy bands. As a conclusion, the time-evolution of the Dirac equation (22) omits the sum over the negative energy bands of the electromagnetic field, without explicitly mentioning it. This omittance is introduced for being consistent with the two bands of conventional Dirac theory. However, the real-and imaginary part of the expansion coefficients E , z k , ( ) w -H y k , ( ) w of the negative energy band of the Maxwell equations can be chosen according to the symmetry scheme of table B1, such that the electric and magnetic field which is transformed from the wave function by the relations(24) becomes exclusively real. Elimination of the imaginary part by the above discussed symmetry procedure and division by a factor of two is equivalent to taking the real part of the wave function. Therefore, the considerations of this appendix are justifying that the real part of the simulated Dirac wave function is proportional to the electric and magnetic field of the underlying Maxwell equations, which can be measured in experiment.  (14) with respect to frequency ω and momentum k, in analogous extension to the symmetry considerations of the expansion (B1). The first column 'coefficient' tells if the real part or the imaginary part of the wave function is considered, the second 'frequency' and third 'momentum' column tells, whether the real/imaginary part of the expansion coefficients are assumed to be fully symmetric or fully anti-symmetric with respect to frequency or momentum, respectively. The last line 'function' tells if the resulting function of this expansion will be exclusively real or exclusively imaginary.