Simulation for formation process of atomic orbitals by the finite difference time domain method based on the eight-element Dirac equation

Hideki Mutoh

doi:10.1515/phys-2022-0262

Article Open Access

Simulation for formation process of atomic orbitals by the finite difference time domain method based on the eight-element Dirac equation

Hideki Mutoh

Published/Copyright: July 13, 2023

Published by

Become an author with De Gruyter Brill

Submit Manuscript Author Information Explore this Subject

From the journal Open Physics Volume 21 Issue 1

Abstract

Using the finite difference time domain (FDTD) method based on the eight-element Dirac equation, we found that a stable Dirac field wave packet with low velocity can be created without explicit consideration of Zitterbewegung (the rapid oscillatory motion of elementary particles), which is difficult in one-dimensional simulations. Furthermore, we successfully simulated the formation process of atomic orbitals for the first time without any physical approximations by calculating the eight-element Dirac field propagation in the central electric force potential. Initially, a small unstable orbital appears, which rapidly grows and results in a large stable orbital with a radius equal to the Bohr radius divided by the atomic number, as given by the solution of the Schrödinger equation. The FDTD calculation based on the conventional four-element Dirac equation cannot produce such reasonable orbitals owing to the spatial asymmetry of the 4 × 4 Dirac matrices. This method has the potential to be used for transient analyses of not only atomic or molecular orbitals but also interactions among elementary particles.

Keywords: Dirac equation; Maxwell’s equations; FDTD method; Zitterbewegung; atomic orbital

1 Introduction

Transient analyses of spatial distributions of atomic and molecular orbitals are crucial for understanding chemical reactions. Experimental investigations of molecular orbitals and their time-dependent spatial distributions have been conducted using Penning ionization electron spectroscopy [1–10], photoelectron spectroscopy with angular distribution [11–15], and electron momentum spectroscopy [16–21]. By contrast, theoretical studies of molecular orbitals have primarily focused on the Schrödinger equation [22–30]. Considering that the orbital spatial distributions change rapidly occurring on a timescale comparable to the orbital length divided by the speed of light, it is essential to solve the Dirac equation instead of the Schrödinger equation for more accurate analyses of the time-dependent spatial distributions of atomic and molecular orbitals. The finite difference time domain (FDTD) method [31,32], which can be used for transient analysis of electromagnetic fields, could also be used for analyzing the Dirac field [33,34] as the Dirac equation is almost equivalent to Maxwell’s equations, except for electron charge and mass [35,36]. In this study, we demonstrate that the FDTD method based on the eight-element Dirac equation successfully calculates the time-dependent Dirac field. Moreover, it reveals the first-ever simulation of the formation process of atomic orbitals, starting from the initial state of a free electron and an atomic nucleus without any physical approximations. This is achieved by calculating the eight-element Dirac field propagation in the central electric force potential.

2 Eight-element Dirac equation

The Dirac equation is given by:

(1) ( i ℏ γ μ ∂ μ − m c ) ψ = 0 ,

where ℏ is the Planck constant, ψ is a wave function vector consisting of four components, m is the mass, c is the light speed in free space, and γ μ are the gamma matrices that satisfy the following equation [37–40]:

(2) 1 2 ( γ μ γ ν + γ ν γ μ ) = η μ ν ≡ 1 : μ = ν = 0 0 : μ ≠ ν − 1 : μ = ν = 1 , 2 , 3 .

For example, gamma matrices are given by:

(3) γ 0 = 1 0 0 0 0 1 0 0 0 0 − 1 0 0 0 0 − 1 , γ 1 = 0 0 0 1 0 0 1 0 0 − 1 0 0 − 1 0 0 0 , γ 2 = 0 0 0 − i 0 0 i 0 0 i 0 0 − i 0 0 0 , γ 3 = 0 0 1 0 0 0 0 − 1 − 1 0 0 0 0 1 0 0 .

These matrices have spatial asymmetry, which means that only one or two of γ 1 , γ 2 , and γ 3 include imaginary numbers. Since it is difficult to calculate the propagation of the Dirac fields by the FDTD method in the aforementioned case as shown later, the Dirac equation should be extended to consist of spatially symmetric matrices as follows.

(4) γ ¯ μ ∂ μ + m c ℏ γ ¯ 4 ψ = 0 ,

where γ ¯ 1 , γ ¯ 2 , and γ ¯ 3 are the real number matrices and γ ¯ μ ( μ = 0 , 1 , 2 , 3 , 4 ) satisfy

(5) 1 2 ( γ ¯ μ γ ¯ ν + γ ¯ ν γ ¯ μ ) = η ¯ μ ν ≡ 1 : μ = ν = 0 , 4 0 : μ ≠ ν − 1 : μ = ν = 1 , 2 , 3 ,

where γ ¯ μ cannot be represented by 4 × 4 matrices. Eqs. (4) and (5) give the Klein-Gordon equation:

(6) γ ¯ μ ∂ μ + m c ℏ γ ¯ 4 2 ψ = □ + m 2 c 2 ℏ 2 ψ = 0 ,

where □ is d’Alembertian defined by □ ≡ ∂ 0 2 − ∇ 2 . To obtain γ ¯ μ , we introduce the following spatially symmetric delta matrices:

