Nothing Special   »   [go: up one dir, main page]
Skip to main content
Springer Nature Link
Log in

Numerical Methods and Comparison for the Dirac Equation in the Nonrelativistic Limit Regime

  • Published:
Journal of Scientific Computing Aims and scope Submit manuscript

Abstract

We analyze rigorously error estimates and compare numerically spatial/temporal resolution of various numerical methods for the discretization of the Dirac equation in the nonrelativistic limit regime, involving a small dimensionless parameter \(0<\varepsilon \ll 1\) which is inversely proportional to the speed of light. In this limit regime, the solution is highly oscillatory in time, i.e. there are propagating waves with wavelength \(O(\varepsilon ^2)\) and O(1) in time and space, respectively. We begin with several frequently used finite difference time domain (FDTD) methods and obtain rigorously their error estimates in the nonrelativistic limit regime by paying particular attention to how error bounds depend explicitly on mesh size h and time step \(\tau \) as well as the small parameter \(\varepsilon \). Based on the error bounds, in order to obtain ‘correct’ numerical solutions in the nonrelativistic limit regime, i.e. \(0<\varepsilon \ll 1\), the FDTD methods share the same \(\varepsilon \)-scalability on time step and mesh size as: \(\tau =O(\varepsilon ^3)\) and \(h=O(\sqrt{\varepsilon })\). Then we propose and analyze two numerical methods for the discretization of the Dirac equation by using the Fourier spectral discretization for spatial derivatives combined with the symmetric exponential wave integrator and time-splitting technique for temporal derivatives, respectively. Rigorous error bounds for the two numerical methods show that their \(\varepsilon \)-scalability is improved to \(\tau =O(\varepsilon ^2)\) and \(h=O(1)\) when \(0<\varepsilon \ll 1\). Extensive numerical results are reported to support our error estimates.

This is a preview of subscription content, log in via an institution to check access.

Access this article

Log in via an institution

Subscribe and save

Springer+ Basic
$34.99 /Month
  • Get 10 units per month
  • Download Article/Chapter or eBook
  • 1 Unit = 1 Article or 1 Chapter
  • Cancel anytime
Subscribe now

Buy Now

Price excludes VAT (USA)
Tax calculation will be finalised during checkout.

Instant access to the full article PDF.

Institutional subscriptions

Fig. 1
Fig. 2
Fig. 3

Similar content being viewed by others

