A Study of Space-Time Discretizations for the Dirac Equation

A Study of Space-Time Discretizations for the Dirac Equation

Robert Vaselaar, Hyun Lim, Jung-Han Kimn Department of Mathematics and Statistics, South Dakota State University, Brookings, SD 57007

We study several numerical discretization techniques for the one-space plus one-time dimensional Dirac equation, including finite difference and space-time finite element methods. Two finite difference schemes and several space-time finite elements function spaces are analyzed with respect to known analytic solutions. Further we propose a finite element discretization along the equations’ characteristic lines, creating diamond-shaped elements in the space-time plane. We show that the diamond shaped elements allow for physically intuitive boundary conditions, improve numerical efficiency, and reduce the overall error of the computed solution as compared to the other finite difference and space-time finite element discretizations studied in this paper.

I Introduction

The Dirac equation governs all spin particles, known as fermions. While solutions to the Dirac equation may be used to derive quantifiable predictions of particle behavior from quantum physics, these solutions are sometimes difficult to find in experimentally interesting scenarios. Numerical methods for the Dirac equation may be able to bridge some of the gap between theoretical and experimental particle physics fillion2012 ; widom1996neutrino .

These include many finite difference based methods for lattice quantum chromodynamics, such as the Dirac Wilson equation PhysRevD.10.2445 which have been used in conjunction with modern numerical methods such as Krylov subspace solvers Sakurai2010113 ; Nakamura201234 and scalable additive Schwarz preconditioners Luscher2004209 . The limitations of these methods are also well known and include the inability to reconcile all limitations simultaneously. This is particularly important when considering the problem of fermion doubling, a condition where the number of particles considered must naturally double for each space-time dimension included on the lattice, and chiral symmetry, which is usually broken by most numerical methods that prevent fermion doubling Nielsen1981219 ; Chandrasekharan2004373 .

Other numerical methods for the Dirac equation include radial formulations created to investigate the energy spectrum of heavy atomic ions almanasrehgconverg ; kullie2004 ; kullie2001 ; Desclaux2003453 . These are based on the Dirac equation in the presence of a coulomb potential resulting a relativistic eigenvalue problem, using both finite difference and finite element numerical methods.

The finite element method has also been used to calculate the propagation of free fermions in space. Analysis of the finite element method combined with Crank-Nicholson time stepping scheme demonstrates that solutions may show inconsistent and impossible physical behavior, such as superluminal propagation, depending on the step size and propagation method used muller1998 . Since using Lagrangian interpolation elements in one dimension are algebraically similar to using finite differences, it is also natural that this choice of function space has the same problems of fermion doubling and numerical instability as its finite difference relativemuller1998 . In this implementation, physically consistent behavior of the particle depended on particle momentum, finite element size, and time step size chosen.

In this paper, several implicit space-time discretizations based on the finite difference and Galerkin methods are presented. This presentation will show that simulation behavior is directly affected by choice of discretization method and function space. The problem domain is then rotated by in the space-time plane, forming diamond-shaped tensor elements, and the solution is recalculated using the rotated domain. This rotated domain shows substantially reduced error and improved performance when compared to the other space-time discretizations listed here. The goal of this research is to create a discrete form of the Dirac equation that shows good agreement with the analytic solution as well as low error and the absence of faster-than-light propagation. Further, we would prefer a solution that does not modify the original Dirac operator in order keep as many of its original physical properties as possible.

This paper is organized as following; First we present the weak form of the gauge-free Dirac Equation in section II. In section III two space-time finite difference and one finite element method are presented along with their numerical results to observe their behavior and performance. Then in section IV three finite element discretizations using space-time tensor elements are presented along with their numerical results. Sources of possible simulation error are also presented and analyzed. Proceeding from the discussion of error we propose our diamond-shaped approach in section V and show how this approach addresses the errors observed and improves simulation efficiency. We conclude by discussing future research opportunities in Section VI.

Ii Space-Time Methods

ii.1 Weak Formulation of the Dirac Equation

The one dimensional Dirac operator may be expressed as follows