(7) δ 0 = i 0 0 0 0 i 0 0 0 0 i 0 0 0 0 i , δ 1 = 0 1 0 0 − 1 0 0 0 0 0 0 − 1 0 0 1 0 , δ 2 = 0 0 1 0 0 0 0 1 − 1 0 0 0 0 − 1 0 0 , δ 3 = 0 0 0 1 0 0 − 1 0 0 1 0 0 − 1 0 0 0 .

Now, we define matrices γ ¯ μ ( μ = 0 , 1 , 2 , 3 ) and γ ¯ 4 as:

(8) γ ¯ μ ≡ 0 δ μ ∗ δ μ 0 , γ ¯ 4 ≡ I 4 0 0 − I 4 ,

where I 4 is the 4 × 4 unit matrix. Since γ ¯ μ ( μ = 0 , 1 , 2 , 3 , 4 ) satisfy Eq. (5), we obtain the eight-element Dirac equation of Eq. (4). When we introduce the four-element wave functions ψ A and ψ B , which satisfy

(9) ψ = ψ A ψ B ,

Eq. (4) is rewritten as:

(10) m c ∕ ℏ δ μ ∗ ∂ μ δ μ ∂ μ − m c ∕ ℏ ψ A ψ B = 0 .

Then, we obtain

(11) i ∂ 0 ψ A + δ k ∂ k ψ A − m c ℏ ψ B = 0 , − i ∂ 0 ψ B + δ k ∂ k ψ B + m c ℏ ψ A = 0 .

When † means the Hermite conjugate, we obtain

(12) − i ∂ 0 ψ A † − ∂ k ψ A † δ k − m c ℏ ψ B † = 0 , i ∂ 0 ψ B † − ∂ k ψ B † δ k + m c ℏ ψ A † = 0 .

Here, we define four current C μ as follows:

(13) C 0 ≡ ψ A † ψ A + ψ B † ψ B = ψ † ψ , C k ≡ − i ( ψ A † δ k ψ A − ψ B † δ k ψ B ) .

Then, Eqs. (11)–(13) give

(14) ∂ μ C μ = 0 .

Now, we define ψ ¯ † as:

(15) ψ ¯ † ≡ ψ † γ ¯ 0 .

Then, we obtain

(16) C μ = ψ ¯ † γ ¯ μ ψ .

Therefore, C 0 and C k can be regarded as probability density and probability current density, respectively.

Next, we consider spin. When we introduce ψ + and ψ − as

(17) ψ + ≡ ψ A + ψ B , ψ − ≡ ψ A − ψ B ,

Eq. (11) is rewritten as:

(18) i ∂ 0 ψ + + δ k ∂ k ψ − − m c ℏ ψ + = 0 , i ∂ 0 ψ − + δ k ∂ k ψ + + m c ℏ ψ − = 0 .

Under vector potential A and scalar one ϕ , we should substitute ∂ k + i e A k ∕ ℏ c for ∂ k and ∂ 0 − i e ϕ ∕ ℏ c for ∂ 0 in the equations. When we assume that E is the total energy and the time dependence of the wave functions is exp ( − i E t ∕ ℏ ) , we could substitute − i E ∕ ℏ c for ∂ 0 as:

(19) E ℏ c + e ℏ c ϕ − m c ℏ ψ + + ∂ k ∂ k + i e ℏ c A k ψ − = 0 , E ℏ c + e ℏ c ϕ + m c ℏ ψ − + ∂ k ∂ k + i e ℏ c A k ψ + = 0 .

Considering the nonrelativistic condition, E is given as:

(20) E = E N R + m c 2 ,

where E N R is the nonrelativistic energy. When we assume ∣ E N R ∣ ≪ m c 2 and ∣ e ϕ ∣ ≪ m c 2 , Eq. (18) gives

(21) ψ − ≈ − ℏ 2 m c δ k ∂ k + i e ℏ c A k ψ + .

Then, Eqs (19) and (21) give

(22) ( E N R + e ϕ ) ψ + = ℏ 2 2 m δ k ∂ k + i e ℏ c A k δ k ∂ k + i e ℏ c A k ψ + = − ℏ 2 2 m ∇ + i e ℏ c A 2 + i e ℏ 2 m c δ k B k ψ + .

Therefore, the second term of the right side of the aforementioned equation shows the magnetic moment − i ℏ ∕ 2 m c . When we define s k as:

(23) s k ≡ − i ℏ 2 δ k ,

we obtain

(24) s 1 s 2 − s 2 s 1 = i ℏ s 3 , s 2 s 3 − s 3 s 2 = i ℏ s 1 , s 3 s 1 − s 1 s 3 = i ℏ s 2 , ( s 1 ) 2 + ( s 2 ) 2 + ( s 3 ) 2 = 3 4 ℏ 2 .

Therefore, s is a spin operator with a spin quantum number of 1/2. Eq. (22) can be regarded as the wave equation of a particle with a spin quantum number of 1/2 under vector potential A and scalar one ϕ .

3 Dirac field propagation analysis by FDTD method