References

  1. Abanin, D.A., Morozov, S.V., Ponomarenko, L.A., Gorbachev, R.V., Mayorov, A.S., Katsnelson, M.I., Watanabe, K., Taniguchi, T., Novoselov, K.S., Levito, L.S., Geim, A.K.: Giant nonlocality near the Dirac point in graphene. Science 332, 328–330 (2011)

    Article  Google Scholar 

  2. Abenda, S.: Solitary waves for the Maxwell–Dirac and Coulomb–Dirac models. Ann. Inst. Henri Poincaré Phys. Theor 68, 229–244 (1998)

    MathSciNet  MATH  Google Scholar 

  3. Ablowitz, M.J., Zhu, Y.: Nonlinear waves in shallow honeycomb lattices. SIAM J. Appl. Math. 72, 240–260 (2012)

    Article  MathSciNet  MATH  Google Scholar 

  4. Anderson, C.D.: The positive electron. Phys. Rev. 43, 491–498 (1933)

    Article  Google Scholar 

  5. Antoine, X., Lorin, E., Sater, J., Fillion-Gourdeau, F., Bandrauk, A.D.: Absorbing boundary conditions for relativistic quantum mechanics equations. J. Comput. Phys. 277, 268–304 (2014)

    Article  MathSciNet  MATH  Google Scholar 

  6. Archilla, B.G., Sanz-Serna, J.M., Skeel, R.D.: Long-time-step methods for oscillatory differential equations. SIAM J. Sci. Comput. 20, 930–963 (1998)

    Article  MathSciNet  MATH  Google Scholar 

  7. Arnold, A., Steinrück, H.: The ‘electromagnetic’ Wigner equation for an electron with spin. ZAMP 40, 793–815 (1989)

    Article  MathSciNet  MATH  Google Scholar 

  8. Bao, W., Cai, Y.: Mathematical theory and numerical methods for Bose–Einstein condensation. Kinet. Relat. Models 6, 1–135 (2013)

    Article  MathSciNet  MATH  Google Scholar 

  9. Bao, W., Cai, Y.: Optimal error estmiates of finite difference methods for the Gross–Pitaevskii equation with angular momentum rotation. Math. Comput. 82, 99–128 (2013)

    Article  MATH  Google Scholar 

  10. Bao, W., Cai, Y.: Uniform and optimal error estimates of an exponential wave integrator sine pseudospectral method for the nonlinear Schrödinger equation with wave operator. SIAM J. Numer. Anal. 52, 1103–1127 (2014)

    Article  MathSciNet  MATH  Google Scholar 

  11. Bao, W., Cai, Y., Zhao, X.: A uniformly accurate multiscale time integrator pseudospectral method for the Klein–Gordon equation in the nonrelativistic limit regime. SIAM J. Numer. Anal. 52, 2488–2511 (2014)

    Article  MathSciNet  MATH  Google Scholar 

  12. Bao, W., Dong, X.: Analysis and comparison of numerical methods for the Klein–Gordon equation in the nonrelativistic limit regime. Numer. Math. 120, 189–229 (2012)

    Article  MathSciNet  MATH  Google Scholar 

  13. Bao, W., Dong, X., Zhao, X.: An exponential wave integrator pseudospectral method for the Klein–Gordon–Zakharov system. SIAM J. Sci. Comput. 35, A2903–A2927 (2013)

    Article  MathSciNet  MATH  Google Scholar 

  14. Bao, W., Dong, X., Zhao, X.: Uniformly correct multiscale time integrators for highly oscillatory second order differention equations. J. Math. Study 47, 111–150 (2014)

    MathSciNet  MATH  Google Scholar 

  15. Bao, W., Shi, J., Markowich, P.A.: On time-splitting spectral approximation for the Schrödinger equation in the semiclassical regime. J. Comput. Phys. 175, 487–524 (2002)

    Article  MathSciNet  MATH  Google Scholar 

  16. Bao, W., Shi, J., Markowich, P.A.: Numerical study of time-splitting spectral discretizations of nonlinear Schrödinger equations in the semi-classical regimes. SIAM J. Sci. Comput. 25, 27–64 (2003)

    Article  MathSciNet  MATH  Google Scholar 

  17. Bao, W., Li, X.: An efficient and stable numerical method for the Maxwell–Dirac system. J. Comput. Phys. 199, 663–687 (2004)

    Article  MathSciNet  MATH  Google Scholar 

  18. Bechouche, P., Mauser, N., Poupaud, F.: (Semi)-nonrelativistic limits of the Dirac eqaution with external time-dependent electromagnetic field. Commun. Math. Phys. 197, 405–425 (1998)

    Article  MATH  Google Scholar 

  19. Bechouche, P., Mauser, N., Selberg, S.: On the asymptotic analysis of the Dirac–Maxwell system in the nonrelativistic limit. J. Hyper. Differ. Equat. 2, 129–182 (2005)

    Article  MathSciNet  MATH  Google Scholar 

  20. Bolte, J., Keppeler, S.: A semiclassical approach to the Dirac equation. Ann. Phys. 274, 125–162 (1999)

    Article  MathSciNet  MATH  Google Scholar 

  21. Booth, H.S., Legg, G., Jarvis, P.D.: Algebraic solution for the vector potential in the Dirac equation. J. Phys. A: Math. Gen. 34, 5667–5677 (2001)

    Article  MathSciNet  MATH  Google Scholar 

  22. Bournaveas, N.: Local existence for the Maxwell–Dirac equations in three space dimensions. Commun. Part. Differ. Equ. 21, 693–720 (1996)

    Article  MathSciNet  MATH  Google Scholar 

  23. Brinkman, D., Heitzinger, C., Markowich, P.A.: A convergent 2D finite-difference scheme for the Dirac–Poisson system and the simulation of graphene. J. Comput. Phys. 257, 318–332 (2014)

    Article  MathSciNet  MATH  Google Scholar 

  24. Cirincione, R.J., Chernoff, P.R.: Dirac and Klein Gordon equations: convergence of solutions in the nonrelativistic limit. Commun. Math. Phys. 79, 33–46 (1981)

    Article  MathSciNet  MATH  Google Scholar 

  25. Das, A.: General solutions of Maxwell–Dirac equations in 1 + 1 dimensional space-time and spatial confined solution. J. Math. Phys. 34, 3986–3999 (1993)

    Article  MathSciNet  MATH  Google Scholar 

  26. Das, A., Kay, D.: A class of exact plane wave solutions of the Maxwell–Dirac equations. J. Math. Phys. 30, 2280–2284 (1989)

    Article  MathSciNet  MATH  Google Scholar 

  27. Deuflhard, P.: A study of extrapolation methods based on multistep schemes without parasitic solutions. ZAMP 30, 177–189 (1979)

    Article  MathSciNet  MATH  Google Scholar 

  28. Dirac, P.A.M.: The quantum theory of the electron. Proc. R. Soc. Lond. A 117, 610–624 (1928)

    Article  MATH  Google Scholar 

  29. Dirac, P.A.M.: A theory of electrons and protons. Proc. R. Soc. Lond. A 126, 360–365 (1930)

    Article  MATH  Google Scholar 

  30. Dirac, P.A.M.: Principles of Quantum Mechanics. Oxford University Press, London (1958)

    MATH  Google Scholar 

  31. Dolbeault, J., Esteban, M.J., Séré, E.: On the eigenvalues of operators with gaps: applications to Dirac operator. J. Funct. Anal. 174, 208–226 (2000)

    Article  MathSciNet  MATH  Google Scholar 

  32. Esteban, M., Séré, E.: Existence and multiplicity of solutions for linear and nonlinear Dirac problems. Partial Differ. Equ. Appl. 12, 107–112 (1997)

    MathSciNet  MATH  Google Scholar 

  33. Esteban, M., Séré, E.: An overview on linear and nonlinear Dirac equations. Discrete Contin. Dyn. Syst. 8, 381–397 (2002)

    Article  MathSciNet  MATH  Google Scholar 

  34. Faou, E., Schratz, K.: Asymptotic preserving schemes for the Klein–Gordon equation in the non-relativistic limit regime. Numer. Math. 126, 441–469 (2014)

    Article  MathSciNet  MATH  Google Scholar 

  35. Fefferman, C.L., Weistein, M.I.: Honeycomb lattice potentials and Dirac points. J. Am. Math. Soc. 25, 1169–1220 (2012)

    Article  MathSciNet  MATH  Google Scholar 

  36. Fefferman, C.L., Weistein, M.I.: Wave packets in honeycomb structures and two-dimensional Dirac equations. Commun. Math. Phys. 326, 251–286 (2014)

    Article  MathSciNet  MATH  Google Scholar 

  37. Ferreira, A., Gomes, J.V., Nilsson, J., Mucciolo, E.R., Peres, N.M.R., Catro Neto, A.H.: Unified description of the dc-conductivity of monolayer and bilayer graphene at finite densities based on resonant scatterers. Phys. Rev. B 83, 165402 (2011)

    Article  Google Scholar 

  38. Fillion-Gourdeau, F., Lorin, E., Bandrauk, A.D.: Resonantly enhanced pair production in a simple diatomic model. Phys. Rev. Lett. 110, 013002 (2013)

    Article  Google Scholar 

  39. Fillion-Gourdeau, F., Lorin, E., Bandrauk, A.D.: A split-step numerical method for the time-dependent Dirac equation in 3-D axisymmetric geometry. J. Comput. Phys. 272, 559–587 (2014)

    Article  MathSciNet  MATH  Google Scholar 

  40. Foldy, L.L., Wouthuysen, S.A.: On the Dirac theory of spin \(1/2\) particles and its nonrelavistic limit. Phys. Rev. 78, 29–36 (1950)

    Article  MATH  Google Scholar 

  41. Fushchich, W.I., Shtelen, W.M.: On some exact solutions of the nonlinear Dirac equation. J. Phys. A: Math. Gen. 16, 271–277 (1983)

    Article  MathSciNet  MATH  Google Scholar 

  42. Gautschi, W.: Numerical integration of ordinary differential equations based on trigonometric polynomials. Numer. Math. 3, 381–397 (1961)

    Article  MathSciNet  MATH  Google Scholar 

  43. Gérad, P., Markowich, P.A., Mauser, N.J., Poupaud, F.: Homogenization limits and Wigner transforms. Commun. Pure Appl. Math. 50, 321–377 (1997)

    MathSciNet  MATH  Google Scholar 

  44. Gesztesy, F., Grosse, H., Thaller, B.: A rigorious approach to relativistic corrections of bound state energies for spin-\(1/2\) particles. Ann. Inst. Henri Poincaré Phys. Theor 40, 159–174 (1984)

    MATH  Google Scholar 

  45. Gosse, L.: A well-balanced and asymptotic-preserving scheme for the one-dimensional linear Dirac equation. BIT Numer. Math. 55, 433–458 (2015)

    Article  MathSciNet  MATH  Google Scholar 

  46. Grigore, D.R., Nenciu, G., Purice, R.: On the nonrelativistic limits of the Dirac Hamiltonian. Ann. Inst. Henri Poincaré Phys. Theor 51, 231–263 (1989)

    MathSciNet  MATH  Google Scholar 

  47. Gross, L.: The Cauchy problem for the coupled Maxwell and Dirac equations. Commun. Pure Appl. Math. 19, 1–15 (1966)

    Article  MathSciNet  MATH  Google Scholar 

  48. Hammer, R., Pötz, W., Arnold, A.: Single-cone real-space finite difference scheme for the time-dependent Dirac equation. J. Comput. Phys. 265, 50–70 (2014)

    Article  MathSciNet  MATH  Google Scholar 

  49. Hammer, R., Pötz, W., Arnold, A.: A dispersion and norm preserving finite difference scheme with transparent boundary conditions for the Dirac equation in (1 + 1)D. J. Comput. Phys. 256, 728–747 (2014)

    Article  MathSciNet  MATH  Google Scholar 

  50. Hairer, E., Lubich, C., Wanner, G.: Geometric Numerical Integration. Springer, Berlin (2002)

    Book  MATH  Google Scholar 

  51. Hochbruck, M., Lubich, C.: A Gautschi-type method for oscillatory second-order differential equations. Numer. Math. 83, 402–426 (1999)

    Article  MathSciNet  MATH  Google Scholar 

  52. Hochbruck, M., Ostermann, A.: Exponential integrators. Acta Numer. 19, 209–286 (2000)

    Article  MathSciNet  MATH  Google Scholar 

  53. Huang, Z., Jin, S., Markowich, P.A., Sparber, C., Zheng, C.: A time-splitting spectral scheme for the Maxwell–Dirac system. J. Comput. Phys. 208, 761–789 (2005)

    Article  MathSciNet  MATH  Google Scholar 

  54. Hunziker, W.: On the nonrelativistic limit of the Dirac theory. Commun. Math. Phys. 40, 215–222 (1975)

    Article  MathSciNet  Google Scholar 

  55. Iserles, A.: A First Course in the Numerical Analysis of Differential Equations. Cambridge University Press, Cambridge (2008)

    Book  MATH  Google Scholar 

  56. Iserles, A., Norsett, S.P.: From high oscillation to rapid approximation I: modified Fourier expansions. IMA J. Numer. Anal. 28, 862–887 (2008)

    Article  MathSciNet  MATH  Google Scholar 

  57. Kammerer, C.F.: Semi-classical analysis of a Dirac equaiton without adiabatic decoupling. Monatsh. Math. 142, 281–313 (2004)

    Article  MathSciNet  MATH  Google Scholar 

  58. Lubich, C.: On splitting methods for Schrödinger–Poisson and cubic nonlinear Schrödinger equations. Math. Comput. 77, 2141–2153 (2008)

    Article  MATH  Google Scholar 

  59. Masmoudi, N., Mauser, N.J.: The selfconsistent Pauli equaiton. Monatsh. Math. 132, 19–24 (2001)

    Article  MathSciNet  MATH  Google Scholar 

  60. Mauser, N.J.: Rigorous derivation of the Pauli equation with time-dependent electromagnetic field. VLSI Design 9, 415–426 (1999)

    Article  Google Scholar 

  61. Najman, B.: The nonrelativistic limit of the nonlinear Dirac equation. Ann. Inst. Henri Poincaré 9, 3–12 (1992)

    Article  MathSciNet  MATH  Google Scholar 

  62. Neto, A.H.C., Guinea, F., Peres, N.M.R., Novoselov, K.S., Geim, A.K.: The electronic properties of the graphene. Rev. Mod. Phys. 81, 109–162 (2009)

    Article  Google Scholar 

  63. Novoselov, K.S., Geim, A.K., Morozov, S.V., Jiang, D., Katsnelson, M.I., Grigorieva, I.V., Dubonos, S.V., Firsov, A.A.: Two-dimensional gas of massless Dirac fermions in graphene. Nature 438, 197–200 (2005)

    Article  Google Scholar 

  64. Novoselov, K.S., Geim, A.K., Morozov, S.V., Jiang, D., Zhang, Y., Dubonos, S.V., Grigorieva, I.V., Firsov, A.A.: Electric filed effect in atomically thin carbon films. Science 306, 666–669 (2004)

    Article  Google Scholar 

  65. Novoselov, K.S., Jiang, Z., Zhang, Y., Morozov, S.V., Stormer, H.L., Zeitler, U., Maan, J.C., Boebinger, G.S., Kim, P., Geim, A.K.: Room-temperature quantum Hall effect in graphene. Science 315, 1379 (2007)

    Article  Google Scholar 

  66. Nraun, J.W., Su, Q., Grobe, R.: Numerical approach to solve the time-dependent Dirac equation. Phys. Rev. A 59, 604–612 (1999)

    Article  Google Scholar 

  67. Schedin, F., Geim, A., Morozov, S., Hill, E., Blake, P., Katsnelson, M., Novoselov, K.: Detection of individual gas molecules absorbed on graphene. Nat. Mater. 6, 652–655 (2007)

    Article  Google Scholar 

  68. Schoene, A.Y.: On the nonrelativistic limits of the Klein–Gordon and Dirac equations. J. Math. Anal. Appl. 71, 36–74 (1979)

    Article  MathSciNet  MATH  Google Scholar 

  69. Shebalin, J.V.: Numerical solution of the coupled Dirac and Maxwell equations. Phys. Lett. A 226, 1–6 (1997)

    Article  MathSciNet  MATH  Google Scholar 

  70. Shen, J., Tang, T.: Spectral and High-Order Methods with Applications. Science Press, Beijing (2006)

    MATH  Google Scholar 

  71. Smith, G.D.: Numerical Solution of Partial Differential Equations: Finite Difference Methods. Clarendon Press, Oxford (1985)

    MATH  Google Scholar 

  72. Spohn, H.: Semiclassical limit of the Dirac equaiton and spin precession. Ann. Phys. 282, 420–431 (2000)

    Article  MATH  Google Scholar 

  73. Strang, G.: On the construction and comparision of difference schemes. SIAM J. Numer. Anal. 5, 505–517 (1968)

    Article  Google Scholar 

  74. Thaller, B.: The Dirac Equation. Springer, New York (1992)

    Book  MATH  Google Scholar 

  75. Veselic, K.: Perturbation of pseudoresolvents and analyticity in \(1/c\) ofrelativistic quantum mechanics. Commun. Math. Phys. 22, 27–43 (1971)

    Article  MATH  Google Scholar 

  76. White, G.B.: Splitting of the Dirac operator in the nonrelativistic limit. Ann. Inst. Henri Poincaré 53, 109–121 (1990)

    MathSciNet  MATH  Google Scholar 

  77. Wu, H., Huang, Z., Jin, S., Yin, D.: Gaussian beam methods for the Dirac equation in the semi-classical regime. Commun. Math. Sci. 10, 1301–1315 (2012)

    Article  MathSciNet  MATH  Google Scholar 

