# Closed-form solutions for the Salpeter equation

Fernando Alonso-Marroquin  CIPR, King Fahd University of Petroleum and
Minerals, Dhahran, 31261, Kingdom of Saudi Arabia
<EMAIL_ADDRESS>Yaoyue Tang  Modelling and Simulation Research
Group, The University of Sydney, Sydney NSW 2006, Australia  Fatemeh Gharari
Department of Statistics and Computer Science, University of Mohaghegh
Ardabili, Ardabil, Iran  M. N. Najafi  Department of Physics, University of
Mohaghegh Ardabili, P.O. Box 179, Ardabil, Iran<EMAIL_ADDRESS>
###### Abstract

We propose integral representations and analytical solutions for the
propagator of the $1+1$ dimensional Salpeter Hamiltonian, describing a
relativistic quantum particle with no spin. We explore the exact Green
function and an exact solution for a given initial condition, and also find
the asymptotic solutions in some limiting cases. The analytical extension of
the Hamiltonian in the complex plane allows us to formulate the equivalent
stochastic problem, namely the Bäumer equation. This equation describes
relativistic stochastic processes with time-changing anomalous diffusion. This
Bäumer propagator corresponds to the Green function of a relativistic
diffusion process that interpolates between Cauchy distributions for small
times and Gaussian diffusion for large times, providing a framework for
stochastic processes where anomalous diffusion is time-dependent.

quantum relativistic mechanics, anomalous diffusion

##  I introduction

The importance of the relativistic quantum mechanics (RQM) is not only due to
the requirements of the standard model of elementary particles [1], but also
as an essential ingredient for describing the new phases of matter, like two-
dimensional electron gas [2], graphene [3], relativistic plasmas [4, 5, 6],
neutron stars [7], Weyl semi-metals [8], relativistic hydrodynamics [9], and
other materials with Dirac cone dispersion in their band structure [10]. RQM
has vast applications in other fields, such as the stochastic phenomena which
is realized using a time Wick’s rotation. The examples are the relativistic
stochastic phenomena [11, 12, 13, 14], the relativistic Ornstein-Uhlenbeck
process [15, 16], the relativistic Brownian motion [17, 18], the relativistic
Fokker-Planck equation [19, 20, 21], the classical and relativistic Gross-
Pitaevskii equation [22] and the relativistic diffusion [23]. Due to its wide
applications, the standard theory for relativistic (classical and quantum)
statistical mechanics is nowadays well-developed [24]. There are different
ways in order to construct a relativistic version of classical quantum
mechanics that respect relativity due to their covariant formulation. The
Klein-Gordon equation [1] containing second time derivative, the Dirac
equations [25] which is of the first order in the time derivative, and the
Salpeter equation [26, 27], which is known as the “square root” of the Klein-
Gordon operator are the famous formulations to this end. The Salpeter equation
governs scalar particles since its relativistic covariance is applied at the
operator level, while the Dirac equation and the Klein-Gordon equation are
capable of describing higher dimensional representations of the Lorenz group.
Although working with a square root of the Klein-Gordon operator makes the
Salpeter equation difficult to handle, the fact that it resolves some
essential problems with the Klein-Gordon equation like the negative energy
particles is a privilege, making it worth putting under intense studies [27].
The square-root operator induces non-locality for Salpeter’s equation, a fact
that makes it troublesome. An analysis for the conserved current, and also
second quantization formulation can be found in [28]. Despite its importance,
there are limited works that focus on its structure and consequences. For
example, the causality structure of this model, its plane wave solutions, the
superposition of quantum states, entanglement structure, second quantized
field theoretical description, and the corresponding particle algebra that
retrieves the causality, the closed form of the green function, and the dual
stochastic theory (after the Wick rotation) is missing in the literature. The
aim of the present paper is to process the later gap, i.e., the closed form of
the propagator of the free particles, and the dual stochastic interpretation.

A motivation for the Salpeter stochastic description is the crossover
phenomena in diffusion processes that takes place for example in the analysis
of stochastic processes, such as the evolution of the probability density
function of the S&P 500 stock market index [29]. This crossover is manifested
in a superdiffusive regime for small times and a Gaussian diffusive regime for
large times. Attempts have been made to use quantum mechanic models to measure
stock return volatility, and reported a power-law relation in the tails of
price return distribution [30]. For small momentum values, the Salpeter
equation is analogous to the Schrödinger equation (where in the dual
stochastic equation, the solution is normal diffusion), while in the large
momentum it corresponds to massless free particles [27]. On the other hand,
the Green function of the stochastic equation is related to the Cauchy
distributions [31]. Therefore, we expect a crossover between these two limits
for intermediate momenta, which induces a crossover in terms of the
spatiotemporal coordinates.

This paper presents analytical solutions of the Salpeter propagator in its
quantum and stochastic versions. The paper is organized as follows: In the
next section, we introduce the Salpeter equation to be analyzed in this paper
and provide the integral representation and analytical solutions of the
propagator. In Sec.III we present the dual stochastic description, namely the
Bäumer equation, and analyze the crossover regimes of the propagator of this
relativistic stochastic process.

##  II Salpeter Hamiltonian

Figure 1: Real and imaginary part of the function $f(z)=\sqrt{1+z^{2}}$, where
$z=p+iq$. The principal branch used in this calculation is obtained from the
branch of the complex square root:
$\sqrt{z}=\sqrt{\rho}(\cos\phi/2)+i*\sin(\phi/2))$, where $z=\rho e^{i\phi}$
and $-\pi\leq\phi\leq\pi$. This branch is extended along the imaginary axis by
imposing $\Im{f(iq)}>0$ in the first quadrant, and $\Im{f(iq)}<0$ in the
fourth quadrant. This guarantee that the defined function is analytical in the
first and fourth quadrant of the complex plane.

Our derivation is based on the classical formulation of the RQM of a free
particle, resulting to an integral representation of the Green (propagator)
function, see Appendix A. The propagator function $G(x,t)$ is used to
calculate the wave function $\psi(x,t)$ in terms of its initial value
$\psi_{0}(x)$ as

$\psi(x,t)=\int_{-\infty}^{\infty}{G(x-x_{0},t)\psi_{0}(x_{0})dx_{0}}.$ (1)

The propagator depends on the Hamiltonian of the particle and is derived in
Appendix A using textbook’s quantum mechanics:

$G(x,t)=\frac{1}{2\pi}\int_{-\infty}^{\infty}{e^{-\frac{i}{\hbar}(H(p)t+px)}dp}$
(2)

Let us assume that the Hamiltonian is an even function of the momentum
$H(-p)=H(p)$, There are two consequences from this property: First, the
propagator will be an event function with respect to $x$, i.e.
$G(-x,t)=G(x,t)$. Second, the propagator can be calculated as the sum
$G=G_{+}+G_{-}$ as

$G_{\pm}(x,t)=\frac{1}{2\pi}\int_{0}^{\infty}{e^{-\frac{i}{\hbar}H(p)t\pm\frac{i}{\hbar}px}dp}$
(3)

Adding $G_{+}(x,t)$ and $G_{-}(x,t)$ we obtain an useful expression of the
propagator $G(x,t)$ as

$\begin{split}G(x,t)&=G_{+}(x,t)+G_{-}(x,t)\\\
&=\frac{1}{2\pi}\int_{0}^{\infty}{e^{-\frac{i}{\hbar}H(p)t}(e^{\frac{i}{\hbar}px}+e^{-\frac{i}{\hbar}px})}dp\\\
&=\frac{1}{\pi}\int_{0}^{\infty}{e^{-\frac{i}{\hbar}H(p)t}\cos(\frac{px}{\hbar})}dp\end{split}$
(4)

Let us consider the Hamiltonian of a free, massive, and spin-less relativistic
particle in one dimension, namely the Salpeter Hamiltonian

$H_{S}(p)=\sqrt{m^{2}c^{4}+c^{2}p^{2}}-mc^{2}$ (5)

Due to the square root, there is no standard differential equation for this
Hamiltonian. In Eq. (5) we subtract the rest energy $mc^{2}$ to recover the
classical Schrödinger equation in the limit of small velocities $|p|<<mc$.
This limit leads to

$H_{S}(p)\approx\frac{p^{2}}{2m}$ (6)

Then the classical Schrödinger equation can be expressed in space coordinates
as a partial differential equation

$i\partial_{t}\psi(x,t)=-\frac{\hbar^{2}}{2m}\partial_{x}^{2}\psi(x,t),$ (7)

A simple transformation of the Schrödinger equation in quantum mechanics can
lead to the diffusion equation in a stochastic equation. This transformation
can be done by either transforming time into imaginary $(t\to it)$ or
transforming the Hamiltonian into imaginary $(H\to iH)$. Either transformation
leads to the normal diffusion in the stochastic process

$\partial_{t}P(x,t)=D\partial_{x}^{2}P(x,t),$ (8)

where $D=\hbar^{2}/(2m)$. For the relativistic case, it is not possible to
express the Schödinger equation as a standard partial differential equation
and the solution of the propagators is not well-known. Dirac was able to
provide a relativistic equation by converting the wave equation into a four-
dimensional function that inherits the spin of the particles. In this paper,
we present an alternative by assuming that the wave function is a scalar so
that it represents a relativistic quantum particle without spin. As shown in
Appendix A, the quantum bra-ket (bracket) formalism allows us to express the
propagator as integral given by Eq. (2) without the need to rely on fractional
differential equations. The analysis relies instead on performing the
analytical continuation of the Salpeter Hamiltonian on the complex plane.
First, the Hamiltonian is written as

$H_{S}(p)=mc^{2}\left(f\left(\frac{p}{mc}\right)-1\right),$ (9)

where $f(p)=\sqrt{1+p^{2}}$. Then we define $f(z)$, where $z=p+iq$, as the
analytical continuation of $f(p)$. Figure 1 shows the real and imaginary part
of the principal analytical branch of the complex function
$f(z)=\sqrt{1+z^{2}}$. The real part of the function has a saddle point at
$z=0$. The function is analytic in some regions containing the real axis
$\imaginary(z)=0$. This is where the real Hamiltonian function is defined. The
complex function is analytical in the complex plane except for the two branch
cuts $[-i\infty,-i)$ and $[i,i\infty)$ where the derivatives do not exist. We
can slightly modify this function by using the function
$f^{*}(p+iq)=\lim_{\epsilon\to 0}{f(p+q+i\epsilon)}$, where $\epsilon>0$. This
will correct the values along the imaginary axis so that the new function
becomes analytical in the first and fourth quadrants, including the imaginary
axis.

For the sake of simplicity, for the rest of the paper we rescale the variables
$x$, $t$, and $p$ in such a way that $\hbar=c=1$. Thus Eq. 5 becomes

$H_{S}(p)=\sqrt{m^{2}+p^{2}}-m$ (10)

Then substituting Eq. 10 into Eq 4, the propagator we want to solve can be
defined as

$G(x,t)=\frac{1}{2\pi}\int_{-\infty}^{\infty}{e^{-iH_{S}(p)t+ipx}dp}$ (11)

We will also use the splitting the integral as $G=G_{+}+G_{-}$ where

$G_{\pm}(x,t)=\frac{1}{2\pi}\int_{0}^{\infty}{e^{-iH_{S}(p)t\pm ipx}dp}$ (12)

The numerical solution to the above integrals is very unstable due to the
heavy oscillatory behaviour of the integrands. We will show that more stable
numerical solutions are possible by using complex variable analysis to obtain
integral representation of the propagators. We will also provide analytical
solutions for these integral representations. In the next sections we present
stable integral representations and analytical solutions in different cases.
We will consider the analytical extension for the integrands and use complex
variable analysis to handle the integrals.