where are the usual Pauli matrices defined as


Here the Pauli matrices are chosen such that the variables and form a Minkowski space-time, which an essential relationship in the Dirac equation.

In this case we will consider the initial value problem given by


In the gage-free case analytic solutions may be computed directly which will give us a basis for comparison for our numerical results.

Using the continuous time Galerkin method the weak form may be expressed as follows. The objective is to find such that


Iii Numerical Results of Different Numerical Approaches

In this section we will show some results from two finite difference approaches and one finite element method to the Dirac equation. These are the central difference method, the staggered finite difference formulation, and the finite element method using triangular finite elements. In each method we observe significant non-physical effects in the space-time boundary value problem introduced previously. For the sake of comparison, we refer the reader to figure 1 which shows the analytic solution to the space-time boundary value problem proposed above.

(a) Real Component of
(b) Imaginary Component of
Figure 1: Analytic Solution of the Massless Initial Value Problem

iii.1 Central Difference Discretization

For an implicit implementation of the two-dimensional Dirac equation using the finite difference method, the integral used in the bilinear form above may be replaced with a double summation


Often referred to as the naive discretization, the matrix in this approach is built according to equation 7. Here we use the following central difference definitions for the partial derivative operators.


The central finite difference discretization was implemented using the bilinear form of the finite difference method shown in equation 7 and the initial value was introduced via a matrix partitioning scheme. The result of the balanced difference discretization when applied to the dimensional Dirac initial value problem is shown by Figure 2 and Table 1.

(a) Real ,
(b) Imaginary ,
(c) Real ,
(d) Imaginary ,
Figure 2: Central Difference Solution. Note that in Area 1 and Area 2, the wave function has shifted rightward, indicating super-luminal propagation, which is forbidden for massless solutions of the Dirac equation

Figure 2 shows that the wave function is similar to the analytic solution when the space and time step sizes are equal. However, when unequal step sizes are used the wave shape deteriorates and is shifted rightward, indicating speeds in excess of the speed of light, which is physically impossible.

Mesh Size Matrix Size BICGSTAB Iterations Error
Table 1: Numerical Performance of Central Difference Discretization

Table 1 shows that norm of the error initially initially improves with a finer mesh, but does not improve uniformly and does not appear to tend toward zero with finer mesh spacings. Further, when the spacing is unequal, , the error is substantially larger, which is expected due to its non-physical behavior.

iii.2 Balanced Difference Discretization

When used to create an explicit propagator, the central difference discretization does not necessarily conserve the probability current of the wave-function. To address this shortcoming the partial derivative stencil in equation 9 is replaced by stecils that are arranged symmetrically with respect to space and time as follows wessels1999 .

Unlike the original paper wessels1999 , where this discretization is used to construct an explicit propagator, our implementation is fully implicit in both time and space.

(a) Real ,
(b) Imaginary ,
(c) Real ,
(d) Imaginary ,
Figure 3: Balanced Difference Solution. Note that in Area 1 and Area 2, the wavefunction has shifted leftward, indicating sub-luminal propagation, which is inconsistent with the expected behavior of massless solutions

Figure 3 shows that while the wave shape is choppy, when tested with equal time and space step sizes, it holds a continuous pattern in the overall shape of the analytic solution and the solution shows the correct propagation speed of . However, with unequal space and time step sizes the propagation speed is visibly slowed to .

Mesh Size Matrix Size BICGSTAB Iterations Residual Error
Table 2: Numerical Performance of Balanced Difference Discretization

iii.3 Triangular Lagrangian Elements

Figure 4: Local Grid Square of a Triangular Finite Element Discretization

Triangular Lagrangian finite elements represent one of the most commonly used shapes in the finite element method. In this case the nodes of the discrete Dirac equation are arranged on a grid pattern, so each square is composed of two triangular elements as shown in figure 4. Assuming the single grid element is a unit square with local variables and , the interpolation polynomials on triangle are given by

similarly, interpolation polynomials for are given by