Download references

Acknowledgements

Part of this work was done when the authors were visiting the Institute for Mathematical Sciences at the National University of Singapore in 2015.

Author information

Authors and Affiliations

  1. Department of Mathematics, National University of Singapore, Singapore, 119076, Singapore

    Weizhu Bao & Xiaowei Jia

  2. Beijing Computational Science Research Center, No. 10 West Dongbeiwang Road, Haidian District, Beijing, 100193, People’s Republic of China

    Yongyong Cai

  3. Department of Mathematics, Purdue University, West Lafayette, IN, 47907, USA

    Yongyong Cai

  4. Institut Elie Cartan de Lorraine, Inria Nancy-Grand Est, Université de Lorraine, 54506, Vandoeuvre-lès-nancy Cedex, France

    Qinglin Tang

Authors
  1. Weizhu Bao
    View author publications

    You can also search for this author in PubMed Google Scholar

  2. Yongyong Cai
    View author publications

    You can also search for this author in PubMed Google Scholar

  3. Xiaowei Jia
    View author publications

    You can also search for this author in PubMed Google Scholar

  4. Qinglin Tang
    View author publications

    You can also search for this author in PubMed Google Scholar

Corresponding author

Correspondence to Yongyong Cai.

Additional information

This work was partially supported by the Ministry of Education of Singapore Grant R-146-000-223-112 (W. Bao and X. Jia), the Natural Science Foundation of China Grant U1530401 and the NSF Grants DMS-1217066 and DMS-1419053 (Y. Cai), and the ANR project BECASIM ANR-12-MONU-0007-02 (Q. Tang).

Appendices

Appendix 1

Proof of Theorem 2.1 for the LFFD method

Define the local truncation error \(\tilde{\xi }^{n}=(\tilde{\xi }_0^n,\tilde{\xi }_1^n,\ldots ,\tilde{\xi }_M^n)^T\in X_{M}\) of the LFFD (2.6) with (2.10) and (2.11) as follows, for \(0\le j\le M-1\) and \(n\ge 1\),

$$\begin{aligned} \tilde{\xi }_{j}^{n}:&=\left[ i\delta _{t}\varPhi + \frac{i}{\varepsilon }\sigma _1\delta _{x}\varPhi -\frac{1}{\varepsilon ^2}\sigma _3\varPhi + \left( A_{1,j}^n\sigma _1-V_j^nI_2\right) \varPhi \right] _{t=t_n,x=x_j}, \end{aligned}$$
(6.1)
$$\begin{aligned} \tilde{\xi }_{j}^{0}:&=i\delta _{t}^+\varPhi (0,x_{j})+\frac{i}{\varepsilon }\sigma _1\delta _x\varPhi _0(x_j)-\left( \frac{1}{\varepsilon ^2}\sigma _3 +V_j^0I_2-A_{1,j}^0\sigma _1\right) \varPhi _0(x_j). \end{aligned}$$
(6.2)

Applying the Taylor expansion in (6.1) and (6.2) we obtain for \( j=0,1,\ldots , M-1\) and \(n\ge 1\),

$$\begin{aligned} \tilde{\xi }_{j}^0=i\tau \partial _{tt}\varPhi (\tau ^{\prime },x_j)+ \frac{i}{\varepsilon }h^2\partial _{xxx}\varPhi _0(x^\prime ),\quad \tilde{\xi }_{j}^n=i\tau ^2\partial _{ttt}\varPhi (t^\prime ,x_j) +\frac{i}{\varepsilon }h^2\partial _{xxx}\varPhi (t_n,x^{\prime \prime }), \end{aligned}$$