###  II.1 Classical Limit: case $m>0$ and $t\gg m$

We start with the integration of Eq. 11 by considering the case when $t$ is
large. First, we perform an analytical continuation fo the integrand of Eq. 11
in the complex plane. We notice that the integrand of this equation is an
analytical function in a region in the complex plane containing the real axis
($z=p+iq$ and $q=0$). We consider $t$ real and positive. When t is large, the
major contribution of this integrand is when $H_{S}(p)$ is imaginary and
minimal, which corresponds to the value $p=0$. Under these conditions, we can
appeal to Debbie’s method based on the saddle point approximation presented in
Section 31 of the Book of Dennery and Krzywicki [32]. The method consists of
replacing the Hamiltonian by its expansion up to the quadratic term around is
minimal point $p=0$, i.e. $H_{S}(p)\approx p^{2}/(2m)$. Thus,

$G(x,t)\sim\frac{1}{2\pi}\int_{-\infty}^{\infty}{e^{-i\frac{p^{2}}{2m}t+ipx}dp}=\sqrt{\frac{m}{2\pi
it}}e^{-\frac{mx^{2}}{2it}}$ (13)

Here we use the following integration formula

$\int_{-\infty}^{\infty}{e^{iap-
bp^{2}}dp}=\sqrt{\frac{\pi}{b}}e^{\frac{a^{2}}{4b}}$ (14)

Note that this solution corresponds to the propagator of classical quantum
mechanics. The evolution of a wave function can be obtained from its initial
condition using Eq. 74. We can also recover the Green function of the
diffusion equation by making the transformation $t\to it$ and $1/(2m)\to D$:

$G(x,t)\to\frac{1}{\sqrt{\pi Dt}}e^{-\frac{x^{2}}{Dt}}$ (15)

The propagator accounts for the probability distribution function of the
particle that was at $x=0$ at $t=0$. which corresponds to the classical
diffusion process.

###  II.2 Integral representation, Case $m>0$ and $0<x<t$

Here we derive an integral representation of the propagator whose integrand
does not pose heavy oscillatory behaviour. The method is based on the Cauchy
Theorem of complex variable analysis, that is used to convert the integral on
the real axis into an integral over imaginary axis that removes heavy
oscillations. The transformed integral are thus well-defined and easily
integrable. The requirements to transform an integral over the real axis into
the imaginary axis relies on two conditions: The complex continuation of the
integrand is analytical in a quadrant of the complex plane; and the integrand
vanishes as the complex argument is large enough in that quadrant. The details
of the use of complex variable are included in two Theorems formulated and
demonstrated in Appendix C. Here we will just show that more stable numerical
solutions are possible by using these Theorems to obtain integral
representation of the propagators. We will also provide analytical solutions
for these integral representations.

First we provide an integral representation for the propagator of relativistic
particles with arbitrary mass $m>0$ inside the light cone ($0<x<t$). We
replace Eq. 10 into Eq. 12 to obtain

$G_{\pm}(x,t)=\frac{1}{2\pi}\int_{0}^{\infty}{e^{\Phi_{\pm}(p)}dp},$ (16)

where

$\Phi_{\pm}(p)=-i\sqrt{m^{2}+p^{2}}t-mt\pm ipx.$ (17)

Now we consider that for $0<x<t$ and $-\pi/2<\theta<=0$ the condition

$\lim_{R\to\infty}{\Re{\Phi(Re^{i\theta})}}=\lim_{R\to\infty}{-R(-t\pm
x)\sin(\theta)}=-\infty.$ (18)

This condition is satisfied in the fourth quadrant in the complex plane. Thus,
applying Theorem 2 (See Appendix C) to the analytical branch shown in Fig. 1
and substituting $p$ by $-iq$ we have

$\begin{split}&G_{\pm}(x,t)=\frac{e^{imt}}{2\pi}\int_{0}^{\infty}{e^{-i\sqrt{m^{2}+p^{2}}t\pm
ipx}dp}\\\
&=\frac{e^{imt}}{2\pi}\int_{0}^{\infty}{e^{-i\sqrt{m^{2}+(-iq)^{2}}t\pm
i(-iq)x}d(-iq)}\\\
&=\frac{e^{imt}}{2i\pi}\left(\int_{0}^{m}+\int_{m}^{\infty}\right){e^{-i\sqrt{m^{2}+(-iq)^{2}}t\pm
qx}dq}\\\
&=\frac{e^{imt}}{2i\pi}\left(\int_{0}^{m}{e^{-i\sqrt{m^{2}-q^{2}}t\pm
qx}dq}+\int_{m}^{\infty}{e^{-\sqrt{q^{2}-m^{2}}t\pm qx}dq}\right).\end{split}$
(19)

The full propagator $G(x,t)=G_{+}(x,t)+G_{-}(x,t)$ becomes

$G(x,t)=\frac{e^{imt}}{i\pi}\int_{0}^{\infty}{e^{-\sqrt{q^{2}-m^{2}}t}\cosh{(qx)}}dq$
(20)

The integrand in this equation does not have singular oscillations and it
easily integrated by standard numerical methods.

###  II.3 Integral representation, Case $m>0$ and $x>t>0$

Now we obtain the expression for the propagator outside of hte light cone
($x>t>0$). Under these conditions $\Phi_{-}$ in Eq. 17 satisfies the condition
Eq. 18 in the four quadrant where $\sin{\theta}<0$. For $\Phi_{+}$ we note
that $-t+x>0$. For this case the integral is calculated in the first quadrant
where $\sin{\theta}>0$. Therefore we apply only the loop integral for the
first quadrant for $\Phi_{+}$ and for the fourth quadrant for $\Phi_{-}$ to
Eq. 12. Applying Theorem 2 on $G_{-}$ for the first quadrant and using the
change of variable $p\to iq$

$\begin{split}&G_{-}(x,t)=\frac{e^{imt}}{2\pi}\int_{0}^{\infty}{e^{-i\sqrt{m^{2}+p^{2}}t-ipx}dp}\\\
&=\frac{e^{imt}}{2\pi}\int_{0}^{\infty}{e^{-i\sqrt{m^{2}+(-iq)^{2}}t-i(-iq)x}d(-iq)}\\\
&=\frac{e^{imt}}{2i\pi}\left(\int_{0}^{m}{e^{-i\sqrt{m^{2}-q^{2}}t-qx}dq}+\int_{m}^{\infty}{e^{-\sqrt{q^{2}-m^{2}}t-qx}dq}\right).\end{split}$
(21)

And applying Theorem 1 on $G_{+}$ for the fourth quadrant and using the change
of variable $p\to-iq$

$\begin{split}&G_{+}(x,t)=\frac{e^{imt}}{2\pi}\int_{0}^{\infty}{e^{-i\sqrt{m^{2}+p^{2}}t+ipx}dp}\\\
&=\frac{e^{imt}}{2\pi}\int_{0}^{\infty}{e^{-i\sqrt{m^{2}+(iq)^{2}}t+i(iq)x}d(iq)}\\\
&=\frac{-e^{imt}}{2i\pi}\left(\int_{0}^{m}{e^{-i\sqrt{m^{2}-q^{2}}t-qx}dq}+\int_{m}^{\infty}{e^{\sqrt{q^{2}-m^{2}}t-qx}dq}\right).\end{split}$
(22)

Adding both parts and using the identity $1/i=-i$, the full propagator becomes

$G(x,t)=\frac{ie^{imt}}{\pi}\int_{m}^{\infty}{\sinh{(\sqrt{q^{2}-m^{2}}t)}e^{-qx}dq}$
(23)

Note that this solution corresponds to non-zero solution for the propagator
outside of this light cone, suggesting faster-than-light propagation. While
the relativistic quantum mechanic formulation allows for the possibility of
particles having speeds higher than light, the green function decays rapidly
outside the light cone. This is associated with certain challenges and
inconsistencies, such as causality violations and the lack of a clear physical
interpretation. Quantum field theories, by treating particles as excitations
of a vacuum state, –i. e. quantum fields– resolve this problem, where the
Green function is expressed explicitly in terms of the expectation value of
quantum fields. More explicitly, one considers the commutation
(anticommutation) relations for bosonic and fermionic creation and destruction
operators. As a result the retarded and advanced Green functions are found to
be zero (non-zero) outside (inside) the light cone, the expression of which is
different for fermions and bosons [33].

###  II.4 Massless particles: Case $m=0$

When the particle has no mass, the propagator for $0<x<t$ is calculated by
taking $m=0$ in the last line of Eq. 19

$G_{\pm}(x,t)=\frac{1}{2i\pi}\int_{0}^{\infty}{e^{-qt\pm
qx}dq}=\frac{1}{2i\pi}\frac{-1}{-t\pm x}.$ (24)

Adding both parts, the full propagator of mass-less particles for $x<t$
becomes

$G(x,t)=\frac{1}{\pi}\frac{it}{x^{2}-t^{2}}$ (25)

For the case $0<t<x$ the propagator is obtained taking $m=0$ in Eq. 21 and 22

$\begin{split}&G_{-}(x,t)=\frac{1}{2i\pi}\int_{0}^{\infty}{e^{-qt-
qx}dq}=\frac{1}{2i\pi}\frac{1}{x+t}\\\
&G_{+}(x,t)=\frac{i}{2\pi}\int_{0}^{\infty}{e^{qt-
qx}dq}=\frac{i}{2\pi}\frac{1}{x-t}\\\ \end{split}$ (26)

Adding both parts, the full propagator of mass-less particles for $x>t$
becomes

$G(x,t)=\frac{1}{\pi}\frac{it}{x^{2}-t^{2}}$ (27)

This result is the same that Eq. 25 for case $x<t$. We notice a singularity at
the light cone $x=t$ and, again, non-zero values at $x>t>0$ indicating faster-
than-light-propagation. As in the case of massive particles, a quantum field
theory formulation is required to avoid the apparent violation of the
causality principle in the classical interpretation of quantum mechanics.

Figure 2: (a) Comparison between the numerical and analytical solution of the
Salpeter propagator $G(x,t)$ concerning $x$, across various $t$ values. (a)
shows the real part and (b) the positive (blue) and negative (red) values of
the imaginary parts in log-log scale. The findings from both methods exhibit
strong consistency and alignment. Notably, the propagator displays a
singularity at the light cone, represented by $x=t$. (b) The phase change
compared between the numerical and analytical solutions, depicted in terms of
$\frac{1}{2}-\frac{\theta(x,t)}{\pi}$ with respect to $x$. For all $t$ values,
the phase consistently converges to $-\pi/2$ in both the numerical and
analytical solutions. (c) and (d) are the Wave function calculated using Eq.
37. The real part is plotted in blue, and the imaginary part is plotted in
green. (c) For a massless particle $m=0$. (d) For a particle $m=1$.

###  II.5 Analytical Solution

Within this subsection, we establish the analytical solution for the Salpeter
propagator. Substitute the Salpeter Hamiltonian (Eq. 10) into the propagator
(Eq. 4) we have

$\begin{split}G(x,t)&=\frac{1}{\pi}\int_{0}^{\infty}{e^{-i(\sqrt{m^{2}+p^{2}}-m)t}\text{cos}(px)}dp\\\
&=\frac{1}{\pi}e^{imt}\int_{0}^{\infty}{e^{-i\sqrt{m^{2}+p^{2}}t}\text{cos}(px)}dp\\\
\end{split}$ (28)

Using the identity [27]