Evaluating the finite element the integral in equation 6 establishes an algebraic relationship between the nodes that for non-boundary elements is equivalent to a finite difference stencil. The finite difference stencil for triangular Lagrangian elements are calculated to be

The columns of the matrices above correspond to the spatial dimension and the rows correspond to the temporal dimension . It is apparent that the stencils above are not symmetric with respect to space and time. This means that the choice of element shape may bias the finite difference stencil along the characteristic line or , depending on which triangle orientation is chosen.

(a) Real ,
(b) Imaginary ,
(c) Real ,
(d) Imaginary ,
Figure 5: Triangular Lagrangian Element Solution. Note that the calculated wave fades out quickly and propagates in the wrong direction, compared with the analytic solution.

Figure 5 shows that instead of instability, first-order Lagrangian space-time finite elements lose wave amplitude very quickly and appear over-damped. From a physical perspective, the particle is disappearing into space. While the precise reason for this disappearance is unknown, it may be due to the finite difference stencils’ bias in the opposite direction of particle propagation. When unequal step sizes in space and time were tested, the wave function begins to move to the left in the opposite direction of the analytic solution and at a speed greater than the speed of light.

Mesh Size Matrix Size BICGSTAB Iterations Error
Table 3: Numerical Performance of Triangular Lagrangian Elements

Iv Tensor Element Based Approaches

In this section we will show three space-time discretizations that use sqaure shaped ”tensor” finite elements along with a selected basis function space to form the weak Dirac boundary value problem Each of these approaches shows overall convergence to the shape of the analytic solution without the presence of superluminal, subluminal, or counter-directional wave functions that were present with the previous approaches.

iv.1 Polynomial Hermite Tensor Elements

These functions are conceptually related to the third-order piecewise Hermite interpolation polynomials given by

Figure 6: Hermite basis functions where .

We use the tensor product to create a composite function that has continuity in a two dimensional plane, which is required for analytic solutions in quantum mechanics. This set also allows the set of second-order mixed partials to be varied independently. However, the continuity of mixed partials has no special physical significance in this case, so these functions are removed in order to reduce the degrees of freedom present in the discrete system.

(a) Real() Hermite Elements Plot
(b) Real() Error Plot
(c) Real , Hermite Elements Plot where
(d) Real , Error Plot where
Figure 7: Solution Curve of Polynomial Hermite Tensor Element

Figure 7 shows that when Hermite tensor elements are used as a function space for the given initial value problem, the overall behavior of the wave-function is consistent with the analytic solution both when and when , although the error function is substantial.

Mesh Size Matrix Size BICGSTAB Iterations Error
Table 4: Numerical Performance of Hermite Tensor Element Discretization

iv.2 Trigonometric Hermite Tensor Elements

Here we will chose our basis functions for the finite element vectors to be trigonometric functions given by

These are conceptually similar to the Hermite polynomials. As previously, we use the tensor product to create a composite function that has continuity.

(a) Real() Trigonometric Elements Plot
(b) Real() Error Plot
(c) Real() Trigonometric Elements Plot where
(d) Real() Error Plot where
Figure 8: Solution Curve of Trigonometric Tensor Element

From viewing the solution and error plot in figure 8 we can see that although the overall shape of the solution is very similar to the analytic solution the error wave is relatively large.

Mesh Size Matrix Size BICGSTAB Iterations Error
Table 5: Numerical Performance of Trigonometric Tensor Element Discretization

iv.3 Linear Lagrangian Tensor Elements

The function space of the linear Lagrangian elements is formed from the tensor product of the first order Lagrangian interpolation polynomials in the and directions. Since they are first order, there is only one degree of freedom per node, or four degrees of freedom per tensor element, making this element type much simpler than the previous elements shown.

The Lagrangian interpolation polynomials are given by the following expression.

(a) Real() Lagrangian Tensor Elements Plot
(b) Real() Error Plot
(c) Real() Lagrangian Tensor Elements Plot where
(d) Real() Error Plot where
Figure 9: Solution Curve of Lagrangian Tensor Element