The FDTD method is one of the simplest methods for transient analysis of field propagation, because it can give field spatial distribution dependence on time by only substituting a pair of field vectors each other to discretized equations starting from a given initial state. Since the Dirac equation is quite similar to Maxwell’s equations, the Dirac field propagation could be calculated by the FDTD method, which is popularly used for propagation analysis of electromagnetic field. We compared the calculation results of the FDTD method based on the 1D-like two-element, the conventional four-element, and the extended eight-element Dirac equations. Figure 1 shows the analyzed structure consisting of a cube with a side length L , where the origin exists at the center of the cube. The Dirac field wave packet is created by the vibration of the top surface of the structure. We adopt the following wave function as the initial wave packet of the Dirac field, which has finite values in the whole region and satisfies the Klein–Gordon equation:

(25) ψ = sin ( k r ′ ) exp ( − i ω t ′ ) k r ′ ,

where r ′ and t ′ satisfy the following equations of Lorentz transformation using the wave packet moving velocity v and β ≡ v ∕ c .

(26) r ′ = x 2 + y 2 + ( z − v t ) 2 1 − β 2 , t ′ = 1 1 − β 2 1 − v z c 2 .

k and ω are wave number and angular frequency, respectively, and satisfy

(27) k 2 = ω 2 c 2 − m 2 c 2 ℏ 2 .

Figure 1

Analyzed structure.

3.1 Discretization for the two- and four-element Dirac equation

The four-element Dirac equation of Eq. (1) is rewritten as:

(28) ∂ 0 + i m c ℏ ψ 0 + ( ∂ 1 − i ∂ 2 ) ψ 3 + ∂ 3 ψ 2 = 0 , ∂ 0 + i m c ℏ ψ 1 + ( ∂ 1 + i ∂ 2 ) ψ 2 − ∂ 3 ψ 3 = 0 , − ∂ 0 + i m c ℏ ψ 2 + ( − ∂ 1 + i ∂ 2 ) ψ 1 − ∂ 3 ψ 0 = 0 , − ∂ 0 + i m c ℏ ψ 3 + ( − ∂ 1 − i ∂ 2 ) ψ 0 + ∂ 3 ψ 1 = 0 .

Here, ψ μ is written by the sum of real and imaginary parts as:

(29) ψ μ = ψ r μ + i ψ i μ ,

where ψ r μ and ψ i μ are the real functions of space-time. The aforementioned equations are discretized by the FDTD method, assuming that i , j , k , and n of ψ r : i , j , k , n μ and ψ i : i , j , k , n μ denote the grid address along x , y , z , and time axes, respectively. The aforementioned functions are defined at the position in the cell of the FDTD method as shown in Figure 2(a), where the real and imaginary parts are defined at different positions and the element positions are spatially asymmetric. The discretized equations of Eq. (28) are

ξ μ η μ ψ r : i , j , k , n + 1 − ζ μ 2 μ − ψ r : i , j , k , n − ζ μ 2 μ c △ t − m c 4 ℏ ∑ l = 0 1 ∑ m = 0 1 ψ i : i − ξ μ l , j − m , k , n + ζ μ μ + ∑ ν = 0 3 ξ μ γ μ ν 1 ψ r : i , j , k , n + ζ ν 2 ν − ψ r : i − ξ μ , j , k , n + ζ ν 2 ν △ x + i γ μ ν 2 ψ i : i , j , k , n + ζ ν 2 ν − ψ i : i , j − 1 , k , n + ζ ν 2 ν △ y + η μ γ μ ν 3 ψ r : i , j , k , n + ζ ν 2 ν − ψ r : i , j , k − η μ , n + ζ ν 2 ν △ z = 0 ,

(30) ξ μ η μ ψ i : i , j , k , n + 1 − ζ μ 2 μ − ψ i : i , j , k , n − ζ μ 2 μ c △ t + m c 4 ℏ ∑ l = 0 1 ∑ m = 0 1 ψ r : i + ξ μ l , j + m , k , n + ζ μ μ + ∑ ν = 0 3 ξ μ γ μ ν 1 ψ i : i + ξ μ , j , k , n + ζ ν 2 ν − ψ i : i , j , k , n + ζ ν 2 ν △ x − i γ μ ν 2 ψ r : i , j + 1 , k , n + ζ ν 2 ν − ψ r : i , j , k , n + ζ ν 2 ν △ y + η μ γ μ ν 3 ψ i : i , j , k , n + ζ ν 2 ν − ψ i : i , j , k − η μ , n + ζ ν 2 ν △ z = 0 ,

where △ x , △ y , △ z , and △ t denote the grid spaces of x , y , z , and time axes, respectively, and

(31) ξ μ ≡ 1 : μ = 0 , 2 − 1 : μ = 1 , 3 , η μ ≡ 1 : μ = 0 , 3 − 1 : μ = 1 , 2 , ζ μ ≡ 0 : μ = 0 , 1 1 : μ = 2 , 3 .

Therefore,