$\begin{split}\int_{0}^{\infty}&{dx\,\,\text{exp}(-\alpha\sqrt{x^{2}+\beta^{2}})\text{cos}\gamma
x}\\\
&=\frac{\alpha\beta}{\sqrt{\alpha^{2}+\gamma^{2}}}K_{1}(\beta\sqrt{\alpha^{2}+\gamma^{2}}),\text{Re}[\alpha]>0,\text{Re}[\beta]>0\end{split}$
(29)

We can rewrite the propagator as

$\begin{split}G(x,t)&=\frac{1}{\pi}e^{imt}\int_{0}^{\infty}{e^{-\alpha\sqrt{\beta^{2}+p^{2}}}\text{cos}(\gamma
p)}dp\\\ \end{split}$ (30)

where $\alpha=it$, $\beta=m$, and $\gamma=x$. Then according to the right-hand
side of Eq. 29, we have

$\begin{split}G(x,t)&=\frac{imte^{imt}}{\pi\sqrt{x^{2}-t^{2}}}K_{1}(m\sqrt{x^{2}-t^{2}}),\
x^{2}>t^{2}\end{split}$ (31)

where $K_{1}$ is the modified Bessel function of the second kind. For
$x^{2}<t^{2}$, we can utilize the property of the modified Bessel functions,
given as:

$K_{\alpha}(x)=\frac{\pi}{2}i^{\alpha+1}H_{\alpha}^{(1)}(ix),\,\,-\pi<\text{arg}\,x\leq\frac{\pi}{2}$
(32)

where $H_{\alpha}^{(1)}$ is the Hankel function of the first kind. Replacing
$x^{2}-t^{2}$ by $-(t^{2}-x^{2})$ in Eq. 31 we have:

$\begin{split}G(x,t)&=\frac{imte^{imt}}{\pi\sqrt{-(t^{2}-x^{2})}}\frac{\pi}{2}i^{2}H_{1}^{(1)}(im\sqrt{-(t^{2}-x^{2})})\\\
&=-\frac{mte^{imt}}{2\sqrt{t^{2}-x^{2}}}H_{1}^{(1)}(-m\sqrt{t^{2}-x^{2}}),\
x^{2}<t^{2}\end{split}$ (33)

In summary, the analytical solution of the Salpeter propagator is

$G(x,t)=\begin{cases}\frac{imte^{imt}}{\pi\sqrt{x^{2}-t^{2}}}K_{1}\left(m\sqrt{x^{2}-t^{2}}\right),&x^{2}>t^{2}\\\
-\frac{mte^{imt}}{2\sqrt{t^{2}-x^{2}}}H_{1}^{(1)}\left(-m\sqrt{t^{2}-x^{2}}\right),&x^{2}<t^{2},\end{cases}$
(34)

where $H_{1}^{(1)}$ is the Hankel function of first kind, and $K_{1}$ is the
modified Bessel function of the second kind.

In Fig. 2, we present the numerical and analytical solutions in terms of the
real and imaginary part of the propagator and the phase at varying $t$. Both
solutions are indistinguishable. We also compare the cases $m=1$ and $m=0$.
The numerical solution inside the light cone $x<t$ is obtained by numerical
integration of Eq 20. The numerical solution outside the cone $x>t$ is the
numerical integration of Eq. 23. We note that the hyperbolic functions of
these integrals are sensitive to large arguments, but this issue can be solved
by expressing them in terms of exponential functions. The analytical solution
is calculated by substituting the corresponding values of $x$ and $t$ into
Equation 34. In this context, we choose $t=2^{n}/100,n=0,1,2...8$, and
$x\in[0,5t]$ as the calculation range. The solutions for $m=1$ and $m=0$ are
asymptotically similar when $d=|x-t|\to 0$. Here $d$ is the distance to the
light cone $|x|=t$. For each time $t$, the propagators present a singularity
at $d=0$. This singularity proves to be of the order $1/d$ so that it produces
a sudden change in the values of the propagator and a discontinuity of the
argument of this function. As we will see in the next section, this
singularity should be carefully treated once we calculate the wave function of
the Salpeter equation.

###  II.6 Salpeter wave function

Using the Eq. 1, one finds the wave function in later times, given its form in
the initial time. After the change of variable $x^{\prime}=x-x_{0}$ this
equation can be written as

$\psi(t,x)=\int_{-\infty}^{\infty}{G(x^{\prime},t)\psi_{0}(x-x^{\prime})dx^{\prime}}.$
(35)

Strictly speaking, this integral does not exist because the propagator has two
first order singularities at $x=\pm t$. To obtain a meaningful solution for
the wave function, we need to circumvent the singularities by choosing an
appropriate integration path in the complex plane. A meaningful solution is
obtained if the integration is performed on the deformed path $C$ shown in Fig
3. We redefine the wave function as the path integral

$\psi^{*}(t,x)=\int_{C}{G(z,t)\psi_{0}(x-z)dz}.$ (36)

This integral can be performed by applying the theorem of residues of the
complex variable [32] over the closed path shown in Fig 3. After calculating
the residues inside the closed path, the integration becomes:

$\begin{split}\psi^{*}(t,x)&=PV\int_{-\infty}^{\infty}{G(x^{\prime},t)\psi_{0}(x-x^{\prime})dx^{\prime}}\\\
&+\frac{1}{2}\psi_{0}(x-t)+\frac{1}{2}\psi_{0}(x+t).\end{split}$ (37)

where $PV$ refers to the principal value of the integral [32]. To see the
proof of Eq. 37, and other details, refer to Appendix D. For the initial
condition $\psi_{0}(x)$ we choose a wave function with compact support in the
interval $[-\delta/2,\delta/2]$. The simplest case is to assume that the wave
function is initially constant and real on the interval, by normalization we
found that the initial wave function is constant in this interval and zero
elsewhere $\psi_{0}(x)=\delta^{-1/2}(|x|<\delta/2)$. Using this initial
condition turns out to be problematic since the wave function is non-zero at
$x\pm\delta$, and it leads to four singularities in the principal value in Eq.
37 at $|x-t|\pm\delta$. To remove these singularities, we need to assume that
the initial wave function is zero at $x\pm\delta$. In other words, we need to
guarantee that the initial wave function is continuous. Our next best option
for the initial condition is an initially harmonic wave function:

$\psi_{0}(x)=\sqrt{\frac{2}{\delta}}\cos{\frac{\pi x}{\delta}}(|x|<\delta/2)$
(38)

For $m>0$, We calculate the wave function based on the numerical integration
of the propagator, and the results are plotted in Figure 2c. In spite of the
singularity of the propagator at $x=\pm t$, the wave functions are shown to be
continuous and differentiable everywhere. The propagation of the wave function
is somewhat similar to the one in the wave equation, where one initial pulse
is split into two equal pulses moving left and right. The main difference here
is that the pulses are complex and include not only propagation but also
dissipation. The probability density is calculated as $\rho=\phi^{*}\phi$ and
shown in Figure 2d. The continuity equation is verified by calculating the
limit of the area under $\rho$ for different values of the value $\epsilon$
used in the calculation of the principal value of the wave functions. It is
verified that for all values of time used in Figure 2, the limit of the
integral is very close to one with an error below $0.004\%$. This is
consistent with the continuity equation for the Salpeter equation, already
demonstrated by Kowalski and Rembielinski [27]

Figure 3: Contour in the complex plane used to calculate the Salpeter wave
function.

##  III Bäumer Hamiltonian

Now we consider the Relativistic Diffusion Equation (RDE), which is obtained
by a Wick rotation, or equivalently, multiplying the Hamiltonian by the
imaginary unit $i$. The corresponding Hamiltonian has been investigated by
Bäumer et al. [31] and will be referred here as the Bäumer Hamiltonian

$H_{B}(p)=i\sqrt{m^{2}+p^{2}}-im.$ (39)

Replacing this Hamiltonian into the Schrödinger equation 57, and using Eq. 69
the RDE is obtained as the following fractional differential equation:

$2D[\sqrt{1+\frac{\partial^{2}}{\partial x^{2}}}-1]P=\frac{\partial
P}{\partial t},$ (40)

where $D\equiv\frac{m}{2}$ and $P$ is the probability density in the
stochastic space, equivalent to the wave function in the quantum space. The
quantum mechanics’ mathematical formalism allows a simple solution of this
fractional differential equation by solving the propagator that in this case
is the Green function of the RDE.

The propagator in Eq. 63 for the Hamiltonian given by Eq. 39 becomes

$G(x,t)=\frac{1}{2\pi}\int_{-\infty}^{\infty}{e^{-\sqrt{m^{2}+p^{2}}t+mt+ipx}dp}$
(41)

whose components are

$G_{\pm}(x,t)=\frac{1}{2\pi}\int_{0}^{\infty}{e^{-\sqrt{m^{2}+p^{2}}t+mt\pm
ipx}dp}$ (42)

We are considering the different cases to discuss the solution below:

###  III.1 Classical Limit: Case $m>0$ and $t\gg m$

The classical limit is the same as the asymptotic solution when $t\to\infty$.
In such case the propagator 41 becomes

$G(x,t)=\frac{1}{2\pi}\int_{-\infty}^{\infty}{e^{-(p^{2}/2m)t+ipx}dp}$ (43)

The solution corresponds to a Gaussian distribution

$G(x,t)=\sqrt{\frac{m}{{2\pi t}}}e^{-\frac{mx^{2}}{2t}}$ (44)

We can also recover the Green function of the diffusion equation by making the
transformation $t\to it$ in the classical Schrödinger equation and taking
$1/(2m)\to D$:

$G(x,t)\to\frac{1}{\sqrt{\pi Dt}}e^{-\frac{x^{2}}{Dt}}$ (45)

###  III.2 Integral representation for $m>0$ and $0\leq x<\infty$ (Inner
solution)

For the general case, the components in Eq.42 can be summed up to obtain;

$G(x,t)=\frac{e^{mt}}{\pi}\int_{0}^{\infty}{e^{-\sqrt{m^{2}+p^{2}}t}\cos{(px)}dp}$
(46)

We note that the integrand of has strongly oscillatory behavior as
$x\to\infty$. This integral is therefore the inner solution of the propagator
that can be applied for moderate values of $x$ but will not effectively
converge at the tails. An outer solution is sought in the next section using
complex variable analysis.

###  III.3 Integral representation for $m>0$ and $x>0$ (Outer solution)

We write the two terms of the propagator as Eq. 16,

$G_{\pm}(x,t)=\frac{1}{2\pi}\int_{0}^{\infty}{e^{\Phi_{\pm}(p)}dp},$ (47)

and the asymptotic behaviour of $\Phi_{\pm}(z)=-iH_{B}(z)t\pm xz$ for
$z=Re^{i}\theta$ as $R\to\infty$ follows

$\lim_{R\to\infty}{\Re{\Phi(Re^{i\theta})}}=\lim_{R\to\infty}{-R(t\cos{\theta}\pm
x\sin{\theta}})$ (48)

The condition given by Eq. 48 is satisfied in the first quadrant of the
complex plane for $G_{+}$ and in the fourth quadrant for $G_{-}$. The
integration of $G_{+}$ over the first quadrant gives

$\begin{split}&G_{+}(x,t)=\frac{1}{2\pi}\int_{0}^{\infty}{e^{-\sqrt{m^{2}+p^{2}}t+mt+ipx}dp}\\\
&=\frac{e^{mt}}{2\pi}\int_{0}^{\infty}{e^{-\sqrt{m^{2}+(iq)^{2}}t+i(iq)x}d(iq)}\\\
&=\frac{ie^{mt}}{2\pi}\left(\int_{0}^{m}{e^{-\sqrt{m^{2}-q^{2}}t-qx}dq}+\int_{m}^{\infty}{e^{-i\sqrt{q^{2}-m^{2}}t-qx}dq}\right)\\\
\end{split}$ (49)

