%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                                                                 %%
%% Please do not use \input{...} to include other tex files.       %%
%% Submit your LaTeX manuscript as one .tex document.              %%
%%                                                                 %%
%% All additional figures and files should be attached             %%
%% separately and not embedded in the \TeX\ document itself.       %%
%%                                                                 %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%\documentclass[referee,sn-basic]{sn-jnl}% referee option is meant for double line spacing

%%=======================================================%%
%% to print line numbers in the margin use lineno option %%
%%=======================================================%%

%%\documentclass[lineno,sn-basic]{sn-jnl}% Basic Springer Nature Reference Style/Chemistry Reference Style

%%======================================================%%
%% to compile with pdflatex/xelatex use pdflatex option %%
%%======================================================%%

%%\documentclass[pdflatex,sn-basic]{sn-jnl}% Basic Springer Nature Reference Style/Chemistry Reference Style

%%\documentclass[sn-basic]{sn-jnl}% Basic Springer Nature Reference Style/Chemistry Reference Style
\documentclass[pdflatex,sn-mathphys]{sn-jnl}% Math and Physical Sciences Reference Style
%%\documentclass[sn-aps]{sn-jnl}% American Physical Society (APS) Reference Style
%%\documentclass[sn-vancouver]{sn-jnl}% Vancouver Reference Style
%%\documentclass[sn-apa]{sn-jnl}% APA Reference Style
%%\documentclass[sn-chicago]{sn-jnl}% Chicago-based Humanities Reference Style
%%\documentclass[sn-standardnature]{sn-jnl}% Standard Nature Portfolio Reference Style
%%\documentclass[default]{sn-jnl}% Default
%%\documentclass[default,iicol]{sn-jnl}% Default with double column layout

%%%% Standard Packages
%%<additional latex packages if required can be included here>
%%%%

%%%%%=============================================================================%%%%
%%%%  Remarks: This template is provided to aid authors with the preparation
%%%%  of original research articles intended for submission to journals published 
%%%%  by Springer Nature. The guidance has been prepared in partnership with 
%%%%  production teams to conform to Springer Nature technical requirements. 
%%%%  Editorial and presentation requirements differ among journal portfolios and 
%%%%  research disciplines. You may find sections in this template are irrelevant 
%%%%  to your work and are empowered to omit any such section if allowed by the 
%%%%  journal you intend to submit to. The submission guidelines and policies 
%%%%  of the journal take precedence. A detailed User Manual is available in the 
%%%%  template package for technical guidance.
%%%%%=============================================================================%%%%

\jyear{2022}%

%% as per the requirement new theorem styles can be included as shown below
\theoremstyle{thmstyleone}%
\newtheorem{theorem}{Theorem}%  meant for continuous numbers
%%\newtheorem{theorem}{Theorem}[section]% meant for sectionwise numbers
%% optional argument [theorem] produces theorem numbering sequence instead of independent numbers for Proposition
\newtheorem{proposition}[theorem]{Proposition}% 
%%\newtheorem{proposition}{Proposition}% to get separate numbers for theorem and proposition etc.

\theoremstyle{thmstyletwo}%
\newtheorem{example}{Example}%
\newtheorem{remark}{Remark}%

\theoremstyle{thmstylethree}%
\newtheorem{definition}{Definition}%

\raggedbottom
%%\unnumbered% uncomment this for unnumbered level heads

\begin{document}

%%\title{Brownian-Huygens propagation: modeling wave functions with discrete particle-antiparticle random walks}

\title[Brownian-Huygens propagation]{Brownian-Huygens propagation: modeling wave functions with discrete particle-antiparticle random walks}

%%=============================================================%%
%% Prefix	-> \pfx{Dr}
%% GivenName	-> \fnm{Joergen W.}
%% Particle	-> \spfx{van der} -> surname prefix
%% FamilyName	-> \sur{Ploeg}
%% Suffix	-> \sfx{IV}
%% NatureName	-> \tanm{Poet Laureate} -> Title after name
%% Degrees	-> \dgr{MSc, PhD}
%% \author*[1,2]{\pfx{Dr} \fnm{Joergen W.} \spfx{van der} \sur{Ploeg} \sfx{IV} \tanm{Poet Laureate} 
%%                 \dgr{MSc, PhD}}\email{iauthor@gmail.com}
%%=============================================================%%