where \(t^\prime \in (0,\tau )\), \(t^{\prime \prime }\in (t_{n-1},t_{n+1})\), \(x^\prime ,x^{\prime \prime }\in (x_{j-1},x_{j+1})\). Noticing (2.1) and the assumptions (A) and (B), we have

$$\begin{aligned} |\tilde{\xi }_{j}^0|\lesssim \frac{\tau }{\varepsilon ^4}+\frac{h^2}{\varepsilon }, \quad |\tilde{\xi }_{j}^n|\lesssim \frac{\tau ^2}{\varepsilon ^6}+\frac{h^2}{\varepsilon }, \quad j=0,1,\ldots , M-1, \quad n\ge 1, \end{aligned}$$
(6.3)

which immediately implies

$$\begin{aligned} \Vert \tilde{\xi }^n\Vert _{l^\infty }=\max _{0\le j\le M-1}|\tilde{\xi }_j^n| \lesssim \frac{\tau ^2}{\varepsilon ^6}+\frac{h^2}{\varepsilon }, \ \Vert \tilde{\xi }^n\Vert _{l^2}\lesssim \Vert \tilde{\xi }^n\Vert _{l^\infty } \lesssim \frac{\tau ^2}{\varepsilon ^6}+\frac{h^2}{\varepsilon }, \ n\ge 1. \end{aligned}$$
(6.4)

Subtracting (2.6) from (6.1), noticing (2.32), we get for \(0\le j\le M-1\) and \( n\ge 1\),

$$\begin{aligned} i\delta _t\mathbf{e} _{j}^{n}=-\frac{i}{\varepsilon }\sigma _1\delta _x\mathbf{e} _{j}^{n}+\frac{1}{\varepsilon ^2}\sigma _3 \mathbf{e} _{j}^{n}+\left( V_j^{n}I_2-A_{1,j}^{n}\sigma _1 \right) \mathbf{e} _{j}^{n}+\tilde{\xi }_j^n, \end{aligned}$$
(6.5)

where the boundary and initial conditions are given as

$$\begin{aligned} \mathbf{e} _0^{n}=\mathbf{e} _{M}^{n}, \quad \mathbf{e} _{-1}^{n}=\mathbf{e} _{M-1}^n, \quad n\ge 0,\quad \mathbf{e} _j^0=\mathbf{0}, \quad j=0,1,\ldots ,M. \end{aligned}$$
(6.6)

For the first step, we have

$$\begin{aligned} \Vert \mathbf{e} ^1\Vert _{l^2}=\tau \Vert \tilde{\xi }^0\Vert _{l^2}\lesssim \frac{\tau ^2}{\varepsilon ^4}+\frac{\tau h^2}{\varepsilon }\lesssim \frac{h^2}{\varepsilon }+\frac{\tau ^2}{\varepsilon ^6}. \end{aligned}$$
(6.7)

Denote \({\mathcal {E}}^{n+1}\) for \(n=0, 1, \ldots \) as

$$\begin{aligned} {\mathcal {E}}^{n+1}=&\Vert \mathbf{e} ^{n+1}\Vert ^2_{l^2}+\Vert \mathbf{e} ^{n}\Vert _{l^2}^2+2\,\text {Re}\left( \tau h\sum \limits _{j=0}^{M-1}(\mathbf{e} _j^{n+1})^*\sigma _1\delta _x\mathbf{e} ^n_j\right) \nonumber \\&-2\,\text {Im}\left( \frac{\tau h}{\varepsilon ^2}\sum \limits _{j=0}^{M-1}(\mathbf{e} _j^{n+1})^*\sigma _3\mathbf{e} _j^n\right) ; \end{aligned}$$
(6.8)

and under the stability condition (2.33), e.g., \(\tau \le \frac{\varepsilon ^2\tau _1 h}{\varepsilon ^2h V_{\mathrm{max}}+\sqrt{h^2+\varepsilon ^2(1+\varepsilon h A_{1,\mathrm{max}})^2}}\) with \(\tau _1=\frac{1}{4}\), which implies \(\frac{\tau }{h}\le \frac{1}{4}\) and \(\frac{\tau }{\varepsilon ^2}\le \frac{1}{4}\), using Cauchy inequality, we can get that

$$\begin{aligned} \frac{1}{2}\left( \Vert \mathbf{e} ^{n+1}\Vert ^2_{l^2}+\Vert \mathbf{e} ^{n}\Vert _{l^2}^2\right) \le {\mathcal {E}}^{n+1}\le \frac{3}{2}\left( \Vert \mathbf{e} ^{n+1}\Vert ^2_{l^2}+\Vert \mathbf{e} ^{n}\Vert _{l^2}^2\right) , \quad n\ge 0. \end{aligned}$$
(6.9)

It follows from (6.7) that

$$\begin{aligned} {\mathcal {E}}^1\lesssim \left( \frac{h^2}{\varepsilon }+\frac{\tau ^2}{\varepsilon ^6}\right) ^2. \end{aligned}$$
(6.10)

Multiplying (6.5) from the left by \(2h\tau (\mathbf{e} _j^{n+1}+\mathbf{e} _j^{n-1})^*\), taking the imaginary part, then summing for \(j=0,1,\ldots ,M-1\), using Cauchy inequality, noting (6.4) and (6.9), we get for \(n\ge 1\),

$$\begin{aligned} {\mathcal {E}}^{n+1}-{\mathcal {E}}^{n}\lesssim&h\tau \sum _{j=0}^{M-1}\left( (A_{1,\mathrm max}+V_\mathrm{max})|\mathbf{e} _j^n|+|\tilde{\xi }_j^n| \right) (|\mathbf{e} _j^{n+1}|+|\mathbf{e} _j^{n-1}|)\\ \lesssim&\tau ({\mathcal {E}}^{n+1}+{\mathcal {E}}^{n})+\tau \left( \frac{h^2}{\varepsilon }+\frac{\tau ^2}{\varepsilon ^6}\right) ^2, \quad n\ge 0. \end{aligned}$$

Summing the above inequality for \(n=1, 2, \ldots , m-1\), we get

$$\begin{aligned} {\mathcal {E}}^m-{\mathcal {E}}^1\lesssim \tau \sum _{k=1}^m{\mathcal {E}}^k+m\tau \left( \frac{h^2}{\varepsilon }+\frac{\tau ^2}{\varepsilon ^6}\right) ^2, \quad 1\le m\le \frac{T}{\tau }. \end{aligned}$$
(6.11)

Taking \(\tau _0\) sufficiently small, using the discrete Gronwall’s inequality and noticing (6.10), we obtain from the above equation that

$$\begin{aligned} {\mathcal {E}}^m \lesssim \left( \frac{h^2}{\varepsilon }+\frac{\tau ^2}{\varepsilon ^6}\right) ^2, \quad 1\le m\le \frac{T}{\tau }, \end{aligned}$$
(6.12)

which immediately implies the error bound (2.34) in view of (6.9). \(\square \)

Appendix 2

Proof of Theorem 3.1 for the sEWI-FP method

Define the error function \(\mathbf{e} ^n(x)\) for \(n=0,1,\ldots \) as

$$\begin{aligned} \mathbf{e} ^n(x)= \begin{pmatrix} e_1^n(x)\\ e_2^n(x) \\ \end{pmatrix} :=P_M\varPhi (t_n,x)-\varPhi _M^n(x)=\sum _{l=-M/2}^{M/2-1}\widehat{\mathbf{e} }_l^n e^{i\mu _l(x-a)},\ a\le x\le b.\quad \end{aligned}$$
(6.13)

Using the triangular inequality and standard interpolation result, we get

$$\begin{aligned} \Vert \varPhi ( t_n,x)-\varPhi _M^n(x)\Vert _{L^2}\le & {} \Vert \varPhi ( t_n,x)-P_M\varPhi (t_n,x)\Vert _{L^2}+\Vert \mathbf{e} ^n(x)\Vert _{L^2}\nonumber \\\le & {} h^{m_0}+\Vert \mathbf{e} ^n(x)\Vert _{L^2}, \end{aligned}$$
(6.14)