The integration of $G_{-}$ over the four quadrants results in

$\begin{split}&G_{-}(x,t)=\frac{1}{2\pi}\int_{0}^{\infty}{e^{-\sqrt{m^{2}+p^{2}}t+mt-
ipx}dp}\\\
&=\frac{e^{mt}}{2\pi}\int_{0}^{\infty}{e^{-\sqrt{m^{2}+(-iq)^{2}}t-i(-iq)x}d(-iq)}\\\
&=\frac{-ie^{mt}}{2\pi}(\int_{0}^{m}{e^{-\sqrt{m^{2}-q^{2}}-qx}dq}+\int_{m}^{\infty}{e^{i\sqrt{q^{2}-m^{2}}-qx}dq})\end{split}$
(50)

Adding both parts and the full propagator becomes

$G(x,t)=\frac{-e^{mt}}{\pi}\int_{m}^{\infty}{\sin{(\sqrt{q^{2}-m^{2}}t)}e^{-qx}dq}.$
(51)

Interestingly, Using the identity $sinh(iz)=isin(z)$ and replacing $t$ by $it$
we can recover the outer-cone solution of the Salpeter equation, Eq. 23.

###  III.4 Case $m=0$

For massless stochastic particles, the components of the propagator in Eq. 42
becomes

$G_{\pm}(x,t)=\frac{1}{2\pi}\int_{0}^{\infty}{e^{-pt\pm ipx}dp}$ (52)

Both integrals can be analytically integrated. Summing up, the solution is

$G_{m=0}(x,t)=\frac{1}{\pi}\frac{t}{t^{2}+x^{2}}$ (53)

This is the Cauchy distribution with scale parameter given by $t$.

###  III.5 Case $x\gg t$

The behavior of the tail of the Bäumer distribution can be obtained by using
$x>>t$ in Eq.51.

$G(x,t)=\frac{-e^{m(t-x)}}{\pi}$ (54)

Interestingly, the Bäumer distribution exhibits exponential-tailed behavior
that is not the case for either the Cauchy distribution or the Gaussian
distribution.

Figure 4: Solution of the Bäumer propagation that describes the relativistic
diffusion propagator. $G(x,t)$ is plotted in blueagainst $x$ for different
values of $t$. The propagator is compared to the Cauchy distribution (green)
and Gaussian distribution (red). The inset shows the exponential law of the
tails. Figure 5: Evolution of the second moment of the Bäumer propagator
$<x^{2}>$ against $t$. The solid line corresponds to power law $y=t$

###  III.6 Analytical solution

The analytical solution of the Bäumer propagator is obtained by substituting
the imaginary Hamiltonian into the Salpeter propagator

$\begin{split}G(x,t)&=\frac{1}{\pi}\int_{0}^{\infty}{e^{-i^{2}(\sqrt{m^{2}+p^{2}}-m)t}\text{cos}(px)}dp\\\
&=\frac{1}{\pi}e^{-mt}\int_{0}^{\infty}{e^{\sqrt{m^{2}+p^{2}}t}\text{cos}(px)}dp\\\
\end{split}$ (55)

Using the identity in Eq. 29, we have $\alpha=-t,\beta=t$ and $\gamma=x$, and
the solution to the integral can be written as:

$G(x,t)=-\frac{mte^{-mt}}{\pi\sqrt{x^{2}+t^{2}}}K_{1}(m\sqrt{x^{2}+t^{2}})$
(56)

Note that the solution of the propagator defines a diffusion process with a
transition from Cauchy to diffusion, shown in Fig. 4. The propagator (i.e.,
The Green function) converges asymptotically to a Cauchy distribution as $t\to
0$, a Gaussian distribution as $t\to\infty$, and an exponential function when
$x\to\infty$.

The Bäumer propagator defines a diffusion process that is more complex than
Gaussian diffusion. In many cases, the diffusivity is characterized by a power
law on the evolution of the peak $G(0,t)\sim t^{-\alpha}$, or the power law of
the evolution of the second moment $\langle x^{2}G(x,t)\rangle_{x}\sim
t^{-\alpha}$. The diffusion process corresponds to normal diffusion when
$\alpha=2$, super-diffusion when $\alpha>2$, and sub-diffusion when
$\alpha<2$. Fig. 5 shows the time evolution of both the peak and second moment
of the propagator. The evolution of the peak shows a crossover between
subdiffusion ($\alpha=1$) for small times, and normal diffusion ($\alpha=2$)
for large times. A different diffusion process results from the second moment
that exhibits normal diffusion at all times. This complexity of the diffusion
process with crossover regions over space-time is also observed in the
temporal evolution of the PDFs of the predicted price return [29]. However,
the values of the exponents are different from the Bäumuer diffusion process.
In particular, the diffusion process in stock markets is not only fractional
but also nonlinear. Thus, the use of quantum mechanics in stock market
prediction is still severely limited due to the current limitations on the
mathematical description of the nonlinearity of quantum mechanics. On the
other hand, the link between stochastic equations and quantum mechanics by the
Wick rotation provides interesting possibilities, like the option of solving
complex non-linear stochastic problems, such as the Porous Media Fractional
Equation [34, 35] and then using the Wick rotation to obtain the corresponding
non-linear quantum problems.

##  IV Conclusions

We presented the numerical, analytical, and asymptotic solutions of the
Salpeter propagator for relativistic quantum mechanics of a free particle with
no spin, and its stochastic counterpart, the Bäumer propagator. The solution
of the Bäumer propagator converges asymptotically to Cauchy distribution for
small times, and to Gaussian distribution with large times. The Salpeter
propagation is recovered by transforming time by imaginary time. The solution
of this quantum mechanical propagator shows important properties, such as non-
zero values outside of the light cone, propagation with diffusion around the
light cone, and singularity of the light cone. Special treatment is required
to handle this singularity in the calculation of the wave function using the
concept of principal value of integrals. A compact-supported initial wave
function with no discontinuities at the edges allows us to remove the
singularity, leading to well-posed wave functions. The conservation of
probability is checked by testing that the integral of the probability density
remains one for all times.

This paper has focused on the use of complex variable analysis for the
derivation of solutions and the careful testing of these solutions using
numerical analysis. The interpretation of these solutions in the context of
relativistic quantum mechanics, quantum field theories, and relativistic
stochastic dynamics is a matter of future publication.

## References

  * Greiner _et al._ [2000] W. Greiner _et al._ , _Relativistic quantum mechanics_ , Springer, 2000, vol. 2.
  * Melrose and Hayes [1984] D. Melrose and L. M. Hayes, _Aust. J. Phys_ , 1984, 37, 639–50.
  * Novoselov _et al._ [2005] K. S. Novoselov, A. K. Geim, S. V. Morozov, D. Jiang, M. I. Katsnelson, I. V. Grigorieva, S. Dubonos and a. Firsov, _nature_ , 2005, 438, 197–200.
  * Umstadter [2001] D. Umstadter, _Physics of Plasmas_ , 2001, 8, 1774–1785.
  * Umstadter _et al._ [1996] D. Umstadter, S.-Y. Chen, A. Maksimchuk, G. Mourou and R. Wagner, _Science_ , 1996, 273, 472–475.
  * Arnold _et al._ [2002] P. Arnold, G. D. Moore and L. G. Yaffe, _Journal of High Energy Physics_ , 2002, 2002, 030.
  * Watts _et al._ [2016] A. L. Watts, N. Andersson, D. Chakrabarty, M. Feroci, K. Hebeler, G. Israel, F. K. Lamb, M. C. Miller, S. Morsink, F. Özel _et al._ , _Reviews of Modern Physics_ , 2016, 88, 021001.
  * Lv _et al._ [2015] B. Lv, H. Weng, B. Fu, X. P. Wang, H. Miao, J. Ma, P. Richard, X. Huang, L. Zhao, G. Chen _et al._ , _Physical Review X_ , 2015, 5, 031013\.
  * Rezzolla and Zanotti [2013] L. Rezzolla and O. Zanotti, _Relativistic hydrodynamics_ , Oxford University Press, 2013.
  * Malko _et al._ [2012] D. Malko, C. Neiss and A. Görling, _Physical Review B_ , 2012, 86, 045443.
  * Dunkel and Hänggi [2009] J. Dunkel and P. Hänggi, _Physics Reports_ , 2009, 471, 1–73.
  * Debbasch and Chevalier [2007] F. Debbasch and C. Chevalier, AIP Conference Proceedings, 2007, pp. 42–48.
  * De Angelis _et al._ [1986] G. De Angelis, G. Jona-Lasinio, M. Serva and N. Zanghi, _Journal of Physics A: Mathematical and General_ , 1986, 19, 865.
  * Deng [2009] W. Deng, _SIAM journal on numerical analysis_ , 2009, 47, 204–226.
  * Debbasch _et al._ [1997] F. Debbasch, K. Mallick and J.-P. Rivet, _Journal of statistical physics_ , 1997, 88, 945–966.
  * Rigotti [2005] M. Rigotti, _Ph.D. thesis_ , Master’s thesis, ETH Zürich, Institute for Theoretical Physics, 2005.
  * Pal and Deffner [2020] P. Pal and S. Deffner, _New Journal of Physics_ , 2020, 22, 073054\.
  * Kakushadze [2017] Z. Kakushadze, _Physica A: Statistical Mechanics and its Applications_ , 2017, 475, 59–76.
  * Hakim [1968] R. Hakim, _Journal of Mathematical Physics_ , 1968, 9, 1805–1818.
  * Kleinert and Zatloukal [2013] H. Kleinert and V. Zatloukal, _Physical Review E_ , 2013, 88, 052106\.
  * Serva [1988] M. Serva, Annales de l’IHP Physique théorique, 1988, pp. 415–432.
  * Kleinert [2012] H. Kleinert, _Europhysics Letters_ , 2012, 100, 10001.
  * Chevalier and Debbasch [2008] C. Chevalier and F. Debbasch, _Journal of mathematical physics_ , 2008, 49, 043303.
  * Mi Hakim [2011] R. Mi Hakim, _Introduction to relativistic statistical mechanics: classical and quantum_ , World scientific, 2011.
  * Thaller [2013] B. Thaller, _The dirac equation_ , Springer Science & Business Media, 2013\.
  * Foldy [1956] L. L. Foldy, _Physical Review_ , 1956, 102, 568.
  * Kowalski and Rembieliński [2011] K. Kowalski and J. Rembieliński, _Physical Review A_ , 2011, 84, 012108.
  * Lämmerzahl [1993] C. Lämmerzahl, _Journal of mathematical physics_ , 1993, 34, 3918–3932.
  * Alonso-Marroquin _et al._ [2019] F. Alonso-Marroquin, K. Arias-Calluari, M. Harré, M. N. Najafi and H. J. Herrmann, _Physical Review E_ , 2019, 99, 062313.
  * Ahn _et al._ [2024] K. Ahn, L. Cong, H. Jang and D. S. Kim, _Financial Innovation_ , 2024, 10, 6.
  * Baeumer _et al._ [2010] B. Baeumer, M. M. Meerschaert and M. Naber, _Physical Review E_ , 2010, 82, 011132.
  * Dennery and Krzywicki [1996] P. Dennery and A. Krzywicki, _Mathematics for physicists_ , Dover Publications, Inc., 1996.
  * Peskin [2018] M. E. Peskin, _An introduction to quantum field theory_ , CRC press, 2018\.
  * Tang _et al._ [2024] Y. Tang, F. Gharari, K. Arias-Calluari, F. Alonso-Marroquin and M. N. Najafi, _Phys. Rev. E_ , 2024, 109, 024310.
  * Gharari _et al._ [2021] F. Gharari, K. Arias-Calluari, F. Alonso-Marroquin and M. N. Najafi, _Physical Review E_ , 2021, 104, 054140.
  * Hildebrand [1962] F. B. Hildebrand, _(No Title)_ , 1962.