\author*[1,2]{\fnm{Hrvoje J.} \sur{Hrgov\v{c}i\'{c}}}\email{hrgovcic@alum.mit.edu}


%%\affil*[1]{\orgdiv{Department}, \orgname{Organization}, \orgaddress{\street{Street}, \city{City}, \postcode{100190}, \state{State}, \country{Country}}}

%%\affil[2]{\orgdiv{Department}, \orgname{Organization}, \orgaddress{\street{Street}, \city{City}, \postcode{10587}, \state{State}, \country{Country}}}

%%\affil[3]{\orgdiv{Department}, \orgname{Organization}, \orgaddress{\street{Street}, \city{City}, \postcode{610101}, \state{State}, \country{Country}}}

%%==================================%%
%% sample for unstructured abstract %%
%%==================================%%
\newcommand{\etal}{\textit{et al}.}
\newcommand{\ie}{\textit{i}.\textit{e}.}
\newcommand{\eg}{\textit{e}.\textit{g}.}
\newcommand{\cf}{\textit{cf}.}

\abstract{We present a simple method of discretely modeling solutions of the classical wave and Klein-Gordon equations using variations of random walks on a graph. Consider a collection of particles executing random walks on an undirected bipartite graph embedded in $\mathbb{R}^D$ at discrete times $\mathbb{Z}$, and assume those walks are “heat-like”, in the sense that a (conserved) density of particles obeys the heat equation in the continuum limit.

If the particles possess a binary degree of freedom (so that they may be said to be either particles or “antiparticles”, \ie, “positive” or “negative), then there exist closely related $\textit{branching}$ random walks on the same graph that are “wave-like”, in the sense that their (also conserved) net density obeys the $D$-dimensional classical wave equation.

Such wave-like paths can be generated even on random graphs.

The transformation by which the heat-like random walks become wave-like branching random walks is as follows: at every time step, any “incoming” particle arriving at any node $X$ of the graph along some edge $e_{X}$ creates a particle-antiparticle pair at that node, with the stipulation that the newly created particle with the opposite sign of the incoming particle must initially step (in “Huygens” back-propagation fashion) along $e_{X}$ in the $\textit{reverse}$ direction, while the other two take a step (in “Brownian” fashion) along any edge of $X$ (including possibly $e_{X}$) with equal likelihood, and chosen independently. 

An additional degree of freedom (resulting in “bra” and “ket” particles) leads to quasi-probability densities proportional to the square of the wave functions.}
%%================================%%
%% Sample for structured abstract %%
%%================================%%

% \abstract{\textbf{Purpose:} The abstract serves both as a general introduction to the topic and as a brief, non-technical summary of the main results and their implications. The abstract must not include subheadings (unless expressly permitted in the journal's Instructions to Authors), equations or citations. As a guide the abstract should not exceed 200 words. Most journals do not set a hard limit however authors are advised to check the author instructions for the journal they are submitting to.


\keywords{ random walk, wave propagation, graph theory, cellular automata, discrete quantum fields}

%%\pacs[JEL Classification]{D8, H51}

%%\pacs[MSC Classification]{35A01, 65L10, 65L12, 65L20, 65L70}

\maketitle



%%\section{Introduction}\label{sec1}
\section{\label{sec:level1}Introduction: Heat-like and wave-like walks} 

Consider a collection of particles executing random walks on an undirected bipartite graph embedded in $\mathbb{R}^{D}$. In other words, at any integer time $t$, every particle is located on some node of the graph, and then, in the subsequent time step, it is located at some neighbor of that node. For the purposes of this paper, we shall assume henceforth that unless otherwise specified, any edge of a node where a particle is located is equally likely to be chosen for the subsequent step (and note that there may be more than one edge connecting any two nodes). We shall say that particles operating under these simple dynamics are executing Brownian motion. (We shall also implicitly assume all particles are initially on the same subgraph of the bipartite graph, and on the other one in the next, and so on.) The random walks of such particles are said to be \textit{ “heat-like”} if, in the continuum limit, the expected density $\rho(\textbf{x}, t)$ of the particles obeys the heat equation \cite[Chapter~1.3]{GFL10}






\begin{equation}
  \frac{\partial}{\partial t} \rho = k \nabla^{2} \rho
  \label{eqn:ContinuousHeat}
\end{equation}


for some positive constant k. Equivalently, a distribution is heat-like if the expected location of a particle that is located at $x_{0}$ at time  $t_{0}$  has, at some later time,  $t_{0}+t$,  the following Gaussian distribution centered on $x_{0}$ with a variance proportional to $t$:


\begin{equation}
\begin{split}
K(t, x, x_{0}) = \frac{1}{(4\pi k t)^{D/2}} e^{-\frac{r^{2}}{4\pi kt}} \\
\textrm{ where } r^{2} = \|x - x_{0}\|^{2} \textrm{ and } (t>0) .
\label{eqn:HeatKernel}
\end{split}
\end{equation}
Throughout this discussion, we shall restrict ourselves to solutions of the heat equation whose integral over $\mathbb{R}^{D}$ is conserved, except where otherwise noted. Assume next that these particles are endowed with an extra binary degree of freedom so that they may be said to be either “positive” or “negative” (or equivalently, can exist either as “particles” or “antiparticles”). The particles and their associated walks are “wave-like” if the expected \textit{net density} $\tilde{\rho}(\textbf{x}, t)$ (\ie, density differential of positive and negative particles) satisfies  the classical wave equation

\begin{equation}
  \frac{\partial^{2}}{\partial t^{2}} \tilde\rho = c^{2} \nabla^{2} \tilde\rho .
  \label{eqn:ContinuousWave}
\end{equation}


We shall show that any set of heat-like paths on a stationary graph becomes a wave-like set of (two-state) particles and paths if, at every node to which a particle travels, it spawns a particle-antiparticle pair (\ie, a pair of particles in opposite states) and if, furthermore, the new particle that has the opposite state of the original incoming particle always initially steps in the reverse direction of the original incoming particle (Figure ~\ref{fig:BrownianHuygensDiagram}). 

\begin{figure}[h]
\centering
\includegraphics[width=0.5\textwidth]{Fig1.tiff}
\caption{Brownian-Huygens propagation: any incoming particle (or antiparticle) spawns a particle-antiparticle pair, and the newly-created particle that is in the opposite state of the incoming particle must initially travel in the opposite direction of the incoming particle; the other outgoing particles randomly and independently choose any edge of the graph for their subsequent step; all outgoing particles become incoming particles in the subsequent time step.}
\label{fig:BrownianHuygensDiagram}
\end{figure}


In other words, for \textit{any} random lattice or graph on which the endpoints of random walks have a probability distribution that is a Gaussian centered on the starting point of that walk with a variance proportional to time (and given the Central Limit Theorem, the variety of walks which satisfy that requirement is wide), there is a related set of “mostly” random walks of particles (and antiparticles) whose distribution satisfies the classical wave equation. (We say “mostly” only to account for the sign- and direction-changing backpropagation steps whose direction is therefore not random). As long as the earlier Brownian motion particles and their implied paths are “heat-like”, these modified particles and their paths are assuredly “wave-like”. That being the case, collections of particles that exhibit wave-like properties are no more paradoxical than, say, particles that exhibit diffusion-like properties. All that is needed to transform the latter into the former is the ability for particles to exist in two states and stability of the graph or lattice to the extent that an edge that was traversed in one time step may be traversed in the reverse direction in the following time step.

Note that it is only for ease of visualization that we say the incoming particle has spawned an anti-particle pair. We could say instead that it is the incoming particle that always changes state (\ie, that it reversed its sign) and also reverses its direction while also spawning two other particles with the opposite sign of the one it now possesses. We prefer the latter memory-aid only because physicists will likely regard that as more familiar and evocative. In either case, the net particle count is strictly conserved. 

The first part of the proof will show that the above-defined particle-antiparticle dynamics do indeed lead to wave-like distributions in the specific case of regular lattices $\mathbb{Z}^D$ (for which the associated regular random walks are indeed heat-like).

Thereafter, we shall show that the normalized net count of the branching random walks in question can be expressed solely as a linear combination of values of the kernel of the heat equation – in other words, the expression does not depend on the particulars of the underlying graph on which the random walks take place (as long as the starting walks are heat-like). Therefore, what was shown to hold for walks on a regular $D$-dimensional lattice must be true for \textit{any} graph configuration on which random walks of a given number of $T$ steps have endpoints that are Gaussianly distributed with variances proportional to $T$. 

In other words, classical wave motion is compatible with a surprisingly wide variety of structures and media.

\section{Regular Lattices in D dimensions}


We shall first consider Brownian particles walking randomly along points of $\mathbb{Z}^{D}$ so that a step is always of length 1. Random walks that always traverse a unit distance along any axis have the property that the sums of the components at two successive times alternate in parity (and in this case, the parity of a point’s components distinguishes the two parts of the bipartite graph). For convenience, we shall restrict ourselves to the subgraph consisting of nodes whose integer components have a sum that is even whenever the time $t$ is even.

The transition matrix that gives the conditional probability that a particle at some point $X$ will in the subsequent time be found at some nearest neighbor of $X$, call it “HeatMatrix”, is an $N\times N$ matrix of the form
%%\usepackage{amsmath} % for 'bmatrix' environment
%%\begin{document}
\begin{equation}
\textrm{HeatMatrix} = \frac{1}{N} \times
 \begin{bmatrix} 
    1 & 1 & \dots \\
    \vdots & \ddots & \\
    1 &        & 1
    \end{bmatrix}
\end{equation}
%%\end{document}
where $N$ is the number of neighbors at a given node (\ie, the degree). For the regular lattice in $\mathbb{Z}^{D}$ we are considering, the degree of any node is $2D$, but so as to keep our notation general, we shall continue to use $N$ to signify the number of nearest neighbors at a node.

Note that HeatMatrix is idempotent and that its rows and columns sum to 1. As noted earlier, we shall refer to random walks whose expected values satisfy the above matrix as Brownian motion, regardless of the lattice in which they propagate.

We can also consider \textit{branching} random walks \cite[Chapter~1.1]{SHI12} for which any incoming particle produces a number of “descendant” particles at any step. For example, a particle that always produces another particle upon arriving at any node will have a related transition matrix, $\textrm{HeatMatrix}_{2}$ that is twice the value of HeatMatrix.

\begin{equation}
\textrm{HeatMatrix}_{2} = 2 \times \textrm{HeatMatrix}
\label{eqn:HeatMatrix2}
\end{equation} 

In that case, the number of particles will of course double at each time. We shall refer to this type of dynamics as Double Branching Brownian motion regardless of the graph in which it takes place. Note that if we choose to count paths, we will typically assume that a particle that was created at some time inherits the past locations of the particle that spawned it, so that every particle’s path, regardless of when that particle was created, has a history with the same number of steps.

We can also choose to give every particle a degree of freedom so that it can be either “positive” or “negative”, in which case, particles in the negative state will be referred to as “antiparticles”. We can likewise define net particle counts on any set of points as the difference between the count of positive and negative particles, and likewise we can define the net particle density. We can even choose to have particles and antiparticles traveling along the same edge at any time annihilate each other, though we will forego doing this whenever we wish to consider the set of all possible paths a particle might take.

Consider also the trivial dynamics wherein every particle always reverses direction and sign at every step. We shall refer to that dynamics as Huygens (back)propagation. In that case, the transition matrix that gives the expected net count of particles on an outgoing edge of some node in terms of the expected particle count along some “incoming” edge of that same node is
\begin{equation}
\textrm{HuygensMatrix} = -I
\label{eqn:Huygens}
\end{equation}

where $I$ is the $N \times N$ identity matrix.

To reiterate, due to the systolic nature of random walks on bipartite graphs where a particle is never included in its set of neighbors (in other words, cycles of length 1 are forbidden), the diagonal matrix elements in the above equation relate motions that are oppositely directed (and, due to the negation, of opposite sign) even though they refer to the same edge – in other words, they are a dynamics exclusively of back-and-forth path (and sign) reversal.

If we define the “momentum” of any particle at any time (on any arbitrary graph embedded in $\mathbb{R}^D$) as a unit vector in the direction of the propagation of a particle traversing across an edge in a single time step, multiplied by the speed required to traverse an edge along with some unit mass (\ie, ignoring any relativistic corrections merely for the sake of simplicity), we see that the momentum is exactly conserved for particles exhibiting Huygens behavior, whereas in the case of Brownian motion, or the doubled branching thereof, the expected value of the momentum is zero on any graph where the expected vector distance (as determined from a graph’s embedding in $\mathbb{R}^D$) from a node to its nearest neighbors has zero magnitude.

And as noted earlier, we can endow the particles with a sign and then calculate the expected net particle count (or net particle density). Using such two-state particles, we can consider a dynamical process that combines Double Branching Brownian motion and Huygens backpropagation, so that any incoming particle at a given node (being in either a positive or negative incoming state) results in an outgoing configuration of two other particles in the same incoming state (in other words a branching random process) and also an additional particle having the opposite sign, and that latter particle must travel in the opposite direction of the incoming particle. Each of the three outgoing particles will in the subsequent step likewise lead to three other outgoing particles in the same manner. (Note, however, that if we allow oppositely signed particles traveling along the same edge to annihilate, then it is possible for one or both of the particles with the same sign as the incoming particle to choose the reverse direction as well, in which case there will only be one outgoing particle in that step).

Even though the outgoing particles possess no properties outside their sign, we will say that the incoming particle executing a normal random walk spawned a particle-antiparticle pair and the newly created particle with the opposite sign of the incoming particle is constrained to traverse along the same edge the incoming particle used (though in the opposite direction). 

We could have also described the dynamics by saying it was the incoming particle that changed state and reversed direction (in other words, followed Huygens backpropagation), in addition to creating two other particles traveling along randomly chosen edges to neighboring points and with both those particles being in the same state that the incoming particle was, and stipulated that both new particles will likewise subsequently follow the Huygens backpropagation rule. However, the former description involving the creation of a particle-antiparticle pair will likely be more evocative and familiar to those who are familiar with modern physics, so that is the convention we shall follow.

This particular kind of propagation – as previously noted, we shall refer to it as “Brownian-Huygens propagation”, regardless of the of graph in which the particles travel -- is of special interest, since it conserves the particle count and has numerous other interesting properties; for example, it can be shown that if the associated Brownian motion paths on the same graph are heat-like (so that their continuum limit distributions are rotationally invariant Gaussians) then the expected momentum for Brownian-Huygens propagation is also conserved. 

Note that as was the case with Double Branching Brownian motion, every particle we create through Brownian-Huygens propagation has, in terms of its ancestry, a well-defined “wave path” that extends all the way back to some initial configuration. However, there is additional structure to consider than was the case with Brownian motion, given that an ancestor might have traveled in the reverse direction of a previous step and changed sign. But purely in terms of a sequence of locations, that takes no account of whether a step was Huygens or Brownian, a given wave path has an obvious correspondence to a specific element in the set of all possible Brownian motion paths. We shall call that element the “Brownian shadow” of that wave path, and because of the multiplicity inherent in any branching process, many wave paths will have the same Brownian shadow.

The matrix that gives the expected net particle count at an outgoing edge, given the presence of an incoming particle along some other edge, which we shall call WaveMatrix, is as follows:

\begin{equation}
\textrm{WaveMatrix} = 
 \begin{bmatrix} 
    (\frac{2}{N} - 1) & \frac{2}{N} & \frac{2}{N} & \dots \\
    \frac{2}{N} & (\frac{2}{N} - 1)  & \frac{2}{N}& \dots   \\
    \frac{2}{N} & \frac{2}{N} & (\frac{2}{N} - 1)  & \\
    \vdots & \vdots & \ddots & \\
   \frac{2}{N} &   \frac{2}{N}   &   & (\frac{2}{N} - 1) 
    \end{bmatrix} .
\end{equation}


%%\begin{equation}
%%\textrm{WaveMatrix} = 
%% \begin{bmatrix} 
%%    (2/N - 1) & 2/N & 2/N & \dots \\
%%    2/N & (2/N - 1)  & 2/N& \dots   \\
%%    2/N & 2/N & (\frac{2}{N} - 1)  & \\
%%    \vdots & \vdots & \ddots & \\
%%   2/N &   2/N   &   & (2/N - 1) 
%%    \end{bmatrix} .
%%\end{equation}

Which may be expressed more compactly as
\begin{equation}
\textrm{WaveMatrix}  = - I + 2 \times \textrm{HeatMatrix} .
\end{equation}

Note that WaveMatrix is a relation between the \textit{edges} of a graph (so that it relates particle propagations on time \text{intervals} between two pairs of successive integer times, as opposed to a relationship among function values at nodes at two specific times), so that the resultant value of the expected number of particles at any point is just the sum of the incoming or outgoing expected particle counts along the edges (or “flows”) of that node. (We shall regard HeatMatrix in the same way, even though, as we shall see, the underlying difference equation is first-order, which means that a particle’s contribution to the overall amplitude at a given node can be determined without reference to the edge along which it arrived.)

As was the case with HeatMatrix, WaveMatrix has rows and columns whose sum is 1, so that the sum of flows along any incoming edges converging on some node is the same as the sum of the particle counts in the node’s outgoing flows. And since HeatMatrix is idempotent, it is easy to prove from the previous relation between WaveMatrix and HeatMatrix that the former is involutory. Therefore, for a collection of particles obeying Brownian-Huygens propagation, the sum of the squares of the expected incoming particle counts or flows along any incoming edges to a node is the same as the sum of the squares of the expected outgoing particle counts (we stress that this equality, as was the case with momentum, only applies to the squares of the \textit{expected} particle counts).

In the case of regular Brownian motion, since a particle has an equal likelihood of traversing to any of its neighbors (and likewise, the expected number of particles at any node is just the average of the expected particle count at its neighbors in the previous step), then considering again a regular lattice $\mathbb{Z}^D$, the expected particle count at any time is related to the previous neighboring particle counts according to the following difference equation

\begin{equation}
  \begin{split}
  \Psi(t+1) = & \frac{1}{N} (\Psi(x+1) + \Psi(x-1) + \\
    & \Psi(y+1) + \Psi(y-1) + \\
    &  \Psi(z+1) + \Psi(z-1) \dots) .
  \end{split}
  \label{eqn:Heat}
\end{equation}

%% $\textrm{HeatMatrix}_{2}$ obeys, naturally, a nearly identical partial difference equation, except that the RHS is double of what it is in the previous equation (though the similarity, given the exponential growth in the latter equation, is deceptive).

%% \begin{equation}
%%  \begin{split}
%%  \Psi(t+1) = & \frac{2}{N} (\Psi(x+1) + \Psi(x-1) + \\
%%    & \Psi(y+1) + \Psi(y-1) + \\
%%    &  \Psi(z+1) + \Psi(z-1) \dots) .
%%  \end{split}
%%  \label{eqn:Heat2}
%%\end{equation}

In all of the partial difference equations we present here, we will, for simplicity, usually specify only the coordinate of a function that deviates from $\textbf{x}$ and $t$, or else omit arguments altogether, \textit{e.g.}:
\begin{equation}
 \psi(x) = \psi(y) = \psi(t) = \psi = \psi(x,y,z,…,t) .
\end{equation}

The Huygens propagation also follows a trivial difference equation due to the fact that particles reverse direction at every subsequent time step:

\begin{equation}
  \Psi(t+1) = -\Psi(t-1)
  \label{eqn:Huygens}
\end{equation}  

In this case, the difference equation applies not just to the expected value of the particle count at t+1, but the actual particle count, since Huygens dynamics are deterministic. 

If we choose to visualize a particle as traveling along an edge from one node of an edge to another before being transformed into some other particle(s), we shall assume that the change-of-state (or particle creation) occurs just before the particle arrives at a neighbor node, and that will specify a particle’s contribution to any count taken along the nodes. (In this paper, however, we shall primarily be concerned with dynamics in which the net count is conserved, so that whether particles change state just before or just after arriving at a node is irrelevant, as long as the convention is consistent for any changes of state.)

Given that the RHS in the equations for \eqref{eqn:Heat} and \eqref{eqn:Huygens}  refer to different times, the partial difference equation in the case of Brownian-Huygens propagation is easily seen to be


\begin{equation}
  \begin{split}
  \Psi(t+1)  + \Psi(t-1) = & \frac{2}{N} (\Psi(x+1) + \Psi(x-1) + \\
    & \Psi(y+1) + \Psi(y-1) + \\
    &  \Psi(z+1) + \Psi(z-1) \dots) .
  \end{split}
  \label{eqn:BrownianHuygens}
\end{equation}


By subtracting the “central” term $2\Psi(x,y,z \dots t)$ from both the LHS and RHS of the above equation (and in the case of the latter, apportioning that central term equally to each pair of terms inside the parentheses), and taking the continuum limit, we see that equation \eqref{eqn:Heat} converges to the familiar heat equation \eqref{eqn:ContinuousHeat}, where $k = 1/N$, whereas equation \eqref{eqn:BrownianHuygens} converges to the $D$-dimensional wave equation \eqref{eqn:ContinuousWave} where $c^{2} = 2/N$ (or, in the case of regular lattices, $c^{2} = 1/D$). 

In particular, with the chosen scaling, the coefficients of HeatMatrix are thereby all seen to be equal to $k$ whereas the off-diagonal and diagonal elements of WaveMatrix are equal to $c^{2}$ and $(c^{2} - 1)$ respectively. And even though the particles in this space move one lattice length per time interval, the speed of any wave pulse is $\sqrt{2/N}$, just as with lattice gases \cite{MAR86}.

As noted earlier, since the heat equation is first-order, calculating the likelihood that a particle at some node will be at a specific neighboring node in the subsequent step requires no earlier knowledge of where the particle had previously been. This is not the case with Brownian-Huygens and the related second-order wave equation. In that case, we must know both the arrival point of a given particle, and also the edge from which the particle came in the earlier step, so that we know along which edge the Huygens backpropagation particle must travel.





Given the nature of \textrm{HeatMatrix}, any sequence of nodes such that each successive member of the sequence is a neighbor of the preceding member constitutes an acceptable path, and has the same probability of being traversed as does any other path of equal length. It is a well-known property of random walks on a square lattice \cite{LUD18} that for any two points $X_{0}$ and $X_{1}$, and any corresponding times  $T_{0}$  and  $T_{1}$  with  $T_{1} > T_{0}$ , the value $K(T_{1} - T_{0}, X_{1}, X_{0})$ of the kernel of the heat equation is equal (in the continuum limit) to the number of all paths of length $t = T_{1} - T_{0}$ , that begin at $X_{0}$ and lead to $X_{1}$ in $t$ steps (divided by the number of all paths of length $t$).

We shall refer to the set of all such acceptable paths as HeatPaths($T$), and for convenience we shall omit mention of $X_{0}$ and $X_{1}$. We can similarly define WavePaths($T$), though as mentioned before, in order to fully specify an element of the latter set, we must specify which of its steps were Huygens (and these can only happen for path steps that are reversals of the immediately previous step); moreover, as noted earlier, given that two Brownian particles are generated at every time step, then even after fully specifying which step of a path, if any, was a Huygens step, there will be multiple identical paths in WavePaths($T$) having the same Brownian shadow in HeatPaths($T$). Finally, in order to be able to precisely determine which of the initial steps of any element of WavePaths($T$) can be Huygens (\ie, if that first step, along with its earlier step, form a 2-cycle on the graph), we shall assume that all the paths in WavePaths($T$) arrived at $X_{0}$ from the same neighbor $X_{\textrm{nbr}}$, though that, too, will be suppressed for ease of notation.

\section{Expected Particle Counts}

The same issue of multiplicity would affect the set of all paths generated by particles obeying Double Branching Brownian motion. That being the case, a kernel for its equation, which we will denote as $K_{2}(T_{1} - T_{0}, X_{0}, X_{1})$, is simply equal the kernel of the heat equation times $2^{T}$. And given that this multiplier does not depend on the nature of the lattice, then \textit{any} heat-like set of random walks on a graph, not just those on a square or even necessarily on a regular lattice, that are transformed into Double Branching Brownian motion, will have probability distributions and associated kernels that obey that same simple relationship to the regular heat kernel.

The same generality, with regard to what kind of graph the particles are propagating on, also holds for Huygens propagation, and however pathological it becomes in the continuum limit, any density of particles will also persist faithfully at every other time step without regard to the particulars of the underlying lattice. 

That being the case, it should not be altogether surprising if the same was true of Brownian-Huygens propagation as well, though the situation is admittedly more complicated upon convolving the double branching and Huygens backpropagation. However, there are some subsets of WavePaths($T$) which can easily be related to Double Branching Brownian motion and Huygens propagation. In particular, consider all elements of WavePaths($T$) that are purely Huygens propagation for some number of steps $T^{*}$, and then afterwards purely Double Branching Brownian motion. For any such path, the expected particle count is just a product of the expected particle count for both sections of the path. (Note that if the number of Huygens steps is odd, then the Brownian section must begin at the nearest neighbor from which the particle reached $X_{0}$.) We shall call this the Compositing Property of Brownian-Huygens paths, and will make use of it later. And even for the subset of paths which begin with a series of strictly Brownian steps followed by a series of strictly Huygens direction- and sign-reversing steps, the expected particle count likewise does not depend on the particulars of the underlying graph (assuming that regular Brownian motion on the same path produces heat-like paths), though in that case, the endpoint of the portion of the path that is Brownian is described by a Gaussian probability distribution that is centered around the point at which the Double Branching Brownian motion began. 

If we could show that this property (\ie, of being independent of lattice particulars) could be shown to hold for all the paths of WavePaths($T$), then again, it would mean that $\textit{any}$ set of heat-like random walks that are transformed according to the rules for Brownian-Huygens produce particle distributions that obey the wave equation. 

%%Because this expected particle count is, in effect, a propagator for the wave equation, we can expect that it will require some reference to two points at two different initial times, and knowledge of how those points are situated will indeed depend on the lattice. Moreover, we would of course still need to know the particulars of the graph in order to assign boundary conditions at two successive times. But as long as the solutions that are being fit

%%We would of course still need to know the particulars of the graph in order to assign boundary conditions at two successive times, but we are explicitly assuming that those particulars vanish in the continuum limit so as to yield solutions of the heat equation.

One way of going about proving that, that will require no more than what is in the rest of this paragraph, would be to keep extending the just-mentioned Compositing Property, given that every path in WavePaths($T$) can be broken up into mutually exclusive subsets for which the paths are either purely Double Branching Brownian motion (along which their contribution to the expected particle count can be expressed in terms of the heat kernel), or else, purely Huygens processes. For example, if the path in question consisted of $m_{1}$ Brownian steps, followed by $m_{2}$ Huygens steps followed by $T_{1} - m_{1} - m_{2}$ Brownian steps, the overall expected particle count is likewise a product of three terms. And since the contribution to the overall expected particle count for each of these subsets can be expressed in terms of the heat kernel without regard to the specifics of the lattice, we could show that the same is true for any set of such three-fold motley paths with identical sublengths.

But we shall instead pursue a modified version of this approach that is more pedantic, but thereby also more explicit (and that does not rely, even implicitly, on complicated and non-intuitive combinatorics that must respect the fact that Huygens steps can only happen on path reversals, and also the many integrations over the beginning and ending of each Brownian segment of any path). 

%% [DEPRECATE] Since the behavior of any particle on the lattice is independent of any other particles arriving or departing from a given node (if we forego allowing annihilation of oppositely signed particles), we only need to prove paths are wave-like for a single particle traveling from $X_{0}$ at  $T_{0}$ (from some specific nearest neighboring node in the previous time step)  to the point $X_{1}$ at  $T_{1}$. 

%% [DEPRECATE] As noted previously, even when extending our dynamics to either branching random walks or Brownian-Huygens propagation, a particle created at a given node will always have, in terms of its ancestry, a well-defined inferred random walk in the set of random walks corresponding to regular Brownian motion that extends back to $X_{0}$ at  $T_{0}$ (call that reconstructed or implied path the Brownian shadow of the inferred walk, but note that in the presence of branching there will be many inferred random walks with the same Brownian shadow.

In the case of standard Brownian motion on any lattice, given that the choice of any subsequent nearest-neighbor destination at a node is completely random, the expected particle count at any point $X_{1}$ and  $T_{1}$,  given that a particle was at $X_{0}$ and some earlier  $T_{0}$, is proportional to the number of all paths of length  $T_{1}$ - $T_{0}$  leading from $X_{0}$ to $X_{1}$, and if the average degree density in a given region of a lattice increases, the node density must decrease so as to leave unchanged the “heat-like” density of endpoints of paths with that given length.  In the case of paths for Double Branching Brownian motion, the expected particle count at $X_{1}$ and  $T_{1}$  is, as noted earlier, simply $2^{T}$ times the expected particle count in the case of Brownian motion. Whereas, in the case of Brownian-Huygens propagation, then regardless of the graph that is being considered, the value of the normalization constant will be exactly the same as it would be in the case of Brownian motion, since the net particle count is always conserved at every step.

Admittedly, the expected value of a given path or set of paths is, in the case of Brownian-Huygens, not immediately apparent, given the sign-reversing Huygens steps that occur at every particle’s step. Let us therefore first consider another subset of Brownian-Huygens paths where the multiplicity is even easer to determine than was the case when discussing the Compositing Property, namely, the subset of paths in which none of the steps are sign-reversing “Huygens” steps. The paths in that case are purely Brownian (or, more precisely, the paths of Double Branching Brownian motion).

%%Let’s call this subset WavePaths($T$), and remember that it not only depends on the same elided arguments as HeatPaths($T$), it also depends on the specific edge that all the associated paths traveled along to get to $X_{0}$ at $t=T_{0}$, for reasons which will be made clear below. (However, since the heat equation is first-order in time, and the wave equation is second-order, it should already be evident that when attempting to express waves in terms of the heat equation, some directionality will need to be added.) 

Note that while Huygens steps can only happen on steps in which a particle reversed direction (relative to the state of the particle in the prior step) and also changed sign, there can also be direction reversals without sign changes that have the same multiplicity associated with them as any other step.

So for this special subset, the number of paths is the same as for those following the dynamics implied by HeatMatrix2, in other words $2^{T}$ times the number of Brownian shadows that terminate at $X_{1}$ at time  $T_{1}$ , so that the expected particle count is likewise $2^{T}$ times the kernel of the heat equation at $X_{1}$ at time  $T_{1}$. (The same is true for paths of length zero, where the expected particle count is trivially a delta function corresponding to the initial conditions.)

Let us continue with this approach by partitioning the set of all acceptable wave paths of length $T$ into $2^{T}$ subsets, with each subset corresponding to a string composed of B’s and H’s that indicate which steps of the path are either Huygens or Brownian. Given that every element of WavePaths($T$) implicitly includes the same “starter step” from a specific neighbor $X_{\textrm{nbr}}$ of the nominal first step, let us agree to preface any such identifying string with an S as a placeholder for this starter step immediately preceding $T_{0}$. For example, if T=5, the subset of paths having just one Huygens step in the second position after the starter step would correspond to the string SBHBBB.

In other words, we will partition the set of wave paths into subsets which have the same number of Huygens steps at exactly the same times, so that their multiplicity (and sign) factor will likewise be identical, \ie,

\begin{equation}
\textrm{MultSgn} = 2^{n_{B}}*(-1)^{n_{H}} ,
\label{eqn:multiplicityandsign}
\end{equation}


and where, in considering the string of B’s and H’s for the subset in question, $ n_{B}$ and $n_{H}$ are the respective counts of each character, so that their sum is $T$.


For convenience, let us also agree to initially defer counting any subsets corresponding to an odd number of initial Huygens steps until we are done with all the other subsets. As we shall see, the paths in this deferred category -- as was already mentioned in our discussion of the Compositing Property -- require special consideration. So again, if the string specifying a given subset of WavePaths($T$) begins (after the S starter step) with an odd number $q$ of H’s and then either terminates (if $q = T$) or else contains a B in the following position, it shall be excluded from our count until subsets corresponding to all other strings have been considered.

\subsection{Subset of paths containing a single Huygens step at some specific time}
With that in mind, consider the subset of paths of length $T > 0$, for which there is a single Huygens step somewhere in the string, say at $T^{*}$, and note this cannot be the initial $t=T_{0}$ step (since the subset of paths that contain a single initial Huygens step is among those we have just agreed to exclude until the end of our count; likewise, for the same reason, $T$ must implicitly be greater than 1 for this subset. Note that we are considering the subset whose Huygens step occurs at some specific time, not the larger subset that is the union of all such distinct single-Huygens-step subsets. So, instead of considering all paths having only a single Huygens step (apart from the ones just excluded), a suitable subset would be all paths that have a single Huygens step specifically at the third position, \ie, SBBH…, \ie, T*=3 (where we have assumed for convenience that $T >= 3$).

Note that for any of these subsets, if we consider the path obtained by removing the single step and also the preceding one (the latter cannot be Huygens since there was only in the path by assumption), then since the Huygens step must be a reversal, the result of that two-step removal from the path is another path that begins at $X_{0}$ and terminates at $X_{1}$, but is of length $T-2$, since we have just excised from the path a cycle of length 2. Furthermore, the path has no other Huygens steps, since we are assuming there was only one to begin with. Apart from the factor of -1 that comes from the presence of the Huygens step in the longer string, and a factor of two from the excised Brownian step, then in accordance with the dynamics of HeatMatrix2, the expected particle count associated with both sets is exactly the same, and in the continuum limit, it reduces to

\begin{equation}
\int K(T^{*}-T^{0}, x’, x_{0})* K(T^{1} - T^{*}, x_{1}, x’) dx’,
\end{equation}

(again note that for this example we are assuming that $T^{*} = 3$). Note also that regardless of how many neighbors a point has at any of these paths at time $t= T^{*}$, its subsequent step must reverse or “undo” it, and therefore the presence of the additional 2-cycle will not affect the expected number of particles. 

In accordance with our supposition that all the paths in this space are heat-like, the expected contribution to the net particle count of any single subset with a Huygens step at $T^{*}$ is $-K(T-2, $X_{1},  $X_{0}$) times the multiplicity-and-sign factor for this subset, which from \eqref{eqn:multiplicityandsign}  is 2^{n_{B}}*(-1)^{n_{H}}
\begin{equation}
\textrm{MultSgn} = -2^{T-1}
\end{equation}

Having computed this quantity for a given value of $T^{*}$, we could take the additional step of grouping all those subsets into one term using simple combinatorics, but our primary purpose, which, again, is merely to show that the contribution to any such subset can be expressed without reference to the particulars of the lattice, has already been satisfied, so for our purposes, that additional step is unnecessary.



\subsection{A specific 1-dimensional example}
To give a specific example, consider a solution of the one-dimensional wave equation whose boundary conditions at $t=-1$ and $t=0$ are zero except at the point -1 at $t = -1$ and 0 at $t= T_{0} =0$, where the value at both points is one; in other words, a delta-function traveling wave moving along the positive direction with a speed of 1. The solution to the heat equation for boundary conditions that are zero everywhere at t=0 except at 0, where the value of the solution is 1, is just Pascal’s triangle, up to a normalization constant, with the value at any time $T_{1}$  and node $X_{1}$ equal to the count of paths of length  $T_{1}$  whose final value is $X_{1}$, divided by $2^{T_{1}}$.

Therefore, if we were to generate 1-dimensional paths as part of a Monte Carlo simulation, then the fraction of runs for which the path endpoint at $T_{1}$ is $X_{1}$ will approach the heat equation at those endpoint coordinates. Likewise, if we consider Double Branching Random motion, the expected particle count in that case would be $2^{T_{1}}$ times that expected value.

Likewise, if we generate a very large number of simulations of Brownian-Huygens propagation starting with a particle that arrives at the origin at $t=0$, having in the prior step been $x=-1$, and consider not only the full path of each particle (that includes ancestral steps prior to its generation, so that all paths similarly start with a step to the origin at $t= T_{0}$ ), but also the string that tracks at which step of a path a particle is in a Huygen’s state and therefore changed sign, then the observed count of particles divided by the number of simulations will approach the wave equation. 

Observe that for $T_{1}=4$, the nonzero values of Pascal’s triangle are $(1, 4, 6, 4, 1)$ all divided by 16, where $x$ is respectively $(-4, -2, 0, 2, 4)$. Let us consider paths of length 4 (plus a starter step) for which the associated strings all have a single H at the step corresponding at $T^{*} = 3$, \ie, SBBHB, and set $X_{1}$ to -2. In that case, there are just two paths to consider, with the coordinates $(0, -1, -2, -1, -2)$ and $(0, -1, 0, -1, -2)$. The value of the heat kernel needed to calculate the contribution for this subset is two steps earlier, and for this value of $X_{1}$ it is $1/4$ whereas the multiplicity and sign factor is -8,  so that the resultant value for this subset is their product, -2. Later on, we shall calculate the rest of the subsets and verify that it gives us the correct value for the wave equation at $x=-2$ and $t=4$, which is zero.






\subsection{Subset of paths containing exactly two Huygens step at two specific times}

Consider next a lattice subset of paths that have only two Huygens steps; let us assume they are at two specific times, $i$ and $j$, with $j > i$. Again, for any subset under consideration, the $i$ and $j$ are fixed; the subset of paths that have only two Huygens steps anywhere along the path is the union of all the subsets we are considering here, but we are considering each of them separately. In order that the subset not be in the excluded category, then if $i=1$, it must be the case that $j=2$; otherwise, any other distinct values of $i$ and $j$ less than or equal to $T_{1}$  would be acceptable. If $j=i+1$, then we can excise the two successive Huygens steps and thereby obtain a path that again has a length of $T-2$ and has no Huygens steps. If $j > i+1$, then (since we excluded the case where $i=1$ and $j>2$, given the resultant odd number of initial Huygens steps) we can excise both Huygens steps, and also the Brownian step preceding each of the Huygens steps, thereby obtaining a path of length $T- 4$ that also has no Huygens steps. In either case, using the same logic that we applied for subsets corresponding to strings with a single Huygens step at some specific time, we are left with a term proportional to either $K(T-2,\dots)$ or $K(T-4,\dots)$ times some power of 2, neither of which have any dependence on the particulars of the graph in which these particles propagate. 

It is also easy to see that this category correctly captures, in the case where  $T_{1} =2$, the contribution of the subset consisting solely (apart from the presumed starter step) of two Huygens steps; since the initial particle arriving at the origin must, according to the rules of Brownian-Huygens propagation, create a Huygens particle, and the same is true for that Huygens particle in the next step, the contribution in that case is 1, which again, does not depend on the particulars of the graph. 

\subsection{Subsets with three or more Huygens steps}
Consider next a lattice subset that has exactly three Huygens steps – let us assume they are at fixed positions $i, j,$ and $k,$ with $i < j < k$. Again, we will exclude from consideration any subset whose initial step is a Huygens step followed by a Brownian step, or else three Huygens steps; otherwise, all other distinct values of $i$ and $j$ and $k$ less than or equal to  $T_{1}$  are acceptable. Topologically, we can consider three cases: 1) three Huygens steps, each separated by one or more Brownian steps; 2) two consecutive Huygens steps with another separated from them by one or more Brownian steps; or 3) three consecutive Huygens steps, where the initial step is Brownian, so as to avoid being in the excluded category. In the first case, we can excise the three Huygens steps along with the three preceding Brownian steps and arrive at a term of the form $K(T-6, \dots)$ times the multiplicity-and-sign factor, which, since $n_{H}$ is odd, is a negative power of 2 for any of these three-Huygens cases. In the second case, we can excise the lone Huygens step and preceding Brownian step, and also the two successive Huygens steps, and and arrive (after multiplying through by the multiplicity-and-sign factor) at a term proportional to $-K(T-4, \dots)$. In the final case, we will excise all three Huygens steps and the preceding Brownian step and arrive again at a term of the form $ -K(T-4, \dots)$ (note the power of two from the multiplicity-and-sign factor will differ for these last two cases, given that differing numbers of Brownian steps were excised). 

In any of these cases, we can again see that the contributions can be written without any reference to the underlying graph. 

Clearly, paths with larger numbers of Huygens steps (assuming the subsets in question are not in the excluded category) can be treated in the same manner. Any grouping of an even number of Huygens steps can be excised together, whereas a grouping of any odd number of Huygens steps (again, we assume they do not occur at the beginning of the path) will likewise be excised along with the preceding Brownian step.

Since each excision trims off at least two path steps, there can be at most $T_{1} /2$ of them. While this approach is wholly impractical as a computation algorithm, given that simply solving the numerical difference equation is much simpler, the very fact that such a decomposition is possible, and that is does not refer any specification of an underlying graph (beyond assuming that regular Brownian motion there is heat-like to the extent that the sum of paths can be computed by way of the heat kernel) is all we needed to establish. 


\subsection{Paths beginning with an odd number of Huygens steps}
Therefore, it remains only to consider those excluded cases for which the paths begin with an odd number $m$ of Huygens steps and then either terminate, or else continue with a Brownian step. Recall yet again that all the wave paths under consideration initially arrived at $X_{0}$ at $t = T_{0}$  from the same edge. Since any excision of a group of odd number of Huygens steps requires an earlier step to be excised as well, we can instead excise that previous starter step along with the initial odd Huygens steps, the dynamics of which is given by pure Huygens propagation, so that the expected particle count  (for that series of Huygens steps considered in isolation) is -1, and the destination is just the neighboring node $X_{\textrm{nbr}}$ from which the particle arrived at the origin. Therefore, the expected particle count of the remaining portion of the path (\ie, apart from the factor of -1) will be the same as for a path whose initial step likewise begins at $X_{\textrm{nbr}}$). Once the remaining portion of that path has been excised of all its Huygens steps, as was done with the earlier subsets, we can use the Compositing Property mentioned earlier to obtain the final value, and see that even for this excluded category of paths, their contribution to the expected particle count is (apart from the multiplicity-and-sign factor) is equal to the value of the propagator $K(t - m, X_{1}, X_{\textrm{nbr}})$. 

Therefore, our task (of showing that the relationship between heat-like walks and the associated expected value of particle counts in the case of Brownian-Huygens propagation does not depend on the particulars of the graph on which the motion takes place) has been accomplished, since the sum of the paths can be expressed as $2^{T}$ terms of the heat kernel at times and locations as specified above, without reference to any graph or lattice whatsoever. Moreover, we know from what was determined by way of regular lattices in $D$ dimensions, the sum of the associated Brownian-Huygens wave paths on that graph gives us a propagator for the wave equation -- in other words, the wave paths in that case are \textit{wave-like} (thereby validating our choice to call them wave paths in the first place). So it therefore follows that \textit{any} set of random walks on a lattice whose endpoints are described by Gaussian distributions with variances proportional to time (and given the Central Limit Theorem, the variety of such walks can be said to be large), can be transformed into wave-like paths of two-state particles whose endpoint densities obey the wave equation (assuming that motion along an edge in any step can be reversed in the subsequent step, as Brownian-Huygens propagation requires). Note that if the degree varies for the different nodes of a graph, then by assumptions the effective scaling terms in the heat equation satisfied by those heat-like paths is no longer simply $1/N$, but rather, some composite term related to the effective number of neighbors in the lattice, and in that case, the same can be said for the resultant speed of light in the case of the wave-like paths.
%%Therefore, our task -- of showing that the relationship between heat-like walks and the associated expected value of particle counts in the case of Brownian-Huygens propagation does not depend on the particulars of the graph on which the motion takes place -- has been accomplished, since the sum of the paths can be expressed in terms of the heat kernel (which makes no reference to any lattice). While we still have to know enough about the particulars of any lattice to be able to assign values to the underlying graph at two successive times (or, for the space of non-trivial solutions that interest us, to assign values along all the edges that connect nodes at two subsequent times), the assumption that the associated random walks are heat-like (to the extent that their expected endpoint locations are representable by the Gaussian kernels of the heat equation) is enough to likewise assure that the corresponding Brownian-Huygens evolution (which is likewise representable by Gaussian kernels) is independent of the particulars of the lattice.

%%Moreover, we know from what was determined by way of regular lattices in $D$ dimensions, the sum of the $2^{T}$ cases we have just enumerated gives us a propagator for the wave equation -- in other words, the wave paths in that case are \textit{wave-like} (thereby validating our choice to call them wave paths in the first place). Moreover, it follows that \textit{any} set of random heat-like random walks on a lattice can be transformed into wave-like paths of two-state particles whose endpoint densities obey the wave equation (assuming that motion along an edge in any step can be reversed in the subsequent step, as Brownian-Huygens propagation requires). Note that if the degree varies for the different nodes of a graph, then the effective scaling terms in the heat equation satisfied by those heat-like paths is no longer simply $1/N$, but rather, some composite term related to the effective number of neighbors in the lattice, and in that case, the same can be said for the resultant speed of light in the case of the wave-like paths.

\subsection{Return to the 1-dimensional example}
Now we can consider the heat kernel weights for the rest of the 1-dimensional paths considered in the earlier example we noted, wherein a solution of the wave equation at $x=-2$ and $t=4$ is calculated directly in terms of values of the heat kernel, as derived by way of Pascal’s triangle.

\begin{table}[]
 \centering
\begin{tabular}{llrrrc} \toprule
\multicolumn{6}{l}{Four-Step 1-D WavePaths($0\to -2$) (starter step at -1)} \\
\hspace{0.5cm} Path   & String & Wt & Krnl & $2^{n_{B}}$ & Krnl(x,t)\\ \toprule
(0,-1,-2,-3,-2) & SHBBH & 2  & $1/2$ & 4  & $(-1,1)$ \\
(0,-1,-2,-3,-2) & SHBBB & -3 & $3/8$ & 8  & $(-1,3)$ \\
(0,-1,-2,-3,-2) & SBBBH & -2 & $1/4$ & 8  & $(-2,2)$ \\
(0,-1,-2,-3,-2) & SBBBB & 4  & $1/4$ & 16 & $(-2,4)$ \\
(0,-1,-2,-1,-2) & SHBHH & -1 & $1/2$ & 2  & $(-1,1)$ \\
(0,-1,-2,-1,-2) & SHBHB & 2  & $1/2$ & 4  & $(-1,1)$ \\
(0,-1,-2,-1,-2) & SHBBH & -  & $1/2$ & 4  & $(-1,1)$ \\
(0,-1,-2,-1,-2) & SHBBB & -  & $3/8$ & 8  & $(-1,3)$ \\
(0,-1,-2,-1,-2) & SBBHH & 1  & $1/4$ & 4  & $(-2,2)$ \\
(0,-1,-2,-1,-2) & SBBHB & -2 & $1/4$ & 8  & $(-2,2)$ \\
(0,-1,-2,-1,-2) & SBBBH & -  & $1/4$ & 8  & $(-2,2)$ \\
(0,-1,-2,-1,-2) & SBBBB & -  & $1/4$ & 16 & $(-2,4)$ \\
(0,-1,0,-1,-2)  & SHHHB & -1 & $1/2$ & 2  & $(-1,1)$ \\
(0,-1,0,-1,-2)  & SHHBB & 1  & $1/4$ & 4  & $(-2,2)$ \\
(0,-1,0,-1,-2)  & SHBHB & -  & $1/2$ & 4  & $(-1,1)$ \\
(0,-1,0,-1,-2)  & SHBBB & -  & $3/8$ & 8  & $(-1,3)$ \\
(0,-1,0,-1,-2)  & SBHHB & 1  & $1/4$ & 4  & $(-2,2)$ \\
(0,-1,0,-1,-2)  & SBHBB & -2 & $1/4$ & 8  & $(-2,2)$ \\
(0,-1,0,-1,-2)  & SBBHB & -  & $1/4$ & 8  & $(-2,2)$ \\
(0,-1,0,-1,-2)  & SBBBB & -  & $1/4$ & 16 & $(-2,4)$ \\
(0,1,0,-1,-2)   & SBHBB & -  & $1/4$ & 8  & $(-2,2)$ \\
(0,1,0,-1,-2)   & SBBBB & -  & $1/4$ & 16 & $(-2,4)$ \\ \toprule
 & \textbf{Total:} & 0 & & &\\ \toprule
\end{tabular}
\caption{\label{tab:TablePascalTriangle}Calculating the 1-dimensional wave equation at $x=-2$ and $t=4$ by tabulating all paths of length 4 (plus an implicit starter step to the origin beginning at x=-1) that begin at zero and terminate at $x=-2$. Only paths with a non-zero weight are shown, and weights for a given subset of paths (characterized by the indicated string) are displayed only for the first listed path of the subset; otherwise, they are suppressed. Since the 1-dimensional positive-direction traveling wave specified here has only a single non-zero value of 1 (at $x=4, t=4$), the sum corresponding to the point $x=-2$ is zero. The remaining columns are explained in the text.}
\end{table}



Table \eqref{tab:TablePascalTriangle} displays the 22 distinct paths (and 12 distinct strings) leading from 0 to -2 in four steps that produce non-zero contributions, \ie, weights, to the solution of the wave equation. The remaining 4 possible strings (for a total of $2^{4}=16$) have zero weights and so have been ignored; for example, the string SHHBH, after excising both the initial pair of H’s as well as the remaining BH string, is just the empty string (apart from the starter step), and that means that the relevant value of the heat kernel is at $t=0$ and $x=-2$, where it is zero. Some strings (\textit{e.g.}, SHBBH) have multiple paths, and though these are both displayed, the weight is specified only for the first instance (and accounts for the entire set of paths corresponding to that string) and is left blank for any subsequent listings of the same string. In each case, the magnitude of the weight (apart from the factor of $-1$ to the power of the number of Huygens steps) is simply of the product of the heat kernel at some (possibly earlier) coordinates (those $x$ and $t$ coordinates are listed in the final column) times two to the power of the number of Brownian steps in the corresponding path. 

As can be seen by inspection, the sum of the weights in this case is indeed zero at $x=-2$, as it should be.

Since we have assumed a starter step that arrives at the origin from the point $x=-1$, a similar tabulation for all space and time produces non-zero values only for the values $x=t$, where the value of the specified traveling wave is 1. For example, at the wave-front value $x=t=4$, only the SBBBB string produces a non-zero value, and the value there is simply the multiplicity-and-sign factor, $2^{4}$. times the value of Pascal’s triangle there, \ie, 1, times the overall normalization constant which is the inverse of the first term, for a result of 1. (Any other strings would, according to the formalism laid out here, involve values of the heat kernel at some earlier space-time coordinates that lie outside the reach of any paths with duration less than 4 – \ie, at earlier times the destination is outside the light cone  -- and so their weight must be zero.)

This can be directly applied to higher dimensions as well (including, for example, hexagonal or triangular regular lattices in two dimensions, arguably the easiest configurations to consider after the two-dimensional square lattice). As noted in the earlier section, \textit{any} heat-like walks can be modified in such a way, though in cases of irregular lattices, closed form solutions of the heat equation (apart from the continuum limit case of a Gaussian) will be harder to obtain.

\section{Additional degrees of freedom, anisotropy, and mass}

\begin{figure}[h]
\centering
\includegraphics[width=0.5\textwidth]{Fig2.tiff}
\caption{(L) A graph with one edge connecting each neighbor along which oppositely signed particles annihilate. (R) A graph with three edges connecting any nearest neighbor, suitable for particles with three binary degrees of freedom.}
\label{fig:MultiplyConnectedGraphs}
\end{figure}



Note that, as mentioned before, there can be multiple edges connecting two nodes. For example, consider a regular $D$-dimensional lattice where the number of edges connecting any two nearest neighbors is some integer $m > 1$ (Figure ~\ref{fig:MultiplyConnectedGraphs}). In that case, HeatMatrix and WaveMatrix have the same general form as before, but $N$, and therefore the size of the matrices, is now multiplied by $m$; otherwise, the idempotency and involutority remain unchanged. Such multi-edge variations are useful if we desire our particles to have additional states – for example, charge. And if we choose, for any reason, to impose a “Pauli exclusion principle” on the particles – in other words, mandate that only particle in a given state can travel along an edge at any time -- then if we apply that exclusion only to particles that are in identical states, such multi-channel graphs can ensure that annihilation or exclusion only involves particles that have identical states.

%%[DEPRECATE] Consider also that there are rare fluctuations in our dynamical rules and that in some rare instances, a particle chooses in rare instances, with a probability proportional to $q\epsilon^{2}$ (with the epsilon being the scale factor that tends to zero in taking the continuum limit), to engage in two reversal steps without any sign change, and without spawning any additional particles. (Equivalently, we can say that all three particles in that case must back-propagate along the same edge that the Huygens step does, for a net count of one.) Perhaps this occurs as a result of some inhomogeneity in the graph or lattice, or because some other kind of particle or event in an edge that allows the Brownian-Huygens motion to be altered at both nodes in that edge. Such “lapses” in Brownian-Huygens propagation still obviously conserve the particle count, and while they temporarily disrupt the momentum, as long as the particle continues after those two steps to return to the usual Brownian-Huygens propagation (with the same notion of previous direction having been observed) then these back-forth directionally-antisymmetric events amount to nothing more than another episode momentary tunneling (and note that even with regular Brownian-Huygens, momentum is only conserved on average, and not at every single step). But to the extent that a collection of particles is proceeding along a wave-front, this lapse amounts to a pause, and dampens the progress of that wavefront, so that it is intuitively a good candidate for modifying the wave equation to account for wavefunctions appropriate for quantum mechanical particles with a non-zero mass.

 %%As with the earlier example of excising Huygens steps from the summation of wave paths that go into a wave amplitude, these two-step lapses can also be excised, and if we expand their effect on the wave amplitude (expressed as a sum over paths) in powers of $q\epsilon^{2}$  for only a few terms (given that the resultant probability is small) we are left with a summation of corrections of the form
%%\begin{equation}
%% (q\epsilon^{2})^{j}  \textrm{WavePaths}(t-2j)  .
%%  \label{eqn:DEPRECATEDKleinGordon}
%%\end{equation}
%%(Had we not chosen \textit{two} successive steps for this lapse, then we would have to excise a preceding Brownian step along with the actual “lapse” step, which means we would also have to treat any single-initial lapse at the beginning of a path in separate fashion, but since the successive lapse steps are assumed to occur in pairs, we can just excise the pair of them together from the path in the same way that was done for pairs of successive Huygens steps.)
%%If the probability of such an event is sufficiently small so that everything beyond the first term can be ignored, the difference equation for this modified Brownian-Huygens becomes

\begin{figure}[h]
\centering
\includegraphics[width=0.5\textwidth]{Fig3.tiff}
\caption{A case where a graph’s anisotropy can be compensated with a scale transformation (in this case, stretching along the $y$ axis or compressing along the $x$ axis), thereby making the associated random walks heat-like.}
\label{fig:AnisotropicGraphs}
\end{figure}


Note that the formalism developed here can also be applied to anisotropic graphs and lattices that require a scale transformation in one or more dimensions in order for the random walks on the graph to become heat-like (Figure ~\ref{fig:AnisotropicGraphs}). Also, the discreteness, and the manner in which the density of particles scales in one dimension or another can be used to our advantage.


\begin{figure}[h]
\centering
\includegraphics[width=0.5\textwidth]{Fig4.tiff}
\caption{Adding a “small dimension”, isomorphic to $\mathbb{Z}/2\mathbb{Z}$, along the $z$ axis, by doubling a two-dimensional $xy$ graph. Note that not all the nodes on the graph are required to possess edges leading along the added dimension.}
\label{fig:AddedDimension}
\end{figure}

In particular, consider a graph where a random subset of particles (such that their density is constant throughout the space) has an added dimension. For convenience, we may assume the new dimension is minimally small, \ie, isomorphic to $\mathbb{Z}/2\mathbb{Z}$). To visualize this in two dimensions, imagine that a random graph in the $xy$ plane is duplicated and then stacked along the $z$ axis, and some subset of points (for now, assume they have a uniform density) have pairs of edges leading up and down the $z$ axis to a node’s stacked duplicate (Figure ~\ref{fig:AddedDimension}).  The evenness of the extra dimension ensures that the graph’s two parts can still be distinguished simply according to the parity of the components of this added dimension (and the same holds true for any regular graphs that are toroidal in one or more dimensions, so that we shall stipulate that for any such $\mathbb{Z}/m\mathbb{Z}$ dimension, that $m$ must likewise be even). In any graph where the edges do not lie along some axis of $\mathbb{R}^{D}$, say a hexagonal two-dimensional lattice, or some irregular $D$-dimensional graph, the correspondence between a given particle motion and an axis is not well-defined as it is in the case of the regular $D$-dimensional lattices we began with. But the fact that this added dimension, by construction, has a distinct axis and a uniform size means that movement in an edge along that dimension -- let us refer to it henceforth as a “mass dimension” -- is distinct from movement along any of the remaining “spatial” edges. Moreover, there is now an eigenvector in our space of functions corresponding to phase rotation along the axis of that new dimension. Let us isolate that mode and impose it onto the dynamics itself, by demanding that any particle that has taken a Brownian-Huygens propagation step along a spatial edge that then takes a step on an edge along the mass dimension (which means it appears to remains stationary with regard to the remaining dimensions, in the same way that a particle in a regular lattice is stationary with respect to any set of axes other than the one it is propagating along) incurs a factor of -1, just as with a Huygens backpropagation step, and likewise, another sign change when it chooses to take a step along a spatial edge after having traversed some non-zero number of steps along the mass dimension – in other words, any transition into or out of motion along the mass dimension involves a change of sign. We shall call this modified Brownian-Huygens propagation. The length of a step along this mass dimension is also arbitrary, apart from inducing a scale factor $\alpha$ corresponding to that dimension length. (We shall further assume that the length is uniform throughout the space.)
%%Therefore, the node locations on the original graph $\textbf{x}$ become $(\textbf{x}, i)$.


Note that on this added dimension, the $2^{\textrm{nd}}$-order partial derivative of any function $f,$ \ie,

\begin{equation}
\frac{f(x, i+1, t) - 2 f(x, i, t) + f(x, i-1, t)}{\epsilon^{2}} ,
 \label{eqn:dscretespatial2ndderiv}
\end{equation}

where $\epsilon$ approaches zero in the continuum limit, has a related function differing only in the sign of the second summand of the numerator and a normalization factor of $1/4$:
\begin{equation}
\frac{f(x, i+1, t) + 2 f(x, i, t) + f(x, i-1, t)}{4\epsilon^{2}} ,
 \label{eqn:dscretetemporal2ndderiv}
\end{equation}

which we shall call, for obvious reasons, a $2^{\textrm{nd}}$-order partial \textit{smoothing} of $f$. Even though $f(x, i+1, t)$ and $f(x, i-1, t)$ are the same point (since the added dimension is isomorphic to $\mathbb{Z}/2\mathbb{Z}$), expressing them in more familiar form will make what follows easier. We shall also use the flexibility afforded to us by our discrete space to consider for the duration of this section only those expected particle counts that are so small in magnitude that the $2^{\textrm{nd}}$-order partial smoothing remains bounded even when $\epsilon$ is close enough to zero for our space to be indistinguishable from the continuum. We shall further limit ourselves to functions that are already smooth enough so that the above smoothing is arbitrarily close to the continuum value of the function at $x$ (at which point $i$ is assumed to be either 0 or 1 so as to properly specify the non-zero part of the bipartite graph). 

Assuming that Brownian walks in this expanded graph is (up to the aforementioned scale transformation along the mass dimension) are heat-like, then Brownian-Huygens propagation in the graph will create walks that are wave-like. As for the modified Brownian-Huygens propagation with the added sign change incurred when switching between a step along some spatial edge and one that is aligned along the mass dimension, and vice versa, we have the following related equation

\begin{equation}
\frac{\partial^2}{\partial t^2} \Psi = \frac{2 }{N’}\left(\nabla^{2}_{N’} \Psi - \frac{\alpha_{c}}{4} \left(\Psi(i+1) +2\Psi  + \Psi(i-1)\right)\right) ,
 \label{eqn:preKleinGordon}
\end{equation}


where $\Psi(i \pm 1)$ signifies the nearest neighbor locations along the newly added dimension (again, given the isomorphism to $\mathbb{Z}/2\mathbb{Z}$, both refer to the same node) and $N’=N+2$, and $\alpha_{c}$ is the resultant value of the scale factor $\alpha$ in the quasi-continuum limit, and where the discrete Laplacian operator  $\nabla^{2}_{N’}$ conforms to the number of neighbors $N’$ of the node in question (and where the discreteness of the $2^{\textrm{nd}}$-order time derivative is understood). Given the smoothing operator we have chosen, the central term $\Psi(\textbf{x},i,t)$ in all the respective $2^{\textrm{nd}}$-order derivatives or smoothings in \eqref{eqn:dscretespatial2ndderiv} and \eqref{eqn:dscretetemporal2ndderiv}  still cancels on both sides of the above equation (though $2/N’$ is of course smaller in magnitude than $2/N$). Therefore, in the near-continuum limit we are considering, the above equation reduces to the familiar Klein-Gordon equation with $\alpha_{c}$ proportional to the square of the associated mass, justifying our usage of “mass dimension”.
%%The added dimension can therefore be referred to as a “mass dimension”.

%%We can define a branching process similar to Brownian-Huygens propagation in that any movement in the original graph that does not change the position along the added dimension is as before (in other words, a particle in that case must still shift from one node of the original graph to a neighbor); however, any particle movement along the new dimension (for which the original coordinates stay the same) incurs a sign shift. 

If we consider the previous graph’s WaveMatrix to have been expanded by two columns and rows (representing forward and reverse step-wise movement along the added “mass dimension”) the coefficients there are negated relative to the rest of the matrix, except for the $2 \times 2$ submatrix at the very bottom of the diagonal, where the coefficients retain the same sign. (All of the coefficients, however will change in the sense that $N$ has now increased to $N+2$.) We shall call this modified matrix KGMatrix.

Although KGMatrix is still involutory, so that the sum of the expected squares of the amplitude in the edges is conserved under this kind of propagation, the sum of the components along any row or column is no longer unity as before, so that the sum of the amplitudes is no longer strictly conserved. However, given that a particle always incurs a change of sign with every transition from a step along the mass dimension to a step along some spatial edge (or vice versa), it is easy to see, even in the case where the length of the mass dimension is some larger even number than 2, that the sum of the amplitudes in the even-indexed slices along the mass dimension axis minus the sum of the amplitudes in the odd-indexed slices is the same absolutely conserved quantity that would be the sum of the amplitudes if the particle were following the usual Brownian-Huygens propagation rules we defined earlier, instead of the modified version. So in that sense, the conservation of amplitudes still holds.

If we consider any paths in this expanded space, we can easily excise any steps taken along the mass dimension (which will always occur in even numbered groups) by repeating the approach described in the previous section (that was used to progressively excise Huygens steps from a set of paths) and use that to expand any values of the wave functions for a massive particle into either a sum of terms consisting of factors of the regular wave equation, or else, a sum of terms of the propagator of the heat equation (with the excised walks taking place on a graph without any added mass dimension).

Depending on the equation we are interested in modeling, we could introduce still other kinds of alternate small dimensions (with possibly longer lengths), and specify more complicated phase factor transitions, rather than the single mass dimension and the simple negation that will suffice for the purpose of this paper, and in that way account for more complicated particle symmetries. In other words, we did not need to actually double and stack the graph – the extra dimensionality is just a way to incorporate an internal structure to the particles. For that matter, we did not need to introduce signs or antiparticles, and could have continued with scalar particles, and instead, accommodated that two-state degree of freedom by similarly modifying the graph.

The approach here is in some ways a discrete and very rudimentary analog of string theory, and the manner in which dimension is treated there, but the goal of this section is simply to demonstrate that the approach we have chosen is not limited to the wave equation, and that discreteness, far from being a disadvantage, allows for a flexibility that is peculiarly well-suited for the path-integral equations of modern physics.
%%, to the extent that the usual approaches, that dictate strict conformity to continuous analytical methods, seem Procrustean in comparison.


\begin{figure}[h]
\centering
\includegraphics[width=0.5\textwidth]{Fig5.tiff}
\caption{Two graphs, each consisting of an irregular 1-dimensional graph along with a mass dimension isomorphic to $\mathbb{Z}/2\mathbb{Z}$. Due to the differences in the density of edges along the mass dimension, the associated mass parameter of the associated Klein-Gordon differs between the two graphs. (Note also that the left-hand side of the top graph has a different density of such edges, so that even within that graph, assuming the density differentials of such edges within different regions persists even in the continuum limit, the mass would then vary across those differing regions.)}
\label{fig:KGgraph}
\end{figure}


Moreover, the discussion on anisotropic graphs applies to the mass dimension as well. It is not necessary that every node of the graph contain an edge leading along the mass dimension. A graph that is “sparse” in such nodes would require a scale transformation, relative to one that is not, in order for the generated waves to have identical densities and properly Gaussian evolution distributions, and in this case, that would mean changing the magnitude of the mass term (Figure ~\ref{fig:KGgraph}). Moreover, if the density of such mass dimension edges were different in one region of the space than another, the mass parameter would likewise vary, so that this formalism easily accommodates the curved spacetimes of general relativity in which mass varies. If we admit exogenous forces and interactions (and any such deviation from our rules that that entails), then any edges along the mass dimension can be conditioned on some exogenous variable or coupled to the presence of another kind of particle in the space; and given that all this can be accomplished with even a random arrangement of nodes and edges, spacetime can expand or shrink more easily than would be the case if the ability to model wave equations required the pristine symmetry of a regular lattice.

%%Finally, we shall leave it to the reader to verify that in some inertial reference frame along, say, the $x$ axis, the length along that axis will appear to shrink (so that the density of these mass dimension edges, and therefore the mass term, will increase) by a factor of $\gamma$, the familiar Lorentz factor of special relativity \cite[p.~63]{FEY63}, as might be expected from any dynamics that, in the continuum limit, models a relativistically invariant equation.

Therefore, this Brownian-Huygens approach is suitable not just the classical wave equation, but also quantum mechanical wavefunctions for particles with a mass. Since the transfer matrix approach already suffices to transform the second-order wave equations into first-order (a discrete analog of Hamilton’s transformation of second-order Lagrangian mechanics into a system of first-order equations), there is no need to further pursue Dirac’s approach of “taking the square root of the d’Alembertian operator”. 

\section{Quasi-probability particle densities}
In the rest of the paper, we shall for simplicity’s sake ignore mass and return to simple Brownian-Huygens propagation. In addition to endowing particles with a sign, let us expand the universe to two distinct types of particles propagating on the same graph, whose randomly chosen path steps are completely independent. (Note that even with this two-fold expansion of particles, the added degrees of freedom still amount to far less complexity than is contained in, say, a complex-valued 4-vector solution of the Dirac equation). Following Dirac’s nomenclature, we shall refer to this expanded universe of particles as being either “bra” or “ket”, respectively. Let us assume that the net density of both types of particles at $t=0$ is equal to some function $\Psi(\textbf{x}, 0)$, for which the integral of the square of the function is unity. Then, as has previously been noted \cite{HRG92}, at any subsequent time $t’$, the product at some point $\textbf{x}$ of the net particle count of the bra particles times the net particle count of the ket particles has a positive expected value that is everywhere equal to the square of the expected value of the wave equation $\Psi(\textbf{x}, t’)$. (Recall that since WaveMatrix was previously shown to be involutory, then the expected sum of the squares propagating along the edges of the graph is conserved, as is the sum.) By the law of large numbers, for a sufficiently fine-grained graph, and a sufficiently large region of that graph, the expected value of this count can be made so that the likelihood of it being negative when summed over this region is arbitrarily small, even if this region is itself a very small region of the entire graph. Admittedly, since Brownian-Huygens is a branching process, so that there is no single particle trajectory to speak of, the resultant wave dynamics here is more closely related to that of a quantum field than a wave-function for a single particle, but given the latter is only a non-relativistic limiting case, that is more a problem for the notion of single-particle wavefunctions than this model. 

Such generalized densities have similarities to quasi-probability distributions by Wigner and Moyal \cite{WIG32} \cite{MOY49} \cite{DIR42}. Given the fact that they can be negative, such quasi-probabilities do not violate Bell’s theorem (\cite{FEY87}, despite their being local, and the same holds true of this approach.

Let us define a measurement as the operation of integrating any such bra-ket product count over some space-time path, so as to include whatever is on (or sufficiently close to) that path. Note that the presence of a negative bra-ket summand reduces the count, and different sequences of measurements along different paths will lead to different perceptions (or “worlds”),  similar in some respects to how the distributions of electrons on an antenna can yield distinctly different programs depending on the frequency with which it is sampled; \cf~the Wigner-von Neumann interpretation of quantum mechanics (\cite{von32}, \cite{WIG61}) in which the subjective perception of an observer depends on, and is intrinsic to, a quantum mechanical measurement. Assuming all observers in such a space are of a sufficiently large macroscopic size, and only capable of registering macroscopic counts (even if those are taken over, say, the retina of an eye that is focused on a photoelectric cell that is designed to register the presence of a quantum event), then most any measurement would inherently involve nonnegative counts, to the extent that exceptions would be so rare as to be regarded as inconceivable.

\section{Discussion}
 

Surprisingly, the fundamental equations of quantum field theory can be modeled to arbitrary precision by extensions of random walks on a random graph, generated with rules so simple that any child who is familiar with the concept of negative numbers can follow and understand them.  

The path-integral approach to quantum field theory by Dirac \cite{DIR33}, and Feynman \cite{FEY65}  is elegant and intuitively appealing; however, it has left physicists -- and in particular, those holding to the notion that physics should be fundamentally simple at some basic level -- with the philosophical quandary of assuming that any subatomic particles, created and annihilated in vast numbers even in tiny regions of the spacetime vacuum, are each able to traverse all possible paths to a given destination in the course of their propagation, including paths implying arbitrarily supraliminal speeds, and integrate at each instant of these paths a Lagrangian. It is a task that the smartest minds could not accomplish, and yet every single field particle is somehow able to execute it. (Admittedly, it was not much better during the reign of Newtonian mechanics, wherein even the tiniest particles were implicitly assumed to have sufficient expertise in vector addition and trigonometry to the extent that they are capable of retaining and transferring energy, momentum, and angular momentum in any collision or interaction, with seemingly infinite precision \cite[Note 5]{DAR17}.)

Numerous physicists have voiced the notion, in the words of Anthony Zee, that at sufficiently small time scales “new physics must appear” \cite[p.~172]{ZEE10}, and various proposals have been offered in anticipation of a need for a more generalized (and ultimately simpler) paradigm of motion. For example, cellular automata have for decades been regarded by some as conceptually a simpler and more plausible foundation for subatomic behavior \cite{TOF03} \cite{THO16} \cite{WOL02} \cite{MAR03} \cite{TOF87} \cite{FRE92}, but there has been little effort devoted to the specifics of how the rules that govern cellular automata could recreate the fundamental  and well-established equations of modern physics. 

With regard to the heat equation, and the summation of paths implicit in its kernel, there is, on the contrary, no particular mystery or vagueness as to how particles acting with simple rules suffice to explain that equation and its many variations, in ways that are accessible to anyone familiar with Pascal's triangle and its continuum limit. This is in sharp contrast to the Schr\"o dinger equation, which is often characterized as formally being simply a heat equation with an imaginary scale factor $k$, without any elaboration as to what, if anything, that might mean in a physical sense.

The formalism presented here builds on that clear connection between what is fundamentally simple and discrete, and the collective phenomena that emerge therefrom, and extends it to the wave equation and other relativistic quantum wavefunctions of field theory, and thereby simplifies Dirac's and Feynman's path-integral formulation (\cf~\cite{TOF03a}). And unlike the continuous-case propagators first derived by Hadamard \cite[Chapter~4.4]{COU89} which have distinctly different solutions according to the dimension of the space, the propagators shown here  are easily generalized to any number of dimensions, as is true of the kernel of the heat equation to which they are related.  Moreover, the simplified approach explains how the interactions \textit{between} fields -- which as Feynman’s diagrams have shown, can be expanded in graph-theoretic terms -- are similar in nature to the interactions \textit{within} a single field that govern its internal propagation even in the absence of outside interactions. And all this without so much as a single complex or floating-point number. Granted, the introduction of quasi-probabilities -- even those that converge to Newtonian values at sufficiently macroscopic scales -- is something that will always be difficult to understand for macroscopic creatures living in a visibly Newtonian world, but there is nothing particularly surprising or unfathomable in recognizing that those who live in a macroscopic realm have powers of perception that conform to those scales. 

The formalism presented here, while poorly suited as a computational scheme (though that may change with the emergence of quantum computing or massively parallel computation), clearly demonstrates how radically simple particles motions can, not by complexity, but by sheer number, recreate the fundamental wave equations of modern physics. 

\section{Conclusion}\label{sec13}
Not applicable





%%\end{document}


\backmatter

\bmhead{Supplementary information}
Not applicable
%%If your article has accompanying supplementary file/s please state so here. 

%%Authors reporting data from electrophoretic gels and blots should supply the full unprocessed scans for key as part of their Supplementary information. This may be requested by the editorial team/s if it is missing.

%%Please refer to Journal-level guidance for %%any specific requirements.

\bmhead{Acknowledgments}

Thanks to Tommasso Toffoli, Norman H. Margolus, and Henry Reich.

\section*{Declarations}

Not applicable

\begin{itemize}
\item Funding
Not applicable
\item Conflict of interest/Competing interests
Not applicable
\item Ethics approval 
Not applicable
\item Consent to participate
Not applicable 
\item Consent for publication
Not applicable
\item Availability of data and materials%\item Code availability
Not applicable
\item Authors' contributions
Not applicable
\end{itemize}

%%\noindent
%%If any of the sections are not relevant to your manuscript, please include the heading and write `Not applicable' for that section. 

%%===================================================%%
%% For presentation purpose, we have included        %%
%% \bigskip command. please ignore this.             %%
%%===================================================%%
%%\bigskip
%%\begin{flushleft}%
%%Editorial Policies for:

%%bigskip\noindent
%%Springer journals and proceedings: %%\url{https://www.springer.com/gp/editorial-policies}

%%\bigskip\noindent
%%Nature Portfolio journals: %%\url{https://www.nature.com/nature-research/editorial-policies}

%%\bigskip\noindent
%%\textit{Scientific Reports}: %%\url{https://www.nature.com/srep/journal-policies/editorial-policies}

%%\bigskip\noindent
%%BMC journals: %%\url{https://www.biomedcentral.com/getpublished/editorial-policies}
%%\end{flushleft}

%%\begin{appendices}

%%\section{Section title of first appendix}\label{secA1}

%%An appendix contains supplementary information that is not an essential part of the text itself but which may be helpful in providing a more comprehensive understanding of the research problem or it is information that is too cumbersome to be included in the body of the paper.

%%=============================================%%
%% For submissions to Nature Portfolio Journals %%
%% please use the heading ``Extended Data''.   %%
%%=============================================%%

%%=============================================================%%
%% Sample for another appendix section			       %%
%%=============================================================%%

%% \section{Example of another appendix section}\label{secA2}%
%% Appendices may be used for helpful, supporting or essential material that would otherwise 
%% clutter, break up or be distracting to the text. Appendices can consist of sections, figures, 
%% tables and equations etc.

%%\end{appendices}

%%===========================================================================================%%
%% If you are submitting to one of the Nature Portfolio journals, using the eJP submission   %%
%% system, please include the references within the manuscript file itself. You may do this  %%
%% by copying the reference list from your .bbl file, paste it into the main manuscript .tex %%
%% file, and delete the associated \verb+\bibliography+ commands.                            %%
%%===========================================================================================%%

\bibliography{sn-bibliography}% common bib file
%% if required, the content of .bbl file can be included here once bbl is generated
%%\input sn-article.bbl

%% Default %%
%%\input sn-sample-bib.tex%

\end{document}
