QuantumFDTD - A computational framework for the relativistic Schrödinger equation

. We extend the publicly available quantumfdtd code. It was originally intended for solving the time-independent three-dimensional Schrödinger equation via the finite-di ff erence time-domain (FDTD) method and for extracting the ground, first, and second excited states. We (a) include the case of the relativistic Schrödinger equation and (b) add two optimized FFT-based kinetic energy terms for the non-relativistic case. All the three new kinetic terms are computed using Fast Fourier Transform (FFT). We release the resulting code as version 3 of quantumfdtd. Finally, the code now supports arbitrary external file-based potentials and the option to project out distinct parity eigenstates from the solutions. Our goal is quark models used for phenomenological descriptions of QCD bound states, described by the three-dimensional Schrödinger equation. However, we target any field where solving either the non-relativistic or the relativistic three-dimensional Schrödinger equation is required.


Introduction
In our work [1] we extended the quantumfdtd software package [2][3][4][5].Besides containing a detailed study about numerical stability of the new quantumfdtd code, Ref. [1] is a reference manual for the code.
The legacy code is a parallelized solver for the three-dimensional non-relativistic Schrödinger equation with a selection of analytical real and complex potentials.The solver is based on the finite-difference time-domain (FDTD) method.Via an iterative method (inherited from the legacy quantumfdtd code [2][3][4][5]), the code is able to extract the ground, first, and second excited states.The excited states are determined from multiple nearby temporal snapshots of the iterative procedure that are consecutively subjected to orthogonal projection.The analytical hard-coded potentials could be easily extended by writing the appropriate code inside the potentials.cppsource file.The chosen FDTD method implies that the boundary conditions in the edge of the 3D lattice cube are Dirichlet ones, with the wavefunction, Ψ, vanishing on the boundary, Ψ(boundary) = 0.The only library that is required to run the code is MPI, for the parallelization.
In our extension of the quantumfdtd package three new kinetic terms are included.These are based on the Fast Fourier Transform (FFT), that requires periodic boundary conditions in the edge of the 3D lattice cube.For the non-relativistic and FFT-based solvers, most of the differences from the FDTD results come from this difference in the boundary conditions, as was extensively discussed on our work [1].Part of such a study will be also reviewed here.
We consider two non-relativistic FFT-based kinetic terms, the difference between them being the usage of the Symanzik effective field theory [6,7] for reducing the discretization errors.The third FFT-based kinetic term is the relativistic one, intended for solving the 3D relativistic Schrödinger equation.For the computation of the FFT we require our code to be linked to the FFTW_MPI version 3 library. 1he target of this software is quark models used for phenomenological description of QCD bound states.Most of these models are described by the three-dimensional Schrödinger equation with different potentials.In particular, these models have successfully described below-threshold charmonium production and bottomonium spectra and have helped to establish confidence in QCD as the first-principles description of hadronic matter [9][10][11].Nonrelativistic effective field theory methods have been used for a first-principles approach to potential-based non-relativistic QCD (pNRQCD) [12,13].In order to describe quarkonium evolution in the quark-gluon plasma, these potential models have been extended to finite temperature [14][15][16] and non-equilibrium [3,4,[17][18][19].In the latter case, the potentials are no longer real-valued or spherically symmetric.In the full non-equilibrium case, a full threedimensional solver, like the one described in this work, is necessary.
Finally, for integration into a lattice field theory workflow (and for the general convenience of the user independently of the research field) we have extended the code by allowing for loading arbitrary complex potentials from ASCII files with numerical values.We now require the GNU Scientific Library (GSL), 2 with the CBLAS link. 3 We also include several post-processing scripts.The most powerful ones aim to aid the extraction of excited states by exploiting their symmetries.For running these scripts, Python 3 is required. 4

New kinetic terms and iterative procedure
The Schrödinger Hamiltonian can be split into a kinetic piece H K and a potential V( r).The non-relativistic kinetic term H nr K , already implemented (in its FDTD discretized version) in the legacy quantumfdtd code, and the relativistic one H rel K are given, respectively, by where p = (p 1 , p 2 , p 3 ) is the spatial three momentum.We provide the following four implemented kinetic terms: the original one H (0) K , based on the FDTD discretization, and three new kinetic terms, H (1)  K , H (2)  K (both non-relativistic) and H (3) K (relativistic), Here, êl denotes the (signed) vector in the l direction (l = ±1, ±2, ±3) of one unit of the lattice spacing A and k l ≡ Ap l ∈ {−π + 2π/N, . . ., 0, 2π/N, 4π/N, . . ., π} are the lattice momentum coordinates.The FFT and IFFT operators denote the Fast Fourier Transform and Inverse Fast Fourier Transform, that are computed via the FFTW_MPI version 3 library [8], as stated in the previous section.Finally, Ψ( r) ≡ Ψ[(Ax 1 , Ax 2 , Ax 3 )] for convenience reasons.
For solving the Schrödinger equation, an iterative method inherited from the legacy quantumfdtd code [2][3][4][5] is used.Since this method has not been modified from the legacy code, the reader is referred to [2] for its description.In brief, it is based on a Wick rotation of the temporal variable i t → τ, so that the Schrödinger equation turns into Its general solution allows for the extraction of both the wavefunctions Ψ n and energies E n of the ground n = 0 and excited n ≥ 1 states, once the intermediate wavefunctions Ψ( r, τ) have been obtained.In order to extract such a solution Ψ( r, τ), we discretize the temporal derivative in Eq. 7.For details see Ref. [1].

Results
Ref. [1] includes a comprehensive study of all the relativistic and non-relativistic kinetic terms, H (i) K (i = 0, . . ., 3).Thus, we summarize here only some of the most important results.For the non-relativistic kinetic terms, we use Coulomb potential (for which the analytical solution is known), a lattice spacing of A = 0.12 fm, and a reduced mass of M = 0.3 GeV.Our goal is demonstrate the capabilities of the code by reproducing the analytically known wavefunctions and eigenenergies.Of course, both discretization errors and finite volume effects will lead to deviations.Finite volume effects, in particular, are expected to produce a discrepancy between H (0 K (based on the FDTD method) and H (1)  K , H (2)  K (based on FFT), due to the different boundary conditions: Dirichlet ones for H (0) K and periodic ones for H (1)  K and H (2) K .K and H (1)  K are almost indistinguishable [1].Right: relativistic kinetic term H (3)  K with the modified harmonic oscillator potential (Ref.[24]) momentum space wavefunctions with the relativistic kinetic term H (3)  K .Here, y(p) = √ 4πpΨ(p) in order to compare to the theoretical predictions following [24].
In Fig. 1 we can see the analytical ground, first, and second excited states Ψ th 1s , Ψ th 2s , Ψ th 2p of the Coulomb potential.In the same plots, we show parity projections of the wavefunctions computed by quantumfdtd with the H (0) K kinetic term (FDTD based, left) and H (2)  K (FFT based, right).Note that, because of the structure of the exact (analytical) solution, P + Ψ ns ( r) = Ψ ns ( r), P − p k Ψ ns ( r) = 0, P + Ψ n p ( r) = 0 and P p k Ψ n p ( r) = Ψ n p ( r), where n = 1, 2, 3, . . .and n = 2, 3, . . . .These parity projectors, that have been implemented [1] in the Python 3 scripts included with quantumfdtd, allow for a more precise extraction of the ground, first, and second excited states.Indeed, the results shown in Fig. 1 are very close to the analytically known solution.
In the left panel of Fig. 1 (FDTD based, Dirichlet boundary conditions) the wavefunctions decay too fast near the edge (expected with Dirichlet boundary conditions).On the right (periodic boundary conditions) the wavefunctions show artifacts near the boundary.These artifacts are due to the breaking of spherical symmetry.This is also expected for periodic boundary conditions with a non-spherically symmetric, cubic, lattice.Finally, on both plots there are artifacts near the origin r/A → 0, where discretization effects are more prominent (and the Coulomb potential has a pole).The wavefunction of H (1)  K is almost indistinguishable from that of H (2)  K (see Ref. [1]).On the left of Fig. 2 we show the non-relativistic wavefunction obtained using a finer lattice N = 128.The wavefunctions of H (0) K , H (1)  K , and H (2)  K are almost indistinguishable and reproduce almost exactly the theoretical solutions Ψ th 1s , Ψ th 2s , Ψ th 2p .The only difference is a small artifact very close to the origin r/A = 0, where discretization effects (finite A) are most prominent.
Finally, as a means to check the relativistic kinetic term H (3) K , we included the right plot of Fig. 2 to check a specially modified version of the relativistic harmonic oscillator potential for which an approximate analytical solution exists [24].In momentum space (chosen by Ref. [24]) our solver is able to reproduce the analytical result of [24] for the ground, 1s, and first excited, 2s, states to sub percent levels.Furthermore, our code is not restricted to s-wave states as is the analytical procedure of Ref. [24].For instance, we can compute the 2p state, as shown in [1].Our goal here is a direct comparison between our code [1] and the analytical procedure [24] where the latter one is applicable, so that we can cross-check the relativistic kinetic term H (3)  K .As shown in Fig. 2, our code behaves correctly at sub percent level in this application of the relativistic solver.

Conclusions
In this work we have extended the legacy quantumfdtd code to support the relativistic Schrödinger equation and we have added two new solvers of the non-relativistic one.The new solvers are based on the Fast Fourier Transform (FFT) and have periodic boundary conditions instead of Dirichlet ones.Furthermore, we have performed a comprehensive study of the stability of the solvers and compared them against analytical (exact) solutions of the nonrelativistic Coulomb potential and a relativistic particle in a harmonic oscillator potential.Furthermore, we have extended the code to increase its usability by the lattice field theory community and by the larger scientific community interested in solving the 3D relativistic and non-relativistic Schrödinger equation.
There are a number of possible extensions to the present version of quantumfdtd that are beyond the scope of this work.From easiest to most difficult: variations of the coupling strength of the hard-coded potentials; new projection operators (for instance, d-wave states); implementing a suitable binary format for the wavefunctions; an automatized algorithm to look for saddle points of E i (τ) (i = 0, 1, 2), that is, automatizing the detection of the stabilization of the wavefunctions and eigenenergies; an implicit Crank-Nicolson method for the solver (difficult for the relativistic kinetic term, which is non-linear); and permitting the (anti)-symmetric solution of a two-particle system in a background potential with an interaction potential.