\section{Introduction}
\label{intro}

The three-dimensional Schr{\"o}dinger equation has a closed-form analytical solution for only a very small class of systems (i.e. a free particle in a box, harmonic oscillator, or some central potentials).
In all other cases, one has to rely on approximations, perturbation theory, or numerical solutions.
In particular, most quark models used for phenomenological descriptions of QCD bound states are described by the three-dimensional Schr{\"o}dinger equation with a wide variety of potentials.
Potential-based quark models have enjoyed a long history of success in describing below-threshold charmonium production and bottomonium spectra and have helped to establish confidence in QCD as the first-principle description of hadronic matter, see e.g. Refs.\cite{Ebert:2002pp, Eichten:2007qx, Segovia:2016xqb}.
In recent years, non-relativistic effective field theory methods have allowed for a first-principles approach to potential-based non-relativistic QCD (pNRQCD)~\cite{Brambilla:1999xf,Brambilla:1999qa}.
Additionally, in order to describe quarkonium evolution in the quark-gluon plasma, such potential models have been extended to finite temperature \cite{Laine:2006ns,Brambilla:2008cx,Brambilla:2010xn} and non-equilibrium \cite{Dumitru:2007hy,Burnier:2009yu,Dumitru:2009fy,Dumitru:2009ni,Guo:2018vwy}, in which case, the potentials are no longer real-valued or spherically symmetric.
In the full non-equilibrium case, it is necessary to use a full three-dimensional Schr{\"o}dinger equation solver that allows for complex potentials

In the past, a parallelized solver called \texttt{quantumfdtd} has been made available for real-valued~\cite{Strickland:2009ft} and complex-valued \cite{Margotta:2011ta} potentials.
The two previous versions of \texttt{quantumfdtd} use local operators in a finite-difference time-domain (FDTD) scheme for solving the non-relativistic three-dimensional Schr{\"o}dinger~equation.
However, if light quarks are involved, such as in heavy-light mesons, e.g. in B mesons, these have to be described as relativistic degrees of freedom.
The light quark's relativistic dispersion results in a much larger average radius of the state, namely with the wave-function's falloff in the case of a linearly rising Cornell potential $V(r) = -A/r+C r$ like $\psi(r) \sim \exp(-\mu r)$ instead of $\psi(r) \sim \exp(-\mu r^{3/2})$ ($\mu$ being the reduced mass of the system)~\cite{Duncan:1994uq}.
While spin-orbit and spin-spin corrections can be accounted for in terms of angular-momentum-dependent contributions beyond a simple central potential, at a fundamental level, the change to relativistic dispersion requires an implementation of a relativistic kinetic term in the Schr{\"o}dinger equation.
For this reason we extend the existing  \texttt{quantumfdtd}  parallelized 3D Schr{\"o}dinger solver~\cite{Strickland:2009ft, Margotta:2011ta} to permit a relativistic kinetic term.

These relativistic quark model calculations have a purpose beyond the purely phenomenological calculations.
It has been clearly demonstrated in the past~\cite{Duncan:1994uq} that heavy-light meson wave-functions from relativistic quark models provide interpolating operators for heavy-light mesons with excellent excited-state suppression in first-principle lattice QCD calculations, see i.e. Fig.~\ref{fig:effMass} for a very recent application.
In this context, the question of heavy-light tetraquarks looms large.
A three-dimensional relativistic Schr{\"o}dinger equation solver may be able to deliver interpolating operators for first-principle lattice QCD studies that are able to enhance overlap factor for certain geometries\footnote{The present version of the \texttt{quantumfdtd} code is not particularly well suited for an extension to systems with two fermionic degrees of freedom interacting with each other in the background of another potential, which would be needed if the two light quarks were not represented as a light diquark.} (i.e. spherical symmetry in a light-diquark heavy-antidiquark system, or cylindric symmetry in a light-diquark heavy-antiquark heavy-antiquark system) of tetraquarks, while effectively suppressing the excited states.
Such an approach might make a resolution of the geometric structure of tetraquarks feasible and distinguish between a wide range of possible scenarios.

\begin{figure}[!ht]
\centering
\igg{effM_fin.pdf}
\caption{\label{fig:effMass}
Effective mass plot for a heavy-strange meson correlation function (static limit for the heavy quark with 1 step of 4D hypercubic smearing~\cite{Hasenfratz:2001hp},
and the highly improved staggered quark discretization~\cite{Follana:2006rc} for the strange quark) with different source and sink operators.
The source operators are either point-type (namely, random Z(2) noise sources with no momentum projection) or wall-type (namely,
corner wall sources, with momenta  $p_i = n\pi/A$).
The sinks are either point sinks, wall sinks or relativistic quark model inspired wave-functions obtained with the code discussed here.
Simulations were conducted in (2+1+1)-flavor QCD with physical sea and lattice spacing $a \approx 0.15$~fm
using MILC ensembles~\cite{Bazavov:2017lyh} and the publicly available MILC code~\cite{MILC}.
The dotted line shows the average of the data obtained using the wave-function sink and the gray band shows the respective standard deviation.
It is clearly visible how much excited state suppression can be gained by using wave-functions at the sink.}
\end{figure}

In this work, we extend the original \texttt{quantumfdtd} code~\cite{Strickland:2009ft, Margotta:2011ta}.
The new version, like the old version, allows one to compute the ground-state and low-lying excited state
three-dimensional (3D) wave-functions and their corresponding energies for both real- and complex-valued potentials.
We implement a new relativistic kinetic term and add two new non-relativistic ones, both of which are Fast Fourier Transform (FFT)
based and implemented using the paralleized \texttt{FFTW} library. In addition, the original code has been extended so that it supports arbitrary external potentials
and the option to project out distinct parity eigenstates from the solutions.  These options also work with the original implementation
of the finite-difference time-domain (FDTD) solver for the non-relativistic Schr{\"o}dinger equation.
Because of the original code design, we use a 3-D Cartesian, homogeneous and isotropic lattice with $N$ points
in each dimension and lattice spacing $A$ (in units of $\text{GeV}^{-1}$).
No adaptive engine has been implemented.  The solver depends on the MPI library to distribute the lattice between the MPI processes.
Finally, we note that by construction, the original FDTD solver for the non-relativistic Schr{\"o}dinger equation used Dirichlet boundary conditions.
For the new kinetic terms, based on the Fast Fourier Transform (FFT), periodic boundary conditions are used.
The potential can be centered either in the lattice volume at the point $(N/2,N/2,N/2)$ or in the origin at the point $(0,0,0)$.
We use the original iterative solver (see Ref.~\cite{Strickland:2009ft} and section~\ref{subsec_iter} where we review the method).

This paper is organized as follows:
In section~\ref{sec_rel_kin} we introduce the newly implemented kinetic terms (subsection~\ref{subsec_kin}),
the iterative procedure (subsection~\ref{subsec_iter}), and the parity projection scripts (subsection~\ref{subsec_parity_proy}).
In section~\ref{sec_examples} we show and compare results obtained using the different kinetic terms.
In section~\ref{sec_compile}, the program requirements are listed and basic instructions for compilation and running are provided.
In section~\ref{sec_params}, we show the basic configuration parameters of the program (subsection~\ref{subsec_basic_params})
and list the different initial conditions (subsection~\ref{subsec_icond}) that are implemented.
We furthermore list the hard-coded potentials in subsection~\ref{subsec_hpot},
and present a newly implemented option that allows the user to load potentials from external files (subsection~\ref{subsec_epot}).
The new post-processing scripts, which include parity projection scripts for wave-functions, are explained in section~\ref{subsec_post_scripts}.
We end with a conclusion in section~\ref{sec_conclus}.

\section{The new relativistic kinetic term}
\label{sec_rel_kin}

The Schr{\"o}dinger Hamiltonian
\begin{equation}
H = H_{K} + V(\vec{r}) \,,
\end{equation}
can be generally split into a kinetic piece, $H_{K}$, and into a potential $V(\vec{r})$.
On the one hand, the non-relativistic kinetic term is given by
\begin{equation}
\label{eq:HKnr}
H_{K}^{nr} = \sum\limits_{i=1,2,3} \frac{p_{i}^{2}}{2m} \,,
\end{equation}
where $\vec{p} = (p_{1},p_{2},p_{3})$ is the spatial three momentum. 
On the other hand, the relativistic term is given by 
\begin{equation}
\label{eq:HKrel}
H_{K}^{rel} = \sqrt{ m^2 + \sum\limits_{i=1,2,3} p_{i}^{2} } \,,
\end{equation}
if neither spin nor antiparticle degrees of freedom are given consideration. 


In the following subsections we describe the newly implemented (relativistic) kinetic terms.

\subsection{Kinetic terms}
\label{subsec_kin}

There are three newly implemented kinetic terms, $H_{K}^{(i)}$, ($i = 1,2,3$), on top of the one based on
the finite-difference time-domain (FDTD) method, $H_{K}^{(0)}$, with Dirichlet boundary conditions ($\Psi(\text{boundary}) \equiv 0$).
Its action on a wave-function $\Psi = \Psi(\vec{r})$, ($\vec{r} = (Ax_{1},Ax_{2},Ax_{3})$ being spatial coordinates and $x_{1},x_{2},x_{3} \in \{0,1,\ldots, N-1\}$) 
is given by

\begin{equation}
\label{eq:HK0}
H_{K}^{(0)} \Psi = - \frac{1}{2M} \sum_{l = \pm 1}^{\pm 3} \frac{1}{2A}\frac{\Psi(r + \hat{e}_l) - \Psi(\vec{r})}{A} \,,
\end{equation}
where $\hat{e}_{l}$ denotes the (signed) unit vector in the $l$ direction ($l = \pm 1,\pm 2, \pm 3$). 
We write $\Psi(\vec{r}) \equiv \Psi[(Ax_{1}, Ax_{2}, Ax_{3})]$ for convenience reasons.


In particular, the FDTD kinetic term $H_{K}^{(0)}$ is based on Taylor expansions on the lattice sites and assumes continuity and differentiability.
In momentum space this term corresponds to 
\begin{equation}
\label{eq:HK0_momentum_space_1}
H_{K}^{(0)} \Psi(\vec{p}) = \frac{1}{2M} \sum\limits_{l=1}^{3} \frac{4}{A^2}\sin^{2}\left(\frac{A p_{l}}{2}\right) \Psi(\vec{p}) \,,
\end{equation}
where$ A p_{l}$ is real-valued and satisfies $-\pi < A p_{l} \le \pi$.  These coefficients are not constrained to a finite set of momenta (due to the absence of periodic boundary conditions). 
In the following we write $k_l \equiv A p_l$ for convenience. 
Note that there are no fundamental limitations for a modified FDTD kinetic term with periodic boundary conditions,
although it would require a substantial addition of structure encoding the communication across this boundary in the present version of the code.

The newly implemented kinetic terms depend on the FFT and, hence, use periodic boundary conditions, and have the momenta constrained to the eigenmodes of the finite box, i.e.
\begin{equation}
\label{eq:eigenmodes}
k_{l} \equiv A p_{l} \in \{ -\pi+2\pi/N, \ldots, 0, 2\pi/N, 4\pi/N, \ldots \pi \} \,,
\end{equation}
for $l = 1,2,3$.
The action of the new kinetic terms on a wave-function $\Psi = \Psi(\vec{r})$ ($\vec{r} = (Ax_{1}, Ax_{2}, Ax_{3})$ being spatial coordinates)
is given in terms of the Fast Fourier Transform (FFT) and Inverse Fast Fourier Transform (IFFT) operators (realized through \texttt{FFTW}) by

\begin{align}
& H_{K}^{(1)} \Psi = \frac{1}{2A^{2}M N^{3}} \cdot \mathop{\text{IFFT}}\left[ \sum\limits_{l=1,2,3} \left(k_{l}\right)^{2} \cdot \mathop{\text{FFT}}[\Psi] \right] \,, \label{eq:HK1} \\
& H_{K}^{(2)} \Psi = \frac{1}{2A^{2}M N^{3}} \cdot \mathop{\text{IFFT}}\left[ 4\sum\limits_{l=1,2,3} \sin^{2}\!\left(\frac{k_{l}}{2}\right) \cdot \mathop{\text{FFT}}[\Psi] \right] \,, \label{eq:HK2} \\
& H_{K}^{(3)} \Psi = \frac{1}{A N^{3}} \cdot \mathop{\text{IFFT}}\left[ \sqrt{4\sum\limits_{l=1,2,3} \sin^{2}\!\left(\frac{k_{l}}{2}\right) + (AM)^{2}} \cdot \mathop{\text{FFT}}[\Psi] \right] \,,\label{eq:HK3}
\end{align}
where the FFT and IFFT operations are defined via 
\begin{align}
& \left(\mathop{\text{FFT}}[\psi]\right)(k_1,k_2,k_3) = \sum\limits_{x_1,x_2,x_3} \Psi \cdot \exp\left(
-\sum_{l=1}^3 k_{l}\cdot x_l \right) \,, \\
& \left(\mathop{\text{IFFT}}[\tilde{\Psi}]\right)(x_1,x_2,x_3) = \!\! \sum\limits_{k_{1},k_{2},k_{3}} \!\! \tilde{\Psi} \cdot \exp\left( 
\sum_{l=1}^3 k_{l}\cdot x_l \right) \,.
\end{align}
Due to the realization of the kinetic terms through the FFT, optimized communications across the inter-node boundaries, as done in the FDTD approach, are not required here.

The term $H_{K}^{(1)}$ is a discretized FFT implementation of the classical non-relativistic kinetic term $H_{K}^{nr}$, Eq.~\eqref{eq:HKnr}.
For this reason, the discretization errors of Eq.~\eqref{eq:HK0_momentum_space_1}, namely,
\begin{equation}
 \label{eq:HK0_momentum_space}
 \left[ H_{K}^{(0)} -H_{K}^{nr} \right] \Psi(\vec{p}) = 
 \frac{1}{2M}\sum\limits_{l=1}^3 \frac{2}{A^2}\sum\limits_{n=2}^\infty (-1)^n \frac{(A p_l)^{2n}}{(2n)!}  \Psi(\vec{p}) \,,
\end{equation}
are completely absent in $H_{K}^{(1)}$.
Note that although $H_{K}^{(1)}$ formally seems to avoid the breaking of the $O(3)$ symmetry to the cubic group, this symmetry breaking is implicit in the FFT.
The $2\pi$ periodicity and restriction to the interval of Eq.~\eqref{eq:eigenmodes} is automatically realized in the FFT implementation.
Using the Symanzik effective field theory~\cite{Symanzik:1983dc, Symanzik:1983gh} the discretization errors of $H_{K}^{(0)}$
could be systematically reduced by including more extended stencils with adjusted coefficients, e.g.
\begin{align}
\label{eq:HK0imp_momentum_space}
&\begin{aligned}
H_{K}^{(0), \text{imp}} \Psi(\vec{p}) &= \frac{1}{2M}\sum\limits_{l=1}^{3} \frac{4}{A^2} \left[ c_0\sin^2\left(\frac{A p_l}{2}\right)+ c_1\sin^2(A p_l) \right] \Psi(\vec{p}), \\
&= \frac{1}{2M}\sum\limits_{l=1}^{3} \frac{2}{A^2} \sum\limits_{n=1}^\infty \frac{c_0+2^{2n}c_1}{(2n)!} (Ap_l)^{2n} \Psi(\vec{p}) \,,
\end{aligned} \\
\label{eq:HK0imp_coordinate_space}
&\begin{aligned}
H_{K}^{(0), \text{imp}} \Psi(\vec{r}) &= \frac{1}{2M}\sum\limits_{l=\pm 1}^{\pm 3} \frac{1}{2A} \frac{c_0\left[ \Psi(\vec{r}+\hat{e}_l)  -\Psi(\vec{r}) \right] 
+ c_1\left[ \Psi(\vec{r}+2\hat{e}_l) -\Psi(\vec{r}) \right]}{A} \,.
\end{aligned}
\end{align}
The choice $c_0 = 4/3$, $c_1 = -1/12$ restores the leading term to its continuum coefficient and cancels the second term in the series.
Hence, there cannot be any local coordinate space representation of $H_{K}^{(1)}$, since infinitely many stencils up to infinite extent
would be required to completely remove the discretization errors.

The term $H_{K}^{(2)}$ is an FFT implementation of the classical non-relativistic kinetic term $H_{K}^{(0)}$ for periodic boundary conditions.
Due to the different boundary conditions there are differences between those terms that vanish in the infinite volume limit.
However, the discretization errors of $H_{K}^{(0)}$ and $H_{K}^{(2)}$ are formally equivalent as they have the same coefficients
in the Symanzik effective theory in the infinite volume limit.

The term $H_{K}^{(3)}$ is the newly implemented relativistic kinetic term, which realizes the most simple discretized version of Eq.~\eqref{eq:HKrel}.
It is inspired by a relativistic quark model (RQM)~\cite{Duncan:1994uq} and can be used by means of the \texttt{KINTERM} parameter (see Section~\ref{sec_params}).
There is no local or continuously differentiable coordinate space representation of the relativistic kinetic term, as is evident
from the momentum dependent term under the square root.
Using Symanzik effective theory~\cite{Symanzik:1983dc, Symanzik:1983gh} the discretization errors of $H_{K}^{(3)}$
could be systematically reduced by including more extended stencils with adjusted coefficients just as in the non-relativistic case.
It is also straightforward to construct an relativistic analog to $H_{K}^{(1)}$.
None of these two obvious extensions have been pursued in this version of the program.

Note that the FDTD kinetic term, $H_{K}^{(0)}$, and each of the FFT kinetic terms, $H_{K}^{(i)}$, $i = 1,2,3$,
are qualitatively different and, in practice, exhibit different numerical behavior.
On the one hand, comparing the results of $H_{K}^{(0)}$ and $H_{K}^{(2)}$ permits a quantitative study of
finite volume errors while using a single volume, due to the different boundary conditions.
On the other hand, comparing the results of $H_{K}^{(1)}$ and $H_{K}^{(2)}$ permits a quantitative study
of discretization errors while using a single lattice spacing, due to different coefficients in the Symanzik effective theory.
Nevertheless, the dependence of the low-energy spectrum on these differences is usually very small.
On the contrary, the relativistic kinetic term encodes qualitatively different physics already in the low-energy spectrum.

\subsection{The iterative procedure}
\label{subsec_iter}

The iterative procedure is unmodified w.r.t. \texttt{quantumfdtd v.2}.
That is, it is the very same procedure as described in section~3 of Ref.~\cite{Strickland:2009ft},
\begin{equation}
\label{iter_proced} \Psi(x_{1},x_{2},x_{3},\tau+\Delta\tau) = \mathcal{A} \Psi(x_{1},x_{2},x_{3},\tau) - \mathcal{B} \Delta\tau H_{K} \Psi(x_{1},x_{2},x_{3},\tau) \,.
\end{equation}
Equation~\ref{iter_proced} is the result of a Wick rotation of the Schr{\"o}dinger Equation, $it \to \tau$.
The evolution with $\tau$ is given by Eq.~(2.5) of Ref.~\cite{Strickland:2009ft},
\begin{equation}
\label{iter_decay} \Psi(x_{1},x_{2},x_{3},\tau) = \sum_{i=0}^{\infty} a_{n}\Psi_{i}(x_{1},x_{2},x_{3}) e^{-E_{i} \tau} \,,
\end{equation}
where $E_{i}$ is the $i$-th eigen-energy of the Hamiltonian and $\{a_0,a_1,\dots\}$ are the decomposition coefficients
of the initial guess $\Psi(x_{1},x_{2},x_{3},0)$ in the basis of eigenvectors.

In order to compute excited states, overlaps between the last wave-function snapshots\footnote{%
See flag \texttt{SNAPUPDATE} on page~\pageref{flag_snapupdate}, section~\ref{flag_snapupdate}.}
are used.  Note that, at some point $e^{-E_{i} \tau}$ will become 
smaller than the machine precision. Hence, for large enough values of $\tau$, the only wave-function
that can be recovered is the ground state and the overlap procedure will not be able to recover any excited states,
as we will see in detail in the next section~\ref{sec_examples}.

We improve the previous version of the code by adding the option of using a kinetic term $H_{K}^{(i)}$, $i = 1,2,3$,
based on the FFT scheme of section~\ref{subsec_kin}, instead of the one based on finite differences, $H_{K}^{(0)}$.
In all cases, the values of $\mathcal{A}$ and $\mathcal{B}$ are taken to be the same as in~\cite{Strickland:2009ft},
\begin{equation}
\mathcal{A} = \frac{1 - \frac{\Delta\tau}{2} V_{b}(x_{1},x_{2},x_{3})}{1 + \frac{\Delta\tau}{2} V_{b}(x_{1},x_{2},x_{3})} \,,\quad \mathcal{B} = \frac{1}{1 + \frac{\Delta\tau}{2} V_{b}(x_{1},x_{2},x_{3})} \,.
\end{equation}
The time step $\Delta\tau$ is an independent parameter that should be chosen as small as necessary
to reduce the discretization errors of the time evolution and achieve good convergence,
while being as large as permissible to minimize the computational expense. 
Both of these limitations depend on the combination of the kinetic term and the potential.  

\subsection{Parity projection}
\label{subsec_parity_proy}

We include a script\footnote{%
See section~\ref{subsec_post_scripts}, page~\pageref{subsec_post_scripts}.}
for computing the positive, $P^{+}$, negative, $P^{-}$, and negative-around-an-axis, $P^{-}_{p_{k}}$, parity projections of the wave-function,
\begin{align}
\label{eq_parity} & P^{\pm} \Psi(\vec{r}) = \frac{1}{2} \left(\Psi(\vec{r}) \pm \Psi(-\vec{r})\right) \,, \\
\label{eq_parity_axis} &\begin{aligned}
P^{-}_{\vec{p}_{k} = \hat{e}_{3}} \Psi(x_{1},x_{2},x_{3}) & = \frac{1}{2} \left[ \Psi(x_{1},x_{2},x_{3}) + \Psi(-x_{1},-x_{2},x_{3}) \right. \\
& \left. - \Psi(x_{1},x_{2},-x_{3}) - \Psi(-x_{1},-x_{2}-x_{3}) \right] \,.
\end{aligned}
\end{align}

Here, we have chosen $\vec{p}_{k} = \hat{e}_{3} = (0,0,1)$, although the script accepts $\hat{e}_{1}$ and $\hat{e}_{2}$ as well.
The aim is to separate the positive parity $s$-waves from the negative-parity $p$-waves.
For the $p$-wave states, the negative-around-an-axis, $P^-_{\vec{p}_k}$ separate the states with different magnetic quantum number $m \in \{0,\ \pm 1\}$. 
It is straightforward to extend these scripts for projection to the higher multipoles, if the need should arise. 

We include a second script that normalizes the wave-function via enforcing that $\int_{V} dV \lvert\Psi\rvert^2 = 1$, within the lattice volume $V$.

\section{Examples}
\label{sec_examples}

\subsection{Initial conditions and studied cases}
\label{subsec_examples}

\begin{figure}[ht]
\centering
\begin{minipage}{.54\textwidth}
\igg{potential_KIN0.png}
\end{minipage}
\begin{minipage}{.44\textwidth}
\caption{\label{fig_pots}
The Coulomb potential used for the following analysis.
The lattice spacing is $A=0.12$~fm. We use $N=64$ points in each spacial direction, and the potential is centered
in the lattice volume at the point $(N/2,N/2,N/2) = (32,32,32)$.}
\end{minipage}
\end{figure}

We discuss the convergence of the iterative procedure (see the previous sections~\ref{subsec_kin} and~\ref{subsec_iter}) in the following.
Furthermore, we study the performance of the overlap method for separating the excited states (see section~\ref{subsec_iter}),
as well as the behavior of the parity projections (see section~\ref{subsec_parity_proy}).

We use, unless stated otherwise, the Coulomb-like potential as our benchmark one,
setting \texttt{POTENTIAL=2} (see subsection~\ref{subsec_hpot}) as it is shown in Fig.~\ref{fig_pots}. 
This lattice implementation of the Coulomb potential has a discretization artifact at its center, leading to suppressed contributions from radii at the order of the lattice spacing, $r \sim A$. This is the dominant discretization artifact at small distances in all of the following examples. 
Along these lines we use the Coulomb-like initial conditions, setting \texttt{INITCONDTYPE=2} (see subsection~\ref{subsec_icond}).
We set the lattice spacing $A=0.12$~fm and use a reduced mass $M=0.3$~GeV. 
Due to $(A\cdot M)^2 \approx 0.03$, we expect mild finite mass discretization errors. 

In all of the cases, we do not plot points that are closer than $0.15N$ to one of the edges\footnote{%
$\texttt{max\_\allowbreak{}sep\_\allowbreak{}edge} = 0.3$ (see section~\ref{subsec_post_python}).}, in order to avoid zones where finite volume effects are prominent.
All of the wave-functions are normalized within the resulting box so that
\begin{equation}
\sum\limits_{r} \lvert\Psi(r)\rvert^2 = 1 \,,
\end{equation}
where the sum over $r$ evaluates over all of the points that are effectively plotted.

We plot the positive parity projection $P^+$ (Eq.~\ref{eq_parity}), and the negative parity around-an-axis $P^-_{\vec{p}_k}$.
We do not plot a projected wave-function if its weight is less than $10^{-5}$.
However, we need to stress that the overlapping procedure used for extracting ground and excited states out of the
iterative procedure (section~\ref{subsec_iter}) can also introduce numerical artefacts affecting excited states.

This section is organized as follows:
In subsection~\ref{subsec_ex_relativ}, we study the convergence of the relativistic $H_{K}^{(3)}$ kinetic term.
A selection of wave-functions at different $\tau$ is shown, as well as the evolution of the energies and the weights of the parity projections for each state.
In subsection~\ref{subsec_ex_nonrelativ}, we compare the non-relativistic kinetic terms $H_{K}^{(0)}$ (FDTD based), $H_{K}^{(1)}$, and $H_{K}^{(2)}$ (both FFT based).
In subsection~\ref{subsec_ex_relativ_vs_nonrelativ_high_mass}, we increase the reduced mass to $M=3$~GeV,
and compare the relativistic kinetic term $H_{K}^{(3)}$ with the non-relativistic one $H_{K}^{(2)}$.
In subsection~\ref{subsec_ex_increase_volume}, we compare the non-relativistic kinetic terms $H_{K}^{(0)}$ (FDTD based)
with $H_{K}^{(2)}$ (FFT based) at an increased lattice volume using $N=128$.
In subsection~\ref{subsec_ex_parity_fixing}, we study the parity fixing machinery of \texttt{quantumfdtd}, for both symmetrization and anti-symmetrization.
We use the non-relativistic kinetic term $H_{K}^{(0)}$ (FDTD) as an example.
Finally, in subsection~\ref{subsec_conv_conditions}, we summarize the convergence conditions and numerical behavior of the program.

\subsection{\texorpdfstring{The relativistic kinetic term $H_{K}^{(3)}$}{The relativistic kinetic term H\_K(3)}}
\label{subsec_ex_relativ}

In this section, the relativistic FFT based kinetic term $H_{K}^{(3)}$ (Eq.~\ref{eq:HK3}, section~\ref{subsec_kin}), is studied.
We use $N=64$ points in each direction and plot the convergence of the energies with $\tau$ (left plot in Fig.~\ref{fig_energ_KIN3}),
the relative weights of the parity projection for the 3 lowest states vs. $\tau$ (right plot in Fig.~\ref{fig_energ_KIN3} and Fig.~\ref{fig_proj_KIN3}).
Finally we plot the parity projections of the wave-functions for 4 selected values of $\tau$ (Figs.~\ref{fig_fig_KIN3_1} and~\ref{fig_fig_KIN3_2}).

Note that, for $\tau > 125$~GeV$^{-1}$, the weights of the excited states fall below the computational accuracy, $\epsilon$, of the computer.
For a 64 bit architecture, $\epsilon \sim 10^{-16}$.
Below this threshold we expect numerical artefacts which can be seen in the right plot of Fig.~\ref{fig_energ_KIN3},
where the exponential decay breaks for a weight of the minor component of the parity projection $\sim 10^{-21}$.
See also, on the left side of Fig.~\ref{fig_energ_KIN3}, that the energy of the first excited state decays to that of
the ground state because the overlap method is no longer able to separate the first excited state from the ground state once numerical artefacts show up.

For large values of $\tau$, the energy of the \emph{excited states} decays to that of the ground state, whose energy eigenvalue remains stable throughout the evolution.
The relative weights of the parity projections of the excited states decay exponentially, as expected from Eq.~\ref{iter_decay},
but their behavior starts to deviate from the exponential decay for high enough values of $\tau$, where numerical artefacts are expected.

Note also a first crucial point at $\tau \sim 60-65$~GeV$^{-1}$, where the energy, $E_2$, of the second excited state degenerates with $E_1$.
This can also be seen in the relative weights of the parity projections of the second excited state shown in the right plot of Fig.~\ref{fig_proj_KIN3}.

For high values of $\tau$, but not high enough to hit the aforementioned cancellation issues, the parity projected wave-functions are stable with $\tau$.
The negative parity projections around an axis, $P^-_{\vec{p}_k}$, of different separated states are compatible (they actually come from the same $2p$ state).
Additionally, there are 2 positive parity projections, $P^+$. The $1s$ ground state and the $2s$ excited state. This can be seen in Fig.~\ref{fig_fig_KIN3_2}.

In the captions of Figs.~\ref{fig_fig_KIN3_1} and~\ref{fig_fig_KIN3_2} we also list the average radii~in the box normalization, $\langle r \rangle^{\Box}$,
for each of the plotted wave-functions, computed after normalizing the wave functions in the effective lattice volume.
Their values lie well within the considered lattice volume such that we do not expect considerable finite volume effects.
Note also that, for $\tau > 100$~GeV$^{-1}$ (compare $\tau=104.72$ and $123.3$~GeV$^{-1}$ on Fig.~\ref{fig_fig_KIN3_2}),
these average radii become stable with $\tau$~(see the table associated with Fig.~\ref{fig_fig_KIN3_2}).
The average radii of the parity projections, $P^{+}\Psi_1$ and $P^{+}\Psi_2$, are compatible.
Similarly for the parity projections, $P^{-}_{\vec{p}_k}\Psi_1$ and $P^{-}_{\vec{p}_k}\Psi_2$.
That is, the overlapping procedure of section~\ref{subsec_iter} can separate the first excited states of each parity,
but cannot identify the higher states. 

Note that the typical radii in the box normalization are substantially smaller for the relativistic kinetic term than for the non-relativistic kinetic term (see subsection~\ref{subsec_ex_nonrelativ}). 
This is somewhat counter-intuitive given that the more compact wave-function with the Cornell potential is known to exhibit a slower decay for the relativistic kinetic term compared to the non-relativistic kinetic term (see section~\ref{intro}).

\begin{figure}[ht]
\centering
\ig{energ_KIN3.png} \hfill
\ig{proj_KIN3_0.png} 
\caption{\label{fig_energ_KIN3}
Left: Energies obtained as a function of $\tau$ using the relativistic kinetic term $H_{K}^{(3)}$.\\
Right: Weights of the parity projection of the corresponding ground state as a function of $\tau$.}
\end{figure}

\begin{figure}[ht]
\centering
\ig{proj_KIN3_1.png} \hfill
\ig{proj_KIN3_2.png}
\caption{\label{fig_proj_KIN3}
Weights of the parity projections of the first and second excited state as in Fig~\ref{fig_energ_KIN3}.}
\end{figure}

\begin{figure}[ht]
\centering
\ig{figure_KIN3_500.png} \hfill
\ig{figure_KIN3_4500.png} 
\caption{\label{fig_fig_KIN3_1}
Parity-projected wave-functions as a function of $r/A$ using the relativistic kinetic term $H_{K}^{(3)}$ at two different snapshot time steps.
The average radii in units of $A$ are given by}
\small
\begin{tabular}{r|c|c|c|c|c|c}
& $P_{\vec{p}_{k}}^{-} \Psi_{0}$ & $P^{+} \Psi_{0}$ & $P_{\vec{p}_{k}}^{-} \Psi_{1}$ & $P^{+} \Psi_{1}$ & $P_{\vec{p}_{k}}^{-} \Psi_{2}$ & $P^{+} \Psi_{2}$ \\
\hline
$\langle r \rangle^{\Box}_{\text{left}}$ & $18.17$ & $6.11$ & $19.58$ & $4.97$ & $10.26$ & $8.78$ \\
\hline
$\langle r \rangle^{\Box}_{\text{right}}$ & $-$ & $3.36$ & $15.95$ & $16.49$ & $16.28$ & $11.83$
\end{tabular}
\end{figure}

\begin{figure}[ht]
\centering
\ig{figure_KIN3_8500.png} \hfill
\ig{figure_KIN3_10000.png} 
\caption{\label{fig_fig_KIN3_2}
Parity-projected wave-functions as a function of $r/A$ using the relativistic kinetic term $H_{K}^{(3)}$ at two different snapshot time steps as in Fig.~\ref{fig_fig_KIN3_1}.
The average radii in units of $A$ are given by}
\small
\begin{tabular}{r|c|c|c|c|c|c}
& $P_{\vec{p}_{k}}^{-} \Psi_{0}$ & $P^{+} \Psi_{0}$ & $P_{\vec{p}_{k}}^{-} \Psi_{1}$ & $P^{+} \Psi_{1}$ & $P_{\vec{p}_{k}}^{-} \Psi_{2}$ & $P^{+} \Psi_{2}$ \\
\hline
$\langle r \rangle^{\Box}_{\text{left}}$ & $-$ & $3.36$ & $15.79$ & $15.86$ & $15.79$ & $15.86$ \\
\hline
$\langle r \rangle^{\Box}_{\text{right}}$ & $-$ & $3.36$ & $15.78$ & $15.87$ & $15.78$ & $15.87$
\end{tabular}
\end{figure}

\FloatBarrier

\subsection{Comparison between the different non-relativistic kinetic terms}
\label{subsec_ex_nonrelativ}

We now compare the three non-relativistic kinetic terms introduced in section~\ref{subsec_kin}:
$H_{K}^{(0)}$, (Eq.~\ref{eq:HK0}), based on the FDTD with Dirichlet boundary conditions; $H_{K}^{(1)}$, (Eq.~\ref{eq:HK1}), based on a discretized FFT implementation of the continuum kinetic term $H_{K}^{nr}$, (Eq.~\ref{eq:HKnr}); and $H_{K}^{(2)}$, (Eq.~\ref{eq:HK2}), based on the discretized FFT implementation of the FDTD kinetic term. Both $H_{K}^{(1)}$ and $H_{K}^{(2)}$ use periodic boundary conditions.

In the wave-function plots shown below\footnote{%
Right plots of Figs.~\ref{fig_comp_KIN_0} to~\ref{fig_comp_KIN_2} and all of the plots of Figs.~\ref{fig_fi_comp_KIN_0} to~\ref{fig_fi_comp_KIN_2}.}, we show with dashed lines the analytic solutions using the \emph{non-relativistic} kinetic term with a Coulomb potential, whereas the solid lines are those computed via the \texttt{quantumfdtd} program.
We use $N=64$ points in each direction.
The analytic solutions are the well known eigen-solutions of the Hydrogen atom,
\begin{equation}
\label{eq:non-rel-theo}
R_{1s}(r) = N_{1s} e^{-\mu r} \,, \quad R_{2s}(r) = N_{2s} (2-\mu r) e^{-\mu r/2} \,, \quad R_{2p}(r) = N_{2p} \mu r e^{-\mu r/2} \,,
\end{equation}
where $\mu$ is the reduced mass and the continuum normalization is given by
\begin{equation}
N_{1s}^{\infty} = 2\mu^{3/2} \,, \quad N_{2s}^{\infty} = \frac{\mu^{3/2}}{2\sqrt{2}} \,, \quad N_{2p}^{\infty} = \frac{\mu^{3/2}}{2\sqrt{6}}\,.
\end{equation}
However, for the sake of comparability, we use the box normalization
\begin{equation}
\int_{V} dV |R(r)|^{2} = \sum_{x \in V} |R(x)|^{2} = 1 \,.
\end{equation}
The corresponding average radii in the continuum are given by
\begin{align}
\label{eq:analytic_avg_radii}
\langle r \rangle^{\infty}_{1s} = \frac{3}{2\mu} \,, \quad \langle r \rangle^{\infty}_{2s} = \frac{6}{\mu} \,, \quad \langle r \rangle^{\infty}_{2p} = \frac{5}{\mu} \,.
\end{align}
We list the average radii obtained using the box normalization as $\langle r \rangle^{\Box}$.

Although the $1s$, $2p$, and $2s$ states can be recognized, and the wave-functions that can be extracted via parity projection from different \texttt{quantumfdtd} states (see section~\ref{subsec_parity_proy}) are compatible, there are some differences between the 3 kinetic terms, due to finite volume errors or due to the discretization errors.
A study of the former is performed in the next subsection~\ref{subsec_ex_increase_volume} where we increase the lattice volume to $N=128$.

For the FDTD kinetic term, $H_{K}^{(0)}$, the most prominent finite volume effect is a faster than expected decay of the $2p$ and $2s$ wave-functions close to the borders of the lattice volume which can be seen in Figs.~\ref{fig_comp_KIN_0} and~\ref{fig_fi_comp_KIN_0}.
This is expected because of the Dirichlet boundary conditions.
The $1s$ wave-function decays exponentially, so that its theoretical value at the borders of the finite volume is also low and therefore, this wave-function is not as severely affected by the Dirichlet boundary conditions.

For the FFT kinetic terms, $H_{K}^{(1)}$ and $H_{K}^{(2)}$, the most prominent finite size effects can be seen in Figs.~\ref{fig_comp_KIN_1}, \ref{fig_comp_KIN_2}, \ref{fig_fi_comp_KIN_1}, and~\ref{fig_fi_comp_KIN_2} and are located at values of $r$ close to the lattice spacing $A$. 
We remind that this artifact is due to the lattice implementation of the Coulomb potential.
The periodic boundary conditions smoothen the behavior of the wave-function close to the borders, but at the expense of being less able to minimize the discretization artifacts for values close to the origin, where the continuum potential is actually divergent.
This effect is not present for the $2p$ wave-function, that vanishes at the origin.
But it is especially prominent for the $2s$ wave-function.

The continuum average radii in units of $A$ are given by $\langle r \rangle^{\infty}_{1s} = 8.21$, $\langle r \rangle^{\infty}_{2s} = 32.83$, and $\langle r \rangle^{\infty}_{2p} = 27.36$.
However, in the considered effective lattice volumes, the average radii of the continuum solutions of Eq.~\ref{eq:non-rel-theo} (in units of $A$) are $\langle r \rangle^{\Box}_{1s} = 8.14$, $\langle r \rangle^{\Box}_{2s} = 20.30$, and $\langle r \rangle^{\Box}_{2p} = 19.06$.
The lower values are linked to finite volume effects, that make larger values of $r$ being discarded.
These box normalized theoretical average radii should be compared with those listed in the captions of Figs.~\ref{fig_fi_comp_KIN_0} to~\ref{fig_fi_comp_KIN_2}.
The wave-functions of figures~\ref{fig_comp_KIN_0} to~\ref{fig_comp_KIN_2}, with $\tau=6.16$~GeV$^{-1}$, are still not stabilized on $\tau$.
Hence, there are big differences between theoretical and simulated values.
This is especially prominent for the $1s$ state, that still has non-vanishing components of positive parity higher states with larger radii like the $2s$.

Things become more interesting with figures~\ref{fig_fi_comp_KIN_0} to~\ref{fig_fi_comp_KIN_2}, where the left plot is $\tau=61.6$~GeV$^{-1}$ and the right one, $\tau=123.2$~GeV$^{-1}$.
The values of $\langle r \rangle^{\Box}_{P^+\Psi_0}$ get closer to $\langle r \rangle^{\Box}_{1s}$, the ones of $\langle r \rangle^{\Box}_{P^-_{\vec{p}_k}\Psi_{1,2}}$ get closer to $\langle r \rangle^{\Box}_{2p}$, and the ones of $\langle r \rangle^{\Box}_{P^+\Psi_{1,2}}$ get closer to $\langle r \rangle^{\Box}_{2s}$.
However, note that $\langle r \rangle^{\Box}_{2p} \approx 20.30$ is pretty close to $\langle r \rangle^{\Box}_{2s}\approx 19.06$.
Indeed, the only difference between these states is the angular quantum number $l$, while $n=1$ is equal.
This fact makes the comparisons harder.
Anyhow, the average radii of the wave-functions of the FFT kinetic terms $H_{K}^{(1)}$ (Figs.~\ref{fig_fi_comp_KIN_1}) and $H_{K}^{(2)}$ (Fig.~\ref{fig_fi_comp_KIN_2}) are pretty similar.
But the FDTD kinetic term $H_{K}^{(0)}$ (Fig.~\ref{fig_fi_comp_KIN_0}) shows smaller average radii than the FFT ones.
This is especially prominent on the $2s$ and $2p$ states.

Looking at the tail of the wave-functions, the decay of the 
excited state wave-functions is faster in the FDTD case.
This fast decay is explicitly enforced by the Dirichlet boundary conditions.
For the FFT kinetic terms, periodic boundary conditions are used.
Such boundary conditions do not enforce $\Psi(\text{boundary}) \equiv 0$, and the tail of the corresponding wave-functions
are closer to the continuum ones as can be seen in Figs.~\ref{fig_fi_comp_KIN_1} and~\ref{fig_fi_comp_KIN_2}.
The effect of the Dirichlet boundary conditions would be an unphysically low value of the average radii, an effect that is shown in the table in Fig.~\ref{fig_fi_comp_KIN_0}.
The $1s$ state is less affected by finite volume effects because its theoretical average radius $\langle r \rangle^{\Box}_{1s}\approx 8.14$ is well within the finite lattice volume.

\begin{figure}[ht]
\centering
\ig{energ_KIN0.png} \hfill
\ig{figure_KIN0_500.png}
\caption{\label{fig_comp_KIN_0}
Left: Energies obtained as a function of $\tau$ using the non-relativistic kinetic term $H_{K}^{(0)}$.\\
Right: Parity-projected wave-functions as a function of $r/A$ using the non-relativistic kinetic term $H_{K}^{(0)}$.
The broad line of the $P^-_{\vec{p}_k}\Psi_2$ wave-function denotes lack of spherical symmetry, which is explicitly broken by the discretization.
For low values of $\tau$ (as in this plot, $\tau=6.16$~GeV$^{-1}$), higher states with quantum number $l>0$ could show up as broad lines on the plots.
Such states decay exponentially with $\tau$, according to Eq.~\ref{iter_decay}, so that they do not show up for higher values of $\tau$.
In particular, the negative parity-projected ground state $P^-_{\vec{p}_k}\Psi_0$ is dominated by the $\Psi_{2p}$ contribution, and the positive parity-projected first excited state $P^+\Psi_1$ is actually dominated by the true ground state $\Psi_{1s}$.
The average radii in units of $A$ are given by}
\small
\begin{tabular}{r|c|c|c|c|c|c|c|c|c}
& $P_{\vec{p}_{k}}^{-} \Psi_{0}$ & $P^{+} \Psi_{0}$ & $P_{\vec{p}_{k}}^{-} \Psi_{1}$ & $P^{+} \Psi_{1}$ & $P_{\vec{p}_{k}}^{-} \Psi_{2}$ & $P^{+} \Psi_{2}$ & $\Psi_{1s}^{\text{th}}$ & $\Psi_{2s}^{\text{th}}$ & $\Psi_{2p}^{\text{th}}$ \\
\hline
$\langle r \rangle^{\Box}$ & $18.74$ & $15.58$ & $17.25$ & $9.60$ & $21.41$ & $13.94$ & $8.14$ & $20.30$ & $19.06$
\end{tabular}
\end{figure}

\begin{figure}[ht]
\centering
\ig{energ_KIN1.png} \hfill
\ig{figure_KIN1_500.png}
\caption{\label{fig_comp_KIN_1}
Energies and parity projected wave-functions as in Fig.~\ref{fig_comp_KIN_0} using the non-relativistic kinetic term $H_{K}^{(1)}$. The average radii in units of $A$ are given by}
\small
\begin{tabular}{r|c|c|c|c|c|c|c|c|c}
& $P_{\vec{p}_{k}}^{-} \Psi_{0}$ & $P^{+} \Psi_{0}$ & $P_{\vec{p}_{k}}^{-} \Psi_{1}$ & $P^{+} \Psi_{1}$ & $P_{\vec{p}_{k}}^{-} \Psi_{2}$ & $P^{+} \Psi_{2}$ & $\Psi_{1s}^{\text{th}}$ & $\Psi_{2s}^{\text{th}}$ & $\Psi_{2p}^{\text{th}}$ \\
\hline
$\langle r \rangle^{\Box}$ & $18.85$ & $16.66$ & $24.50$ & $9.60$ & $20.11$ & $9.26$ & $8.14$ & $20.30$ & $19.06$
\end{tabular}
\end{figure}

\begin{figure}[ht]
\centering
\ig{energ_KIN2.png} \hfill
\ig{figure_KIN2_500.png}
\caption{\label{fig_comp_KIN_2}
Energies and parity projected wave-functions as in Fig.~\ref{fig_comp_KIN_0} using the non-relativistic kinetic term $H_{K}^{(2)}$. The average radii in units of $A$ are given by}
\small
\begin{tabular}{r|c|c|c|c|c|c|c|c|c}
& $P_{\vec{p}_{k}}^{-} \Psi_{0}$ & $P^{+} \Psi_{0}$ & $P_{\vec{p}_{k}}^{-} \Psi_{1}$ & $P^{+} \Psi_{1}$ & $P_{\vec{p}_{k}}^{-} \Psi_{2}$ & $P^{+} \Psi_{2}$ & $\Psi_{1s}^{\text{th}}$ & $\Psi_{2s}^{\text{th}}$ & $\Psi_{2p}^{\text{th}}$ \\
\hline
$\langle r \rangle^{\Box}$ & $18.84$ & $16.55$ & $24.47$ & $9.40$ & $20.07$ & $9.29$ & $8.14$ & $20.30$ & $19.06$
\end{tabular}
\end{figure}

\begin{figure}[ht]
\centering
\ig{figure_KIN0_5000.png} \hfill
\ig{figure_KIN0_10000.png}
\caption{\label{fig_fi_comp_KIN_0}
Parity projected wave-functions as in Fig.~\ref{fig_comp_KIN_0} using the non-relativistic kinetic term $H_{K}^{(0)}$ at two later snapshot time steps. The average radii in units of $A$ are given by}
\small
\begin{tabular}{r|c|c|c|c|c|c|c|c|c}
& $P_{\vec{p}_{k}}^{-} \Psi_{0}$ & $P^{+} \Psi_{0}$ & $P_{\vec{p}_{k}}^{-} \Psi_{1}$ & $P^{+} \Psi_{1}$ & $P_{\vec{p}_{k}}^{-} \Psi_{2}$ & $P^{+} \Psi_{2}$ & $\Psi_{1s}^{\text{th}}$ & $\Psi_{2s}^{\text{th}}$ & $\Psi_{2p}^{\text{th}}$ \\
\hline
$\langle r \rangle^{\Box}_{\text{left}}$ & $17.64$ & $8.36$ & $17.64$ & $19.09$ & $16.87$ & $18.96$ & \multirow{2}{*}{$8.14$} & \multirow{2}{*}{$20.30$} & \multirow{2}{*}{$19.06$} \\
\cline{1-7}
$\langle r \rangle^{\Box}_{\text{right}}$ & $-$ & $8.38$ & $17.63$ & $18.90$ & $17.63$ & $18.90$ &  &  & 
\end{tabular}
\end{figure}

\begin{figure}[ht]
\centering
\ig{figure_KIN1_5000.png} \hfill
\ig{figure_KIN1_10000.png}
\caption{\label{fig_fi_comp_KIN_1}
Parity projected wave-functions as in Fig.~\ref{fig_comp_KIN_1} using the non-relativistic kinetic term $H_{K}^{(1)}$ at two later snapshot time steps. The average radii in units of $A$ are given by}
\small
\begin{tabular}{r|c|c|c|c|c|c|c|c|c}
& $P_{\vec{p}_{k}}^{-} \Psi_{0}$ & $P^{+} \Psi_{0}$ & $P_{\vec{p}_{k}}^{-} \Psi_{1}$ & $P^{+} \Psi_{1}$ & $P_{\vec{p}_{k}}^{-} \Psi_{2}$ & $P^{+} \Psi_{2}$ & $\Psi_{1s}^{\text{th}}$ & $\Psi_{2s}^{\text{th}}$ & $\Psi_{2p}^{\text{th}}$ \\
\hline
$\langle r \rangle^{\Box}_{\text{left}}$ & $18.46$ & $8.38$ & $18.45$ & $21.32$ & $18.41$ & $21.99$ & \multirow{2}{*}{$8.14$} & \multirow{2}{*}{$20.30$} & \multirow{2}{*}{$19.06$} \\
\cline{1-7}
$\langle r \rangle^{\Box}_{\text{right}}$ & $-$ & $8.49$ & $18.47$ & $20.93$ & $18.47$ & $21.70$ &  &  & 
\end{tabular}
\end{figure}

\begin{figure}[!ht]
\centering
\ig{figure_KIN2_5000.png} \hfill
\ig{figure_KIN2_10000.png}
\caption{\label{fig_fi_comp_KIN_2}
Parity projected wave-functions as in Fig.~\ref{fig_comp_KIN_2} using the non-relativistic kinetic term $H_{K}^{(2)}$ at two later snapshot time steps. The average radii in units of $A$ are given by}
\small
\begin{tabular}{r|c|c|c|c|c|c|c|c|c}
& $P_{\vec{p}_{k}}^{-} \Psi_{0}$ & $P^{+} \Psi_{0}$ & $P_{\vec{p}_{k}}^{-} \Psi_{1}$ & $P^{+} \Psi_{1}$ & $P_{\vec{p}_{k}}^{-} \Psi_{2}$ & $P^{+} \Psi_{2}$ & $\Psi_{1s}^{\text{th}}$ & $\Psi_{2s}^{\text{th}}$ & $\Psi_{2p}^{\text{th}}$ \\
\hline
$\langle r \rangle^{\Box}_{\text{left}}$ & $18.44$ & $8.33$ & $18.44$ & $21.29$ & $18.39$ & $21.97$ &\multirow{2}{*}{$8.14$} & \multirow{2}{*}{$20.30$} & \multirow{2}{*}{$19.06$} \\
\cline{1-7}
$\langle r \rangle^{\Box}_{\text{right}}$ & $-$ & $8.43$ & $18.45$ & $20.93$ & $18.45$ & $20.46$ &  &  & 
\end{tabular}
\end{figure}

\FloatBarrier

\subsection{\texorpdfstring{Comparison of the relativistic $H_{K}^{(3)}$ and the non-relativistic $H_{K}^{(2)}$ kinetic terms at high mass}{Comparison of the relativistic H\_K(3) and the non-relativistic H\_K(2) kinetic terms at high mass}}
\label{subsec_ex_relativ_vs_nonrelativ_high_mass}

In this subsection, we increase the reduced mass up to $M=3$~GeV~for demonstration purposes.
In all of the following figures, Figs.~\ref{fig_energ_HEAVYMASS} to~\ref{fig_figure_HEAVYMASS_2}, the relativistic kinetic term,
$H_{K}^{(3)}$, (Eq.~\ref{eq:HK3}) is being used in the left plots, whereas the non-relativistic one, $H_{K}^{(2)}$,
(Eq.~\ref{eq:HK2}) is being used in the right plots, respectively.

Our next goal is to compare the numerical behavior of the relativistic and the non-relativistic kinetic terms
in the regime of high mass, where they should become more similar.
Comparing these results with the results from the previous subsections~\ref{subsec_ex_relativ} and~\ref{subsec_ex_nonrelativ}
shows that the increase of the reduced mass from $M=0.3$~GeV to $M=3$~GeV leads to the expected decrease of the characteristic average radius of the wave-functions.
Due to $(AM)^2 \approx 3.3$ we expect large finite mass discretization errors. For this reason, this setup is not appropriate for production purposes, where one should reduce the lattice spacing by (approximately) the inverse factor from $A=0.12$ fm to $A \lesssim 0.02$ fm. 
The energies and parity projections of the wave-functions evolve faster with $\tau$ (see Figs.~\ref{fig_energ_HEAVYMASS} to \ref{fig_proj_HEAVYMASS_2}),
such that the energy of the second excited state, $E_2$, degenerates with $E_{1}$ at $\tau \sim 25$~GeV$^{-1}$
and $E_1$ degenerates with $E_{0}$ at $\tau \sim 70$~GeV$^{-1}$.
For $\tau > 60-65$~GeV$^{-1}$, numerical cancellation effects start showing up (see Figs.~\ref{fig_proj_HEAVYMASS_0} to~\ref{fig_proj_HEAVYMASS_2}).

As expected, the wave-functions (see Fig.~\ref{fig_figure_HEAVYMASS_1} and~\ref{fig_figure_HEAVYMASS_2}) of the relativistic and
the non-relativistic cases are closer to one another due to the increased reduced mass.
However, the results are visibly non-smooth due to discretization artifacts.
Note that the increase on the reduced mass $M$ by an order of magnitude has decreased the characteristic average radii of the wave-functions.
For the non-relativistic case, $\langle r \rangle^{\infty}_{1s} = 0.82$ in units of $A$, below a single lattice spacing.
For the $2s$ and $2p$ excited states, the continuum average radii are $\langle r \rangle^{\infty}_{2s} = 3.28$
and $\langle r \rangle^{\infty}_{2p} = 2.74$, respectively. This strongly impacts a precise analysis of this case.
On the right side of Figs.~\ref{fig_figure_HEAVYMASS_2} and~\ref{fig_figure_HEAVYMASS_2}, it can be seen that,
although the wave-function computed via the non-relativistic FFT kinetic term $H_{K}^{(2)}$ are close to the continuum limit, there are clear discrepancies that originate in the discretization of the potential.
This is very prominent for the $1s$ and $2s$ states, that have a non-vanishing value close to the origin. 
Anyhow, in Figs.~\ref{fig_figure_HEAVYMASS_1} to~\ref{fig_figure_HEAVYMASS_3}, it can be seen that the wave-functions computed
via the relativistic kinetic term, $H_{K}^{(3)}$, and the non-relativistic one, $H_{K}^{(2)}$,
(both based on a FFT with periodic boundary conditions) are pretty close.
Compare with the cases of right sides of Figs.~\ref{fig_fi_comp_KIN_2} (right side, non-relativistic $H_{H}^{(2)}$, step=10000, $M=0.3$~GeV)
and~\ref{fig_fig_KIN3_2} (right side, relativistic kinetic term $H_{K}^{(3)}$, step=10000),
where the differences between relativistic and non-relativistic are quite prominent.
The similarities can be also seen in Fig.~\ref{fig_energ_HEAVYMASS}, where the decay curves are quite similar
after taking into account an additive factor of $3$~GeV on the energies due to the rest mass on the relativistic approach.

\begin{figure}[ht]
\centering
\ig{energ_KIN3_HEAVYMASS.png} \hfill
\ig{energ_KIN2_HEAVYMASS.png} 
\caption{\label{fig_energ_HEAVYMASS}
Energies obtained as a function of $\tau$ using the relativistic kinetic term $H_{K}^{(3)}$ (left)
and $H_{K}^{(2)}$ (right) using the increased reduced mass $M = 3$~GeV.}
\end{figure}

\begin{figure}[ht]
\centering
\ig{proj_KIN3_HEAVYMASS_0.png} \hfill
\ig{proj_KIN2_HEAVYMASS_0.png} 
\caption{\label{fig_proj_HEAVYMASS_0}
Weights of the parity projection of the ground state as a function of $\tau$ using
the relativistic kinetic term $H_{K}^{(3)}$ (left) and $H_{K}^{(2)}$ (right) using the increased reduced mass $M = 3$~GeV.}
\end{figure}

\begin{figure}[ht]
\centering
\ig{proj_KIN3_HEAVYMASS_1.png} \hfill
\ig{proj_KIN2_HEAVYMASS_1.png} 
\caption{\label{fig_proj_HEAVYMASS_1}
Weights of the parity projection of the first excited state as in Fig.~\ref{fig_proj_HEAVYMASS_0}. Left: $H_{K}^{(3)}$; right: $H_{K}^{(2)}$.}
\end{figure}

\begin{figure}[ht]
\centering
\ig{proj_KIN3_HEAVYMASS_2.png} \hfill
\ig{proj_KIN2_HEAVYMASS_2.png} 
\caption{\label{fig_proj_HEAVYMASS_2}
Weights of the parity projection of the second excited state as in Fig.~\ref{fig_proj_HEAVYMASS_0}. Left: $H_{K}^{(3)}$; right: $H_{K}^{(2)}$.}
\end{figure}

\begin{figure}[ht]
\centering
\ig{figure_KIN3_HEAVYMASS_500.png} \hfill
\ig{figure_KIN2_HEAVYMASS_500.png}
\caption{\label{fig_figure_HEAVYMASS_1}
Parity-projected wave-functions as a function of $r/A$ using the relativistic kinetic term $H_{K}^{(3)}$ (left)
and the non-relativistic kinetic term $H_{K}^{(2)}$ (right) using the increased reduced mass $M = 3$~GeV.
The average radii in units of $A$ are given by}
\small
\begin{tabular}{r|c|c|c|c|c|c|c|c|c}
& $P_{\vec{p}_{k}}^{-} \Psi_{0}$ & $P^{+} \Psi_{0}$ & $P_{\vec{p}_{k}}^{-} \Psi_{1}$ & $P^{+} \Psi_{1}$ & $P_{\vec{p}_{k}}^{-} \Psi_{2}$ & $P^{+} \Psi_{2}$ & $\Psi_{1s}^{\text{th}}$ & $\Psi_{2s}^{\text{th}}$ & $\Psi_{2p}^{\text{th}}$ \\
\hline
$\langle r \rangle^{\Box}_{\text{left}}$ & $2.24$ & $1.40$ & $2.35$ & $2.20$ & $1.63$ & $3.28$ & $-$ & $-$ & $-$ \\
\hline
$\langle r \rangle^{\Box}_{\text{right}}$ & $2.60$ & $1.49$ & $2.63$ & $2.54$ & $1.78$ & $3.14$ & 1.02 & 3.39 & 2.74 
\end{tabular}
\end{figure}

\begin{figure}[ht]
\centering
\ig{figure_KIN3_HEAVYMASS_1000.png} \hfill
\ig{figure_KIN2_HEAVYMASS_1000.png} 
\caption{\label{fig_figure_HEAVYMASS_2}
Parity-projected wave-functions as in Fig.~\ref{fig_figure_HEAVYMASS_1} at a later snapshot time step. Left: $H_{K}^{(3)}$; right: $H_{K}^{(2)}$.
The average radii in units of $A$ are given by}
\small
\begin{tabular}{r|c|c|c|c|c|c|c|c|c}
& $P_{\vec{p}_{k}}^{-} \Psi_{0}$ & $P^{+} \Psi_{0}$ & $P_{\vec{p}_{k}}^{-} \Psi_{1}$ & $P^{+} \Psi_{1}$ & $P_{\vec{p}_{k}}^{-} \Psi_{2}$ & $P^{+} \Psi_{2}$ & $\Psi_{1s}^{\text{th}}$ & $\Psi_{2s}^{\text{th}}$ & $\Psi_{2p}^{\text{th}}$ \\
\hline
$\langle r \rangle^{\Box}_{\text{left}}$ & $2.13$ & $1.27$ & $2.15$ & $3.52$ & $2.11$ & $3.52$ & $-$ & $-$ & $-$ \\
\hline
$\langle r \rangle^{\Box}_{\text{right}}$ & $2.57$ & $1.37$ & $2.57$ & $3.90$ & $2.11$ & $3.75$ & 1.02 & 3.39 & 2.74 \end{tabular}
\end{figure}

\begin{figure}[ht]
\centering
\ig{figure_KIN3_HEAVYMASS_3000.png} \hfill
\ig{figure_KIN2_HEAVYMASS_3000.png} 
\caption{\label{fig_figure_HEAVYMASS_3}
Parity-projected wave-functions as in Fig.~\ref{fig_figure_HEAVYMASS_1} at a later snapshot time step. Left: $H_{K}^{(3)}$; right: $H_{K}^{(2)}$.
The average radii in units of $A$ are given by}
\begin{tabular}{r|c|c|c|c|c|c|c|c|c}
& $P_{\vec{p}_{k}}^{-} \Psi_{0}$ & $P^{+} \Psi_{0}$ & $P_{\vec{p}_{k}}^{-} \Psi_{1}$ & $P^{+} \Psi_{1}$ & $P_{\vec{p}_{k}}^{-} \Psi_{2}$ & $P^{+} \Psi_{2}$ & $\Psi_{1s}^{\text{th}}$ & $\Psi_{2s}^{\text{th}}$ & $\Psi_{2p}^{\text{th}}$ \\
\hline
$\langle r \rangle^{\Box}_{\text{left}}$ & $-$ & $1.27$ & $2.11$ & $3.65$ & $2.11$ & $3.65$ & $-$ & $-$ & $-$ \\
\hline
$\langle r \rangle^{\Box}_{\text{right}}$ & $-$ & $1.37$ & $2.56$ & $4.07$ & $2.56$ & $4.07$ & 1.02 & 3.39 & 2.74 \end{tabular}
\end{figure}

\FloatBarrier

\subsection{\texorpdfstring{Comparison between the non-relativistic kinetic terms at larger $N$}{Comparison between the non-relativistic kinetic terms at larger N}}
\label{subsec_ex_increase_volume}

In this subsection we compare the FDTD kinetic term, $H_{K}^{(0)}$, (left plots in all figures)
with the FFT kinetic term, $H_{K}^{(2)}$, (right plots in all figures) for a increased lattice volume, using $N=128$ and reduced mass $M=0.3$~GeV.
This should be compared with the $N=64$ case from the previous section~\ref{subsec_ex_nonrelativ}.
The dashed lines again represent the analytic solutions for the \emph{non-relativistic} kinetic term with Coulomb potential.

From Figs.~\ref{fig_energ_comp_KIN_128} to~\ref{fig_figure_comp_KIN_128_2}, we plot the FDTD-based non-relativistic kinetic term,
$H_{K}^{0}$ (Dirichlet boundary conditions), on the left vs. the non-relativistic FFT-based, $H_{K}^{(2)}$ (periodic boundary conditions), on the right.

The continuum average radii in units of $A$ using the box normalization are given by $\langle r \rangle^{\Box}_{1s} = 8.22$, $\langle r \rangle^{\Box}_{2s} = 30.77$, and $\langle r \rangle^{\Box}_{2p} = 26.18$, respectively.
As it happened for $N=64$ in section~\ref{subsec_ex_nonrelativ}, the infinite volume continuum solutions of Eq.~\ref{eq:non-rel-theo} in units of $A$ have larger average radii: $\langle r \rangle^{\infty}_{1s} = 8.21$, $\langle r \rangle^{\infty}_{2s} = 32.83$, and $\langle r \rangle^{\infty}_{2p} = 27.36$, respectively.
However, note that the difference is not as prominent as in the $N=64$ case.
Actually, $r/A=31$ falls well inside the lattice volume of $N^3=128^3$, so that finite volume effects are not so dramatic.

Indeed, when comparing the results with the ones from the previous sections, we see that the finite volume effects are greatly reduced.
Even if we are comparing Dirichlet boundary conditions (FDTD kinetic term) vs. periodic boundary conditions (FFT kinetic term).
Additionally, the FDTD ($H_{K}^{(0)}$) and FFT ($H_{K}^{(2)}$) based solutions agree with one another and with the theoretical predictions (dashed lines).
This can especially be seen in Fig.~\ref{fig_figure_comp_KIN_128_2}, and on the associated table of average radii.
Compare with the previous Figs.~\ref{fig_fi_comp_KIN_0} (FDTD, Dirichlet boundary conditions)
and~\ref{fig_fi_comp_KIN_2} (FFT, periodic boundary conditions) in section~\ref{subsec_ex_nonrelativ} ($N=64$).
In particular, the problem with the tail behavior of Figs.~\ref{fig_fi_comp_KIN_0}, that was attributed to finite volume effects
and the Dirichlet boundary conditions, disappears when the volume is increased so that it can accommodate most of the wave-functions.

The computed values for the average radii in the box normalization, $\langle r \rangle^{\Box}_{P^+\Psi_{0,1,2}}$,
are still systematically higher than the corresponding theoretical values $\langle r\rangle^{\Box}_{1s,2s}$,
with the FDTD (Dirichlet boundary conditions)
being more accurate than the FFT (periodic boundary conditions). This is due to discretization errors close to $r\sim A$,
where $1s$ and $2s$ wave-functions have a maxium at the origin. In the case of $\langle r \rangle^{\Box}_{P^-_{\vec{p}_k}\Psi_{1,2}}$ and the theoretical $\langle r \rangle^{\Box}_{2p}$, the computed value is smaller
than the corresponding theortical value because it has a node at the origin.

On Fig.~\ref{fig_figure_comp_KIN_128_2} it can be seen that, especially for the FFT case (right side),
the behavior of the $2s$ and $2p$ of the computed wave-functions deviate from the continuum limit close to the origin.
The continuum limit curves are higher than the computed ones close to the origin.
This points to an issue with the lattice implementation of the Coulomb potential on the origin as mentioned earlier.
This is to be expected, since the Coulomb potential has a pole on the origin.
Hence, we are facing discretization artifacts originating from the discretization of the potential.
The FFT approach seems slightly worse than the FDTD approach in handling this kind of effect.

\begin{figure}[ht]
\centering
\ig{energ_KIN0_128.png} \hfill
\ig{energ_KIN2_128.png}
\caption{\label{fig_energ_comp_KIN_128}
Energies obtained as a function of $\tau$ using the non-relativistic kinetic terms $H_{K}^{(0)}$ (FDTD-based, left)
and the non-relativistic kinetic term $H_{K}^{(2)}$ (FFT-based, right) using the increased lattice volume ($N = 128$).}
\end{figure}

\begin{figure}[ht]
\centering
\ig{proj_KIN0_128_0.png} \hfill
\ig{proj_KIN2_128_0.png} 
\caption{\label{fig_proj_comp_KIN_128_0}
Weights of the parity projection of the ground state as a function of $\tau$ using the non-relativistic kinetic terms $H_{K}^{(0)}$
(FDTD-based, left) and the non-relativistic kinetic term $H_{K}^{(2)}$ (FFT-based, right) using the increased lattice volume ($N = 128$).}
\end{figure}

\begin{figure}[ht]
\centering
\ig{proj_KIN0_128_1.png} \hfill
\ig{proj_KIN2_128_1.png}
\caption{\label{fig_proj_comp_KIN_128_1}
Weights of the parity projection of the first excited state as in Fig.~\ref{fig_proj_comp_KIN_128_0}.}
\end{figure}

\begin{figure}[ht]
\centering
\ig{proj_KIN0_128_2.png} \hfill
\ig{proj_KIN2_128_2.png}
\caption{\label{fig_proj_comp_KIN_128_2}
Weights of the parity projection of the second excited state as in Fig.~\ref{fig_proj_comp_KIN_128_0}.}
\end{figure}

\begin{figure}[ht]
\centering
\ig{figure_KIN0_128_1500.png} \hfill
\ig{figure_KIN2_128_1500.png}
\caption{\label{fig_figure_comp_KIN_128_1}
Parity-projected wave-functions as a function of $r/A$ using the non-relativistic kinetic term $H_{K}^{(0)}$
(FDTD-based, left) and the non-relativistic kinetic term $H_{K}^{(2)}$ (FFT-based, right) using the increased lattice volume ($N = 128$).
The average radii in units of $A$ are given by}
\small
\begin{tabular}{r|c|c|c|c|c|c|c|c|c}
& $P_{\vec{p}_{k}}^{-} \Psi_{0}$ & $P^{+} \Psi_{0}$ & $P_{\vec{p}_{k}}^{-} \Psi_{1}$ & $P^{+} \Psi_{1}$ & $P_{\vec{p}_{k}}^{-} \Psi_{2}$ & $P^{+} \Psi_{2}$ & $\Psi_{1s}^{\text{th}}$ & $\Psi_{2s}^{\text{th}}$ & $\Psi_{2p}^{\text{th}}$ \\
\hline
$\langle r \rangle^{\Box}_{\text{left}}$ & $26.21$ & $19.28$ & $26.35$ & $11.61$ & $31.17$ & $27.47$ & \multirow{2}{*}{$8.22$} & \multirow{2}{*}{$30.77$} & \multirow{2}{*}{$26.18$} \\
\cline{1-7}
$\langle r \rangle^{\Box}_{\text{right}}$ & $26.22$ & $19.39$ & $26.28$ & $11.54$ & $41.80$ & $27.10$ &  &  & 
\end{tabular}
\end{figure}

\clearpage

\begin{figure}[!ht]
\centering
\ig{figure_KIN0_128_10000.png} \hfill
\ig{figure_KIN2_128_10000.png}
\caption{\label{fig_figure_comp_KIN_128_2}
Parity-projected wave-functions as in Fig.~\ref{fig_figure_comp_KIN_128_1} at a later snapshot time step.
The average radii in units of $A$ are given by}
\small
\begin{tabular}{r|c|c|c|c|c|c|c|c|c}
& $P_{\vec{p}_{k}}^{-} \Psi_{0}$ & $P^{+} \Psi_{0}$ & $P_{\vec{p}_{k}}^{-} \Psi_{1}$ & $P^{+} \Psi_{1}$ & $P_{\vec{p}_{k}}^{-} \Psi_{2}$ & $P^{+} \Psi_{2}$ & $\Psi_{1s}^{\text{th}}$ & $\Psi_{2s}^{\text{th}}$ & $\Psi_{2p}^{\text{th}}$ \\
\hline
$\langle r \rangle^{\Box}_{\text{left}}$ & $-$ & $8.50$ & $26.05$ & $30.77$ & $25.85$ & $30.71$ & \multirow{2}{*}{$8.22$} & \multirow{2}{*}{$30.77$} & \multirow{2}{*}{$26.18$} \\
\cline{1-7}
$\langle r \rangle^{\Box}_{\text{right}}$ & $-$ & $8.50$ & $26.16$ & $31.69$ & $26.14$ & $32.66$ &  &  & 
\end{tabular}
\end{figure}

\subsection{Quantumfdtd internal parity fixing machinery}
\label{subsec_ex_parity_fixing}

In this section, we investigate the ability to specify the initial spatial symmetry of the wave-function using the  \texttt{INITSYMMETRY} parameter (see section~\ref{subsec_icond}).
We set \texttt{INITSYMMETRY=1} for the left plots in all figures, which corresponds to symmetrization along the $z$ axis,
and we set \texttt{INITSYMMETRY=2} for the right plots in all figures, which corresponds to anti-symmetrization.
We again use $N=64$ and $M=0.3$~GeV and the non-relativist kinetic term $H_{K}^{(0)}$.
The dashed lines again represent the analytic solutions for the non-relativistic kinetic term with Coulomb potential.

The results should be compared with the figures of subsection~\ref{subsec_ex_nonrelativ} and with the FDTD based kinetic term, $H_{K}^{(0)}$, therein.
That is, with Figs.~\ref{fig_comp_KIN_0} and~\ref{fig_fi_comp_KIN_0}.

It can be seen that the choice \texttt{INITSYMMETRY=1} removes the $2p_z$ state completely.
However, there is still a residual negative parity component, most likely due to higher lying states or due to $2p$ states that are symmetric around the $z$ axis.
On the other hand, \texttt{INITSYMMETRY=2} removes all $s$ states, leaving only $2p$ and higher lying states. 
In particular, there are $3d$ states, which are odd under $P^-_{\vec{p}_k}$ for two directions while being even under parity.
These states are not separated by the post-processing parity projections (see section~\ref{subsec_parity_proy}).
Hence, these positive parity wave-functions mix contributions from states with different geometry that show as the broad bands in the Figures even for late time steps.
Note that, on the right energy plot of Fig.~\ref{fig_energ_comp_MIKE}, the ground state energy no longer corresponds to the ground state energy of the left plot.
It now corresponds to the first excited state energy of the left plot which is the $2p$ state, since it is the lowest lying state compatible with the anti-symmetrization.

The continuum average radii in units of $A$ are given by $\langle r \rangle^{\infty}_{1s} = 8.21$, $\langle r \rangle^{\infty}_{2s} = 32.83$, and $\langle r \rangle^{\infty}_{2p} = 27.36$.
As explained in the previous sections~\ref{subsec_ex_nonrelativ} and~\ref{subsec_ex_increase_volume}, the values are lower when computed using the box normalization: $\langle r \rangle^{\Box}_{1s} = 8.14$, $\langle r \rangle^{\Box}_{2s} = 20.30$, and $\langle r \rangle^{\Box}_{2p} = 19.06$, respectively.
At sufficiently large times, the average radii agree reasonably well with the theory expectations.

\begin{figure}[ht]
\centering
\ig{energ_KIN0_MIKE_PARITY.png} \hfill
\ig{energ_KIN0_MIKE_PARITY_ANTISYM.png}
\caption{\label{fig_energ_comp_MIKE}
Energies obtained as a function of $\tau$ using the non-relativistic kinetic terms $H_{K}^{(0)}$ with \texttt{INITSYMMETRY=1}
(left, symmetric along the $z$ axis) and \texttt{INITSYMMETRY=2} (right, antisymmetric along the $z$ axis).}
\end{figure}

\begin{figure}[ht]
\centering
\ig{proj_KIN0_MIKE_PARITY_0.png} \hfill
\ig{proj_KIN0_MIKE_PARITY_ANTISYM_0.png}
\caption{\label{fig_proj_comp_MIKE_0}
Weights of the parity projection of the ground state as a function of $\tau$ using the non-relativistic kinetic terms $H_{K}^{(0)}$
with \texttt{INITSYMMETRY=1} (left, symmetric along the $z$ axis) and \texttt{INITSYMMETRY=2} (right, antisymmetric along the $z$ axis).}
\end{figure}

\begin{figure}[ht]
\centering
\ig{proj_KIN0_MIKE_PARITY_1.png} \hfill
\ig{proj_KIN0_MIKE_PARITY_ANTISYM_1.png}
\caption{\label{fig_proj_comp_MIKE_1}
Weights of the parity projection of the first excited state as in Fig.~\ref{fig_proj_comp_MIKE_0}.}
\end{figure}

\begin{figure}[ht]
\centering
\ig{proj_KIN0_MIKE_PARITY_2.png} \hfill
\ig{proj_KIN0_MIKE_PARITY_ANTISYM_2.png}
\caption{\label{fig_proj_comp_MIKE_2}
Weights of the parity projection of the second excited state as in Fig.~\ref{fig_proj_comp_MIKE_0}.}
\end{figure}

\begin{figure}[ht]
\centering
\ig{figure_KIN0_MIKE_PARITY_1500.png} \hfill
\ig{figure_KIN0_MIKE_PARITY_ANTISYM_1500.png}
\caption{\label{fig_comp_MIKE_1}
Parity-projected wave-functions as a function of $r/A$ using the non-relativistic kinetic term $H_{K}^{(0)}$
with \texttt{INITSYMMETRY=1} (left, symmetric along the $z$ axis) and \texttt{INITSYMMETRY=2} (right, antisymmetric along the $z$ axis).
The average radii in units of $A$ are given by}
\small
\begin{tabular}{r|c|c|c|c|c|c|c|c|c}
& $P_{\vec{p}_{k}}^{-} \Psi_{0}$ & $P^{+} \Psi_{0}$ & $P_{\vec{p}_{k}}^{-} \Psi_{1}$ & $P^{+} \Psi_{1}$ & $P_{\vec{p}_{k}}^{-} \Psi_{2}$ & $P^{+} \Psi_{2}$ & $\Psi_{1s}^{\text{th}}$ & $\Psi_{2s}^{\text{th}}$ & $\Psi_{2p}^{\text{th}}$ \\
\hline
$\langle r \rangle^{\Box}_{\text{left}}$ & $-$ & $12.44$ & $-$ & $11.66$ & $-$ & $21.71$ & \multirow{2}{*}{$8.14$} & \multirow{2}{*}{$20.30$} & $-$ \\
\cline{1-7}\cline{10-10}
$\langle r \rangle^{\Box}_{\text{right}}$ & $18.64$ & $22.09$ & $15.09$ & $21.92$ & $17.48$ & $23.76$ &  &  & 19.06
\end{tabular}
\end{figure}

\clearpage

\begin{figure}[!ht]
\centering
\ig{figure_KIN0_MIKE_PARITY_10000.png} \hfill
\ig{figure_KIN0_MIKE_PARITY_ANTISYM_10000.png}
\caption{\label{fig_comp_MIKE_2}
Parity-projected wave-functions as in Fig.~\ref{fig_comp_MIKE_1} at a later snapshot time step.
The average radii in units of $A$ are given by}
\small
\begin{tabular}{r|c|c|c|c|c|c|c|c|c}
& $P_{\vec{p}_{k}}^{-} \Psi_{0}$ & $P^{+} \Psi_{0}$ & $P_{\vec{p}_{k}}^{-} \Psi_{1}$ & $P^{+} \Psi_{1}$ & $P_{\vec{p}_{k}}^{-} \Psi_{2}$ & $P^{+} \Psi_{2}$ & $\Psi_{1s}^{\text{th}}$ & $\Psi_{2s}^{\text{th}}$ & $\Psi_{2p}^{\text{th}}$ \\
\hline
$\langle r \rangle^{\Box}_{\text{left}}$ & $-$ & $8.38$ & $-$ & $18.90$ & $-$ & $18.90$ & \multirow{2}{*}{$8.14$} & \multirow{2}{*}{$20.30$} & $-$ \\
\cline{1-7}\cline{10-10}
$\langle r \rangle^{\Box}_{\text{right}}$ & $17.63$ & $-$& $18.28$ & $22.13$ & $19.64$ & $22.01$ &  &  & 19.06
\end{tabular}
\end{figure}

\subsection{Convergence conditions}
\label{subsec_conv_conditions}

The convergence of the wave-functions shown in Figs.~\ref{fig_energ_KIN3} to~\ref{fig_comp_MIKE_2} is pretty clean.
Discretization errors that are predominantly due to the lattice implementation of the Coulomb potential are visible on the coarser lattices and excited states are affected by finite volume effects in small lattice volumes.
The wave-functions generated by the \texttt{quantumfdtd} can be considered reliable, only after the following conditions are met:
\begin{itemize}
\item Appropriate parity projection and normalization has been performed using the provided post-processing scripts for separating the wave-functions (see section~\ref{subsec_parity_proy}).
\item Using the snapshot facility of the code for multiple values of $\tau$ the wave-function has been shown to be stable within small changes of $\tau$ and the energy $E_{s,i}(\tau)$ should be in a plateau, i.e., flat or, at least, in a saddle point.
\item $A < M^{-1}$, such that the typical radius of the ground state is substantially larger than the lattice spacing $A$, $\langle r \rangle^{\Box}_{P^+\Psi_0} \gg A$. This condition is violated in section~\ref{subsec_ex_relativ_vs_nonrelativ_high_mass} in the example.
\item $2\langle r \rangle^{\Box}_{P^{\pm}\Psi_{1,2}} \ll N\cdot A$, such that typical radii of the excited states are substantially smaller than half of the box size $N\cdot A$.
\end{itemize}

\clearpage

\section{Compilation and running}
\label{sec_compile}

We require $\text{\texttt{\$\{NCORES\}}}$ MPI processes, with $\text{\texttt{\$\{NCORES\}}}$ being a divisor of $N$.
The following external libraries are used:
\begin{itemize}
\item MPI library
\item \texttt{FFTW\_MPI}, version 3 \footnote{
\href{http://fftw.org/}{http://fftw.org/} and Ref.~\cite{fftw3}}
\item GNU Scientific Library (GSL) \footnote{
\href{https://www.gnu.org/software/gsl/}{https://www.gnu.org/software/gsl/} and Ref.~\cite{gsl}},
with the CBLAS link \footnote{
\href{https://www.netlib.org/blas/}{https://www.netlib.org/blas/} and Refs.~\cite{blas1,blas2}}
\item For some of the post-processing scripts, Python~3 is required\footnote{
\href{https://www.python.org/downloads/}{https://www.python.org/downloads/}}
\end{itemize}
There are both, open source (Open MPI\footnote{
\href{https://www.open-mpi.org/}{https://www.open-mpi.org/}})
and proprietary implementations of the MPI library.
The open source Open~MPI~2.1.1 and the proprietary
IBM MPI~1.4 versions of the MPI library have been used
for testing purposes.
Concerning \texttt{FFTW\_MPI} and GSL, both are open source.
Versions \texttt{FFTW\_MPI}~3.5.7 and GSL~2.5 have been used
for testing.
Hence, the usage of this code does not require any proprietary software.

The directory structure is as follow:
\begin{description}
\item[\texttt{build/src}:] Source code.
\item[\texttt{build/include}:] Header files.
\item[\texttt{build/obj}:] Folder for intermediate object files.
\item[\texttt{data}:] Folder for default output files.
\item[\texttt{input}:] Folder for default input files. The file \texttt{input/params.txt} is necessary for running the program.
\item[\texttt{license}:] GNU Public License.
\item[\texttt{archive}:] Several versions of the \texttt{Makefile}, for different machine configurations (architecture).
\item[\texttt{scripts}:] Examples of scripts for running the program and post-processing of the output data.
This includes the Python3~\cite{python} module \texttt{quantumfdtd.py}, that offers a Python class (\texttt{quantumfdtd}) for the post-analysis inside the Pandas~\cite{pandas, reback2020pandas, mckinney-proc-scipy-2010} framework as well as the scripts \texttt{symmetrize.sh}, for computing parity projections, and several scripts for extracting plots\footnote{%
\texttt{postprocessing.sh} for the potential, energies and wave-functions and \texttt{postprocessing\_\allowbreak{}sym.sh} for the parity projections.
In both cases, plots of the full 3D wave-functions snapshots are obtained if they have been computed.}.
\end{description} 

The root directory contains a \texttt{Makefile} that allows one to build
the program provided that all of the external libraries are installed.
If the libraries are installed in non-standard locations,
this \texttt{Makefile} should be modified accordingly.

The following commands, executed in the root folder of the provided tarball,
are required for building and executing the program:
\begin{verbatim}
$ make
$ mpirun -N ${NCORES} ./mpisolve
\end{verbatim}
Here, \texttt{\$\{NCORES\}} is the number of cores to be used which should be
a divisor of the number of the spatial grid points.\\[1em]
\noindent
The source code is available on GitHub: \url{https://github.com/quantumfdtd/quantumfdtd_v3}

\section{Parameters of the program}
\label{sec_params}

The parameters listed in the following subsections should be defined in the standard parameter file located at \texttt{input/params.txt}.
The \texttt{input} directory should be located in the directory where the program is also run from.
The configuration file accepts full line comments, that are marked with an initial double backslash, \texttt{\textbackslash\textbackslash}, on the commented line.

For convenience, the program accepts command line parameters as well.
They can be parsed as follows:
\begin{verbatim}
$ mpirun -N ${NCORES} ./mpisolve -KEY1 ARG1 -KEY2 ARG2 ...
\end{verbatim}
Unless stated differently, all of the numerical values are separated by tabs in all of the output files.

\subsection{Basic configuration}
\label{subsec_basic_params}

This program is based on an iterative method for finding both the ground state energy and the two first excited state energies.
The corresponding wave-functions are computed as well.
Hence, most of the basic configuration flags that are required for running the program
are related both to the setup of the lattice (discretization) and the termination conditions of the iterative procedure.
The following parameters are required for running the program:
\begin{description}
\item[\texttt{NUM}:] Number of spatial grid points.
Note that they should be divisible by the number of computational nodes, i.e., $\text{\texttt{NUM}} \mod N_{\text{MPI cores}} = 0$.
\item[\texttt{A}:] Spatial grid spacing. Default units: 1/GeV.
\item[\texttt{MASS}:] Reduced mass of the simulated system. Default units: GeV.
\item[\texttt{KINTERM}:] Selection of the kinetic term to be used (see subsection~\ref{subsec_kin}):
\begin{itemize}
\item The non-relativistic kinetic term based on finite difference time domain (\texttt{0}).
\item The relativistic kinetic term (\texttt{3}).
\item 2 other non-relativistic kinetic terms based on FFT (\texttt{1},\texttt{2}).
\end{itemize}
\item[\texttt{POTENTIAL}:] Potential to simulate.
There are 6 hard-coded potentials described in subsection~\ref{subsec_hpot} and 4 different options to read the potential from file described in subsection~\ref{subsec_epot}.
\item[\texttt{EPS}:] Temporal grid spacing. Should be $\leq \text{\texttt{A}}^{2}/3$.
\item[\texttt{UPDATE}:]\label{explan_bare} How many steps should be taken before checking the \texttt{SNAP\-UPDATE}, \texttt{SNAP\-DUMP} and \texttt{STEPS} variables, and computing observables\footnote{%
Ground state energy and the expectation values $\braket{r^{2}}$, $\braket{x}$, $\braket{y}$, and $\braket{z}$.
Note that here $\braket{r^{2}}$ is w.r.t. the center of lattice, i.e., w.r.t. the point $(x_{1},x_{2},x_{3}) = (N/2,N/2,N/2)$.}
and printing them to command line.\\
If the \texttt{SAVEDECAY} flag is enabled, observables are also saved to the file \emph{decay.dat}, as an ASCII table: $T$, $E_{b}$.
Here, $T$ stands for the total time (\texttt{EPS}$\times$\texttt{\#steps}) and $E_{b}$ is the bare ground state energy.\\
Note that some hard-coded potentials split the potential in two parts\footnote{%
There is an option for using this splitting for externally provided potentials as well, see sec.~\ref{subsec_epot}.}:
A \emph{bare} ($V_{b}$) and a \emph{subtracted} ($V_{s}$) potential.
$V_{b}$ is the one used inside the iterative procedure, for computing the wave-functions $\Psi_i$ and energies $E_{b}$.
Afterwards, the \emph{subtracted} value of the energy, $E_{s}$, is computed via
\begin{equation}
E_{s,i} = E_{b,i} - \int_{V} |\Psi(r)|^{2} V_{s}(r) dV \,,
\end{equation}
where the integral extends over all of the lattice positions.
This is correct for the potential $V(r)$ in the particular case $V(r) = V_{b}(r) - V_{s}(r)$, where $V_{s}(r) = V_{s}$, a scalar constant.
\item[\texttt{SNAPUPDATE}:]\label{flag_snapupdate} How many steps should be taken before enforcing both
normalization of the ground state wave-function ($\int_{V} dV |\psi_{0}|^{2} = 1$) and its required symmetries (see the flag \texttt{INITCONDTYPE}).
It triggers the check of the tolerance condition (see \texttt{TOLERANCE} flag). It also triggers
the storage of the 2 snapshots that are required for the extraction of the first and second excited states
via overlapping (see section~\ref{subsec_iter}).
\item[\texttt{SNAPDUMP}:] In case \texttt{DUMPSNAPS} (see below) is set to a non-null value,
how many steps should be taken before storing a wave-function snapshots to disk.
\item[\texttt{TOLERANCE}:] Convergence tolerance w.r.t. the ground state energy.
The stop condition being\footnote{%
$E_{b,i} - E_{b,i-1}$ stands for the difference between the current bare ground state energy and the one obtained when \texttt{SNAPUPDATE} was triggered previously.} $E_{b,i} - E_{b,i-1} < \texttt{TOLERANCE}$. It can be disabled by setting a negative value.
\item[\texttt{STEPS}:] Maximum number of steps to take.
This limit is enforced each time a number of steps given by the variable \texttt{UPDATE} is evaluated.
Hence, the actual number of steps to take is exactly \texttt{STEPS} only if \texttt{STEPS} is divisible by \texttt{UPDATE}.
Otherwise, it will be the smallest multiple of \texttt{UPDATE} bigger than \texttt{STEPS}. This is the actual number of steps
to take if the \texttt{TOLERANCE} parameter is set to a negative value.
\item[\texttt{SAVEWAVEFNCS}:] Save a full snapshot of the 3d wave-functions at the end of the run.
Set to \texttt{1} for activation.
\item[\texttt{SAVEDECAY}:] Save a table of energy vs. time.
This is particularly useful for checking the convergence of the iterative procedure to the ground state, especially, when using the FFT-based kinetic terms.
\item[\texttt{SAVEPOT}:] Name of the external file where the potential that is used will be saved to.
The file format is an ASCII list containing:
\begin{equation*}
x_{1}/A \quad x_{2}/A \quad x_{3}/A \quad r^2/A^2 \quad \mathop{\text{Re}} V \quad \mathop{\text{Im}} V \,,
\end{equation*}
where $r^2/A^2=(x_{1}/A)^2 + (x_{2}/A)^2 + (x_{3}/A)^2$, centered at the point $(x_{1},x_{2},x_{3})=(0,0,0)$.
\item[\texttt{DATAFOLD}:] Name of the folder where the full output, wave-functions, and snapshots will be saved to.
By default, is set to \texttt{data}.
The file format is an ASCII list containing:
\begin{equation*}
x_{1}/A \quad x_{2}/A \quad x_{3}/A \quad r^2/A^2 \quad \mathop{\text{Re}}\psi_i \quad \mathop{\text{Im}}\psi_i \,,
\end{equation*}
where $r^2/A^2=(x_{1}/A)^2 + (x_{2}/A)^2 + (x_{3}/A)^2$, centered at the point $(x_{1},x_{2},x_{3})=(0,0,0)$, and $i=0$ (ground state), $1$ (first excited state), or $2$ (second excited state).
\item[\texttt{DUMPSNAPS}:] Dump debugging files containing either "snapshot slices" or full 3D wave-functions:
\begin{itemize}
\item Disables snapshot recording (\texttt{0}).
\item Store "snapshot slices" of the ground state (\texttt{1}). Only valid for potentials centered on $(N/2,N/2,N/2)$. Do not use when centering the potential on $(0,0,0)$.
\item Store full 3D wave-functions of the ground state (\texttt{2}).
\item Store "snapshot slices" of the ground state, and of the first and second excited states (\texttt{3}). Only valid for potentials centered on $(N/2,N/2,N/2)$. Do not use when centering the potential on $(0,0,0)$.
\item Stores full 3D wave-functions of the groundstate, and of the first and second excited states (\texttt{4}).
\end{itemize}
The full 3D wave-functions are stored in the same way as the ones triggered by \texttt{DATAFOLD}.\\
The file format for the "slices" consists of three concatenated ASCII lists, separated by a row containing "\&\&", respectively.
These lists are, respectively:
\begin{align*}
x_{1}/A && \mathop{\text{Re}} V[x_{1}, N/(2A), N/(2A)] && \mathop{\text{Im}} V[x_{1}, N/(2A), N/(2A)] \\
x_{2}/A && \mathop{\text{Re}} V[N/(2A), x_{2}, N/(2A)] && \mathop{\text{Im}} V[N/(2A), x_{2}, N/(2A)] \\
x_{3}/A && \mathop{\text{Re}} V[N/(2A), N/(2A), x_{3}] && \mathop{\text{Im}} V[N/(2A), N/(2A), x_{3}]
\end{align*}
In both cases, the output is stored in the folder \texttt{data/snapshot} and such a folder must exist prior to running the program.
The file names is \texttt{data/\allowbreak{}snapshot\_\allowbreak{}wave-function\_\allowbreak{}step\_\allowbreak{}i.dat}, where \texttt{step} is the step associated to the snapshot ($\tau=\text{EPS}\cdot\text{step}$) and \texttt{i} is an index running over the computational nodes.
\end{description}


\subsection{Initial conditions}
\label{subsec_icond}

The following parameters allow changing the initial state condition. The Coulomb-type initial conditions
are intended to be a sensible initial guess for spherically symmetric confining central potentials.

\begin{description}
\item[\texttt{INITCONDTYPE}:] Initial condition to use.
Valid values:
\begin{description}
\item[\texttt{0}:] Read initial condition from \texttt{wave-function\_\allowbreak{}0\_\allowbreak{}\#.dat} files in the \texttt{data} directory.
\item[\texttt{1}:] Random Gaussian noise with standard deviation given by the \texttt{SIG} parameter.
\item[\texttt{2}:] Non-relativistic Coulomb-type initial condition, centered at the same point where the potential is centered.
I.e., either in the center of the lattice volume at the point $(N/2,N/2,N/2)$ (this applies to all of the hard-coded potentials with code <\texttt{100}), or in the origin of the lattice at the point $(x_{1},x_{2},x_{3}) = (0,0,0)$, otherwise.
By default, we introduce a linear combination of the states $1s$, $2s$, $2p(m=0)$ and the real part of the $2p(l=1)$.
Note that the selection of the $z$ axis of the $2p$ states is chosen through the flag \texttt{INITCONDAXIS}.
\item[\texttt{3}:] Constant wave-function with value $0.1$.
\item[\texttt{4}:] Boolean test grid, $V = (x_{1} \mod 2) \cdot (x_{2} \mod 2) \cdot (x_{3} \mod 2)$.
\item[\texttt{5}:] This was an intermediate step during the development of the new version of \texttt{quantumfdtd}.
The idea was to define an initial wave-function containing a sensible selection of Fourier components.
We ended up with the function shown in Listing~\ref{listing} on page~\pageref{listing}.
If the need arises to modify this function in order to include own guesses for particular potentials, we advise looking into the function\newline
\texttt{void setInitialConditions(int seedMult)}, contained in the file \texttt{build/\allowbreak{}src/\allowbreak{}initialconditions.cpp}.
\end{description}
\item[\texttt{INITCONDAXIS}:] For \texttt{INITCONDTYPE} set to \texttt{2}, this flag sets the $z$ axis of the $2p_z$ components of the wave-function.
For the choices \texttt{0}, \texttt{1}, or \texttt{2} the axis is set to $x_{1}$, $x_{2}$, or $x_{3}$, respectively.
\item[\texttt{INITSYMMETRY}:] Symmetrize the wave-function.
Valid values are:
\begin{description}
\item[\texttt{0}:] Disable option.
\item[\texttt{1}:] Symmetric w.r.t. the $x_{3}$ direction.
\item[\texttt{2}:] Antisymmetric w.r.t. the $x_{3}$ direction.
\item[\texttt{3}:] Symmetric w.r.t. the $x_{2}$ direction.
\item[\texttt{4}:] Antisymmetric w.r.t. the $x_{2}$ direction.
\end{description}
Note that because of the internal memory layout, it is not possible to offer options in the $x_{1}$ direction without a high MPI overhead.
\item[\texttt{SIG}:] Standard deviation - used when \texttt{INITCONDTYPE} is set to \texttt{1}.
\end{description}

\clearpage

{
\captionsetup[table]{name=Listing}
\begin{table}[p]
\centering
\begin{minipage}[t]{\textwidth}
\centering
\fbox{\lstinputlisting{initcondtype_5.cpp}}
\end{minipage}
\caption{\label{listing}
Implementation of \texttt{INITCONDTYPE=5}: a selection of Fourier components as initial guess.}
\end{table}
}

\clearpage

\subsection{Hard-coded potentials}
\label{subsec_hpot}

The program comes bundled with several hard-coded potentials, that were already present in the original code~\cite{Strickland:2009ft, Dumitru:2009ni, Dumitru:2009fy}.
The user should be aware that some of these potentials may be unphysical depending on the boundary conditions.
\begin{description}
\item[\texttt{0}:] No potential ($V(\vec{r}) = 0$).
This is equivalent to an infinitely deep 3d well due to boundary conditions - provided that $\text{\texttt{KINTERM}} = 0$ is used.
\item[\texttt{1}:] 3d square well in the center of the lattice volume:
\begin{equation}
V(x_{1},x_{2},x_{3}) = -10 \, \Theta(L/4 - |x_{1}|) \Theta(L/4 - |x_{2}|) \Theta(L/4 - |x_{3}|) \,,
\end{equation}
where $\Theta(x)$ is the Heaviside Theta function.
\item[\texttt{2}:] Radial Coulomb potential:
\begin{equation}
V(x_{1},x_{2},x_{3}) =
\begin{cases}
-\frac{1}{r} \,, & \text{for } r > \texttt{A} \,, \\
-\frac{1}{\texttt{A}} \,, & \text{for } r \leq \texttt{A} \,,
\end{cases}
\end{equation}
where $r = \sqrt{x_{1}^2 + x_{2}^2 + x_{3}^2}$.
\item[\texttt{3}:] Elliptical Coulomb potential:
\begin{equation}
V(x_{1},x_{2},x_{3}) =
\begin{cases}
-\frac{1}{\tilde{r}}  \,, & \text{for } \tilde{r} >  \texttt{A} \,, \\
-\frac{1}{\texttt{A}} \,, & \text{for } \tilde{r} \leq \texttt{A} \,,
\end{cases}
\end{equation}
where $\tilde{r} = \sqrt{x_{1}^2 + x_{2}^2 + 4 x_{3}^2}$.
\item[\texttt{4}:] 3d harmonic oscillator:
\begin{equation}
V(x_{1},x_{2},x_{3}) = -\frac{1}{2} (x_{1}^2 + x_{2}^2 + x_{3}^2) \,.
\end{equation}
\item[\texttt{5}:] Complex 3d harmonic oscillator:
\begin{equation}
V(x_{1},x_{2},x_{3}) = -\frac{1}{2} (1 + \text{i}) (x_{1}^2 + x_{2}^2 + x_{3}^2) \,.
\end{equation}
\item[\texttt{6}:] Cornell potential with string breaking:\\
We take the string-breaking scale to be $r_{\text{sb}} = 5.5745$~GeV$^{-1}$ and define
\begin{equation}
\hat{r} =
\begin{cases}
r \,, & \text{for } r < r_{\text{sb}} \,, \\
r_{\text{sb}} \,, & \text{for } r \geq r_{\text{sb}} \,,
\end{cases}
\end{equation}
where $r = \sqrt{x_{1}^2 + x_{2}^2 + x_{3}^2}$, and
\begin{equation}
V(x_{1},x_{2},x_{3}) = -\frac{0.385}{\hat{r}} + \texttt{SIGMA} \, \hat{r} + 4 \, \texttt{MASS} \,,
\end{equation}
where \texttt{MASS} is the \texttt{quantumfdtd} reduced mass and \texttt{SIGMA} is an independent parameter.
\item[>\texttt{99}:] If $100$ is added to any of the previous potential codes, the corresponding potential will be loaded with its origin centered at the point $(x_{1},x_{2},x_{3}) = (0,0,0)$ instead of at the point $(x_{1},x_{2},x_{3}) = (N/2,N/2,N/2)$.
This is especially useful for periodic boundary conditions.
\end{description}


\subsection{External potentials}
\label{subsec_epot}
The following values of \texttt{POTENTIAL} offer different options to read external potentials from files:
\begin{description}
\item[\texttt{90}:] Read external potential as a table containing
\begin{equation*}
x_{1} \text{\textvisiblespace} x_{2} \text{\textvisiblespace} x_{3} \text{\textvisiblespace} \Re(V) \text{\textvisiblespace} \Im(V) \,,
\end{equation*}
separated by spaces.
\item[\texttt{91}:] Read external potential as a table containing
\begin{equation*}
R^{2} \text{\textvisiblespace} \Re(V) \text{\textvisiblespace} \Im(V) \,,
\end{equation*}
separated by spaces.
\item[>\texttt{189}:] As in the previous section~\ref{subsec_hpot}, \texttt{190} and \texttt{191} are, respectively, the same that \texttt{90} and \texttt{91} but centered on $(0,0,0)$.
\end{description}

The following configuration flags are only valid when loading potentials from external files\footnote{Values of the \texttt{POTENTIAL} parameter: \texttt{90}, \texttt{91}, \texttt{190} and \texttt{191}}:
\begin{description}
\item[\texttt{EXPOT}:] Name of the file containing the external potential.
The file format depends on the selected kind of external potential as described above.
\item[\texttt{POTCRITR}:] Defines a point from which on (in units of the lattice spacing), the potential will be a $\chi^{2}$ adjustment to the following function:
\begin{equation}
V_{A,B,\sigma}(r) = A + B/r + \sigma\cdot r
\end{equation}
The parameters $A$, $B$, and $\sigma$ are written to console. Only valid for values of the \texttt{POTENTIAL} parameter: \texttt{90} and \texttt{190}.
\item[\texttt{POTFLATR}:] Defines a point from which on (in units of the lattice spacing), the potential is flat taking the value of the point previous to it. Only valid for \texttt{POTENTIAL} values \texttt{90} and \texttt{190}.
\end{description}

\subsection{Post-processing scripts}
\label{subsec_post_scripts}

The following \texttt{bash} and \texttt{Python~3} scripts have been included for the post-processing of the wave-functions.
They can be found in the \texttt{scripts} folder.
Unless otherwise stated, they should be run from the parent folder where both the \texttt{data} and \texttt{input} folders are located:
\begin{description}
\item[\texttt{symmetrize.sh}:] Computes the projections and normalizations of all of the wave-functions,
using the scripts \texttt{symmetrize\_\allowbreak{}wf.py} and \texttt{normalize.py}.
It accepts direct output stored in the \texttt{data} folder or it accepts wave-functions
which are stored in \texttt{gz} compressed \texttt{tar} archives named \texttt{wave-function.tgz} which allows for batch evaluation of results from many runs.
Non-normalized output files named \texttt{wave-function\_\allowbreak{}\%i\_\allowbreak{}\%p.dat}
and normalized ones named \texttt{wave-function\_\allowbreak{}\%i\_\allowbreak{}\%p\_\allowbreak{}norm.dat}
are generated. \texttt{\%i}=0,1,2 denotes the ground, first and second excited states, respectively.
\texttt{\%p} stands for the projection of the wave-function, as follows:
\begin{description}
\item[\texttt{wave-function\_\allowbreak{}\%i\_\allowbreak{}all.dat.gz}:]
The wave-function as computed by \texttt{quantumfdtd}, but in a single file instead of one file per computational node.
\item[\texttt{wave-function\_\allowbreak{}\%i\_\allowbreak{}p.dat.gz}:] Positive parity projection $P^+$ (Eq.~\ref{eq_parity}).
\item[\texttt{wave-function\_\allowbreak{}\%i\_\allowbreak{}m.dat.gz}:] Negative parity projection $P^-$ (Eq.~\ref{eq_parity}).
\item[\texttt{wave-function\_\allowbreak{}\%i\_\allowbreak{}pk.dat.gz}:] Negative parity projection around-an-axis $P^-_{\vec{p}_k}$ (Eq.~\ref{eq_parity_axis}).
\item[\texttt{wave-function\_\allowbreak{}\%i\_\allowbreak{}all\_\allowbreak{}norm.dat.gz}:] Full normalized wave-function.
\item[\texttt{wave-function\_\allowbreak{}\%i\_\allowbreak{}p\_\allowbreak{}norm.dat.gz}:] Normalized $P^+$ projection.
\item[\texttt{wave-function\_\allowbreak{}\%i\_\allowbreak{}m\_\allowbreak{}norm.dat.gz}:] Normalized $P^-$ projection.
\item[\texttt{wave-function\_\allowbreak{}\%i\_\allowbreak{}pk\_\allowbreak{}norm.dat.gz}:] Normalized $P^-_{\vec{p}_k}$ projection.
\item[\texttt{wave-function\_\allowbreak{}\%i\_\allowbreak{}proj.dat}:] Weights of the parity projections.
\end{description}
For all of the cases but the normalized projection $P^-_{\vec{p}_k}$, the output is a gzipped ASCII file containing:
\begin{equation*}
x_{1}/A \quad x_{2}/A \quad x_{3}/A \quad r^2/A^2 \quad \mathop{\text{Re}}\Psi_i \quad \mathop{\text{Im}}\Psi_i \,,
\end{equation*}
For the normalized projection $P^-_{\vec{p}_k}$, the output follows the format:
\begin{equation*}
x_{1}/A \quad x_{2}/A \quad x_{3}/A \quad r/A \quad \cos\theta \quad \mathop{\text{Re}}\Psi_i \quad \mathop{\text{Im}}\Psi_i \,,
\end{equation*}
\item[\texttt{postprocess.sh}:] Generates plots of the potential and the wave-functions.
\item[\texttt{postprocess\_\allowbreak{}sym.sh}:] Generates plots of the parity projected wave-functions obtained with \texttt{symmetrize.sh}.
\item[\texttt{allpostprocess.sh}:] Runs \texttt{postprocess.sh}, \texttt{symmetrize.sh} and \texttt{postprocess\_\allowbreak{}sym.sh}, sequentially.
\item[\texttt{cleandatafiles.sh}:] Runs \texttt{rm -rf data/*.* \allowbreak{}data/snapshot/*.* \allowbreak{}debug/*.*} for cleaning up a previous run.
\item[\texttt{symmetrize\_\allowbreak{}wf.py}:] Computes the positive and negative parity projections of the given wave-function.
It is called as:\\
\texttt{../scripts/symmetrize\_\allowbreak{}wf.py [params.txt] wave-function\_\allowbreak{}\%i\_\allowbreak{}\%d.dat}\\
The optional parameter \texttt{[params.txt]} encode the path to the config file,
from which the simulation parameters like $N$, $A$, or the centering of the lattice are taken. If not provided,
the script expects to be called from the \texttt{data} folder and tries to load the config file from the relative path \texttt{../input/params.txt}.\\
\texttt{\%i} should be substituted by the numbers 0, 1 and 2 for the ground, first and second excited states, respectively.\\
\texttt{\%d} should be left on the actual call. It means that the script will look for files where \texttt{\%d} is a numerical index starting on 0.
This procedure has been implemented because \texttt{quantumfdtd} splits the output data between as many files as computational nodes are used.
This \texttt{\%d} is the numerical index of these files.
In practice, for the ground state wave-function, the call (from the \texttt{data} folder) is\\
\texttt{../scripts/symmetrize\_\allowbreak{}wf.py wave-function\_\allowbreak{}0\_\allowbreak{}\%d.dat}
\item[\texttt{normalize.py}:] Normalizes the wave-funtions to unity over the lattice volume.
It is expected to be run from inside the \texttt{data} folder, but it can be run from any folder where wave-functions are stored.
It is called as:\\
\texttt{../scripts/symmetrize\_\allowbreak{}wf.py wave-function.dat out.dat}.
\end{description}

\subsection{Post-processing Python module}
\label{subsec_post_python}

We provide a post-processing Python~3 module in \texttt{scripts/quantumfdtd.py},
whose class \texttt{quantumfdtd} aims at helping with the analysis via Python~3~\cite{python}, Pandas~\cite{pandas, reback2020pandas, mckinney-proc-scipy-2010}, NumPy~\cite{numpy, oliphant2006guide, 5725236}, and Matplotlib~\cite{matplotlib, 4160265}.
It accepts as input the output of either \texttt{symmetrize.sh} or \texttt{allpostprocess.sh} from the previous section~\ref{subsec_post_scripts}.
The class functions are:
\begin{description}
\item[\texttt{\_\_init\_\_(self, case, config=None)}:] Constructor. It requires a \texttt{case} parameter that is a parent folder where \texttt{input}
and \texttt{data} folders of an individual run of \texttt{quantumfdtd} and \texttt{allpostprocess.sh} are stored.
There is an optional parameter, \texttt{config}, to provide an absolute or relative path to a \texttt{params.txt} (config file).
The function \texttt{load\_\allowbreak{}energies}, that loads the energies of the final wave-function
(and the intermediate snapshots, if they are taken), is also called by the constructor.
\item[\texttt{load\_case(self, case, config=None)}:] In case the default constructor is called, this function accepts the same arguments as the constructor in order to perform a proper initialization of the class.
It also allows to load a different \emph{case}\footnote{%
Path where a run of \texttt{quantumfdtd} is stored.} without destroying and creating a new class.
\item[\texttt{load\_energies(self)}:] This function is implicitly called by the constructor and by the \texttt{load\_\allowbreak{}case} function, so that it does not need to be called explicitly.
It loads the eigen-energies as a function of $\tau$ into the \texttt{self.e\_wf} variable, as a Pandas \texttt{DataFrame}.
It also loads the \texttt{decay.dat} table, that encodes the eigen-energy of the ground state with a higher frequency on $\tau$, and returns it as a \texttt{DataFrame} via the \texttt{self.e\_dec} variable.
\item[\texttt{load\_potential(self)}:] This function loads the potential and returns it as a \texttt{DataFrame} with the columns $x$, $y$, $z$ (integers, coordinates in lattice space in units of $A$), $r$ (in units of $A$), \texttt{re\_v} and \texttt{im\_v} (real and imaginary part of the potential, in GeV).
\item[\texttt{load\_wf(self, snap, state, figure, min\_sep\_edge=-1., normalization=True)}:]\hfill\\
This function loads a wave function and returns it as a \texttt{DataFrame}.
If \texttt{snap=-1}, the final wave-function is loaded.
Otherwise, an intermediate 3D snapshot is loaded, whose \texttt{STEP} value is given by the \texttt{snap} parameter.
Remember that $\tau = \text{\texttt{STEP}}\cdot\text{\texttt{EPS}}$.
The argument \texttt{state} is an integer value from 0 to 2, corresponding to the ground state, or to the first or second excited state.
The \texttt{figure} variable is \texttt{m}, \texttt{p}, \texttt{pk} or \texttt{all}, corresponding to the projection that is going to be loaded: negative parity projection $P^{-}$, positive parity projection $P^{+}$, negative-around-an-axis $P^{-}_{p_{k}}$ and non-projected wave-function, respectively.
\item[\texttt{fix\_positive\_wf(self,wf,lim=0.05)}:] This function returns the wave-function, inverting its sign if the integration of its real values for the points $r < A \cdot N \cdot \text{\texttt{lim}}$ is negative.
\item[\texttt{get\_avg\_r(self,wf)}:] This function returns the average radius, $\langle r \rangle$, of the wave-function passed by the \texttt{wf} parameter.
\item[\texttt{load\_proj\_components(self)}:] This function adds additional rows to \texttt{self.e\_wf}.
In these rows, the weights of the parity projections of each excited state are stored.
If snapshots are taken, this option is very useful for measuring the decay of the different excited states through the iterative procedure explained in section~\ref{subsec_iter}.
\item[\texttt{get\_snaps(self)}:] This function returns an \texttt{array} of snapshots.
For each snapshot, the value of \texttt{STEP} is stored on the array, where $\tau = \text{\texttt{STEP}}\cdot\text{\texttt{EPS}}$.
\end{description}

\section{Conclusions}
\label{sec_conclus}

In this work, the \texttt{quantumfdtd} code version~2.1~\cite{Strickland:2009ft,Margotta:2011ta}
has been extended to include the case of a relativistic Schr{\"o}dinger equation.
This is accomplished via a FFT implementation for the RQM inspired kinetic Hamiltonian and working with the iterative solver in imaginary time.
Additionally, we have extended the program to have support for arbitrary external potentials, introduced as ASCII tables, and included several analysis scripts.
One of them allows to extract parity projections of the extracted wave-functions.
This greatly improves the flexibility and usability of the original program
and makes this version~3 valuable even if the original finite-difference time-domain (FDTD) solver
(for the non-relativistic Schr{\"o}dinger equation) is being used.

We have seen in section~\ref{sec_examples} that the ground-state wave-function
is stable for high values of the time evolution parameter $\tau$.
The output wave-functions are linear combinations of the first excited states.
The overlap method is an approximation to increase the weight of the first
and second excited state in the corresponding wave-functions.
But the evolution Eq.~\ref{iter_decay}, and a periodic internal normalization to one of the evolved wave-function,
makes the weight of the excited states in the wave-function decay exponentially.
These weights will eventually decay below the machine precision,
eliminating the information about excited states from the runs with the highest values of $\tau$.
The parity projection scripts are a useful tool for separating excited states
when they are nearly degenerated in energy which is especially noticeable for the $2s$ and the $2p$-states.

To sum up: the ground-state wave-function should be easy to recover in all of the situations,
provided that the corresponding energy is in a plateau and the wave-function is stable w.r.t. $\tau$.
For increased precision, the positive parity projection can be used.
The $2p$-state can be separated by means of the negative-around-an-axis parity projection operator.
The $2s$-state can be also separated, but is more challenging due to the behavior of the tail after the first node on a periodic potential.
A large lattice volume can help here.
In the examples presented in this paper, we focused on relatively small lattices ($N \leq 128$) {that are most relevant for interfacing with lattice gauge theory simulations} in order to assess the numerical accuracy of the newly introduced relativistic and non-relativistic kinetic energy terms.
In practice, due to the parallelized implementation of both the FDTD and FFT-based algorithms it is possible to use \texttt{quantumfdtd} on quite large lattices.
For example, there have been past studies with lattices as large as $N=1024$~\cite{Strickland:2009ft}.

Possible future improvements include:
\begin{itemize}
\item Implementing an automatized algorithm in order to look for saddle points of $E_{s,i}(\tau)$ for $i=0,\,1,\,2$.
\item Implementing new projection operators on the post-processing stage (see subsection~\ref{subsec_parity_proy}). For instance, separation of $d$ states. 
For larger overlaps this may involve including contributions from higher states in the initial wave-function (see subsection~\ref{subsec_icond}).
\item Implementing a suitable binary format for the wave-functions.
\end{itemize}
Such improvements are beyond the scope of the present work.
However, it would be useful to consider them as a guidance for further improvements of the \texttt{quantumfdtd} program.

\section*{Acknowledgments}
R.L.D. was financially supported by the Ram{\'o}n Areces Foundation, the INFN post-doctoral fellowship AAOODGF-2019-0000329 and the Spanish grant MICINN: PID2019-108655GB-I00.
M.S. was supported by the U.S. Department of Energy, Office of Science, Office of Nuclear Physics Award No.~DE-SC0013470.
J.H.W. was supported by the U.S. Department of Energy, Office of Science, Office of Nuclear Physics and Office of Advanced Scientific Computing Research within the framework of Scientific Discovery through Advance Computing (SciDAC) award Computing the Properties of Matter with Leadership Computing Resources.
We acknowledge the support by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany's Excellence Strategy -- EXC-2094 -- 390783311 via the Excellence Cluster "ORIGINS".
The simulations have been carried out on the computing facilities of the Computational Center for Particle and Astrophysics (C2PAP) and the Leibniz Supercomputing Center (SuperMUC), on the local theory cluster (T30 cluster) of the Physics Department of the Technische Universit{\"a}t M{\"u}nchen (TUM), and on local computing facilities at the INFN-Firenze.
