]thmTheorem[section]

# Martingale drift of Langevin dynamics and classical canonical spin
statistics

Ken Sekimoto  Laboratoire Matière et Systèmes Complexes, UMR CNRS 7057,
Université Paris Cité,
10 Rue Alice Domon et Léonie Duquet, 75013, Paris, France  Laboratoire
Gulliver, UMR CNRS 7083, ESPCI Paris, Université PSL
10 rue Vauquelin, 75005, Paris, France<EMAIL_ADDRESS>
###### Abstract

A martingale is a stochastic process that encodes a kind of fairness or
unbias, which is associated with a reference process. Here we show that, if
the reference process $x_{t}$ evolves according to the Langevin equation with
drift $a(x),$ and if $a(x_{t})$ is martingale, then its amplitude is the
Langevin function, which originally described the canonical response of a
single classical Heisenberg spin under static field. Furthermore, the
asymptotic limit of $x_{t}/t$ obeys the ensemble statistics of such Heisenberg
spin.

##  I Backgrounds

Paul Langevin derived in 1905 what we call nowadays the Langevin function [1]
that gives the canonical equilibrium response of paramagnetism under static
magnetic field. It was based on the Boltzmann-Einstein statistics and
explained the Curie’s law, which Pierre Curie has established experimentally
ten years before [2]. In the modern language, a classical 3D Heisenberg spin
under a uniform magnetic field along an axis undergoes the thermal fluctuation
according to the Boltzmann weight, or the relative probability,
$e^{x\cos\theta},$ where $\theta$ is the polar angle of the spin relative to
the field axis and $x$ is the strength of the magnetic field scaled by the
temperature and magnetic moment of the spin. The mean polarisation of the spin
is given by the average of $\cos\theta$ weighted by the above probability and
integrated over the solid angles. The result is the Langevin function,
$\coth{x}-{\frac{1}{x}}$ [1]. That it behaves like ${\frac{1}{3}}x$ for small
$|x|$ explained the Curie’s law.

Three years after Albert Einstein [3] had formulated Guoy’s qualitative idea
of Brownian motion [4, 5], P. Langevin introduced in 1908 the first stochastic
differential equation that is nowadays called the Langevin equation [6]. In
the modern language a simple and generic form of Langevin equation in
$d$-dimension reads

$\frac{d{\vec{x}_{t}}}{dt}={\vec{a}}({\vec{x}_{t}})+\vec{\xi}_{t},$ (1)

which contains the drift, ${\vec{a}}({\vec{x}}_{t}),$ and the white Gaussian
noise, $\vec{\xi},$ with the zero mean and the unit diagonal covariance,
$\langle\vec{\xi}_{t},\vec{\xi}_{s}\rangle=\bm{1}\delta(t-s),$ where $\bm{1}$
is the unit tensor. While we nowadays know how this equation admits the first
[7, 8] and second [9, 10] laws of thermodynamics, it was Langevin who
revolutionarized the notion of the evolution equation as a mapping between the
path ensembles. In the 1900’s, however, no link between the Langevin function
and Langevin equation has been known, to the author’s knowledge.

It is around 1940 that Jean Viller and Doob modernized the historically old
notion of martingale as a powerful concept of the modern probability theory,
see for example the Chap.1 of [11] for a historical review. The martingale is
the mathematical expression of the idea of the fairness for the future: The
continuous-time stochastic process $\vec{Y}_{t}$ is said to be martingale
associated with the stochastic process $\\{\vec{X}_{s}\\}$ if (i)
$\vec{Y}_{t}$ is causally determined by $\\{\vec{X}_{s}\\}_{0\leq s\leq t}$
and (ii) its conditional expectation for future shows the fairness:

$\langle\vec{Y}_{t}|\\{\vec{X}_{u}\\}_{0\leq u\leq
s}\rangle=\vec{Y}_{s}\quad\mbox{for $\forall t\geq s$},$ (2)

let alone the other rigorous mathematical conditions. Here
$\\{\vec{X}_{u}\\}_{0\leq u\leq s}$ denotes the history of $\vec{X}_{u}$ over
the peiode $0\leq u\leq s$ and $\langle R|{\it cond}\rangle$ means the
expectation of the random variable $R$ under the condition, cond. For the
later convenience we introduce the differential version which follows from (2)
applied to $Y_{t+dt}$ as well as to $Y_{t}$:

$\langle d\vec{Y}_{t}|\\{\vec{X}_{u}\\}_{0\leq u\leq s}\rangle=0\mbox{ for
}\forall t\geq s,$ (3)

As the probability theory provides with many useful theorems about martingale,
efforts have often been made to convert a reference process of interest into a
martingale process to draw new properties of $X_{t}.$ The population dynamics
or mathematical finances have developed that approach since longtime [11]. In
the stochastic thermodynamics people recently recognized that the
exponentiated entropy production is the martingale process of the category
called path-probability ratio or Radon-Nikodym density process [12, 13].

A distinct type of martingale other than the path-probability ratio has also
been found in what we call Progressive Quenching [14, 15, 16, 17]. Starting
from an equilibrium many spin system, we fix, one after another, the spins at
the instantaneous state they take. Through such an unbiased protocol, the
martingale process was found not in the progressively fixed spins, but in the
mean spin to be fixed. Those already fixed spins influence on the latter
through the molecular field. We have coined a name “hidden martingale” for it
because, in the continuous-time counterpart, it is not ${\vec{x}}_{t}$ in (1)
but the drift ${\vec{a}}({\vec{x}}_{t})$ that should be martingale — This is
the very starting point of the present study. Below, in Sec.II we show that,
if the drift ${\vec{a}}({\vec{x}}_{t})$ in (1) is martingale associated with
the process (1) itself, then the drift is the Langevin function111The contents
essentially up to Sec.II is accessible as arXiv:2305.04976v1.. The emergence
of a spin-related function is surprising because the Langevin equation (1)
works in a non-compact space. However, a further surprise in Sec.III is that
the canonical statistics of such spin emerges in the ensemble of the long-time
asymptotes of $x_{t}/t.$

##  II Langevin function as self-harmonic drift

### IIA. Family of harmonic functions associated with a drift

When the process $\vec{Y}_{t}=\vec{h}(\vec{X}_{t})$ is martingale associated
with the process $\\{\vec{X}_{s}\\}$, the function $\vec{h}$ is said to be
harmonic [11]. In particular, if the process
$\\{\vec{X}_{s}\\}=\\{{\vec{x}}_{t}\\}$ is generated by (1) and that
$\vec{h}({\vec{x}}_{t})$ takes the value in the same space as ${\vec{x}}_{t}$
does, the condition for the harmonicity reads

$({\vec{a}}\cdot\nabla)\vec{h}+{\frac{1}{2}}\Delta\vec{h}=0.$ (4)

To see this it suffices to develop $d\vec{h}({\vec{x}}_{t})$ up to the order
$\mathcal{O}(dt)$, that is,

$d\vec{h}({\vec{x}}_{t})=(d{\vec{x}}_{t}\bullet\nabla)\vec{h}+{\frac{1}{2}}(d{\vec{x}}_{t}\bullet\nabla)^{2}\vec{h},$
(5)

where $\bullet$ is the Itô-type product [11]. We substitute into it the
stochastic differential version of (1),

$d{\vec{x}}_{t}={\vec{a}}({\vec{x}}_{t})dt+d\vec{W}_{t},$ (6)

where $\vec{W}_{t}$ is the $d$-dimensional standardized Brownian motion, i.e.,
the Wiener process, such that $d\vec{W}_{t}d\vec{W}_{t}=\bm{1}dt$ with
$\bm{1}$ being the unit tensor. The imposition of the martingale condition (3)
on $\vec{Y}_{t}=\vec{h}({\vec{x}}_{t})$ then leads to (4). The form (4)
explains the denomination ‘harmonic’ because, in the absence of the drift,
${\vec{a}}\equiv 0,$ the Eq.(4) is the vector Laplace equation.

Given a drift ${\vec{a}}$ we can conceive the family of harmonic functions
associated with the process obeying (1), We shall denote this family by
$\mathcal{H}_{{\vec{a}}}:$

$\mathcal{H}_{{\vec{a}}}={\left\\{{\vec{h}\,\,;\,({\vec{a}}\cdot\nabla)\vec{h}+{\frac{1}{2}}\Delta\vec{h}=0}\right\\}}.$
(7)

### IIB. Langevin function as fixed point

Among the families associated with different drift, the self-referential
condition

${{\vec{a}}^{\,*}}\in\mathcal{H}_{{\vec{a}}^{\,*}}$ (8)

defines the special drift ${\vec{a}}^{\,*}$ as a kind of fixed point, which we
shall call the self-harmonic drift when it exists. More concretely,

$({{\vec{a}}^{\,*}}\\!\cdot\\!\nabla){{\vec{a}}^{\,*}}+{\frac{1}{2}}\Delta{{\vec{a}}^{\,*}}=\vec{0}.$
(9)

This is the key equation of the present paper.

There is possibilities of asymmetric solutions of (9). Apparently the
cylindrical lift-up of solution for the lower-dimensions is an anisotropic
solution. Nevertheless, we here seek for the solutions which are isotropic
with respect to the origin, ${\vec{x}}=0.$ (Note: The process ${\vec{x}}_{t}$
can nevertheless start from any ${\vec{x}}$ off the origin.)

We then assume the form222We use the notation ${\vec{x}}$ rather than
$\vec{r}$ because it is generally not be a spatial position and, moreover,
$L(x)$ is more common than $L(r)$ in the literature.:

${\vec{a}}^{\,*}({\vec{x}})=L_{d}({\|\vec{{x}}\|})\hat{x},$ (10)

where $\hat{x}\equiv{{\vec{x}}}/{{\|\vec{{x}}\|}}$ is the unit vector along
${\vec{x}}$ and the suffix $d$ stands for the space dimension. Then (9)
implies

$L_{d}^{\prime\prime}(x)+2L_{d}(x)L_{d}^{\prime}(x)+\frac{d-1}{x}{\left({L_{d}^{\prime}(x)-\frac{L_{d}(x)}{x}}\right)}=0.$
(11)

In the two-dimensional family of solutions of (11) the invariance under the
similarity transformation, $L_{d}(x)\to\alpha L_{d}(\alpha x)$ provides with
one parameter. If we write the first integral of (11) as

$L_{d}^{\prime}(x)+(L_{d}(x))^{2}+\frac{d-1}{x}L_{d}(x)=\alpha^{2},$ (12)

then we can reduce the problem to the finding of $L_{d}(x)$ for (12) with
$\alpha=1.$ Then in the remaining one-parameter family of solutions for (12),
some analysis using Mathematica${}^{\mbox{\textregistered}}$ tells that, at
least for $d=$2,3 and 4, there is a unique solution which does not diverge at
$x=0.$ For example in $d=3,$ the solution
$L_{3}^{(\beta)}(x)\equiv\coth(x+\beta)-{\frac{1}{x}}$ is regular only with
$\beta=0.$ For $d=1$ we can solve (9) directly so that $\alpha=1$ in (12) and
choose appropriately the origin. Altogether, the solutions thus identified are
as below:

$\displaystyle L_{1}(x)$ $\displaystyle:=$ $\displaystyle\tanh(x)$ (13)
$\displaystyle L_{2}(x)$ $\displaystyle:=$
$\displaystyle\frac{I_{1}(x)}{I_{0}(x)}$ (14) $\displaystyle L_{3}(x)$
$\displaystyle:=$ $\displaystyle\coth{x}-{\frac{1}{x}}$ (15) $\displaystyle
L_{4}(x)$ $\displaystyle:=$
$\displaystyle\frac{x[I_{0}(x)+I_{2}(x)]-2I_{1}(x)}{2xI_{1}(x)},$ (16)

for $x\neq 0$ and $L_{d}(0)=0$ for any dimension, $d,$ where $I_{n}(x)$ are
the $n$-th modified Bessel functions of the first kind. A standard but
inspiring approach to obtain these solution from (12) with $\alpha=1$ is to
use the (Riccati) transformation,

$L_{d}\equiv\frac{Z^{\prime}_{d}}{Z_{d}}=\frac{d}{dx}[\ln Z_{d}(x)],$ (17)

which renders (12) into the linear equation,

$\frac{d^{2}Z_{d}}{dx^{2}}+\frac{d-1}{x}\frac{dZ_{d}}{dx}=Z_{d}.$ (18)

By noticing the radial part of the Laplacian operator on the l.h.s. of (18),
we find

$Z_{d}(\|{\vec{x}}\|)\propto\oint_{\|\hat{S}\|=1}e^{\hat{S}\cdot{\vec{x}}}d\Omega_{S},$
(19)

that is, the partition function for a single classical Heisenberg spin under
the non-dimensionalized external field, ${\vec{x}}.$ Through (17) the drift
is, therefore, the canonical average of the spin $\hat{S}$ under this field.
All $L_{d}(x)$’s are odd in $x$ and their graphs look similar to the simplest
one, $L_{1}(x),$ or the original Langevin function, $L_{3}(x).$ They also have
the unique limit, $\lim_{x\to\infty}L_{d}(x)=1.$ We may call $L_{d}(x)$ the
$d$-dimensional Langevin function. In any case the Langevin equation meets
here with the Langevin functions in the context of the martingale. For
completeness we list below $Z_{d}(x)$ with the numerical coefficient so that
$Z_{d}(0)$ is the surface area of $d$-dimensional unit hyper-sphere:

$\displaystyle(Z_{1},Z_{2},Z_{3},Z_{4})$ (20)
$\displaystyle={\left({2\cosh(x),2\pi
I_{0}(x),4\pi\frac{\sinh(x)}{x},4\pi^{2}\frac{I_{1}(x)}{x}}\right)}.$ (21)

#### Remark on the solution family

— The change of the units of length and time causes the rewriting of (6). If
we introduce $\vec{y}$ and $\tau$ through ${\vec{x}}=\alpha\vec{y}$ and
$t=\alpha^{2}C\tau,$ then (1) with (10) for ${\vec{a}}({\vec{x}}_{t})$ leads
to

$d\vec{y}_{\tau}=C\alpha\,L_{d}(\alpha{\|\vec{y}_{\tau}\|}{})\hat{y}_{\tau}d\tau+\sqrt{C}d\vec{W}_{\tau},$
(22)

where $\hat{y}_{\tau}\equiv\vec{y}_{\tau}/\|\vec{y}_{\tau}\|$ and we have used
the statistical equivalence, $d\vec{W}_{\kappa
t}\simeq\sqrt{\kappa}\,d\vec{W}_{t},$ for the standard $d$-dimensional Wiener
process. In (22) we see that, apart from the flexibility of the self-harmonic
function, $\alpha L_{d}(\alpha{\|\vec{y}_{\tau}\|}{}),$ there is a specific
scale relationship between the drift and diffusion for the drift to be
martingale. See Appendix A for more argument.

##  III Statistical property of the asymptotic limits : “microscope”

In general the Langevin functions $L_{d}(x)$ can appear in being unrelated
with the canonical statistics of a Heisenberg spin. For example, $L_{2}(x)$
has appeared as the 1+1-dimensional nucleation-controlled polymer crystal
growth rate [20, 21], with no relevance to the XY-spin. In this section we
show that the statistics of a Heisenberg spin indeed appears in the asymptotic
behavior of the Langevin dynamics with self-harmonic drift.

### IIIA. Convergence of the self-harmonic drift and asymptotic behavior of
Langevin process

