# Single-cone real-space finite difference schemes

for the Dirac von Neumann equation

###### Abstract

Two finite difference schemes for the numerical treatment of the von Neumann equation for the (2+1)D Dirac Hamiltonian are presented. Both utilize a single-cone staggered space-time grid which ensures a single-cone energy dispersion to formulate a numerical treatment of the mixed–state dynamics within the von Neumann equation. The first scheme executes the time-derivative according to the product rule for ”bra” and ”ket” indices of the density operator. It therefore directly inherits all the favorable properties of the difference scheme for the pure–state Dirac equation and conserves positivity. The second scheme proposed here performs the time-derivative in one sweep. This direct scheme is investigated regarding stability and convergence. Both schemes are tested numerically for elementary simulations using parameters which pertain to topological insulator surface states. Application of the schemes to a Dirac Lindblad equation and real–space–time Green’s function formulations are discussed.

###### keywords:

Dirac equation, vonNeumann equation, finite-difference scheme, staggered grid, fermion doubling, topological insulator## 1 Introduction - preliminaries and definitions

### 1.1 The Dirac equation and numerical schemes

In context with graphene and, more recently, topological insulator surface states, the Dirac equation has received renewed interest within the physics community. Introduced by P. A. M. Dirac in 1928, the (3+1)D Dirac equation traditionally has been a model in high-energy physics to describe relativistic spin-1/2 particles and, as a field theory, is a key ingredient to the standard model of elementary particle physics.dirac ; sakurai ; ryder ; itzykson ; srednicki In condensed matter and atomic physics its importance lies in the non-relativistc limit providing the spin-orbit interaction, which is instrumental to an understanding of atomic spectra and represents the foundation for the field of spintronics.spintr In electronic structure calculations the full Dirac equation has been used to describe inner shell electrons.wien1 ; wien2 While the (1+1)D and (2+1)D Dirac equation, allowing a two-dimensional representation of the Clifford algebra, in the early days of quantum physics was used for simplicity in formal model studies, condensed matter physics has recently identified physical realizations of the (2+1)D Dirac equation in form of low energy electronic excitations in graphene and topological insulators (TI).neto ; qi ; hsieh ; hasan ; xia ; analytis

With the interest in graphene and TIs as components in electronic and spintronic devices efficient schemes for the numerical solution of the (2+1)D Dirac equation have become desirable. Numerical approaches for solving the Dirac equation have taken several approaches. For investigations of relativistic electrons in atomic physics, a (2+1)D FFT-split-operator code was usedpiazza ; mocken . In such an approach, fast Fourier transformation between the space and momentum representation is used to compute the Dirac propagator. The computational effort scales like where is the number of grid-points. An efficient code using operator splitting in real space was introduced for the (3+1)D case.gourdeau It leads to the highly efficient operations count of . Alternatively, approaches using finite-difference approximations to the first-order space-time derivatives have been used and have foremostly been applied in lattice quantum chromo dynamics (QCD).QCD Traditionally they were plagued by the fermion doubling problem associated with centered difference approximations to first-order space- and time-derivatives.Stacey We note that such a 2d finite-difference scheme for the Dirac–Poisson system was presented recently.brinkman

Recently we have identified a class of staggered grids on which fermion doubling can be avoided in arbitrary space dimensions, with explicit schemes given for the (1+1)D, (2+1)D, and (3+1)D case.hammer1d ; hammer3d This class of schemes scales linearly in the number of grid points and allows a gauge-invariant formulation of electromagnetic fields on the lattice. Moreover, this grid can be used in any space-time formulation involving the Dirac operator, such as the density matrix and Green’s function approach, while yielding a single Dirac cone.

In this work we develop two numerical schemes for the mixed-state time evolution under the Dirac Hamiltonian in (2+1)D , which we term the Dirac-von-Neumann equation

(1) |