##  Appendix A Quantum Propagator

We can write the Schrödinger equation for a particle in the braket notation:

$\hat{H}\ket{\psi}=i\hbar\partial_{t}\ket{\psi};\quad\ket{\psi(t_{0})}=\ket{\psi_{0}}.$
(57)

Here the $\hat{H}$ is the Hamiltonian of the particle, $\ket{\psi}$ the state
of the particle at time $t$, and $\ket{\psi_{0}}$ the initial condition at
time $t_{0}$. Let us define the temporal evolution operator by

$\hat{U}(t_{0},t)\ket{\psi(t_{0})}=\ket{\psi(t)}.$ (58)

Using Eq. 57, the evolution equation of $\hat{U}$ be written as

$i\hbar\partial_{t}\hat{U}=\hat{H}\hat{U};\quad\hat{U}(t_{0})=\hat{1}.$ (59)

Since the particle is in a state $\ket{\psi(t_{0})}$ at $t=t_{0}$, the
operation of temporal evolution brings the particle to a new state
$\ket{\psi(t)}$. Then, the probability of finding the particle in a position
$x$ is calculated by a braket operation $\psi(t,x)=\braket{x}{\psi(t)}$. This
is the wave function of the particle, which can be calculated by multiplying
Eq. 58 by $\bra{x}$

$\psi(t,x)=\bra{x}\hat{U}(t_{0},t)\ket{\psi(t_{0})}.$ (60)

Let us introduce the completeness operator for the observable $\hat{x}$

$\hat{1}=\int{dx_{0}\ket{x_{0}}\bra{x_{0}}}$ (61)

Inserting this operator into Eq. 60 we obtain

$\begin{split}\psi(t,x)&=\bra{x}\hat{U}(t_{0},t)\hat{1}\ket{\psi(t_{0})}\\\
&=\int{dx_{0}\bra{x}\hat{U}(t_{0},t)\ket{x_{0}}\braket{x_{0}}{\psi(t_{0})}}\\\
&=\int{dx_{0}G(x,t,x_{0},t_{0})\psi(t_{0},x_{0})}.\end{split}$ (62)

Here we define the propagator by the function

$G(x,t,x_{0},t_{0})=\bra{x}\hat{U}(t_{0},t)\ket{x_{0}}$ (63)

According to Eq. 74, it is noted that the propagation behaves similarly to the
Green functions of the partial differential equation. This function can be
calculated for some specific cases. If the Hamiltonian is independent of time,
the evolution operator of Eq. 59 can be easily solved to obtain

$\hat{U}(t,t_{0})=e^{-\frac{i}{\hbar}\hat{H}(t-t_{0})}.$ (64)

Eq. 64 can be replaced with Eq. 63 to obtain the explicit expression of the
propagator. This requires to introduce the monochromatic wave functions that
result from the eigenstates of the momentum observable $\hat{p}$

$\ket{p}=\frac{1}{\sqrt{2\pi}}\int{dx_{0}e^{\frac{i}{\hbar}px_{0}\ket{x_{0}}}}$
(65)

The wave functions of the eigenstates $\ket{x}$ of the observables $\hat{x}$
is given by the Dirac’s delta function

$\psi_{x_{0}}(x)=\braket{x}{x_{0}}=\delta(x-x_{0})$ (66)

In agreement with Eq. 61, The property of the Dirac function hold for any
function $f(x)$

$\int{dxf(x)\delta(x-x_{0})}=f(x_{0}).$ (67)

The corresponding wave function of the eigenstates $\ket{p_{0}}$ of the
observable $\hat{p}$ is given by

$\psi_{p_{0}}(x)=\braket{x}{p_{0}}=\frac{1}{\sqrt{2\pi}}\int{dx_{0}e^{\frac{i}{\hbar}px_{0}}\braket{x}{x_{0}}}=\frac{1}{\sqrt{2\pi}}e^{\frac{i}{\hbar}p_{0}x}$
(68)

Let us note that $i\hbar\partial_{x}\psi_{p_{0}}(x)=p\psi_{p_{0}}(x)$ which
indicates that

$\bra{x}\hat{p}\ket{p_{0}}=-i\hbar\partial_{x}\braket{x}{k}$ (69)

The completeness relation for the $\hat{p}$ observable also holds

$\hat{1}=\int{dp_{0}\ket{p_{0}}\bra{p_{0}}}$ (70)

Let us assume a Hamiltonian of the free particle

$\hat{H}=\hat{H}(p)$ (71)

Using Eq. 70 in the propagator (Eq. 63) gives

$\begin{split}&G(x,t,x_{0},t_{0})=\bra{x}e^{-\frac{i}{\hbar}\hat{H}(t-t_{0})}\ket{x_{0}}\\\
&=\int{dp\bra{x}e^{-\frac{i}{\hbar}\hat{H}{(p)}(t-t_{0})}\ket{p}\braket{p}{x_{0}}}\\\
&=\int{dpe^{-\frac{i}{\hbar}H{(p)}(t-t_{0})}\braket{x}{p}\braket{p}{x_{0}}}\\\
&=\frac{1}{2\pi}\int{dp\exp{-\frac{i}{\hbar}H{(p)}(t-t_{0})+\frac{i}{\hbar}p(x-x_{0})}}\end{split}$
(72)

Hence the propagator satisfies

$G(x,t,x_{0},t_{0})=G(x-x_{0},t-t_{0}).$ (73)

Using this property, and assuming that $t_{0}=0$, we can simplify the
calculation of the wave function by Eq. 62 to

$\psi(t,x)=\int_{-\infty}^{\infty}{G(x-x_{0},t)\psi_{0}(x_{0})dx_{0}}.$ (74)

Where $\psi_{0}(x)$ is the initial condition and

$G(x,t)=\frac{1}{2\pi}\int_{-\infty}^{\infty}{e^{-\frac{i}{\hbar}H(p)t+i\hbar
px}dp}$ (75)

##  Appendix B Scaling Relations

In this appendix, we focus on the scaling properties of the Salpeter equation,
which helps us to recognize the general features of the model. Note that the
full functionality of the green function is $G(x,t;x^{\prime},t^{\prime};m)$,
but without loss of generality we set $x^{\prime}=t^{\prime}\equiv 0$, and
$G(x,t,m)\equiv G(x,t;0,0;m)$, keeping in mind that eventually we should
replace $x\to x-x^{\prime}$ and $t\to t-t^{\prime}$. We have the following
integral expression

$G(x,t,m)\equiv\frac{1}{2\pi}e^{imt}\int_{0}^{\infty}e^{-i\sqrt{m^{2}+p^{2}}t}\cos{(px)}\text{d}p.$
(76)

One can easily show by inspection that for any positive scaling parameter
$\lambda>0$ we have

$G(\lambda x,\lambda t,m)=\lambda^{-1}G(x,t,\lambda m),$ (77)

so that under the transformation

$x\rightarrow\lambda x\ ,\ t\rightarrow\lambda t\ ,\
m\rightarrow\lambda^{-1}m$ (78)

we have $G\rightarrow\lambda^{-1}G$. Therefore, we expect the following
identity:

$G(x,t)=m\tilde{G}\left[m(ax+bt)\right],$ (79)

where $\tilde{G}$ is a general function, that should be determined using other
methods, and $a$ and $b$ are some constant complex numbers. The argument of
$\tilde{G}$ guarantees the invariance of $\tilde{G}$ under scaling
transformation. The limit $m\to 0$ helps to identify this function. To this
end, we add a vanishingly small quantity to time, i.e. $t\to t-i0^{+}$ (which
guarantees that the integral is well-defined in the UV limit $p\to\infty$), so
that

$G(x,t,m=0)=\frac{1}{2\pi}\left[\frac{i}{x-t+i0^{+}}-\frac{i}{x+t-i0^{+}}\right].$
(80)

Using the identity

$\frac{1}{X+i0^{\pm}}=\text{P.V.}\frac{1}{X}\mp i\pi\delta(X),$ (81)

where P.V. stands for principal value, we obtain

$G(x,t,m=0)=\frac{it}{\pi(x^{2}-t^{2})}+\frac{1}{2}\delta(x-t)+\frac{1}{2}\delta(x+t),$
(82)

where $\delta(x)$ is a Dirac delta function, and the first term on the right
is P.V. One may compare this result with the orthogonality identity for plane
waves

$\frac{1}{2\pi}\int_{-\infty}^{\infty}e^{ip(x\pm
t)}\text{d}p=\frac{1}{\pi}\lim_{p\to\infty}\frac{\sin p(x\pm t)}{x\pm
t}=\delta(x\pm t),$ (83)

from which we realize that the first term in the r.h.s. of Eq. 82 is due to
the presence of absolute value in the dispersion relation $E=\left|p\right|$
in the $m\to 0$ limit. From this limit we realize that $a=1$ and $b=\pm 1$. We
additionally find that for very small $m$ we have

$\tilde{G}(x,t,\text{small}\
m)=\frac{1}{2\pi}\left[\frac{i}{m(x-t)+i0^{+}}-\frac{i}{m(x+t)-i0^{+}}\right].$
(84)

The other quantity of interest is $\left\langle x^{2}\right\rangle$. From the
symmetry Eq. 77 one finds that

$\left\langle x^{2}\right\rangle=m^{-2}\tilde{f}(mt)=t^{2}\tilde{g}(mt),$ (85)

where $\tilde{f}(y)=y^{2}\tilde{g}(y)$ are general smooth functions. For the
massless free particle limit, we know

$\left\langle x^{2}\right\rangle_{\text{massless}}=t^{2},$ (86)

from which we see

$\left.\tilde{g}(y)\right|_{m\to 0}\to 1.$ (87)

##  Appendix C Integral Transformation Theorem

Analytical integration of improper functions is available only in simple
cases. If the integrand poses strong oscillatory behavior, numerical
integration is also inconvenient and often lacks accuracy. Using the Cauchy
theorem of complex variable analysis, it is possible to transform integrals
with oscillatory behavior to integrands with smooth exponential decreases that
are easier to integrate.

######  Theorem 1.

If $f(z)=e^{\Phi(z)}$ is is a function that is analytical in the first quarter
of the complex plane:

$D_{1}=\\{z=\rho e^{i\theta}=p+iq;0\leq\theta\leq\pi/2\\}$ (88)

and $\Phi(z)$, with $z=\rho e^{i\theta}=p+iq$ satisfies

$\lim_{\rho\to\infty}{\Re{\Phi(\rho e^{i\theta})}}=-\infty\quad z\in D_{1}$
(89)

its integral over the semi-real axis $\mathbb{R}\cap D_{n}$ can be transformed
into an integral over the semi-imaginary axis $\mathbb{I}\cap D_{1}$ as

$\int_{0}^{\infty}{e^{\Phi(p)}dp}=\int_{0}^{\infty}{e^{\Phi(iq)}idq}.$ (90)

###### Proof.

Since $f(z)$ is analytical in the first quadrant, its integral over any of the
circular sectors on the first quadrant is zero

$\oint_{D_{1}}{f(z)dz}=\oint_{D_{1}}{e^{\Phi(z)}dz}=0$ (91)

Applying this property to a circular sector of infinite radius,

$\int_{0}^{\infty}{e^{\Phi(p)}dp}+\lim_{\rho\to\infty}\int_{0}^{\frac{\pi}{2}}{e^{\Phi(\rho
e^{i\theta})}d\theta}+\int_{\infty}^{0}{e^{\Phi(iq)}idq}=0$ (92)

