A Statistical Dirac Field Model on Minkowski Spacetime

Brenden McDearmon

doi:10.20944/preprints202504.1510.v1

Submitted:

17 April 2025

Posted:

18 April 2025

You are already at the latest version

Abstract

A statistical field model is constructed on Minkowski spacetime for a self-interacting Dirac field that fluctuates with respect to an auxiliary parameter. The fluctuating Dirac field undergoes a variational Hamiltonian flow while exchanging action with an action bath. Assuming ergodicity, the partition function is calculated and used to provide a perturbative expansion for the field's correlation functions. Non-perturbative calculations are carried out and the resulting 2-point correlation functions have a distinct light cone structure as expected for a relativistic field theory.

Keywords:

Relativistic statistical field theory

Subject:

Physical Sciences - Theoretical Physics

Introduction

Traditionally, non-perturbative quantum field theory calculations are performed on Euclidean lattices [1,2,3]. Physics, however, is relativistic. Euclidean lattice models lack a unique time-like direction and instead rely on analytical continuation to connect the Euclidean correlation functions to observables on Minkowski spacetime [3,4,5]. While several methods are known for generalizing these reconstructions to some curved spactimes, it would still be desirable to have models formulated on spacetime directly [6,7,8,9]. For this reason, various “real-time” methods have been developed to non-perturbatively calculate correlation functions directly on Minkowski spacetime. Many of these real-time models, however, face convergence issues because the actions of quantum field theories are usually not bounded from below and because fields are often weighted by a complex phase factor [10,11]. Other authors have addressed some of these issues using various relativistic statistical field theory constructions [12,13,14,15,16].

Recently, however, the present author has developed a statistical field model with two key features: (1) a relativistic action that is bounded from below and (2) an extended variational phase space whereby the fluctuating fields exchange action with an action bath [17]. Assuming ergodicity, the combination of these two features provides a relativistic partition function that samples a canonical ensemble of field configurations. The present article refines this approach and applies in to an interacting Dirac field model.

Preliminaries

The action of the free Dirac equation is typically expressed as

S_{f r e e}^{D i r a c} : = \int \bar{ψ} (x) (i \partial - m I) ψ (x) d^{4} x

where I is the 4x4 identity matrix,

\bar{ψ} : = ψ^{†} γ^{0}

, and

\partial : = \sum_{μ = 0}^{3} γ^{μ} \frac{\partial}{\partial x^{(μ)}}

with

γ^{μ}

denoting the Dirac representation gamma matrices and

x^{(μ)}

providing coordinates on Minkowski spacetime. To get the classical equations of motion for

ψ

, one usually sets the variational derivative of the action with respect to

\bar{ψ}

equal to zero; however, one could equivalently set the variational derivative with respect to

ψ^{†}

equal to zero to afford

γ^{0} (i \partial - m I) ψ (x) = 0

. Often, it is more efficient to deal with the momentum space version of the Dirac equation where, after Fourier transformation, the partial differential equations become the matrix equation

γ^{0} (k - m I) ψ (k) = 0

where now

k = γ^{0} k^{(0)} - \sum_{μ = 1}^{3} γ^{μ} k^{(μ)}

. The matrix

γ^{0} (k - m I)

is diagonalizable. Diagonalizing this matrix gives

γ^{0} (k - m) = P (k) \tilde{D} (k) P^{- 1} (k)

where

\tilde{D} (k)

is the diagonal matrix defined in equation (1) with

ω_{0} : = \sqrt{{(k^{(1)})}^{2} + {(k^{(2)})}^{2} + {(k^{(3)})}^{2} + m^{2}}

and

P (k)

is an invertible matrix.

\begin{matrix} \tilde{D} (k) : = (\begin{matrix} k^{(0)} - ω_{0} & 0 & 0 & 0 \\ 0 & k^{(0)} - ω_{0} & 0 & 0 \\ 0 & 0 & k^{(0)} + ω_{0} & 0 \\ 0 & 0 & 0 & k^{(0)} + ω_{0} \end{matrix}) \end{matrix}

(1)

The diagonalization of a matrix describes, in a different basis, the same operator as the original matrix. Redefining the Dirac field in this new basis as