Since the velocity $d{\vec{x}}_{t}/dt$ is on the average oriented along
${\vec{x}}_{t}$ (see (10)) and that $L_{d}(\|{\vec{x}}_{t}\|)$ is nonnegative
for any nonzero ${\vec{x}}_{t},$ we may expect that ${\vec{x}}_{t}$ most
probably grows unboundedly for a long time. However, its further detail is not
clear at first glance. We, therefore, focus first on the evolution of the
drift, ${\vec{a}}^{\,*}({\vec{x}}_{t}).$ Because of the martingale condition
(9) the development of $d{\vec{a}}^{\,*}({\vec{x}}_{t})$ contains no (drift)
term for $\mathcal{O}(dt)$ (see (5)), and the result reads

$d{\vec{a}}^{\,*}({\vec{x}}_{t})=d\vec{W}_{t}\bullet(\nabla{\vec{a}}^{\,*}),$
(23)
$\nabla{\vec{a}}^{\,*}=\frac{dL_{d}(\chi)}{d\chi}\hat{a}^{*}\hat{a}^{*}+\frac{L_{d}(\chi)}{\chi}(\bm{1}-\hat{a}^{*}\hat{a}^{*}),$
(24)

where $\hat{a}^{*}\equiv{\vec{a}}^{\,*}/\|{\vec{a}}^{\,*}\|$ and $\chi$ is
inversely determined so that $L_{d}(\chi)=\|{\vec{a}}^{\,*}\|.$ Thus
$(\nabla{\vec{a}}^{\,*})$ is a function of ${\vec{a}}^{\,*}.$ The process
generated by $d{\vec{a}}_{t}=d\vec{W}_{t}\bullet{\sf M}({\vec{a}}_{t})$ with
${\sf M}({\vec{a}})$ being any rank-2 tensor as function of ${\vec{a}}$ is
martingale associated with the Wiener process, $\vec{W}_{s}$ $(0\leq
s<\infty).$ The self-harmonic drift ${\vec{a}}^{\,*}$ defined above, however,
has several particular properties on top of it, which we will discuss below.

We notice that ${\vec{a}}^{\,*}({\vec{x}}_{t})$ of each realisation converges
for $t\to\infty$ and the limit
${\vec{a}}^{\,*}_{\infty}\equiv\lim_{t\to\infty}{\vec{a}}^{\,*}({\vec{x}}_{t})$
is on the hypersphere, $\|{\vec{a}}^{\,*}_{\infty}\|=1.$ The case of $d=1$ is
particularly simple: Eqs.(23) and (24) reads

$da^{*}_{t}=(1-{a^{*}_{t}}^{2})\bullet dW_{t}$ (25)

because $(\tanh x)^{\prime}=1-\tanh^{2}x.$ The Fokker-Planck (FP) equation
equivalent to (25) has the potential, $\ln[(1-{a^{*}}^{2})^{2}]$ and the
diffusion coefficient ${\frac{1}{2}}(1-{a^{*}}^{2})^{2}.$ The variable
${\vec{a}}^{\,*}_{t}$ should be pushed towards $\pm 1$ and stuck there. In
fact the singular densities $\delta(a_{\infty}^{*}\mp 1)$ or their linear
combinations are all stationary solutions of that FP equation.

While the behavior of ${\vec{a}}^{\,*}({\vec{x}}_{t})$ for $d>1$ is less clear
than for $d=1,$ the fact that $\nabla a^{*}$ vanishes on the hyper-sphere
$\|{\vec{a}}^{\,*}\|=1$ (i.e. at $\chi=+\infty$) indicates that the hyper-
sphere is the natural absorbing boundary for the evolution of
${\vec{a}}^{\,*}({\vec{x}}_{t}).$ More formal and general argument for the
existence of limit is provided by the convergence theorem of (sub)martingale
(see, for e.g., §4.1.5 of [11]). According to this theorem, because
${\vec{a}}^{\,*}({\vec{x}})$ is bounded, i.e., $\|{\vec{a}}^{\,*}\|\leq 1,$
the martingale process ${\vec{a}}^{\,*}({\vec{x}}_{t})$ should have a limit;

${\vec{a}}^{\,*}_{\infty}\equiv\lim_{t\to\infty}{\vec{a}}^{\,*}({\vec{x}}_{t}).$
(26)

On the other hand, if the limit ${\vec{a}}^{\,*}_{\infty}$ were not on the
natural boundary, the noise for ${\vec{x}}_{t}$ could still relocate it.
Therefore, $\|{\vec{a}}^{\,*}_{\infty}\|=1$ is concluded.

Next we consider the trajectory of ${\vec{x}}_{t}$ that leads to
$\lim_{t\to\infty}{\vec{a}}^{\,*}({\vec{x}}_{t})={\vec{a}}^{\,*}_{\infty}$ ?
Fig.1 shows three samples in $d=2,$ having started from the origin,
${\vec{x}}_{0}=0.$ The left figure is the close-up view of the initial part of
the trajectories and the right figure shows the entire trajectory up to $t=16$
including the left figure333The $0_{+}$ is a numerical artificial
infinitesimal, being used to avoid the removable singularity of $L_{d}(x)$ at
$x=0$. In the early stage with $\|{\vec{x}}_{t}\|\lesssim 1,$ the noise term
$\vec{\xi}_{t}$ dominates over the drift, then there is a slow cross-over to
the long-time ballistic behavior. This observation together with (1) or (6),
we understand that ${\vec{x}}_{t}$ becomes asymptotically ballistic in the
sense that

$\lim_{t\to\infty}\frac{{\vec{x}}_{t}}{t}={\vec{a}}^{\,*}_{\infty},$ (27)

which also imply $\lim_{t\to\infty}\frac{\|{\vec{x}}_{t}\|}{t}=1.$ In Appendix
B we show (27) in more detail. Note that each realisation of trajectory ends
up with a particular orientation of ${\vec{a}}^{\,*}_{\infty}.$

Figure 1:  2D trajectories from ${\vec{x}}_{0}=(0_{+},0_{+}).$ (Left) Close-up
view near the origin. (Right) Whole trajectories up to $t=16.$

### IIIB. Distribution of trajectories at $t\to\infty$ as “spin” microscope

Each trajectory of ${\vec{x}}_{t}$ or of ${\vec{a}}^{\,*}({\vec{x}}_{t})$ has
a the limiting value ${\vec{a}}^{\,*}_{\infty}$ as random variable and, by the
definition of self-harmonicity, we have the martingality:

$\langle{\vec{a}}^{\,*}_{\infty}|{\vec{x}}_{0}\rangle={\vec{a}}^{\,*}({\vec{x}}_{0}).$
(28)

This is a constraint on the statistics of ${\vec{a}}^{\,*}_{\infty}.$ However,
it is only for $d=1$ that the probability of realizing $a_{\infty}^{*}=\pm 1$
and the initial data $a^{*}(x_{0})$ are trivially related by (28),

$Prob(a_{\infty}^{*}=s)=\frac{1+s\,a^{*}(x_{0})}{2},\quad s=\pm 1.$

By contrast, the statistics of ${\vec{a}}^{\,*}_{\infty}$ on the
$d$($>1$)-dimensional hyper-sphere surface is a by no means trivial. Fig.2
gives an qualitative idea of how the trajectories of ${\vec{x}}_{t}$ up to
$t=40$ depends the initial data, ${\vec{x}}_{0}.$ Roughly, the larger is
$\|x_{0}\|,$ the more polarized is the orientation of ${\vec{x}}_{t}.$

(a) $x_{0}=0$

(b) $x_{0}=1$

(c) $x_{0}=2$

(d) $x_{0}=4$

Figure 2: 2D trajectories from ${\vec{x}}_{0}=(x_{0},0_{+})$ with (a)
$x_{0}=0$, (b) $x_{0}=1$, (c) $x_{0}=2$ and (d) $x_{0}=4.$ In all cases the
duration is up to $t=40$ and we used the identical set of noise histories for
$\vec{\xi}_{t}$ in (1)

We shall parameterize the orientation of ${\vec{a}}^{\,*}_{\infty}$ by its
orthogonal projection,
$\cos\theta\equiv{\vec{a}}^{\,*}_{\infty}\cdot({\vec{x}}_{0}/\|{\vec{x}}_{0}\|),$
onto the axis along ${\vec{x}}_{0},$ where $0\leq\theta\leq\pi.$ When
${\vec{x}}_{0}=0,$ the distribution of ${\vec{a}}^{\,*}_{\infty}$ through (27)
must be isotropic and the cumulative probability distribution of $\cos\theta$
should rigorously obey
$Prob^{\rm(can)}(\cos\theta<\chi)={\frac{1}{\pi}}\int_{-1}^{\chi}\frac{d\xi}{\sqrt{1-\xi^{2}}}=1-{\frac{1}{\pi}}\arccos\chi.$

Figure 3: The 2D orientational distributions of
$\lim_{t\to\infty}(\hbox{}{{\vec{x}}_{t}}/{t})={\vec{a}}^{\,*}_{\infty}$ for
the different starting points, $x_{0}=0,1,2,$ and $4,$ from top to bottom. We
represent the orientation of ${\vec{a}}^{\,*}_{\infty}$ using $\cos\theta,$
where $\theta$ is the angle between ${\vec{a}}^{\,*}_{\infty}$ and
${\vec{x}}_{0}=(x_{0},0),$ and we show the distributions by the empirical
cumulative probability. Numerically the ensemble of ${\vec{a}}^{\,*}_{\infty}$
is approximated by 3000 realisations of ${\vec{x}}_{t}/t|_{t=40}.$ The red
dashed curves represent the canonical equilibrium distribution of a single
unitary spin under the non-dimensionalized field ${\vec{x}}_{0},$ see (29) in
the main text.

In Fig.3 the top dashed curve (in red) represents this formula. The train of
dots (in blue) along this curve are the results obtained from the numerical
data over 3000 trajectories, where
${\vec{a}}^{\,*}_{\infty}=\lim_{t\to\infty}\frac{{\vec{x}}_{t}}{t}$ is
approximated by $\left.\frac{{\vec{x}}_{t}}{t}\right|_{t=40}.$ The deviations
from the theoretical curve shows the errors due to the finiteness of sampling.
The other three trains of blue dots in Fig.3 represent the numerically
obtained cumulative probabilities $\cos\theta$ for $x_{0}=1,2$ and $4,$
respectively, in the descending order. (Note: We took the cumulative
probability because it can be empirically reconstructed just by plotting the
normalized rank $r/3000$ vs the corresponding value of $\cos\theta_{r},$
without any binning or smoothing.)

Quite surprisingly, the results almost surely obey the canonical statistics of
a Heisenberg spin under the non-dimensonalized field,
${\vec{x}}_{0}=(x_{0},0),$ which leads to the cumulated probability,

$Prob^{\rm(can)}(\cos\theta<\chi|{\vec{x}}_{0})={\frac{1}{\pi
I_{0}(x_{0})}}\int_{-1}^{\chi}\frac{e^{x_{0}\xi}d\xi}{\sqrt{1-\xi^{2}}}.$ (29)

This formula is shown by the dashed curves (in red) in Fig.3. Thus the thermal
distribution of the (fictitious) spin orientations that would yield the
expectation ${\vec{a}}^{\,*}({\vec{x}}_{0})$ has been mapped into the static
distribution of the “spin,” ${\vec{a}}^{\,*}_{\infty},$ characterizing the
trajectory of ${\vec{x}}_{t}$ in the limit of $t\to\infty.$ Being consistent
with (28), our proposition for the probability density
$\rho{\left({{\vec{a}}^{\,*}_{\infty}}\right)}$ is, therefore,

$\rho{\left({{\vec{a}}^{\,*}_{\infty}}\right)}=\frac{e^{{\vec{x}}_{0}\cdot{\vec{a}}^{\,*}_{\infty}}}{Z_{d}(\|{\vec{x}}_{0}\|)},\qquad{\vec{a}}^{\,*}_{\infty}\in
S^{d-1},$ (30)

where ${\vec{a}}^{\,*}_{\infty}$ means
$\lim_{t\to\infty}\frac{{\vec{x}}_{t}}{t}$ and $S^{d-1}$ is the
$(d-1)$-surface of the $d$-dimensional unit sphere. This mapping from
${\vec{x}}_{0}$ to the distribution of
$\lim_{t\to\infty}\frac{{\vec{x}}_{t}}{t}$ serves as a microscope, which
allows to assess ${\vec{x}}_{0}$ through the repeated measurements of
${\vec{x}}_{t}$ at sufficiently large $t.$ Notice that the mapping,
${\vec{x}}_{0}\mapsto Prob^{\rm(can)}(\cos\theta<\chi|{\vec{x}}_{0}),$ is
asymptotically independent of the micro-macro ratios such as
$\|{\vec{x}}_{0}\|/\|{\vec{x}}_{t}\|.$ This is in contrast to the case of
assessing the initial position ${\vec{x}}_{0}$ of a Brownian particle from its
position at a later time, ${\vec{x}}_{t}.$ In the latter case the signal-to-
background ratio lessens with time. See Appendix C for more explanation.

##  IV Discussion

The present work may bring two things. First, we have found a connection
between the Langevin equation and the Langevin functions through the
martingale process. Usually the relationship between the Langevin equation and
canonical distribution is through the Einstein’s relation by which the
stationary state of the Langevin equation become canonical (if the drift has a
potential). In the present case, however, the evolution of neither
${\vec{x}}_{t}$ nor ${\vec{a}}({\vec{x}}_{t})$ has regular stationary density;
${\vec{x}}_{t}$ grows unboundedly while ${\vec{a}}({\vec{x}}_{t})$ is trapped
asymptotically at a point on the surface of unit hyper-sphere. Nevertheless,
the canonical spin statistics, which could give the Langevin function as
response, emerges in the distribution of asymptotic limit of
${\vec{x}}_{t}/t.$ In Appendix D we present a qualitative parable to give a
glimpse of our results.

Secondly, our report is another example in which the self-referential
condition brings us a non-trivial outcome, like the celebrated diagonal
arguments by Cantor and by Gödel, the fractals, to mention but a few. Usually
the martingale process is constructed from a reference process so that we can
apply useful theorems of the martingale theory to the former, and eventually,
we gain insights about the reference process. The present setup asks, instead,
what type of reference process posesses the self-referential martingale
property.

The underlying physics by which the martingale brought the Langevin function
is unknown. That the linearising transformation for the Riccati equation,
$L_{d}(x)=Z_{d}^{\prime}(x)/Z_{d}(x),$ takes the form of the canonical
response might be a clue to this question. Also unknown is how the martingale
constraint transmits the initial molecular field ${\vec{x}}_{0}$ to the
asymptotic distribution of “spin,” ${\vec{a}}^{\,*}_{\infty}.$ In short the
present findings bring us questions, rather than answers, about the relevance
of martingale in physics beyond being a mere mathematical tool. We, therefore,
would like to conclude our report by several open questions: Can the space-
time harmonic function (see, for e.g., §3.2.5.3 of [11]) be also a fixed
point? Can the microscope idea developed below (30) be generalizable to non-
martingale drift ? Are there links to the stochastic thermodynamics [23, 24,
25]? Can there be a parallel framework for the state space having topologies
other than Euclidean one? Will there be a quantum counterpart? Does the self-
harmonic drift optimize some physical entity or information ?

###### Acknowledgements.

The author acknowledges Charles Moslonka for the regular discussion on the
martingale and the Progressive Quenching. The author gratefully appreciates
the helpful comments on the early draft by Matteo Polettini and the insightful
discussion with Felix Ritort.

##  Appendix A Family of self-harmonic drifts obtained by rescaling