From Figure 9 we see that the Lagrangian tensor elements also produce a numerical result very close to the analytic solution. This behavior is consistent both when and when .

Mesh Size Matrix Size BICGSTAB Iterations Error
Table 6: Numerical Performance of Lagrangian Tensor Elements

Table 6 shows that the Lagrangian tensor elements have substantially lower error, smaller matrix size, greater mesh refinement and more efficient convergence than either the Hermite or trigonometric tensor elements. Unfortunately, the norm of the error actually increases with greater mesh refinement. Possible sources of this remaining error will be analyzed in the following section.

iv.4 Error Analysis of Linear Lagrangian Tensor Elements

Figure 10: Error Analysis of Lagrangian Tensor Element. In Area 1 we see boundary error that propagates inward from the natural boundary conditions on the right, left, and rear edges of the domain. In Area 2 the error wave is composed of closely spaced peaks of period 2h, which are effectively invisible to the partial derivative operator for non-boundary nodes

From looking at the error wave in figure 9 and the analysis in figure 10 we make the following observations about the sources of error. The first source of error for Lagrangian tensor elements may lie in the stencil of the and operators. For non-boundary nodes, linear Lagrangian tensor elements introduce the following algebraic relationship between nodes.

It is apparent that the value of the partial derivative as calculated by these operators would approach zero as the period of the wave approaches . The error wave across the right-hand side appears to have a period of exactly , making it ”invisible” to the discrete form of our partial derivative operator.

A second source of error may come from the boundary conditions across both sides and . Here an error wave springs ex-nihilo from the side and propagates parallel to the solution. It may be possible to eliminate such waves by choosing Dirichlet boundary conditions. However, these conditions would imply knowledge solution before the solution is calculated. The source of this problem is that the domain sides are not completely contained by the light-cone of our initial condition.

If the domain were sufficiently wide as to preclude the wave from reaching the boundary, it would then be appropriate to apply Dirichlet boundary conditions to the sides of the experiment. However, this would also add siginificant empty space to the domain and computational cost to the experiment.

From the error observations above we draw the following conclusions. One, boundary conditions should utilize light cone causality to ensure a unique solution. Two, momentum and energy operators should be able to detect tightly spaced, erroneous wave patterns and prevent them from appearing in the solution..

V Diamond Shapend Tensor Elements

In order to reduce the error waves observed in the Lagrangian tensor element discretization, we propose the following element shape for discretizing the 1+1 Dirac equation, shown in Figure 11. The plane is rotated to create two new axis that will will name ”right” and ”left”. This shape gives us two important advantages.

One, we may impose Dirichlet boundary conditions across both the right and left axis. The entire domain is then contained within the light cone of the ”initial” conditions, meaning that the solution should be unique, at least from the physical perspective, since no new information can enter the domain. This is shown graphically in figure 11(a).

(a) Diamond Shaped Space-Time Domain
(b) Diamond Shaped Space-Time Element
Figure 11: Diamond Shaped Domain and Single Element Composition

Linear interpolation polynomials for the diamond tensor elements may be defined as follows.

Two, the partial derivative stencils now become more complex and should be better able to detect the closely chopped error waves that were present in the Lagrangian tensor element solution. For non-boundary elements, these linear interpolation polynomials introduce the following algebraic relationship between nodes for the two partial derivative operators of the Dirac equation.

Using this domain, element shape, and interpolation polynomial set with the weak form defined in equation  6 generates the following solution shape shown in Figure 12. The magnitude of the error wave is extremely small when compared to the solution, and shows that the finite element solution is nearly exact when one considers the values of the wave-function at the node points. The scale of the error at the nodes is around twelve orders of magnitude lower than the error at the nodes for other methods considered.

(a) Real() Solution Plot
(b) Real() Error Plot
Figure 12: Solution Curve of Diamond Lagrangian Tensor Element. Note the reduction in error scale compared with previous methods tested
Mesh Size Matrix Size BICGSTAB Iterations Error
Table 7: Numerical Performance of Diamond Tensor Elements