ψ^{'} (k) : = P^{- 1} (k) ψ (k)

, the Dirac equation is now expressed as

\tilde{D} (k) ψ^{'} (k) = 0

. Using the inverse Fourier transformation then gives the spacetime representation of this diagonalized Dirac equation as the pseudodifferential equation

(D ψ^{'}) (x) : = \frac{1}{4 π^{2}} \int \tilde{D} (k) ψ^{'} (k) e^{- i k \cdot x} d^{4} k = 0

where

k \cdot x : = k^{(0)} x^{(0)} - k^{(1)} x^{(1)} - k^{(2)} x^{(2)} - k^{(3)} x^{(3)}

. This diagonalized Dirac operator is conceptually similar to the pseudodifferential operator square root of the Klein-Gordon equation described previously by Lämmerzahl [18]. Indeed, one can define an adjoint operator

D^{a d j}

that composes with D to give the Klein-Gordon equation. To see this, return to the free Dirac action and this time set the variational derivative with respect to

ψ

to zero to get

ψ^{†} (x) γ^{0} (i \overset{\leftarrow}{\partial} - m I) = 0

. Taking the Hermitian conjugate of both sides gives

- (i \partial + m I) γ^{0} ψ (x) = 0

. The momentum space representation of the matrix

- (i \partial + m I) γ^{0}

is again diagonalizable and affords the diagonal matrix defined in Equation (2).

\begin{matrix} {\tilde{D}}^{a d j} (k) : = (\begin{matrix} - k^{(0)} - ω_{0} & 0 & 0 & 0 \\ 0 & - k^{(0)} - ω_{0} & 0 & 0 \\ 0 & 0 & - k^{(0)} + ω_{0} & 0 \\ 0 & 0 & 0 & - k^{(0)} + ω_{0} \end{matrix}) \end{matrix}

(2)

Composing

{\tilde{D}}^{a d j} (k)

with

\tilde{D} (k)

gives

- ({(k^{(0)})}^{2} - {(k^{(1)})}^{2} - {(k^{(2)})}^{2} - {(k^{(3)})}^{2} - m^{2}) I

. It follows that

(D^{a d j} D ψ^{'}) (x) : = \frac{1}{4 π^{2}} \int {\tilde{D}}^{a d j} (k) \tilde{D} (k) ψ^{'} (k) e^{- i k \cdot x} d^{4} k = - (□ + m^{2}) ψ^{'} (x) = 0

as desired. Because the original basis for the Dirac field is unimportant, the tick mark, ′ , will be dropped from the notation going forward and

ψ (x)

will be understood to be a collection of four complex-valued fields

ψ^{(α)} (x) \in C

,

α = 0, 1, 2, or 3

.

These diagonalized free Dirac equations are interesting because one can define smooth, positive-definite integral kernels that are approximate solutions. For example, define

Δ_{m}^{D} (x; ϵ)

according to:

\begin{matrix} Δ_{m}^{D} (x; ϵ) : = (\begin{matrix} Δ_{m}^{+} (x; ϵ) & 0 & 0 & 0 \\ 0 & Δ_{m}^{+} (x; ϵ) & 0 & 0 \\ 0 & 0 & Δ_{m}^{-} (x; ϵ) & 0 \\ 0 & 0 & 0 & Δ_{m}^{-} (x; ϵ) \end{matrix}), \end{matrix}

(3)

where

ϵ

is a small positive number and

Δ_{m}^{\pm} (x; ϵ)

are defined as follows:

\begin{matrix} Δ_{m}^{\pm} (x; ϵ) : = \int \frac{1}{2 ω_{0} \sqrt{ϵ}} e x p (- \frac{π {(k^{(0)} \mp ω_{0})}^{2}}{ϵ} - ϵ ω_{0}) e^{- i k \cdot x} d^{4} k . \end{matrix}

(4)

The integral kernel

Δ_{m}^{D} (x; ϵ)

is an approximate solution of the diagonalized free Dirac operator D in the sense that

{lim}_{ϵ \to 0^{+}} (D Δ_{m}^{D}) (x; ϵ) = 0

because:

\begin{matrix} lim_{ϵ \to 0^{+}} \int \frac{k^{(0)} \mp ω_{0}}{2 ω_{0} \sqrt{ϵ}} e x p (- \frac{π {(k^{(0)} \mp ω_{0})}^{2}}{ϵ} - ϵ ω_{0}) e^{- i k \cdot x} d^{4} k = 0 \end{matrix}

(5)

since

\frac{1}{\sqrt{ϵ}} e x p (- \frac{π {(k^{(0)} \mp ω_{0})}^{2}}{ϵ})

approaches a Dirac delta distribution centered on

(k^{(0)} \mp ω_{0})

as

ϵ \to 0^{+}

.

Moreover,

Δ_{m}^{D} (x; ϵ)

is positive-definite because each of

Δ_{m}^{\pm} (x; ϵ)

is. That is, for any complex-valued function h, with complex conjugate

h^{*}

, that is not identically zero, these integral kernels have the following property:

\begin{matrix} positive - definitness for m, ϵ > 0 : \int \int h^{*} (x_{2}) Δ_{m}^{\pm} (x_{2} - x_{1}; ϵ) h (x_{1}) d^{4} x_{2} d^{4} x_{1} > 0 \end{matrix}

(6)

This property follows from the Plancherel theorem [19] whereby

\int \int h^{*} (x_{2}) Δ_{m}^{\pm} (x_{2} - x_{1}; ϵ) h (x_{1}) d^{4} x_{2} d^{4} x_{1}

is equal to

\int {\tilde{Δ}}_{m}^{\pm} (k; ϵ) {| \tilde{h} (k) |}^{2} d^{4} k

and

{\tilde{Δ}}_{m}^{\pm} (k; ϵ) > 0

everywhere when

m, ϵ > 0

.

The integral kernels also have smooth, closed form spacetime representations which can be found by explicitly evaluating the Fourier integrals. First, evaluate the integral with respect to

k^{(0)}

to get:

\begin{matrix} Δ_{m}^{\pm} (x; ϵ) = \frac{1}{{(2 π)}^{2}} e x p (- \frac{ϵ {(x^{(0)})}^{2}}{4 π}) \int \frac{1}{2 ω_{0}} e x p (- (ϵ \pm i x^{(0)}) ω_{0}) e^{i \sum_{j = 1}^{3} k^{(j)} x^{(j)}} d^{3} k . \end{matrix}

(7)

Switching to spherical coordinates with

ρ^{2} : = \sum_{j = 1}^{3} {(k^{(j)})}^{2}

,

t : = x^{(0)}

, and

r^{2} : = \sum_{j = 1}^{3} {(x^{(j)})}^{2}

gives:

\begin{matrix} Δ_{m}^{\pm} (t, r; ϵ) = \\ \frac{1}{4 π} e x p (- \frac{ϵ t^{2}}{4 π}) \int_{0}^{\infty} \frac{ρ^{2}}{\sqrt{ρ^{2} + m^{2}}} e x p (- (ϵ \pm i t) \sqrt{ρ^{2} + m^{2}}) \int_{- 1}^{1} e^{i ρ r \cos (θ)} d \cos (θ) d ρ . \end{matrix}

(8)

The integral with respect to the angle

θ

can be performed to give1:

\begin{matrix} Δ_{m}^{\pm} (t, r; ϵ) = \frac{1}{2 π r} e x p (- \frac{ϵ t^{2}}{4 π}) \int_{0}^{\infty} \frac{ρ}{\sqrt{ρ^{2} + m^{2}}} e x p (- (ϵ \pm i t) \sqrt{ρ^{2} + m^{2}}) \sin (ρ r) d ρ . \end{matrix}

(9)

According to the Bateman Table of Integral Transforms at page 75, equation 36 [21], this integral has the following closed-form solution:

\begin{matrix} Δ_{m}^{\pm} (t, r; ϵ) = \frac{m e^{- ϵ t^{2} / (4 π)}}{2 π \sqrt{r^{2} + {(ϵ \pm i t)}^{2}}} K_{1} [m \sqrt{r^{2} + {(ϵ \pm i t)}^{2}}], \end{matrix}

(10)

where

K_{1}

is the modified cylindrical Bessel function of the second kind and order one.

With these properties,

Δ_{m}^{D} (x; ϵ)

encodes the physics of the free Dirac equation into a smooth, positive-definite integral kernel with a closed-form spacetime representation.

Constructing the Statistical Dirac Field Model

Let X denote a simply connected region of Minkowski spacetime having a boundary,

\partial X

, and finite 4-volume,

{Vol}_{X} : = \int_{X} d^{4} x

. The Dirac field,

ψ

, is a collection of four complex-valued scalar fields,

ψ^{(α)}

, on X that is typically arranged in column vector notation as:

\begin{matrix} ψ (x) : = (\begin{matrix} ψ^{(0)} (x) \\ ψ^{(1)} (x) \\ ψ^{(2)} (x) \\ ψ^{(3)} (x) \end{matrix}) . \end{matrix}

(11)

Suppose that this Dirac field “fluctuates” with respect to some parameter

λ \in [0, \infty)

. The variation of

ψ

for any given value of

(x, λ) \in (X, [0, \infty))

is defined as

\dot{ψ} (x, λ) : = {\frac{d ψ (x, λ + κ)}{d κ}|}_{κ = 0}

. To prescribe physics to these fluctuations, it is useful to couple

ψ

to an additional fluctuating field defined as:

\begin{matrix} π^{(ψ)} (x) : = (\begin{matrix} π_{(ψ)}^{(0)} (x) \\ π_{(ψ)}^{(1)} (x) \\ π_{(ψ)}^{(2)} (x) \\ π_{(ψ)}^{(3)} (x) \end{matrix}), \end{matrix}

(12)

where each

π_{(ψ)}^{(α)} (x)

is a complex-valued scalar field. Fluctuations of

π^{(ψ)}

are defined by

{\dot{π}}^{(ψ)} (x, λ) : = {\frac{d π^{(ψ)} (x, λ + κ)}{d κ}|}_{κ = 0}

. As these fields fluctuate, they can exchange action with an action bath by interacting with the action bath parameter

s (λ)

which is a global, real-valued parameter with

\dot{s} (λ) : = {\frac{d s (λ + κ)}{d κ}|}_{κ = 0}

. The action bath parameter

s (λ)

will likewise be coupled to an additional parameter denoted by

π^{(s)} (λ)

with

{\dot{π}}^{(s)} (λ) : = {\frac{d π^{(s)} (λ + κ)}{d κ}|}_{κ = 0}

. Together, these fields define the degrees of freedom of the system. To define how these fields fluctuate consider the following action:

\begin{matrix} S : = s (\frac{{(π^{(s)})}^{2}}{2 M_{s}} + 8 \frac{{Vol}_{X}}{β} \ln (s) + \int_{X_{2}} \frac{| π^{(ψ)} |^{2} (x_{2})}{s^{2}} + L [ψ (x_{2})] d^{4} x_{2} - S_{0}), \end{matrix}

(13)

\begin{matrix} with L [ψ (x_{2})] : = ψ^{†} (x_{2}) \int_{X_{1}} {[Δ_{m}^{D}]}^{- 1} (x_{2} - x_{1}; ϵ) ψ (x_{1}) + q^{2} V [ψ (x_{2}), ψ (x_{1})] d^{4} x_{1} . \end{matrix}

(14)

Here,

S_{0}

,

M_{s}

,

β

, and q are real-valued constants that do not depend on

λ

and the explicit

λ

-dependence of the degrees of freedom is omitted in favor of a more compact notation. The constant

S_{0}

in particular is defined such that

{S|}_{λ = 0} = 0

. In

L

, the term

{[Δ_{m}^{D}]}^{- 1} (x_{2} - x_{1}; ϵ)

is the operator inverse of

Δ_{m}^{D} (x_{2} - x_{1}; ϵ)

and necessarily exists and is positive-definite because

Δ_{m}^{D} (x_{2} - x_{1}; ϵ)

is positive-definite. The operator inverse is used here so that the 2-point correlation function of the free field is proportional to

Δ_{m}^{D} (x_{2} - x_{1}; ϵ)

as shown in one of the sections below. Additionally, the term

V [ψ (x_{2}), ψ (x_{1})]

is an interaction term for the Dirac field, an example for which is provided in the examples section below. Overall, the term

L [ψ (x_{2})]

must be bounded from below, and this limitation constrains the set of acceptable interaction terms

V [ψ (x_{2}), ψ (x_{1})]

.

The development of the degrees of freedom with respect to

λ

is then defined by the following variational analog of Hamilton’s differential equations:

\begin{matrix} {\dot{π}}^{(ψ)} (x, λ) = - {\frac{δ S}{δ ψ^{†}}|}_{x, λ}, \end{matrix}

(15)

\begin{matrix} \dot{ψ} (x, λ) = {\frac{δ S}{δ {π^{(ψ)}}^{†}}|}_{x, λ}, \end{matrix}

(16)

\begin{matrix} {\dot{π}}^{(s)} (λ) = - {\frac{\partial S}{\partial s}|}_{λ}, and \end{matrix}

(17)

\begin{matrix} \dot{s} (λ) = {\frac{\partial S}{\partial π^{(s)}}|}_{λ} . \end{matrix}

(18)

Here,

\frac{δ S}{δ (\cdot)}

denotes the variational derivative. These equations define a symplectic flow on the variational phase space consisting of the degrees of freedom. This flow conserves the total action of the system, assuming appropriate boundary conditions, because the total action acts as a variational Hamiltonain [17,22,23,24]. When the fields are discretized for numerical simulation, these variational equations become a finite system of partial differential equations that can be integrated numerical with respect to

λ

as demonstrated in the examples section below. Thus, assuming a suitable prescription for all of the degrees of freedom at

λ = 0

, one can obtain a 1-parameter family of system configurations from which observables can be calculated.

Observables, Expectation Values, and the Partition Function

The observables of interest for this system are functionals of Dirac field, e.g.,

O [ψ^{(α)}] \in C

. The expectation value of an observable is given by its average with respect to

λ

. Specifically,

{〈O [ψ^{(α)}]〉}_{λ} : = {lim}_{λ \to \infty} \frac{1}{Λ} \int_{0}^{Λ} O [ψ^{(α)} (λ)] d λ

. These expectation values can be computed during numerical simulations. Assuming ergodicity of the system’s dynamics with respect to

λ

, one can also calculate expectation values using the partition function.

The partition function is defined by the ergodic average over the variational phase space using the phase space measure given by:

\begin{matrix} d Γ : = C δ (S) d s d π^{(s)} D [R e {π^{(ψ)}}] D [I m {π^{(ψ)}}] D [R e {ψ}] D [I m {ψ}], \end{matrix}

(19)

where

D [R e {ψ}] : = \prod_{x \in {1, 2 . . . N}, α \in {0, 1, 2, 3}} d R e {ψ^{(α)}} (x)

is the functional integral measure over N discretized points, the other functional integral measures are defined analogously, and C is a normalization constant. Define the following change of variables:

{\underset{̲}{π}}^{(ψ)} : = \frac{π^{(ψ)}}{s}

and

{\underset{̲}{π}}^{(ψ) †} : = \frac{π^{(ψ) †}}{s}

. This change of variables is equivalent to a re-scaling of

λ

. See the article by Bond, Leimkuhler, and Laird for an analogous construction used in classical statistical mechanics [25]. This change of variables changes the integration measure in the partition function by a factor of

s^{8 N}

to give:

\begin{matrix} Z & = C \int s^{8 N} δ (S) d s d π^{(s)} D [R e {{\underset{̲}{π}}^{(ψ)}}] D [I m {{\underset{̲}{π}}^{(ψ)}}] D [R e {ψ}] D [I m {ψ}] . \end{matrix}

(20)

Recall the identity

\frac{d}{d c} δ [f (c)] = \frac{δ [f (c - c^{'})]}{(\frac{d f}{d c} |_{(c^{'})})}

, where

c^{'}

is the isolated zero of

f (c)

. For the systems of interest here,

c^{'}

is given by the following equation:

\begin{matrix} c^{'} = e x p [\frac{- β}{8 V o l_{X}} (\frac{{(π^{(s)})}^{2}}{2 M_{s}} + \int_{X_{2}} {| {\underset{̲}{π}}^{(ψ)} |}^{2} (x_{2}) + L [ψ (x_{2})] d^{4} x_{2} - S_{0})] . \end{matrix}

(21)

Assigning a unit volume to each of the N points in the functional integral measure gives

V o l_{X} = N

. Performing the integral with respect to

d s

affords:

\begin{matrix} Z \propto \int e x p [- β (\frac{{(π^{(s)})}^{2}}{2 M_{s}} + & \int_{X_{2}} {| {\underset{̲}{π}}^{(ψ)} |}^{2} (x_{2}) + L [ψ (x_{2})] d^{4} x_{2} - S_{0})] \\ d π^{(s)} D [R e {{\underset{̲}{π}}^{(ψ)}}] D [I m {{\underset{̲}{π}}^{(ψ)}}] D [R e {ψ}] D [I m {ψ}] . \end{matrix}

(22)

Now, the integrals with respect to

d π^{(s)}

,

D [R e {{\underset{̲}{π}}^{(ψ)}}]

,

D [I m {{\underset{̲}{π}}^{(ψ)}}]

are all Gaussian and can be evaluated and, after collecting constant terms, affords the partition function:

\begin{matrix} Z = Z_{0} \int e x p [- β \int_{X_{2}} L [ψ (x_{2})] d^{4} x_{2}] D [R e {ψ}] D [I m {ψ}] . \end{matrix}

(23)

The partition function is helpful for calculating expectation values because, by the assumption of ergodicity, the expectation value of some observable

O

taken with respect to

λ

agrees with its ensemble average. In equations:

{〈O [ψ^{(α)}]〉}_{λ} : = lim_{Λ \to \infty} \frac{1}{Λ} \int_{0}^{Λ} O [ψ^{(α)} (λ)] d λ,

(24)

\begin{matrix} {〈O [ψ^{(α)}]〉}_{Γ} : = Z_{0} \int O [ψ^{(α)}] e x p [- β \int_{X_{2}} L [ψ (x_{2}))] d^{4} x_{2}] D [R e {ψ}] D [I m {ψ}], and \end{matrix}

(25)

\begin{matrix} {〈O [ψ^{(α)}]〉}_{λ} = {〈O [ψ^{(α)}]〉}_{Γ} . \end{matrix}

(26)

Because the expectation value with respect to

λ

agrees with the ensemble average, let

〈O [ψ^{(α)}]〉

denote the expectation value determined by either method.

Correlation Functions

The correlation functions encode the physics of the system. These correlation functions can be calculated perturbatively using the partition function, as demonstrated in this section, or non-perturbatively by numerically integrating the system’s evolution with respect to

λ

, as demonstrated in the next section. To begin, consider the cumulant generating function:

\begin{matrix} \ln (Z [η]) : = \\ \ln (Z_{0} \int e x p [\int_{X_{2}} i η^{†} (x_{2}) ψ (x_{2}) - i ψ^{†} (x_{2}) η (x_{2}) - β L [ψ (x_{2})] d^{4} x_{2}] D [R e {ψ}] D [I m {ψ}]), \end{matrix}

(27)

where

η

is a source term, and recall the definition of

L [ψ (x_{2})]

:

\begin{matrix} L [ψ (x_{2})] : = ψ^{†} (x_{2}) \int_{X_{1}} {[Δ_{m}^{D}]}^{- 1} (x_{2} - x_{1}; ϵ) ψ (x_{1}) + q^{2} V [ψ (x_{2}), ψ (x_{1})] d^{4} x_{1} . \end{matrix}

(28)

The correlation functions are calculated in terms of the cumulant generating function by taking variational derivatives with respect to

η^{†}

and

η

followed by setting the source terms to zero:

\begin{matrix} {\frac{δ^{a + b} \ln (Z [η])}{δ η^{†} (x_{1}) . . . δ η^{†} (x_{a}) δ η (x_{a + 1}) . . . δ η (x_{a + b})}|}_{η^{†}, η = 0} = 〈{(- 1)}^{b} i^{a + b} ψ (x_{1}) . . . ψ (x_{a}) ψ^{†} (x_{a + 1}) . . . ψ^{†} (x_{a + b})〉 . \end{matrix}

(29)

These correlation functions can be expressed perturbatively in the interaction term,

q^{2} V [ψ (x_{2}), ψ (x_{1})]

, by expanding them as a series in the parameter q. The generating function can be compactly expressed by making use of the following indentity terms of a variational operator acting on the free partition function as follows:

\begin{matrix} Z [η] = e x p [- i β q^{2} \int_{X_{2}} \int_{X_{1}} V [- i \frac{δ}{δ η (x_{2})}, - i \frac{δ}{δ η (x_{1})}] d^{4} x_{2} d^{4} x_{1}] Z_{f r e e} [η] \end{matrix}

(30)

where the details of this expression can be found in Bailin and Love at Chaper 5 [26]. The free partition function, defined as:

\begin{matrix} Z_{f r e e} [η] : = Z_{0} \int & e x p [\int_{X_{2}} i η^{†} (x_{2}) ψ (x_{2}) - i ψ^{†} (x_{2}) η (x_{2}) \\ - β ψ^{†} (x_{2}) \int_{X_{1}} {[Δ_{m}^{D}]}^{- 1} (x_{2} - x_{1}; ϵ) ψ (x_{1}) d^{4} x_{2} d^{4} x_{1}] D [R e {ψ}] D [I m {ψ}], \end{matrix}

(31)

is a standard Gaussian integral which can be evaluated to give:

\begin{matrix} Z_{f r e e} [η] \propto e x p [- \frac{1}{β} \int_{X_{2}} \int_{X_{1}} η^{†} (x_{2}) Δ_{m}^{D} (x_{2} - x_{1}; ϵ) η (x_{1}) d^{4} x_{2} d^{4} x_{1}] . \end{matrix}

(32)

Combining the above equations, the correlation functions can be expressed in terms of various contractions of the integral kernel

Δ_{m}^{D} (x_{2} - x_{1}; ϵ)

, and one could use a diagrammatic technique similar to the Feynman diagrams to keep track of these terms. In particular, this shows that the 2-point correlation function of the free field is given by:

\begin{matrix} {〈ψ (x_{2}) ψ^{†} (x_{1})〉}_{f r e e} = \frac{1}{β} Δ_{m}^{D} (x_{2} - x_{1}; ϵ) . \end{matrix}

(33)

Examples

An advantage of the model developed in the previous sections is that it is easy to discretize and simulate numerically. In each of the following examples, the field

ψ

was evaluated on a set of 4096 points selected at random from the region X of Minkowski spacetime that satisfies

\sqrt{{(x_{0})}^{2}} < 1

and

\sqrt{{(x_{1})}^{2} + {(x_{2})}^{2} + {(x_{3})}^{2}} + \sqrt{{(x_{0})}^{2}} < 1

for each

x \in X

with a uniform distribution of points in the (t, r) plane. Each component of

ψ

and

π^{(ψ)}

was initialized with a random phase and magnitude ranging from -1.75 to 1.75. The variables s and

π^{(s)}

where initialized to 1 and 0, respectively. The parameter values

m = 1

and

β = 1

were used. The value of

ϵ

was set to 0.01 for the integral kernel,

Δ_{m}^{D} (x_{2} - x_{1}; ϵ)

, and the integral kernel was divided by its value at

x_{2} - x_{1} = 0

to standardize the results.

Once an initial configuration was established, each example was numerically integrated with respect to

λ

using the explicit leap-frog algorithm provided by Bond, Leimkuhler, and Laird [25] with a step size of

Δ λ = 0.1

. Each system was “equilibrated" by stepping the system forward for 1,000,000 steps. After this “equilibration” period, each system was further numerically integrated with respect to

λ

for an additional 5,000,000 steps during which data was collected. The collected data included expectation values of the form

〈\bar{ψ} (f_{i}) ψ (f_{0})〉

where

\bar{ψ} (f_{i}) : = \sum_{x \in X} \bar{ψ} (x) f_{i} (x)

and

ψ (f_{0}) : = \sum_{x \in X} ψ (x) f_{0} (x)

with

{f_{i} (x)}

being a set of 225 equally sized, real-valued functions that tesselate X by taking the value 1 in a triangularly supported region of the (t, r) plane and taking the value 0 outside of each function’s respective region of support. Further,

f_{0} (x)

covers the origin where (t, r)=(0, 0).

Example 1: Free Field

In this example, the free field model was used (i.e., q=0). The real and imaginary parts of the expectation value

〈\bar{ψ} (f_{i}) ψ (f_{0})〉

were each normalized to a maximum absolute value of 1 and are ploted in Figure 1A and 1B, respectively. From the calculation of the 2-point correlation function in the previous section, one would expect

〈\bar{ψ} (f_{i}) ψ (f_{0})〉

to approximate

〈{(ψ^{(0)} (x_{2}))}^{*} ψ^{(0)} (x_{1})〉 + 〈{(ψ^{(1)} (x_{2}))}^{*} ψ^{(1)} (x_{1})〉 - 〈{(ψ^{(2)} (x_{2}))}^{*} ψ^{(2)} (x_{1})〉 - 〈{(ψ^{(3)} (x_{2}))}^{*} ψ^{(3)} (x_{1})〉

which is to proportional to

I m {Δ_{m}^{+} (x_{2} - x_{1}; ϵ)}

. Moreover,

l i m_{ϵ \to 0^{+}} I m {Δ_{m}^{+} (x_{2} - x_{1}; ϵ)} = Δ_{m} (x_{2} - x_{1})

which is the Pauli-Jordon propagator.

As seen in the plot of the imaginary part (Figure 1B), the correlation function has three distinct regions: (I) the interior of the past-directed light cone, (II) the acausal region, and (III) the interior of the future-direct light cone. Regions, (I) and (III) have a small oscillatory character whereas region (II) is substantially zero. Moreover, region (II) is separated from regions (I) and (III) by a sharp light cone with the past-directed light cone and future-direct light cone being antisymmetric with respect to time inversion. All of these properties are consistent with the expectation that the result should be similar to the Pauli-Jordon propagator.

The real-part of the correlation (Figure 1A), by contrast, is more noisy with a large negative feature near the origin and some smaller oscillatory behavior in the interiors of the future and past directed light cones. These effects may be due to the non-interacting nature of the separate

ψ^{(α)}

modes which do not cancel each other out to give zero. This hypothesis is supported in the following interacting model wherein the real part of the correlation function is substantially feature-less noise.

Example 2: Interacting Model

In this model, the following interaction term was used with q=0.1:

\begin{matrix} q^{2} V [ψ (x_{2}), ψ (x_{1})] = q^{2} {(j^{μ} (x_{2}))}^{†} Δ^{S} (x_{2} - x_{1}; ϵ) j^{μ} (x_{1}), where \end{matrix}

(34)

\begin{matrix} j^{μ} (x_{2}) : = \bar{ψ} (x_{2}) γ^{μ} ψ (x_{2}), and \end{matrix}

(35)

\begin{matrix} Δ^{S} (x_{2} - x_{1}; ϵ) : = d i a g (Δ^{T} (x_{2} - x_{1}; ϵ), Δ^{T} (x_{2} - x_{1}; ϵ), {Δ^{T}}^{*} (x_{2} - x_{1}; ϵ), {Δ^{T}}^{*} (x_{2} - x_{1}; ϵ)), \end{matrix}

(36)

\begin{matrix} Δ^{T} (x_{2} - x_{1}; ϵ) : = (Δ_{ϵ}^{+} (x_{2}^{(0)} - x_{1}^{(0)}, {\vec{x}}_{2} - {\vec{x}}_{1}; ϵ) + Δ_{ϵ}^{+} (x_{1}^{(0)} - x_{2}^{(0)}, {\vec{x}}_{2} - {\vec{x}}_{1}; ϵ)) / Δ_{m}^{+} (0, ϵ) . \end{matrix}

(37)

Here, the

j^{μ} (x_{2})

and

j^{μ} (x_{1})

terms are currents generated by the fluctuating Dirac field and they interact through the time-symmetrized integral kernel

Δ^{T} (x_{2} - x_{1}; ϵ)

and its complex conjugate

{Δ^{T}}^{*} (x_{2} - x_{1}; ϵ)

where

ϵ

additionally parameterizes a small but positive mass term that, in the limit as

ϵ

goes to zero, may mimic the exchange of a photon. The interaction kernel

Δ^{T} (x_{2} - x_{1}; ϵ)

was found to be positive-definite by Cholesky decomposition for the system under consideration; so, other interaction models could consider using an inverse of this kernel. The factor

Δ_{m}^{+} (0, ϵ)

was included to standardize the magnitude.

As seen in the plot of the imaginary part (Figure 2B), the correlation function has three distinct regions: (I) the interior of the past-directed light cone, (II) the acausal region, and (III) the interior of the future-direct light cone. Regions, (I) and (III) have a small oscillatory character whereas region (II) is substantially zero. Moreover, region (II) is separated from regions (I) and (III) by a sharp light cone with the past-directed light cone and future-direct light cone being antisymmetric with respect to time inversion. These properties suggest that this interacting model may obey causality as would be required for a physically meaningful model.

The real-part of the correlation (Figure 2A), by contrast, looks like feature-less noise. While both plots were normalized to have a greatest absolute value of 1, before normalization the maximum absolute value for the real part was about two orders of magnitude smaller than that of the imaginary part. This suggests that the real part may be effectively zero.

Discussion

These results demonstrate an implementation of a statistical Dirac field model on Minkowski spacetime. The 2-point correlation functions of both the free field model and the interacting model have light cone structures similar to what one would require from a relativistic field theory. Future work may involve attempting to construct a fermionic algebra of observables with the appropriate anti-commutation relations.

References

Roepstorff, G. Path integral approach to quantum physics: an introduction; Springer Science & Business Media, 2012.
Montvay, I.; Münster, G. Quantum fields on a lattice; Cambridge University Press, 1994.
Glimm, J.; Jaffe, A. Quantum physics: a functional integral point of view; Springer Science & Business Media, 2012.
Osterwalder, K.; Schrader, R. Axioms for Euclidean Green’s functions. Communications in mathematical physics 1973, 31, 83–112. [CrossRef]
Osterwalder, K.; Schrader, R. Axioms for Euclidean Green’s functions II. Communications in mathematical physics 1975, 42, 281–305. [CrossRef]
Ashtekar, A.; Marolf, D.; Mourao, J.; Thiemann, T. Osterwalder-Schrader reconstruction and diffeomorphism invariance. arXiv preprint quant-ph/9904094 1999.
Visser, M. How to Wick rotate generic curved spacetime. arXiv preprint arXiv:1702.05572 2017.
Baldazzi, A.; Percacci, R.; Skrinjar, V. Wicked metrics. Classical and Quantum Gravity 2019, 36, 105008. [CrossRef]
Schrader, R. Reflection positivity in simplicial gravity. Journal of Physics A: Mathematical and Theoretical 2016, 49, 215202. [CrossRef]
Berges, C.E.; Rammelmuller, L.; Loheac, A.C.; Ehmann, F.; Braun, J.; Drut, J.E. Complex Langevin and other approaches to the sign problem in quantum many-body physics. Physics Reports 2021, 892, 1–54. [CrossRef]
Berges, J.; Borsanyi, S.; Sexty, D.; Stamatescu, I.O. Lattice simulations of real-time quantum fields. Physical Review D 2007, 75, 045007. [CrossRef]
Giachello, M.; Gradenigo, G. Symplectic quantization and the Feynman propagator: a new real-time numerical approach to lattice field theory. arXiv preprint arXiv:2403.17149 2024.
Tilloy, A. Interacting quantum field theories as relativistic statistical field theories of local beables. arXiv preprint arXiv:1702.06325 2017.
Floerchinger, S. Quantum field theory as a bilocal statistical field theory. arXiv preprint arXiv:1004.2250 2010.
Hitschfeld, B.S. Probability Density Functionals in Perturbation Theory within the QFT in-in Formalism. arXiv preprint arXiv:1901.09734 2019.
Kent, A. Path integrals and reality. arXiv preprint arXiv:1305.6565 2013.
McDearmon, B. A relativistic statistical field theories on a discretized Minkoski lattice with quantum-like observables. arXiv preprint arXiv:2405.2038 2024.
Lämmerzahl, C. The pseudodifferential operator square root of the Klein-Gordon equation. Journal of mathematical physics 1993, 34, 3918–3932. [CrossRef]
Rudin, W. Functional Analysis; McGraw-Hill, 1973.
Scharf, G. Finite quantum electrodynamics: the causal approach; Courier Corporation, 2014.
Bateman, H.; Erdélyi, A.; Magnus, W.; Oberhettinger, F. Tables of integral transforms 1954. 1.
McDearmon, B. A momentum space formulation for some relativistic statistical field theories with quantum-like observables. arXiv preprint arXiv:2406.04365 2024.
McDearmon, B. Euclidean Quantum Field Theory from Variational Dynamics. arXiv arXiv:2303.12666 2023.
McDearmon, B. Euclidean Quantum Gravity from Variational Dynamics. arXiv arXiv:2302.05713 2023.
Bond, S.D.; Leimkuhler, B.J.; Laird, B.B. The Nose-Poincare method for constant temperature molecular dynamics. Journal of Computational Physics 1999, 151, 114–134. [CrossRef]
Bailin, D.; Love, A. Introduction to Gauge Field Theory Revised Edition; Taylor & Francis, 1993.

1	A similar integral is evaluated in spherical coordinates in the book by Sharf at equations (2.3.8) and (2.3.9) [20].

Figure 1. Normalized real part (A) and imaginary part (B) of the correlation function

〈\bar{ψ} (f_{i}) ψ (f_{0})〉

of the free field model in the (t, r) plane. Before normalization, the maximum absolute value for the real part was 0.8 times that of the imaginary part.

Figure 1. Normalized real part (A) and imaginary part (B) of the correlation function

〈\bar{ψ} (f_{i}) ψ (f_{0})〉

of the free field model in the (t, r) plane. Before normalization, the maximum absolute value for the real part was 0.8 times that of the imaginary part.

Figure 2. Normalized real part (A) and imaginary part (B) of the correlation function

〈\bar{ψ} (f_{i}) ψ (f_{0})〉

of the interacting field model in the (t, r) plane. Before normalization, the maximum absolute value for the real part was 0.018 times that of the imaginary part.

Figure 2. Normalized real part (A) and imaginary part (B) of the correlation function

〈\bar{ψ} (f_{i}) ψ (f_{0})〉

of the interacting field model in the (t, r) plane. Before normalization, the maximum absolute value for the real part was 0.018 times that of the imaginary part.

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

© 2025 by the author. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

Copyright: This open access article is published under a Creative Commons CC BY 4.0 license, which permit the free download, distribution, and reuse, provided that the author and preprint are cited in any reuse.

A Statistical Dirac Field Model on Minkowski Spacetime

Abstract

Keywords:

Subject:

Introduction

Preliminaries

Constructing the Statistical Dirac Field Model

Observables, Expectation Values, and the Partition Function

Correlation Functions

Examples

Example 1: Free Field

Example 2: Interacting Model

Discussion

References

MDPI Initiatives

Important Links

Subscribe