(32) ψ r : i , j , k , n + 1 − ζ μ 2 μ = ψ r : i , j , k , n − ζ μ 2 μ + ξ μ η μ α 0 4 ∑ l = 0 1 ∑ m = 0 1 ψ i : i − l , j − m , k , n − ζ μ 2 μ − ∑ ν = 0 3 η μ γ μ ν 1 α 1 ψ r : i , j , k , n + ζ ν 2 ν − ψ r : i − 1 , j , k , n + ζ ν 2 ν + α 2 ψ i : i , j , k , n + ζ ν 2 ν − ψ i : i , j − 1 , k , n + ζ ν 2 ν + α 3 ξ μ γ μ ν 3 ψ r : i , j , k , n + ζ ν 2 ν − ψ r : i , j , k − 1 , n + ζ ν 2 ν , ψ i : i , j , k , n + 1 − ζ μ 2 ν = ψ i : i , j , k , n − ζ μ 2 ν − ξ μ η μ α 0 4 ∑ l = 0 1 ∑ m = 0 1 ψ i : i + l , j + m , k , n − ζ μ 2 0 − ∑ ν = 0 3 η μ γ μ ν 1 α 1 ψ i : i + 1 , j , k , n + ζ ν 2 ν − ψ i : i , j , k , n + ζ ν 2 ν − α 2 ψ r : i , j + 1 , k , n + ζ ν 2 ν − ψ r : i , j , k , n + ζ ν 2 ν + α 3 ξ μ γ μ ν 3 ψ i : i , j , k , n + ζ ν 2 ν − ψ i : i , j , k − 1 , n + ζ ν 2 ν ,

where

(33) α 0 ≡ m c 2 △ t ℏ , α 1 ≡ c △ t △ x , α 2 ≡ c △ t △ y , α 3 ≡ c △ t △ z .

$Figure 2 Definition position of the Dirac field elements in the cell of the FDTD method. (a) The four-element Dirac field and (b) the eight-element Dirac field, where the white, red, blue, and green circles denote ψ 0 {\psi }^{0} , ψ 1 {\psi }^{1} , ψ 2 {\psi }^{2} , and ψ 3 {\psi }^{3} , respectively.$

Figure 2

Definition position of the Dirac field elements in the cell of the FDTD method. (a) The four-element Dirac field and (b) the eight-element Dirac field, where the white, red, blue, and green circles denote ψ 0 , ψ 1 , ψ 2 , and ψ 3 , respectively.

When we consider the one-dimensional analysis of the Dirac field propagation along z -axis, assuming ∂ 1 ψ μ = ∂ 2 ψ μ = 0 , Eq. (28) is rewritten as:

(34) ∂ 0 + i m c ℏ ψ 0 + ∂ 3 ψ 2 = 0 , − ∂ 0 + i m c ℏ ψ 2 − ∂ 3 ψ 0 = 0 , ∂ 0 + i m c ℏ ψ 1 − ∂ 3 ψ 3 = 0 , − ∂ 0 + i m c ℏ ψ 3 + ∂ 3 ψ 1 = 0 .

Since we can define the real and imaginary parts at the same position in this case, we obtain the discretized one-dimensional Dirac field as follows:

(35) ψ i , j , k , n + 1 0 = ( 1 − i α 0 ) ψ i , j , k , n 0 − α 3 ψ i , j , k , n + 1 2 2 − ψ i , j , k − 1 , n + 1 2 2 , ψ i , j , k , n + 1 2 2 = ( 1 + i α 0 ) ψ i , j , k , n − 1 2 2 − α 3 ( ψ i , j , k + 1 , n 0 − ψ i , j , k , n 0 ) , ψ i , j , k , n + 1 3 = ( 1 + i α 0 ) ψ i , j , k , n 3 + α 3 ψ i , j , k , n + 1 2 1 − ψ i , j , k − 1 , n + 1 2 1 , ψ i , j , k , n + 1 2 1 = ( 1 − i α 0 ) ψ i , j , k , n − 1 2 1 + α 3 ( ψ i , j , k + 1 , n 3 − ψ i , j , k , n 3 ) .

The one-dimensional Dirac field propagation along z -axis can be calculated by the two elements, the pair of ψ 0 and ψ 2 or that of ψ 1 and ψ 3 , because the upper two and lower two equations of Eqs. (34) and (35) are independent.

3.2 Discretization for the eight-element Dirac equation

As same as the four-element case, we can obtain discretized equations for the eight-element Dirac field ψ − : i , j , k , n μ and ψ + : i , j , k , n μ defined at the positions in Figure 2(b), where the real and imaginary parts can be defined at the same position and the element positions are spatially symmetric. The discretized equations of Eq. (18) are

(36) ψ + : i , j , k , n + 1 μ − ψ + : i , j , k , n μ c △ t + i m c ℏ ψ + : i , j , k , n μ − i λ μ ∑ ν = 0 3 λ ν δ μ ν 1 ψ − : i + λ μ λ ν , j , k , n + 1 2 ν − ψ − : i , j , k , n + 1 2 ν △ x + δ μ ν 2 ψ − : i , j + λ μ λ ν , k , n + 1 2 ν − ψ − : i , j , k , n + 1 2 ν △ y + δ μ ν 3 ψ − : i , j , k + λ μ λ ν , n + 1 2 ν − ψ − : i , j , k , n + 1 2 ν △ z = 0 , ψ − : i , j , k , n + 1 2 μ − ψ − : i , j , k , n − 1 2 μ c △ t − i m c ℏ ψ − : i , j , k , n − 1 2 μ − i λ μ ∑ ν = 0 3 λ ν δ μ ν 1 ψ + : i , j , k , n ν − ψ + : i − λ μ λ ν , j , k , n ν △ x + δ μ ν 2 ψ + : i , j , k , n ν − ψ + : i , j − λ μ λ ν , k , n ν △ y + δ μ ν 3 ψ + : i , j , k , n ν − ψ + : i , j , k − λ μ λ ν , n ν △ z = 0 ,