Thus,

$\int_{0}^{\infty}{e^{\Phi(p)}dp}=\int_{0}^{\infty}{e^{\Phi(iq)}idq}-\lim_{\rho\to\infty}\int_{0}^{\frac{\pi}{2}}{e^{\Phi(\rho
e^{i\theta})}d\theta}$ (93)

The second term is expected to vanish under the condition Eq. 90 ∎

Now we provide a similar theorem for the four quadrants of the complex plane.

######  Theorem 2.

If $f(z)=e^{\Phi(z)}$ is is a function that is analytical in the four quarter
of the complex plane:

$D_{4}=\\{z=\rho e^{i\theta}=p+iq;-\pi/2\leq\theta\leq 0\\}$ (94)

and $\Phi(z)$, with $z=\rho e^{i\theta}=p+iq$ satisfies

$\lim_{\rho\to\infty}{\Re{\Phi(\rho e^{i\theta})}}=-\infty\quad z\in D_{4}$
(95)

the integral over the semi-real axis $\mathbb{R}\cap D_{4}$ can be transformed
into an integral over the semi-imaginary axis $\mathbb{I}\cap D_{4}$ as

$\int_{0}^{\infty}{e^{\Phi(p)}dp}=\int_{0}^{\infty}{e^{\Phi(-iq)}(-i)dq}$ (96)

###### Proof.

The theorem will be proof for the fourth quadrant. Since $f(z)$ is analytical
in the four quadrants, its integral over any of the circular sectors on the
first quadrant is zero

$\oint_{D_{4}}{f(z)dz}=\oint_{D_{4}}{e^{\Phi(z)}dz}=0$ (97)

Applying this property to a circular sector of infinite radius,

$\int_{0}^{\infty}{e^{\Phi(p)}dp}+\lim_{\rho\to\infty}\int_{0}^{-\frac{\pi}{2}}{e^{\Phi(\rho
e^{i\theta})}d\theta}+\int_{-\infty}^{0}{e^{\Phi(iq)}idq}=0$ (98)

Thus, replacing $q$ by $-q$ in the last integral,

$\int_{0}^{\infty}{e^{\Phi(p)}dp}=\int_{0}^{\infty}{e^{\Phi(-iq)}(-idq)}-\lim_{\rho\to\infty}\int_{0}^{-\frac{\pi}{2}}{e^{\Phi(\rho
e^{i\theta})}d\theta}.$ (99)

The second term is expected to vanish under the condition 95 ∎

##  Appendix D Theorem of Residues

The main focus is to calculate the following integral

$\psi^{*}(t,x)=\int_{C}{G(z,t)\psi_{0}(x-z)dz}.$ (100)

where $C$ is the contour in the complex plane shown in Fig. 3. The function
$G(z,t)$ is the analytical continuation of the propagator (Eq. 11) in the
complex plane:

$G(z,t)=\frac{1}{2\pi}\int_{-\infty}^{\infty}{e^{-iH_{S}(p)t+ipz}dp}$ (101)

and the function $\psi_{0}(z)$ is the initial wave function (Eq. 38)

$\psi_{0}(z)=\sqrt{\frac{2}{\delta}}\cos{\frac{\pi z}{\delta}}(|z|<\delta/2).$
(102)

In this appendix, we provide a proof of Eq. 37. The integral 37 diverges
because of the non-integrable poles on the integration contour. Since these
poles are of the first order, the Cauchy principal value $(PV)$ of the
integral exists and corresponds to avoiding the poles with semi-circular arcs.
To prove, we use the Cauchy residue theorem which states that if the integrand
is analytic on a simply connected domain except for a finite number of
isolated singularities, then the integral of that around a closed curve is
equal to $2\pi i$ times the sum of the residues of the integrand at the
singular points inside the closed curve. First, we consider the case $m=0$ and
show that the $PV$ is:

$PV\int_{-\infty}^{\infty}{G(y,t)\psi_{0}(x-y)dy}=\frac{1}{2}\left(\psi_{0}(x-t)+\psi_{0}(x+t)\right).$
(103)

The contour in Fig. 3 is chosen so that there are two poles inside it and so
that the little circles around each of the poles are so small. The contour is

$C=C^{R}+C^{\epsilon}_{1}+C^{\epsilon}_{2}+C^{\infty}_{-\infty},$ (104)

where $C^{\infty}_{-\infty}$ is the integral in Eq. 103,
$C^{\epsilon}_{i},\,i=1,2,$ are the contributions due to integrating over the
poles at $\pm t$, and $C^{R}$ is the (vanishing) contribution due to the large
arc.
To integrate around the pole $"-t"$, let us parameterize $z=-t+\epsilon
e^{i\theta}$, then $dz=i\epsilon e^{i\theta}d\theta$. For the case $m=0,$ by
considering Eq. 27, integrating around the associated semi-circle yields

$\displaystyle\int_{C^{1}_{\epsilon}}$ $\displaystyle{G(z,t)\psi_{0}(x-z)dz}$
$\displaystyle=\lim_{\epsilon\to 0}\int_{\pi}^{2\pi}{G(-t+\epsilon
e^{i\theta},t)\psi_{0}(x+t-\epsilon e^{i\theta})i\epsilon e^{i\theta}d\theta}$
$\displaystyle=\lim_{\epsilon\to
0}\frac{-t}{\pi}\int_{\pi}^{2\pi}{\dfrac{\psi_{0}(x+t-\epsilon
e^{i\theta})}{\epsilon e^{i\theta}(-2t+\epsilon e^{i\theta})}\epsilon
e^{i\theta}d\theta}$ $\displaystyle=\frac{1}{2}\psi_{0}(x+t).$ (105)

Similarly, to integrate around the pole $t$, we parameterize $z=t+\epsilon
e^{i\theta}$, then $dz=i\epsilon e^{i\theta}d\theta$. By considering Eq. 27,
integrating around the other semi-circle yields

$\int_{C^{2}_{\epsilon}}{G(z,t)\psi_{0}(x-z)dz}=\frac{1}{2}\psi_{0}(x-t).$
(106)

Also, the integration along the semi-circular portion $C^{R}$ of the contour
$C$ vanishes as the radius of integration becomes very large, i.e., as
$R\to\infty,$

$\int_{C^{R}}{G(z,t)\psi_{0}(x-z)dz}=0,$ (107)

and by considering the counter 104, this completes the proof of Eq. 37.
To prove Eq. 103, since there are two poles inside of the contour, we have by
Cauchy’s residue theorem,

$\int_{C}{G(z,t)\psi_{0}(x-z)dz}=2i\pi\left(Res(-t)+Res(t)\right).$ (108)

By using Theorem IV(p. 557) from [36] and Equations 105, and 106,

$Res(\pm t)=\frac{1}{2i\pi}\psi_{0}(x\mp t),$ (109)

and by considering the counter 104, this completes the proof of Eq. 103.
Now, we consider the case $m\neq 0.$ We follow a similar scenario on the same
contour. Therefore, under these conditions, Equations 104, 107, and 108 work
for our case and we should investigate Eq. 103. To achieve this, we apply the
analytical solution of the Salpeter propagator, as described in Eq. 34, for
integrals 105 and 106. We emphasize the case $m\neq 0$ by denoting

$G^{m}(z,t)=\begin{cases}\frac{imte^{imt}}{\pi\sqrt{z^{2}-t^{2}}}K_{1}\left(m\sqrt{z^{2}-t^{2}}\right),&|z|>t\\\
\\\
-\frac{mte^{imt}}{2\sqrt{t^{2}-z^{2}}}H_{1}^{(1)}\left(-m\sqrt{t^{2}-z^{2}}\right),&|z|<t,\end{cases}$
(110)

where $H_{1}^{(1)}$ is the Hankel function of first kind, and $K_{1}$ is the
modified Bessel function of the second kind.
Integrating around the first semi-circle yields

$\displaystyle\int_{C^{1}_{\epsilon}}$
$\displaystyle{G^{m}(z,t)\psi_{0}(x-z)dz}$ $\displaystyle=\lim_{\epsilon\to
0}\int_{\pi}^{2\pi}{G^{m}(-t+\epsilon e^{i\theta},t)\psi_{0}(x+t-\epsilon
e^{i\theta})i\epsilon e^{i\theta}d\theta}$ $\displaystyle=\lim_{\epsilon\to
0}\frac{-imte^{imt}}{\pi}\int_{\pi}^{2\pi}\frac{K_{1}\left(m\epsilon^{1/2}e^{i\theta/2}\sqrt{-2t+\epsilon
e^{i\theta}}\right)}{\epsilon^{1/2}e^{i\theta/2}\sqrt{-2t+\epsilon
e^{i\theta}}}...$
$\displaystyle\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,...\times\psi_{0}(x+t-\epsilon
e^{i\theta})i\epsilon e^{i\theta}d\theta$

Applying the fact $K_{1}(y)\to\frac{1}{y}$ as $y\to o,$ we get

$\int_{C^{\epsilon}_{1}}{G^{m}(z,t)\psi_{0}(x-z)dz}=\frac{1}{2}\psi_{0}(x+t)e^{imt}.$
(111)

Similarly, integrating around the second semi-circle yields

$\int_{C^{\epsilon}_{2}}{G^{m}(z,t)\psi_{0}(x-z)dz}=\frac{1}{2}\psi_{0}(x-t)e^{imt}.$
(112)

We can easily acquire

$Res(\pm t)=\frac{e^{imt}}{2i\pi}\psi_{0}(x\mp t),$ (113)

then, the $PV$ is:

$\displaystyle PV\int_{-\infty}^{\infty}$
$\displaystyle{G(y,t)\psi_{0}(x-y)dy}$
$\displaystyle=\frac{1}{2}\left(\psi_{0}(x-t)+\psi_{0}(x+t)\right)e^{imt}.$
(114)

##  Appendix E Analytical expressions based on perturbative expansion

In this section we present an analytical expansion of the salpeter integral,
which serves as a perturbative expansion. It helps much to realize the
analytical form of the function around e.g. the light cone, and other regions
of interest, and also shows the internal structure and consistency of the
analytical expressions. We follow the expansion scheme presented according to
an ad hoc weight function. In this method we multiply the integrand by a
weight function $f(p,p_{0})$ which makes the approximate integration feasible
in an appropriate limit, and $p_{0}$ is a characteristic cutoff momentum
dividing the integral domain to two intervals, each being treated differently
in the perturbation expansion. In our case, $p_{0}\equiv m$ is the
characteristic momentum that separates “large” and “short” $p$ values. The
choice of $f(p,p_{0})$ is not unique, and can be any form depending on the
form of the integral, for example it can be either ($p_{0}\equiv m$)

$\begin{split}f(p,p_{0})&=\exp\left[-\left(\frac{p}{p_{0}}\right)^{2}\right]\\\
&\text{or}\\\ f(p,p_{0})&=\Theta(p-p_{0}),\end{split}$ (115)

where $\Theta(x)$ is a step function. In this paper choose the second “cut”
function, and $p_{0}\equiv m$. We then cast the formula Eq. 2 to the
following:

$\begin{split}G(x,t,0,0)=&\int_{0}^{\infty}\text{d}p\left[f(p,p_{0})+1-f(p,p_{0})\right]\times\\\
&\exp[-itH(p)]\cos\left(px\right)\\\ &\equiv
G_{1}(x,t)+G_{2}(x,t).\end{split}$ (116)

The integrand of function

$G_{1}(x,t)=\int_{0}^{\infty}\text{d}pf(p,p_{0})\exp[-itH(p)]\cos\left(px\right)$
(117)