It describes the coherent mixed state dynamics of Dirac fermions which is useful for simulations of quantum transport on TI surfaces near the Dirac point. Both schemes to be presented in the following sections utilize the staggered grid of Ref. hammer3d to define centered differences over a single lattice spacing, thereby eliminating the very source for fermion doubling from the start. The first scheme (”bra-ket”) scheme, introduced in Sect. 2, treats the time derivative of the density operator within the product rule as indicated by the commutator on the r.h.s. of Eq. 1: apply from the left and therewith propagate the ket in time, apply from the right and therewith propagate the bra in time, and then form the difference. The second (”direct”) scheme treats the time-derivative of in one step and is introduced in Sect. 3, with further details given in the Appendix. Numerical examples, one for each scheme, are given in Sect. 4. Finally, we briefly discuss the extension of this approach to the Lindblad equation and space-time formulations of the Green’s function method for the Dirac Hamiltonian in Sect. 5 and close with Summary and Outlook.

### 1.2 The continuum formulation of the von Neumann equation for the (2+1)D Dirac Hamiltonian

We consider the effective model Hamiltonian which accounts for the energy spectrum of TI surface states near the Dirac pointqi ^{1}^{1}1Note that the TI representation differs from the standard representation used in Ref. hammer3d , however, shifting and converts the latter to the former.

(2) |

denote the Pauli matrices.
Using the abbreviation^{2}^{2}2For now, the electromagnetic vector potential is set equal to zero.

and (omitting space–time arguments)

this Hamilton operator takes the form of a matrix

Inserted into the von Neumann equation (1) one obtains a set of first-order partial differential equation in space and time for the density operator elements , conveniently written as the matrix identity

Note that, as in the original von Neumann equation, Hamilton matrix operator elements to the right of the density operator matrix elements act to the left (and vice versa). The two-component nature of the spin- Dirac fermion suggests this form. Choosing a continuous space representation for the orbital degrees of freedom one arrives at

(3) |

(4) |

Here we omit the time variable , common to all terms, for brevity, and use the real-space versions of the abbreviations defined above, i.e.,

and (including arguments)

Note that in the real-space representation we have

### 1.3 Placement onto a space-time grid

The task at hand now is to develop a space-time finite-difference approximation to Eq. 4 which avoids fermion doubling. This is achieved by using the staggered grid introduced in an earlier paper to accommodate the real-space density matrix in which upper and lower spinor component(s) are placed on two adjacent time-sheets.hammer3d The proper implementation is found by inspection of the density operator for a pure state spinor. In (2+1)D one has a pure-state ket and a ket-bra projector for the density operator

(5) |

This shows that the Pauli indices of and respectively, take the face-centered rectangular -grid () and -grid () position. For a single-cone representation the latter arehammer3d

Here the time index is labeled . Note that, for given time, and are placed onto adjacent time sheets, respectively, and . Each time sheet contains a rectangular face-centered spatial grid, symmetrically staggered relative to the ones on the two adjacent time sheets, such that symmetric difference quotients replace the respective partial derivatives on the grid. The grid spacings are denoted by , , and . We also introduce the associated grids

These grids are obtained from their non-barred counterparts by a forward time-shift by and thus share the spatial positions with the former. These will be the grids on which we place mass, scalar potential, and all partial derivatives. Common to all four grids is their symmetry in space-time, as required in a covariant theory. Since most dynamic problems are formulated as an initial–value problem in time, however, it is useful to view each of the two grids as a set of time sheets. Each time sheet in turn consists of a rectangular-, for , or cubic-face-centered spatial grid, for , which in units is characterized by the set of relative positions .

The matrix elements (and the vector potential to be introduced later) live on the grids

Note that leads to the need for two adjacent time sheets (e.g., and ) for the representation of all matrix elements of on the lattice: while and each are placed on a single time sheet, and require the use of two adjacent time sheets. This leads one to the definition of the trace of using two adjacent time sheets of the grid:

###### Definition 1.

The trace of at a given time () on the grid,

(6) | |||||

is defined as the sum of the trace of with respect to the lattice sites and the trace of with respect to the lattice sites . The summation runs over all lattice sites of the respective time sheet, .

We extend this definition to lattice sums for non-diagonal density matrix elements

(7) |

with constant relative indices and .

###### Definition 2.

Centered spatial difference operators are defined on the barred grids as follows

and

Furthermore, we will use the time-average operation

Formally these operators are defined on the barred grids , however, the execution involves density matrix elements on the unbarred grids. Note that when applied to the density matrix, these operators may act on the first (bra) or second set of indices (ket).

Partial derivatives, scalar potential and mass terms, respectively, are implementing into Eq. (3) by the substitutions