For different interpretations of ${\vec{x}}_{t}$ the self-harmonicity of the
drift implies different constraints between the drift and noise terms. For
example, when $\vec{y}_{t}$ represents a position of a particle in a heat bath
of temperature $T$ and friction coefficient $\gamma,$ the choice
$(\alpha,C)=(1,{2\mathcal{D}})$ with $\mathcal{D}=\frac{{k_{\rm B}T}}{\gamma}$
in (22) leads to

$d\vec{y}_{\tau}=\frac{2{k_{\rm
B}T}}{\gamma}L_{d}({\|\vec{y}_{\tau}\|}{})\hat{y}_{\tau}d\tau+\sqrt{\frac{2{k_{\rm
B}T}}{\gamma}}\,d\vec{W}_{\tau},$ (31)

Since $L_{d}({\|\vec{y}_{\tau}\|}{})\hat{y}_{\tau}=\nabla\ln
Z_{d}(\|\vec{y}\|),$ this equation describes the motion in the potential,
$(-2{k_{\rm B}T}\ln Z_{d}(\|\vec{y}\|)).$ To assure the self-harmonicity the
energy scale of the potential is thus constrained.

Another case to apply (22) is through an analogy with the Progressive
Quenching of $d$-dimensional Heisenberg spins on a complete network [14, 15,
16, 17], which we have mentioned in the last paragraph of Sec.I. In this
context the appropriate choice is $(\alpha,C)=(\frac{J}{{k_{\rm
B}T}},\frac{{k_{\rm B}T}}{J}),$ which brings (22) into the following:

$d\vec{y}_{\tau}=L_{d}{\left({\frac{J\|\vec{y}_{\tau}\|}{{k_{\rm
B}T}}}\right)}\hat{y}_{\tau}d\tau+\sqrt{\frac{{k_{\rm
B}T}}{J}}\,d\vec{W}_{\tau}.$ (32)

In this picture we regard $L_{d}{\left({\hbox{}{J\|\vec{y}_{\tau}\|}/{{k_{\rm
B}T}}}\right)}\hat{y}_{\tau}$ on the r.h.s. of (32) as the equilibrium mean of
a Heisenberg spin under the molecular field $J\vec{y}_{\tau}$ with $J$ being a
spin-spin coupling constant and ${k_{\rm B}T}$ being the temperature. Then a
new fragment of spin $d\vec{y}_{\tau}$ appears at every interval $d\tau$ and
joins the quenched part, $\vec{y}_{\tau}.$ While the conditional mean of
$d\vec{y}_{\tau}$ follows the canonical statistics, there is the Wiener noise
$\propto d\vec{W}_{\tau}$ on top of it. To assure the self-harmonicity the
amplitude of the noise must be of $\sqrt{\hbox{}{{k_{\rm B}T}}/{J}}.$

##  Appendix B Long-time convergence of ${\vec{x}}_{t}/t$

In the main text we see that, for every particular trajectory of
${\vec{a}}^{\,*}({\vec{x}}_{t})$ starting from $x_{t}=x_{0},$ there is a
limit, ${\vec{a}}^{\,*}_{\infty}$ with $\|{\vec{a}}^{\,*}_{\infty}\|=1.$ Below
we will see that the ratio, ${{\vec{x}}_{t}}/{t},$ also converges to
${\vec{a}}^{\,*}_{\infty}.$

Once we admit the existence of the limit ${\vec{a}}^{\,*}_{\infty},$ there
exists a $u$ such that
$\|{\vec{a}}^{\,*}(x_{t})-{\vec{a}}^{\,*}_{\infty}\|<\epsilon$ for $\forall
t\geq u$ for a given constant $\epsilon(>0).$ When we integrate (6) over the
interval $[u,t]$ as

${\vec{x}}_{t}={\vec{x}}_{u}+\int_{u}^{t}{\vec{a}}^{\,*}({\vec{x}}_{s})ds+\int_{s=u}^{s=t}dW_{s},$
(33)

we see that ${\vec{x}}_{u}/t\sim t^{-1}$ and $\int_{s=u}^{s=t}dW_{s}/t\sim
t^{-1/2}$ for $t\to\infty.$ We, therefore, focus on
$\int_{u}^{t}{\vec{a}}^{\,*}({\vec{x}}_{s})ds,$ which reads

$\int_{u}^{t}{\vec{a}}^{\,*}({\vec{x}}_{s})ds=(t-u){\vec{a}}^{\,*}_{\infty}+\int_{u}^{t}[{\vec{a}}^{\,*}({\vec{x}}_{s})-{\vec{a}}^{\,*}_{\infty}]ds.$

While the magnitude of the integrand on the r.h.s. is already bounded by
$\epsilon,$ we expect the integrand to decay almost like $s^{-1}$ because for
large ${\vec{x}}_{s}$ the Langevin equation (1) means roughly
$d{\vec{x}}_{s}/ds\simeq(1-\mathcal{O}(\|{\vec{x}}_{s}\|))\hat{x}_{s}+d\vec{W}_{t},$
where we noticed $L_{d}(x)=1-\mathcal{O}(x^{-1})$ for $x\to\infty.$ Therefore,
the time integral on the r.h.s. is dominated by $(t-u)\epsilon.$ Dividing each
term in (33) by $t,$ we reach the claimed result;

$\lim_{t\to\infty}\frac{{\vec{x}}_{t}}{t}={\vec{a}}^{\,*}_{\infty},\mbox{ when
}\lim_{t\to\infty}{\vec{a}}^{\,*}({\vec{x}}_{t})={\vec{a}}^{\,*}_{\infty}.$

##  Appendix C Inferring the initial position of a Brownian particle a
posteriori

When a Brownian particle starts from ${\vec{x}}_{0},$ its position at time
$t,$ that is ${\vec{x}}_{t}={\vec{x}}_{0}+\vec{W}_{t},$ is martingale
associated with the Wiener process $\vec{W}_{t}.$ While the mean
$\langle{\vec{x}}_{t}|{\vec{x}}_{0}\rangle$ remains ${\vec{x}}_{0},$ the
variance $\langle({\vec{x}}_{t}-{\vec{x}}_{0})^{2}|{\vec{x}}_{0}\rangle$ grows
linearly in time. Therefore, unlike the case of self-harmonic drift, the
inference of ${\vec{x}}_{0}$ through the observations of ${\vec{x}}_{t}$ at a
“macroscopic” distance $R(\gg\|{\vec{x}}_{0}\|\geq 0)$ is hard to realize
under the signal-to-background ratio decaying with $R.$ In fact, through the
analysis of the exit problem [26], the detected position on the circle,
$\|{\vec{x}}\|=R,$ obeys the cumulative probability,

$\displaystyle Prob(\cos\theta<\chi)$ $\displaystyle=$
$\displaystyle\frac{1-\epsilon^{2}}{2}{\left({{\frac{1}{\sqrt{1-2\epsilon\chi+\epsilon^{2}}}}-{\frac{1}{1+\epsilon}}}\right)}$
(34) $\displaystyle=$
$\displaystyle\frac{1+\chi}{2}{\left[{1-\frac{3(1-\chi)}{2}\epsilon+\mathcal{O}(\epsilon^{2})}\right]},$
(35)

where $\epsilon\equiv\|{\vec{x}}_{0}\|/R.$ We see that the signal, the term
$-\frac{3(1-\chi)}{2}\epsilon,$ vanishes in the macroscopic limit,
$\epsilon\to 0.$

##  Appendix D Parable of the contents of present work

Here we give a qualitative parable to give a glimpse of the system that we
found and analyzed, in disregarding the logical order and some notations in
the main text.

Suppose there is a $d$-dimensional space with Euclidean metric in which the
“field” is radially oriented outwards and the strength at the position
${\vec{x}}$ is simply ${\vec{x}}.$ A space-ship that travels in this space has
a “compass” made of the $d$-dimensional spin $\hat{S}$ with $\|\hat{S}\|=1$
that feels the field ${\vec{x}}$ in which the space-ship finds itself. The
thermal fluctuation of this spin is, however, so rapid that the navigator of
the space-ship can only observe the thermal average,
$\langle\hat{S}\rangle=L_{d}(\|{\vec{x}}\|)\hat{x},$ where $\hat{x}$ is the
unit vector along ${\vec{x}}$ and $L_{d}(x)$ is the Langevin function that
saturates to unity for $x\to\infty.$ Starting from ${\vec{x}}_{0}$ at time
$t=0,$ the instantaneous speed of the space-ship at position ${\vec{x}}_{t}$
is the sum of the drift which is simply the mean spin,
$\langle\hat{S}_{t}\rangle\,(=L_{d}(\|{\vec{x}}_{t}\|)\hat{x_{t}}),$ and the
standard white-Gaussian noise, $\xi_{t},$ i.e., the time derivative of a
Wiener process $\vec{W}_{t}.$ After a large enough time, the space-ship is far
from the origin and the mean speed ${\vec{x}}_{t}/t$ converges to a unit
vector, which we denote by $\hat{S}_{\infty}.$ The first finding (the
martingale property) is that the ensemble average of the asymptotic mean speed
$\hat{S}_{\infty}$ over the trajectories from ${\vec{x}}_{0}$ is identical to
the initial speed, that is,
$\langle\hat{S}_{\infty}\rangle_{{\vec{x}}_{0}}=\langle\hat{S}_{0}\rangle_{{\vec{x}}_{0}},$
where the r.h.s. is the thermal average of the spin at ${\vec{x}}_{0}.$ The
second finding is that the ensemble distribution of $\hat{S}_{\infty}$ is
identical to the thermal distribution of $\hat{S}_{0}$ at the starting point
${\vec{x}}_{0},$ the level of information which the navigator did could not
access.

## References

  * Langevin [1905] P. Langevin, J. Phys. Theor. Appl. 4, 678 (1905).
  * Curie [1895] P. Curie,  _Propriétés magnétiques des corps à diverses températures_ (Thesis: Gauthier-Villars, Paris (France), 1895).
  * Einstein [1905] A. Einstein, Ann. Phys. (Leipzig) 17, 549 (1905).
  * Gouy [1888] Gouy, Journ. de Phys., 2e série 7, 561 (1888).
  * Gouy [1889] Gouy, Comptes rendus 109, 102 (1889).
  * Langevin [1908] P. Langevin, C.R. Acad. Sci. (Paris) 146, 530 (1908).
  * Sekimoto [1997] K. Sekimoto, J. Phys. Soc. Jpn. 66, 1234 (1997).
  * Sekimoto [1998] K. Sekimoto, Prog. Theor. Phys. Suppl. 130, 17 (1998).
  * Seifert [2005] U. Seifert, Phys. Rev. Lett. 95, 040602 (2005).
  * Seifert [2012] U. Seifert, Reports on Progress in Physics 75, 126001 (2012).
  * Édgar Roldán  _et al._ [2023] Édgar Roldán, I. Neri, R. Chetrite, S. Gupta, S. Pigolotti, F. Jülicher, and K. Sekimoto, Martingales for physicists (2023), arXiv:2210.09983 [cond-mat.stat-mech] .
  * Chetrite and Gupta [2011] R. Chetrite and S. Gupta, J. Stat. Phys. 143, 543 (2011).
  * Neri  _et al._ [2017] I. Neri, E. Roldán, and F. Jülicher, Phys. Rev. X 7, 011019 (2017).
  * Ventéjou and Sekimoto [2018] B. Ventéjou and K. Sekimoto, Phys. Rev. E 97, 062150 (2018).
  * Etienne and Sekimoto [2018] M. Etienne and K. Sekimoto, ACTA PHYSICA POLONICA B 49, 883 (2018).
  * Moslonka and Sekimoto [2020] C. Moslonka and K. Sekimoto, Phys. Rev. E 101, 062139 (2020).
  * Moslonka and Sekimoto [2022] C. Moslonka and K. Sekimoto, Phys. Rev. E 105, 044146 (2022).
  * Note [1] The contents essentially up to Sec.II is accessible as arXiv:2305.04976v1.
  * Note [2] We use the notation ${\mathaccentV{vec}17E{x}}$ rather than $\mathaccentV{vec}17E{r}$ because it is generally not be a spatial position and, moreover, $L(x)$ is more common than $L(r)$ in the literature.
  * Bennett  _et al._ [1981] C. H. Bennett, M. Büttiker, R. Landauer, and H. Thomas, Journal of Statistical Physics 24, 419 (1981).
  * Goldenfeld [1984] N. Goldenfeld, Journal of Physics A: Mathematical and General 17, 2807 (1984).
  * Note [3] The $0_{+}$ is a numerical artificial infinitesimal, being used to avoid the removable singularity of $L_{d}(x)$ at $x=0$.
  * Sekimoto [2010] K. Sekimoto,  _Stochastic Energetics (Lecture Notes in Physics, vol. 799)_ (Springer, 2010).
  * Peliti and Pigolotti [2021] L. Peliti and S. Pigolotti,  _Stochastic Thermodynamics: An Introduction_ (Princeton University Press, 2021).
  * Shiraishi [2023] N. Shiraishi,  _An Introduction to Stochastic Thermodynamics, From Basic to Advanced_ (Springer, 2023).
  * Sekimoto [2021] K. Sekimoto, Derivation of the first passage time distribution for markovian process on discrete network (2021).


∎

11institutetext: Anadijiban Das 22institutetext: Department of Mathematics,
Simon Fraser University, Burnaby, British Columbia, V5A 1S6, Canada
22email<EMAIL_ADDRESS>33institutetext: Rupak Chatterjee 44institutetext: Center
for Quantum Science and Engineering
Department of Physics, Stevens Institute of Technology, Castle Point on the
Hudson, Hoboken, NJ 07030, USA
44email<EMAIL_ADDRESS>
# Discrete Phase Space-Continuous Time Relativistic Klein-Gordon and Dirac
Equations, and a New Non-Singular Yukawa Potential

Anadijiban Das  Rupak Chatterjee

(Received: date / Accepted: date)

###### Abstract