can be expanded for small $p$ values ($p\ll p_{0}$), while

$G_{2}(x,t)=\int_{0}^{\infty}\text{d}p\left[1-f(p,p_{0})\right]\exp[-itH(p)]\cos\left(px\right)$
(118)

should be expanded around small $\frac{1}{p}\ll\frac{1}{p_{0}}$. We end up to
($\epsilon_{1}\equiv\frac{p}{m}$, and
$g_{1}(mt;\epsilon_{1})\equiv\exp\left[-imt\left(\sqrt{1+\epsilon_{1}^{2}}\right]-1\right)$,
$\rho\equiv mx$)

$\begin{split}G_{1}(x,t)&=m\int_{0}^{1}\text{d}\epsilon_{1}e^{-imt}g_{1}(mt;\epsilon_{1})\cos(\epsilon_{1}mx)\\\
&=m\sum_{n=0}^{\infty}e^{-imt}\frac{g_{1}^{(n)}\left(mt\right)}{n!}\int_{0}^{1}\text{d}\epsilon_{1}\epsilon_{1}^{n}\cos(\epsilon_{1}mx)\\\
&=me^{-imt}\sum_{n=0}^{\infty}\frac{g_{1}^{(n)}\left(mt\right)}{n!}\text{Re}\left(-i\frac{\partial}{\partial\rho}\right)^{n}\int_{0}^{1}\text{d}\epsilon_{1}e^{i\rho\epsilon_{1}}\\\
&=me^{-imt}\sum_{n=0}^{\infty}\frac{g_{1}^{(n)}\left(mt\right)}{n!}f_{n}(mx)\end{split}$
(119)

where

$g_{1}^{(n)}(mt)\equiv\frac{\partial^{n}}{\partial\epsilon_{1}^{n}}\left.g_{1}(mt;\epsilon_{1})\right|_{\epsilon_{1}=0},$
(120)

and