We note that the solution shows substantially lower error when compared to any of the methods previously presented. This is further confirmed by our results in  7. Here we see that the norm of the error is much lower than for the other methods tested, and that the numerical simulation converged more quickly as well. In the case of a the element matrix the diamond tensor element solution converged in iterations versus iterations for the Lagrangian tensor elements, and the norm of the error was (Table  7) for the diamond tensor elements versus (Table  6) for the Lagrangian tensor elements.

Finally, as with the other space-time tensor element approaches, no superluminal or subluminal behavior was observed when tested with unequal space and time spacings.

v.1 Rotation Tests

To test the effectiveness of other possible domain rotation angles, the domain was rotated about the origin counter-clockwise from to . This is shown conceptually in Figure  13 demonstrating how the domain rotates about the origin of the space-time plain.

Figure 13: Conceptual Diagram of a Domain Rotation in Space-Time

After performing this rotation on the domain, a similar initial value problem as the one given in  6 was then solved on the new rotated domain. The algebraic formulation rotated domain is given below.

In this test, and , where is the angle of rotation. The -axis is scaled by a factor of so that the wave function remains centered in the domain throughout the rotation. This rotation has the added advantage of following the path of the solution more closely, since high-energy, low mass solutions to the Dirac equation tend to move along the characteristic lines of equation; which is to say, particles that have high energy and low mass move at nearly the speed of light.

This test was also conducted with non-zero masses, and results compared to a solution calculated using numerical Fourier transformation. Due to stability concerns in the massive case, a Dirichlet boundary condition was added to the side and the the center of the wave function was moved from to to keep the wave function from colliding with the Dirichlet boundary condition on the wall. The domain shape was slightly altered to .

The number of GMRES iterations and the norm of the error were recorded and plotted against the rotation angle used. The results are shown in Figure  15 and  14.

Figure 14: , the Error Norm of the Computed Wave Function vs the Angle of Rotation
Figure 15: GMRES iterations to reach a residual of vs the Angle of Rotation

In Figure  14 the norm of the error is shown to decrease steadily as the angle is increased from to and reaches its lowest point at . This relationship is demonstrated in each case tested, regardless of particle mass.

Likewise, in Figure  15 the number of GMRES iterations is shown to decrease steadily as the angle is increased from to and reaches its lowest point at . This relationship is demonstrated for each domain size tested but only in the massless case.

Interestingly, this relationship is changed somewhat as the mass increases. At a mass of and , the fewest iterations are used at around . When , the fewest iterations are required at . Finding the source of this small off-angle efficiency improvement requires further investigation.

These results demonstrate that a rotation in space-time shows the lowest error for all angles and mass levels tested. Further, they also demonstrate that a rotation is either optimal or nearly optimal from a GMRES iterations perspective; however, this relationship is less straightforward than the correlation with error and will require further research to establish the relationship between the massive particle and the optimal rotation angle for algorithmic efficiency.

Vi Conclusion

From the data above we have shown several space-time approaches that may be useful in numerical calculations of the Dirac equation in a rectangular space-time domain. We have also shown that a physically-motivated selected of element and domain shape can substantially improve performance and reduce error for the numerical experiments considered above. Since this improvement was shown for an equation dominated by first-order operators, it may be possible to use a similar approach for other equations with unstable first-order operators as well.

The results above also show that problems with super-luminal and physically inconsistent propagation may be addressed by choice of discretization and using a fully implicit method. This is then corrected without reference to the problem of Fermion doubling as was suggested by Müller et al. in muller1998 .

To expand the usefulness of this numerical approach, further research should be conducted in several areas. One, the sample problem should be developed for and dimensional settings. Two, a more formal error analysis should be conducted to understand the root cause of the simulation behavior above.Three, scalable preconditioners should be investigated for new numerical solutions to the Dirac equation, especially given the size of Dirac-based problems in dimensions (or more). Finally, this model should be tested for suitability and performance in more realistic, inhomogeneous or nonlinear settings.