This paper deals with the second quantization of interacting relativistic
Fermionic and Bosonic fields in the arena of discrete phase space and
continuous time. The mathematical formulation involves partial difference
equations. The corresponding Feynman diagrams and a new $S^{\\#}$-matrix
theory is developed. In the special case of proton-proton Møller scattering
via an exchange of a neutral meson, the explicit second order element $\langle
f|S^{\\#}_{(2)}|i\rangle$ is deduced. In the approximation of very low
external three-momenta, a new Yukawa potential is explicitly derived from
$\langle f|S^{\\#}_{(2)}|i\rangle$. Moreover, it is rigorously proved that
this new Yukawa potential is divergence-free. The mass parameter of the
exchanged meson may be set to zero to obtain a type of scalar Boson exchange
between hypothetical Fermions. This provides a limiting case of a new Coulomb
type potential directly from the new singularity free Yukawa potential. A
divergence-free Coulomb potential between two Fermions at two discrete points
is shown to be proportional to the Euler beta function. Within this
relativistic discrete phase space continuous time, a single quanta is shown to
occupy the hyper-tori $S^{1}_{n^{1}}\times S^{1}_{n^{3}}\times S^{1}_{n^{3}}$
where $S^{1}_{n}$ is a circle of radius $\sqrt{2n+1}$.

###### Keywords:

Discrete phase space partial difference-differential equations non-singular
Yukawa potential

###### pacs:

11.10Ef 11.10Qr 11.15Ha 02.30Em 03.65Fd

##  1 Introduction

Quantum mechanics has been exactly represented by the standard partial
differential equations of Schrödinger, the matrix mechanics of Heisenberg, and
the phase space continuous time continuum of Weyl and Wigner Weyl ; Wigner .
In recent years, an exact representation of quantum mechanics has been
introduced in the discrete phase space and continuous time arena DasI ; DasII
; DasIII ; DasIV ; DasV  with the use of partial difference-differential
equations.

The second quantization of relativistic interacting fields in a background
space-time continuum has given rise to the usual $S$-matrix theory Jauch ;
Peskin ; Weinberg  and its well known problem of divergences and subsequent
renormalization attempts. On the other hand, the analogous $S^{\\#}$-matrix in
the arena of discrete phase space and continuous time has been shown to soften
the the degree of divergences DasIII ; DasIV ; DasV ; DasVI . In fact, in
DasVI , it was shown the the $S^{\\#}$-matrix formulation for a second order
electron-electron scattering (or Møller scattering) led, in the low momenta
approximation of two external electrons, to a new Coulomb potential that was
completely divergence-free.

In this paper, we investigate the problem of a new Yukawa style potential
arising from the discrete phase space-continuous time representation of
relativistic quantum field theory and the corresponding $S^{\\#}$-matrix. We
show this new Yukawa style potential is completely divergence-free. In section
2, we summarize briefly the notations used in the present paper. Section 3
defines various partial difference operators DasIII ; DasIV ; DasV ; Jordan
and the corresponding basis of Hermite polynomials Gradshteyn . The following
three sections describe the second quantization of free relativistic scalar
and Fermionic fields and thereafter, interacting fields in the discrete phase
space-continuous time $S^{\\#}$-matrix formulation DasIII ; DasIV ; DasV ;
Jauch ; Peskin ; Weinberg ; DasVI . The specific example of Møller scattering
is investigated.

Section 7 describes the new discrete phase space Yukawa-style potential. The
Green’s function of the corresponding partial difference equation in such a
static scalar field with Yukawa mass parameter $\mu$ is shown to have a zero
discrete phase space distance limit of
$\mu\;\exp(\mu^{2})\;\Gamma(-1/2,\mu^{2})$ where $\Gamma(-1/2,\mu^{2})$ is a
non-singular incomplete gamma function Gradshteyn  for $\mu>0$. This
divergence free quality is the main attraction of our new Yukawa potential.

In section 8, we take the $\mu=0$ limit of our new Yukawa potential to find a
new non-singular Coulomb potential. The Green’s function of the corresponding
partial difference equation is shown to be proportional to the Euler beta
function similar to that found in string theory Green . The coincidence
discrete space limit is shown to produce a non-singular result similar to our
previous work DasVI . We point out that a single field quanta in this
relativistic discrete phase space formulation occupy hyper-tori
$S^{1}_{n^{1}}\times S^{1}_{n^{3}}\times S^{1}_{n^{3}}$ where $S^{1}_{n}$ is a
circle of radius $\sqrt{2n+1}$ as discussed extensively in DasRCI ; DasRCII .
Finally, section 9 summarizes some of the key results of the paper while
illustrating the divergence free nature of our new Yukawa style potential.

##  2 Notations and preliminary definitions

There exists a characteristic length $l>0$ in this theory, which is
conjectured to be the Planck length. We employ fundamental units characterized
by $\hbar=c=l=1$. Thus, all mathematical expressions involving physical
phenomena appear as physically dimensionless numbers. Greek indices take
values from $\\{1,2,3,4\\}$ whereas the Roman indices take values from
$\\{1,2,3\\}$. Einstein’s summation convention is adopted in both cases. We
denote the flat space-time metric by $\eta_{\mu\nu}$ with the corresponding
diagonal matrix $[\eta_{\mu\nu}]:=diag[1,1,1,-1]$. Therefore, we use a
signature of +2 in this paper. An element of our discrete space phase and
continuous time is expressed as
$(\mathbf{n},x^{4})\equiv(n^{1},n^{2},n^{3},t)\in\mathbb{N}^{3}\times\mathbb{R},n^{j}\in\mathbb{N}$
for $j\in{1,2,3}$ and $x^{4}\equiv t\in\mathbb{R}$.

Let a real or complex-valued function $f$ from
$\mathbb{N}^{3}\times\mathbb{R}$ into $\mathbb{R}$ or $\mathbb{C}$ be denoted
as $f(\mathbf{n},t)=f(n^{1},n^{2},n^{3},t)$. We denote the first quantized
wave function $\phi(\mathbf{n},t)$ and the second quantized wave function by
the same symbol. The context should indicate the quantization order.

##  3 Partial Difference and Differential Operators

We denote various partial difference operators and the partial differential
operators as follows DasIII ; DasIV ; DasV ; DasVI

$\begin{array}[]{c}\Delta_{j}f(\mathbf{n},t):=f(...,n^{j}+1,...,t)-f(...,n^{j},...,t)\\\
\\\
\Delta_{j}^{{}^{\prime}}f(\mathbf{n},t):=f(...,n^{j},...,t)-f(...,n^{j}-1,...,t)\\\
\\\
\Delta_{j}^{\\#}f(\mathbf{n},t):=\dfrac{1}{\sqrt{2}}\left[\sqrt{n^{j}+1}f(...,n^{j}+1,...,t)-\sqrt{n^{j}}f(...,n^{j}-1,...,t)\right]\\\
\\\
\Delta_{j}^{o}f(\mathbf{n},t):=\dfrac{1}{\sqrt{2}}\left[\sqrt{n^{j}+1}f(...,n^{j}+1,...,t)+\sqrt{n^{j}}f(...,n^{j}-1,...,t)\right]\\\
\\\ \partial_{t}f(\mathbf{n},t):=\dfrac{\partial}{\partial
t}[f(\mathbf{n},t)]\end{array}$ (1)

We now introduce Hermite polynomials and some pertinent properties by the
following equations DasIII ; DasIV ; DasV ; Gradshteyn ,

$\begin{array}[]{c}H_{n^{j}}(k_{j}):=(-1)^{n^{j}}e^{(k_{j})^{2}}\dfrac{d^{n^{j}}}{(dk_{j})^{n^{j}}}[e^{-(k_{j})^{2}}]\\\
\\\
\dfrac{d^{2}}{(dk_{j})^{2}}[H_{n^{j}}(k_{j})]-2k_{j}\dfrac{d}{dk_{j}}[H_{n^{j}}(k_{j})]+2n^{j}[H_{n^{j}}(k_{j})]=0\\\
\\\ \dfrac{d}{dk_{j}}[H_{n^{j}}(k_{j})]=2n^{j}[H_{n^{j}-1}(k_{j})],\;n^{j}\geq
1\\\ \\\
H_{n^{j}+1}(k_{j})=2k_{j}[H_{n^{j}}(k_{j})]-2n^{j}[H_{n^{j}-1}(k_{j})],\;\;n^{j}\geq
1\end{array}$ (2)

and the scaled Hermite function

$\begin{array}[]{c}\xi_{n^{j}}(k_{j}):=\dfrac{(i)^{n^{j}}e^{-(k_{j})^{2}/2}H_{n^{j}}(k_{j})}{(\pi)^{1/4}2^{(n^{j}/2)}\sqrt{(n^{j})!}}\\\
\\\
\displaystyle\int\limits_{{\mathbb{R}}^{3}}\left\\{\prod_{j=1}^{3}\left[\xi_{n^{j}}(k_{j})\overline{\xi_{\hat{n}^{j}}(k_{j})}\right]\right\\}dk_{1}dk_{2}dk_{3}=\delta_{n^{1}\hat{n}^{1}}\delta_{n^{2}\hat{n}^{2}}\delta_{n^{3}\hat{n}^{3}}=:\delta^{3}_{\mathbf{n}\hat{\mathbf{n}}}\\\
\\\ -i\Delta^{\\#}\xi_{n^{j}}(k_{j})=k_{j}\xi_{n^{j}}(k_{j})\\\ \\\
-i\Delta^{\\#}\overline{\xi_{n^{j}}(k_{j})}=-k_{j}\overline{\xi_{n^{j}}(k_{j})}\end{array}$
(3)

##  4 The second quantization of the free relativistic partial difference-
differential scalar wave field equation

The second quantized Klein-Gordon scalar wave field over the discrete phase
space-continuous time arena will be denoted as
$\phi(\mathbf{n},t)=\phi^{\dagger}(\mathbf{n},t)$. It is Hermitian linear
operator acting on a Hilbert space bundle. The linear operator
$\phi(\mathbf{n},t))$ satisfies the partial difference-differential equations
DasIII ; DasIV ; DasV :

$\delta^{ab}\Delta_{a}^{\\#}\Delta_{b}^{\\#}\phi(\mathbf{n},t)-(\partial_{t})^{2}\phi(\mathbf{n},t)-\mu^{2}\phi(\mathbf{n},t)=0$
(4)

where $\mu>0$ denotes the mass of a neutral scalar Boson. A class of exact
plane wave solutions of the above equation is furnished by

$\begin{array}[]{c}\phi^{(-)}(\mathbf{n},t)=\displaystyle\int\limits_{{\mathbb{R}}^{3}}d^{3}\mathbf{k}\,\,\,[2\omega(\mathbf{k})]^{-1/2}\left\\{a(\mathbf{k})\left[\displaystyle{\prod_{j=1}^{3}}\xi_{n^{j}}(k_{j})\right]e^{-i\omega
t}\right\\}\\\ \\\
\phi^{(+)}(\mathbf{n},t)=\displaystyle\int\limits_{{\mathbb{R}}^{3}}d^{3}\mathbf{k}\,\,\,[2\omega(\mathbf{k})]^{-1/2}\left\\{a^{\dagger}(\mathbf{k})\left[\displaystyle{\prod_{j=1}^{3}}\overline{\xi_{n^{j}}(k_{j})}\right]e^{i\omega
t}\right\\}\\\ \\\
\phi(\mathbf{n},t)=\phi^{(-)}(\mathbf{n},t)+\phi^{(+)}(\mathbf{n},t)=\phi^{\dagger}(\mathbf{n},t)\\\
\\\
\omega\equiv\omega(\mathbf{k})=k^{4}=-k_{4}=+\sqrt{\mathbf{k}\cdot\mathbf{k}+\mu^{2}}>0\end{array}$
(5)

where these improper integrals are supposed to converge uniformly DasVII ;
Whittaker . The linear operator $\phi(\mathbf{n},t)$ representing Bosonic
quanta are assumed to possess the relativistic four-momentum
$(k^{1},k^{2},k^{3},k^{4})=(k_{1},k_{2},k_{3},+\omega)=(\mathbf{k},+\omega)$.

The canonical quantization rules for operators that act on an infinite
dimensional Hilbert space $a_{\mu}(\mathbf{k})$ and
$a^{\dagger}_{\mu}(\mathbf{k})$ are assumed to be the commutators:

$\begin{array}[]{c}[a(\mathbf{k}),a^{\dagger}(\hat{\mathbf{k}})]=-[a^{\dagger}(\hat{\mathbf{k}}),a(\mathbf{k})]=\delta^{3}(\mathbf{k}-\hat{\mathbf{k}})\mathbf{I}(\mathbf{k})\\\
\\\
\;[a(\mathbf{k}),a(\hat{\mathbf{k}})]=[a^{\dagger}(\mathbf{k}),a^{\dagger}((\hat{\mathbf{k}})]=\mathbf{0}(\mathbf{k})\end{array}$
(6)

where $\mathbf{I}(\mathbf{k})$ and $\mathbf{0}(\mathbf{k})$ are the identity
and zero operators respectively and $\delta^{3}(\mathbf{k}-\hat{\mathbf{k}})$
is the Dirac distribution function DasI ; DasII . These commutation relations
imply the following commutators for our second quantized Bosonic field

$\begin{array}[]{c}[\phi^{(+)}(\mathbf{n},t),\phi^{(+)}(\hat{\mathbf{n}},\hat{t})]=[\phi^{(-)}(\mathbf{n},t),\phi^{(-)}(\hat{\mathbf{n}},\hat{t})]=\mathbf{0}\\\
\\\
\;[\phi^{(-)}(\mathbf{n},t),\phi^{(+)}(\hat{\mathbf{n}},\hat{t})]=-i\Delta_{(+)}(\mathbf{n},t;\hat{\mathbf{n}},\hat{t};\mu)\mathbf{I}\\\
\\\
\;[\phi^{(+)}(\mathbf{n},t),\phi^{(-)}(\hat{\mathbf{n}},\hat{t})]=i\Delta_{(-)}(\mathbf{n},t;\hat{\mathbf{n}},\hat{t};\mu)\mathbf{I}\\\
\\\
\;[\phi(\mathbf{n},t),\phi^{\dagger}(\hat{\mathbf{n}},\hat{t})]=-i\Delta(\mathbf{n},t;\hat{\mathbf{n}},\hat{t};\mu)\mathbf{I}\\\
\\\ \end{array}$ (7)

Here, $\Delta_{(\pm)}(\mathbf{n},t;\hat{\mathbf{n}},\hat{t};\mu)$ and
$\Delta(\mathbf{n},t;\hat{\mathbf{n}},\hat{t};\mu)$ are non-singular Green’s
functions for equation (4) as discussed in the Appendix.

##  5 Second quantization of the free relativistic Fermionic anti-Fermionic
partial difference-differential wave equation

Consider the following irreducible representation of $4\times 4$ Dirac
matrices with real and complex entries DasVI ; DasVII

$\begin{array}[]{l}\gamma^{1}_{(4\times 4)}:=\begin{bmatrix}0&0&0&1\\\
0&0&1&0\\\ 0&1&0&0\\\ 1&0&0&0\\\ \end{bmatrix},\;\;\gamma^{2}_{(4\times
4)}:=\begin{bmatrix}0&0&0&-i\\\ 0&0&i&0\\\ 0&-i&0&0\\\ i&0&0&0\\\
\end{bmatrix}\\\ \\\ \gamma^{3}_{(4\times 4)}:=\begin{bmatrix}0&0&1&0\\\
0&0&0&-1\\\ 1&0&0&0\\\ 0&-1&0&0\\\ \end{bmatrix},\;\;\gamma^{4}_{(4\times
4)}:=\begin{bmatrix}-i&0&0&0\\\ 0&-i&0&0\\\ 0&0&i&0\\\ 0&0&0&i\\\
\end{bmatrix}\end{array}$ (8)

and their important properties

$\begin{array}[]{c}\gamma^{a\dagger}=\gamma^{a},\;\gamma^{4\dagger}=-\gamma^{4}\\\
\\\
\gamma^{\mu}\gamma^{\nu}+\gamma^{\nu}\gamma^{\mu}=2\eta^{\mu\nu}I_{(4\times
4)}\end{array}$ (9)

The Dirac bispinor wave field
$\psi(\mathbf{n},t)=[\psi(\mathbf{n},t)]_{(4\times 1)}$ is a $4\times 1$
column vector field with entries of complex numbers. It physically represents
a relativistic Fermion-anti-Fermion wave field defined over discrete phase
space and continuous time. We introduce a corresponding conjugate $1\times 4$
row wave field as follows,

$\tilde{\psi}(\mathbf{n},t):=i\psi^{\dagger}(\mathbf{n},t)\gamma^{4}$ (10)

The discrete phase space continuous time Fermionic-anti-Fermionic wave
equations are furnished by DasVI ; DasVII

$\begin{array}[]{c}\gamma^{a}\Delta^{\\#}_{a}\psi(\mathbf{n},t)+\gamma^{4}\partial_{t}\psi(\mathbf{n},t)+m\psi(\mathbf{n},t)=0_{(4\times
1)}\\\ \\\
\;[\Delta^{\\#}_{a}\tilde{\psi}(\mathbf{n},t)]\gamma^{a}+[\partial_{t}\tilde{\psi}(\mathbf{n},t)]\gamma^{4}-m\tilde{\psi}(\mathbf{n},t)=0_{(1\times
4)}\end{array}$ (11)

Here, the positive constant $m>0$ represents the mass of the Fermionic-anti-
Fermionic particle. These equations are essentially the partial difference-
differential versions of the Dirac equation. We now explore a class of exact
plane wave solutions by the following trial solution

$\psi(\mathbf{n},t)=\zeta(\mathbf{p},p_{4})\left[\prod_{j=1}^{3}\xi_{n^{j}}(p_{j})\right]e^{ip_{4}t}$
(12)

where $\zeta(\mathbf{p},p_{4})$ is a $(4\times 1)$-column vector function of
four-momentum variables $(\mathbf{p},p_{4})$. By substituting this trial
wavefunction into the partial difference-differential Dirac equation above, we
arrive at the following algebraic equations

$\begin{array}[]{c}\eta^{\mu\nu}p_{\mu}p_{\nu}+m^{2}=0,\\\ \\\
p^{4}=-p_{4}=\pm\sqrt{\delta^{ab}p_{a}p_{b}+m^{2}}=\pm\sqrt{||\mathbf{p}||^{2}+m^{2}}\\\
\\\ E\equiv
E(\mathbf{p}):=+\sqrt{\mathbf{p}\cdot\mathbf{p}+m^{2}}=+\sqrt{||\mathbf{p}||^{2}+m^{2}}>0\end{array}$
(13)

and the following four linearly independent solutions

$\begin{array}[]{c}\zeta_{(r)}(\mathbf{p},p^{4})=\zeta_{(r)}(\mathbf{p},E)=:u_{(r)}(\mathbf{p}),\;\;E=E(\mathbf{p})>0\\\
\\\
\zeta_{(r)}(\mathbf{p},p_{4})=\zeta_{(r)}(\mathbf{p},-E)=:v_{(r)}(\mathbf{p}),\;\;-E=-E(\mathbf{p})<0\\\
\\\ r\in\\{1,2\\}\end{array}$ (14)

where $r\in\\{1,2\\}$ physically represents the spin-up and down cases for the
Fermionic-anti-Fermionic quantas. the four explicit $(4\times 1)$ column
vector solutions are listed below DasVI ; DasVII

$\begin{array}[]{c}u_{(1)}(\mathbf{p})=[(m+E)/2E]^{1/2}\begin{bmatrix}1\\\
0\\\ -i(m+E)^{-1}p_{3}\\\ -i(m+E)^{-1}(p_{1}+ip_{2})\\\ \end{bmatrix},\\\ \\\
u_{(2)}(\mathbf{p})=[(m+E)/2E]^{1/2}\begin{bmatrix}0\\\ 1\\\
-i(m+E)^{-1}(p_{1}-ip_{2})\\\ i(m+E)^{-1}p_{3}\\\ \end{bmatrix},\\\ \\\
v_{(1)}(\mathbf{p})=[(m+E)/2E]^{1/2}\begin{bmatrix}i(m+E)^{-1}p_{3}\\\
-i(m+E)^{-1}(p_{1}+ip_{2})\\\ 1\\\ 0\\\ \end{bmatrix},\\\ \\\
v_{(2)}(\mathbf{p})=[(m+E)/2E]^{1/2}\begin{bmatrix}i(m+E)^{-1}(p_{1}-ip_{2})\\\
-i(m+E)^{-1}p_{3}\\\ 0\\\ 1\\\ \end{bmatrix}.\end{array}$ (15)

Here, $u_{(1)}(\mathbf{p})$ and $u_{(2)}(\mathbf{p})$ represent Fermionic
quanta and $v_{(1)}(\mathbf{p})$ and $v_{(2)}(\mathbf{p})$ represent anti-
Fermionic quanta solutions. The above solutions also satisfy the
orthonormality conditions

$\begin{array}[]{c}\tilde{u}_{(r)}(\mathbf{p})\cdot
u_{(s)}(\mathbf{p})=-\tilde{v}_{(r)}(\mathbf{p})\cdot
v_{(s)}(\mathbf{p})=\delta_{(rs)}\\\ \\\ \tilde{u}_{(r)}(\mathbf{p})\cdot
v_{(s)}(\mathbf{p})=\tilde{v}_{(r)}(\mathbf{p})\cdot
u_{(s)}(\mathbf{p})=0\end{array}$ (16)

In later sections, we will have to consider the case of very low values of
three-momentum $\mathbf{p}$. Therefore, it useful to derive approximations of
low external momenta $||\mathbf{p}||$ to the above solutions. The following
expansions to order $O\left(||\mathbf{p}||^{4}\right)$ are

$\begin{array}[]{c}E(\mathbf{p})=m+\left(||\mathbf{p}||^{2}/2m\right)+O\left(||\mathbf{p}||^{4}\right),\\\
\\\
\;[m+E(\mathbf{p})/2E(\mathbf{p})]^{1/2}=[1-(1/2)\left(||\mathbf{p}||/2m\right)^{2}]+O\left(||\mathbf{p}||^{4}\right)\\\
\\\ u_{(1)}=\begin{bmatrix}1\\\ 0\\\ 0\\\
0\end{bmatrix}+\begin{bmatrix}-(1/2)\left(||\mathbf{p}||^{2}/2m\right)\\\ 0\\\
-i(p_{3}/2m)[1-(3/2)\left(||\mathbf{p}||/2m\right)^{2}]\\\
-i\left(\dfrac{p_{1}+ip_{2}}{2m}\right)[1-(3/2)\left(||\mathbf{p}||/2m\right)^{2}]\\\
\end{bmatrix}+O\left(||\mathbf{p}||^{4}\right)\\\ \\\
u_{(2)}=\begin{bmatrix}0\\\ 1\\\ 0\\\ 0\end{bmatrix}+\begin{bmatrix}0\\\
(1/2)\left(||\mathbf{p}||^{2}/2m\right)\\\
-i\left(\dfrac{p_{1}-ip_{2}}{2m}\right)[1-(3/2)\left(||\mathbf{p}||/2m\right)^{2}]\\\
+i(p_{3}/2m)[1-(3/2)\left(||\mathbf{p}||/2m\right)^{2}]\\\
\end{bmatrix}+O\left(||\mathbf{p}||^{4}\right)\end{array}$ (17)

A class of exact plane wave solutions of the Dirac partial difference-
differential equations above is furnished by the following Fourier-Hermite
integrals DasIII ; DasIV ; DasV ,

$\begin{array}[]{c}\psi_{(4\times
1)}^{(-)}(\mathbf{n},t)=\displaystyle\int\limits_{{\mathbb{R}}^{3}}d^{3}\mathbf{p}\,\,\,[m/E(\mathbf{p})]^{1/2}\left\\{\sum_{r=1}^{2}\alpha_{(r)}(\mathbf{p})u_{(r)}(\mathbf{p})\left(\prod_{j=1}^{3}\xi_{n^{j}}(p_{j})\right)e^{-iEt}\right\\}\\\
\\\ \psi_{(4\times
1)}^{(+)}(\mathbf{n},t)=\displaystyle\int\limits_{{\mathbb{R}}^{3}}d^{3}\mathbf{p}\,\,\,[m/E(\mathbf{p})]^{1/2}\left\\{\sum_{r=1}^{2}\beta^{\dagger}_{(r)}(\mathbf{p})v_{(r)}(\mathbf{p})\left(\prod_{j=1}^{3}\overline{\xi_{n^{j}}(p_{j})}\right)e^{iEt}\right\\}\\\
\\\ \psi_{(4\times 1)}(\mathbf{n},t)=\psi_{(4\times
1)}^{(-)}(\mathbf{n},t)+\psi_{(4\times 1)}^{(+)}(\mathbf{n},t)\\\ \\\
\tilde{\psi}_{(1\times
4)}^{(-)}(\mathbf{n},t)=\displaystyle\int\limits_{{\mathbb{R}}^{3}}d^{3}\mathbf{p}\,\,\,[m/E(\mathbf{p})]^{1/2}\left\\{\displaystyle\sum_{r=1}^{2}\alpha^{\dagger}_{(r)}(\mathbf{p})\tilde{u}_{(r)}(\mathbf{p})\left(\prod_{j=1}^{3}\overline{\xi_{n^{j}}(p_{j})}\right)e^{iEt}\right\\}\\\
\\\ \tilde{\psi}_{(1\times
4)}^{(+)}(\mathbf{n},t)=\displaystyle\int\limits_{{\mathbb{R}}^{3}}d^{3}\mathbf{p}\,\,\,[m/E(\mathbf{p})]^{1/2}\left\\{\displaystyle\sum_{r=1}^{2}\beta_{(r)}(\mathbf{p})\tilde{v}_{(r)}(\mathbf{p})\left(\prod_{j=1}^{3}\xi_{n^{j}}(p_{j})\right)e^{-iEt}\right\\}\\\
\\\ \tilde{\psi}_{(1\times 4)}(\mathbf{n},t)=\tilde{\psi}_{(1\times
4)}^{(-)}(\mathbf{n},t)+\tilde{\psi}_{(1\times
4)}^{(+)}(\mathbf{n},t)\end{array}$ (18)

The 4-component column vector $\psi_{(4\times 1)}^{(-)}(\mathbf{n},t)$ and the
4-component row vector $\tilde{\psi}_{(1\times 4)}^{(-)}(\mathbf{n},t)$ are
associated with the Fermionic wave field whereas the vectors $\psi_{(4\times
1)}^{(+)}(\mathbf{n},t)$ and $\tilde{\psi}_{(1\times 4)}^{(+)}(\mathbf{n},t)$
are the anti-Fermionic wave fields.

Now, we shall introduce the canonical or second quantization of the free
Fermionic-anti-Fermionic wave fields. We adopt a two-dimensional pre-Hilbert
space bundle for this purpose. We postulate that the Dirac fields above act as
vector bundles in a two-dimensional pre-Hilbert space bundle. Moreover, the
Fourier-Hermite coefficients
$\alpha_{(r)}(\mathbf{p}),\alpha^{\dagger}_{(r)}(\mathbf{p}),\beta_{(r)}(\mathbf{p}),\beta^{\dagger}_{(s)}(\hat{\mathbf{p}})$
act as linear operators on this pre-Hilbert space bundle. We assume that these
linear operators satisfy anti-commutation rules Jauch ,

$\begin{array}[]{c}\;[A,B]_{+}\equiv AB+BA=[B,A]_{+}\\\ \\\
\;[\alpha_{(r)}(\mathbf{p}),\alpha_{(s)}(\hat{\mathbf{p}})]_{+}=[\beta_{(r)}(\mathbf{p}),\beta_{(s)}(\hat{\mathbf{p}})]_{+}=\mathbf{0}\\\
\\\
\;[\alpha^{\dagger}_{(r)}(\mathbf{p}),\alpha^{\dagger}_{(s)}(\hat{\mathbf{p}})]_{+}=[\beta^{\dagger}_{(r)}(\mathbf{p}),\beta^{\dagger}_{(s)}(\hat{\mathbf{p}})]_{+}=\mathbf{0}\\\
\\\
\;[\alpha_{(r)}(\mathbf{p}),\beta_{(s)}(\hat{\mathbf{p}})]_{+}=[\alpha^{\dagger}_{(r)}(\mathbf{p}),\beta^{\dagger}_{(s)}(\hat{\mathbf{p}})]_{+}=\mathbf{0}\\\
\\\
\;[\alpha_{(r)}(\mathbf{p}),\beta^{\dagger}_{(s)}(\hat{\mathbf{p}})]_{+}=[\alpha^{\dagger}_{(r)}(\mathbf{p}),\beta_{(s)}(\hat{\mathbf{p}})]_{+}=\mathbf{0}\\\
\\\
\;[\alpha_{(r)}(\mathbf{p}),\alpha^{\dagger}_{(s)}(\hat{\mathbf{p}})]_{+}=[\beta_{(r)}(\mathbf{p}),\beta^{\dagger}_{(s)}(\hat{\mathbf{p}})]_{+}=\delta_{(rs)}\delta^{3}(\mathbf{p}-\hat{\mathbf{p}})\mathbf{I}\end{array}$
(19)

There are key physical motivations for choosing pre-Hilbert space bundle
operators
$\alpha_{(r)}(\mathbf{p}),\alpha^{\dagger}_{(s)}(\hat{\mathbf{p}}),\beta_{(r)}(\mathbf{p}),\beta^{\dagger}_{(s)}(\hat{\mathbf{p}})$
satisfying anti-commutation relations in contrast to Hilbert space operators
$a_{\mu}(\mathbf{k}),a^{\dagger}_{\nu}(\hat{\mathbf{k}})$ satisfying
commutation relations and acting on an infinite dimensional Hilbert space
bundle. These differences are due to the fact that Bosonic fields obey Bose-
Einstein statistics, whereas the Fermionic fields obey Fermi-Dirac statistics.

Now, we present various anti-commutation relations between various quantized
Fermionic and anti-Fermionic wave fields,

$\begin{array}[]{c}\;[\psi^{(-)}(\mathbf{n},t),{\psi}^{(+)}(\hat{\mathbf{n}},\hat{t})]_{+}=[\psi^{(-)}(\mathbf{n},t),\tilde{\psi}^{(+)}(\hat{\mathbf{n}},\hat{t})]_{+}=[0]_{(4\times
4)}\\\ \\\
\;[\tilde{\psi}^{(+)}(\mathbf{n},t),{\psi}^{(-)}(\hat{\mathbf{n}},\hat{t})]_{+}=[\tilde{\psi}^{(+)}(\mathbf{n},t),\tilde{\psi}^{(-)}(\hat{\mathbf{n}},\hat{t})]_{+}=[0]_{(4\times
4)}\\\ \\\ \\\
\;[\psi^{(-)}(\mathbf{n},t),\tilde{\psi}^{(-)}(\hat{\mathbf{n}},\hat{t})]_{+}=i[S_{(+)}^{\\#}(\mathbf{n},t;\mathbf{\hat{n}},\hat{t};m)]_{(4\times
4)}\\\ \\\
\;[\tilde{\psi}^{(+)}(\mathbf{n},t),\psi^{(+)}(\hat{\mathbf{n}},\hat{t})]_{+}=i[S_{(-)}^{\\#}(\mathbf{n},t;\mathbf{\hat{n}},\hat{t};m)]_{(4\times
4)}\\\ \\\ \end{array}$ (20)

The various non-singular Green’s functions
$S_{(a)}^{\\#}(\mathbf{n},t;\mathbf{\hat{n}},\hat{t};m)$ will be elaborated
further in the Appendix.

Let us prove the first non-trivial anti-commutation relation above.

$\begin{array}[]{c}\;[\psi^{(-)}(\mathbf{n},t),\tilde{\psi}^{(-)}(\hat{\mathbf{n}},\hat{t})]_{+}=\displaystyle\sum_{r=1}^{2}\displaystyle\sum_{s=1}^{2}\displaystyle\int\limits_{{\mathbb{R}}^{3}}\displaystyle\int\limits_{{\mathbb{R}}^{3}}d^{3}\mathbf{p}d^{3}\mathbf{\hat{p}}\,\,\,\left[\dfrac{m}{\sqrt{E(\mathbf{p})\hat{E}(\mathbf{\hat{p}})}}\right]\\\
\\\
\;[\alpha_{(r)}(\mathbf{p})\alpha^{\dagger}_{(s)}(\mathbf{\hat{p}})]_{+}[u_{(r)}(\mathbf{p})\tilde{u}_{(s)}(\mathbf{\hat{p}})]\displaystyle\left(\prod_{j=1}^{3}\xi_{n^{j}}(p_{j})\overline{\xi_{\hat{n}^{j}}(\hat{p}_{j})}\right)e^{-iEt+i\hat{E}\hat{t}}\\\
\\\ or\\\ \\\
\;[\psi^{(-)}(\mathbf{n},t),\tilde{\psi}^{(-)}(\hat{\mathbf{n}},\hat{t})]_{+}=\displaystyle\sum_{r=1}^{2}\displaystyle\int\limits_{{\mathbb{R}}^{3}}\left[\dfrac{m}{E(\mathbf{p})}\right][u_{(r)}(\mathbf{p})\tilde{u}_{(r)}(\mathbf{p})]\\\
\displaystyle\prod_{j=1}^{3}\xi_{n^{j}}(p_{j})\overline{\xi_{\hat{n}^{j}}(p_{j})}e^{-iE(t-\hat{t})}d^{3}\mathbf{p}\end{array}$
(21)

Now, it can be shown that Jauch

$\left[\dfrac{m}{E(\mathbf{p})}\right]\displaystyle\sum_{r=1}^{2}[u_{(r)}(\mathbf{p})\tilde{u}_{(r)}(\mathbf{p})]=\left[\dfrac{-i\gamma^{j}p_{j}+i\gamma^{4}E+mI}{2E}\right]_{(4\times
4)}$ (22)

and therefore

$\begin{array}[]{c}\;[\psi^{(-)}(\mathbf{n},t),\tilde{\psi}^{(-)}(\hat{\mathbf{n}},\hat{t})]_{+}=\\\
\\\
-\displaystyle\int\limits_{{\mathbb{R}}^{3}}\left[\dfrac{i\gamma^{j}p_{j}-i\gamma^{4}E-mI}{2E}\right]\displaystyle\prod_{j=1}^{3}\xi_{n^{j}}(p_{j})\overline{\xi_{\hat{n}^{j}}(p_{j})}e^{-iE(t-\hat{t})}d^{3}\mathbf{p}\end{array}$
(23)

where in the appendix, it is shown the the Fermionic Green’s function is

$\begin{array}[]{c}[S_{(+)}(\mathbf{n},t;\mathbf{\hat{n}},\hat{t};m)]_{(4\times
4)}=\\\ \\\
i\displaystyle\int\limits_{{\mathbb{R}}^{3}}\left[\dfrac{i\gamma^{j}p_{j}-i\gamma^{4}E-mI}{2E}\right]\displaystyle\prod_{j=1}^{3}\xi_{n^{j}}(p_{j})\overline{\xi_{\hat{n}^{j}}(p_{j})}e^{-iE(t-\hat{t})}d^{3}\mathbf{p}\end{array}$
(24)

Therefore, our result follows.

##  6 Interactions of relativistic fields, the $S^{\\#}$-matrix in discrete
phase space-continuous time, and Møller Scattering

The relativistic Lagrangian of the second quantized interacting Fermionic and
Bosonic fields is assumed to given by DasIII ; DasIV ; DasV ; Weinberg

$L_{int}(\mathbf{n},t):=-igN[\tilde{\psi}(\mathbf{n},t)\psi(\mathbf{n},t)\phi(\mathbf{n},t)].$
(25)

Here, $g$ represents the strong coupling constant between the Fermionic-anti-
Fermionic field and a Bosonic field. Furthermore, $N[\cdots]$ stands for
normal ordering. The scattering matrix in the background of discrete phase
space and continuous time, denoted by $S^{\\#}$-matrix, is defined by the
operator-valued infinite series DasIII ; DasIV ; DasV ; Jauch ; Peskin ;
Weinberg

$\begin{array}[]{c}S^{\\#}=I+\sum\limits_{j=1}^{\infty}\dfrac{(g)^{j}}{j!}\sum\limits_{\mathbf{n}^{1}=(0)}\cdots\sum\limits_{\mathbf{n}^{j}=(0)}\displaystyle\int\limits_{\mathbb{R}}dt^{1}\\\
\\\
\cdots\displaystyle\int\limits_{\mathbb{R}}dt^{j}T\left\\{N[\tilde{\psi}(\mathbf{n}^{1},t_{1})\psi(\mathbf{n}^{1},t^{1})\phi(\mathbf{n}^{1},t^{1})]\right.\\\
\\\
\cdots\left.N[\tilde{\psi}(\mathbf{n}^{j},t^{j})\psi(\mathbf{n}^{1},t^{j})\phi(\mathbf{n}^{j},t^{j})]\right\\}=:I+\displaystyle\sum_{j=1}^{\infty}S^{\\#}_{(j)}\end{array}$
(26)

Here, $T$ denotes Wick’s time ordering operation. We distinguish the
scattering matrix by the notation $S^{\\#}$-matrix from the usual notation of
S-matrix in continuous space-time because the physics in the discrete phase
space and continuous time is different from the physics in the space-time
continuum.

Now, let us consider a particular physical process characterized by the
initial Hilbert space state vector $|i\rangle$ and the final Hilbert space
state vector $|f\rangle$. The $S^{\\#}$-matrix elements for such a process is
provided by

$\langle f|S^{\\#}|i\rangle=\langle
f|i\rangle+\sum\limits_{j=1}^{\infty}\langle f|S^{\\#}_{(j)}|i\rangle$ (27)

Furthermore, we need to define the following vertex distribution functions

$\begin{array}[]{c}\delta^{\\#}_{(3)}(\mathbf{p},\mathbf{q},\mathbf{k}):=\displaystyle\sum_{n^{1}=0}^{\infty}\sum_{n^{2}=0}^{\infty}\sum_{n^{3}=0}^{\infty}\left[\prod_{j=1}^{3}\xi_{n^{j}}(p_{j})\xi_{n^{j}}(q_{j})\xi_{n^{j}}(k_{j})\right]\\\
\\\
\delta^{\\#}_{(3)}(\mathbf{p},-\mathbf{q},-\mathbf{k}):=\displaystyle\sum_{n^{1}=0}^{\infty}\sum_{n^{2}=0}^{\infty}\sum_{n^{3}=0}^{\infty}\left[\prod_{j=1}^{3}\xi_{n^{j}}(p_{j})\overline{\xi_{n^{j}}(q_{j})}\overline{\xi_{n^{j}}(k_{j})}\right]\end{array}$
(28)

The function $\delta^{\\#}_{(3)}(\mathbf{p},-\mathbf{q},-\mathbf{k})$ above is
different from the standard delta function factor
$(2\pi)^{3}\delta^{3}(\mathbf{p},-\mathbf{q},-\mathbf{k})$ found in the usual
theory in space-time continuum. Also, the problem of reducing equation (28)
into much simpler functions is an unsolved problem.

We provide Feynman graph rules to evaluate succinctly each term of the
$S^{\\#}$-matrix series in Table 1.

Table 1: Feynman Graphs in four-momentum space

Note that the relativistic Feynman prescriptions are identical between our
$S^{\\#}$-matrix and space-time continuum $S$ matrix except on the vertices.
At a vertex in our formalism, one has the relativistic term

$g\delta^{\\#}_{(3)}(\mathbf{p},-\mathbf{q},-\mathbf{k})\delta(p^{4}-q^{4}-k^{4}).$
(29)

As an illustration of computing explicitly an element of$\langle
f|S^{\\#}_{(j)}-I|i\rangle$ in (27), using the corresponding component of the
Feynman rules in Table 1, we consider the case of Møller scattering of two
Fermions by exchange of one Boson as depicted in figure 1.

Figure 1: The Møller scattering of two Fermions by exchange of one Boson

The element $\langle
p^{{}^{\prime}}_{2},p^{{}^{\prime}}_{1}|S^{\\#}_{(2)}|p_{1},p_{2}\rangle$
corresponding to figure 1 will be evaluated using the very low momentum
approximation for the external Fermion lines as discussed above in (17).
Following this, we obtain

$\begin{array}[]{c}\langle
p^{{}^{\prime}}_{2},p^{{}^{\prime}}_{1}|S^{\\#}_{(2)}|p_{1},p_{2}\rangle\approx\left(\dfrac{g^{2}}{4\pi}\right)\left[\dfrac{m^{2}}{\sqrt{E^{{}^{\prime}}_{1}E^{{}^{\prime}}_{2}E_{1}E_{2}}}\right]\delta(E^{{}^{\prime}}_{2}+E^{{}^{\prime}}_{1}-E_{2}-E_{1})\\\
\\\
\;[\delta_{(r^{\prime}_{1},r_{1})}\delta_{(r^{\prime}_{2},r_{2})}]\displaystyle\int\limits_{{\mathbb{R}}^{3}}\left\\{\delta^{\\#}_{(3)}(\mathbf{p_{1}},-\mathbf{p^{{}^{\prime}}_{1}},-\mathbf{k})[\mathbf{k}\cdot\mathbf{k}+(\mu)^{2}]^{-1}\delta^{\\#}_{(3)}(\mathbf{p_{2}},-\mathbf{p^{{}^{\prime}}_{2}},+\mathbf{k})\right\\}d^{3}\mathbf{k}\\\
\\\
=\left(\dfrac{g^{2}}{4\pi}\right)\left[\dfrac{m^{2}}{\sqrt{E^{{}^{\prime}}_{1}E^{{}^{\prime}}_{2}E_{1}E_{2}}}\right]\delta(E^{{}^{\prime}}_{2}+E^{{}^{\prime}}_{1}-E_{2}-E_{1})[\delta_{(r^{\prime}_{1},r_{1})}\delta_{(r^{\prime}_{2},r_{2})}]\\\
\\\
\displaystyle\int\limits_{{\mathbb{R}}^{3}}\biggl{\\{}[\mathbf{k}\cdot\mathbf{k}+\mu^{2}]^{-1}\left[\displaystyle\sum_{\mathbf{n}=0}^{\infty(3)}\prod_{a=1}^{3}\xi_{n^{a}}(p_{1a})\overline{\xi_{n^{a}}(p^{\prime}_{1a})}\overline{\xi_{n^{a}}(k_{a})}\right]\\\
\\\
\left[\displaystyle\sum_{\mathbf{\hat{n}}=0}^{\infty(3)}\prod_{b=1}^{3}\xi_{\hat{n}^{b}}(p_{2b})\overline{\xi_{\hat{n}^{b}}(p^{\prime}_{2b})}\xi_{\hat{n}^{b}}(k_{b})\right]\biggr{\\}}d^{3}\mathbf{k}\end{array}$
(30)

Inside the above equation, one may find the double Hermite transform for the
Green’s function defined as

$G^{\\#}(\mathbf{n},\mathbf{\hat{n}};\mu):=\displaystyle\int\limits_{{\mathbb{R}}^{3}}[\mathbf{k}\cdot\mathbf{k}+\mu^{2}]^{-1}\left[\displaystyle\prod_{j=1}^{3}\overline{\xi_{n^{j}}(k_{j})}\xi_{\hat{n}^{j}}(k_{j})\right]d^{3}\mathbf{k}$
(31)

To compare and contrast the above result (31) with the usual second order
S-matrix elements for Møller scattering over space-time continuum, we furnish
from Jauch ; DasVI

$\begin{array}[]{c}\langle
p^{{}^{\prime}}_{2},p^{{}^{\prime}}_{1}|S_{(2)}|p_{1},p_{2}\rangle\approx\left(\dfrac{g^{2}}{4\pi}\right)\left[\dfrac{m^{2}}{\sqrt{E^{{}^{\prime}}_{1}E^{{}^{\prime}}_{2}E_{1}E_{2}}}\right]\delta(E^{{}^{\prime}}_{2}+E^{{}^{\prime}}_{1}-E_{2}-E_{1})\\\
\\\
\;[\delta_{(r^{\prime}_{1},r_{1})}\delta_{(r^{\prime}_{2},r_{2})}]\displaystyle\int\limits_{{\mathbb{R}}^{3}}\left\\{\delta^{3}(\mathbf{p_{1}},-\mathbf{p^{{}^{\prime}}_{1}},-\mathbf{k})[\mathbf{k}\cdot\mathbf{k}+(\mu)^{2}]^{-1}\delta^{3}(\mathbf{p_{2}},-\mathbf{p^{{}^{\prime}}_{2}},+\mathbf{k})\right\\}d^{3}\mathbf{k}\\\
\\\
=\left(\dfrac{g^{2}}{4\pi}\right)\left[\dfrac{m^{2}}{\sqrt{E^{{}^{\prime}}_{1}E^{{}^{\prime}}_{2}E_{1}E_{2}}}\right]\delta(E^{{}^{\prime}}_{2}+E^{{}^{\prime}}_{1}-E_{2}-E_{1})[\delta_{(r^{\prime}_{1},r_{1})}\delta_{(r^{\prime}_{2},r_{2})}]\\\
\\\
\displaystyle\int\limits_{{\mathbb{R}}^{3}}\biggl{\\{}[\mathbf{k}\cdot\mathbf{k}+\mu^{2}]^{-1}\left[\dfrac{1}{(2\pi)^{3/2}}\displaystyle\int\limits_{{\mathbb{R}}^{3}}\exp\\{-i(\mathbf{p_{1}}-\mathbf{p^{{}^{\prime}}_{1}}-\mathbf{k})\cdot\mathbf{x_{1}}\\}d^{3}\mathbf{x_{1}}\right]\\\
\\\
\left[\dfrac{1}{(2\pi)^{3/2}}\displaystyle\int\limits_{{\mathbb{R}}^{3}}\exp\\{-i(\mathbf{p_{2}}-\mathbf{p^{{}^{\prime}}_{2}}-\mathbf{k})\cdot\mathbf{x_{2}}\\}d^{3}\mathbf{x_{2}}\right]\biggr{\\}}d^{3}\mathbf{k}\end{array}$
(32)

The above equation may be regarded as the double Fourier transform for the
Green’s function

$G(\mathbf{x_{1}},\mathbf{x_{2}};\mu)=G(\mathbf{x_{1}}-\mathbf{x_{2}};\mu):=\displaystyle\dfrac{1}{(2\pi)^{3}}\int\limits_{{\mathbb{R}}^{3}}[\mathbf{k}\cdot\mathbf{k}+\mu^{2}]^{-1}\exp\\{-i(\mathbf{x_{1}}-\mathbf{x_{2}})\cdot\mathbf{k}\\}d^{3}\mathbf{k}$
(33)

with the limit

$\lim_{\mathbf{x_{1}}\rightarrow\mathbf{x_{2}}}|G(\mathbf{x_{1}}-\mathbf{x_{2}};\mu)|\rightarrow\infty$
(34)

This equation implies the following static partial differential equation

$\begin{array}[]{c}\delta^{ab}\partial_{(1)^{a}}\partial_{(1)^{b}}G(\mathbf{x_{1}}-\mathbf{x_{2}};\mu)-\mu^{2}G(\mathbf{x_{1}}-\mathbf{x_{2}};\mu)\\\
\\\
=\delta^{ab}\partial_{(2)^{a}}\partial_{(2)^{b}}G(\mathbf{x_{1}}-\mathbf{x_{2}};\mu)-\mu^{2}G(\mathbf{x_{1}}-\mathbf{x_{2}};\mu)=-\delta^{(3)}(\mathbf{x_{1}}-\mathbf{x_{2}})\end{array}$
(35)

Now, the potential between two Fermions exchanging one Boson is furnished by,

$V(\mathbf{x_{1}}-\mathbf{x_{2}})=g^{2}G(\mathbf{x_{1}}-\mathbf{x_{2}};\mu^{2})=\dfrac{g^{2}}{(2\pi)^{2}}\int\limits_{{\mathbb{R}}^{3}}[\mathbf{k}\cdot\mathbf{k}+\mu^{2}]^{-1}\exp\\{-i(\mathbf{x_{1}}-\mathbf{x_{2}})\cdot\mathbf{k}\\}d^{3}\mathbf{k}$
(36)

To obtain a closed form solution, we transform to spherical polar coordinates

$\begin{array}[]{c}\mathbf{x}=\mathbf{x_{1}}-\mathbf{x_{2}},\;\;\;r=||\mathbf{x}||,\;\;\;k=||\mathbf{k}||\\\
\\\ \mathbf{x}=(r\sin\theta\cos\phi,r\sin\theta\sin\phi,r\cos\theta)\\\ \\\
\mathbf{k}=(k\sin\hat{\theta}\cos\hat{\phi},k\sin\hat{\theta}\sin\hat{\phi},k\cos\hat{\theta})\\\
\\\
\cos\gamma=\sin\theta\sin\hat{\theta}\cos\phi\cos\hat{\phi}+\sin\theta\sin\hat{\theta}\sin\phi\sin\hat{\phi}+\cos\theta\cos\hat{\theta},\;\;\;0\leq\gamma\leq\pi\\\
\\\ \mathbf{x}\cdot\mathbf{k}=rk\cos\gamma=:rky\end{array}$ (37)

thereby obtaining

$V(r)=\dfrac{g^{2}}{(2\pi)^{2}}\int\limits_{0}^{\infty}\int\limits_{0}^{\pi}\int\limits_{-\pi}^{\pi}[\mathbf{k}\cdot\mathbf{k}+\mu^{2}]^{-1}e^{irk\cos\gamma}k^{2}\sin\gamma
dkd\gamma d\phi$ (38)

or

$\begin{array}[]{c}V(r)=\dfrac{g^{2}}{(2\pi)^{2}}\displaystyle\int\limits_{0}^{\infty}\int\limits\limits_{-1}^{1}[\mathbf{k}\cdot\mathbf{k}+\mu^{2}]^{-1}e^{irky}k^{2}dkdy\\\
\\\
=-i\dfrac{g^{2}}{(2\pi)^{2}}\left(\dfrac{1}{r}\right)\displaystyle\int\limits_{0}^{\infty}[\mathbf{k}\cdot\mathbf{k}+\mu^{2}]^{-1}[e^{irk}-e^{-irk}]kdk\\\
\\\
=-i\dfrac{g^{2}}{(2\pi)^{2}}\left(\dfrac{1}{r}\right)\displaystyle\int\limits_{0}^{\infty}[\mathbf{k}\cdot\mathbf{k}+\mu^{2}]^{-1}[2i\sin(rk)]kdk\\\
\\\ \end{array}$ (39)

Note that $[\mathbf{k}\cdot\mathbf{k}+\mu^{2}]^{-1}[2i\sin(rk)]k$ is an even
function whereas $[\mathbf{k}\cdot\mathbf{k}+\mu^{2}]^{-1}[2i\cos(rk)]k$ is an
odd function of $k$. Therefore we may write the above as

$V(r)=-i\dfrac{g^{2}}{(2\pi)^{2}}\left(\dfrac{1}{r}\right)\displaystyle\int\limits_{-\infty}^{\infty}k[\mathbf{k}\cdot\mathbf{k}+\mu^{2}]^{-1}e^{irk}dk$
(40)

The integrand has two simple poles at $k\pm i\mu$. Using a standard counter-
clockwise contour $C$ around the pole $k=i\mu$ in the upper half plane, the
integral above using the Cauchy residue theorem is

$\displaystyle\int\limits_{-\infty}^{\infty}k[\mathbf{k}\cdot\mathbf{k}+\mu^{2}]^{-1}e^{irk}dk=\oint\limits_{C}k[\mathbf{k}\cdot\mathbf{k}+\mu^{2}]^{-1}e^{irk}dk=-i\pi
e^{-\mu r}$ (41)

Therefore, our potential is

$V(r)=-\left(\dfrac{g^{2}}{4\pi}\right)\dfrac{e^{-\mu r}}{r}$ (42)

which is the well know Yukawa potential for the strong interaction between two
protons exchanging a neutral meson. It is regular for $r>0$ but has a
singularity at $r=0$. in Cartesian coordinates, the Yukawa potential for two
interacting Fermions is

$V(\mathbf{x_{1}},\mathbf{x_{2}};\mu)=-\left(\dfrac{g^{2}}{4\pi}\right)\dfrac{\exp(-\mu||\mathbf{x_{1}}-\mathbf{x_{2}}||)}{||\mathbf{x_{1}}-\mathbf{x_{2}}||}$
(43)

In the case of one Fermion situated at the origin $\mathbf{x_{2}}=(0,0,0)$ and
the other situated at $\mathbf{x_{1}}=\mathbf{x}$. the corresponding Yukawa
potential is given by,

$\begin{array}[]{c}V(\mathbf{x_{1}},\mathbf{x_{2}};\mu)=-\left(\dfrac{g^{2}}{4\pi}\right)\dfrac{\exp(-\mu||\mathbf{x}||)}{||\mathbf{x}||}\\\
\\\
=-\left(\dfrac{g^{2}}{4\pi}\right)\dfrac{\exp(-\mu\sqrt{(x^{1})^{2}+(x^{3})^{2}+(x^{3})^{2}})}{\sqrt{(x^{1})^{2}+(x^{3})^{2}+(x^{3})^{2}}}\end{array}$
(44)

##  7 Discrete phase space and a new non-singular Yukawa potential

The static approximation to the scalar field equation (4) in the background of
discrete phase space is furnished by the following partial difference equation

$\delta^{ab}\Delta_{a}^{\\#}\Delta_{b}^{\\#}\phi(\mathbf{n})-\mu^{2}\phi(\mathbf{n})=0$
(45)

with a Green’s function given by

$\begin{array}[]{c}G^{\\#}(\mathbf{n},\mathbf{\hat{n}};\mu)=G^{\\#}(n^{1},n^{2},n^{3},\hat{n}^{1},\hat{n}^{2},\hat{n}^{3};\mu)\\\
\\\
:=\displaystyle\int\limits_{{\mathbb{R}}^{3}}[\mathbf{k}\cdot\mathbf{k}+(\mu)^{2}]^{-1}\left[\displaystyle\prod_{j=1}^{3}\xi_{n^{j}}(k_{j})\overline{\xi_{\hat{n}^{j}}(k_{j})}\right]d^{3}\mathbf{k}\end{array}$
(46)

By direct substitution of this Green’s function into the partial difference
equation (45), one finds, using (2) and (3)

$\begin{array}[]{c}\delta^{ab}\Delta_{a}^{\\#}\Delta_{b}^{\\#}G^{\\#}(\mathbf{n},\mathbf{\hat{n}};\mu)-\mu^{2}G^{\\#}(\mathbf{n},\mathbf{\hat{n}};\mu)\\\
\\\
=\delta^{ab}\hat{\Delta}_{a}^{\\#}\hat{\Delta}_{b}^{\\#}G^{\\#}(\mathbf{n},\mathbf{\hat{n}};\mu)-\mu^{2}G^{\\#}(\mathbf{n},\mathbf{\hat{n}};\mu)\\\
\\\
=-\delta_{(n^{1}\hat{n}^{1})}\delta_{(n^{2}\hat{n}^{2})}\delta_{(n^{3}\hat{n}^{3})}=:-\delta^{3}_{(\mathbf{n}\mathbf{\hat{n}})}\end{array}$
(47)

Furthermore, from using (2) and (3), one has from (46)

$\begin{array}[]{c}G^{\\#}(\mathbf{n},\mathbf{0};\mu)=\left[\dfrac{i^{n^{1}+n^{2}+n^{3}}}{\pi^{3/2}2^{(n^{1}+n^{2}+n^{3})/2}\sqrt{n^{1}!n^{2}!n^{3}!}}\right]\times\\\
\displaystyle\int\limits_{{\mathbb{R}}^{3}}e^{-\mathbf{k}\cdot\mathbf{k}}[\mathbf{k}\cdot\mathbf{k}+\mu^{2}]^{-1}\left[\displaystyle\prod_{j=1}^{3}H_{n^{j}}(k_{j})\right]d^{3}\mathbf{k}\end{array}$
(48)

and

$G^{\\#}(n^{1},0,0,0,0,0;\mu)=\dfrac{i^{n^{1}}}{\pi^{3/2}2^{\frac{n^{1}}{2}}\sqrt{n^{1}!}}\displaystyle\int\limits_{{\mathbb{R}}^{3}}e^{-\mathbf{k}\cdot\mathbf{k}}[\mathbf{k}\cdot\mathbf{k}+\mu^{2}]^{-1}H_{n^{1}}(k_{1})d^{3}\mathbf{k}$
(49)

Now, we introduce spherical polar coordinates in three-momentum space

$\mathbf{k}=(k\cos\theta,k\sin\theta\cos\phi,k\sin\theta\sin\phi),\;\;y=\cos\theta$
(50)

to transform (49) into

$\begin{array}[]{c}G^{\\#}(n^{1},0,0,0,0,0;\mu)=\dfrac{i^{n^{1}}}{\pi^{3/2}2^{\frac{n^{1}}{2}}\sqrt{n^{1}!}}\times\\\
\\\
\displaystyle\int\limits_{0}^{\infty}\int\limits_{0}^{\pi}\int\limits_{-\pi}^{\pi}e^{-k^{2}}\;[k^{2}+\mu^{2}]^{-1}H_{n^{1}}(k\cos\theta)k^{2}\sin\theta
dkd\theta d\phi\\\ \\\
=\dfrac{i^{n^{1}}}{\sqrt{\pi}2^{\frac{n^{1}}{2}-1}\sqrt{n^{1}!}}\displaystyle\int\limits_{0}^{\infty}\int\limits_{-1}^{1}e^{-k^{2}}\;[k^{2}+\mu^{2}]^{-1}H_{n^{1}}(ky)k^{2}dkdy\end{array}$
(51)

Letting

$x=k^{2}>0,\;\;k=+\sqrt{x}>0$ (52)

(51) becomes

$G^{\\#}(n^{1},0,0,0,0,0;\mu)=\dfrac{i^{n^{1}}}{\sqrt{\pi}2^{\frac{n^{1}}{2}}\sqrt{n^{1}!}}\displaystyle\int\limits_{0}^{\infty}\int\limits_{-1}^{1}e^{-x}\;[x+\mu^{2}]^{-1}H_{n^{1}}(\sqrt{x}y)\sqrt{x}dxdy$
(53)

The coincidence limit of $n^{1}\rightarrow 0$ of (53) yields Gradshteyn

$G^{\\#}(0,0,0,0,0,0;\mu)=\mu\;e^{\mu^{2}}\Gamma\left(-\dfrac{1}{2},\mu^{2}\right)$
(54)

Here, the incomplete gamma function is defined by Gradshteyn

$\begin{array}[]{c}\Gamma\left(-\dfrac{1}{2},\mu^{2}\right):=\displaystyle\int\limits_{\mu^{2}}^{\infty}w^{-3/2}e^{-w}dw\\\
\\\
\dfrac{\partial}{\partial\mu^{2}}\Gamma\left(-\dfrac{1}{2},\mu^{2}\right)=-\mu^{3}e^{-\mu^{2}}\end{array}$
(55)

By the above three equations, it has been clearly demonstrated that the
function $G^{\\#}(0,0,0,0,0,0;\mu)$ is divergence-free.

##  8 Discrete phase space, a new non-singular Coulomb potential, and Beta
functions

The Coulomb potential may be obtained by setting $\mu=0$ in the previous
section. The Green’s function from (46) becomes

$\begin{array}[]{c}G^{\\#}(\mathbf{n},\mathbf{\hat{n}};0)=G^{\\#}(n^{1},n^{2},n^{3},\hat{n}^{1},\hat{n}^{2},\hat{n}^{3};0)\\\
\\\
:=\displaystyle\int\limits_{{\mathbb{R}}^{3}}[\mathbf{k}\cdot\mathbf{k}]^{-1}\left[\displaystyle\prod_{j=1}^{3}\xi_{n^{j}}(k_{j})\overline{\xi_{\hat{n}^{j}}(k_{j})}\right]d^{3}\mathbf{k}\end{array}$
(56)

or using Hermite polynomials

$\begin{array}[]{c}G^{\\#}(\mathbf{n},\mathbf{0};0)=\left[\dfrac{i^{n^{1}+n^{2}+n^{3}}}{\pi^{3/2}2^{(n^{1}+n^{2}+n^{3})/2}\sqrt{n^{1}!n^{2}!n^{3}!}}\right]\times\\\
\displaystyle\int\limits_{{\mathbb{R}}^{3}}e^{-\mathbf{k}\cdot\mathbf{k}}[\mathbf{k}\cdot\mathbf{k}]^{-1}\left[\displaystyle\prod_{j=1}^{3}H_{n^{j}}(k_{j})\right]d^{3}\mathbf{k}\end{array}$
(57)

and

$G^{\\#}(n^{1},0,0,0,0,0;0)=\dfrac{i^{n^{1}}}{\pi^{3/2}2^{\frac{n^{1}}{2}}\sqrt{n^{1}!}}\displaystyle\int\limits_{{\mathbb{R}}^{3}}e^{-\mathbf{k}\cdot\mathbf{k}}[\mathbf{k}\cdot\mathbf{k}]^{-1}H_{n^{1}}(k_{1})d^{3}\mathbf{k}$
(58)

Once again, using the spherical coordinates of (50), we deduce that

$G^{\\#}(n^{1},0,0,0,0,0;0)=\dfrac{i^{n^{1}}}{\sqrt{\pi}2^{(\frac{n^{1}}{2}-1)}\sqrt{n^{1}!}}\displaystyle\int\limits_{0}^{\infty}\int\limits_{-1}^{1}e^{-k^{2}}H_{n^{1}}(ky)dkdy$
(59)

and

$G^{\\#}(2n^{1},0,0,0,0,0;0)=\dfrac{(-1)^{n^{1}}}{\sqrt{\pi}2^{(n^{1}-1)}\sqrt{(2n^{1})!}}\displaystyle\int\limits_{0}^{\infty}\int\limits_{-1}^{1}e^{-k^{2}}H_{2n^{1}}(ky)dkdy$
(60)

and

$\begin{array}[]{c}G^{\\#}(2n^{1}+1,0,0,0,0,0;0)=\dfrac{i(-1)^{n^{1}}}{\sqrt{\pi}2^{(n^{1}-1/2)}\sqrt{(2n^{1}+1)!}}\times\\\
\displaystyle\int\limits_{0}^{\infty}\int\limits_{-1}^{1}e^{-k^{2}}H_{2n^{1}+1}(ky)dkdy\end{array}$
(61)

As $H_{2n^{1}}(ky)$ and $H_{2n^{1}+1}(ky)$ are even and odd functions
respectively of the variable $ky$, it is clear that

$\begin{array}[]{c}\displaystyle\int\limits_{0}^{\infty}\int\limits_{-1}^{1}e^{-k^{2}}H_{2n^{1}}(ky)dkdy=\dfrac{1}{2}\displaystyle\int\limits_{-\infty}^{\infty}\int\limits_{-1}^{1}e^{-k^{2}}H_{2n^{1}}(ky)dkdy\\\
\\\ and\\\ \\\
\displaystyle\int\limits_{0}^{\infty}\int\limits_{-1}^{1}e^{-k^{2}}H_{2n^{1}+1}(ky)(ky)dkdy=0\end{array}$
(62)

Furthermore, from Gradshteyn , using the relation

$\displaystyle\int\limits_{-\infty}^{\infty}e^{-k^{2}}H_{2n^{1}}(ky)dk=(-1)^{n^{1}}\sqrt{\pi}\dfrac{(2n^{1})!}{(n^{1}!)}(1-y^{2})^{n^{1}}\\\
\\\ $ (63)

our Green’s functions reduce to

$\begin{array}[]{c}G^{\\#}(2n^{1},0,0,0,0,0;0)=\left[\dfrac{\sqrt{(2n^{1})!}}{2^{n^{1}}(n^{1}!)}\right]\displaystyle\int\limits_{-1}^{1}(1-y^{2})^{n^{1}}dy\\\
\\\ and\\\ \\\ G^{\\#}(2n^{1}+1,0,0,0,0,0;0)=0\end{array}$ (64)

The integral above is the Euler beta function Gradshteyn

$\begin{array}[]{c}\displaystyle\int\limits_{-1}^{1}(1-y^{2})^{n^{1}}dy=2^{(2n^{1}+1)}B(n^{1}+1,n^{1}+1)\\\
\\\
=2^{(2n^{1}+1)}\dfrac{[\Gamma(n^{1}+1)]^{2}}{\Gamma(2n^{1}+2)}=2^{(2n^{1}+1)}\dfrac{[(n^{1})!]^{2}}{[(2n^{1}+1)!]}\end{array}$
(65)

leaving us with

$G^{\\#}(2n^{1},0,0,0,0,0;0)=\dfrac{2^{(n^{1}+1)}\;(n^{1})!\;\sqrt{(2n^{1})!}}{(2n^{1}+1)!}=\dfrac{2^{(n^{1}+1)}\;(n^{1})!}{(2n^{1}+1)\sqrt{(2n^{1})!}}$
(66)

as was initially derived in DasVI . In the coincidence limit $n^{1}\rightarrow
0_{+}$, we have

$G^{\\#}(0,0,0,0,0,0;0)=+2$ (67)

indicating that the discrete phase space Coulomb potential has no singularity
as was derived in a different manner in DasVI .

##  9 A summary of various Discrete phase space non-singular potentials and
Green’s functions

Here, we summarize the various discrete phase space and continuous time
Green’s functions and potentials in a special way. From (46), we have

$G^{\\#}(\mathbf{n},\mathbf{\hat{n}};\mu)=\displaystyle\int\limits_{{\mathbb{R}}^{3}}[\mathbf{k}\cdot\mathbf{k}+(\mu)^{2}]^{-1}\left[\displaystyle\prod_{j=1}^{3}\xi_{n^{j}}(k_{j})\overline{\xi_{\hat{n}^{j}}(k_{j})}\right]d^{3}\mathbf{k}$
(68)

and it’s various special cases

$\begin{array}[]{l}(i)\;G^{\\#}(\mathbf{n},\mathbf{0};\mu)=\left[\dfrac{i^{n^{1}+n^{2}+n^{3}}}{\pi^{3/2}2^{(n^{1}+n^{2}+n^{3})/2}\sqrt{n^{1}!n^{2}!n^{3}!}}\right]\times\\\
\displaystyle\int\limits_{{\mathbb{R}}^{3}}e^{-\mathbf{k}\cdot\mathbf{k}}[\mathbf{k}\cdot\mathbf{k}+\mu^{2}]^{-1}\left[\displaystyle\prod_{j=1}^{3}H_{n^{j}}(k_{j})\right]d^{3}\mathbf{k}\\\
\\\
(ii)\;G^{\\#}(n^{1},0,0,0,0,0;\mu)=\dfrac{i^{n^{1}}}{\pi^{3/2}2^{\frac{n^{1}}{2}}\sqrt{n^{1}!}}\displaystyle\int\limits_{{\mathbb{R}}^{3}}e^{-\mathbf{k}\cdot\mathbf{k}}[\mathbf{k}\cdot\mathbf{k}+\mu^{2}]^{-1}H_{n^{1}}(k_{1})d^{3}\mathbf{k}\\\
\\\
(iii)\;G^{\\#}(2n^{1},0,0,0,0,0;\mu)=\dfrac{(-1)^{n^{1}}}{\pi^{3/2}2^{n^{1}}\sqrt{(2n^{1})!}}\displaystyle\int\limits_{{\mathbb{R}}^{3}}e^{-\mathbf{k}\cdot\mathbf{k}}[\mathbf{k}\cdot\mathbf{k}+\mu^{2}]^{-1}H_{2n^{1}}(k_{1})d^{3}\mathbf{k}\\\
\\\ (iv)\;G^{\\#}(2n^{1}+1,0,0,0,0,0;\mu)=0\\\ \\\
(v)\;G^{\\#}(\mathbf{n},\mathbf{0};0)=\left[\dfrac{i^{n^{1}+n^{2}+n^{3}}}{\pi^{3/2}2^{(n^{1}+n^{2}+n^{3})/2}\sqrt{n^{1}!n^{2}!n^{3}!}}\right]\times\\\
\displaystyle\int\limits_{{\mathbb{R}}^{3}}e^{-\mathbf{k}\cdot\mathbf{k}}[\mathbf{k}\cdot\mathbf{k}]^{-1}\left[\displaystyle\prod_{j=1}^{3}H_{n^{j}}(k_{j})\right]d^{3}\mathbf{k}\\\
\\\
(vi)\;G^{\\#}(n^{1},0,0,0,0,0;0)=\dfrac{i^{n^{1}}}{\pi^{3/2}2^{\frac{n^{1}}{2}}\sqrt{n^{1}!}}\displaystyle\int\limits_{{\mathbb{R}}^{3}}e^{-\mathbf{k}\cdot\mathbf{k}}[\mathbf{k}\cdot\mathbf{k}]^{-1}H_{n^{1}}(k_{1})d^{3}\mathbf{k}\\\
\\\
(vii)\;G^{\\#}(2n^{1},0,0,0,0,0;0)=\dfrac{(-1)^{n^{1}}}{\pi^{3/2}2^{n^{1}}\sqrt{(2n^{1})!}}\displaystyle\int\limits_{{\mathbb{R}}^{3}}e^{-\mathbf{k}\cdot\mathbf{k}}[\mathbf{k}\cdot\mathbf{k}]^{-1}H_{2n^{1}}(k_{1})d^{3}\mathbf{k}\\\
\\\ (viii)\;G^{\\#}(2n^{1}+1,0,0,0,0,0;0)=0\end{array}$ (69)

Figure 2: The negative of the usual singular Coulomb potential function
$W(x^{1};0)$ versus its non-singular counterpart $W^{\\#}(m;0)$

To reduce the notation, let us define
$W^{\\#}(n,\mu):=G^{\\#}(n,0,0;0,0,0;\mu)$, i.e.
$W^{\\#}(2n^{1},\mu)=G^{\\#}(2n^{1},0,0;0,0,0;\mu)$ and compare the discrete
formalism to the space-time continuum function implied by equation (42)

$W(x^{1};\mu)=G(x^{1},0,0;0,0,0;\mu)=\dfrac{1}{4\pi}\;\dfrac{e^{-\mu|x^{1}|}}{|x^{1}|}$
(70)

and the Coulomb potential

$W(x^{1};0)=G(x^{1},0,0;0,0,0;0)=\dfrac{1}{4\pi}\;\dfrac{1}{|x^{1}|}$ (71)

Note that the actual potential function $V$ between two Fermions interacting
via a Boson must include the coupling coefficient $g$ and $e$, i.e.

$V(x^{1};\mu)=W(x^{1};\mu)=\dfrac{e^{2}}{4\pi}\;\dfrac{e^{-\mu|x^{1}|}}{|x^{1}|}$
(72)

and

$V^{\\#}(n^{1};\mu)=-g^{2}W^{\\#}(n^{1},\mu)=\dfrac{-g^{2}\;i^{n^{1}}}{\pi^{3/2}2^{\frac{n^{1}}{2}}\sqrt{n^{1}!}}\displaystyle\int\limits_{{\mathbb{R}}^{3}}e^{-\mathbf{k}\cdot\mathbf{k}}[\mathbf{k}\cdot\mathbf{k}+\mu^{2}]^{-1}H_{n^{1}}(k_{1})d^{3}\mathbf{k}$
(73)

It is clear that continuum space-time function
$W(x^{1};\mu)=\frac{1}{4\pi}\;\frac{e^{-\mu|x^{1}|}}{|x^{1}|}$ from above has
a singularity at $x^{1}=0$ whereas the discrete phase space version

$W^{\\#}(2n^{1};0)=\dfrac{2^{(n^{1}+1)}\;(n^{1})!}{(2n^{1}+1)\sqrt{(2n^{1})!}}$
(74)

has no singularities even at $n^{1}\rightarrow 0_{+}$ as shown in figure 2. We
have already shown that $W^{\\#}(n^{1};\mu)$ is non-singular in equations (54)
and (55).

## Appendix: Discrete phase space-continuous time non-singular Green’s
functions for Fermionic-anti-Fermionic relativistic field equations

Recall that the relativistic Boson field equation in the discrete phase space-
continuous time formalism is given by

$\delta^{ab}\Delta_{(a)}^{\\#}\Delta_{(b)}^{\\#}\phi(\mathbf{n},t)-(\partial_{t})^{2}\phi(\mathbf{n},t)-\mu^{2}\phi(\mathbf{n},t)=0$
(75)

where $\mu>0$ is the mass of the Boson particle. The associated Green’s
function are given by DasIII ; DasIV ; DasV ,

$\begin{array}[]{c}\Delta_{(a)}^{\\#}(\mathbf{n},t;\mathbf{\hat{n}},\hat{t};\mu)=(2\pi)^{-1}\displaystyle\int_{{\mathbb{R}}^{3}}\left\\{\left[\prod_{j=1}^{3}\xi_{n^{j}}(k_{j})\overline{\xi_{\hat{n}^{j}}(k_{j})}\right]\right.\\\
\left.\displaystyle\int_{C_{(a)}}\left[\left(\eta^{\alpha\beta}k_{\alpha}k_{\beta}+\mu^{2}\right)^{-1}\mathbf{exp}(-ik_{4}(t-\hat{t})]dk^{4}\right]\right\\}d^{3}\mathbf{k}\\\
\\\
D_{(a)}(\mathbf{n},t;\hat{\mathbf{n}},\hat{t}):=\Delta_{(a)}^{\\#}(\mathbf{n},t;\hat{\mathbf{n}},\hat{t}:0)=\dfrac{1}{2\pi}\displaystyle\int_{{\mathbb{R}}^{3}}\left\\{\left[\prod_{j=1}^{3}\xi_{n^{j}}(k_{j})\overline{\xi_{\hat{n}^{j}}(k_{j})}\right]\right.\\\
\left.\displaystyle\int_{C_{(a)}}\left[\left(\eta^{\alpha\beta}k_{\alpha}k_{\beta}\right)^{-1}\mathbf{exp}(-ik_{4}(t-\hat{t})]dk^{4}\right]\right\\}d^{3}\mathbf{k}\end{array}$
(76)

The Green’s functions above involve nine contours in the complex $k^{4}$-plane
as exhibited explicitly in figure 3 with
$w=w(\mathbf{k}):=+\sqrt{\mathbf{k}\cdot\mathbf{k}+\mu^{2}}>0$.

Figure 3: Various contours $C_{(a)}$ in the complex $k^{4}$ plane. Figure 4:
Various contours $C_{(a)}$ in the complex $p^{4}$ plane.

Now, recall the Fermion-anti-Fermion wave equations in the discrete phase
space-continuous time formalism along with its Green’s function is furnished
by

$\begin{array}[]{c}\gamma^{a}\Delta^{\\#}_{(a)}\psi(\mathbf{n};t)+\gamma^{4}\partial_{t}\psi(\mathbf{n};t)+m\psi(\mathbf{n};t)=[0]_{(4\times
1)}\\\ \\\
\;[\Delta^{\\#}_{(a)}\tilde{\psi}(\mathbf{n};t)]\gamma^{a}+[\partial_{t}\tilde{\psi}(\mathbf{n};t)]\gamma^{4}-m\psi(\mathbf{n};t)=[0]_{(1\times
4)}\\\ \\\ S^{\\#}_{(a)}(\mathbf{n},t;\hat{\mathbf{n}},\hat{t};m)_{(4\times
4)}=(2\pi)^{-1}\displaystyle\int_{{\mathbb{R}}^{3}}\left\\{\left[\prod_{j=1}^{3}\xi_{n^{j}}(p_{j})\overline{\xi_{\hat{n}^{j}}(p_{j})}\right]\right.\\\
\left.\displaystyle\int_{C{(a)}}\dfrac{(i\gamma^{\mu}p_{\mu}-mI)e^{-ip^{4}(t-\hat{t})}}{\eta^{\alpha\beta}p_{\alpha}p_{\beta}+m^{2}}dp^{4}\right\\}d^{3}\mathbf{p}\end{array}$

The nine $4\times 4$ matrix Green’s functions above satisfy DasIII ; DasIV ;
DasV

$\begin{array}[]{c}(\gamma^{j}\Delta^{\\#}_{j}+\gamma^{4}\partial_{t}+mI)S^{\\#}_{(a)}(\mathbf{n},t;\hat{\mathbf{n}},\hat{t};m)_{(4\times
4)}=[0]_{(4\times 4)}\\\ \textbf{for contours}\,\,\,C,C_{(\pm)},C_{(1)}\\\ \\\
and\\\ \\\
=-\delta_{n^{1}\hat{n}^{1}}\delta_{n^{2}\hat{n}^{2}}\delta_{n^{3}\hat{n}^{3}}\delta(t-\hat{t})[I]_{(4\times
4)}\\\ \textbf{for
contours}\,\,\,C_{(R)},C_{(A)},C_{(P)},C_{(F+)},C_{(F-)}\end{array}$ (77)

with the contours depicted in figure 4.

## References

  * (1) H. Weyl,  Z. Physik. 46, 1, (1927).
  * (2) E. Wigner,  Phys. Rev. 40, 749, (1932).
  * (3) A. Das and P. Smozynski,  Found. Phys. Letter 7, 21 (1994).
  * (4) A. Das and P. Smozynski,  Found. Phys. Letter 7, 127 (1994).
  * (5) A. Das,  Can. J. Phys. 88, 73 (2010).
  * (6) A. Das,  Can. J. Phys. 88, 93 (2010).
  * (7) A. Das,  Can. J. Phys. 88, 111 (2010).
  * (8) J.M. Jauch and F. Rohrlich,  The Theory of Photons and Electrons,  Addison-Wesley, Cambridge, MA, (1955).
  * (9) M.E. Peskin and D.V. Schroeder,  An Introduction to Quantum Field Theory,  Addison-Wesley, Cambridge, MA, (1995).
  * (10) S. Weinberg  The Quantum Theory of Fields, I & II,  Cambridge University Press, Cambridge, MA, (1995).
  * (11) A. Das, R. Chatterjee, and T. Yu,  Mod. Phys. Lett. A 35, 24 (2020).
  * (12) C. Jordan,  Calculus of Finite Differences,  Chelsea Publishing, New York, (1965).
  * (13) I.S. Gradshteyn and I.M. Ryzhik  Table of Integrals, Series and Products,  Academic Press, San Diego, (1980).
  * (14) M.B. Green, J.H. Schwarz, and E. Witten,  Superstring Theory Vol 1 & 2,  Cambridge University Pres, Cambridge, (1987).
  * (15) A. Das, and R. Chatterjee,  Mod. Phys. Lett. A 36, 20 (2021).
  * (16) A. Das, and R. Chatterjee,  arXiv:2101.09191 (2021) (to appear in Mod. Phys. Lett. A)
  * (17) A. Das,  The Special Theory of Relativity (A mathematical exposition),  Springer-Verlag, Berlin, (1996).
  * (18) E.T. Whittaker, and G.N. Watson,  A Course of Modern Analysis,  Cambridge University Pres, Cambridge, (1961).