$\begin{split}f_{n}(\rho)&\equiv\text{Re}\left[\left(-i\frac{\partial}{\partial\rho}\right)^{n}\left(\frac{e^{i\rho}-1}{i\rho}\right)\right]\\\
&=\left\\{\begin{matrix}(-1)^{\frac{n+1}{2}}\left(\frac{\partial}{\partial\rho}\right)^{n}\left(\frac{\cos\rho-1}{\rho}\right)&n=\text{odd}\\\
(-1)^{\frac{n}{2}}\left(\frac{\partial}{\partial\rho}\right)^{n}\left(\frac{\sin\rho}{\rho}\right)&n=\text{even}\end{matrix}\right..\end{split}$
(121)

Note that $g_{1}^{(n)}=0$ for all odd $n$ values, while for the even ones we
have

$\begin{split}&g_{1}^{(0)}(y)=1,\\\ &g_{1}^{(2)}(y)=-iy,\\\
&g_{1}^{(4)}(y)=3y(i-y),\\\ &g_{1}^{(6)}(y)=15y(-3i+3y+iy^{2}),\\\
&g_{1}^{(8)}(y)=15y(15i+15y-6iy^{2}+y^{3}),\\\ &...\end{split}$ (122)

The asymptotic behavior of $g_{1}^{(2n)}$ is

$g_{1}^{(2n)}(y)\to\left\\{\begin{matrix}A_{n}^{(1)}y^{n}&y\gg 1\\\
B_{n}^{(1)}y&y\ll 1\end{matrix}\right.$ (123)

where $A_{n}^{(1)}$ and $B_{n}^{(1)}$ are two expansion coefficients that are
given in Eq. 122.
$f_{n}(\rho)$ can also be written in terms of generalized hypergeometric
function as follows

$f_{n}(\rho)=\frac{{{}_{1}F_{2}}\left[\frac{n+1}{2},\left\\{\frac{1}{2},\frac{n+3}{2}\right\\},-\frac{\rho^{2}}{4}\right]}{n+1}.$
(124)

The explicit form of $f_{n}$ is

$\begin{split}&f_{0}(\rho)=\frac{\sin\rho}{\rho}\\\
&f_{1}(\rho)=\frac{\cos\rho-1+\rho\sin\rho}{\rho^{2}}\\\
&f_{2}(\rho)=\frac{2\rho\cos\rho+(\rho^{2}-2)\sin\rho}{\rho^{3}}\\\
&f_{3}(\rho)=\frac{6+3(\rho^{2}-2)\cos\rho+\rho(\rho^{2}-6)\sin\rho}{\rho^{4}}\\\
&f_{4}(\rho)=\frac{4\rho(\rho^{2}-6)\cos\rho+(24-12\rho^{2}+\rho^{4})\sin\rho}{\rho^{5}}\\\
&...\end{split}$ (125)

Note that in the limit $\rho\to 0$ we have

$\begin{split}f_{n}(\rho)&=\text{Re}\left[\left(-i\frac{\partial}{\partial\rho}\right)^{n}\sum_{m=1}^{\infty}\frac{1}{m!}(i\rho)^{m-1}\right]\\\
&=\text{Re}\left[\frac{n!}{(n+1)!}+i\frac{(n+1)!}{(n+2)!}\rho+\frac{(n+2)...4\times
3}{(n+3)!}\rho^{2}+O(\rho^{3})\right]\\\
&=\frac{1}{n+1}-\frac{\rho^{2}}{2(n+3)}+O(\rho^{3}).\end{split}$ (126)

For the large $\rho$ limit one can see by inspection that

$\begin{split}\left.f_{n}(\rho)\right|_{\rho\gg 1}\to
a_{n}\frac{\sin\rho}{\rho}\end{split}$ (127)

where $a_{n}=1$ for all $n$ values.
For $G_{2}$ we have:

$\begin{split}G_{2}(x,t)&=\int_{m}^{\infty}\text{d}pe^{-ipt\sqrt{1+\left(\frac{m}{p}\right)^{2}}}\cos
px\\\ &=\int_{m}^{\infty}\text{d}pg_{2}(mt,\epsilon_{2})e^{-ipt}\cos px\\\
&=m\sum_{n=0}^{\infty}\frac{g_{2}^{(n)}(mt)}{n!}\int_{1}^{\infty}\text{d}\tilde{p}\tilde{p}^{-n}e^{-i\tilde{p}mt}\cos(\tilde{p}mx)\end{split}$
(128)

where $\tilde{p}\equiv p/m$,
$\epsilon_{2}\equiv\frac{1}{\epsilon_{1}}=\frac{m}{p}$, and

$\begin{split}g_{2}(mt,\epsilon_{2})&\equiv\exp\left[-i\frac{mt}{\epsilon_{2}}\left(\sqrt{1+\epsilon_{2}^{2}}-1\right)\right]\\\
&=\sum_{n=0}^{\infty}\frac{1}{n!}g_{2}^{(n)}(mt)\epsilon_{2}^{n}=\sum_{n=0}^{\infty}\frac{1}{n!}g_{2}^{(n)}(mt)\tilde{p}^{-n},\end{split}$
(129)

and

$g_{2}^{(n)}(mt)\equiv\left.\frac{\partial^{n}}{\partial\epsilon_{2}^{n}}g_{2}(mt,\epsilon_{2})\right|_{\epsilon_{2}=0}$
(130)

is the coefficients of the Taylor expansion. Notice that when
$\epsilon_{2}\rightarrow 0$,
$\epsilon_{2}^{-1}\left(\sqrt{1+\epsilon_{2}^{2}}-1\right)\rightarrow 0$.
The coefficients are easily obtained by inspection as follows:

$\begin{split}&g_{2}^{(0)}(y)=(-i)^{0},\\\
&g_{2}^{(1)}(y)=(-i)^{1}\frac{y}{2},\\\
&g_{2}^{(2)}(y)=(-i)^{2}\frac{y^{2}}{2^{2}},\\\
&g_{2}^{(3)}(y)=(-i)^{3}\frac{y}{2^{3}}(3!+y^{2}),\\\
&g_{2}^{(4)}(y)=(-i)^{4}\frac{y^{2}}{2^{4}}(4!+y^{2}),\\\
&g_{2}^{(5)}(y)=(-i)^{5}\frac{y}{2^{5}}(2\times
5!+\frac{1}{2}5!y^{2}+y^{4}),\\\
&g_{2}^{(6)}(y)=(-i)^{6}\frac{y^{2}}{2^{6}}(\frac{5}{2}\times
6!+5!y^{2}+y^{4}),\\\ &g_{2}^{(7)}(y)=(-i)^{7}\frac{y}{2^{7}}(5\times
7!+\frac{3}{2}7!y^{2}+\frac{7}{4}5!y^{4}+y^{6}),\\\
&g_{2}^{(8)}(y)=(-i)^{8}\frac{y^{2}}{2^{8}}(7\times
8!+\frac{14}{3}7!y^{2}+\frac{7}{15}6!y^{4}+y^{6})...\end{split}$ (131)

Note that

$g_{2}^{(n)}\to\left\\{\begin{matrix}(-i)^{n}\frac{y^{n}}{2^{n}}&y\gg 1\\\
(-i)^{n}\frac{n!}{2^{n}}(A_{n}^{(2)})y^{m_{n}}&y\ll 1\end{matrix}\right.$
(132)

where $m_{n}=2$ and $1$ for $n$ even and odd respectively, and $A_{n}^{(2)}$
is an expansion coefficient given in Eq. 131. To process the Eq. 128 we use
the fact that

$\begin{split}I_{n}(mt,mx)&\equiv\int_{1}^{\infty}\text{d}\tilde{p}\tilde{p}^{-n}e^{-i\tilde{p}mt}\cos(\tilde{p}mx)\\\
&=\frac{1}{2}\text{Ei}_{n}\left[im(t-x)\right]+\frac{1}{2}\text{Ei}_{n}\left[im(t+x)\right],\end{split}$
(133)

where $\text{Ei}_{n}(y)\equiv\int_{1}^{\infty}z^{-n}e^{-yz}\text{d}z$ is the
exponential integral function. Therefore, we find that

$\begin{split}G_{2}(x,t)=\frac{m}{2}\sum_{n=0}^{\infty}\frac{g_{2}^{(n)}(mt)}{n!}&\left(\text{Ei}_{n}\left[im(t-x)\right]\right.\\\
&\left.+\text{Ei}_{n}\left[im(t+x)\right]\right)\end{split}$ (134)

Note that the singularities of the propagator are due to the presence of the
exponential integral functions in $G_{2}(x,t)$. Given that

$\begin{split}\text{Ei}_{0}(z)&=\frac{e^{-z}}{z},\\\ \lim_{z\to
0}\text{Ei}_{1}(z)&=-\gamma-\ln z,\\\ \lim_{z\to
0}\text{Ei}_{n}(z)&=\frac{1}{n-1}\ \ ,\ n>1\end{split}$ (135)

where $\gamma$ is an Euler–Mascheroni constant, we find that in the limit
$x\to t$, $G_{2}$ becomes singular according to (the same happens for $x=-t$)

$G_{2}(x,t)\to\frac{1}{2i(t-x)}+\frac{im^{2}t}{4}\ln[im(t-x)].$ (136)

Finally, we find

$G(x,t;0,0)=m\sum_{n=0}^{\infty}\frac{g_{1}^{(n)}\left(mt\right)}{n!(n+1)}{{}_{1}F_{2}}\left[\frac{n+1}{2},\left\\{\frac{1}{2},\frac{n+3}{2}\right\\},-\frac{m^{2}x^{2}}{4}\right]+\frac{m}{2}\sum_{n=0}^{\infty}\frac{g_{2}^{(n)}(mt)}{n!}\left(\text{Ei}_{n}\left[im(t-x)\right]+\text{Ei}_{n}\left[im(t+x)\right]\right).$
(137)

This expansion is singular at the Dirac cone, i.e. $x=t$ and $x=-t$ due to the
presence of the exponential integral functions. More precisely, from Eq. 136
we have:

$\begin{split}\left.G(x,t;0,0)\right|_{x\to\pm t}&=\frac{1}{2i(t\mp
x)}+\frac{im^{2}t}{4}\ln[im(t\mp x)]\\\
&=\frac{1}{2z}+\frac{m^{2}z}{8}\ln(mz),\end{split}$ (138)

where we have defined $z\equiv\tau+ix$, $\bar{z}\equiv\tau-ix$, and
$\tau\equiv it$ is the Wick-rotated (imaginary) time. In the above equation we
used $z=\tau+ix\approx 2it$ in the limit $x\to t$ (the same holds for $x\to-t$
with $z$ replaced by $\bar{z}$). All in all we find that

$\begin{split}&\left.G(x,t;0,0)\right|_{x\to
t}+\left.G(x,t;0,0)\right|_{x\to-t}=\\\
&G(x,t,m=0)+\frac{m^{2}}{8}\left[z\ln(mz)+\bar{z}\ln(m\bar{z})\right].\end{split}$
(139)

In Appendix B we discuss the scaling properties of the propagator. Especially
we show that the integral is proportional to the modified Bessel functions of
the second type. More precisely,

$\begin{split}G(x,t)&=\frac{imte^{imt}}{\pi\sqrt{x^{2}-t^{2}}}K_{1}\left(m\sqrt{x^{2}-t^{2}}\right),\
x^{2}\geq t^{2}\\\
&=-\frac{mte^{imt}}{2\sqrt{t^{2}-x^{2}}}H_{1}^{(1)}\left(-m\sqrt{t^{2}-x^{2}}\right),\
x^{2}<t^{2},\end{split}$ (140)

where $K_{1}(y)$ is the modified Bessel function of second type, and
$H_{1}^{(1)}$ is a Hankel function of first type. This can also be written in
terms of the imaginary time $\tau$ as follows

$\begin{split}G(x,\tau)&=\frac{m}{\pi}\cos\theta
e^{mr\cos\theta}K_{1}\left(mr\right),\ x^{2}\geq t^{2}\\\
&=\frac{im}{2}\cos\theta
e^{mr^{\prime}\cos\theta}H_{1}^{(1)}\left(mr^{\prime}\right),\
x^{2}<t^{2},\end{split}$ (141)

where $r\equiv\sqrt{x^{2}+\tau^{2}}$, $r^{\prime}\equiv\sqrt{t^{2}-x^{2}}=ir$,
and $\theta\equiv\tan^{-1}\left(\frac{x}{\tau}\right)$. For the stochastic
process (SP) applications, one may use the first branch, as follows

$G_{\text{SP}}(r,\theta)=\frac{m}{\pi}\cos\theta
e^{mr\cos\theta}K_{1}\left(mr\right).$ (142)

This function has been shown in Fig. 6

Figure 6: A plot for Eq. 142 for $\tau=2^{n}$, $n=0,1,2,3,4,5,6,7$ (from top
to bottom).

##  Appendix F Differential Equation

Here we present the connection of the Salpeter equation to the Klein-Gordon
differential equation. By differentiating the Eq. 11 with respect to $x$ and
$t$ One easily finds ($x,t\neq 0$)

$\left(-\partial_{x}^{2}+\partial_{t}^{2}+m^{2}\right)G\left(x,t,m\right)=0,$
(143)

which can be understood using the fact that $H\equiv i\partial_{t}$ and
$p\equiv-i\partial_{x}$. Above, we considered a complex time $t+i0^{+}$ for
control the integral limits. In an extreme situation, one may Wick-rotate the
time, so that $\tau=it$ describes a new coordination which is treated like $x$
as follows ($x,\tau\neq 0$)

$\left(\nabla^{2}-m^{2}\right)G^{(\tau)}=0,$ (144)

which is a homogeneous Helmholtz equation in two dimensions for $x,\tau\neq
0$, and $G^{(\tau)}(x,\tau,m)\equiv G(x,-i\tau,m)$, and
$\nabla^{2}\equiv\partial_{x}^{2}+\partial_{\tau}^{2}$. Note that this
equation respects the scaling relation given in Eq.77. The resulting integral
after the Wick rotation describes also the relativistic probability density of
stochastic processes. In terms of this quantity, we have

$G^{(\tau)}(x,\tau,m=0)=\frac{\tau}{\pi(x^{2}+\tau^{2})}+\frac{1}{2}\delta(x-i\tau)+\frac{1}{2}\delta(x+i\tau),$
(145)

i.e. when $x$ is real, only the first term remains, which is the Cauchy
distribution. Using the identity

$\lim_{\epsilon\to
0}\frac{1}{\pi}\frac{\epsilon}{x^{2}+\epsilon^{2}}=\delta(x),$ (146)

we find that

$\begin{split}\nabla^{2}G^{(\tau)}(x,\tau,m=0)&=-\frac{8\tau}{x^{2}+\tau^{2}}\lim_{\epsilon\to
0}\frac{1}{\pi}\frac{\epsilon^{2}}{(x^{2}+\tau^{2}+\epsilon^{2})^{2}}\\\
&=-4f(x,\tau)\delta(x)\delta(\tau),\end{split}$ (147)

where

$\begin{split}f(x,\tau)\equiv&\frac{2\tau}{(x^{2}+\tau^{2})}=\frac{1}{\tau+ix}+\frac{1}{\tau-
ix}\\\ &=2\pi G^{(\tau)}(x,\tau,m=0)\end{split}$ (148)

We immediately conclude that the same happens for general $m$

$\begin{split}\left(\nabla^{2}-m^{2}\right)G^{(\tau)}(x,\tau,m)=-4f(x,\tau)\delta(x)\delta(\tau)\end{split}$
(149)

Defining $r\equiv\sqrt{x^{2}+\tau^{2}}$ and
$\theta\equiv\arctan\frac{x}{\tau}$, we rewrite it as ($r,\theta\neq 0$)

$G_{rr}^{(\tau)}+\frac{1}{r}G_{r}^{(\tau)}+\frac{1}{r^{2}}G_{\theta\theta}^{(\tau)}-m^{2}G^{(\tau)}=0,$
(150)

where $G_{rr}^{(\tau)}$ means double derivative with respect to $r$ and etc.
Then, using the method of separation of variables $G^{(\tau)}\equiv
R(r)\Theta(\theta)$, we find

$\begin{split}&\Theta^{\prime\prime}+n^{2}\Theta=0\\\
&r^{2}R^{\prime\prime}+rR^{\prime}-\left((mr)^{2}+n^{2}\right)R=0,\end{split}$
(151)

from which one finds

$\Theta=a_{n}\sin(n\theta)+b_{n}\cos(n\theta)$ (152)

and $R$ is a modified Bessel function of the second kind

$R(r)=K_{n}(mr)=\int_{0}^{\infty}\exp[-mr\cosh t]\cosh(nt)\text{d}t.$ (153)

Note that we have two independent radial solutions and the only one which is
infinity at $r=0$ and regular at $r\to\infty$ is this term. We find eventually

$G^{(\tau)}=\frac{1}{2\pi}\sum_{n}K_{n}(mr)\left(a_{n}\sin(n\theta)+b_{n}\cos(n\theta)\right)$
(154)

where $a_{n}$ and $b_{n}$ should be found using Eq. 145 for $m=0$. Given that

$\begin{split}&K_{n}(y)\to\sqrt{\frac{\pi}{2}}\frac{\exp(-y)}{\sqrt{y}},\ \
y\to\infty\\\ &K_{0}(y)\to-\gamma-\ln\frac{y}{2}+O(y^{2})\ \ y\to 0\\\
&K_{1}(y)\to\frac{1}{y}+\frac{y}{2}\ln\frac{y}{2}+O(y),\ \ y\to 0\\\
&K_{n}(y)\to(-1)^{n-1}\left(\frac{y}{2}\right)^{n}\ln\frac{y}{2}\left(\frac{1}{n!}+O(y^{2})\right),\
n\geq 2,y\to 0.\end{split}$ (155)

Comparing this result with Eq. 138, we notice that the singularity of the
integral for $r\to 0$ (equivalent to $x=\pm t$) is by $1/r$ and $r\ln r$, so
that the only term that contributes in the summation is $n=1$. More precisely,
Eq. 138 gives us

$\begin{split}G\to&\frac{m}{2}\left[\frac{1}{mr}e^{-i\theta}+\frac{1}{2}\frac{mre^{i\theta}}{2}\ln(mr)\right]\\\
+&\frac{m}{2}\left[\frac{1}{mr}e^{i\theta}+\frac{1}{2}\frac{mre^{-i\theta}}{2}\ln(mr)\right]\\\
&=m\left[\frac{1}{mr}+\frac{1}{2}\frac{mr}{2}\ln(mr)\right]\cos\theta\end{split}$
(156)

Then, we find that $a_{1}=m$, and $b_{1}=0$, and all the other $a_{n},b_{n}$
($n\neq 1$) are zero. The following integral identity confirms the above
result

$\begin{split}\int_{0}^{\infty}&{dx\,\,\text{exp}(-\alpha\sqrt{x^{2}+\beta^{2}})\cos\gamma
x}\\\
&=\frac{\alpha\beta}{\sqrt{\alpha^{2}+\gamma^{2}}}K_{1}(\beta\sqrt{\alpha^{2}+\gamma^{2}}),\text{Re}[\alpha]>0,\text{Re}[\beta]>0.\end{split}$
(157)

Using this, we find that

$\begin{split}G(x,t)&=\frac{imte^{imt}}{\pi\sqrt{x^{2}-t^{2}}}K_{1}\left(m\sqrt{x^{2}-t^{2}}\right),\
x^{2}>t^{2}\\\
&=-\frac{mte^{imt}}{2\sqrt{t^{2}-x^{2}}}H_{1}^{(1)}\left(m\sqrt{t^{2}-x^{2}}\right),\
x^{2}<t^{2},\end{split}$ (158)

where $H_{1}^{(1)}$ is a Hankel function of first kind. Using new complex
variables $z\equiv\tau+ix$ and $\bar{z}\equiv\tau-ix$, we can simplify the
notation. First we note that
$4\partial\bar{\partial}=\partial_{\tau}^{2}+\partial_{x}^{2}$, where
$\partial\equiv\partial_{z}$, and $\bar{\partial}\equiv\partial_{\bar{z}}$.
Equation 149 gives rise to

$\left(4\partial\bar{\partial}-m^{2}\right)G^{(\tau)}(z,\bar{z},m)=-2\left(\frac{1}{z}+\frac{1}{\bar{z}}\right)\delta^{2}(z)$
(159)

where $\delta^{2}(z)\equiv\delta(x)\delta(t)$, and we have used
$\partial\frac{1}{\bar{z}}=\bar{\partial}\frac{1}{z}=\delta^{2}(z)$. For the
$m=0$ case, $G^{(\tau)}$ is a summation of the holomorphic and anti-
holomorphic terms, i.e.

$G^{(\tau)}(z,\bar{z},m=0)=G_{1}(z)+G_{2}(\bar{z})=\frac{1}{\pi}\left[\frac{1}{z}+\frac{1}{\bar{z}}\right]$
(160)

When $m\neq 0$ however, it is not decomposable to these two parts. In this
case, Eq. 157 tells us that

$G^{(\tau)}(z,\bar{z},m)\propto\frac{1}{\sqrt{z\bar{z}}}K_{1}(m\sqrt{z\bar{z}}).$
(161)