where

(37) λ μ ≡ 1 : μ = 1 , 2 , 3 − 1 : μ = 0 .

Therefore,

(38) ψ + : i , j , k , n + 1 μ = ( 1 − i α 0 ) ψ + : i , j , k , n μ − i λ μ ∑ ν = 0 3 λ ν α 1 δ μ ν 1 ψ − : i + λ μ λ ν , j , k , n + 1 2 ν − ψ − : i , j , k , n + 1 2 ν + α 2 δ μ ν 2 ψ − : i , j + λ μ λ ν , k , n + 1 2 ν − ψ − : i , j , k , n + 1 2 ν + α 3 δ μ ν 3 ψ − : i , j , k + λ μ λ ν , n + 1 2 ν − ψ − : i , j , k , n + 1 2 ν , ψ − : i , j , k , n + 1 2 μ = ( 1 + i α 0 ) ψ − : i , j , k , n − 1 2 μ + i λ μ ∑ ν = 0 3 λ ν { α 1 δ μ ν 1 ( ψ + : i , j , k , n ν − ψ + : i − λ μ λ ν , j , k , n ν ) + α 2 δ μ ν 2 ( ψ + : i , j , k , n ν − ψ + : i , j − λ μ λ ν , k , n ν ) + α 3 δ μ ν 3 ( ψ + : i , j , k , n ν − ψ + : i , j , k − λ μ λ ν , n ν ) } .

3.3 Comparison among the two-, four-, and eight-element Dirac field propagation

Figure 3 shows the wave packet shape dependence on the propagation time for the two-, four-, and eight-element fields in the cube of L = 10 pm . Figure 4(a) shows that the moving distance of the wave packet peak is equal to the product of the wave packet velocity and the propagation time in the four- and eight-element cases, although the wave packet speed by 1D-like calculation in the two-element case is equal to the light speed in free space. We found that explicit consideration of Zitterbewegung [41,42] is unnecessary for the three-dimensional calculation. As shown in Figure 4(b), the wave packet peak intensity decreases with increasing propagation time in the four-element case, although it is almost constant except early period in the eight-element case. The instability of the four-element Dirac field propagation seems to be caused by the asymmetric definition position of the field elements shown in Figure 2(a).

Figure 3

Time dependence of the Dirac field intensity distribution. (a), (b), (c), and (d) are the two-element field at t = 0.01 , 0.02, 0.03, and 0.04 as, (e), (f), (g), and (h) are the four-element field at t = 0.01 , 0.02, 0.03, and 0.04 as, (i), (j), (k), and (l) are the eight-element field at t = 0.01 , 0.02, 0.03, and 0.04 as, respectively.

$Figure 4 Comparison among the two-, four-, and eight-element Dirac field propagation with β = 0.5 \beta =0.5 , (a) peak z z -coordinate dependence on time and (b) peak intensity dependence on time, where the blue, green, and red lines show the two-, four-, and eight-element fields, respectively, and the broken black line denotes the theoretical value.$

Figure 4

Comparison among the two-, four-, and eight-element Dirac field propagation with β = 0.5 , (a) peak z -coordinate dependence on time and (b) peak intensity dependence on time, where the blue, green, and red lines show the two-, four-, and eight-element fields, respectively, and the broken black line denotes the theoretical value.

4 Simulation for formation process of atomic orbitals

In the electric central force potential ϕ = Z e ∕ 4 π ε 0 r generated by an atomic nucleus with its atomic number of Z , the eight-element Dirac equations are given by:

(39) i ∂ 0 + e ϕ ℏ c + δ k ∂ k ψ A − m c ℏ ψ B = 0 , − i ∂ 0 − e ϕ ℏ c + δ k ∂ k ψ B + m c ℏ ψ A = 0 ,

and

(40) i ∂ 0 + e ϕ ℏ c ψ + + δ k ∂ k ψ − − m c ℏ ψ + = 0 , i ∂ 0 + e ϕ ℏ c ψ − + δ k ∂ k ψ + + m c ℏ ψ − = 0 .

In the discretized equations of Eqs. (36) and (38), m c is replaced with m c ± e ϕ ∕ c . We calculated the eight-element Dirac field propagation in the cube of L = 10 pm based on the aforementioned equations and assuming an injected electron velocity parameter β of 0.5. Figure 5(a)–(c) show the log-scale intensity on x z -plane at early propagation times of 0.03 as, 0.0325 as and 0.035 as. We found that after t = 0.03 as, a cone-shaped peak appears at the origin and rapidly grows with increasing propagation time, while the injected wave packet peak intensity simultaneously decreases. During this phase, the cone-shaped wave function surrounding the origin is proportional to exp ( − r ∕ r init ) , where r init denotes the initial orbital radius. As shown in Figure 5(d), r init is inversely proportional to the atomic number Z and considerably smaller than the analytical value of r B ∕ Z provided by the Schrödinger equation, where r B denotes the Bohr radius. This is explained by the following equations as derived from Eq. (39):

(41) □ − 2 i e ϕ ℏ c ∂ 0 − e 2 ϕ 2 ℏ 2 c 2 + m 2 c 2 ℏ 2 + e ℏ c ( δ k ∂ k ϕ ) ψ A = 0 , □ − 2 i e ϕ ℏ c ∂ 0 − e 2 ϕ 2 ℏ 2 c 2 + m 2 c 2 ℏ 2 − e ℏ c ( δ k ∂ k ϕ ) ψ B = 0 .