where \(0\le n\le \frac{T}{\tau }\), and the above result means that we only need estimate \(\Vert \mathbf{e} ^n(x)\Vert _{L^2}\).

Define the local truncation error \(\xi ^n(x)=\sum _{l=-M/2}^{M/2-1}\widehat{\xi }_l^ne^{i\mu _l(x-a)}\in Y_M\) of the sEWI-FP (3.17) for \(n\ge 1\) as

$$\begin{aligned} \widehat{\xi }_l^n= \widehat{(\varPhi (t_{n+1}))}_l +2i\sin (\frac{\tau \varGamma _l}{\varepsilon ^2})\,\widehat{(\varPhi (t_n))}_l-\widehat{(\varPhi (t_{n-1})}_l +2i\frac{\varepsilon ^2}{\delta _l}\sin (\frac{\tau \delta _l}{\varepsilon ^2})\widehat{(G(t_n)\varPhi (t_n))}_l, \end{aligned}$$
(6.15)

and for \(n=0\) as

$$\begin{aligned} \widehat{\xi }_l^0=\widehat{(\varPhi (\tau ))}_l -e^{-i\tau \varGamma _l/\varepsilon ^2}\widehat{(\varPhi (0))}_l +\varepsilon ^2\varGamma _l^{-1}\left[ I_2-e^{-\frac{i\tau }{\varepsilon ^2}\varGamma _l}\right] \widehat{(G(0)\varPhi (0))}_l, \end{aligned}$$
(6.16)

where we write \(\varPhi (t)\) and G(t) in short for \(\varPhi (t,x)\) and G(t, x), respectively.

Firstly, we estimate the local truncation error \(\xi ^n(x)\). Multiplying both sides of the Dirac equation (2.1) by \(e^{i\mu _l(x-a)}\) and integrating over the interval (a, b), we easily recover the equations for \((\widehat{\varPhi (t)})_l\), which are exactly the same as (3.6) with \(\varPhi _M\) being replaced by \(\varPhi (t,x)\). Replacing \(\varPhi _M\) with \(\varPhi (t,x)\), we use the same notations \(\widehat{F}_l^n(s)\) as in (3.7) and the time derivatives of \(\widehat{F}_l^n(s)\) enjoy the same properties of time derivatives of \(\varPhi (t,x)\). Thus, the same representation (3.12) holds for \((\widehat{\varPhi (t_n)})_l\) with \(n\ge 1\). From the derivation of the EWI method, it is clear that the error \(\xi ^n(x)\) comes from the approximations for the integrals in (3.13) and (3.14), and we have

$$\begin{aligned} \widehat{\xi }^0_l=-i\int _0^{\tau }e^{\frac{i(s-\tau )}{\varepsilon ^2}\varGamma _l}(\widehat{F}^0_l(s) -\widehat{F}^0_l(0))ds=-i\int _0^{\tau }\int _0^s e^{\frac{i(s-\tau )}{\varepsilon ^2}\varGamma _l} \partial _{s_1}\widehat{F}^0_l(s_1)\,ds_1ds, \end{aligned}$$

and for \(n\ge 1\)

$$\begin{aligned} \widehat{\xi }^n_l=&-i\int _0^{\tau }\cos ((s-\tau )\delta _l/\varepsilon ^2) \int _0^s\int _{-s_1}^{s_1}\partial _{s_2s_2} \widehat{F}^{n}_l(s_2)\,ds_2ds_1ds \nonumber \\&+\int _0^\tau \sin \left( \frac{(s-\tau )\varGamma _l}{\varepsilon ^2}\right) \int _{-s}^s\partial _{s_1}\widehat{F}_l^n(s_1) \,ds_1\,ds. \end{aligned}$$
(6.17)

For \(n=0\), the above equalities imply \(|\widehat{\xi }^0_l|\lesssim \int _0^{\tau }\int _0^s|\partial _{s_1}\widehat{F}^0_l(s_1)|ds_1ds\) and by the Bessel inequality and assumptions (C) and (D), we find

$$\begin{aligned} \Vert \xi ^0(x)\Vert _{L^2}^2=&(b-a)\sum \limits _{l=-M/2}^{M/2-1}|\widehat{\xi }_l^0|^2 \lesssim (b-a)\tau ^2\int _0^\tau \int _0^s\sum \limits _{l=-M/2}^{M/2-1}|\partial _{s_1} \widehat{F}^0_l(s_1)|^2\,ds_1ds\\ \lesssim&\tau ^2\int _0^\tau \int _0^s\Vert \partial _{s_1}(G(s_1)\varPhi (s_1))\Vert _{L^2}^2\,ds_1ds \lesssim \frac{\tau ^4}{\varepsilon ^4}. \end{aligned}$$

Similarly, for \(n\ge 1\), we obtain

$$\begin{aligned}&\Vert \xi ^n(x)\Vert _{L^2}^2=(b-a)\sum \limits _{l=-M/2}^{M/2-1}|\widehat{\xi }_l^n|^2\\&\lesssim \tau ^3\int _0^{\tau }\int _0^s\int _{-s_1}^{s_1} \sum \limits _{l=-\frac{M}{2}}^{\frac{M}{2}-1}|\partial _{s_2s_2}\widehat{F}_l^n(s_2)|^2\,ds_2ds_1ds\\&\quad +\tau ^2\int _0^\tau \int _{-s}^{s}\frac{(\tau -s)^2}{\varepsilon ^4}\sum \limits _{l=-\frac{M}{2}}^{\frac{M}{2}-1} |\partial _{\theta _1}\widehat{F}_l^{n-1}(\theta _1)|^2\,d\theta _1\,d\theta \,ds\\&\lesssim \tau ^6\Vert \partial _{tt}(G(t)\varPhi (t))\Vert _{L^\infty ([0,T]; (L^2)^2)}^2+\frac{\tau ^6}{\varepsilon ^4}\Vert \partial _{t}(G(t)\varPhi (t))\Vert _{L^\infty ([0,T]; (L^2)^2)}^2 \lesssim \frac{\tau ^6}{\varepsilon ^8}, \end{aligned}$$

where we have used the assumptions (C) and (D). Hence, we derive that

$$\begin{aligned} \Vert \xi ^0(x)\Vert _{L^2}\lesssim \frac{\tau ^2}{\varepsilon ^2},\quad \Vert \xi ^n(x)\Vert _{L^2}\lesssim \frac{\tau ^3}{\varepsilon ^4},\quad n\ge 1. \end{aligned}$$
(6.18)

Now, we look at the error equations. For each fixed \(l=-M/2,\ldots ,M/2-1\), subtracting (3.17) from (6.15), we obtain the equation for the error vector function as

$$\begin{aligned} \widehat{\mathbf{e} }^0_l=\mathbf{0},\quad \widehat{\mathbf{e} }^1_l=\widehat{\xi }^0_l; \quad \widehat{\mathbf{e} }^{n+1}_l-\widehat{\mathbf{e} }^{n-1}_l= -2i\sin (\tau \varGamma _l/\varepsilon ^2)\widehat{\mathbf{e} }^n_l +\widehat{R}^n_l +\widehat{\xi }^n_l, \end{aligned}$$
(6.19)

where \(1\le n\le \frac{T}{\tau }-1,\) and \(R^n(x)=\sum \limits _{l=-M/2}^{M/2-1}\widehat{R}_l^ne^{i\mu _l(x-a)}\in Y_M\) for \(n\ge 1\) is given by

$$\begin{aligned} \widehat{R}^n_l=2i\varepsilon ^2\delta _l^{-1}\sin (\tau \delta _l/\varepsilon ^2)\left( \widehat{(G(t_n)\varPhi (t_n))}_l-\widehat{(G(t_n)\varPhi ^n_M)}_l\right) . \end{aligned}$$
(6.20)

Since \(|\varepsilon ^2\delta _l^{-1}\sin (\tau \delta _l/\varepsilon ^2)|\le \tau \), from (6.20) and the assumption (D), we get

$$\begin{aligned} \Vert R^n(x)\Vert _{L^2}^2=&(b-a)\sum \limits _{l=-M/2}^{M/2-1}|\widehat{R}_l^n|^2\nonumber \\ \lesssim&(b-a)\tau ^2\sum \limits _{l=-M/2}^{M/2-1}\left| \widehat{(G(t_n)\varPhi (t_n))}_l-(\widehat{G(t_n)\varPhi _M^n})_l\right| ^2 \nonumber \\ \lesssim&\tau ^2\Vert G(t_n)\varPhi ( t_n,x)-G(t_n)\varPhi _M^n(x)\Vert _{L^2}^2\lesssim \tau ^2\Vert \varPhi ( t_n,x)-\varPhi _M^n(x)\Vert _{L^2}^2\nonumber \\ \lesssim&\tau ^2 h^{2m_0} +\tau ^2\Vert \mathbf{e} ^{n}(x)\Vert _{L^2}^2. \end{aligned}$$
(6.21)

Multiplying both sides of (6.19) by \(\left( \widehat{\mathbf{e} }^n_l\right) ^*\) from left, taking the real parts, we obtain

$$\begin{aligned} \text {Re}\left( (\widehat{\mathbf{e} }^n_l)^*\widehat{\mathbf{e} }^{n+1}_l\right) -\text {Re}\left( (\widehat{\mathbf{e} }^n_l)^*\widehat{\mathbf{e} }^{n-1}_l\right) =\text {Re}\left( (\widehat{\mathbf{e} }^n_l)^*(\widehat{R}^n_l +\widehat{\xi }^n_l)\right) , \end{aligned}$$

which implies

$$\begin{aligned} |\widehat{\mathbf{e} }^{n+1}_l|^2+|\widehat{\mathbf{e} }^{n}_l|^2 -|\widehat{\mathbf{e} }^{n+1}_l-\widehat{\mathbf{e} }^{n}_l|^2= & {} |\widehat{\mathbf{e} }^{n}_l|^2+|\widehat{\mathbf{e} }^{n-1}_l|^2 -|\widehat{\mathbf{e} }^{n}_l-\widehat{\mathbf{e} }^{n-1}_l|^2\nonumber \\&+2\,\text {Re} \left( (\widehat{\mathbf{e} }^n_l)^*(\widehat{R}^n_l +\widehat{\xi }^n_l)\right) . \end{aligned}$$
(6.22)

Multiplying both sides of (6.19) by \(\left( \widehat{\mathbf{e} }_l^{n+1}-2\widehat{\mathbf{e} }_l^{n}+\widehat{\mathbf{e} }_l^{n-1}\right) ^*\) from left, taking the real parts, we have

$$\begin{aligned}&\left| \widehat{\mathbf{e} }^{n+1}_l- \widehat{\mathbf{e} }^n_l\right| ^2-\left| \widehat{\mathbf{e} }^{n}_l- \widehat{\mathbf{e} }^{n-1}_l\right| ^2 \nonumber \\ =&2\,\text {Im}\left( (\widehat{\mathbf{e} }^{n+1}_l)^*\sin (\tau \varGamma _l/\varepsilon ^2)\widehat{\mathbf{e} }^{n}_l\right) - 2\,\text {Im}\left( (\widehat{\mathbf{e} }^{n}_l)^*\sin (\tau \varGamma _l/\varepsilon ^2)\widehat{\mathbf{e} }^{n-1}_l\right) \nonumber \\&+\text {Re}\left( (\widehat{\mathbf{e} }^{n+1}_l- 2\widehat{\mathbf{e} }^{n}_l+\widehat{\mathbf{e} }^{n-1}_l)^*(\widehat{R}^n_l+\widehat{\xi }_l^n)\right) . \end{aligned}$$
(6.23)

Summing (6.22) and (6.23), then applying Cauchy inequality and triangle inequality, we get

$$\begin{aligned}&|\widehat{\mathbf{e} }^{n+1}_l|^2+|\widehat{\mathbf{e} }^{n}_l|^2-2\, \text {Im}\left( (\widehat{\mathbf{e} }^{n+1}_l)^*\sin (\tau \varGamma _l/\varepsilon ^2) \widehat{\mathbf{e} }^{n}_l\right) \nonumber \\ \le&|\widehat{\mathbf{e} }^{n}_l|^2+|\widehat{\mathbf{e} }^{n-1}_l|^2 -2\text {Im}\left( (\widehat{\mathbf{e} }^{n}_l)^*\sin (\tau \varGamma _l/\varepsilon ^2) \widehat{\mathbf{e} }^{n-1}_l\right) \nonumber \\&+\tau (|\widehat{\mathbf{e} }^{n+1}_l|^2+|\widehat{\mathbf{e} }^{n-1}_l|^2) +\frac{1}{\tau }(|\widehat{R}^n_l|^2+|\widehat{\xi }^n_l|^2). \end{aligned}$$
(6.24)

Denote

$$\begin{aligned} \mathcal{E}^n=\left\| \mathbf{e} ^{n+1}(x)\right\| _{L^2}^2+\left\| \mathbf{e} ^{n}(x)\right\| _{L^2}^2-2(b-a)\sum \limits _{l=-M/2}^{M/2-1} \text {Im}\left( (\widehat{\mathbf{e} }^{n+1}_l)^*\sin (\tau \varGamma _l/\varepsilon ^2)\widehat{\mathbf{e} }^{n}_l\right) , \end{aligned}$$
(6.25)

and it follows from the stability constraint (3.25) that the matrix \(l^2\) norm satisfies \(\Vert \sin (\frac{\tau \varGamma _l}{\varepsilon ^2})\Vert _{l^2}\le \sin (\tau \delta _l/\varepsilon ^2)\le \sin (\pi /3)=\sqrt{3}/2\), which yield the following conclusion

$$\begin{aligned} \mathcal{E}^n\ge&\sum _{k=n}^{n+1}\left\| \mathbf{e} ^{k}(x)\right\| _{L^2}^2- \frac{\sqrt{3}}{2}(b-a)\sum \limits _{l=-M/2}^{M/2-1} \left( \left| \widehat{\mathbf{e} }^{n+1}_l\right| ^2+\left| \widehat{\mathbf{e} }^{n}_l\right| ^2\right) \nonumber \\ =&\frac{2-\sqrt{3}}{2}(\left\| \mathbf{e} ^{n+1}(x)\right\| _{L^2}^2+ \left\| \mathbf{e} ^{n}(x)\right\| _{L^2}^2). \end{aligned}$$
(6.26)

Multiplying (6.24) by \(b-a\) and summing together for \(l=-M/2, \ldots , M/2-1\), in view of the Bessel inequality, we obtain

$$\begin{aligned} \mathcal{E}^{n}-\mathcal{E}^{n-1}\lesssim&\tau (\left\| \mathbf{e} ^{n+1}(x)\right\| _{L^2}^2+\left\| \mathbf{e} ^{n}(x)\right\| _{L^2}^2+\left\| \mathbf{e} ^{n-1}(x)\right\| _{L^2}^2)\nonumber \\&+\frac{1}{\tau }\Vert R^n(x)\Vert ^2_{L^2}+\frac{1}{\tau }\Vert \xi ^n(x)\Vert ^2_{L^2},\quad n\ge 1. \end{aligned}$$
(6.27)

Summing (6.27) for \(n=1,\ldots ,m-1\), using (6.21) and (6.18), we derive

$$\begin{aligned} \mathcal{E}^{m-1}-\mathcal{E}^0\lesssim \tau \sum \limits _{k=1}^{m}\left\| \mathbf{e} ^{k}(x)\right\| _{L^2}^2+\frac{m\tau ^5}{\varepsilon ^8}+m\tau h^{2m_0},\quad 1\le m\le \frac{T}{\tau }. \end{aligned}$$
(6.28)

Since \(\mathbf{e} ^0(x)=\mathbf{0}\) and \(\mathcal{E}^{m-1}\) is bounded from below (6.26), we have for \(1\le m\le \frac{T}{\tau }\),

$$\begin{aligned} \frac{2-\sqrt{3}}{2}(\Vert \mathbf{e} ^{m}(x)\Vert _{L^2}^2+\Vert \mathbf{e} ^{m-1}(x)\Vert _{L^2}^2)-\Vert \mathbf{e} ^{1}(x)\Vert _{L^2}^2\lesssim \tau \sum \limits _{k=1}^{m}\Vert \mathbf{e} ^{k}(x)\Vert _{L^2}^2+\frac{m\tau ^5}{\varepsilon ^8}+m\tau h^{2m_0}. \end{aligned}$$
(6.29)

Noticing \(\Vert \mathbf{e} ^1(x)\Vert _{L^2}\lesssim \frac{\tau ^2}{\varepsilon ^2}\lesssim \frac{\tau ^2}{\varepsilon ^4}\), the discrete Gronwall’s inequality will imply that for sufficiently small \(\tau \),

$$\begin{aligned} \left\| \mathbf{e} ^{m}(x)\right\| _{L^2}^2\lesssim h^{2m_0}+\frac{\tau ^4}{\varepsilon ^8},\quad 1\le m\le \frac{T}{\tau }. \end{aligned}$$
(6.30)

Combining (6.14) and (6.30), we draw the conclusion (3.26). \(\square \)

Appendix 3

Extensions of the sEWI-FS (3.16–3.17) and TSFP (4.4) in 2D and 3D

The sEWI-FS (3.16–3.17), sEWI-FP (3.18–3.19) and TSFP (4.4) can be easily extended to 2D and 3D with tensor grids by modifying the matrices \(\varGamma _l\) in (3.8) and G(t, x) in (4.5) in the TSFP case. For the reader’s convenience, we present the modifications of \(\varGamma _l\) in (3.8) and G(t, x) in (4.5) in 2D and 3D as follows.

For the Dirac equation (1.21) in 2D, i.e. we take \(d=2\) in (1.21). The problem is truncated on \(\varOmega =(a_1, b_1)\times (a_2, b_2)\) with mesh sizes \(h_1=(b_1-a_1)/M_1\) and \(h_2=(b_2-a_2)/M_2\) (\(M_1,M_2\) two even positive integers) in the x- and y-direction, respectively. The wave function \(\varPhi \) is a two-component vector, and the matrix \(\varGamma _l\) in (3.8) will be replaced by

$$\begin{aligned} \varGamma _{jk}=\begin{pmatrix} 1 &{} \varepsilon \mu _j^{(1)}-i\varepsilon \mu _k^{(2)} \\ \varepsilon \mu _j^{(1)}+i\varepsilon \mu _k^{(2)} &{} -1\\ \end{pmatrix},\ \mu _j^{(1)}=\frac{2j\pi }{b_1-a_1},\ \mu _k^{(2)}=\frac{2k\pi }{b_2-a_2}, \end{aligned}$$
(6.31)

where \(-\frac{M_1}{2}\le j\le \frac{M_1}{2}-1\), \(-\frac{M_2}{2}\le k\le \frac{M_2}{2}-1\), and the Schur decomposition \(\varGamma _{jk}=Q_{jk}D_{jk}Q_{jk}^*\) is given as

$$\begin{aligned} Q_{jk}=\begin{pmatrix} \frac{1+\delta _{jk}}{\sqrt{2\delta _{jk}(1+\delta _{jk})}} &{}\frac{-\varepsilon \mu _{j}^{(1)}+i\varepsilon \mu _k^{(2)}}{\sqrt{2\delta _{jk}(1+\delta _{jk})}}\\ \frac{\varepsilon \mu _{j}^{(1)}+i\varepsilon \mu _k^{(2)}}{\sqrt{2\delta _{jk}(1+\delta _{jk})}} &{}\frac{1+\delta _{jk}}{\sqrt{2\delta _{jk}(1+\delta _{jk})}} \end{pmatrix}, \quad D_{jk}=\begin{pmatrix} \delta _{jk} &{}0\\ 0 &{}-\delta _{jk}\\ \end{pmatrix}, \end{aligned}$$
(6.32)

where

$$\begin{aligned} \delta _{jk}=\sqrt{1+\varepsilon ^2(\mu _{j}^{(1)})^2+\varepsilon ^2(\mu _k^{(2)})^2}. \end{aligned}$$
(6.33)

The matrix \(\int _{t_n}^{t_{n+1}}G(t,\mathbf{x} )dt\) in (4.5) becomes \(\int _{t_n}^{t_{n+1}}G(t,\mathbf{x} )dt\) and the Schur decomposition \(\int _{t_n}^{t_{n+1}}G(t,\mathbf{x} )dt=P_\mathbf{x} \varLambda _\mathbf{x} P_{\mathbf{x} }^*\) with \(V_\mathbf{x} ^{(1)}=\int _{t_n}^{t_{n+1}}V(t,\mathbf{x} )dt\), \(A_{l,\mathbf{x} }^{(1)}=\int _{t_n}^{t_{n+1}}A_l(t,\mathbf{x} )dt\) for \(l=1,2\), \(\lambda _\mathbf{x} ^{(1)}=\sqrt{|A_{1,\mathbf{x} }^{(1)}|^2+|A_{2,\mathbf{x} }^{(1)}|^2}\), \(\varLambda _\mathbf{x} =\mathrm{diag}(\varLambda _{\mathbf{x} ,+},\varLambda _{\mathbf{x} ,-})\), \(\varLambda _{\mathbf{x} ,\pm }=V_\mathbf{x} ^{(1)}\pm \lambda _\mathbf{x} ^{(1)}\), and \(P_\mathbf{x} =I_2\) if \(\lambda _\mathbf{x} ^{(1)}=0\) and otherwise

$$\begin{aligned} P_\mathbf{x} =\begin{pmatrix}\frac{1}{\sqrt{2}}&{}\frac{A_{1,\mathbf{x} }^{(1)}-iA_{2,\mathbf{x} }^{(1)}}{\sqrt{2}\lambda _\mathbf{x} ^{(1)}}\\ -\frac{A_{1,\mathbf{x} }^{(1)}+iA_{2,\mathbf{x} }^{(1)}}{\sqrt{2}\lambda _\mathbf{x} ^{(1)}}&{}\frac{1}{\sqrt{2}} \end{pmatrix}. \end{aligned}$$
(6.34)

For the Dirac equation (1.9) in 3D, i.e. we take \(d=3\) in (1.9). The problem is truncated on \(\varOmega =(a_1, b_1)\times (a_2, b_2)\times (a_3, b_3)\) with mesh sizes \(h_1=(b_1-a_1)/M_1\), \(h_2=(b_2-a_2)/M_2\) and \(h_3=(b_3-a_3)/M_3\) (\(M_1,M_2,M_3\) three even positive integers) in x-, y- and z-direction, respectively. The wave function \(\varPsi \) is a four-component vector, and the matrix \(\varGamma _l\) in (3.8) will be replaced by \(\varGamma _{jkl}\) as:

$$\begin{aligned} \varGamma _{jkl}=\begin{pmatrix} 1 &{} 0 &{} \varepsilon \mu _l^{(3)} &{} \varepsilon \mu _j^{(1)}-i\varepsilon \mu _k^{(2)} \\ 0 &{} 1 &{} \varepsilon \mu _j^{(1)}+i\varepsilon \mu _k^{(2)} &{} -\varepsilon \mu _l^{(3)} \\ \varepsilon \mu _l^{(3)} &{} \varepsilon \mu _j^{(1)}-i\varepsilon \mu _k^{(2)} &{} -1 &{} 0 \\ \varepsilon \mu _j^{(1)}+i\varepsilon \mu _k^{(2)} &{} -\varepsilon \mu _l^{(3)} &{} 0 &{} -1 \\ \end{pmatrix}, \end{aligned}$$
(6.35)

where \(-\frac{M_1}{2}\le j\le \frac{M_1}{2}-1,-\frac{M_2}{2}\le k\le \frac{M_2}{2}-1,-\frac{M_3}{2}\le l\le \frac{M_3}{2}-1\) and

$$\begin{aligned} \mu _j^{(1)}=\frac{2j\pi }{b_1-a_1},\quad \mu _k^{(2)}=\frac{2k\pi }{b_2-a_2},\quad \mu _l^{(3)}=\frac{2l\pi }{b_3-a_3}. \end{aligned}$$
(6.36)

The eigenvalues of \(\varGamma _{jkl}\) are

$$\begin{aligned} \delta _{jkl}, \delta _{jkl}, -\delta _{jkl}, -\delta _{jkl},\quad \text {with}\quad \delta _{jkl}=\sqrt{1+\varepsilon ^2\left| \mu _j^{(1)}\right| ^2+\varepsilon ^2\left| \mu _{k}^{(2)}\right| ^2+ \varepsilon ^2\left| \mu _l^{(3)}\right| ^2}. \end{aligned}$$

The corresponding eigenvectors are

$$\begin{aligned} \mathbf {v}^{(1)}_{jkl}=\begin{pmatrix} 1+\delta _{jkl}\\ 0\\ \varepsilon \mu _l^{(3)}\\ \varepsilon \mu _j^{(1)}+i\varepsilon \mu _{k}^{(2)}\end{pmatrix},\quad \mathbf {v}^{(2)}_{jkl}=\begin{pmatrix} 0\\ 1+\delta _{jkl}\\ \varepsilon \mu _j^{(1)}-i\varepsilon \mu _{k}^{(2)}\\ -\varepsilon \mu _l^{(3)}\end{pmatrix}, \end{aligned}$$

and

$$\begin{aligned} \mathbf {v}^{(3)}_{jkl}=\begin{pmatrix} -\varepsilon \mu _l^{(3)}\\ -\varepsilon \mu _j^{(1)}-i\varepsilon \mu _{k}^{(2)}\\ 1+\delta _{jkl}\\ 0\end{pmatrix},\quad \mathbf {v}^{(4)}_{jkl}=\begin{pmatrix} -\varepsilon \mu _j^{(1)}+i\varepsilon \mu _{k}^{(2)}\\ \varepsilon \mu _l^{(3)}\\ 0\\ 1+\delta _{jkl}\end{pmatrix}. \end{aligned}$$

Then the Schur decomposition \(\varGamma _{jkl}=Q_{jkl}D_{jkl}Q^*_{jkl}\) is given as

$$\begin{aligned} D_{jkl}=\mathrm {diag}(\delta _{jkl},\delta _{jkl},-\delta _{jkl},-\delta _{jkl}),\quad Q_{jkl}=\frac{1}{\sqrt{2\delta _{jkl}(1+\delta _{jkl})}}\left( \mathbf {v}^{(1)}_{jkl}, \mathbf {v}^{(2)}_{jkl},\mathbf {v}^{(3)}_{jkl},\mathbf {v}^{(4)}_{jkl}\right) . \end{aligned}$$

The matrix \(\int _{t_n}^{t_{n+1}}G(t,\mathbf{x} )dt\) in (4.5) becomes \(\int _{t_n}^{t_{n+1}}G(t,\mathbf{x} )dt\) and the Schur decomposition \(\int _{t_n}^{t_{n+1}}G(t,\mathbf{x} )dt=P_\mathbf{x} \varLambda _\mathbf{x} P_{\mathbf{x} }^*\) with \(V_\mathbf{x} ^{(1)}=\int _{t_n}^{t_{n+1}}V(t,\mathbf{x} )dt\), \(A_{l,\mathbf{x} }^{(1)}=\int _{t_n}^{t_{n+1}}A_l(t,\mathbf{x} )dt\) for \(l=1,2,3\), \(\lambda _\mathbf{x} ^{(1)}=\sqrt{|A_{1,\mathbf{x} }^{(1)}|^2+|A_{2,\mathbf{x} }^{(1)}|^2+|A_{3,\mathbf{x} }^{(1)}|^2}\), \(\varLambda _\mathbf{x} =\mathrm{diag}(\varLambda _{\mathbf{x} ,+},\varLambda _{\mathbf{x} ,+},\varLambda _{\mathbf{x} ,-},\varLambda _{\mathbf{x} ,-})\), \(\varLambda _{\mathbf{x} ,\pm }=V_\mathbf{x} ^{(1)}\pm \lambda _\mathbf{x} ^{(1)}\), and \(P_\mathbf{x} =I_4\) if \(\lambda _\mathbf{x} ^{(1)}=0\) and otherwise \(P_{\mathbf{x} }=(\mathbf{u}_\mathbf{x} ^{(1)},\mathbf{u}_\mathbf{x} ^{(2)},\mathbf{u}_\mathbf{x} ^{(3)},\mathbf{u}_\mathbf{x} ^{(4)})\),

$$\begin{aligned} \mathbf{u}_\mathbf{x} ^{(1)}=\begin{pmatrix}\frac{-A_{1,\mathbf{x} }^{(1)}+iA_{2,\mathbf{x} }^{(1)}}{\sqrt{2}\lambda _\mathbf{x} ^{(1)}}\\ \frac{A_{3,\mathbf{x} }^{(1)}}{\sqrt{2}\lambda _\mathbf{x} ^{(1)}}\\ 0\\ \frac{1}{\sqrt{2}}\end{pmatrix},\quad \mathbf{u}_\mathbf{x} ^{(2)}=\begin{pmatrix}\frac{-A_{3,\mathbf{x} }^{(1)}}{\sqrt{2}\lambda _\mathbf{x} ^{(1)}}\\ \frac{-A_{1,\mathbf{x} }^{(1)}-iA_{2,\mathbf{x} }^{(1)}}{\sqrt{2}\lambda _\mathbf{x} ^{(1)}}\\ \frac{1}{\sqrt{2}}\\ 0\end{pmatrix}, \end{aligned}$$

and

$$\begin{aligned} \mathbf{u}_\mathbf{x} ^{(3)}=\begin{pmatrix}0\\ \frac{1}{\sqrt{2}}\\ \frac{A_{1,\mathbf{x} }^{(1)}-iA_{2,\mathbf{x} }^{(1)}}{\sqrt{2}\lambda _\mathbf{x} ^{(1)}}\\ \frac{-A_{3,\mathbf{x} }}{\sqrt{2}\lambda _\mathbf{x} ^{(1)}}\end{pmatrix},\quad \mathbf{u}_\mathbf{x} ^{(4)}=\begin{pmatrix}\frac{1}{\sqrt{2}}\\ 0\\ \frac{A_{3,\mathbf{x} }^{(1)}}{\sqrt{2}\lambda _\mathbf{x} ^{(1)}}\\ \frac{A_{1,\mathbf{x} }^{(1)}+iA_{2,\mathbf{x} }^{(1)}}{\sqrt{2}\lambda _\mathbf{x} ^{(1)}}\end{pmatrix} . \end{aligned}$$

For the Dirac equation (1.9) in 2D, we simply let \(\mu _{l}^{(3)}=0\), \(A_3(t,\mathbf{x} )\equiv 0\) in the above 3D case; and for the Dirac equation (1.9) in 1D, we let \(\mu _{k}^{(2)}=\mu _{l}^{(3)}=0\), \(A_2(t,\mathbf{x} )=A_3(t,\mathbf{x} )\equiv 0\) in the above 3D case. Then the sEWI-FP (3.18–3.19) and TSFP (4.4) can be designed accordingly for the Dirac equation (1.9) in 2D and 1D.

Rights and permissions

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Bao, W., Cai, Y., Jia, X. et al. Numerical Methods and Comparison for the Dirac Equation in the Nonrelativistic Limit Regime. J Sci Comput 71, 1094–1134 (2017). https://doi.org/10.1007/s10915-016-0333-3

Download citation

  • Received:

  • Revised:

  • Accepted:

  • Published:

  • Issue Date:

  • DOI: https://doi.org/10.1007/s10915-016-0333-3

Keywords

Search

Navigation