\section{Introduction}

The numerical simulation of plasmas using macroparticle models\footnote{
We distinguish between a macroparticle model -- which represents large numbers of point particles with a single particle of finite spatial extent -- and a particle-in-cell model -- which places those finite spatial extent particles into a spatial grid and carries out charge deposition by interpolating to the grid. Thus, particle-in-cell models are a subset of macroparticle models. In spectral algorithms, both options are viable, although in this paper we consider only gridless methods.
}~\cite{hockney_eastwood:89, birdsall_langdon:85} is a common tool for applications ranging from plasma processing to particle accelerators, fusion to laser-plasma interactions. The traditional approach to such algorithms starts with discretizing the Lorentz force law in time, assuming the fields are given everywhere in space. The Poisson or Maxwell equations are then discretized in space and time, and updated synchronous to when the Lorentz force integrator requires them. The fields used for the Lorentz force law are interpolated from the discrete spatial points at the proper time. Depositing the source terms is then chosen to carefully avoid self-forces on a single particle and conserve the local charge density. The update sequence is also frequently chosen to conserve momentum, energy, or both.

This process is prone to unphysical heating due to various instabilities having to do with the grid, the time discretization, or both. They also impose conservation laws which do not arise naturally from the discretization. Essentially, the problem is that when discretizing the equations of motion there are too many degrees of freedom left to the algorithm designer which may lead to unphysical behavior.

This can be remedied by deriving the algorithms from a Hamiltonian least-action principle. Such multisymplectic algorithms, described most recently by Shadwick, Stamm, and Evstatiev~\cite{shadwick:14a, shadwick:14b}, respect the fact that the dynamics of both particles and fields~\cite{marsdenPatrickShkoller:98,marsdenPekarskyShkollerWest:01} have a symplectic geometry. The benefits of symplectic integration for single-particle mechanics are well-known in the accelerator physics community (see, for example, \cite{forest:06} for an historical overview of the subject and its applications). The idea of multisymplectic particle-in-cell algorithms has generated considerable interest, and a number of one-dimensional examples on a discrete spatial grid have recently been published~\cite{kraus:13,qin:15a,squire_etal:12,xiao_etal:15}

The chief benefit of multisymplectic algorithms is that their solutions must be a solution for some Lagrangian which approximates the continuous, ``real world'' Lagrangian. The numerical solutions are exact solutions to an approximate Lagrangian. This contributes to long-time stability over many oscillations of the system. In the case of plasmas, such algorithms promise to make simulations extremely reliable over hundreds or even thousands of plasma periods.

By making the discretizing approximations in the action integral, and then deriving the resulting equations of motion, we remove many of the undetermined components that can appear when approximations are made to the equations of motion directly. For example, charge deposition and force interpolation arise from a common term in the Lagrangian, and thus any approximations made to the particles and fields will relate these two constructively. The textbook method is to insist on a deposition/interpolation scheme which does not produce self-forces on the particles. This may or may not be consistent with the underlying physics which produces the equations of motion.

What is important to understand here is that the familiar physical principles -- conservation of momentum, conservation of energy, conservation of angular momentum, and on -- are encoded not in the equations of motion but in the symmetries of the action. They appear in the equations of motion, but arise because of symmetries under translation in space and time, rotations, and on, per Noether's Theorem. If we want our algorithms to contain these conservation laws consistently, then our algorithms should come from some action principle which has the appropriate symmetries.

In this paper we introduce an electrostatic canonical gridless macroparticle algorithm using a spectral decomposition of the electrostatic fields. We express the Lagrangian in terms of the electrostatic potential and particle phase space density (Section~\ref{lowLagrangian}), then discretize the Lagrangian into Fourier modes for the fields and macro-particles for the phase space density (Section~\ref{macroparticles}). In Section~\ref{algorithm}, the resulting action integral is discretized into a second order Riemann sum, and we use the discrete Euler-Lagrange equations to arrive at a second order in time spectral algorithm. The linear numerical dispersion relation is derived in Section~\ref{numerical_dispersion}. We demonstrate the momentum and energy conservation properties for a two-dimensional two-body problem in Section~\ref{twobody} and for a thermal plasma in Section~\ref{thermal_plasma}. Finally, we discuss some optimization strategies for implementation in Section~\ref{implementation}.

\section{The electrostatic Lagrangian} \label{lowLagrangian}

Single particle symplectic algorithms can be derived either from a sequence of canonical transformations~\cite{ruth:83}, splitting methods or so-called ``kick codes''~\cite{schachinger_talman:87}, or from a discrete action minimization~\cite[and references therein]{marsden_west:01}. Regardless of the formalism used to derive the symplectic algorithm, the numerical solutions are an exact solution to the equations of motion derived from an approximate action. This is a much tighter constraint than schemes which start with the equations of motion, where any arbitrary discretization scheme likely fails to respect the geometric properties required by a system described by a Hamiltonian least-action principle. This is what gives them their long-term stability.

In this paper, we consider the Lagrangian treatment of the least-action principle. This contains the same geometric structures as a Hamiltonian approach, but will be more convenient for future work with electromagnetics, where the canonical momentum's relation to the vector potential makes explicit integration schemes much more difficult.