Assuming that the time-dependent wave function is proportional to exp ( − r ∕ r init − i ω t ) , Eq. (41) gives the terms of r − n exp ( − r ∕ r init − i ω t ) for n = 0 , 1 , 2 . The term corresponding to n = 0 provides

(42) − ω 2 c 2 − 1 r init 2 + m 2 c 2 ℏ 2 = 0 .

Therefore, if r init is smaller than ℏ ∕ mc ≈ 0.39 pm, ω 2 becomes negative, resulting in an imaginary ω value. This implies that the absolute value of the wave function exp ( − r ∕ r init − i ω t ) will exponentially increase over time, eventually leading to an overflow of peak intensity. To prevent this overflow, we imposed a limit of 1 0 3 for the absolute value of each wave function element within the spherical region of r < 0.5 pm.

$Figure 5 Time dependence of the eight-element Dirac field intensity in the potential of Z e ∕ 4 π ε 0 r Ze/4\pi {\varepsilon }_{0}r with Z = 50 Z=50 for (a) t t = 0.03 as, (b) t t = 0.0325 as, and (c) t t = 0.035 as, and (d) the initial orbital radius dependence on 1 ∕ Z 1/Z .$

Figure 5

Time dependence of the eight-element Dirac field intensity in the potential of Z e ∕ 4 π ε 0 r with Z = 50 for (a) t = 0.03 as, (b) t = 0.0325 as, and (c) t = 0.035 as, and (d) the initial orbital radius dependence on 1 ∕ Z .

Figure 6 shows the log-scale Dirac field intensity as a function of the propagation time on x z -plane in the cube of L = 20 pm, considering the aforementioned limit. At t = 0.05 as, a small cone-shaped orbital appears; similar to that observed in Figure 5, the radius of the orbital expands with time after 0.055 as. Beyond t = 0.09 as, the high-intensity region of the wave function extends across the entire calculation structure, with the orbital shape stabilizing after t = 0.17 as in the region of r > 3 pm. At the final moment of t = 0.2 as, the outer part of the wave function appears to be proportional to exp ( − r ∕ r final ) , where r final denotes the radius of the final orbital. In Figure 7, the solid line illustrates r final dependence on 1 ∕ Z , while the broken line represents the r final = r B ∕ Z . Considering that 1s atomic orbitals are known to be proportional to exp ( − Z r ∕ r B ) [43], our calculations can be considered accurate, as r final closely aligns with the theoretical value provided by the Schrödinger equation.

$Figure 6 Time dependence of the eight-element Dirac field intensity in the potential of Z e ∕ 4 π ε 0 r Ze/4\pi {\varepsilon }_{0}r , with Z = 50 Z=50 . (a) t t = 0.045 as, (b) t t = 0.05 as, (c) t t = 0.055 as, (d) t t = 0.06 as, (e) t t = 0.065 as, (f) t t = 0.07 as, (g) t t = 0.08 as, (h) t t = 0.09 as, (i) t t = 0.12 as, (j) t t = 0.145 as, (k) t t = 0.17 as, and (l) t t = 0.2 as.$

Figure 6

Time dependence of the eight-element Dirac field intensity in the potential of Z e ∕ 4 π ε 0 r , with Z = 50 . (a) t = 0.045 as, (b) t = 0.05 as, (c) t = 0.055 as, (d) t = 0.06 as, (e) t = 0.065 as, (f) t = 0.07 as, (g) t = 0.08 as, (h) t = 0.09 as, (i) t = 0.12 as, (j) t = 0.145 as, (k) t = 0.17 as, and (l) t = 0.2 as.

$Figure 7 Final orbital radius dependence on 1 ∕ Z 1/Z . The solid line shows the calculation value by the FDTD method based on the eight-element Dirac equation, and the broken line denotes r final = r B ∕ Z {r}_{{\rm{final}}}={r}_{B}/Z obtained by the analytical solution of the Schrödinger equation.$

Figure 7

Final orbital radius dependence on 1 ∕ Z . The solid line shows the calculation value by the FDTD method based on the eight-element Dirac equation, and the broken line denotes r final = r B ∕ Z obtained by the analytical solution of the Schrödinger equation.

5 Conclusion

The transient analysis of the Dirac fields has been successfully implemented using the FDTD method based on the eight-element Dirac equation, which includes dual four-element wave functions and five spatially symmetric 8 × 8 matrices. We found that a stable three-dimensional Dirac field wave packet with low velocity can be created without explicitly considering Zitterbewegung, which is difficult in one-dimensional analysis. Furthermore, we achieved the first-ever simulation of the formation process of atomic orbitals starting from the initial state of a free electron and an atomic nucleus, using the eight-element Dirac equation without any physical approximations. Initially, a small unstable orbital appears, which rapidly grows and results in a large stable orbital with a radius equal to the Bohr radius divided by the atomic number, as given by the solution of the Schrödinger equation. The FDTD calculation based on the conventional four-element Dirac equation cannot produce such reasonable orbitals because of the spatial asymmetry of the 4 × 4 Dirac matrices. This method has the potential to be used for transient analyses of not only atomic or molecular orbitals but also interactions among elementary particles.