Acknowledgements We thank you to Professor Dongming Mei from the Department of Physics of University of South Dakota for his discussion on physics interpretation on this results.


  • (1) Andre D. Bandrauk Francois Fillion-Gourdeau, Emmanuel Lorin. Numerical solution of the time-dependent dirac equation in coordinate space without fermion-doubling. arXiv:1107.4650v2, 2012.
  • (2) A Widom and YN Srivastava. Neutrino flavor oscillations using the dirac equation. arXiv preprint hep-ph/9608476, 1996.
  • (3) Kenneth G. Wilson. Confinement of quarks. Phys. Rev. D, 10:2445–2459, Oct 1974.
  • (4) T. Sakurai, H. Tadano, and Y. Kuramashi. Application of block krylov subspace algorithms to the wilson–dirac equation with multiple right-hand sides in lattice qcd. Computer Physics Communications, 181(1):113 – 117, 2010.
  • (5) Y. Nakamura, K.-I. Ishikawa, Y. Kuramashi, T. Sakurai, and H. Tadano. Modified block bicgstab for lattice qcd. Computer Physics Communications, 183(1):34 – 37, 2012.
  • (6) Martin Lüscher. Solution of the dirac equation in lattice qcd using a domain decomposition method. Computer Physics Communications, 156(3):209 – 220, 2004.
  • (7) H.B. Nielsen and M. Ninomiya. A no-go theorem for regularizing chiral fermions. Physics Letters B, 105(2–3):219 – 223, 1981.
  • (8) S. Chandrasekharan and U.-J. Wiese. An introduction to chiral symmetry on the lattice. Progress in Particle and Nuclear Physics, 53(2):373 – 418, 2004.
  • (9) Hasan Almanasreh and Nils Svanstedt. G-convergence of dirac operators. Journal of Function Spaces and Applications, 2012.
  • (10) O. Kullie, D. Kolb, and A. Rutkowski. Two-spinor fully relativistic finite-element (fem) solution of the two-center coulomb problem. Chemical Physics Letters, 383(3–4):215 – 221, 2004.
  • (11) O. Kullie and D. Kolb. High accuracy dirac-finite-element () calculations for and . The European Physical Journal D - Atomic, Molecular, Optical and Plasma Physics, 17:167–173, 2001.
  • (12) J.P. Desclaux, J. Dolbeault, M.J. Esteban, P. Indelicato, and E. Séré. Computational approaches of relativistic models in quantum chemistry. In C. Le Bris, editor, Special Volume, Computational Chemistry, volume 10 of Handbook of Numerical Analysis, pages 453 – 483. Elsevier, 2003.
  • (13) N. Grün C. Müller and W. Scheid. Finite element formulation of the dirac equation and the problem of fermion doubling. Physics Letters A, 242(4–5):245 – 250, 1998.
  • (14) P. P. F. Wessels, W. J. Caspers, and F. W. Wiegel. Discretizing the one-dimensional dirac equation. EPL (Europhysics Letters), 46(2):123, 1999.

Appendix A Solution of the Initial Value Problem

This discussion follows closely to the derivation presented in the appendices of [1]. In order to compare our results with known solutions of the Dirac equation, we will first consider the case of the massless Dirac equation. Since we are interested in the behavior of particles whose mass is very close to zero, this should give us some indication of the fitness of our approach for real-world problems.

Removing the mass term from equation 1 and multiplying both sides by the matrix gives us the following equation.

We then make the following substitutions

and rearrange the equation as follows.


We may further simplify this equation into a first order ODE by taking the Fourier transform with respect to which expresses the massless Dirac operator in momentum space.


Integrating directly from to gives us the general solution to the massless initial value problem in momentum space.

Where we may then apply Euler’s Identity in order to remove the matrix from the exponential


Where is the identitity matrix. If we take our initial function to be a Gaussian wave of the form which may be expressed in momentum space as . Inverting the Fourier transform from equation 13 with the given initial value results in the general initial value solution:


where the values , , , and are define as

Equation 14 may then be used to calculate the analytic solution to any combination of massless Gaussian wave packets with the packet width given by and the momentum set by .

Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description