We begin with the electrostatic Lagrangian for a system of particles and fields, due to Low\footnote{
We use C.G.S. units throughout this paper.
}~\cite{low:58}
\begin{eqnarray} \label{lowlagrangian}
\nonumber\mathcal{L} =& \int d\mathbf{x}_0 d \mathbf{v}_0 \left [ \frac{1}{2} m \left ( \frac{\partial \mathbf{x}}{\partial t}(\mathbf{x}_0, \mathbf{v}_0) \right )^2 - e \varphi(\mathbf{x})\right ] f(\mathbf{x}_0, \mathbf{v}_0) + \dots \\
&\dots + \frac{1}{8 \pi} \int d\mathbf{x} \nabla \varphi \cdot \nabla \varphi .
\end{eqnarray}
An auxiliary requirement is that the phase space volume is conserved, \emph{i.e.} $f(\mathbf{x}_0, \mathbf{v}_0, 0) = f(\mathbf{x}, \mathbf{v}, t)$, which is a statement of the Vlasov equation. The continuum equations of motion come from minimizing the action $\mathcal{S} = \int dt \mathcal{L}$ with respect to $\mathbf{x}$ and $\varphi$, and the resulting equations are the familiar Newton's Second Law for the particle characteristics, and the Poisson equation:
\numparts
\begin{eqnarray}
\frac{\partial}{\partial t} \frac{\delta \mathcal{L}}{\delta \mathbf{x}'} - \frac{\delta \mathcal{L}}{\delta \mathbf{x}} = 0 \rightarrow m \frac{\partial^2 \mathbf{x}}{\partial t^2} + e\nabla \varphi(\mathbf{x}) = 0 \\
\frac{\partial}{\partial_\mu} \frac{\delta \mathcal{L}}{\delta \varphi_\mu} - \frac{\delta \mathcal{L}}{\delta \varphi} = 0 \rightarrow\nabla^2 \varphi = 4 \pi e \int d \mathbf{v} f(\mathbf{x}, \mathbf{v}, t)
\end{eqnarray}
\endnumparts

To derive our discrete equations of motion from an action, we must first decompose $\mathcal{L}$ into a discrete set of coefficients the basis functions of our phase space density and scalar potential, and then approximate $\mathcal{S}$ as some discrete approximation of the action integral. The former is a mode decomposition, while the later is an approximation of the action integral using Riemann sums. Once this is done, we will use the discrete Euler-Lagrange equations from~\cite{marsden_west:01} to get the discrete equations of motion.

\section{Spatial discretization of the Lagrangian}\label{macroparticles}

Two fields must be discretized -- the phase space density and the scalar potential. Discretizing $\varphi$ will specify the Poisson solving algorithm, while discretizing the phase space density will give us a definition of macroparticles.

We discretize the scalar potential in a discrete Fourier-type basis, allowing for shape functions in Fourier space:
\begin{equation}
\varphi(\mathbf{x}) = \frac{1}{\sqrt{2 \pi}^D}\int d\mathbf{k} \sum_{\sigma} e^{i \mathbf{k} \cdot \mathbf{x}} \Psi(\mathbf{k} - \mathbf{k}_\sigma) \tilde{\varphi}_\sigma
\end{equation}
where $\sigma$ is a collective index for the discrete modes and $\Psi$ is a generic shape function for the Fourier modes.

The choice of spectral shape function can be made to minimize ringing. For example, a delta function will yield the unattenuated $\cos (k_\sigma x)$ oscillation. By using tent functions separated by $\Delta k$, on the other hand, it is possible to accurately represent the Fourier transform of the fields in a piecewise linear fashion. This puts an envelope on the oscillations, and damps out ringing. Other options are also possible in this formalism, although we must consider their specific properties, as we do in \S~\ref{algorithm} below. In this way, the delta function option may prove the fastest, although it has some issues with ringing. An example of the tent function shape in real space is given in figure (\ref{tentfunction}).

\begin{figure}
\centering
\includegraphics[scale=0.5]{psiofx.png}
\caption{The real and imaginary parts in the $x$ space for a tent shape function $\Psi(k - k_\sigma)$. The envelope drops off as $(x ~\Delta k)^{-2}$ for the tent function.}
\label{tentfunction}
\end{figure}

To discretize the phase space density, we take shape functions in real space for the positions, and use delta functions in velocity:
\begin{equation}
f(\mathbf{x}, \mathbf{v}, t) = \sum_j w_j \Lambda(\mathbf{x} - \mathbf{x}_j(t)) \delta(\mathbf{v} - \mathbf{v}_j(t))
\end{equation}
where $\mathbf{x}_j$ and $\mathbf{v}_j$ are the co\"{o}rdinates and velocities for the $j^{th}$ macroparticle, which has a weight $w_j$. Here, we use $\mathbf{v}$ as shorthand for $\partial_t \mathbf{x}$. We can understand macroparticles in this context as a finite element Lagrangian picture of the flow of the phase space fluid.

Inserting these expressions into the action, eqn. (\ref{lowlagrangian}), yields a discrete Lagrangian whose independent dynamical variables are the Fourier mode coefficients of the scalar potential and the co\"ordinates and velocities of the macroparticles:
\begin{eqnarray}
\nonumber \mathbf{L}_D =& \sum_j w_j \left \{ \frac{1}{2} m \left ( \frac{\partial \mathbf{x}_j}{\partial t} \right )^2 - e\sum_\sigma \int d\mathbf{k} \tilde{\Lambda}^*(\mathbf{k})  e^{i \mathbf{k} \cdot \mathbf{x}_j} \Psi(\mathbf{k} - \mathbf{k}_\sigma) \tilde{\varphi}_\sigma \right \} \\
&+ \frac{1}{8 \pi} \sum_\sigma \sum_{\sigma'} \tilde{\varphi}_\sigma \tilde{\varphi}_{\sigma'} \int d \mathbf{k} ~ \mathbf{k} \cdot \mathbf{k} \Psi(\mathbf{k} - \mathbf{k}_{\sigma}) \Psi(-\mathbf{k} - \mathbf{k}_{\sigma'}).
\end{eqnarray}
We have used the orthogonality of the Fourier modes when integrating over $\mathbf{x}$ in the scalar potential energy term, and the particle shape function Fourier transform, $\tilde{\Lambda}(\mathbf{k})$, is defined the same as the Fourier transform for $\varphi$ above. Specifically, this requires that
\begin{equation}
\Lambda(\mathbf{x}) = \frac{1}{\sqrt{2 \pi}^D}\int d\mathbf{k} e^{i \mathbf{k} \cdot \mathbf{x}} \tilde{\Lambda}(\mathbf{k}) \rightarrow \tilde{\Lambda}(\mathbf{k}) = \frac{1}{\sqrt{2 \pi}^D} \int d\mathbf{x} e^{-i \mathbf{k} \cdot \mathbf{x}} \Lambda(\mathbf{x})
\end{equation}
from the Dirac delta function identity $\delta(\mathbf{x}) = (2 \pi)^D \int d\mathbf{k} e^{i \mathbf{k} \cdot \mathbf{x}}$. For convenience we will define the $k$-squared matrix as
\begin{equation}
\mathcal{K}_{\sigma \sigma'} \equiv \int d \mathbf{k} ~ \mathbf{k} \cdot \mathbf{k} \Psi(\mathbf{k} - \mathbf{k}_\sigma) \Psi(-\mathbf{k} - \mathbf{k}_{\sigma'}).
\end{equation}

If we consider the simple discrete Fourier transform, with $\Psi(\mathbf{k} - \mathbf{k}_\sigma) = \delta (\mathbf{k} - \mathbf{k}_\sigma)$, this becomes an uncoupled Poisson equation, where modes are coupled only to their complex conjugates. This is the functional form of the continuum Lagrangian, which is invariant under translations and therefore conserves momentum. For a more general shape function, $\mathcal{K}$ has off-diagonal elements that couple different modes. This matrix must be inverted once at the beginning of the simulation to compute the fields. It also requires evaluating various convolution integrals between the particle and field shapes, which may prove computationally prohibitive if an analytical expression is not available.

\section{Discretizing the action and the discrete equations of motion}\label{algorithm}

The continuous-time action integral is now approximated by a discrete Lagrangian which is still a function of continuous time:
\begin{equation}
\mathcal{S} \approx \int dt \mathbf{L}_D [\mathbf{x}_j, \mathbf{v}_j; \tilde{\varphi}_\sigma].
\end{equation}
To find the discrete equations of motion, we must approximate the integral for $\mathcal{S}$ using discrete steps in time, using Riemann sums or higher order approximations. The most direct approach is to use composition methods to make a self-adjoint approximation of the action (see \S 2.4 and \S2.5 of~\cite{marsden_west:01}).

A first order approximation of the action is the Riemann sum
\begin{equation}
\mathbf{S}_D^{(n)} = h \mathbf{L}_D \left [ \mathbf{x}_j^{(n+1)}, \frac{\mathbf{x}_j^{(n+1)} - \mathbf{x}_j^{(n)}}{h}; \tilde{\varphi}_{\sigma}^{(n+1)} \right ]
\end{equation}
for a step size $h$. This is the action for the $n^{th}$ step, and the total action is approximated by $\mathcal{S} \approx \sum_n \mathbf{S}_D^{(n)} + \mathcal{O}(h^2)$. This is a first order approximation. To reach second order accuracy in time, we use a composite scheme of a step with its adjoint. The adjoint of the discrete action is given by
\begin{equation}
\mathbf{S}_D^*(\mathbf{x}^{(n+1)}, \mathbf{x}^{(n)}, h) = -\mathbf{S}_D (\mathbf{x}^{(n)}, \mathbf{x}^{(n+1)}, -h)
\end{equation}
and it is straightforward to show that a self-adjoint action is of an even order of accuracy -- thus by being self-adjoint we are assured at least second order accuracy in time.

By taking the composite action for the $n^{th}$ step
\begin{eqnarray}
\nonumber \mathbf{S}_D = & \underbrace{\sum_i \frac{1}{2} m w_i \left \{  \frac{\left (\mathbf{x}_i^{(n+1)} - \mathbf{x}_i^{(n+\nicefrac{1}{2})}\right )^2}{\nicefrac{h}{2}} + \frac{\left (\mathbf{x}_i^{(n+\nicefrac{1}{2})} - \mathbf{x}_i^{(n)}\right )^2}{\nicefrac{h}{2}} \right \}}_{\textrm{ptcl. move}} \\
\nonumber& - \underbrace{h q\sum_i \sum_\sigma w_i \tilde{\varphi}_{\sigma}^{(n+\nicefrac{1}{2})} \int d\mathbf{k}~ \tilde{\Lambda}^* (\mathbf{k}) \Psi(\mathbf{k} - \mathbf{k}_\sigma) e^{i \mathbf{k} \cdot \mathbf{x}_i^{(n+\nicefrac{1}{2})}}}_{\textrm{interpolation/deposition}}
\\
&+ \underbrace{\frac{h}{8 \pi} \sum_{\sigma, \sigma'} \tilde{\varphi}^{(n+\nicefrac{1}{2})}_{\sigma} \mathcal{K}_{\sigma \sigma'} \tilde{\varphi}^{(n+\nicefrac{1}{2})}_{\sigma'}}_{\textrm{Poisson equation}}
\end{eqnarray}
and minimizing the action using the discrete Euler-Lagrange equations (see \S 2.5 of \cite{marsden_west:01}) we get the discrete equations of motion
\numparts
\begin{eqnarray} \label{del_1}
\nonumber m \frac{\mathbf{x}^{(n+1)}_j - 2 \mathbf{x}^{(n+\nicefrac{1}{2})}_j + \mathbf{x}^{(n)}_j}{\nicefrac{h}{2}} = \dots\\
\dots - hei \sum_{\sigma} \int d\mathbf{k} ~ \mathbf{k}  \tilde{\Lambda}^*(\mathbf{k}) \Psi(\mathbf{k} - \mathbf{k}_\sigma) \tilde{\varphi}_\sigma e^{i \mathbf{k} \cdot \mathbf{x}^{(n+\nicefrac{1}{2})}_j} \\ \label{del_2}
\sum_{\sigma'}\mathcal{K}_{\sigma \sigma'} \tilde{\varphi}_{\sigma'} = 8 \pi e \sum_j w_j  \int d\mathbf{k} \tilde{\Lambda}^*(\mathbf{k})  e^{i \mathbf{k} \cdot \mathbf{x}_j^{(n+\nicefrac{1}{2})}} \Psi(\mathbf{k} - \mathbf{k}_\sigma)\\
\frac{\mathbf{x}^{(n+\nicefrac{3}{2})} - \mathbf{x}^{(n+1)}}{\nicefrac{h}{2}} = \frac{\mathbf{x}^{(n+1)} - \mathbf{x}^{(n+\nicefrac{1}{2})}}{\nicefrac{h}{2}}
\end{eqnarray}
\endnumparts
where $\mathcal{K}_{\sigma \sigma'}$ is the generalization of the usual $\mathbf{k}^2$ in the Poisson equation. The first equation is simply $m \ddot{\mathbf{x}} = q \mathbf{E}$, while the second is the Poisson equation for computing the potential in the middle of the step. The third equation tells us how the boundaries between steps are related. In our case it is the continuity of the ``velocity'', although different discretizations may have different boundary terms.

Note that there is no velocity co\"{o}rdinate -- as discussed by Marsden and West~\cite{marsden_west:01} there is no tangent space in the discrete Lagrangian picture. We may, however, introduce the velocity as a convenient intermediate computational variable. If we define
\begin{equation}
\mathbf{v}^{(n)} = \frac{\mathbf{x}^{(n+\nicefrac{1}{2})} - \mathbf{x}^{(n)}}{\nicefrac{h}{2}}
\end{equation}
then the above equations become
\numparts
\begin{eqnarray}
\mathbf{x}^{(n+\nicefrac{1}{2})}_j = \mathbf{x}^{(n)}_j + \frac{h}{2} \mathbf{v}^{(n)}_j \\
\sum_{\sigma'}\mathcal{K}_{\sigma \sigma'} \tilde{\varphi}_{\sigma'} = 8 \pi e \sum_j w_j  \int d\mathbf{k} \tilde{\Lambda}^*(\mathbf{k})  e^{i \mathbf{k} \cdot \mathbf{x}_j^{(n+\nicefrac{1}{2})}} \Psi(\mathbf{k} - \mathbf{k}_\sigma)\\
\mathbf{v}^{(n+\nicefrac{1}{2})}_j - \mathbf{v}^{(n)}_j = -\frac{h e}{m} i \sum_{\sigma} \int d\mathbf{k} ~ \mathbf{k}  \tilde{\Lambda}^*(\mathbf{k}) \Psi(\mathbf{k} - \mathbf{k}_\sigma) \tilde{\varphi}_\sigma e^{i \mathbf{k} \cdot \mathbf{x}^{(n+\nicefrac{1}{2})}_j} \\
\mathbf{x}^{(n+1)}_j = \mathbf{x}^{(n+\nicefrac{1}{2})}_j + \frac{h}{2} \mathbf{v}^{(n+\nicefrac{1}{2})}_j\\
\mathbf{v}^{(n+1)}_j = \mathbf{v}^{(n+\nicefrac{1}{2})}_j
\end{eqnarray}
\endnumparts
which is the split-step move and accelerate familiar from most particle-based algorithms. We have avoided using a grid to deposit charge to preserve the translational invariance. This means that every particle has an exponential, $e^{i \mathbf{k} \cdot \mathbf{x}}$, that must be calculated.

We have thus far left $\Psi$ almost completely general. To understand conditions on momentum conservation, we derive the total force between two particles. First, let us assume that we have chosen our discrete $\mathbf{k}_{\sigma}$ such that there exists a $\mathbf{k}_{-\sigma} = - \mathbf{k}_{\sigma}$. This is a reasonable choice if we want to make the shape functions $\Lambda$ and $\Psi$ even functions and still return real values for $\varphi(\mathbf{x})$.

If we place a single particle at position $\mathbf{x}_j$, it will create the potential
\begin{eqnarray}
\nonumber \tilde{\varphi}_{\sigma} &=& 8 \pi e w_j \sum_{\sigma'} \mathcal{K}_{\sigma' \sigma}^{-1} \int d\mathbf{k} \tilde{\Lambda}^*(\mathbf{k}) e^{i \mathbf{k} \cdot \mathbf{x}_j} \Psi(\mathbf{k} - \mathbf{k}_{\sigma'})\\
 &=& 8 \pi e w_j \sum_{\sigma'} \mathcal{K}_{\sigma' \sigma}^{-1} f_{\sigma'}(\mathbf{x}_j)
\end{eqnarray}
where
\begin{equation}
f_{\sigma}(\mathbf{x}) = \int d \mathbf{k} \tilde{\Lambda}^* (\mathbf{k}) \Psi(\mathbf{k} - \mathbf{k}_\sigma) e^{i \mathbf{k} \cdot \mathbf{x}_j}
\end{equation}
The force on particle $\ell$ due to particle $j$ is then given by
\begin{equation}
\Delta \mathbf{p}_{j \rightarrow \ell} = 8 \pi h e^2 w_\ell w_j \nabla_{\mathbf{x}_\ell} \sum_{\sigma \sigma'} f_{\sigma'}(\mathbf{x}_j) \mathcal{K}_{\sigma \sigma'}^{-1} f_{\sigma}(\mathbf{x}_\ell)
\end{equation}

If we define
\begin{equation}
F \left (\mathbf{k}, \mathbf{k}' \right ) = \tilde{\Lambda}^*(\mathbf{k}) \tilde{\Lambda}^*(\mathbf{k}') i (\mathbf{k} + \mathbf{k}') e^{i (\mathbf{k} \cdot \mathbf{x}_j + \mathbf{k}' \cdot \mathbf{x}_\ell)}
\end{equation}
then we can then say that the total change of the momentum for this two-body pair is
\begin{eqnarray}
\nonumber \Delta \mathbf{p}_{j \rightarrow \ell} + \Delta \mathbf{p}_{\ell \rightarrow j} & = A \int_{-\infty}^{\infty} d \mathbf{k} d \mathbf{k}' ~ F(\mathbf{k}, \mathbf{k}') \sum_{\sigma, \sigma'} \Psi(\mathbf{k} - \mathbf{k}_\sigma) \mathcal{K}_{\sigma \sigma'}^{-1} \Psi(\mathbf{k}' - \mathbf{k}_{\sigma'})\\
\nonumber & = A \int_{-\infty}^{\infty} d \mathbf{k} d \mathbf{k}' ~ F(\mathbf{k}, \mathbf{k}') \sum_{\sigma, \sigma'} \Psi(\mathbf{k} - \mathbf{k}_{-\sigma}) \mathcal{K}_{-\sigma -\sigma'}^{-1} \Psi(\mathbf{k}' - \mathbf{k}_{-\sigma'}) \\
& = A \int_{-\infty}^{\infty} d \mathbf{k} d \mathbf{k}' ~ F(\mathbf{k}, \mathbf{k}') \sum_{\sigma, \sigma'} \Psi(\mathbf{k} + \mathbf{k}_{\sigma}) \mathcal{K}_{\sigma \sigma'}^{-1} \Psi(\mathbf{k}' + \mathbf{k}_{\sigma'})
\end{eqnarray}
where we have used to symmetry of $\mathcal{K}_{-\sigma, -\sigma'} = \mathcal{K}_{\sigma, \sigma'}$ and the fact that $\sigma$ and $\sigma'$ are dummy indices. We can then flip the sign on $\mathbf{k}$ and $\mathbf{k}'$ and exploit the fact that $\Psi$ is an even function to get that
\begin{eqnarray}
\nonumber \Delta \mathbf{p}_{j \rightarrow \ell} + \Delta \mathbf{p}_{\ell \rightarrow j} &= A \int_{-\infty}^{\infty} d \mathbf{k} d \mathbf{k}' ~ F(-\mathbf{k}, -\mathbf{k}') \sum_{\sigma, \sigma'} \Psi(\mathbf{k} - \mathbf{k}_{\sigma}) \mathcal{K}_{\sigma \sigma'}^{-1} \Psi(\mathbf{k}' - \mathbf{k}_{\sigma'}) \\
&= (\Delta \mathbf{p}_{j \rightarrow \ell} + \Delta \mathbf{p}_{\ell \rightarrow j})^*
\end{eqnarray}
This proves that the total change in momentum is a real number if we choose even functions for our particle and field shapes, and have negative pairs of Fourier modes. However, this will not generally be zero, and we can have some change in the total momentum. This is because the resulting discrete Lagrangian is not, in general, unchanged by translations, so it does not have the associated Noether conserved quantity.

For the special case of a delta function, the Lagrangian is unchanged under a translation, and we can show immediately that momentum is conserved to machine precision. Here, $\mathcal{K}_{\sigma \sigma'} = \mathbf{k}_{\sigma} \cdot \mathbf{k}_{\sigma'} \delta_{\sigma, -\sigma'}$. There are no self-forces
\begin{equation}
\Delta \mathbf{p}_{j \rightarrow j} = 4 \pi h e^2 w i \sum_{\sigma} \frac{\mathbf{k}_{\sigma}}{|\mathbf{k}_\sigma|^2} e^{i \mathbf{k}_\sigma \cdot \mathbf{x}_j} e^{-i \mathbf{k}_{\sigma} \cdot \mathbf{x}_j} = 0
\end{equation}
since every $\mathbf{k}_\sigma$ is paired with a $-\mathbf{k}_{\sigma} = \mathbf{k}_{-\sigma}$, a requirement on the expansion to keep $\varphi(\mathbf{x})$ real. Because there are no self-forces, we can write the force on particle $\ell$ from particle $j$ as
\begin{equation}
\Delta \mathbf{p}_{\ell \rightarrow j} = 4 \pi h e^2 w_\ell w_j i \sum_\sigma \mathbf{k}_{\sigma} e^{i \mathbf{k}_{\sigma} \cdot \mathbf{x}_j} e^{-i \mathbf{k}_{\sigma} \cdot \mathbf{x}_{\ell}} \frac{1}{|\mathbf{k}_\sigma|^2} = - \Delta \mathbf{p}_{j \rightarrow \ell}.
\end{equation}
We must therefore balance our desire for preventing the ringing from using discrete Fourier modes with how much our problem requires exact or nearly-exact momentum conservation. It remains to see whether such systems will grow monotonically in momentum, or whether the total momentum will simply oscillate around the initial total momentum, and how these depend on our shape functions. For now, we will focus on the discrete modes with delta-function shape, as they are the least computationally demanding due to the lack of any numerical integrals in the update sequence.

\section{Numerical Dispersion Relation}\label{numerical_dispersion}

We now derive a linearized dispersion relation for the finite time step spectral algorithm, extending the arguments of \S 9-2 of~\cite{birdsall_langdon:85} for our spectral algorithms. We begin by defining $\mathbf{x} = \mathbf{X} + \delta \mathbf{x}$, where $\mathbf{X}$ is the unperturbed ballistic motion, $\mathbf{X}^{(n)} = \mathbf{x}^{(0)} + \mathbf{v}^{(0)} n h$ and $|\mathbf{k} \cdot \delta \mathbf{x}| \ll 1$ is a small perturbation.

To include the boundary terms, we must include the update sequence from the end of the last ``accelerate'' to the beginning of the next. For convenience, we will denote these with integer indices, since they are all offset by a half-step. Inserting this approximation into eqn.~(\ref{del_1}), along with our re-indexing, gives
\begin{equation}
\frac{\delta \mathbf{x}_j^{(n+1)} - 2 \delta \mathbf{x}_j^{(n)} + \delta \mathbf{x}_j^{(n-1)}}{h} = - \frac{hei}{m} \sum_{\sigma} \int d\mathbf{k} ~ \mathbf{k}  \tilde{\Lambda}^*(\mathbf{k}) \Psi(\mathbf{k} - \mathbf{k}_\sigma) \tilde{\varphi}_\sigma e^{i \mathbf{k} \cdot \mathbf{X}^{(n)}_j}
\end{equation}

 It is interesting to look at a finite width in $\Psi$. This behaves like a finite width in $\mathbf{v}$ in the particle distribution, creating numerical Landau damping even for a cold plasma, with a damping time $\tau^{-1} \sim \sigma_ \mathbf{k} \cdot \mathbf{v}_0$ for a typical width $\sigma_ \mathbf{k}$ to the shape function $\Psi$. We can intuitively understand this as follow: The finite width Fourier shapes represent a distribution of phase velocities. Thus, by introducing these shapes, we create a spread of phase velocities that a given mode couples to the distribution of particles, which is an unphysical but identical configuration to the usual mechanism for physical Landau damping. This may be useful in some scenarios. However, from here on we will consider only the discrete Fourier modes, $\Psi(\mathbf{k} - \mathbf{k}_\sigma) = \delta(\mathbf{k} - \mathbf{k}_\sigma)$.

  Let us assume that the electrostatic field varies as $\tilde \varphi_\sigma = \tilde{\Phi}_\sigma e^{-i \omega n h}$. The integrand of the right hand side of this equation varies as $e^{i \mathbf{k}_\sigma \cdot (\mathbf{x}_j^{(0)} - \mathbf{v}_j^{(0)} n h) - i \omega n h}$, and so too must $\delta \mathbf{x}$. Under this assumption, we have the dynamics for $\delta \mathbf{x}$ given by
\begin{eqnarray}
\nonumber \frac{\delta \mathbf{x}_j^{(n+1)} - 2 \delta \mathbf{x}_j^{(n)} + \delta \mathbf{x}_j^{(n-1)}}{h} = \\ - \frac{hei}{m} \sum_{\sigma}  \mathbf{k}_\sigma  \tilde{\Lambda}^*(\mathbf{k}_\sigma) \tilde{\Phi}_\sigma e^{i \mathbf{k}_\sigma \cdot (\mathbf{x}^{(0)}_j + \mathbf{v}_j^{(0)}nh)}e^{-i \omega n h}
\end{eqnarray}
It is convenient to define
\begin{equation}
\delta \mathbf{x}^{(n)} = \sum_\sigma \delta \tilde \mathbf{x}^\sigma e^{-i \omega_d^\sigma n h}
\end{equation}
which in turn gives that
\begin{equation}
\delta \tilde \mathbf{x}_j^\sigma = -\frac{1}{2 (\cos \omega_d^\sigma h - 1)} \frac{h^2 e i}{m}  \mathbf{k}_\sigma \tilde \Lambda^*(\mathbf{k}_\sigma) \tilde \Phi_\sigma e^{i \mathbf{k}_\sigma \mathbf{x}_j^{(0)}}
\end{equation}

It remains to compute $\tilde \Phi_\sigma$. Looking to eqn.~(\ref{del_2}) and Taylor expanding the right hand side gives
\begin{equation}
|\mathbf{k}_\sigma|^2 \tilde \Phi^*_{\sigma} = 4 \pi e \sum_j w_j \tilde{\Lambda}^*(\mathbf{k}_\sigma) e^{i \mathbf{k}_\sigma \cdot (\mathbf{x}_j^{(0)} + \mathbf{v}_j^{(0)}n h)} e^{i \omega n h} i \mathbf{k}_\sigma \cdot \delta \mathbf{x}_j
\end{equation}
where we have assumed that the unperturbed quantity
\begin{equation}
\sum_j w_j \tilde{\Lambda}^*(\mathbf{k}_\sigma) e^{i \mathbf{k}_\sigma \cdot (\mathbf{x}_j^{(0)} + \mathbf{v}_j^{(0)}n h)} \approx 0
\end{equation}
which is to say that the unperturbed system has approximate charge neutrality.

Combining the solution for the perturbed orbit with the Poisson equation yields
\begin{eqnarray}
\nonumber |\mathbf{k}_\sigma|^2 \tilde{\Phi}_\sigma = & 4 \pi \frac{h^2 e^2}{m} \sum_{j, \sigma'} w_j \tilde{\Lambda}^*(\mathbf{k}_\sigma) e^{-i (\mathbf{k}_\sigma-\mathbf{k}_{\sigma'}) \cdot \mathbf{x}_j^{(0)}} \times \\
 & e^{i (\omega_d^\sigma - \omega_d^{\sigma'})n h} \mathbf{k}_\sigma \cdot \mathbf{k}_{\sigma'} \tilde \Lambda^*(\mathbf{k}_{\sigma'}) \tilde \Phi_{\sigma'} \frac{1}{2 (\cos \omega_d^{\sigma'} h - 1)}
\end{eqnarray}
We note that the sum over $j$ can be approximated as an integral over the initial conditions weighted by the local phase space density
\begin{equation}
\sum_j w_j * \mapsto \int d\mathbf{x}^{(0)} d \mathbf{v}^{(0)} f(\mathbf{v}^{(0)}) *
\end{equation}
where we have assumed the initial system has no charge perturbations. Noting the orthogonality relation with Fourier modes, we get that
\begin{equation}
\tilde \Phi_\sigma = \frac{4 \pi e^2 n_0}{m} h^2|\tilde \Lambda (\mathbf{k}_\sigma)|^2 \int d\mathbf{v} f(\mathbf{v}) \frac{1}{2(\cos \omega_d^\sigma h - 1)} \tilde \Phi_\sigma
\end{equation}
where $n_0$ is the unperturbed initial charge density. This yields the numerical dispersion relation
\begin{equation}
1 = \omega_p^2h^2 |\tilde \Lambda(\mathbf{k}_\sigma)|^2 \int d \mathbf{v} f(\mathbf{v}) \frac{1}{4 \sin^2 \left (\frac{ (\omega - \mathbf{k}_\sigma \cdot \mathbf{v}) h}{2} \right )}
\end{equation}
for each mode $\mathbf{k}_\sigma$. In the limit of $h \rightarrow 0$ and $\tilde \Lambda = 1$ (point particles) this returns the continuum limit. This is as in~\cite{birdsall_langdon:85} but for discrete time steps but continuous space, although the discrete Fourier representation picks specific $\mathbf{k}_\sigma$, and with the particle shape functions appearing to modify the dispersion relation. We can carry out additional analysis on this dispersion relation.

Immediately we see that this algorithm modifies the plasma frequency for each mode to
\begin{equation}
\omega_p^\sigma = \omega_p |\tilde{\Lambda}(\mathbf{k}_\sigma)|.
\end{equation}
This requires that $|\tilde \Lambda(\mathbf{k}_\sigma)|$ stay close to unity until the Debye wavenumber, which is equivalent to requiring that the particle shape function $\Lambda(\mathbf{x})$ be narrower than the Debye length.

Assuming a cold coasting plasma, $f(\mathbf{v}) = \delta(\mathbf{v} - \mathbf{v}_0)$, the dispersion relation becomes
\begin{equation}
4 \sin^2 \left ( \frac{(\omega - \mathbf{k}_\sigma \cdot \mathbf{v}_0)h}{2} \right ) = \omega_p^2 h^2 |\tilde \Lambda(\mathbf{k}_\sigma)|^2
\end{equation}
which in turn requires that
\begin{equation}
\frac{\omega_p h}{2} |\tilde \Lambda(\mathbf{k}_\sigma)| < 1
\end{equation}
which is the modified $\omega_p-\Delta t$ condition for our electrostatic algorithm. This also gives a guideline for accurate spectral fidelity for the linear plasma dispersion relation (must be much smaller than $1$), and the restriction is most stringent for long wavelength modes, assuming $|\tilde \Lambda(\mathbf{k}_\sigma)|$ has a decreasing envelope with $|\mathbf{k}_\sigma|$.

We can furthermore use Nyquist diagrams to evaluate the numerical Penrose criterion, since the dispersion relation has poles at $\omega - \mathbf{k}_\sigma \cdot \mathbf{v} = n \pi$. However, only the $n=0$ pole will contribute, as the residues of all the others are zero. Since dynamics is only parallel to $\mathbf{k}_\sigma$, let us limit ourselves to a one-dimensional dispersion relation parallel to $\mathbf{k}_\sigma$ for simplicity:
\begin{equation}
1 = \omega_p^2h^2 |\tilde \Lambda(\mathbf{k}_\sigma)|^2 \int d v \frac{f(v)}{(\omega - k_\sigma v)^2} \frac{1}{4 \left ( \frac{\sin \left (\frac{ (\omega - k_\sigma v) h}{2} \right )}{\omega - k_\sigma v} \right )^2}.
\end{equation}
This integral is given by the Cauchy principle value integral
\begin{equation}
1 = \omega_p^2 h^2 |\tilde{\Lambda}(\mathbf{k}_\sigma)|^2 \left (\mathcal{P} \int dv \dots + \pi i \frac{f'(\nicefrac{\omega}{k_\sigma})}{h^2 k_\sigma^2} \right ).
\end{equation}
The imaginary part is of interest here. Because it has the identical functional form to the continuum limit with a modified plasma frequency, we can conclude that there are no additional plasma instabilities introduced by the algorithm so long as the $\omega_p-\Delta t$ condition is satisfied.

That instability arises in the $\mathcal{P}$ part of the integral, which must be able to equal $1$ for real values of $\omega$. The sinusoid in the denominator may prevent that due to the repeated pole structure, although the details will depend on $f(\mathbf{v})$. If we assume a cold but not infinitely cold plasma, and a gaussian distribution $f(v) = \frac{\exp(-\nicefrac{v^2}{2 \sigma_v^2})}{\sqrt{2 \pi \sigma_v^2}}$, with $\sigma_v \ll \nicefrac{1}{h k_\sigma}$, then we can approximate the Cauchy principal value asymptotically as
\begin{equation}
\mathcal{P} \int dv \dots \approx \frac{1}{4 \sin^2 \frac{\omega h}{2}} + \mathcal{O}(h k_\sigma \sigma_v)
\end{equation}
which means thermal effects offer a small modification of to the cold $\omega_p-\Delta t$ condition. It also offers an additional guideline for accurate modeling, which is that for an initially thermal distribution $|h k_\sigma \sigma_v| \ll 1$ for the largest values of $k_\sigma$ or else we could introduce a numerical plasma instability with a sufficiently large time step.

\section{Numerical Example -- Two-Body Problem}\label{twobody}

As a first test of the angular and linear momentum conserving properties of this algorithm, we consider the two dimensional Kepler problem. In two dimensions the particles are line charges, and the electrostatic potential is given by $V(r) \sim \ln (\nicefrac{r}{r_0})$, so this differs from the standard $\nicefrac{1}{r}$ Kepler problem. Our test problem is a single electron $(m_e, e)$ and a ``double-positron'' $(2 m_e, -2 e)$. There is some initial total linear momentum. We use tent functions for the particle shapes and discrete Fourier modes, with $\Psi(\mathbf{k}) = \delta(\mathbf{k})$, with a time step of $\unit[5. \times 10^{-9}]{sec}$ for $2 \times 10^6$ time steps. For the given initial conditions, this is approximately 9-10 steps per period, which reasonably resolves the dynamics of the problem.

\begin{figure}
\begin{center}
\includegraphics[scale=0.5]{energy.png}
\end{center}
\caption{Total energy relative to the initial energy of the two-body system.}
\label{energy}
\end{figure}

As is well-known, if one calculates the energy of a conservative Hamiltonian system from simulations using a symplectic integrator, that energy is bounded for exponential time. We see this behavior in fig. (\ref{energy}). After some initial transient increase in the energy, it is conserved at the \unit[1]{\%} level over an extended time. We can also look at the individual energy contributions -- kinetic, capacitive, and field energies -- to see that each is very accurately bounded, so that there is no drift in the individual components that happens to cancel. The larger scale oscillations are a result of computing the energy every 250 time steps, i.e. once every 25 periods or so. This creates much of the rippling we see in fig. (\ref{energy_envelopes}).

\begin{figure}
\begin{center}
\includegraphics[scale=0.5]{energy_envelopes.png}
\end{center}
\caption{The individual contributions of the kinetic, capacitive, and field energies to the total energy of the Kepler problem.}\label{energy_envelopes}
\end{figure}

Fig. (\ref{energy_envelopes}) demonstrates that the energy conservation does not arise from spurious growth in kinetic energy and decay in potential energy, or \emph{vice versa}. As one would expect, the individual contributions oscillate and are well-behaved. The energy conservation then does not arise from any spurious trends in the individual contributions which happen to cancel out on the whole.

Furthermore, the robustness of the behavior is independent of the timing of the co\"ordinates during the time step when the energy is computed. This means that the energy evaluated at the end of a full time step has the same level of precision (although will differ numerically) as the energy evaluated if a half-drift is used to put $\mathbf{x}$ and $\mathbf{v}$ at the same time, and then the fields are computed using $\mathbf{x}$. This is another property familiar from single particle symplectic integrators.

\begin{figure}
\begin{center}
\includegraphics[scale=0.5]{momentum.png}
\end{center}
\caption{Total $p_x$ and $p_y$ for the system, in a $\log$ scale.}
\label{momentum}
\end{figure}

A major feature of this discrete Fourier spectral algorithm is that it is translationally invariant, therefore it should conserve the total momentum. In fig. (\ref{momentum}) we see that the initial total linear momentum is conserved to machine precision. There is an initial nonzero total momentum, and for the first five hundred thousand time steps the momentum change is machine precision zero. The upward trend after this point is likely due to roundoff error.

The algorithm also purports to be second order accurate, meaning that the correction to the continuous action integral is given by $\mathbf{S}_D \approx \mathcal{S} + \mathcal{O}(h^3)$. Based on scaling, where $\mathcal{S} \sim \int E dt$, this would imply that the correction to the energy would be $E_D \approx E + \mathcal{O}(h^2)$. For the two-body problem, we can clearly see this in fig.~(\ref{energyscaling}).

As can be seen, the energy deviation scales with $h^2$ for this particular problem from around \unit[$10^{-3}$]{sec.} where the failure to resolve the physics becomes catastrophic\footnote{At this point, the relative energy error is $\sim 10^3-10^4$.}, down to around \unit[$10^{-8}$]{sec} where machine precision issues come into play. This clearly demonstrates the second order accuracy of the spectral algorithm. It should be noted that this algorithm was implemented using NumPy (v. 1.9.2) in Python (v. 2.7.10) on a 64-bit machine -- handling of machine precision will vary across architectures and programming languages.

\begin{figure}
\begin{center}
\includegraphics[scale=0.5]{energy_error.png}
\end{center}
\caption{Energy error over a single time step versus $h$.}
\label{energyscaling}
\end{figure}

We have thus shown that we can accurately simulate the two-body Kepler problem for a hundred thousand periods of the system with no change in the total momentum and with very well-behaved energy which does not spuriously grow. Such a problem is very difficult to do in conventional gridded PIC, where the total energy grows and the total momentum is changing, which breaks the careful conservation that allows for a closed orbit two-body problem. We see that this multisymplectic integrator preserves a number of key geometric properties of the system -- the angular and linear momentum are conserved exactly, with good behavior for the total energy due to the nature of symplectic algorithms.

\section{Numerical Example -- Thermal Plasma}~\label{thermal_plasma}

Having demonstrated the energy and momentum conserving properties of this algorithm for a two-body problem, we now consider a thermal plasma. Again, we consider periodic boundaries to keep the system closed and make the total energy a quantity that is conserved in the physical system. In conventional PIC, such a system would show spurious energy growth over many plasma periods.

We show simulation results run for $\sim 10^6$ plasma oscillations. We consider a system with the parameters in fig. (\ref{simparams}), which is a very low resolution simulation. We are resolving Debye length physics here -- in practice the spectral method introduces an ultraviolet cutoff on the resolved physics. Any length shorter than $(k_{\textrm{max}})^{-1}$ is removed entirely, as can be verified by failing to resolve the Debye length.

\begin{figure}
\begin{center}
\begin{tabular}{|c | c|}
\hline
Quantity & Value [units] \\
\hline
electron density & $\unit[10^{17}]{cm^{-3}}$\\
$\lambda_D$ & $\unit[3.96\times 10^{-11}]{cm}$ \\
$\omega_p$ & $\unit[8.14 \times 10^{17}]{sec^{-1}}$\\
\hline
$L_{\textrm{domain}}$ & $4 \times \lambda_D$\\
$n_{\textrm{macro}}$ & 16\\
macro. weight & $9.9 \times 10^5$\\
macro. width & $0.1 \lambda_D$ \\
$\delta k$ & $\nicefrac{\pi}{2} (L_{\textrm{domain}})^{-1}$\\
$n_{\textrm{modes}}$ & 80\\
$h^{-1}$ & $5 \times \nicefrac{\omega_p}{2\pi}$\\
\hline
\end{tabular}
\caption{First simulation parameters.}\label{simparams}
\end{center}
\end{figure}

\begin{figure}
\begin{center}
\includegraphics[scale=0.5]{thermal_energy.png}
\end{center}
\caption{Energy change over time for thermal plasma simulation.}
\label{thermal1energy}
\end{figure}

\begin{figure}
\begin{center}
\includegraphics[scale=0.5]{thermal_momentum.png}
\end{center}
\caption{Momentum change over time for thermal plasma simulation.}
\label{thermal1momentum}
\end{figure}

In this simulation, there is no discernible change in the total energy of the system (fig. (\ref{thermal1energy})). This is to be expected, as the fields should be zero in the physical system, and we would expect all the motion to be ballistic. In two or three dimensions, we do not anticipate this level of energy conservation, which likely arises from a lack of degrees of freedom to cause the energy oscillations observed in the two-dimensional two-body problem. We also see the same momentum conservation, fig. (\ref{thermal1momentum}), which oscillates near machine precision throughout the simulation.

In conventional particle-in-cell simulations, even in one dimension, it is known that the plasma will heat up~\cite{ueda_etal:94, mardahl_thesis}. As noted in \cite{mardahl_thesis}, various techniques in conventional algorithms -- such as higher order particle shapes or increasing the grid resolution of the Debye length -- can reduce the numerical heating up to a point. This algorithm removes the heating entirely, as well as dramatically reducing the statistical fluctuations intrinsic in conventional particle-in-cell.

The reduction in statistical fluctuations can be understood by looking at how the variational method would replicate a conventional algorithm. Conventional algorithms have no intrinsic convolution integrals, so we can understand the ``particle shape'' as either a particle shape with delta function fields, or a ``field shape'' with delta function particles. Both options are extremely noisy. By introducing the convolution integrals and using a Fourier basis, we smooth the particle shapes while simultaneously removing high frequency content without aliasing. This has been demonstrated in finite-difference versions of variational algorithms in, e.g., \cite{shadwick:14a}, where a one dimensional electromagnetic finite differencing shows extremely smooth charge distributions for the applied problem. This comes from using a combination of field and particle shape functions, instead of just one or the other.

\section{Implementation and Optimization}~\label{implementation}

This algorithm has a number of properties that differ from conventional particle-in-cell methods that requires some consideration. This includes the discrete mode properties, parallel implementation, the global nature of the field solve, and evaluation of the deposition/interpolation terms. Also, the spectral nature of the algorithm requires some care to minimize the computational costs.

The first is the span of $k$-space. Implemented with periodic boundaries, the longest wavelength is $L$, the simulation domain length. Furthermore, we must resolve adequately the particle width, $\ell$. Thus, our $k$-space modes must span $k \in [2 \pi/L, \sim 2\pi/\ell]$. Furthermore, because the system has periodic boundaries, we must have as our increment in $k$-space $\Delta k = 2 \pi/L$ to maintain periodicity. Thus, we must use $\mathcal{O}(L/\ell)^D$ Fourier modes. This is comparable to the number of cells required for a traditional gridded method.

Because the algorithm is truly spectral (there is no grid in real space) conventional domain decomposition is irrelevant. This simplifies parallel implementations. Each process must have all the field data, and therefore the field solve can be reduced to a single \texttt{MPI\_AllReduce} call on the $\tilde{\varphi}_\sigma$ after the deposition phase. All other operations can be performed locally. This also makes load balancing trivial -- every process will receive an equal number of macro-particles and store all the field information. However, this could present a challenge for very large simulations with substantial amounts of field data. Because the simulation requires no overhead for messaging particle data, and is perfectly load balanced, such spectral algorithms may be ideal for simulations of devices such as magnetrons which have fields everywhere but particles which are isolated to a small region of space throughout much of the simulation.

Computing the deposition and interpolation can be computationally intensive due to the exponential evaluations. This can be optimized in two ways for the discrete Fourier treatment (Fourier shape functions cannot avoid this problem). The first is by recognizing that we are always taking the real part of the fields, so we only have to keep the $\vec{k}_\sigma$ with all components greater than or equal to zero. This is a factor of $2^D$ savings in the amount of information that must be stored. We can also reduce the number of exponentials per particle from $N_{modes}$ to $D$ if we use regularly spaced Fourier modes spaced by $\nicefrac{2 \pi}{L}$ by tracking $\nicefrac{2 \pi x}{L}$. Then we must evaluate the exponential once, and each coefficient is simply a positive integer power of that exponential. Compared to the the na\"{i}ve algorithm, this can represent a substantial savings in computation and memory footprint.

\section{Conclusion}

We have derived an electrostatic spectral macroparticle algorithm from a variational approach. The resulting algorithm closely resembles conventional leapfrog algorithms, although we have bypassed using a grid -- hence the choice to refer to it as a ``spectral macroparticle algorithm'' instead of a ``particle in cell algorithm''. A broad class of algorithms are available with finite widths in $k$-space for the field modes. However, our analysis has shown that discrete Fourier modes (rather than modes with finite width) are exactly momentum conserving, and furthermore observed that finite width modes introduce a numerical Landau damping term to the dispersion relation. For these reasons, we consider the discrete Fourier algorithm as the ideal implementation. This algorithm has demonstrated remarkable energy conservation properties, as one would expect from symplectic algorithms.

We have demonstrated the discrete Fourier algorithm prevents all spatial discretization instabilities, and is only prone to the $\omega_p-\Delta t$ instability that arises from time discretization. This is therefore a remedy for the finite grid instability. We have also seen these algorithms demonstrated in two two-dimensional simulations: a two-body problem and a thermal plasma, both with periodic boundaries. These are, to the author's best knowledge, the first computational examples of two-dimensional variational algorithms. They demonstrate proper second-order energy error, as one would expect.

Because we considered Cartesian electrostatics, we were not presented with a number of issues important for electromagnetic algorithms, particularly gauge invariance and the associated charge conservation, as well as difficulties in non-Cartesian co\"{o}rdinate systems such as the pole at the origin of spherical or cylindrical co\"{o}rdinates. This represents future work.

Variational algorithms have a number of advantages over conventional plasma algorithms. The primary one is that the equations of motion come from a discrete Euler-Lagrange equation. Because of this, all approximations (spectral decomposition, co\"{o}rdinate transformations for particles and fields, time discretization) are made at the level of the action -- the discrete Euler-Lagrange equations then construct the algorithm from these approximations. This gives more freedom to the algorithm developer, as the approximations are guaranteed to be self-consistent.

\section{Acknowledgements}

The author would like to thank B. Shadwick (UNL) for fruitful discussions on variational algorithms. The author would also like to thank D. Bruhwiler and N. Cook (RadiaSoft) and J.-L. Vay and R. Lehe (LBL) for invaluable discussions concerning the content of this manuscript.

This work was sponsored by the Air Force Office of Scientific Research, Young Investigator Program, under contract no. FA9550-15-C-0031. Distribution Statement A. Approved for public release; distribution is unlimited.


\section{The Basic Model}
In a previous article, the author introduced a free wave equation for discrete quantum gravity \cite{gud15}. According to this model, the geometry of a universe is determined by nature (through a quantum formalism) and it is this geometry that we call gravity. We say that a particle is free if it is not acted on by any forces except gravity and gravity is not a force anyway, it is geometry. The propagation amplitude of a free particle is governed by the free wave equation. We showed in \cite{gud15} that solutions to this equation predict that particles tend to move along geodesics given by the geometry This implies that particles are attracted toward regions of large curvature. In this sense, the mass-energy distribution is determined by the geometry which is opposite to the postulates of classical general relativity theory.

In the present article, we extend the work in \cite{gud15} to include other interactions. Unlike the free wave equation which describes the motion of a single particle, the equations we introduce here describe two interacting particles. The main difference is that the free wave equation involves functions of a single variable, while the interaction wave equations involve functions of two variables. We now briefly review the basic model for this theory. For more details, we refer the reader to \cite{gud13,gud14,gud15}.

We shall model the structure of a universe by a causal set or causet \cite{bdghs, hen09,sor03,sur11,wc15}. Mathematically, a
\textit{causet} is a partially ordered set $(x,<)$. For $a,b\in x$, we interpret $a<b$ as meaning that $b$ is in the causal future of $a$. If $a<b$ and there is no $c\in x$ with $a<c<b$ we say that $a$ is a \textit{parent} of $b$ and write $a\prec b$. Denoting the cardinality of $x$ by
$\ab{x}$, a \textit{labeling} of $x$ is a map $\ell\colon x\to\brac{1,2,\ldots ,\ab{x}}$ such that $a<b$ implies $\ell (a)<\ell (b)$. A labeling of $x$ may be considered a ``birth order'' of the vertices of $x$. A \textit{covariant causet} ($c$-\textit{causet}) is a causet that has a unique labeling \cite{gud13,gud14,gud15}. In this article, we shall only model possible universes by $c$-causets.

A \textit{path} in a $c$-causet $x$ is a finite sequence $a_1a_2\cdots a_n$ with $a_1\prec a_2\prec\cdots\prec a_n$. The
\textit{height} $h(a)$ of $a\in x$ is the cardinality, minus one, of a longest path in $x$ that ends with $a$. Two vertices $a,b$ are
\textit{comparable} if $a<b$ or $b<a$. It is shown in \cite{gud14} that a causet $x$ is a $c$-causet if and only if $a,b\in x$ are comparable whenever $h(a)\ne h(b)$. We call the set
\begin{equation*}
S_j(x)=\brac{a\in x\colon h(a)=j}
\end{equation*}
the $j$th \textit{shell} of $x$. Letting $s_j(x)=\ab{S_j(x)}$, $j=0,1,2,\ldots ,k$, we call $\paren{s_0(x),s_1(x),\ldots ,s_k(x)}$ the
\textit{shell sequence} of $x$. A $c$-causet is uniquely determined by its shell sequence. Conversely, any finite sequence of positive integers is the shell sequence for a unique $c$-causet.

Let $\omega =a_1a_2\cdots a_n$ be a path in $x$ where $a_j\in{\mathbb N}$ are the labels of the vertices. The \textit{length} of $\omega$ is
\begin{equation*}
{\mathcal L} (\omega )=\sqbrac{\sum _{j=1}^{n-1}(a_{j+1}-a_j)^2}^{1/2}
\end{equation*}
A \textit{geodesic} from $a$ to $b$ where $a<b$, is a path from $a$ to $b$ of smallest length. For $a<b$, we define the \textit{distance from}
$a$ \textit{to} $b$ to be $d(a,b)={\mathcal L} (\omega )$ where $\omega$ is a geodesic from $a$ to $b$. If $a$ and $b$ are incomparable then
$a$ and $b$ are in the same shell and we define $d(a,b)=\ab{\ell (a)-\ell (b)}$. It can be shown that if $a<b<c$ then
$d(a,c)\le d(a,b)+d(b,c)$ \cite{gud14}. Moreover, it is clear that if $a,b,c\in S_j(x)$ for some $j$, then again, $d(a,c)\le d(a,b)+d(b,c)$. Even if we extend the definition of $d$ by defining $d(b,a)=d(a,b)$ when $a<b$, the triangle inequality does not hold in general.

The shell sequence determines the ``shape'' or geometry of a $c$-causet $x$. We view $x$ as a framework or scaffolding of a possible universe. The vertices represent tiny cells that may or may not be occupied by a particle. This geometry gives the kinematics of the system. The dynamics is described in terms of paths in $x$. A free particle tends to move along a geodesic \cite{gud15}, while interacting particles may be forced to propagate along other paths. The remainder of this article describes how this propagation is governed by discrete wave equations. We have already considered gravitational wave equations in \cite{gud15}. Section~2 presents a force-free wave equation and Section~3 studies electric force wave equations. Finally, Section~4 considers the strong nuclear force. Many examples that illustrate this theory are given for the simplest nontrivial $c$-causet.

\section{Difference Operators}
Let $\omega =a_1a_2\cdots a_n$ and $\omega '=b_1b_2\cdots b_n$ be two paths with the same length in a $c$-causet $x$. We say that
$\omega ,\omega '$ are \textit{equal-time} paths if $h(a_1)=h(b_1)$ from which it follows that $h(a_j)=h(b_j)$, $j=1,2,\ldots ,n$.
Two equal-time paths are \textit{noncrossing} if $a_j\ne b_j$, $j=1,\ldots ,n$. The reason that we shall only consider noncrossing paths when they are equal-time is that we need the condition $d(a_j,b_j)\ne 0$. Notice that this condition automatically holds if $\omega ,\omega '$ are not equal-time. If $v$ is a complex-valued function of two variables whose domain includes $\brac{(a_j,b_j)\colon j=1,2,\ldots ,n}$, we define the
\textit{covariant difference operator} on $v$ by
\begin{equation}
\label{eq21}
\nabla _{\omega ,\omega '}v(a_j,b_j)=d(a_{j-1},b_{j-1})v(a_j,b_j)-d(a_j,b_j)v(a_{j-1},b_{j-1})
\end{equation}
$j=2,3,\ldots ,n$.

We only consider $\nabla _{\omega ,\omega '}$ for noncrossing paths when $\omega ,\omega '$ are equal-time. Notice that
$\nabla _{\omega ,\omega '}d(a_j,b_j)=0$, $j=2,3,\ldots ,n$, which is why we call $\nabla _{\omega ,\omega '}$, the
\textit{covariant} difference operator. We think of the function $v$ in \eqref{eq21} as describing amplitudes of a pair of particles moving along the paths $\omega ,\omega '$. We shall usually assume that the two particles are indistinguishable because it makes the theory slightly simpler, but this assumption can be dropped. We call $v$ a \textit{force-free wave function} for $(\omega ,\omega ')$ if
\begin{equation}
\label{eq22}
\nabla _{\omega ,\omega '}v(a_j,b_j)=0
\end{equation}
$j=2,3,\ldots ,n$. We use this terminology because such a $v$ is supposed to describe a pair of particles that are not subject to forces. This includes gravity which is not really a force because gravity is given by the geometry of the $c$-causet. But all particles are subject to gravity so this is actually an idealized approximation in which gravity is neglected. More realistic situations are treated in the next two sections where we include force terms.

\begin{thm}
\label{thm21}
A function $v(a_j,b_j)$, $j=1,2,\ldots ,n$, satisfies \eqref{eq22} if and only if $v(a_j,b_j)=cd(a_j,b_j)$, $j=2,3,\ldots ,n$ for some fixed
$c\in{\mathbb C}$.
\end{thm}
\begin{proof}
If $v$ satisfies \eqref{eq22} then by \eqref{eq21} we have
\begin{equation*}
v(a_j,b_j)=\frac{d(a_j,b_j)}{d(a_{j-1},b_{j-1})}\,v(a_{j-1},b_{j-1})\quad j=2,3,\ldots ,n
\end{equation*}
If $j\ge 3$ we can continue to obtain
\begin{align*}
v(a_j,b_j)&=\frac{d(a_j,b_j)}{d(a_{j-1},b_{j-1})}\,\frac{d(a_{j-1},b_{j-1})}{d(a_{j-2},b_{j-2})}\,v(a_{j-2},b_{j-2})\\\noalign{\medskip}
  &=\frac{d(a_j,b_j)}{d(a_{j-2},b_{j-2})}\,v(a_{j-2},b_{j-2})
\end{align*}
Continue this process to obtain
\begin{equation*}
v(a_j,b_j)=\frac{d(a_j,b_j)}{d(a_1,b_1)}\,v(a_1,b_1)\end{equation*}
for $j=1,2,\ldots ,n$. Letting $c=v(a_1,b_1)/d(a_1,b_1)$ gives our result. The converse clearly holds.
\end{proof}

We shall mainly consider equal-time noncrossing paths $\omega ,\omega '$ in the sequel because they give a simpler quantum theory. Moreover, in this case $d(a_j,b_j)=\ab{a_j-b_j}$, where we represent a vertex by its label, which is a further simplification. In this situation, let $v$ be a force-free wave function for $(\omega ,\omega ')$ and assume that $v(a_1,b_1)=d(a_1,b_1)$. Then $c=1$ in Theorem~\ref{thm21} so that
\begin{equation*}
v(a_j,b_j)=d(a_j,b_j)=\ab{a_j-b_j}
\end{equation*}
We do not lose generality with this assumption because we shall later normalize our wave functions so the theory is independent of a multiplicative constant. Since $v$ depends on $\omega$ and $\omega '$ we sometimes write
\begin{equation*}
v(a_n,b_n)=v(\omega ,\omega ',a_n,b_n)
\end{equation*}
and think of $v$ as the amplitude that the particles move along $\omega$ and $\omega '$ starting at $(a_1,b_1)$ and ending at $(a_n,b_n)$. According to the usual quantum formalism, the amplitude that the particles move to positions $(a_n,b_n)$ becomes
\begin{equation}
\label{eq23}
\widehat{v} (a_n,b_n)=\sum _{\omega ,\omega '}v(\omega ,\omega ',a_n,b_n)=\sum _{\omega ,\omega '}\ab{a_n-b_n}
\end{equation}
where $(\omega ,\omega ')$ are all pairs of noncrossing paths from $(a_1,b_1)$ to $(a_n,b_n)$. Letting $\tau$ be the number of such path pairs \eqref{eq23} gives
\begin{equation}
\label{eq24}
\widehat{v} (a_n,b_n)=\tau\ab{a_n-b_n}
\end{equation}

Suppose $h(a_1)=h(b_1)=p$ and $h(a_n)=h(b_n)=q$, $p<q$. We then have that
\begin{equation}
\label{eq25}
\tau =s_{p+1}(x)\sqbrac{s_{p+1}(x)-1}s_{p+2}(x)\sqbrac{s_{p+2}(x)-1}\cdots s_q(x)\sqbrac{s_q(x)-1}
\end{equation}
By \eqref{eq24} and \eqref{eq25}, the position amplitude $\widehat{v} (a_n,b_n)$ does not depend on $a_1$ and $b_1$ and only depends on
$a_n$ and $b_n$ through the term $\ab{a_n-b_n}$. We next normalize $\widehat{v} (a_n,b_n)$ to obtain the \textit{position state} $\psi (a_n,b_n)$. Thus, if $S_q(x)=\brac{c_1,c_2,\ldots ,c_r}$ where the $c_j$ are vertex labels in increasing order, then the normalization constant becomes:
\begin{equation*}
N^2=\sum\brac{\ab{\widehat{v} (c_j,c_k)}^2\colon c_j<c_k}
\end{equation*}
and $\psi (c_j,c_k)=\widehat{v} (c_j,c_k)/N$. Since $\tau$ is a constant that appears in every term of $N$ and also in $\widehat{v} (c_j,c_k)$ it cancels out and we can omit it. We then obtain the normalization constant
\begin{equation*}
N_1^2=\sum _{j=1}^{r-1}(r-j)j^2
\end{equation*}
and $\psi (c_j,c_k)=\ab{j-k}/N_1$.

As usual in the quantum formalism, we define the probability that the two particles end at positions $(c_j,c_k)$ to be
\begin{equation*}
p(c_j,c_k)=\ab{\psi (c_j,c_k)}^2=\frac{\ab{j-k}^2}{N_1^2}
\end{equation*}
We call $p(c_j,c_k)$ the \textit{force-free distribution}. The constant $N_1$ is not so important. What is important is the relative sizes of the probabilities $p(c_j,c_k)$. For example, if $c_k-c_j=r-1$ and $c_t-c_s=1$, then
\begin{equation*}
\frac{p(c_s,c_t)}{p(c_j,c_k)}=\frac{1}{(r-1)^2}
\end{equation*}
which can be quite small. However, we must remember that there are $r-1$ vertices with $c_t-c_s=1$ and only one vertex with $c_k-c_j=r-1$ so letting $A$ be the set of vertices satisfying $c_t-c_s=1$ and $B$ be the set satisfying $c_k-c_j=r-1$ we have
\begin{equation*}
\frac{\mathrm{prob} (A)}{\mathrm{prob} (B)}=\frac{1}{r-1}
\end{equation*}
This can still be small but it is not as small as $1/(r-1)^2$.

\begin{exam}{1}
Suppose $S_q=\brac{c_1,c_2,c_3,c_4,c_5}$. We then have
\begin{equation*}
N_1^2=\sum _{j=1}^4(5-j)j^2=4+3\mathrel{\mathlarger\cdot} 2^2+2\mathrel{\mathlarger\cdot} 3^2+4^2=50
\end{equation*}
The probabilities become:
\begin{align*}
p(c_1,c_2)&=p(c_2,c_3)=p(c_3,c_4)=p(c_4,c_5)=\tfrac{1}{50}\\\noalign{\medskip}
p(c_1,c_3)&=p(c_2,c_4)=p(c_3,c_5)=\tfrac{4}{50}=\tfrac{2}{25}\\\noalign{\medskip}
p(c_1,c_4)&=p(c_2,c_5)=\tfrac{9}{50}\\\noalign{\medskip}
p(c_1,c_5)&=\tfrac{16}{50}=\tfrac{8}{25}
\end{align*}
We conclude that being three apart is more probable than any other configuration.\qquad\qedsymbol
\end{exam}

\section{Electric Forces}
We now discuss an electric force wave equation. As in Section~2, let $\omega =a_1a_2\cdots a_n$, $\omega '=b_1b_2\cdots b_n$ be two paths in $x$. If $\omega ,\omega '$ are equal-time paths, we assume they are noncrossing. Suppose we have two charged particles with electric charges $e_1$ and $e_2$. We again let $v$ be a complex-valued function of two variables whose domain includes
$\brac{(a_j,b_j)\colon j=1,2,\ldots ,n}$. If $e_1e_2\ge 0$, the \textit{repulsive electric-force wave equation} is defined to be
\begin{equation}
\label{eq31}
i\nabla _{\omega ,\omega '}v(a_j,b_j)=\frac{e_1e_2}{d(a_j,b_j)}\,v(a_j,b_j)
\end{equation}
$j=2,3,\ldots ,n$. If $e_1e_2\le 0$, the \textit{attractive electric-force wave equation} is defined to be
\begin{equation}
\label{eq32}
i\nabla _{\omega ,\omega '}v(a_j,b_j)=\frac{e_1e_2}{d(a_j,b_j)}\,v(a_{j-1},b_{j-1})
\end{equation}
$j=2,3,\ldots ,n$. Of course, if $e_1=0$ or $e_2=0$, then \eqref{eq31} and \eqref{eq32} reduce to the force-free wave equation \eqref{eq22}. If $v$ satisfies \eqref{eq31} or \eqref{eq32} for $j=2,3,\ldots ,n$, then $v$ is called an \textit{electric force wave function}. The next result shows that \eqref{eq31} or \eqref{eq32} and the \textit{initial condition} $v(a_1,b_1)$ uniquely determine $v$.

\begin{thm}
\label{thm31}
A function $v(a_j,b_j)$, $j=1,2,\ldots ,n$, satisfies \eqref{eq31} if and only if
\begin{equation}
\label{eq33}
v(a_j,b_j)=\frac{i^{j_1}}{\sqbrac{\frac{id(a_{j-1},b_{j-1})}{d(a_j,b_j)}\,-\,\frac{e_1e_2}{d(a_j,b_j)^2}}\cdots
 \sqbrac{\frac{id(a_1,b_1)}{d(a_2,b_2)}\,-\,\frac{e_1e_2}{d(a_2,b_2)^2}}}\,v(a_1,b_1)
\end{equation}
and $v(a_j,b_j)$, $j=1,2,\ldots ,n$, satisfies \eqref{eq32} if and only if
\begin{equation}
\label{eq34}
v(a_j,b_j)=\frac{d(a_j,b_j)^2-ie_1e_2}{d(a_{j-1},b_{j-1})d(a_j,b_j)}\cdots
 \frac{d(a_2,b_2)^2-ie_1e_2}{d(a_1,b_1)d(a_2,b_2)}\,v(a_1,b_1)
\end{equation}
\end{thm}
\begin{proof}
If $v(a_j,b_j)$ satisfies \eqref{eq31}, we have
\begin{align*}
i\sqbrac{d(a_{j-1},b_{j-1}v(a_j,b_j)-d(a_j,b_j)v(a_{j-1},b_{j-1}}=\frac{e_1e_2}{d(a_j,b_j)}\,v(a_j,b_j)\\
\intertext{Hence,}
v(a_j,b_j)=\frac{id(a_j,b_j)}{id(a_{j-1},b_{j-1})-\frac{e_1e_2}{d(a_j,b_j)}}\,v(a_{j-1},b_{j-1})
\end{align*}
Continuing this process we obtain
\begin{align*}
v(a_j,b_j)&=\frac{id(a_j,b_j)}{id(a_{j-1},b_{j-1})-\frac{e_1e_2}{d(a_j,b_j)}}
  \mathrel{\mathlarger\cdot}\frac{id(a_{j-1},b_{j-1})}{id(a_{j-2},b_{j-2})-\frac{e_1e_2}{d(a_{j-1},b_{j-1})}}\,v(a_{j-2},b_{j-2})\\
  &\ \vdots\\
  &=\frac{id(a_j,b_j)}{id(a_{j-1},b_{j-1})-\frac{e_1e_2}{d(a_j,b_j)}}\cdots\frac{id(a_2,b_2)}{id(a_1,b_1)-\frac{e_1e_2}{d(a_2,b_2)}}
  \,v(a_1,b_1)
\end{align*}
from which the result follows. If $v(a_j,b_j)$ satisfies \eqref{eq32} we have
\begin{equation*}
i\sqbrac{d(a_{j-1},b_{j-1}v(a_j,b_j)-d(a_j,b_j)v(a_{j-1},b_{j-1}}=\frac{e_1e_2}{d(a_j,b_j)}\,v(a_{j-1},b_{j-1})
\end{equation*}
Hence,
\begin{align*}
v(a_j,b_j)&=\frac{id(a_j,b_j)^2+e_1e_2}{id(a_{j-1},b_{j-1})d(a_j,b_j}\,v(a_{j-1},b_{j-1})\\\noalign{\medskip}
  &=\frac{d(a_j,b_j)^2-ie_1e_2}{d(a_{j-1},b_{j-1})d(a_j,b_j)}\,v(a_{j-1},b_{j-1})\\
    &\ \vdots\\
  &=\frac{d(a_j,b_j)^2-ie_1e_2}{d(a_{j-1},b_{j-1})d(a_j,b_j)}\cdots\frac{d(a_2,b_2)^2-ie_1e_2}{id(a_1,b_1)d(a_2,b_2)}\,v(a_1,b_1)
\end{align*}
The converses are straightforward.
\end{proof}

\begin{exam}{2}
One of the simplest nontrivial $c$-causets has shell sequence $(1,2,2,3)$. Delineating shells by semicolons, we can label the vertices as $(1;2,3;4,5;6,7,8)$. We now examine two indistinguishable particles both having charge $e$ that move from the first shell $(2,3)$ to the top shell $(6,7,8)$. The paths between these shells are:
$\omega _1=2-4-6$, $\omega _2=3-4-6$, $\omega _3=2-5-6$, $\omega _4=3-5-6$, $\omega _5=2-4-7$, $\omega _6=2-5-7$,
$\omega _7=3-4-7$, $\omega _8=3-5-7$, $\omega _9=2-4-8$, $\omega _{10}=2-5-8$, $\omega _{11}=3-4-8$, $\omega _{12}=3-5-8$. Applying \eqref{eq33} we obtain
\begin{align*}
v(\omega _1,\omega _8,6,7)&=v(\omega _2,\omega _6,6,7)=v(\omega _3,\omega _7,6,7)=v(\omega _4,\omega _5,6,7)\\
  &=-\frac{1}{(i-e^2)^2}=\frac{(1-e^4)-2e^2i}{(1-e^4)^2+4e^4}
\end{align*}
Moreover,
\begin{align*}
v(\omega _1,\omega _{12},6,8)&=v(\omega _1,\omega _{10},6,8)=v(\omega _3,\omega _{11},6,8)=v(\omega _8,\omega _9,6,8)\\
  &=-\frac{-1}{\paren{\frac{i}{2}-\frac{e^2}{4}}(i-e^2)}=\frac{4(2-e^4)-12e^2i}{(2-e^4)^2+9e^4}
\end{align*}
Summing over noncrossing path pairs as in \eqref{eq23} we obtain the position wave functions
\begin{align*}
\widehat{v} (6,7)&=4\sqbrac{\frac{(1-e^4)-2e^2i}{(1-e^4)^2+4e^4}}\\\noalign{\medskip}
\widehat{v} (6,8)&=4\sqbrac{\frac{4(2-e^4)-12e^2i}{(2-e^4)^2+9e^4}}
\end{align*}
and by symmetry $\widehat{v} (7,8)=\widehat{v} (6,7)$. The normalization constant becomes
\begin{equation*}
N^2=2\ab{\widehat{v} (6,7)}^2+\ab{\widehat{v} (6,8)}^2=\frac{32}{(1+e^4)^2}+\frac{256}{(1+e^4)(4+e^4)}
\end{equation*}
We compute the probabilities as
\begin{align*}
p(6,7)=p(7,8)=\frac{\ab{\widehat{v} (6,7)}^2}{N^2}=\frac{4+e^4}{6(4+3e^4)}\\
\intertext{and}
p(6,8)=1-2p(6,7)=\frac{8(1+e^4)}{3(4+3e^4)}
\end{align*}
Now $p(6,7)$ and $p(6,8)$ are functions of $e$ and we write $p(6,7)=p(6,7)(e)$ and $p(6,8)=p(6,8)(e)$. Notice that $p(6,7)(0)=1/6$ and
$p(6,8)(0)=2/3$ which agree with the force-free probabilities of Section ~2. Moreover, for large $e$ we have
\begin{equation*}
\lim _{e\to\infty}p(6,7)(e)=\tfrac{1}{18},\quad\lim _{e\to\infty}p(6,8)(e)=\tfrac{8}{9}
\end{equation*}
Thus, under this repulsive electric force, the particles prefer to be farther apart. The function $p(6,7)(e)$ decreases monotonically from $1/6$ at $e=0$ to $1/18$ at $e\approx 3$. The function $p(6,8)(e)$ increases monotonically from $2/3$ at $e=0$ to $8/9$ at $e\approx 3$.
\qquad\qedsymbol
\end{exam}

Although we shall not use this in the sequel, we now briefly discuss how to place this theory within the tensor product formalism. In Example~2, we think of $S_3(x)=\brac{6,7,8}$ as a 3-dimensional Hilbert space with orthonormal basis
$\widehat{6} ,\widehat{7} ,\widehat{8}$. We can extend the states $\psi (j,k)=\widehat{v} (j,k)/N$ by defining $\psi (k,j)=-\psi(j,k)$. In fact, we see from
\eqref{eq33} and \eqref{eq34} that if $v(a_1,b_1)=-v(b_1,a_1)$ then $v(a_j,b_j)$ and hence $\psi (a_j,b_j)$ are automatically antisymmetric. (We could also make them symmetric in this way.) If we do this antisymmetrization, the state we have constructed in Example~2 can be written as
\begin{align*}
\psi&=\frac{1}{\sqrt{2}}\left[\psi (6,7)\widehat{6}\otimes\widehat{7} +\psi (7,6)\widehat{7}\otimes\widehat{6} +\psi (7,8)\widehat{7}\otimes\widehat{8}
  +\psi(8,7)\widehat{8}\otimes\widehat{7} \right.\\
  &\quad\left.+\psi (6,8)\widehat{6}\otimes\widehat{8} +\phi (8,6)\widehat{8}\otimes\widehat{6}\right]\\
  &=\psi (6,7)\frac{(\,\widehat{6}\otimes\widehat{7} -\widehat{7}\otimes\widehat{6}\,)}{\sqrt{2}}
  +\psi (7,8)\frac{(\,\widehat{7}\otimes\widehat{8} -\widehat{8}\otimes\widehat{7}\,)}{\sqrt{2}}
  +\psi (6,8)\frac{(\,\widehat{6}\otimes\widehat{8} -\widehat{8}\otimes\widehat{6}\,)}{\sqrt{2}}
\end{align*}
Then the probability that the particle starting at vertex~2 ends at vertex~6 while the particle starting at vertex~3 ends at vertex~7 is
\begin{equation*}
\ab{\elbows{\psi ,\widehat{6}\otimes\widehat{7}}}^2=\tfrac{1}{2}\ab{\psi (6,7)}^2
\end{equation*}
Similarly, if we reverse vertices 2 and 3 in the above we obtain
\begin{equation*}
\ab{\elbows{\psi ,\widehat{7}\otimes\widehat{6}}}^2=\tfrac{1}{2}\ab{\psi (6,7)}^2
\end{equation*}
Similar results hold for vertex pairs $(7,8)$ and $(6,8)$.

\begin{exam}{3}
We can adjoin a vertex labeled 9 to the top shell in Example~2 to obtain the shell sequence $(1,2,2,4)$. This $c$-causet can be analyzed by constructing the new paths $\omega _{13}=2-4-9$, $\omega _{14}=2-5-9$, $\omega _{15}=3-4-9$, $\omega _{16}=3-5-9$. We then have
\begin{align*}
v(\omega _1,\omega _{16},6,9)&=v(\omega _2,\omega _{14},6,9)=v(\omega _3,\omega _{15},6,9)=v(\omega _4,\omega _{13},6,9)\\
  &=-\frac{-1}{\paren{\frac{i}{3}-\frac{e^2}{9}}(i-e^2)}=\frac{9(3-e^4)-36e^2i}{(3-e^4)^2+16e^4}
\end{align*}
This gives
\begin{align*}
\widehat{v} (6,9)=36\sqbrac{\frac{(3-e^4)-4e^2i}{(3-e^4)^2+16e^4}}\\
\intertext{Hence,}
\ab{\widehat{v} (6,9)}^2=\frac{36^2}{(9+e^4)(1+e^4)}
\end{align*}
Using the values from Example~2 we obtain
\begin{align*}
p(6,7)&=\frac{\ab{\widehat{v} (6,7)}^2}{3\ab{\widehat{v} (6,7)}^2+2\ab{\widehat{v} (6,8)}^2+\ab{\widehat{v} (6,9)}^2}\\
&=\frac{1}{3+2\frac{\ab{\widehat{v} (6,8)}^2}{\ab{\widehat{v} (6,7)}^2}+\frac{\ab{\widehat{v} (6,9)}^2}{\ab{\widehat{v} (6,7)}^2}}\\
&=\frac{1}{3+\frac{32(1+e^4)}{4+e^4}+\frac{81(1+e^4)}{9+e^4}}
\end{align*}
In a similar way we obtain
\begin{align*}
p(6,8)&=\frac{1}{2+\frac{3(4+e^4)}{16(1+e^4)}+\frac{81(4+e^4)}{16(9+e^4)}}\\
p(6,9)&=\frac{1}{1+\frac{9+e^4}{27(1+e^4)}+\frac{32(9+e^4)}{81(4+e^4)}}
\end{align*}
As in Section~2 we obtain the force-free distribution
\begin{align*}
p(6,7)(0)&=p(7,8)(0)=p(8,9)(0)=\tfrac{1}{20}\\
p(6,8)(0)&=p(7,9)(0)=\tfrac{1}{5}\\
p(6,9)(0)&=\tfrac{9}{20}
\end{align*}
For large values of $e$ (even for $e\approx 3$) we obtain
\begin{equation*}
\lim _{e\to\infty}p(6,7)(e)=\tfrac{1}{116},\quad\lim _{e\to\infty}p(6,8)(e)=\tfrac{16}{116},\quad\lim _{e\to\infty}p(6,9)(e)=\tfrac{81}{116}
\end{equation*}
The functions $p(6,7)(e)$ and $p(6,8)(e)$ decrease monotonically, while $P(6,9)(e)$ increases monotonically. Again, the two particles prefer to move farther apart. This phenomenon becomes more pronounced as more vertices are adjoined to the top shell.\qquad\qedsymbol
\end{exam}

\begin{exam}{4}
This example considers the attractive wave equation \eqref{eq32} with $e_1=-e_2=e$ for the $c$-causet of Example~2. Applying \eqref{eq34} we have
\begin{align*}
v(\omega _1,\omega _8,6,7)&=v(\omega _2,\omega _6,6,7)=v(\omega _3,\omega _7,6,7)=v(\omega _4,\omega _5,6,7)\\
&=(1+ie^2)^2\\
v(\omega _1,\omega _{12},6,8)&=v(\omega _2,\omega _{10},6,8)=v(\omega _3,\omega _{11},6,8)=v(\omega _4,\omega _9,6,8)\\
&=\tfrac{1}{2}(4+ie^2)(1+ie^2)
\end{align*}
We conclude that
\begin{align*}
\widehat{v} (6,7)&=\widehat{v} (7,8)=4(1-e^4)+8e^2i\\
\widehat{v} (6,8)&=2(4-e^4)+10e^2i
\end{align*}
Hence,
\begin{align*}
\ab{\widehat{v} (6,7)}^2&=\ab{\widehat{v} (7,8)}^2=16(1+e^4)^2\\
\ab{\widehat{v} (6,8)}^2&=4(1+e^4)(16+e^4)
\end{align*}
Computing the probabilities gives
\begin{align*}
p(6,7)&=p(7,8)=\frac{\ab{\widehat{v} (6,7)}^2}{2\ab{\widehat{v} (6,7)}^2+\ab{\widehat{v} (6,8)}^2}=\frac{1}{2+\frac{\ab{\widehat{v} (6,8)}^2}{\ab{\widehat{v} (6,7)}^2}}\\
  &=\frac{4(1+e^4)}{3(8+3e^4)}\\
  p(6,8)&=\frac{1}{1+\frac{2\ab{\widehat{v} (6,7)}^2}{\ab{\widehat{v} (6,8)}^2}}=\frac{16+e^4}{3(8+3e^4)}
\end{align*}
Of course, we also have that $p(6,8)=1-2(p(6,7)$. As before, we obtain
\begin{equation*}
p(6,7)(0)=p(7,8)(0)=1/16,\quad p(6,8)(0)=2/3
\end{equation*}
However, in contrast with Examples~2 and 3 we have
\begin{align*}
\lim _{e\to\infty}p(6,7)(e)&=\lim _{e\to\infty}p(7,8)(e)=\tfrac{4}{9}\\
\lim _{e\to\infty}p(6,8)(e)&=\tfrac{1}{9}
\end{align*}
Thus, the particles prefer to be close together under the attractive force. Notice that
\begin{equation*}
p(6,7)(1)=p(7,8)(1)=\tfrac{8}{33},\quad p(6,8)(1)=\tfrac{17}{33}
\end{equation*}
so that for $e\approx 1$, the probabilities that the two particles are within one unit and two units are about the same.\qquad\qedsymbol
\end{exam}

\begin{exam}{5}
We now compute probabilities for the $c$-causet in Example~3 with an attractive electric force as in Example~4. We again have the new paths $\omega _{13},\omega _{14},\omega _{15},\omega _{16}$ for which we obtain
\begin{align*}
v(\omega _1,\omega _{16},6,9)&=v(\omega _2,\omega _{14},6,9)=v(\omega _3,\omega _{15},6,9)=v(\omega _4,\omega _{13},6,9)\\
  &=\frac{9+ie^2}{3}\,(1+ie^2)
\end{align*}
This gives
\begin{align*}
\widehat{v} (6,9)=\tfrac{4}{3}(9-e^4+10e^2i)\\
\intertext{Hence,}
\ab{\widehat{v} (6,9)}^2=\tfrac{16}{9}\,(1+e^4)(81+e^4)
\end{align*}
Using the values from Example~4 we obtain
\begin{align*}
p(6,7)&=\frac{1}{3+\frac{2\ab{\widehat{v} (6,8)}^2}{\ab{\widehat{v} (6,7)}^2}+\frac{\ab{\widehat{v} (6,9)}^2}{\ab{\widehat{v} (6,7)}^2}}
  =\frac{1}{3+\frac{16+e^4}{2(1+e^4)}+\frac{81+e^4}{9(1+e^4)}}\\
  &=\frac{18(1+e^4)}{5(72+13e^4)}
\end{align*}
In a similar way we have
\begin{align*}
p(6,8)&=\frac{9(16+e^4)}{10(72+13e^4)}\\
p(6,9)&=\frac{2(81+e^4)}{5(72+13e^4)}
\end{align*}
As before
\begin{equation*}
p(6,7)(0)=\tfrac{1}{20},\quad p(6,8)(0)=\tfrac{1}{5},\quad p(6,9)(0)=\frac{9}{20}
\end{equation*}
But for large $e$ we have
\begin{equation*}
\lim _{e\to\infty}p(6,7)(e)=\tfrac{18}{65},\quad\lim _{e\to\infty}p(6,8)(e)=\tfrac{9}{130},\quad\lim _{e\to\infty}p(6,9)(e)=\tfrac{2}{65}
\end{equation*}
Again, for $e\approx 3$ and larger, the particles prefer to be together.\qquad\qedsymbol
\end{exam}

\section{Strong Nuclear Force}
We propose that the strong nuclear force acts like an attractive spring force with wave equation given by
\begin{equation}
\label{eq41}
i\nabla _{\omega ,\omega '}v(a_j,b_j)=-gd(a_{j-1},b_{j-1})v(a_{j-1},b_{j-1})
\end{equation}
where $g>0$ is a constant. As in Theorem~\ref{thm31}, a function $v(a_j,b_j)$ satisfies \eqref{eq41} for $j=2,3,\ldots ,n$ if and only if $v$ as the form
\begin{equation}
\label{eq42}
v(a_j,b_j)=\sqbrac{\frac{d(a_j,b_j)}{d(a_{j-1},b_{j-1})}+ig}\cdots\sqbrac{\frac{d(a_2,b_2)}{d(a_1,b_1)}+ig}v(a_1,b_1)
\end{equation}
We now repeat Examples~2 and 3 for the strong nuclear force.

\begin{exam}{6}
As in Example~2, we consider the $c$-causet with shell sequence $(1,2,2,3)$. Applying \eqref{eq42} we obtain
\begin{align*}
v(\omega _1,\omega _8,6,7)&=v(\omega _2,\omega _6,6,7)=v(\omega _3,\omega _7,6,7)=v(\omega _4,\omega _5,6,7)\\
  &=(1+ig)^2\\
v(\omega _1,\omega _{12},6,8)&=v(\omega _1,\omega _{10},6,8)=v(\omega _3,\omega _{11},6,8)=v(\omega _8,\omega _9,6,8)\\
  &=(2+ig)(1+ig)
\end{align*}
Hence,
\begin{align*}
\widehat{v} (6,7)&=4(1+ig)^2=4(1-g^2)+8gi\\
\widehat{v} (6,8)&=4(2+ig)(1+ig)=4(2-g^2)+12gi
\end{align*}
It follows that
\begin{equation*}
\ab{\widehat{v} (6,7)}^2=16(1+g^2)^2,\quad\ab{\widehat{v} (6,8)}^2=16(1+g^2)(4+g^2)
\end{equation*}
We then obtain
\begin{align*}
p(6,7)&=p(7,8)=\frac{1}{2+\frac{\ab{\widehat{v} (6,8)}^2}{\ab{\widehat{v} (6,7)}^2}}=\frac{1}{2+\frac{4+g^2}{1+g^2}}=\frac{1+g^2}{3(2+g^2)}\\
p(6,8)&=1-2p(6,7)=\frac{4+g^2}{3(2+g^2)}
\end{align*}
For the larger $g$ we have
\begin{equation*}
\lim _{g\to\infty}p(6,7)(g)=\lim _{g\to\infty}p(7,8)(g)=\lim _{g\to\infty}p(6,8)(g)=\tfrac{1}{3}
\end{equation*}
Since $g$ is a constant, we should not let $g$ vary. However, if $g$ about 3 or higher, then $p(6,7)$, $p(7,8)$ and $p(6,8)$ are all close to $1/3$.\qquad\qedsymbol
\end{exam}

\begin{exam}{7}
We next consider the strong nuclear force for the $c$-causet with shell sequence $(1,2,2,4)$. As in Example~3 we have the new paths
$\omega _{13},\omega _{14}, \omega _{15},\omega _{16}$. Each of the noncrossing path pairs to $(6,9)$ have the wave function
\begin{align*}
v(\omega _1,\omega _{16},6,9)&=(3+ig)(1+ig)\\
\intertext{Hence,}
\widehat{v} (6,9)&=4(3+ig)(1+ig)=4(3-g^2)+16gi\\
\intertext{and}
\ab{\widehat{v} (6,9)}^2&=16(1+g^2)(9+g^2)
\end{align*}
We conclude that
\begin{align*}
p(6,7)&=\frac{1}{3+\frac{2\ab{\widehat{v} (6,8)}^2}{\ab{\widehat{v} (6,7)}^2}+\frac{\ab{\widehat{v} (6,9)}^2}{\ab{\widehat{v} (6,7)}^2}}
  =\frac{1}{3+\frac{2(4+g^2)}{1+g^2}+\frac{9+g^2}{1+g^2}}\\
  &=\frac{1+g^2}{2(10+3g^2)}
\end{align*}
In a similar way
\begin{equation*}
p(6,8)=\frac{4+g^2}{2(10+3g^2)}\,,\quad p(6,9)=\frac{9+g^2}{2(10+3g^2)}
\end{equation*}
We then obtain
\begin{equation*}
\lim _{g\to\infty}p(6,7)(g)=\lim _{g\to\infty}p(6,8)(g)=\lim _{g\to\infty}p(6,9)(g)=\tfrac{1}{6}
\end{equation*}
which is a bit surprising.\qquad\qedsymbol
\end{exam}

\begin{exam}{8}
The previous example compels us to consider the strong nuclear force for the $c$-causet with shell sequence $(1,2,2,5)$ to see if there is a pattern. We now have the new paths $\omega _{17}=2-4-10$, $\omega _{18}=2-5-10$, $\omega _{19}=3-4-10$,
$\omega _{20}=3-5-10$. Each of the noncrossing path pairs to $(6,10)$ have the wave function
\begin{equation*}
v(\omega _2,\omega _{17},6,10)=(4+ig)(1+ig)
\end{equation*}
Hence,
\begin{align*}
\widehat{v} (6,10)&=4(4+ig)(1+ig)=4(4-g^2)+20gi\\
\intertext{and}
\ab{\widehat{v} (6,10)}^2&=16(1+g^2)(16+g^2)
\end{align*}
We conclude that
\begin{align*}
p(6,7)&=\frac{\ab{\widehat{v} (6,7)}^2}{4\ab{\widehat{v} (6,7)}^2+3\ab{\widehat{v} (6,8)}^2+2\ab{\widehat{v} (6,9)}^2+\ab{\widehat{v} (6,10)}^2}\\\noalign{\medskip}
  &=\frac{1}
  {4+\frac{3\ab{\widehat{v} (6,8)}^2}{\ab{\widehat{v} (6,7)}^2}+\frac{2\ab{\widehat{v} (6,9)}^2}{\ab{\widehat{v} (6,7)}^2}+\frac{\ab{\widehat{v} (6,10)}^2}{\ab{\widehat{v} (6,7)}^2}}
  \\\noalign{\medskip}
  &=\frac{1}{4+\frac{3(4+g^2)}{1+g^2}+\frac{2(9+g^2)}{1+g^2}+\frac{16+g^2}{1+g^2}}=\frac{1+g^2}{10(5+g^2)}
\end{align*}
In a similar way we obtain
\begin{equation*}
p(6,8)=\frac{4+g^2}{10(5+g^2)}\,,\quad p(6,9)=\frac{9+g^2}{10(5+g^2)}\,,\quad p(6,10)=\frac{16+g^2}{10(5+g^2)}
\end{equation*}
We then have
\begin{equation*}
p(6,7)(0)=\tfrac{1}{50},\quad p(6,8)(0)=\tfrac{2}{25},\quad p(6,9)(0)=\tfrac{9}{50},\quad (6,10)(0)=\tfrac{8}{25}
\end{equation*}
as predicted in Section~2. Moreover, we obtain
\begin{equation*}
\lim _{g\to\infty}p(6,7)(g)=\lim _{g\to\infty}p(6,8)(g)=\lim _{g\to\infty}p(6,9)(g)=\lim _{g\to\infty}p(6,10)(g)=\tfrac{1}{10}
\end{equation*}
It is not hard to show that this pattern continues and for a $c$-causet with shell sequence $(1,2,2,n)$ we have
\begin{equation*}
\lim _{g\to\infty}p(6,7)(g)=\lim _{g\to\infty}p(6,8)(g)=\cdots =\lim _{g\to\infty}p(6,n+5)(g)=\tfrac{n(n-1)}{2}
\end{equation*}
We suspect that this strong regularity is due to the fact that we are considering small $c$-causets. For larger $c$-causets we believe that the situation will be much more complex.\qquad\qedsymbol
\end{exam}


\section*{Introduction}

The Chow variety, introduced in 1937 by Chow and van der Waerden~\cite{CvW},
parameterizes algebraic cycles of any fixed dimension and degree in a
projective space,
each given by its Chow form.
 The case of curves in~$\bbbp^3$ goes back to an 1848 paper by
Cayley~\cite{cayley}.  A fundamental problem, addressed by Green and
Morrison~\cite{green_morrison} as well as Gel'fand, Kapranov and
Zelevinsky~\cite[\S4.3]{gkz}, is to describe the equations defining Chow
varieties.  We present a definitive computational solution for the smallest
non-trivial case, namely for cycles of dimension~$1$ and degree~$2$
in~$\bbbp^3$.

The \emph{Chow form} of a  cycle of degree $2$ is a quadratic form in the
Pl\"ucker coordinates of the Grassmannian~$G(2,4)$ of lines
in~$\bbbp^3$. Such a quadric in~$G(2,4)$  represents the set
 of all lines that intersect the given cycle. Quadratic forms in Pl\"ucker
coordinates form a projective space~$\bbbp^{19}$. The Chow variety we are
interested in, denoted~$G(2,2,4)$, is the set of all Chow forms in
that~$\bbbp^{19}$.  The aim of this note is to make  the concepts
in~\cite{cayley,CvW,green_morrison} and \cite[\S 4.3]{gkz} completely explicit.

We start with the~$9$-dimensional subvariety of~$\bbbp^{19}$ whose
points are the {\em coisotropic quadrics} in~$G(2,4)$.  By \cite[\S4.3, Theorem~3.14]{gkz},
this  decomposes as the Chow variety and the variety of Hurwitz
forms~\cite{Stu}, representing lines that are tangent to a
quadric surface in~$\bbbp^3$.  Section~\ref{sec:coiso}
studies the ideal generated by the coisotropy conditions.
We work in a polynomial ring in~$20$ variables,
one for each quadratic Pl\"ucker monomial on $G(2,4)$
minus one for the Pl\"ucker relation.
We derive the coisotropic ideal from the differential characterization of
coisotropy. Proposition~\ref{prop:coisotropic} exhibits the decomposition
of this ideal into three minimal primes.
In particular, this shows that the coisotropic ideal is radical,
and  it hence resolves the degree $2$ case of
 a problem posed in 1986 by Green and Morrison \cite{green_morrison}.
 They wrote:  `We do not know whether [the differential characterization of
coisotropy] generates the full ideal of these Chow variables.'

Section~\ref{sec:chowvar} derives the radical ideal of
the Chow variety $G(2,2,4)$ in~$\bbbp^{19}$.
Its two minimal primes represent Chow forms of plane conics and Chow forms of pairs of
lines. We also study the characterization of Chow forms among all coisotropic
quadrics by the vanishing of certain differential forms.
These represent the integrability of the $\alpha$-distribution
in \cite[\S4.3, Theorem 3.22]{gkz}.
After saturation by the irrelevant ideal, the integrability ideal is found to be radical.

\section{Coisotropic Quadrics}
\label{sec:coiso}

The Grassmannian $G(2,4)$ is a quadric in $\bbbp^5$.
Its points are lines in $\bbbp^3$. We represent these lines using
dual Pl\"ucker coordinates $\vec{p} = (p_{01}, p_{02}, p_{03}, p_{12}, p_{13}, p_{23})$
subject to the Pl\"ucker relation~$p_{01}p_{23}-p_{02}p_{13} + p_{03}p_{12}$.
Following \cite[\S 2]{Stu}, by
{\em dual} coordinates we mean that
$p_{ij}$ is the $ij$-minor of a $2 \times 4$-matrix
whose rows span the line.
The generic quadric in $G(2,4)$   %
is written as a generic quadratic form
\begin{equation}
\label{general_quadric}
Q(\vec{p}) \ =\ \,  \vec{p} \cdot \left(
\begin{array}{cccccc}
c_0 & c_1 & c_2 & c_3 & c_4 & c_5 \\
c_1 & c_6 & c_7 & c_8 & c_9 & c_{10} \\
c_2 & c_7 & c_{11} & c_{12} & c_{13} & c_{14} \\
c_3 & c_8 & c_{12} & c_{15} & c_{16} & c_{17} \\
c_4 & c_9 & c_{13} & c_{16} & c_{18} & c_{19} \\
c_5 & c_{10} & c_{14} & c_{17} & c_{19} & c_{20}
\end{array} \right) \cdot \vec{p}^T \; .
\end{equation}
The quadric
$Q(\vec p)$ is an element in~$V := \bbbc [\vec{p}]_2 / \bbbc  \lbrace
p_{01}p_{23}-p_{02}p_{13} + p_{03}p_{12} \rbrace \simeq \bbbc^{21} / \bbbc$.
Hence, $\vec{c} =  (c_0, c_1, \ldots, c_{20})$ serves as homogeneous coordinates on $\bbbp^{19} = \bbbp(V)$, which~-- due to the Pl\"ucker relation~-- need to be understood modulo
\begin{equation}\label{eq:coordinate-indet}
c_5 \mapsto c_5 + \lambda, \,\,
c_9 \mapsto c_9 - \lambda, \,\,
c_{12} \mapsto c_{12} + \lambda \; .
\end{equation}
The coordinate ring $\bbbq[V]$ is a subring of~$\bbbq[c_0,c_1,\ldots,c_{20}]$,
namely it is the invariant ring of the additive group action~(\ref{eq:coordinate-indet}).
Hence $\bbbq[V]$ is the polynomial ring in  $20$ variables
$\,c_0,c_1,c_2,c_3,c_4,c_5-c_{12},c_6,c_7,c_8, c_9+c_{12},c_{10},c_{11},c_{13},\ldots,c_{20}  $.

We are interested in the $\vec{c}$'s that lead to \emph{coisotropic}
hypersurfaces of~$G(2,4)$. For these, the tangent space
at any point~$\ell$, considered as a subspace of~$\,T_\ell G(2,4) $ $ =
\mathrm{Hom}(\ell, \bbbc^4/\ell)$, has the form $\left\{ \varphi  \
\middle|\  \varphi(a) = 0  \right\} + \left\{\varphi \ \middle|\
\mathrm{im}(\varphi) \subset M  \right\}$, for some~$a\in\ell \backslash \{0\}$
 and some plane~$M$ in~$\bbbc^4/\ell$.
By~\cite[\S 4.3, (3.24)]{gkz}, the quadric hypersurface $\{Q(\vec p) = 0\}$
in~$G(2,4)$  is coisotropic if and only if there exist~$s,t \in \bbbc$ such that
\begin{equation}\label{eq:coisotropy_differential}
\frac{\partial Q}{\partial p_{01}} \cdot \frac{\partial Q}{\partial p_{23}} - \frac{\partial Q}{\partial p_{02}} \cdot \frac{\partial Q}{\partial p_{13}} + \frac{\partial Q}{\partial p_{03}} \cdot \frac{\partial Q}{\partial p_{12}}
= s \cdot Q + t \cdot \left( p_{01}p_{23}-p_{02}p_{13} + p_{03}p_{12} \right) \; .
\end{equation}
Equivalently, the vector $(t,s,-1)^T$ is in the  kernel of the $21 {\times}
3$ matrix in Figure~\ref{fig:bigmat}. The $3\times 3$ minors of this matrix are
all in the subring~$\bbbq[V]$.  The {\em coisotropic ideal}~$I$ is the ideal
of~$\bbbq[V]$ generated by these minors.  The  subscheme~$V(I)$
of~$\bbbp^{19}=\bbbp(V)$ represents all coisotropic hypersurfaces~$\lbrace Q =
0 \rbrace$ of degree two in~$G(2,4)$.
Using computations with {Maple} and {Macaulay2}~\cite{M2}, we found that~$I$
has codimension~10, degree~92 and is minimally generated by~175 cubics.
Besides,~$V(I)$ is the reduced union of three  components, of dimensions nine,
eight and five.
\looseness=-1

\begin{proposition} \label{prop:coisotropic}
The coisotropic ideal is  the intersection of three prime ideals:
\begin{equation}
\label{coisotropy_ideal_decomposition}
I \,\,= \,\,P_{\mathrm{Hurwitz}} \,\cap \, P_{\mathrm{ChowLines}} \,\cap \,P_{\mathrm{Squares}} \; .
\end{equation}
So, $I$ is radical.  The prime $P_{\mathrm{Hurwitz}}$ has codimension 10 and
degree~92, it is minimally generated by~20 quadrics, and its variety $V (
P_{\mathrm{Hurwitz}} )$ consists of Hurwitz forms of quadric surfaces in
$\bbbp^3$. The prime $P_{\mathrm{ChowLines}}$ has codimension~11 and degree~140,
it is minimally generated by~265 cubics, and $V (P_{\mathrm{ChowLines}} )$
consists of Chow forms of pairs of lines in $\bbbp^3$. The prime
$P_{\mathrm{Squares}}$ has codimension~14 and degree~32, it is minimally
generated by~84 quadrics, and $V ( P_{\mathrm{Squares}})$ consists of all
quadrics
 $Q(\vec{p})$ that are squares modulo the Pl\"ucker relation.
\end{proposition}

\begin{figure}
\centering
\[
 \left( \begin{array}{ccc}
    0 & {c}_{0} & 2  {c}_{0}  {c}_{5}-2  {c}_{1}  {c}_{4}+2 {c}_{2}  {c}_{3}\\
    0 & {c}_{1} & {c}_{0}  {c}_{10}-{c}_{1} {c}_{9}+{c}_{2}  {c}_{8}+{c}_{3}  {c}_{7}-{c}_{4} {c}_{6}+{c}_{1}  {c}_{5}\\
    0 & {c}_{2} & {c}_{0} {c}_{14}-{c}_{1}  {c}_{13}+{c}_{2}  {c}_{12}+{c}_{3} {c}_{11}-{c}_{4}  {c}_{7}+{c}_{2}  {c}_{5}\\
    0 & {c}_{3} & {c}_{0}  {c}_{17}-{c}_{1}  {c}_{16}+{c}_{2}  {c}_{15}+{c}_{3} {c}_{12}-{c}_{4}  {c}_{8}+{c}_{3}  {c}_{5}\\
    0 & {c}_{4} & {c}_{0}  {c}_{19}-{c}_{1}  {c}_{18}+{c}_{2}  {c}_{16}+{c}_{3} {c}_{13}-{c}_{4}  {c}_{9}+{c}_{4}  {c}_{5}\\
    1 & {c}_{5} & {c}_{0}  {c}_{20}-{c}_{1}  {c}_{19}+{c}_{2}  {c}_{17}+{c}_{3} {c}_{14}-{c}_{4}  {c}_{10}+{c}_{5}^{2}\\ 0 & {c}_{6} & 2 {c}_{1}  {c}_{10}-2  {c}_{6}  {c}_{9}+2  {c}_{7}  {c}_{8}\\
    0 & {c}_{7} & {c}_{1}  {c}_{14}-{c}_{6}  {c}_{13}+{c}_{7} {c}_{12}+{c}_{8}  {c}_{11}+{c}_{2}  {c}_{10}-{c}_{7} {c}_{9}\\
    0 & {c}_{8} & {c}_{1}  {c}_{17}-{c}_{6} {c}_{16}+{c}_{7}  {c}_{15}+{c}_{8}  {c}_{12}+{c}_{3} {c}_{10}-{c}_{8}  {c}_{9}\\
    -1 & {c}_{9} & {c}_{1} {c}_{19}-{c}_{6}  {c}_{18}+{c}_{7}  {c}_{16}+{c}_{8} {c}_{13}+{c}_{4}  {c}_{10}-{c}_{9}^{2}\\
    0 & {c}_{10} & {c}_{1}  {c}_{20}-{c}_{6}  {c}_{19}+{c}_{7}  {c}_{17}+{c}_{8} {c}_{14}-{c}_{9}  {c}_{10}+{c}_{5}  {c}_{10}\\
    0 & {c}_{11} & 2  {c}_{2}  {c}_{14}-2  {c}_{7}  {c}_{13}+2  {c}_{11} {c}_{12}\\
    1 & {c}_{12} & {c}_{2}  {c}_{17}-{c}_{7} {c}_{16}+{c}_{11}  {c}_{15}+{c}_{3}  {c}_{14}-{c}_{8} {c}_{13}+{c}_{12}^{2}\\
    0 & {c}_{13} & {c}_{2} {c}_{19}-{c}_{7}  {c}_{18}+{c}_{11}  {c}_{16}+{c}_{4} {c}_{14}+{c}_{12}  {c}_{13}-{c}_{9}  {c}_{13}\\
    0 & {c}_{14} & {c}_{2}  {c}_{20}-{c}_{7}  {c}_{19}+{c}_{11} {c}_{17}+{c}_{12}  {c}_{14}+{c}_{5}  {c}_{14}-{c}_{10} {c}_{13}\\
    0 & {c}_{15} & 2  {c}_{3}  {c}_{17}-2  {c}_{8} {c}_{16}+2  {c}_{12}  {c}_{15}\\
    0 & {c}_{16} & {c}_{3} {c}_{19}-{c}_{8}  {c}_{18}+{c}_{4}  {c}_{17}+{c}_{12} {c}_{16}-{c}_{9}  {c}_{16}+{c}_{13}  {c}_{15}\\
    0 & {c}_{17} & {c}_{3}  {c}_{20}-{c}_{8}  {c}_{19}+{c}_{12} {c}_{17}+{c}_{5}  {c}_{17}-{c}_{10}  {c}_{16}+{c}_{14} {c}_{15}\\
    0 & {c}_{18} & 2  {c}_{4}  {c}_{19}-2  {c}_{9} {c}_{18}+2  {c}_{13}  {c}_{16}\\
    0 & {c}_{19} & {c}_{4} {c}_{20}-{c}_{9}  {c}_{19}+{c}_{5}  {c}_{19}-{c}_{10} {c}_{18}+{c}_{13}  {c}_{17}+{c}_{14}  {c}_{16}\\
    0 & {c}_{20} & 2  {c}_{5}  {c}_{20}-2  {c}_{10}  {c}_{19}+2  {c}_{14} {c}_{17}
  \end{array} \right)
  \]
  \label{fig:bigmat}
  %
  \caption{This matrix has rank  $\leq 2$ if and only
  if the quadric given by~$\vec c$ is coisotropic.}
\end{figure}

This proposition answers a question due to Green and Morrison, who had asked in
 \cite{green_morrison} whether $I$ is radical.
 To derive the prime decomposition~(\ref{coisotropy_ideal_decomposition}), we computed the three prime ideals as kernels of homomorphisms of polynomial rings, each expressing the relevant
 geometric condition. This construction ensures that the ideals are prime. We then
verified that their intersection equals $I$. For details, check our computations,
 using the link given at the end of this article.

From the geometric perspective of \cite{gkz}, the third
prime $P_{\mathrm{Squares}}$ is extraneous,
because nonreduced hypersurfaces in $G(2,4)$ are excluded
by Gel'fand, Kapranov
and Zelevinsky.
 Theorem 3.14 in \cite[\S4.3]{gkz} concerns irreducible hypersurfaces, and the identification of Chow forms within the coisotropic hypersurfaces~\cite[\S4.3, Theorem 3.22]{gkz} assumes
   the corresponding polynomial to be squarefree.
With this, the following would be the correct ideal for the coisotropic variety in $\bbbp^{19}$:
\begin{equation}
\label{correct_coisotropic_ideal}
P_{\mathrm{Hurwitz}} \,\cap \,P_{\mathrm{ChowLines}}
\,\,= \,\, \left( I : P_{\mathrm{Squares}} \right) \; .
\end{equation}
This means that the reduced coisotropic quadrics in~$G(2,4)$ are either Chow forms of curves or Hurwitz forms of surfaces.
The ideal in~(\ref{correct_coisotropic_ideal}) has codimension~10, degree~92,
and is minimally generated by~175 cubics and~20 quartics in $\bbbq[V]$.

A slightly different point of view on the coisotropic ideal is presented in a recent
paper of Catanese \cite{catanese}.  He derives a variety in $\mathbb{P}^{20} =
\mathbb{P}(\mathbb{C}[{\bf p}]_2)$ which projects isomorphically onto our variety $V(I) \subset \mathbb{P}^{19}$.
The center of projection is the Pl\"ucker quadric.
To be precise,  Proposition 4.1 in \cite{catanese} states the following:
   For every $Q \in \bbbc[\vec{p}]_2 \backslash
   \bbbc  \left( p_{01}p_{23}-p_{02}p_{13} + p_{03}p_{12} \right)
$    satisfying~(\ref{eq:coisotropy_differential}) there is a unique $\lambda \in \mathbb{C}$ such that
the quadric $Q_\lambda := Q + \lambda \cdot \left( p_{01}p_{23}-p_{02}p_{13} + p_{03}p_{12} \right)$
satisfies
\begin{equation}\label{eq:coisotropy_catanese}
\frac{\partial Q_\lambda}{\partial p_{01}} \cdot \frac{\partial Q_\lambda}{\partial p_{23}} - \frac{\partial Q_\lambda}{\partial p_{02}} \cdot \frac{\partial Q_\lambda}{\partial p_{13}} + \frac{\partial Q_\lambda}{\partial p_{03}} \cdot \frac{\partial Q_\lambda}{\partial p_{12}}
=  t \cdot \left( p_{01}p_{23}-p_{02}p_{13} + p_{03}p_{12} \right)
\end{equation}
for some $t \in \mathbb{C}$.
This implies that $V(I)$ is isomorphic to the variety of all $Q \in \bbbp (\bbbc[\vec{p}]_2) \setminus \lbrace p_{01}p_{23}-p_{02}p_{13} + p_{03}p_{12} \rbrace$ satisfying~(\ref{eq:coisotropy_catanese}). Let $I_2$ be generated by the $2 \times 2$ minors of the $21 \times 2$ matrix that is obtained by deleting the middle column of the matrix in Figure~\ref{fig:bigmat}. Then $V(I_2)$ contains exactly those $Q \in \bbbp (\bbbc[\vec{p}]_2)$ satisfying~(\ref{eq:coisotropy_catanese}), and $V(I)$ is the projection of $V(I_2)$
from the center $(p_{01}p_{23}-p_{02}p_{13} + p_{03}p_{12})$. The ideal $I_2$ has codimension 11, degree 92, and is minimally generated by 20 quadrics. Interestingly, Catanese shows furthermore in \cite[Theorem 3.3]{catanese} that a hypersurface in $G(2,4)$ is coisotropic if and only if it is selfdual
in $\mathbb{P}^5$ with respect to the inner product given by the Pl\"ucker quadric.

\section{The Chow Variety}
\label{sec:chowvar}

In this section we study the Chow variety $G(2,2,4)$ of one-dimensional algebraic cycles of degree two in $\bbbp^3$. By~\cite[\S4.1, Ex.~1.3]{gkz}, the Chow variety $G(2,2,4)$ is the
union of two irreducible components of dimension eight in $\bbbp^{19}$, one
corresponding to planar quadrics and the other to pairs of lines. Formally, this means
that~$G(2,2,4) = V (P_{\mathrm{ChowConic}} ) \cup V (P_{\mathrm{ChowLines}})$,
where $P_{\mathrm{ChowConic}}$ is the homogeneous prime ideal in $\bbbq[V]$ whose variety comprises the Chow forms of irreducible curves of degree two in $\bbbp^3$. The ideal $P_{\mathrm{ChowConic}}$ has codimension 11 and degree 92, and it is minimally generated by 21 quadrics and 35 cubics.
The radical ideal $P_{\mathrm{ChowConic}} \cap P_{\mathrm{ChowLines}}$ has codimension 11, degree $232=92+140$, and  it is minimally generated by 230 cubics.

Since $G(2,2,4)$ should be contained in the coisotropic variety
$V(I)$, it seems that $P_{\mathrm{ChowConic}}$ is missing from the decomposition (\ref{coisotropy_ideal_decomposition}).  Here is the explanation:

\begin{proposition} \label{prop:chowhurw}
Every Chow form of a plane conic in $\bbbp^3$ is also a Hurwitz form.
In symbols,
$\,P_{\mathrm{Hurwitz}} \,\subset \,P_{\mathrm{ChowConic}} \,$ and thus
$\,V ( P_{\mathrm{ChowConic}} ) \,\subset V ( P_{\mathrm{Hurwitz}}) $.
\end{proposition}

Our first proof is by computer: just check the inclusion of ideals in Macaulay2.
For a conceptual proof, we consider a $4 \times 4$-symmetric matrix
 $M = M_0  + \epsilon M_1$, where ${\rm rank}(M_0) =1$.
By \cite[eqn.~(1)]{Stu}, the Hurwitz form
 of the corresponding quadric surface in $\bbbp^3\,$ is
 $\,Q(\vec{p}) = \vec{p} (\wedge_2 M) \vec{p}^T$.
 Divide by $\epsilon$ and  let $\epsilon \rightarrow 0$.
 The limit is  the Chow form of the plane
  conic defined by restricting $M_1$ to ${\rm ker}(M_0) \simeq \bbbp^2$.
  This type of degeneration is familiar from the study of complete quadrics
  \cite{CGMP}.    Proposition \ref{prop:chowhurw}
  explains why the locus of irreducible curves
  is  not visible in~(\ref{coisotropy_ideal_decomposition}).

Gel'fand, Kapranov and Zelevinsky \cite[\S 4.3]{gkz} introduce
a class of differential forms in order to  discriminate Chow forms
among all coisotropic hypersurfaces.
In their setup, these forms represent the integrability of the $\alpha$-distribution
$\mathcal{E}_{\alpha,Z}$.  We shall apply the tools of computational
commutative algebra to shed some light on the
characterization of Chow forms via integrability of $\alpha$-distributions.

For this, we use local affine coordinates instead of Pl\"ucker coordinates.
A point in the Grassmannian $G(2,4)$ is represented as the row space of
the matrix
\begin{equation}
\label{eq:twobyfour}
\left( \begin{array}{cccc}
1 &\, 0 &\, a_2 & a_3 \\
0 &\, 1 & \, b_2 & b_3
\end{array} \right) \; .
\end{equation}
We express the quadrics $Q$ in (\ref{general_quadric}) in terms of the local coordinates
$a_2,a_3,b_2,b_3$, by substituting the Pl\"ucker coordinates with the minors of the
matrix (\ref{eq:twobyfour}), i.e.,
\begin{equation}
p_{01} = 1, \,\,
p_{02} = b_2, \,\,
p_{03} = b_3, \,\,
p_{12} = -a_2, \,\,
p_{13} = -a_3, \,\,
p_{23} = a_2 b_3 - b_2 a_3 \; .
\end{equation}
We consider the following differential $1$-forms on affine $4$-space:
\[
\alpha_1^1 := \frac{\partial Q}{\partial a_2} d a_2 + \frac{\partial Q}{\partial a_3} d a_3,\quad
 \alpha_2^1 := \frac{\partial Q}{\partial a_2} d b_2 + \frac{\partial Q}{\partial a_3} d b_3,
\]
\[
\alpha_1^2 := \frac{\partial Q}{\partial b_2} d a_2 + \frac{\partial Q}{\partial b_3} d a_3,\quad
 \alpha_2^2 := \frac{\partial Q}{\partial b_2} d b_2 + \frac{\partial Q}{\partial b_3} d b_3 \; .
\]
By taking wedge products, we derive the $16$ differential $4$-forms
\begin{equation}
d Q \wedge d \alpha^i_j \wedge \alpha^k_l
\,\,=\,\, q_{ijkl} \cdot d a_2 \wedge d a_3 \wedge d b_2 \wedge d b_3
\,\,\quad \hbox{for $i,j,k,l \in \lbrace 1,2 \rbrace$} \; .
\end{equation}
Here the expressions $q_{ijkl} $
are certain polynomials in $ \bbbq[V][a_2,a_3,b_2,b_3]$.

Theorems~3.19 and~3.22 in~\cite[\S4.3]{gkz} state that a squarefree coisotropic
quadric $Q$ is a Chow form if and only if all~16 coefficients~$q_{ijkl}$ are
multiples of~$Q$.  By taking normal forms of the polynomials $q_{ijkl}$ modulo
the principal ideal $\langle Q \rangle$, we obtain a collection of $720$
homogeneous polynomials in $\vec{c}$.  Among these,~$58$ have degree three,~$340$
have degree four, and~$322$ have degree five.  The aforementioned result
implies that these~$720$ polynomials cut out~$G(2,2,4)$ as a subset of~$\bbbp^{19}$.

The {\em integrability ideal} $J \subset  \bbbq[V]$
is generated by these $720$ polynomials and their analogues from
other affine charts of the Grassmannian, obtained by permuting
columns in~(\ref{eq:twobyfour}).  We know that $V(J)$ equals the union of $G(2,2,4)$ with all double hyperplanes in $G(2,4)$ (corresponding to $P_{\mathrm{Squares}}$)
set-theoretically.   Maple, Macaulay2 and Magma verified for us that it holds
scheme-theoretically:

     \begin{proposition}
The integrability ideal~$J$ is minimally generated by~210 cubics.
Writing $\mathfrak m$ for the irrelevant ideal~$\langle c_0,c_1,\ldots,c_{20} \rangle$
of~$\bbbq[V]$, we have
\begin{equation}
\sqrt{J}\,\, =\,\, (J : \mathfrak m)\, \,=\,\, P_{\mathrm{ChowConic}} \,\cap\, P_{\mathrm{ChowLines}} \,\cap\, P_{\mathrm{Squares}} \; .
\end{equation}
   \end{proposition}

\section*{Conclusion}

We reported on computational experiments with hypersurfaces in the Grassmannian
$G(2,4)$ that are associated to curves and surfaces in $\bbbp^3$.  For degree
$2$, all relevant parameter spaces were described by explicit polynomials in
$20$ variables.
All ideals and computations discussed in this note can be obtained at
\begin{center}
  \url{www3.math.tu-berlin.de/algebra/static/pluecker/}
\end{center}

Many possibilities exist for future work.  Obvious next
milestones are the ideals for the Chow varieties of degree $3$ cycles in
$\bbbp^3$,  and degree $2$ cycles in~$\bbbp^4$.  Methods from  representation
theory promise  a compact  encoding of their generators, in terms of
irreducible ${\rm GL}(4)$-modules.  Another question we aim to pursue is
motivated by the geometry  of condition numbers \cite{condition}: express the
volume of a tubular neighborhood of a coisotropic quadric in   ${\rm G}(2,4)$
as a function   of~$\vec{c}$.

\bibliographystyle{splncs}