Funding information: The author states no funding involved.
Author contributions: The author has accepted responsibility for the entire content of this manuscript and approved its submission.
Conflict of interest: The author states no conflict of interest.

References

[1] Ohno K, Mutoh H, Harada Y. Application of Penning ionization electron spectroscopy to the study of chemical reactions on the solid surface; Photooxidation of naphthacene and rubrene. Surface Sci. 1982;115(3):L128–32. 10.1016/0039-6028(82)90373-9Search in Google Scholar

[2] Ohno K, Fujisawa S, Mutoh H, Harada Y. Application of Penning ionization electron spectroscopy to stereochemistry. Steric shielding effect of methyl groups on Penning ionization in substituted anilines. J Phys Chem. 1982;86(4):440–1. 10.1021/j100393a004Search in Google Scholar

[3] Ohno K, Mutoh H, Harada Y. Study of electron distributions of molecular orbitals by Penning ionization electron spectroscopy. J Am Chem Soc. 1983;105(14):4555–61. 10.1021/ja00352a009Search in Google Scholar

[4] Harada Y, Ohno K, Mutoh H. Penning ionization electron spectroscopy of CO and Fe (CO) 5. Study of electronic structure of Fe (CO) 5 from electron distribution of individual molecular orbitals. J Chem Phys. 1983;79:3251. 10.1063/1.446218Search in Google Scholar

[5] Veszprémi T, Harada Y, Ohno K, Mutoh H. Photoelectron and Penning ionization electron spectroscopic investigation of trimethylsilyl- and t-butyl-thiophenes. J Organometallic Chem. 1983;252(2):121–5. 10.1016/0022-328X(83)80074-6Search in Google Scholar

[6] Veszprémi T, Harada Y, Ohno K, Mutoh H. Photoelectron and Penning ionization electron spectroscopic investigation of trimethylphenylsilane. J Organometallic Chem. 1983;244(2):115–8. 10.1016/S0022-328X(00)98590-5Search in Google Scholar

[7] Veszprémi T, Harada Y, Ohno K, Mutoh H. Photoelectron and Penning electron spectroscopic investigation of phenylhalosilanes. J Organometallic Chem. 1984;266(1):9–16. 10.1016/0022-328X(84)80104-7Search in Google Scholar

[8] Veszprémi T, Harada Y, Ohno K, Mutoh H. Photoelectron and Penning ionization electron spectroscopic investigation of some silazanes. J Organometallic Chem. 1985;280(1):39–43. 10.1016/0022-328X(85)87060-1Search in Google Scholar

[9] Mutoh H, Masuda S. Spatial distribution of valence electrons in metallocenes studied by Penning ionization electron spectroscopy. J Chem Soc. Dalton Trans. 2002;9:1875–81. 10.1039/b111486gSearch in Google Scholar

[10] Yamazaki M, Horio T, Kishimoto N, Ohno K. Determination of outer molecular orbitals by collisional ionization experiments and comparison with Hartree-Fock, Kohn-Sham, and Dyson orbitals. Phys Rev A. 2007;75:032721. 10.1103/PhysRevA.75.032721Search in Google Scholar

[11] Eland JHD. Photoelectron spectroscopy. 2nd ed. London: Butterworth; 1984. 10.1016/B978-0-408-71057-2.50005-3Search in Google Scholar

[12] Paul PM, Toma ES, Breger P, Mullot G, Agostini P. Observation of a train of attosecond pulses from high harmonic generation. Science. 2001;292:1689. 10.1126/science.1059413Search in Google Scholar PubMed

[13] Itatani J, Levesque J, Zeidler D, Niikura H, Pepin H, Kieffer JC, et al. Tomographic imaging of molecular orbitals. Nature. 2004;432:867–71. 10.1038/nature03183Search in Google Scholar PubMed

[14] Yagishita A, Hosaka K, Adachi J. Photoelectron angular distributions from fixed-in-space molecules. J Electron Spectrosc Relat Phenom. 2005;142:295–312. 10.1016/j.elspec.2004.09.005Search in Google Scholar

[15] Niikura H, Wörner HJ, Villeneuve DM, Corkum PB. Probing the spatial structure of a molecular attosecond electron wave packet using shaped recollision trajectories. Phys Rev Lett. 2011;107:093004. 10.1103/PhysRevLett.107.093004Search in Google Scholar PubMed

[16] Weigold E, Hood ST, Teubner PJO. Energy and angular correlations of the scattered and ejected electrons in the electron-impact ionization of argon. Phys Rev Lett. 1973;30:475–8. 10.1103/PhysRevLett.30.475Search in Google Scholar

[17] McCarthy IE, Weigold E. (e, 2e) spectroscopy. Phys Rep. 1976;27:275–371. 10.1016/0370-1573(76)90005-3Search in Google Scholar

[18] Brion CE. Looking at orbitals in the laboratory: the experimental investigation of molecular wavefunctions and binding energies by electron momentum spectroscopy. Int J Quant Chem. 1986;29:1397–428. 10.1002/qua.560290534Search in Google Scholar

[19] McCarthy IE, Weigold E. Electron momentum spectroscopy of atoms and molecules. Rep Prog Phys. 1991;54:789–879. 10.1088/0034-4885/54/6/001Search in Google Scholar