Within a sum over all lattice sites a spatial difference operation (, k=x,y) performed on the second set of indices of the density matrix is equivalent to minus the derivative () performed on the first set .

###### Corollary 1.

###### Proof.

We give the proof for

(9) | |||||

We have used that and that, under the trace and are constant. Note that a combined shift (carried out in the third step) of and is required to stay on the proper sub-grid . This mixes corner and face-center positions of a given time sheet. However, since the sum is to be taken over all grid points this ”integration by parts” rule holds under the trace. Note also how the barred grids used for the definition of center position of the spatial derivative are related to the index of the density matrix elements. Other cases for , , and can be shown in the same fashion. ∎

## 2 The bra-ket scheme

The bra-ket time-propagation scheme is obtained directly by adapting the scheme for the pure-state Dirac equation, applied to the bra- and the ket-side of the density operator. This corresponds to the interpretation

(10) |

###### Definition 3.

Within the representation Eq. (2) the pure state time–evolution of Ref. hammer3d from initial time to final time may be written as

(11) |

(12) |

with the coefficients

living on integer time sheet of grid , and

living on half-integer time sheet of grid . We have used which holds for real-valued mass and scalar potential .

Using the short-hand notation , the progression by one full time step is executed as follows (initial () and final () time is indicated by the respective superscript):

This set of operations is followed by

Executed in this order the scheme is explicit and follows exactly the single-cone time-propagation of Ref. hammer3d applied to both sets of indices of the density matrix, with unprimed and primed operators acting on, respectively, bra and ket. As such all properties of the scheme for the pure-state propagation, such as stability, convergence, and spectral properties, are carried over to the scheme for the density matrix.

###### Definition 4.

With the abbreviations defined above, the single time-step progression under the bra–ket scheme is defined as

(13) |

(14) |

(15) |

(16) |

Def. 4 extends the single-cone pure-state time-propagation scheme in Def. 3 to the one for mixed-states. For pure states, the latter reduces to the former.

The following lemmas and propositions pertaining to the bra-ket scheme are a direct consequence of the properties of the pure state scheme Ref. hammer3d and the observation that the density operator at initial time, without loss of generality, can be written as as sum of pure-state contributions of the form Eq. 5.

###### Lemma 2.

The bra-ket scheme conserves positivity and Hermiticity of the density matrix.

###### Proof.

Both properties are a direct result of the bra–ket scheme: Let, at initial time the density matrix be positive definite, i.e.,

for an arbitrary element of the Hilbert space. Abbreviating the one-step time evolution of a pure state Eqs. (11) and (12) by ,

the bra-ket scheme propagates

since positivity holds for . Hence under the bra-ket scheme positivity is preserved for arbitrary time step. Preservation of Hermiticity is seen from . ∎

Note that the standard norm is not conserved within the pure-state scheme since is not unitary.hammer3d Hence, although the time-evolution superoperator can be cast in Kraus form with a single Kraus operator (living on two adjacent time-sheets), we do not have .choi ; kraus

###### Lemma 3.

###### Proof.

The proof of this lemma as well as the following two propositions rest on the fact that the density matrix at initial time may be cast into diagonal form and then, under the scheme, is placed onto the grid according to

(18) |

where are normalized to one, real, and . Since each term in this sum is propagated individually within this scheme, all the stability properties of the pure-state dynamics apply. The conserved functional for the individual pure-state contribution, , and the definition of the trace on the grid in Eq. (6) as a sum over all lattice sites of the respective time-sheet render the conserved functional for the density operator . ∎

###### Proposition 4.

###### Proof.

Using the decomposition Eq. (18) stability can be shown term by term in the sum. ∎

Comment: Conservation of on the grid corresponds to the conservation of the trace of in the continuum limit. At the same time the conservation of the former within this scheme implies non-conservation of the trace when the latter is defined on the grid according to Eq. (6). This non-conservation is of order and thus can be adjusted by this ratio.

From the analysis of the pure-state scheme we know that Prop. 4 is too restrictive. In fact, the stability condition holds for constant mass and potential. For the special case this less restrictive stability condition reads . This stability condition also holds for an arbitrary space- and time-dependent mass and potential.

###### Definition 5.

We define the functionals^{3}^{3}3For compactness of notation the diagonal time-index is placed onto as a superscript.

Again, the summation runs over all lattice sites of time sheet , i.e., , with and and .