[20] Coplan MA, Moore JH, Doering JP. (e, 2e) spectroscopy. Rev Mod Phys. 1994;66:985–1014. 10.1103/RevModPhys.66.985Search in Google Scholar

[21] Yamazaki M, Oishi K, Nakazawa H, Zhu C, Takahashi M. Molecular orbital imaging of the acetone S 2 excited state using time-resolved (e, 2e) electron momentum spectroscopy. Phys Rev Lett. 2015;114:103005. 10.1364/UP.2016.UTh4A.4Search in Google Scholar

[22] Mulliken RS. The assignment of quantum numbers for electrons in molecules II. Phys Rev. 1928;32:761. 10.1103/PhysRev.32.761Search in Google Scholar

[23] Hückel E. Quantum theory of double linkings. Z Phys. 1930;60:423. 10.1007/BF01341254Search in Google Scholar

[24] Hartree DR. The wave mechanics of an atom with a non-Coulomb central field. Part I. Theory and methods. Proc Camb Phil Soc. 1928;24:89. 10.1017/S0305004100011919Search in Google Scholar

[25] Fock V. Näherungsmethode zur Lösung des quantenmechanischen Mehrkörperproblems. Z Phys. 1930;61:126. 10.1007/BF01340294Search in Google Scholar

[26] Slater JC. Atomic shielding constants. Phys Rev. 1930;36:57. 10.1103/PhysRev.36.57Search in Google Scholar

[27] Hoffman R. An extended Hückel theory. I. Hydrocarbons. J Chem Phys. 1963;39:1397. 10.1063/1.1734456Search in Google Scholar

[28] Pople JA, Segal GA. Approximate self-consistent molecular orbital theory. III. CNDO results for AB2 and AB3 systems. J Chem Phys. 1966;44:3289. 10.1063/1.1727227Search in Google Scholar

[29] Roothaan CCJ. Self-consistent field theory for open shells of electronic systems. Rev Mod Phys. 1960;32:179. 10.1103/RevModPhys.32.179Search in Google Scholar

[30] Hay PJ. Gaussian basis sets for molecular calculations. The representation of 3d orbitals in transition‐metal atoms. J Chem Phys. 1977;66:4377. 10.1063/1.433731Search in Google Scholar

[31] Yee KS. Numerical solution of initial boundary value problems involving Maxwell’s equations in isotropic media. IEEE Trans Antennas Propagat. 1966;14:302. 10.1109/TAP.1966.1138693Search in Google Scholar

[32] Taflove A, Hagness SC. Computational electrodynamics: the finite-difference time-domain method. third edition. Norwood, MA: Artech House; 1995. Search in Google Scholar

[33] Simicevic N. Three-dimensional finite difference-time domain solution of Dirac equation. 2008. arXiv:physics.comp-ph:0812.1807v1. Search in Google Scholar

[34] Simicevic N. Finite difference-time domain solution of Dirac equation and the Klein paradox. 2009. arXiv:quant-ph:0901.3765v1. Search in Google Scholar

[35] Mutoh H. Extension of Maxwell’s equations for charge creation-annihilation and its applications. In: Proceedings on 2018 Progress in Electromagnetics Research Symposium. 2018. p. 2537–45. 10.23919/PIERS.2018.8597609Search in Google Scholar

[36] Mutoh H. Extended Maxwell’s Diamond equations to unify electromagnetism, weak gravitation, and classical and quantum mechanics and their applications for semiconductor devices. In: Proceedings on 2019 Progress in Electromagnetics Research Symposium. 2019. p. 1749–56. 10.1109/PIERS-Spring46901.2019.9017460Search in Google Scholar

[37] Dirac PAM. The principles of quantum mechanics, 4th edition. Oxford: Oxford University Press; 1958. Search in Google Scholar

[38] Peskin ME, Schroeder DV. An introduction to quantum field theory. Boulder, Colorad: Westview; 1995. Search in Google Scholar

[39] Mandl F, Shaw G. Quantum field theory. 2nd edition. Chichester, West Sussex, UK: Wiley; 2010. Search in Google Scholar

[40] Sakurai JJ. Modern quantum mechanics. 2nd edition. San Francisco, CA: Addison-Wesley; 2011. Search in Google Scholar

[41] Schrödinger E. Zur Quantendynamik des Elektrons. Sitzungsberichte der Preußischen Akademie der Wissenschaften. Physikalisch-mathematische Klasse. In: Proceedings of the Prussian Academy of Sciences, Physical-mathematical class. Berlin: Deutsche Akademie der Wissenschaften (German Academy of Sciences); 1931. p. 63–72. Search in Google Scholar

[42] Sakurai JJ. Advanced quantum mechanics. San Francisco, CA: Addison-Wesley; 1967. Search in Google Scholar

[43] Schiff LI. Quantum mechanics. 3rd edition. New York: McGraw-Hill; 1970. Search in Google Scholar

Received: 2023-03-11

Revised: 2023-05-31

Accepted: 2023-06-09

Published Online: 2023-07-13

This work is licensed under the Creative Commons Attribution 4.0 International License.

Articles in the same Issue

https://doi.org/10.1515/phys-2022-0262

Keywords for this article

Dirac equation; Maxwell’s equations; FDTD method; Zitterbewegung; atomic orbital

Creative Commons

BY 4.0