[]Gabriel Nowaskie

[]Department of Physics and Astronomy, The Gatton Academy, Western Kentucky
University, 1906 College Heights Blvd., Bowling Green, 42101, KY, United
States

# Operator Space Manifold Theory: Modeling Quantum Operators with a Riemannian
Manifold

<EMAIL_ADDRESS>*

###### Abstract

The Half-Transform Ansatz (HTA) is a proposed method to solve hyper-geometric
equations in Quantum Phase Space by transforming a differential operator to an
algebraic variable and including a specific exponential factor in the wave
function, but the mechanism which provides this solution scheme is not known.
Analysis of the HTA’s application to the Hydrogen atom suggests an underlying
mechanism which the HTA is a part of. Observations of exponential factors that
act on the wave function naturally suggest modeling quantum operator
definitions as a point on a Riemannian manifold in the 4D Operator Space, a
novel idea we call the Operator Space Manifold Theory. Expanding on this
concept, we find the true nature of the HTA and how Operator Space Manifold
Theory can be used to describe and solve quantum systems by manipulating how a
quantum state perceives position and momentum.

###### keywords:

Half Transform Ansatz, HTA, Quantum Phase Space, Operator Space Manifold,
Hydrogen Phase Space, Operator Space, Schrodinger Equation

##  1 Introduction

With the discovery of the QPSR by Torres-Vega and Frederick [1], a framework
with the capabilities to describe quantum systems in phase space, and the
Heaviside Operational Ansatz (HOA) was created and published in 2004 by Simpao
to generate exact analytical wave functions in the QPSR [2, 3]. The HOA yields
implicit solutions in Quantum Phase Space. By applying a Laplace transform on
the $p$ differential operator, we developed a general scheme that seemed to
work for finding explicit solutions to the Schrodinger Equation or any similar
hyper-geometric equation with the use of the Nikiforov-Uvarov (NU) Method [4].
This was called the Half-Transform Ansatz (HTA), but there is a lack of
theoretical development justifying this method or explaining its nature.

This work aims to produce a theoretical backing to the HTA as continuation to
historical developments in [4], where we develop our Operator Space Manifold
Theory. We develop this theory in the following sections and this paper is
organized as follows. Section 2 is a recap of the Quantum Phase Space
formulation and the HOA. Section 3 is a brief summary of the NU method.
Section 4 provides the original HTA as in [4]. Section 5 shows the application
of the HTA to the Hydrogen atom and its discrepancies with the configuration
space results. Section 6 presents the Operator Space Manifold Theory. Section
7 presents how to traverse the Operator Space Manifold by the use of
transforming one of the differential operators to an algebraic variable.
Section 7 concludes our results and presents a full description of the HTA
with the insight gained from the Operator Space Manifold Theory.

##  2 Recap of the QPSR and Hamiltonian Construction

To begin, we must revisit the Heaviside Operational Ansatz. Only the brief
methods relevant to this work are presented. Historical developments can be
seen in [2, 3]. The variable relations: $x,p,t$ are the configuration space
position, momentum, and time respectively. When transforming to the QPSR we
recall that (where $\hat{}$ denotes operators and $\alpha,\gamma$ are
otherwise free parameters as in [1]):

$\displaystyle H\left(x,\ p,\ t\right)\rightarrow\ \hat{H}\left(\hat{x},\
\hat{p},\ t\right)=\hat{H}\left(i\hbar\partial_{p}+\alpha x,\
-i\hbar\partial_{x}+\gamma p,\ t\right),\ \ni\alpha+\gamma=1$ $\displaystyle
x\rightarrow\ \hat{x}\equiv\ i\hbar\partial_{p}+\alpha\ x,\
p\rightarrow\hat{p}\equiv-i\hbar\partial_{x}+\gamma\ p,t\rightarrow\ t=t$
$\displaystyle\left(x_{1},\ldots,x_{n}\right)\rightarrow\left(\widehat{x_{1}},\ldots,\widehat{x_{n}}\right)=\left(i\hbar\partial_{p_{1}}+\alpha_{1}x_{1},\ldots,i\hbar\partial_{p_{n}}+\alpha_{n}x_{n}\right),\ni\alpha_{j}+\gamma_{j}=1,\
j=1,\ldots,n$ $\displaystyle
H\left(x_{1},\ldots,x_{n};p_{1},\ldots,p_{n};t\right)\rightarrow\hat{H}\left({\hat{x}}_{1},\ldots,{\hat{x}}_{n};{\hat{p}}_{1},\ldots,{\hat{p}}_{n};t\right)$
$\displaystyle\equiv\hat{H}\left(i\hbar\partial_{p_{1}}+\alpha_{1}x_{1},\ldots,i\hbar\partial_{p_{n}}+\alpha_{n}x_{n};-i\hbar\partial_{x_{1}}+\gamma_{1}p_{1},\ldots,-i\hbar\partial_{x_{n}}+\gamma_{n}p_{n};t\right).$
(1)

For our nucleus-electron system, our frame of reference is the nucleus. We
thus transform to a spherical coordinate system and separate the angular
variables, where $L(L+1)$ is the angular coupling term, ${p_{r}}$ is the
radial momentum, and $r$ is the radial distance from the electron to the
nucleus. For a 1D system, the Hamiltonian is:

$\displaystyle\hat{H}=\frac{{\hat{{p_{r}}}}^{2}}{2m}+\hat{V}(r).$ (2)

For a 3D system, the Hamiltonian is:

$\displaystyle\hat{H}=\sum^{N}{\frac{-i\hbar{\hat{p}}_{r}}{m\hat{r}}+\frac{{\hat{p}}_{r}^{2}}{2m}+\frac{\hbar^{2}L(L+1)}{2m{\hat{r}}^{2}}}+\hat{V}(r).$
(3)

Typically, $\hat{r}=-i\hbar\frac{\partial}{\partial{p_{r}}}+\alpha r,\
\hat{p}_{r}=-i\hbar\frac{\partial}{\partial r}+\gamma{p_{r}}$ by convention of
the QPSR. The more general form presented in [3] is:

$\displaystyle\hat{r}=\alpha r+i\hbar\beta\frac{\partial}{\partial{p_{r}}},\
\hat{{p_{r}}}=\gamma{p_{r}}+i\hbar\delta\frac{\partial}{\partial r}\ ,$ (4)

where $\hat{r},\hat{p}_{r}$ are bound by the condition
$\beta\gamma-\alpha\delta=1$ in order to keep consistency of $\hat{r}$ and
$\hat{p}_{r}$ with the Heisenberg uncertainty principle,
$[\hat{r},\hat{p}_{r}]=1$.

##  3 Nikiforov-Uvarov Method

The Nikiforov-Uvarov method is based on solving the hyper-geometric type
second-order differential equation of form

$\displaystyle{\
\Psi}^{\prime\prime}\left(s\right)+\frac{\widetilde{\tau}}{\sigma\left(s\right)}\Psi^{\prime}\left(s\right)+\frac{\
\widetilde{\sigma}\left(s\right)}{\sigma^{2}\left(s\right)}\Psi\left(s\right)=0$
(5)

Where $\sigma(s)$ and $\widetilde{\sigma}\left(s\right)$ are at most, second
degree polynomials, and $\widetilde{\tau}\left(s\right)$ is at most a first-
degree polynomial. $\Psi\left(s\right)$ is a function of hyper-geometric type.
The solution of (5) takes the form

$\displaystyle\Psi\left(s\right)=\phi\left(s\right)y\left(s\right)$ (6)

Substituting (6) into (5) yields

$\displaystyle\sigma\left(s\right)y^{\prime\prime}\left(s\right)+\tau\left(s\right)y^{\prime}\left(s\right)+\lambda\left(s\right)y\left(s\right)=0,$
(7)

where $\phi\left(s\right)$ satisfies the following relation

$\displaystyle\frac{\phi^{\prime}\left(s\right)}{\phi(s)}=\frac{\pi\left(s\right)}{\sigma(s)}.$
(8)

and $y\left(s\right)$ is a hyper-geometric type function, whose polynomial
solutions are obtained from Rodrigues’ relation

$\displaystyle\
y\left(s\right)=y_{n}\left(s\right)=\frac{B_{n}}{\rho\left(s\right)}\frac{d^{n}}{ds^{n}}\left[\sigma^{n}\left(s\right)\rho\left(s\right)\right].$
(9)

where $B_{n}$ is the normalization constant and $\rho\left(s\right)$ is a
weight function satisfying the equation

$\displaystyle\left[\sigma\left(s\right)\rho\left(s\right)\right]^{\prime}=\tau\left(s\right)\rho\left(s\right).$
(10)

The function $\pi\left(s\right)$ is defined as

$\displaystyle\pi\left(s\right)=\left(\frac{\sigma^{\prime}-\widetilde{\tau}}{2}\right)\pm\sqrt{\left(\frac{\sigma^{\prime}-\widetilde{\tau}}{2}\right)^{2}-\widetilde{\sigma}+K\sigma},$
(11)

and $\lambda$ is defined as

$\displaystyle\lambda=K+\pi\prime.$ (12)

The value of $K$ can be calculated under the condition that the square root in
(11) must be the square of a polynomial. Thus, the equation of eigenvalues can
be given as:

$\displaystyle\lambda=\lambda_{n}=-n\tau^{\prime}-\frac{n\left(n-1\right)}{2}\sigma^{\prime\prime},$
(13)

where

$\displaystyle\tau\left(s\right)=\widetilde{\tau}\left(s\right)+2\pi\left(s\right).$
(14)

##  4 The Half-Transform Ansatz

Here we present the HTA. Its original statement can be seen in [4]. Using the
general QPSR definitions from (1), the time independent wave equation becomes

$\displaystyle\hat{H}\binom{i\hbar\partial_{p_{r_{1}}}+\alpha_{1}x_{1},\ldots,i\hbar\partial_{p_{r_{n}}}+\alpha_{n}x_{n};}{-i\hbar\partial_{x_{1}}+\gamma_{1}p_{r_{1}},\ldots,-i\hbar\partial_{x_{n}}+\gamma_{n}p_{r_{n}}}\psi\left(x_{1},\ldots,x_{n};p_{r_{1}},\ldots,p_{r_{n}}\right)=E_{n}\psi(x_{1},\ldots,x_{n};p_{r_{1}},\ldots,p_{r_{n}}).$
(15)

Applying the alternative to the convolution and multi-variable inverse
transform (3A), the convolution of ${p_{r}}$ between (4A) and relation (5A)
yields

$\begin{aligned}
\left[L_{\binom{\left({\partial_{p_{r_{1}}},\ldots,\partial_{p_{r_{n}}}}\right)}{\rightarrow\left(p_{r_{1}},\ldots.,p_{r_{n}}\right)}}^{-1}\left[\hat{H}\binom{i\hbar\partial_{p_{r_{1}}}+\alpha_{1}x_{1},\ldots,i\hbar\partial_{p_{r_{n}}}+\alpha_{n}x_{n};}{-i\hbar\partial_{x_{1}}+\gamma_{1}p_{r_{1}},\ldots,-i\hbar\partial_{x_{n}}+\gamma_{n}p_{r_{n}}}\right]\right]\\\
\ast\ \ \ \ \
\psi\left(x_{1},\ldots,x_{n};p_{r_{1}},\ldots,p_{r_{n}}\right)\\\
\end{aligned}\leavevmode\hbox to23.16pt{\vbox
to5.4pt{\pgfpicture\makeatletter\hbox{\hskip 5.89046pt\lower-5.2pt\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ } {}{{}}{} {}{}
{}{{}{}}{}{}{}{{}}{{}}{{}{}}{{}{}}{{{{}{}{{}} }}{{}} {} {}{}{} { {{}} {}
{}{}{} {}{}{} } { {{}} {} {}{}{} } }{{}{}}{{}{}}{{{{}{}{{}} }}{{}}}
{}\pgfsys@moveto{17.07182pt}{0.0pt}\pgfsys@moveto{17.07182pt}{0.0pt}\pgfsys@curveto{16.32185pt}{-1.50002pt}{14.57182pt}{-2.5pt}{12.07182pt}{-2.5pt}\pgfsys@lineto{10.69069pt}{-2.5pt}\pgfsys@curveto{8.19069pt}{-2.5pt}{6.44066pt}{-3.49998pt}{5.69069pt}{-5.0pt}\pgfsys@curveto{4.94072pt}{-3.49998pt}{3.19069pt}{-2.5pt}{0.69069pt}{-2.5pt}\pgfsys@lineto{-0.69046pt}{-2.5pt}\pgfsys@curveto{-3.19046pt}{-2.5pt}{-4.94049pt}{-1.50002pt}{-5.69046pt}{0.0pt}\pgfsys@stroke\pgfsys@invoke{
} \pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}\leavevmode\hbox
to48.11pt{\vbox to16.67pt{\pgfpicture\makeatletter\hbox{\hskip
24.05667pt\lower-8.33301pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }
{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}}
{{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{-20.72366pt}{-2.5pt}\pgfsys@invoke{
}\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\hbox{{$\left(p_{r_{1}},\ldots,p_{r_{n}}\right)$}}
}}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}$
$\displaystyle\equiv\hat{H}\binom{i\hbar\partial_{p_{r_{1}}}+\alpha_{1}x_{1},\ldots,i\hbar\partial_{p_{r_{n}}}+\alpha_{n}x_{n};}{-i\hbar\partial_{x_{1}}+\gamma_{1}p_{r_{1}},\ldots,-i\hbar\partial_{x_{n}}+\gamma_{n}p_{r_{n}}}\psi\left(x_{1},\ldots,x_{n};p_{r_{1}},\ldots,p_{r_{n}}\right).$
(16)

Applying (16)-(15) with the convolution identities in (3A), (16) is
transformed into

$\begin{aligned}
L_{\binom{\left({p_{r_{1}},\ldots,p_{r_{n}}}\right)}{\rightarrow\left({\overline{p}_{r}}_{1},\ldots.,{\overline{p}_{r}}_{n}\right)}}\
&\left[\hat{H}\binom{i\hbar\partial_{p_{r_{1}}}+\alpha_{1}x_{1},\ldots,i\hbar\partial_{p_{r_{n}}}+\alpha_{n}x_{n};}{-i\hbar\partial_{x_{1}}+\gamma_{1}p_{r_{1}},\ldots,-i\hbar\partial_{x_{n}}+\gamma_{n}p_{r_{n}}}\psi\left(x_{1},\ldots,x_{n};p_{r_{1}},\ldots,p_{r_{n}}\right)\right]\\\
&\equiv\
L_{\binom{\left({p_{r_{1}},\ldots,{p_{r}}}_{n}\right)}{\rightarrow\left({\overline{p}_{r}}_{1},\ldots.,{\overline{p}_{r}}_{n}\right)}}\begin{bmatrix}L_{\binom{\left({\partial_{p_{r_{1}}},\ldots,\partial}_{p_{r_{n}}}\right)}{\rightarrow\left(p_{r_{1}},\ldots.,p_{r_{n}}\right)}}^{-1}\left[\hat{H}\binom{i\hbar\partial_{p_{r_{1}}}+\alpha_{1}x_{1},\ldots,i\hbar\partial_{p_{r_{n}}}+\alpha_{n}x_{n};}{-i\hbar\partial_{x_{1}}+\gamma_{1}p_{r_{1}},\ldots,-i\hbar\partial_{x_{n}}+\gamma_{n}p_{r_{n}}}\right]\\\
\ast\ \ \ \ \ \ \ \
\psi\left(x_{1},\ldots,x_{n};p_{r_{1}},\ldots,p_{r_{n}}\right)\end{bmatrix}\\\
&=E_{n}\check{\psi}\left(x_{1},\ldots,x_{n};{\overline{p}_{r}}_{1},\ldots,{\overline{p}_{r}}_{n}\right)\end{aligned}\leavevmode\hbox
to23.16pt{\vbox to5.4pt{\pgfpicture\makeatletter\hbox{\hskip
5.89046pt\lower-5.2pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ } \par {}{{}}{} {}{}
{}{{}{}}{}{}{}{{}}{{}}{{}{}}{{}{}}{{{{}{}{{}} }}{{}} {} {}{}{} { {{}} {}
{}{}{} {}{}{} } { {{}} {} {}{}{} } }{{}{}}{{}{}}{{{{}{}{{}} }}{{}}}
{}\pgfsys@moveto{17.07182pt}{0.0pt}\pgfsys@moveto{17.07182pt}{0.0pt}\pgfsys@curveto{16.32185pt}{-1.50002pt}{14.57182pt}{-2.5pt}{12.07182pt}{-2.5pt}\pgfsys@lineto{10.69069pt}{-2.5pt}\pgfsys@curveto{8.19069pt}{-2.5pt}{6.44066pt}{-3.49998pt}{5.69069pt}{-5.0pt}\pgfsys@curveto{4.94072pt}{-3.49998pt}{3.19069pt}{-2.5pt}{0.69069pt}{-2.5pt}\pgfsys@lineto{-0.69046pt}{-2.5pt}\pgfsys@curveto{-3.19046pt}{-2.5pt}{-4.94049pt}{-1.50002pt}{-5.69046pt}{0.0pt}\pgfsys@stroke\pgfsys@invoke{
} \pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}\leavevmode\hbox
to48.11pt{\vbox to16.67pt{\pgfpicture\makeatletter\hbox{\hskip
24.05667pt\lower-8.33301pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }
{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}}
{{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{-20.72366pt}{-2.5pt}\pgfsys@invoke{
}\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\hbox{{$\left(p_{r_{1}},\ldots,p_{r_{n}}\right)$}}
}}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}$ (17)
$\displaystyle\equiv\hat{H}\binom{i\hbar{\overline{p}_{r}}_{1}+\alpha_{1}x_{1},\ldots,i\hbar{\overline{p}_{r}}_{n}+\alpha_{n}x_{n};}{-i\hbar\partial_{x_{1}}+\gamma_{1}p_{r_{1}},\ldots,-i\hbar\partial_{x_{n}}+\gamma_{n}p_{r_{n}}}\check{\psi}\left(x_{1},\ldots,x_{n};{\overline{p}_{r}}_{1},\ldots,{\overline{p}_{r}}_{n}\right)={E}_{n}\check{\psi}\left(x_{1},\ldots,x_{n};{\overline{p}_{r}}_{1},\ldots,{\overline{p}_{r}}_{n}\right).$
(18)

Thus, the wave function in phase space can be expressed as the inverse
transform of the solution
$\breve{\psi}\left(x_{1},\ldots,x_{n};{\overline{p}_{r}}_{1},\ldots,{\overline{p}_{r}}_{n};t\right)$.
The radial Schrodinger equation using the 3D full Cornell potential
Hamiltonian is

$\displaystyle\left(-\frac{i\hbar\hat{{p_{r}}}}{m\hat{r}}+\frac{{\hat{{p_{r}}}}^{2}}{2m}+\frac{\hbar^{2}L\left(L+1\right)}{2m\hat{r}}+V(\hat{r})\right)\psi\left(r,\overline{p}_{r}\right)=E_{n}\psi\left(r,\overline{p}_{r}\right).$
(19)

We expand the ansatz previously stated in [4] by stating that the formulation:

$\displaystyle\psi\left(r,{p_{r}}\right)=e^{\frac{i{p_{r}}r}{\hbar}\left(\frac{\gamma}{\delta}\right)}\
\Omega\left(r,{p_{r}}\right),$ (20)

where $\delta$ and $\gamma$ are defined as the coefficients in the momentum
operator, degenerates the Schrodinger Equation in Quantum Phase Space to a
form similar to the configuration space. By making the substitution $A=\alpha
r+i\hbar\beta\overline{p}_{r}$, the Schrodinger equation can be written as an
ordinary differential equation (ODE). This ODE is hyper-geometric, and thus
can be solved using the Nikiforov Uvarov method. To then get the solution in
the domain of $\left(r,p_{r}\right)$, an inverse Laplace transform from
$(\overline{p}_{r}\to p_{r})$ is done. The effects of Heisenberg’s uncertainty
principle and the nature of the Quantum Phase Space can be seen using the HTA.
To present this, we present the application of the HOA to the Hydrogen Atom.

##  5 Application of the HTA to Hydrogen

In this section, we apply the HTA to the 3D TISE in Quantum Phase Space.

Using the general QPSR definitions from (1) along with the Coulomb potential,
the radial, time independent wave equation becomes

$\displaystyle\left(-\frac{i\hbar\hat{{p_{r}}}}{m\hat{r}}+\frac{{\hat{{p_{r}}}}^{2}}{2m}+\frac{\hbar^{2}L(L+1)}{2m\hat{r}}-\frac{ke^{2}}{\hat{r}}\right)\psi\left(r,{\overline{p}_{r}}\right)=E_{n}\psi\left(r,\overline{p}_{r}\right).$
(21)

Using the ansatz
$\psi\left(r,{p_{r}}\right)=e^{\frac{i{p_{r}}r}{\hbar}\left(\frac{\gamma}{\delta}\right)}\Omega\left(r,{p_{r}}\right)$
along with the general operator definitions as in (4) and mapping
$(\partial_{{p_{r}}}\rightarrow\overline{p}_{r})$:

$\displaystyle\begin{pmatrix}-\frac{2e^{2}k\Omega\left(r,\overline{p}_{r}\right)e^{\frac{i\gamma
rp_{r}}{\delta\hbar}}}{\alpha
r+i\beta\overline{p}_{r}\hbar}+\frac{L^{2}\hbar^{2}\Omega\left(r,\overline{p}_{r}\right)e^{\frac{i\gamma
rp_{r}}{\delta\hbar}}}{m(\alpha r+i\beta\overline{p}_{r}\hbar)^{2}}\\\
+\frac{L\hbar^{2}\Omega\left(r,\overline{p}_{r}\right)e^{\frac{i\gamma
rp_{r}}{\delta\hbar}}}{m(\alpha
r+i\beta\overline{p}_{r}\hbar)^{2}}+\frac{2\delta\hbar^{2}e^{\frac{i\gamma
rp_{r}}{\delta\hbar}}}{m(\alpha
r+i\beta\overline{p}_{r}\hbar)}\frac{\partial\Omega\left(r,\overline{p}_{r}\right)}{\partial
r}-\frac{\delta^{2}\hbar^{2}e^{\frac{i\gamma
rp_{r}}{\delta\hbar}}}{m}\frac{\partial^{2}\Omega\left(r,\overline{p}_{r}\right)}{\partial
r^{2}}\end{pmatrix}=2E_{n}\Omega\left(r,\overline{p}_{r}\right)e^{\frac{i\gamma
rp_{r}}{\delta\hbar}}$ (22)

which becomes

$\displaystyle\begin{pmatrix}-\frac{2e^{2}k\Omega\left(r,\overline{p}_{r}\right)}{\alpha
r+i\beta\overline{p}_{r}\hbar}+\frac{L^{2}\hbar^{2}\Omega\left(r,\overline{p}_{r}\right)}{m(\alpha
r+i\beta\overline{p}_{r}\hbar)^{2}}\\\
+\frac{L\hbar^{2}\Omega\left(r,\overline{p}_{r}\right)}{m(\alpha
r+i\beta\overline{p}_{r}\hbar)^{2}}+\frac{2\delta\hbar^{2}}{m(\alpha
r+i\beta\overline{p}_{r}\hbar)}\frac{\partial\Omega\left(r,\overline{p}_{r}\right)}{\partial
r}-\frac{\delta^{2}\hbar^{2}}{m}\frac{\partial^{2}\Omega\left(r,\overline{p}_{r}\right)}{\partial
r^{2}}\end{pmatrix}=2E_{n}\Omega\left(r,\overline{p}_{r}\right)$ (23)

We proceed with the substitution $A=\alpha r+i\hbar\beta\overline{p}_{r}$. The
chain rule yields

$\displaystyle\ \frac{\partial\psi\left(r,\overline{p}_{r}\right)}{\partial
r}=\frac{\partial\Omega(A)}{\partial A}\frac{\partial A}{\partial
r}=\alpha\frac{\partial\Omega(A)}{\partial r},\ \
\frac{\partial^{2}\Omega\left(r,\overline{p}_{r}\right)}{\partial
r^{2}}=\frac{\partial^{2}\Omega(A)}{\partial A^{2}}\frac{\partial
A^{2}}{\partial r^{2}}=\alpha^{2}\frac{\partial^{2}\Omega(A)}{\partial
r^{2}}.$ (24)

(23) transforms into

$\displaystyle\Omega(A)\left(\frac{L^{2}\hbar^{2}}{A^{2}m}+\frac{L\hbar^{2}}{A^{2}m}-\frac{2e^{2}k}{A}\right)-2E_{n}\Omega(A)+\frac{2\alpha\delta\hbar^{2}\Omega^{\prime}(A)}{Am}-\frac{\alpha^{2}\delta^{2}\hbar^{2}\Omega^{\prime\prime}(A)}{m}=0$
(25)

$\displaystyle=\frac{\left(\frac{2A^{2}E_{n}m}{\hbar^{2}}+\frac{2Ae^{2}km}{\hbar^{2}}-L(L+1)\right)}{A^{2}\alpha^{2}\delta^{2}}\Omega(A)-\frac{2\Omega^{\prime}(A)}{A\delta\alpha}+\Omega^{\prime\prime}(A)=0$
(26)

For convenience, we define three constants, $\omega,\zeta,\kappa$, to be

$\displaystyle\omega\equiv L(L+1),\ \ \zeta\equiv\
\frac{2e^{2}km}{\hbar^{2}},\ \ \kappa\equiv-\frac{2E_{n}m}{\hbar^{2}}.$ (27)

Now, we can rewrite (26) as

$\displaystyle\frac{\left(-\omega+A\zeta-A^{2}\kappa\right)}{\alpha^{2}A^{2}\delta^{2}}\Omega(A)-\frac{2\Omega^{\prime}(A)}{\alpha
A\delta}+\Omega^{\prime\prime}(A)=0,$ (28)

the exact form as in (5) to apply the Nikiforov-Uvarov method. Thus, we let

$\displaystyle\sigma\left(A\right)=-\alpha\delta
A,\quad\widetilde{\sigma}\left(A\right)=-\omega+A\zeta-A^{2}\kappa,\quad\widetilde{\tau}\left(A\right)=2.$
(29)

From (11)

$\displaystyle\pi\left(A\right)=\frac{1}{2}(-\alpha\delta-2)+\pm\sqrt{\frac{1}{4}(-\alpha\delta-2)^{2}+A^{2}\kappa-A\zeta-\alpha
A\delta K+\omega}$ (30)

Since the polynomial under the radical must a square of a polynomial, the
discriminant must equal zero.

$\displaystyle(-4\zeta-4\alpha\delta
K)^{2}-4(4\kappa)\left(\alpha^{2}\delta^{2}+4\alpha\delta+4\omega+4\right)=0$
(31)

Solving for $K$,

$\displaystyle
K=-\frac{\alpha\delta\zeta+|\alpha\delta|\sqrt{\kappa\left((\alpha\delta+2)^{2}+4\omega\right)}}{\alpha^{2}\delta^{2}},\quad\frac{-\alpha\delta\zeta+|\alpha\delta|\sqrt{\kappa\left((\alpha\delta+2)^{2}+4\omega\right)}}{\alpha^{2}\delta^{2}}.$
(32)

In an attempt to get a result similar to the configuration space Schrodinger
equation, we want the expression inside the square inside the radical to be a
perfect square, creating the condition: $\alpha\delta+2=1$. Thus,

$\displaystyle\alpha=\frac{-1}{\delta},\quad\alpha=\frac{-3}{\delta}.$ (33)

These options are still bound by Heisenberg’s uncertainty principle,
$\beta\gamma-\alpha\delta=1$. We realize that the option
$\alpha=\frac{-1}{\delta}$ will result in the configuration space operator
definitions, and so our Schrodinger equation solution would become the
configuration space solution. Thus, due to the constraint of Heisenberg’s
uncertainty principle and keeping this equation in phase space, we chose the
option $\alpha=\frac{-3}{\delta}$, which makes the condition:

$\displaystyle\beta\gamma=-2.$ (34)

Due to the exponential factor substitution earlier, it can be seen that the
ansatz,
$\psi(r,p_{r})=e^{\left(\frac{ip_{r}r\gamma}{\hbar\delta}\right)}\Omega(r,p_{r})$,
truly is a degeneration of the PDE to the configuration space. However, this
degeneration is restricted by Heisenberg’s uncertainty principle, and the
factor of 3 in (33) is what embeds the intrinsic uncertainty between the
position and momentum in our phase space equation. We chose the value of $K$
that results in a negative derivative of $\tau$. Thus, we chose the negative
sign. Simplifying K, we have

$\displaystyle K=\frac{1}{3}\left(\zeta-\sqrt{4\kappa\omega+\kappa}\right).$
(35)

This results in

$\begin{aligned}
\pi\left(A\right)=\frac{1}{2}\pm\frac{1}{2}\left(2\sqrt{\kappa}A-\sqrt{4\omega+1}\right)\end{aligned}.$
(36)

The choice of the sign in (36) must once again generate a negative first
derivative in $\tau$; thus, we take the negative sign again:

$\begin{aligned}
\pi\left(A\right)=\frac{1}{2}-\frac{1}{2}\left(2\sqrt{\kappa}A-\sqrt{4\omega+1}\right)\end{aligned}.$
(37)

From (8), $\phi(A)$ can be found by the differential equation

$\displaystyle\frac{\phi^{\prime}\left(A\right)}{\phi(A)}=\frac{\pi\left(A\right)}{\sigma(A)}=-\frac{-A\sqrt{\kappa}+\frac{1}{2}\sqrt{4\omega+1}+\frac{1}{2}}{\alpha
A\delta}=\frac{-A\sqrt{\kappa}+\frac{1}{2}\sqrt{4\omega+1}+\frac{1}{2}}{3A},$
(38)

which yields the exponential solution

$\displaystyle\phi\left(A\right)=e^{\frac{A\sqrt{\kappa}}{\alpha\delta}}A^{-\frac{\sqrt{4\omega+1}+1}{2\alpha\delta}}=e^{-\frac{1}{3}\left(A\sqrt{\kappa}\right)}A^{\frac{1}{6}\left(\sqrt{4\omega+1}+1\right)}.$
(39)

By (14), $\tau$ is

$\displaystyle\tau\left(A\right)=-2A\sqrt{\kappa}+\sqrt{4\omega+1}+3.$ (40)

(10) is used to find the weight function $\rho\left(A\right)$, where

$\displaystyle\left[\sigma\left(A\right)\rho\left(A\right)\right]^{\prime}=\tau\left(A\right)\rho\left(A\right),$
(41)
$\displaystyle\rho(A)\left(\alpha\delta-2A\sqrt{\kappa}+\sqrt{4\omega+1}+3\right)+\alpha
A\delta\rho^{\prime}(A)=0,$ (42)

which gives the solution

$\displaystyle\rho\left(A\right)=e^{\frac{2A\sqrt{\kappa}}{\alpha\delta}}A^{-\frac{\alpha\delta+\sqrt{4\omega+1}+3}{\alpha\delta}}=e^{-\frac{1}{3}\left(2A\sqrt{\kappa}\right)}A^{\frac{1}{3}\sqrt{4\omega+1}}.$
(43)

With $\rho(A)$ and $\sigma\left(A\right)$, the function $y$ can be formulated
as

$\displaystyle
y_{n}\left(A\right)=\frac{B}{\rho\left(A\right)}\frac{d^{n}}{dA^{n}}\left(\sigma^{n}\left(A\right)\rho\left(A\right)\right).$
(44)

where $B$ is a normalization constant. Recalling (6), the wave function must
be

$\displaystyle\psi_{n}\left(A\right)=e^{\frac{i{p_{r}}r}{\hbar}\left(\frac{\gamma}{\delta}\right)}\Omega(A)=e^{\frac{i{p_{r}}r}{\hbar}\left(\frac{\gamma}{\delta}\right)}y_{n}\left(A\right)\phi\left(A\right)=e^{\frac{i{p_{r}}r}{\hbar}\left(\frac{\gamma}{\delta}\right)}\phi\left(A\right)\frac{B}{\rho\left(A\right)}\
\frac{d^{n}}{dA^{n}}\left(\sigma^{n}\left(A\right)\rho\left(A\right)\right)$
$\displaystyle=Be^{\frac{i{p_{r}}r}{\hbar}\left(\frac{\gamma}{\delta}\right)}e^{\frac{A\sqrt{\kappa}}{3}}A^{\frac{1}{6}-\frac{1}{6}\sqrt{4\omega+1}}\frac{d^{n}}{dA^{n}}\left(e^{-\frac{1}{3}\left(2A\sqrt{\kappa}\right)}A^{\frac{1}{3}\sqrt{4\omega+1}}(-\alpha
A\delta)^{n}\right).$ (45)

Thus, to find the wave function with respect to $(r,p_{r})$, (45) can be
inversely Laplace transformed from $\overline{p_{r}}\to p_{r}$ after back-
substituting $A$. Now, to cross check our results with the coordinate space
hydrogen atom, we look at the energy eigenvalues using equations (12) and
(13):

$\displaystyle\lambda=K+\pi^{\prime}\left(A\right)=\frac{1}{3}\left(\zeta-\sqrt{4\kappa\omega+\kappa}\right)-\sqrt{\kappa},$
(46) $\displaystyle\lambda_{n}=\
-n\tau^{\prime}\left(A\right)-\frac{n\left(n-1\right)\sigma^{\prime\prime}\left(A\right)}{2}=2n\sqrt{\kappa}.$
(47)

With condition $\lambda=\lambda_{n}$,

$\displaystyle\frac{1}{3}\left(\zeta-\sqrt{4\kappa\omega+\kappa}\right)-\sqrt{\kappa}=2\sqrt{\kappa}n$
(48)

Solving this equation for $\kappa$, we have

$\displaystyle\kappa=\frac{\zeta^{2}}{\left(6n+\sqrt{4\omega+1}+3\right)^{2}}$
(49)

Back-substituting $\omega\equiv L(L+1),\ \zeta\equiv\
\frac{2e^{2}km}{\hbar^{2}},\ \kappa\equiv-\frac{2\text{En}m}{\hbar^{2}}$, we
find the energy eigenvalue equation to be

$\displaystyle E_{n}=-\frac{e^{4}k^{2}m}{2\hbar^{2}(L+3n+2)^{2}}.$ (50)

This result is not the same as the configuration space solution, and is a
direct consequence of the relationship $\alpha=\frac{-3}{\delta}.$ If we were
to take (48) and (45) and have the variables $\alpha$ and $\delta$ not
simplified, we would see that for $\alpha\to\frac{-1}{\delta}$ and $\gamma=0$
(which conforms with the momentum operator definition in the configuration
space), these expressions would converge to the exact Hydrogen atom wave
function and energy eigenvalues. Because we chose to have a wave function
including both $p$ and $r$, the Heisenberg uncertainty principle embedded a
factor of $3$ that represents the error in measuring $r$ and $p$
simultaneously. This error uniquely spread across the energy eigenvalues and
the wave function depending on what path you take in solving (28). For the
values of $\alpha,\beta,\gamma,\delta$ that define the coordinate space or the
momentum space, the inverse Laplace transform isn’t needed since the
coefficient in front of $p$ or $r$ will then be $0$. Though, this requires a
mechanism that exists which dictates the nature of the operator coefficients
$\alpha,\beta,\gamma,\delta$ the Schrodinger Equation, or any quantum
equation, before this inverse transform. In other words, there must then exist
a framework that defines the nature of the quantum operator coefficients that
is independent of the $\overline{p}$ space previously used. We also need this
framework to describe the nature of the
$e^{\frac{i{p_{r}}r}{\hbar}\frac{\gamma}{\delta}}$ factor.

##  6 Operator Space Manifold Theory

In this section we present a mechanism to describe all of the concerns in the
HTA that generates a overarching theory to describe the Quantum Phase Space.

Consider a 4D space which has 4 degrees of freedom:
$\alpha,\beta,\gamma,\delta.$ In this space there is a curve that is
represented by $\beta\gamma-\alpha\delta=1$, such that Heisenberg’s
uncertainty principle lives on this curve and every point on this curve
represents a possible definition of phase space. Looking at (23),
$e^{\frac{i{p_{r}}r}{\hbar}\frac{\gamma}{\delta}}$ seems to have the the
effect of forcing $\gamma$ to converge to $0$. We interpret this exponent as a
rotation of the form $e^{i\phi}$ and assume that this exponential factor is an
operator that translates a point on the curve $\beta\gamma-\alpha\delta=1$ to
have $\gamma=0$; therefore, each point on our curve $\beta\gamma-\alpha\delta$
must have defined angles. This naturally leads us to use a Riemannian manifold
as the machinery for our framework.

More formally stated, let us consider a Riemannian manifold defined by
$\beta\gamma-\alpha\delta=1$ in the $\alpha,\beta,\gamma,\delta$ 4D Operator
Space. A ”specialized” phase space is the space with operator definitions
defined by any point on the manifold; we can think of these specialized phase
spaces as the image of a point on the manifold. We consider the ”entire” phase
space to be the complete set of specialized phase spaces, or the set of all
images from all points on the manifold. By definition, this manifold also
includes the coordinate space and the momentum space operator definitions,
which are seen as subsets of the entire phase space.

We define each point in the Operator Space by:

$\begin{pmatrix}\alpha\\\ \beta\\\ \gamma\\\ \delta\end{pmatrix}.$

$e^{\frac{i{p_{r}}r}{\hbar}\frac{\gamma}{\delta}}$ is seen as an operator that
converges $\gamma\to 0$, and thus is defined to be a rotational transform on
this manifold such that angle
$\phi=\frac{{p_{r}}r}{\hbar}\frac{\gamma}{\delta}$ is the angle between the
tangent space of a particular point and the 3D space where $\gamma$ is one
less than initially ($1$ in this scenario). To mathematically represent this
we postulate that

$\displaystyle e^{i\phi}\psi=g_{\eta}\begin{pmatrix}\alpha\\\ \beta\\\
\gamma\\\ \delta\end{pmatrix}$ (51)

where $g_{\eta}$ is a $4\times 4$ matrix that acts as the linear transform in
Operator Space analog to a phase shift by $\phi$. These rotational transforms
can also be seen as a translation to a specific point
$(\alpha,\beta,\gamma,\delta)$. For simplicity, we only consider this
translational definition continuing forward. (51) is also a mathematical way
to represent the idea that a quantum system is analog to a point on the
Operator Space Manifold. The effects of a specific $\phi$ on the Operator
Space can be seen when applying $e^{i\phi}$ onto $\psi$. We present a table of
the most basic transformations:

Table 1: 4 fundamental $g_{\eta}$ matrices for different $\phi$. $\phi$ | $g_{\eta}$
---|---
$\displaystyle\phi_{1}=\frac{p_{r}r}{\hbar}\left(\frac{\alpha}{\beta}\right)$ | $\displaystyle\begin{pmatrix}0&0&0&0&\\\ 0&1&0&0&\\\ 0&0&1&0&\\\ 0&0&0&1&\\\ \end{pmatrix}={g_{\eta}}_{1}$
$\displaystyle\phi_{3}=\frac{p_{r}r}{\hbar}\left(\frac{\gamma}{\delta}\right)$ | $\displaystyle\begin{pmatrix}1&0&0&0&\\\ 0&1&0&0&\\\ 0&0&0&0&\\\ 0&0&0&1&\\\ \end{pmatrix}={g_{\eta}}_{3}$
$\displaystyle\phi_{2}=i\log\left(\psi\left(r,p_{r}\right)\right)+C\times r$ | $\displaystyle\begin{pmatrix}1&0&0&0&\\\ 0&0&0&0&\\\ 0&0&1&0&\\\ 0&0&0&1&\\\ \end{pmatrix}={g_{\eta}}_{2}$
$\displaystyle\phi_{4}=i\log\left(\psi\left(r,p_{r}\right)\right)+C\times p_{r}$ | $\displaystyle\begin{pmatrix}1&0&0&0&\\\ 0&1&0&0&\\\ 0&0&1&0&\\\ 0&0&0&0&\\\ \end{pmatrix}={g_{\eta}}_{4}$
|

where $C$ is an arbitrary complex constant. The indices $n$ for
${g_{\eta}}_{n}$ are simply decided by which $1$ across the diagonal becomes
$0$. For notation simplicity when developing the algebra of these matrices, we
also introduce the idea of the a complement matrix, ${g_{\eta}}_{C}$:

$\displaystyle{g_{\eta}}_{C}=I-{g_{\eta}}$ (52)

where $I$ is the identity matrix. The complements of the $g_{\eta}$ matrices
(${g_{\eta}}_{C}$) represent the manipulation on a specific operator
coefficient, or the causality, while the $g_{\eta}$ matrices represent the
result of that manipulation on the entire space, or the effect.

When applying $e^{\frac{i{p_{r}}r}{\hbar}\frac{\gamma}{\delta}}$ twice to
$\psi(r,p_{r})$ we have

$\displaystyle i\delta\hbar\psi^{(1,0)}\left(r,p_{r}\right)e^{\frac{2i\gamma
rp_{r}}{\delta\hbar}}-\gamma p_{r}\psi\left(r,p_{r}\right)e^{\frac{2i\gamma
rp_{r}}{\delta\hbar}},$ (53)

a result corresponding to

$\displaystyle g_{\eta}=\begin{pmatrix}1&0&0&0&\\\ 0&1&0&0&\\\ 0&0&-1&0&\\\
0&0&0&1&\\\ \end{pmatrix}.$ (54)

Applying this transform $n$ times, we see a pattern: each transform shifts the
index on the diagonal affecting $\gamma$ by $-1$ each time, starting
originally from the identity matrix $I$. In other words, the original
$g_{\eta}$ matrix is shifted by ${g_{\eta}}_{C}$ n times, where n is the
number of times you apply the specific operator. This same result holds true
for the other transforms and their respective operator coefficients they
manipulate, and we present a mathematical representation of this below:

$\displaystyle{g^{\prime}_{\eta}}=-n{g_{\eta}}_{C}+{g_{0}},$ (55)

where $g^{\prime}_{\eta}$ is the resultant $g_{\eta}$ vector, ${g_{0}}$ is the
initial $g_{\eta}$ matrix you are working with (typically the Identity matrix
$I$ to represent the entire Operator Space), and ${g_{\eta}}_{C}$ is the
complement of the $g_{\eta}$ matrix for the specific operator you are applying
to $\psi$. For example, if we were to begin with the entire Operator Space
($g_{0}=I$), and applied $e^{\frac{i{p_{r}}r}{\hbar}\frac{\gamma}{\delta}}$ 2
times, we would have the resultant $g^{\prime}_{\eta}$ matrix:

$\displaystyle g^{\prime}_{\eta}=-2\begin{pmatrix}0&0&0&0&\\\ 0&0&0&0&\\\
0&0&1&0&\\\ 0&0&0&0&\end{pmatrix}+\begin{pmatrix}1&0&0&0&\\\ 0&1&0&0&\\\
0&0&1&0&\\\ 0&0&0&1&\end{pmatrix}=\begin{pmatrix}1&0&0&0&\\\ 0&1&0&0&\\\
0&0&-1&0&\\\ 0&0&0&1&\end{pmatrix},$ (56)

the same result as shown in (53). Thus,
${g_{\eta}}_{1},{g_{\eta}}_{2},{g_{\eta}}_{3},{g_{\eta}}_{4}$ are unit
transforms that shift $\alpha,\beta,\gamma,$ or $\delta$ by -1. We define the
inverse of these $e^{i\phi}$ transforms to be when you make the substitution
$\phi\to-\phi$, which then results in a shift of $+1$ instead of $-1$, or
${g_{\eta}}_{C}=-{g_{\eta}}_{C}$.

After defining these transforms that allow us to traverse the Operator Space
Manifold, it is tempting to assume that they all must commute since the basis
created from ${g_{\eta}}_{1},{g_{\eta}}_{2},{g_{\eta}}_{3},{g_{\eta}}_{4}$
covers the entire Operator Space. However, these matrices have a restriction
that appears due to the interaction of $\phi$ with differential operators
creating an additional cross-term. Thus, we can only apply operators
manipulating $\alpha$ and $\beta$ (${g_{\eta}}_{1}$ and ${g_{\eta}}_{2}$) or
$\gamma$ and $\delta$ (${g_{\eta}}_{3}$ and ${g_{\eta}}_{4}$) simultaneously.
The operators that can be applied simultaneously also commute. Therefore, we
can only traverse the Operator Space by the position or the momentum operator
coefficients at a time. This seems to be an expression of Heisenberg’s
uncertainty principle in the Operator Space.

This result, Heisenberg’s uncertainty principle, cannot be worked around, as
assuming there is a translational operator that is able to traverse the entire
Operator Space at once, a cross term appears after the application of either
the momentum or position operator which cannot be canceled due to the
condition resolving to $e^{ix}=0$.

While we have stated that we cannot apply operators that will affect
$\alpha,\beta,\gamma$, and $\delta$ simultaneously, what about these operators
for equations or formulations that include both the momentum and position
operators? Up until now, every formulation presented has been from applying
single operators to $\psi(r,p_{r})$ or choosing which operators can be
applied. But the momentum and position operators will exist in almost every
phase space formulation. For these scenarios, any of these operators on the
Operator Space Manifold would inherently affect the entire Operator Space
since its being applied to both the position and momentum operators. We will
address this problem in the next section, which will also explain the
importance of the HTA.

When applying multiple operators that are allowed to be used at the same time
(those affecting only either the position or momentum operator coefficients at
the same time), it is seen that the same principle from (55) follows, except
${g_{\eta}}_{C}$ can is the addition of any of the ${g_{\eta}}_{C}$ from the
operators being applied:

$\displaystyle{g^{\prime}_{\eta}}=-{g_{\eta}}_{C1}-{g_{\eta}}_{C2}-....-{g_{\eta}}_{Cn}+{g_{0}},$
(57)

where ${g_{\eta}}_{Cn}$ is simply the ${g_{\eta}}_{C}$ for the nth operator
that is being applied simultaneously.

We can summarize the algebra of these matrices as follows: only operators that
act on $\alpha$ and $\beta$ or $\gamma$ and $\delta$ can be used
simultaneously. The operators that can be used simultaneously commute, and a
resultant ${g_{\eta}}_{C}$ matrix is found by simply summing the matrices.
This same result holds true for applying the same operators multiple times.
Thus, the Operator Space Manifold can only be simultaneously traversed using
operators on either $\alpha$ and $\beta$ or $\gamma$ and $\delta$. The
complement matrix represents the causality of operators, and can be seen as
the shift from an original set of operator coefficient definitions to a new
one. This can be extended to any quantum system by stating that a system can
only traverse space exactly by either $r$ or $p$. This result, known as
Heisenberg’s uncertainty principle, can thus find origins in the Operator
Space.

We also briefly introduce the idea of ”position-like” or ”configuration-like”
phase spaces and ”momentum-like” phase spaces. We define ”position-like” phase
spaces to be specialized phase spaces with a $g_{\eta}$ of

$\displaystyle g_{\eta}=\begin{pmatrix}1&0&0&0&\\\ 0&0&0&0&\\\ 0&0&0&0&\\\
0&0&0&1&\end{pmatrix},$ (58)

and ”momentum-like” phase spaces are specialized phase spaces with a
$g_{\eta}$ of

$\displaystyle g_{\eta}=\begin{pmatrix}0&0&0&0&\\\ 0&1&0&0&\\\ 0&0&1&0&\\\
0&0&0&0&\end{pmatrix}.$ (59)

These two 2D spaces represent the analysis of a quantum system by its position
or momentum. It is worth noting that the position and momentum space are
subsets of these spaces respectively, and represent the most simple solutions
to $\beta\gamma-\alpha\delta=1.$ This is possibly a representation of the
Universe taking ”least action” that is reflected within the Operator Space
Manifold.

##  7 Traversing the Operator Space for Formulations with both Position and
Momentum Operators

In this section we will address the problem of traversing the Operator Space
considering the restrictions due to Heisenberg’s uncertainty principle suggest
that it is not allowed.

By the principles of the Operator Space from the last section, we know that
the Operator Space Manifold can only be traversed by operators affecting
$\alpha,\beta$ or $\gamma,\delta$ simultaneously because the the cross-
interaction of the differential operators and the operator coefficients. We
also expressed how for equations or formulations that include both the
momentum and position operators, any of the operators presented would
inherently affect the entire Operator Space since its being applied to both
the position and momentum operators. So, how are we supposed to traverse the
Operator Space Manifold at all?

In order to bypass this restriction, we can take the Laplace transform one of
the differential operators to an algebraic variable. This will redefine one of
our quantum operators such that it will act like a variable. This move is the
vital stepping stone of the HTA, and allows us to traverse the Operator Space
Manifold by either $\alpha,\beta$ or $\gamma,\delta$ simultaneously. The
restriction on only being able to traverse the Operator Space by either
$\alpha,\beta$ or $\gamma,\delta$ still exists because of the cross-
interaction between the angle $\phi$ in $e^{i\phi}$ and the derivative of the
non-transformed quantum operator. Then, after solving the equation or doing
the calculations, the inverse Laplace transform can be used to recover the
information about $p$ or $r$.

But, if we apply a transform that makes either $\alpha,\gamma,\beta,\gamma$
become $0$ for the equation, then the inverse Laplace transform will simply
return the original input. This is because the half transform used in the HTA
can be changed to include the operator coefficient when taking the transform
using the same logic in Section 4. If the coefficient operator is transformed
to become 0, then the transform will be the input multiplied by the Dirac
Delta function of $0$ which will return back the original input is not a
function of either $\alpha r,\gamma p,\delta\overline{r},$ or
$\beta\overline{p}$ anymore. This fully explains how to recover the Hydrogen
atom configuration space results from the phase space Schrodinger Equation,
and can be used to recover any momentum or configuration space result from its
phase space formulation.

The last thing to cover is the residual exponentials which are left behind
after applying one of these Operator Space transformations such as in (45)
during the final formulation of the wave function. When applying these
Operator Space transformations on a state, it is changing the way that state
perceives its quantum operators. When such an ansatz of $\psi=e^{i\phi}\Omega$
is done in an equation, it is modeling that equation by assuming that there is
a function $\Omega$ that is undergoing a transform from one point on the
Operator Space Manifold to another. This changes how the equation $\psi$ is in
will perceive position and momentum.

Because this transform is happening within $\psi$ and not on $\psi$, it does
not change the way $\psi$ itself perceives position and momentum. Thus, while
the equation of $\Omega$ perceives some transform in of
($\alpha,\beta,\gamma,\delta$), $\psi$ should not, and should still contain
information about the entire Operator Space ($\alpha,\beta,\gamma,\delta$).

This conservation of information can be seen in the residual exponential
factor when doing the inverse Laplace transform at the end of evaluating
$\psi$, in which the residual exponential should factor out of the transform
and still exist in the final wave function result. In summary, by using the
HTA with the Operator Space Manifold, we can model a state, $\psi$, in
equations by an Operator Space transform acting on a function, $\omega$, such
that that $\omega$ has specific quantum operator definitions. These definition
are isolated to only affect $\omega$ without actually losing that information
in $\psi$. This information on the transformed operator coefficients is held
within in these residual exponentials.

##  8 Conclusion: Operator Space Manifold Theory and the Nature of the HTA

By applying the HTA as in [4] to the hydrogen atom and considering the
differences between the phase space and the configuration space results, it is
seen that there must be an underlying mechanism dictating the nature of
$\alpha,\beta,\gamma$, and $\delta$ such that these values can be transformed
to new values. Through observation of specific factor acting on $\psi$ and
interpreting them as translations on some space, we develop the Operator Space
Manifold Theory, which models each possible definition of the quantum
operators as a point on a Riemannian manifold with specific restriction and
behaviors as discussed in Sections 6 and 7.

From section 7, we know that the Operator Space Manifold cannot be traversed
for equations with $\hat{r}$ and $\hat{p}$ unless a differential operator is
transformed into an algebraic variable and that corresponding operator is
represented as a collective variable. Thus, we present the HTA as:
transforming a differential operator of a phase space formulation and then
treating the corresponding quantum operator whose differential operator was
transformed as a collective variable. An example for this method can be seen
when. As posited in section 7, the HTA by definition is the method that allows
us to traverse the Operator Space Manifold for equations that have $\hat{r}$
and $\hat{p}$.

Equations with $\hat{r}$ and $\hat{p}$ constitutes almost every formulation in
phase space, and thus the HTA is the primary method to utilize the Operator
Space Manifold to solve equations or generate Phase Space formulations. For
equations or formulations that do not include both quantum operators, the HTA
is not needed to traverse the Operator Space Manifold. However, both scenarios
are still restricted by only traversing the Operator Space with
transformations on only $\alpha,\beta$ or $\gamma,\delta$ simultaneously.

Traversing the Operator Space Manifold can be used to uniquely express
equations and solve equations. One of these methods is the ansatz
$\psi=e^{i\phi}\Omega$, where $e^{i\phi}$ is an Operator Space transformation.
The simplest scenario is using one of the four fundamental transformations as
in table 1. These unit transforms translate a point on the Operator Space
manifold. When stating $\psi=e^{i\phi}\Omega$, the Operator Space transform
$e^{i\phi}$ only acts on $\Omega$. In short, we can use this ansatz to model
$\psi$ by a function which resides in a specific part of the Operator Space
Manifold (has restricted or specific quantum operator definitions). Thus,
$\Omega$ is perceiving its state by a different metric than $\psi$, but the
information changed in $\Omega$ is only localized and is not reflected in
$\psi$ because the residual factor $e^{i\phi}$ left after the transformation
holds the information about the change in the Operator Space basis.

From the developments presented in the paper, suggest not only the validity of
the HTA, but the powerful utilization of Operator Space Manifold Theory, which
the HTA is apart of. It provides a unique toolbox to solve equations and
represent formulations in Quantum Phase Space not seen elsewhere.

##  9 Future Developments and Acknowledgements

### Future Developments

The main path of this work for the future is to further develop the Operator
Space Manifold Theory to provide solutions to the Quantum Phase Space Dirac
Equation. We also plan to use the developments seen here to develop
configuration space formulations in the Quantum Phase Space.

### Acknowledgements

I would like to express my gratitude to the Gatton Academy and the Department
of Physics and Astronomy at Western Kentucky University for their support in
this research.

## References

  * [1] Torres-Vega, G.; Frederick, J.H., _A quantum-mechanical representation in phase space._ , _J. Chem. Phys._ 98 (1993) pg.3103-3120 https://doi.org/10.1063/1.464085
  * [2] Simpao, V. A., _Real wavefunction from Generalized Hamiltonian Schrodinger Equation in quantum phase space via HOA (Heaviside Operational Ansatz): Exact analytical results._ , _J. of Math. Chem._ 52 (2014) pg.1136-1155 https://doi.org/10.1007/s10910-014-0332-2
  * [3] Simpao, V. A., _Toward chemical applications of Heaviside operational Ansatz: exact solution of radial Schrodinger equation for nonrelativistic N-particle system with pairwise 1/rij radial potential in quantum phase space._ , _J. Math. Chem._ 45 (2009) pg.129-140 https://doi.org/10.1007/s10910-008-9372-9
  * [4] Nowaskie, G., _The Half Transform Ansatz: Quarkonium Dynamics in Quantum Phase Space._ (Version 2.0), _ArXiv_ (2023) https://arxiv.org/abs/2303.16356

## Statements and Declarations

The authors declare that no funds, grants, or other support were received
during the preparation of this manuscript.

Financial Interests: The authors have no relevant financial or non-financial
interests to disclose.

The independent author, Gabriel Nowaskie, did all work concerning this paper
from the writing of the manuscript to the calculations and data analysis. All
authors read and approved of this final manuscript.


# Exact polarization energy for clusters of contacting dielectrics

Huada Lian  Department of Materials Science & Engineering, Stanford University
Jian Qin  Department of Chemical Engineering, Stanford University
Corresponding author<EMAIL_ADDRESS>
###### Abstract

The induced surface charges appear to diverge when dielectric particles form
close contacts. Resolving this singularity numerically is prohibitively
expensive because high spatial resolution is needed. We show that the strength
of this singularity is logarithmic in both inter-particle separation and
dielectric permittivity. A regularization scheme is proposed to isolate this
singularity, and to calculate the exact cohesive energy for clusters of
contacting dielectric particles. The results indicate that polarization energy
stabilizes clusters of open configurations when permittivity is high, in
agreement with the behavior of conducting particles, but stabilizes the
compact configurations when permittivity is low.

For particle aggregates or clusters stabilized by electrostatic
interactions,1, 2, 3 the cohesive energy depends on the dielectric
permittivities of both the particles and the medium. Permittivity quantifies
the density of dipoles induced by externally applied electric fields, which is
proportional to polarization in the linear regime. When the permittivity
contrast between the particles and the medium is high, as is often the case,
the polarizations from the two sides of the interface do not fully compensate
each other, resulting in the accumulation of induced surface charges.

Resolving the surface charges is needed to evaluating the electrostatic
interactions among particles in close proximity, and is challenging because
polarization is intrinsically a many-body effect, depending on the positions
of all particles. For instance, careful measurements in colloidal suspensions
have shown that the inter-particle force is non-additive, which is at least
partially due to the polarization effect.4 In another set of experiments on
metallic nanoparticles, it has been found that the aggregation of multiple
particles surrounding a charged particle can be stabilized by the polarization
effect alone.3, 5 More dramatic demonstration of such polarization effects is
found in the so-called like-charge attraction caused by the strong
polarization effect, for particles of large size or permittivity ratios.6, 7

Computational methods, such as the boundary element method,8 the spectral
methods 9, 10 and the image method,11 have been developed to account for this
polarization effect. In practice, these methods all need to evaluate the
induced surface charge densities in one form or another, and have been
successfully applied to study a wide range of problems involving the
aggregation of polarizable particles.12, 13, 14 However, when particles are in
close proximity, the surface charge densities apparently diverge because the
electric field in the narrow gap region is strong even for a small difference
in the electrostatic potentials of particles. This is analogous to the
divergent lubrication force between approaching solid particles in the
Stokesian regime.15 Resolving these apparently diverging charge densities
requires a high spatial resolution that becomes computationally prohibitive
for nearly touching particles.8 For conductors like metalic particles, the
functional form of the divergent charge density has been identified and used
to isolate the singularity obscuring numerical calculations, and to obtain the
exact energy and force for contacting particles.5 For dielectric particles,
the question remains unsolved.

Figure 1:  Surface charges between two particles separated by a distance $h$.
The separation of two surfaces at a radial distance $r$ from the contact
points can be approximated by $h+r^{2}/a$, where $a\equiv
2/(a_{1}^{-1}+a_{2}^{-1})$. The permittivity of the medium is $\epsilon_{\rm
out}$. Surface of conducting particles are equipotential with potential set to
$V$ and $0$. Surfaces of dielectric particles, with permittivities
$\epsilon_{1}$ and $\epsilon_{2}$, have surface potential distributions
$\phi_{1}(r)$ and $\phi_{2}(r)$, respectively.

One way to reveal these singular surface charges is to consider two conducting
particles separated by a small gap, as sketched in Fig. 1. The points on the
two surfaces separated by the minimum distance, or gap distance $h$, are
referred to as contact points. The surface charge density is proportional to
the strength of the normal component of electric field on surfaces, according
to Gauss’s law. When $h$ is small, the electric field lines in the gap region
are nearly parallel to the line connecting the two contact points. The
strength of the electrical field is, in turn, proportional to the difference
in the surface electrical potentials. This relation between the surface
potential and the electrical field is analogous to that between the
longitudinal velocities of converging particles and the transverse velocity of
squzzed fluid velocity in the lubrication theory.15

To calculate the surface charge density, both the surface potentials and the
vertical separations are needed. Since conductors are equipotential, the
surface potentials can be set to $V$ and $0$ respectively. The vertical
separation depends on the radial distance $r$, and can be expressed within the
Derjaguin approximation as $h+r^{2}/a$, where $a$ is the harmonic average of
radii of curvature $a_{1}$ and $a_{2}$ at the two contact points, i.e.,
$a^{-1}=(a_{1}^{-1}+a_{2}^{-1})/2$. Here, the spherical apexes are assumed for
simplicity, but more general curvatures can be treated similarly.
Consequently, the field strength at radius $r$ is $E(r)=V/(h+r^{2}/a)$, and
the surface charge density $\sigma(r)=\frac{1}{4\pi}\epsilon_{\rm
out}V/(h+r^{2}/a)$ on the top surface, where $\epsilon_{\rm out}$ is the
medium permittivity. Integrating $\sigma(r)$ for $r$ from $0$ to $R_{0}$,
where $R_{0}$ is a cutoff of order $a$, results in the singular part of the
surface charge

$\displaystyle Q_{\rm s}$ $\displaystyle=\frac{V\epsilon_{\rm
out}}{4\pi}\int_{0}^{R_{0}}\\!{\rm d}r\frac{2\pi
r}{h+r^{2}/a}\simeq\frac{V\epsilon_{\rm out}a}{4}\ln\left(\frac{a}{h}\right).$
(1)

The unity in the logarithmic term is dropped because $R_{0}\simeq a\gg h$.
Further, the difference between $R_{0}$ and $a$ is neglected because it only
leads to a constant shift. Equation (1) shows how surface charges of the upper
surface become singular as the gap distance decreases. That of the lower
surface is also singular, but of opposite sign. In the limit $h\to 0$, the
energy does not blow up because the potential difference $V$ vanishes once
particles form contact.

The ratio of the singular charge $Q_{\rm s}$ and the potential difference $V$
gives the singular part of the capacitance $c_{\rm s}=\epsilon_{\rm
out}a\ln(a/h)/4$, which has been found previous for spherical dimers.6, 7 The
above analysis shows that this singular capacitance is local. Thus, the same
singularity applies to non-spherical particles, and for aggregates of multiple
particles. This fact has been employed to find the exact energy for clusters
of contacting, conducting particles.5 In this work, the full capacitance array
is first numerically calculated for an ensemble of conducting particles at
finite but small separations. The singular contribution is then subtracted,
leaving a regular part that can be extrapolated to $h=0$. Finally, when the
singular and regular parts are pieced together, the variation of energy with
separation is obtained.

For the dielectric case, a straightforward generalization of the above
treatment fails, because the particles are not equipotential. The potential
difference $\Delta V(r)$ needed to evaluate the charge density in eq. (1) can
not be fully specified by the potential difference between the contact points
$\Delta V(0)$. The variation of $\Delta V(r)$ with $r$ in the gap region is
expected to be quadratic, but the curvature is unknown a priori. In an early,
and analogous work on thermal conduction of composite materials, Batchelor et
al.16 noticed that the potential distribution, which determines the surface
charges, is itself dominated by the contribution from the surface charges
nearby, so that a self-consistent treatment is needed.

To illustrate this point, we consider the dielectric case in Fig. 1. Let the
surface potentials of the two particles be $\phi_{1}(r)$ and $\phi_{2}(r)$,
and the surface charges be $\sigma_{1}(r)$ and $\sigma_{2}(r)$. The surface
potential $\phi_{i}(r)$, with $i=1,2$, can be calculated by integrating the
coulombic potential of the respective surface charges. Near the contact point,
we have

$\phi_{i}(r)=\tilde{V_{i}}+\frac{1}{2\pi\epsilon_{i}}\int_{0}^{\infty}\\!\\!{\rm
d}r^{\prime}\\!\\!\int_{0}^{2\pi}\\!\\!{\rm
d}\theta^{\prime}\frac{r^{\prime}\sigma_{i}(r^{\prime})}{\sqrt{r^{2}+r^{\prime
2}-2rr^{\prime}\cos\theta^{\prime}}}.$ (2)

Here, $\tilde{V}_{i}$ are the contributions to the surface potential from the
surface charge outside the contact region. The integral are those from the
surface charges in the contact region. The prefactor is $\frac{1}{2\pi}$
(instead of $\frac{1}{4\pi}$) because of the well-known jump condition for
surface potentials.17 The distance at denominator is approximated using that
for the flat surface, which leads to negligible error because only the contact
region is of concern. The upper bound is set to infinity for convenience; as
we shall see below, the singular surface charge density dies off rapidly
outside the contact region.

The surface charge density $\sigma_{i}(r)$ in eq. (2) are related to the
potential difference, $\Delta V(r)\equiv\phi_{1}(r)-\phi_{2}(r)$. By analogy
to the conductor case, the electric field is approximately vertical and its
magnitude is $\Delta V(r)/(h+r^{2}/a)$. Then applying Gauss’s law, we get the
charge density on the top surface

$\sigma_{1}(r)=\epsilon_{\rm out}\left(1-\epsilon_{{\rm
r},1}^{-1}\right)\frac{\Delta V(r)}{h+r^{2}/a},$ (3)

where $\epsilon_{{\rm r},1}\equiv\epsilon_{1}/\epsilon_{\rm out}$. The charge
density $\sigma_{2}$ is given similarly, but with a negative sign.
Substituting the charge densities to eq. (2) and taking the difference gives
an integral equation for $\Delta V(r)$. The dependence on the unknown,
$\tilde{V}_{1}-\tilde{V}_{2}$, can be factored out by introducing an
auxiliary, $f(r)\equiv 1-\Delta V(r)/(\tilde{V}_{1}-\tilde{V}_{2})$, which
satisfies

$f(r)=\frac{\epsilon_{{\rm r},1}^{-1}+\epsilon_{{\rm
r},2}^{-1}}{2\pi}\int_{0}^{\infty}\\!\\!{\rm d}r^{\prime}\
\frac{1-f(r^{\prime})}{h+r^{\prime 2}/a}\frac{4r^{\prime}}{r+r^{\prime}}K(x).$
(4)

Here the integral over the azimuthal angle has been replaced with the complete
elliptic function of the first kind $K(x)$, where
$x\equiv\frac{4rr^{\prime}}{(r+r^{\prime})^{2}}$. When
$\epsilon_{1}=\epsilon_{2}$, eq. (4) reduces to Batchelor’s original result,
eq. (4.5) in ref.16. Since $f(r)$ is proportional to the contribution from the
surface charges in the contact region, we see that it is dominated by
particles of lower permittivity. When both permittivities approach infinity,
we have $f(r)=0$ and $\Delta V(r)=\tilde{V}_{1}-\tilde{V}_{2}$, which is
identical to the expression for the conductors discussed above. More general
cases are discussed below.

The solution to eq. (4) is uniquely determined by the normalized distance
$h/a$ and the average permittivity, $\epsilon_{\rm r}\equiv 2/(\epsilon_{\rm
r,1}^{-1}+\epsilon_{\rm r,2}^{-1})$. As shown by Batchelor,16 eq. (4) can be
non-dimensionalized to

$f(\rho)=\frac{1}{\pi}\int_{0}^{\infty}\\!{\rm
d}\rho^{\prime}\,\frac{1-f(\rho^{\prime})}{\lambda+\rho^{\prime
2}}\frac{4\rho^{\prime}}{\rho+\rho^{\prime}}K(x)$ (5)

in which $\rho\equiv r\epsilon_{\rm r}/a$ and $\lambda\equiv h\epsilon_{\rm
r}^{2}/a$. The numerically solved $f(\rho)$ for a few representative $\lambda$
values are shown in Fig. 2a. For large gap distance, with $\lambda\gg 1$,
$f(\rho)$ is nearly uniform, as expected. For smaller $\lambda$ values,
$f(\rho)$ decreases from $f(0)$ with $r$ monotonically. The difference in
electric potential at the contact points is proportional to $1-f(0)$. The
value of $f(0)$ increases as $h$ decreases, and reaches unity at $h=0$,
ensuring that the surface potential is continuous at the contact point. The
variation of $f(\rho)$ obtained from the above local analysis is confirmed by
directly solving the full potential distribution for dielectric dimers (inset,
Fig. 2a).

Similar to eq. (1) for the conductor case, the singular part of surface
charges on particle $1$ is given, with $R_{0}$ being the regularizing cutoff
of order $a$, by

$Q_{\rm s,1}=\frac{(\tilde{V}_{1}-\tilde{V}_{2})\epsilon_{\rm
out}}{4\pi}\left(1-\epsilon_{{\rm r},1}^{-1}\right)\int_{0}^{R_{0}}\\!\\!{\rm
d}r\frac{2\pi r\left[1-f(r)\right]}{h+r^{2}/a}$ (6)

The dependence on the dielectric permittivity appears in the prefactor
$1-1/\epsilon_{{\rm r},1}$ and in $f(r)$. The singular surface charge $Q_{\rm
s,2}$ on the particle 2 is given analogously, but with a negative sign.
However, because of the dependence on the dielectric permittivity, $Q_{{\rm
s},2}$ and $Q_{{\rm s},1}$ do not add up to zero. The first term in the square
bracket gives the same $\ln(a/h)$ singular contributions as eq. (1). The
second term represents the correction due to dielectric screening.

Figure 2: (a) Variation of surface potential $f(\rho)$ along the radial
direction for $\lambda=0$, $0.5$, $1.0$, $10.0$ and $1000$. The inset shows
the normalized potential difference for spherical dimer with a minimum
separation $h=0.1a$ and parameters $(Q_{i},\epsilon_{i},a_{i})$: $(1,100,1)$
and $(-1,10,2)$. Dots are numerically solved by a spectral method in
bispherical coordinate 10 and the solid line is proportional to $1-f(r)$. (b)
The singular capacitance $c_{1}$ is presented as a function of normalized
surface separation $h/a$ and relative permittivity $\epsilon_{\rm r}$.

The singular capacitance defined by $c_{1}(h)\equiv Q_{{\rm
s},1}/(\tilde{V}_{1}-\tilde{V}_{2})$ also has these two contributions. In the
non-dimensionalized form, it reads

$\displaystyle c_{1}(h)$ $\displaystyle=\frac{1-\epsilon_{{\rm
r},1}^{-1}}{4}a\epsilon_{\rm
out}\left(\ln\frac{a}{h}-P\left(\lambda\right)\right),\quad$ (7)
$\displaystyle P(\lambda)$ $\displaystyle\equiv\int_{0}^{\infty}\\!{\rm
d}\rho\,\frac{2\rho}{\lambda+\rho^{2}}f(\rho).$

In the last term, the upper bound is set to infinity because $f(\rho)$ decays
rapidly (as $\ln\rho/\rho$, ref.16). The dielectric correction is contained in
the term $P(\lambda)$, which depends on permittivity $\epsilon_{\rm r}$ and
relative gap distance $h/a$ through the combination $\lambda=h\epsilon_{\rm
r}^{2}/a$. The behavior of $P(\lambda)$ is derived from that of $f(\rho)$. For
$\lambda\gg 1$, $P(\lambda)$ is vanishingly small, so the capacitance is
dominated by $\ln(a/h)$, the characteristic behavior of conductors. On the
other hand, when $\lambda\to 0$, the leading contribution to $P(\lambda)$ is
$-\ln\lambda$, which cancels exactly the $\ln(a/h)$ dependence, leaving a term
$2\ln\epsilon_{\rm r}$ that diverges instead with the average permittivity.
Further, we note that $c_{1}(h)$ approaches the conductor limit for
$\epsilon_{\rm r}\gg 1$.

The difference in contact charges between dielectric and conductor cases is
solely contained in $P(\lambda)$. For intermediate separation, $c_{1}(h)$
shows a singular $\ln(a/h)$ dependence that is similar to the conductor
behavior. However, as long as $\epsilon_{\rm r}$ is finite, this logarithmic
behavior will eventually be cut off by contribution from $P(\lambda)$ at
sufficiently small separation. Therefore, unlike the conductor case, the
contact capacitance for dielectrics is finite at $h=0$. Instead, it approaches
a constant proportional to $2\ln\epsilon_{\rm r}$. Physically, the logarithmic
$h$-dependence originates from the accumulated polarization charges at the
interface. It is cut off for dielectrics because the potential difference,
which gives rise to the polarization charge, is self-consistently determined
by the latter. The weaker polarization effect for dielectrics eventually can
not keep up with the needed potential difference for producing the
polarization charge. The variation of $c_{1}$ on $h$ and $\epsilon_{\rm r}$
are shown in Fig. 2b. The crossover can be estimated by setting
$\lambda=h\epsilon_{\rm r}^{2}/a=1$. The logarithmic divergence, although
being cut at small gap separation, still plagues numerical calculations in
practice.18

This type of crossover behavior has been confirmed in a study on dielectric
spherical dimers using the image method.19 The limiting value of the
capacitance $2\ln\epsilon_{\rm r}$ as well as the crossover is also consistent
with the analytically known result between a conducting sphere and a
dielectric plane.20 However, unlike these two work on dimer particles,
Batchelor’s result is strictly local, suggesting that the same type of
singularity as represented by eq. (7) is applicable to all contact region in
cluster of multiple dielectric particles, irrespective of the particle shape.
In the following, we show that, the approach for isolating the singularity for
contact between conductors can be generalized to treat the dielectric
clusters, yielding the exact contact energy and distance-dependence with
modest numerical cost.

We consider the cluster consisting of $n$ dielectric spheres. Given the free
charge distribution $\rho\lower 3.0pt\hbox{${}_{\rm f}$}({\bf r})$, the
potential $\phi$ is governed by Poisson’s equation $\nabla\cdot\epsilon({\bf
r})\nabla\phi=-{\rho\lower 3.0pt\hbox{${}_{\rm f}$}({\bf r})}/{\epsilon_{\rm
out}}$, where $\epsilon_{\rm out}$ is the medium permittivity. The (relative)
material permittivity $\epsilon({\bf r})$ is set to $\epsilon_{\rm
r}\equiv\epsilon_{\rm in}/\epsilon_{\rm out}$ inside particles and unity in
the medium. Although the general charge distribution poses no added
difficulty, for simplicity, we focus on the case when only the free surface
charge $\sigma_{\rm f}$ are present. In such case, the system energy can be
written as surface integrals of the product of the free surface charge and the
surface potential over all surfaces. Furthermore, when the surface charge
distribution is uniform, the energy reduces to the sum of the product of the
total charge and the average surface potential. Therefore, we denote the set
of net charges on particles by ${{\bf Q}}$, and the mean surface potentials
${{\bf V}}$. The net charges and the mean potentials are linearly related,
i.e.,

${\bf Q}={\bf C}{\bf V},$ (8)

where ${\bf C}$ is an $n{\times}n$ “capacitance array”. In this notation, the
total energy can be expressed as $\mathcal{E}=\frac{1}{2}{{\bf Q}}\cdot{\bf
C}^{-1}\cdot{\bf Q}$. For convenience, the total energy $\mathcal{E}$
presented below is always normalized by the self-energy of a single isolated
sphere with the radius $a$ and a net charge $q$, $q^{2}/(8\pi\epsilon_{\rm
out}a)$. We note that this formulation applies to arbitrary particle shapes
and cluster configurations. If the surface charge distribution is non-uniform
or any external excitation exists, a multipole expansion of surface charge in
terms of dipole, quadrupole etc. is needed. The constitutive relation eq. (8)
and the quadratic expansion to energy can be generalized to include the
contribution from these higher multipoles and external electric fields. For
the purpose of demonstrating how the contact singularity is isolated, we focus
on the case of the uniform surface charge distribution.

Our method for dielectrics is an extension to our earlier work on conductors.5
Because the conductors are equipotential, the mean potential is also the
surface potential at every point on the surface. Since the capacitance $\bf C$
in eq. (8) contains two types of contributions, one type dominated by close
contacts between neighboring particles, and the other type from the remaining
interactions from all particles, we decompose the capacitance $\bf C$ as
follows

${\bf Q}=\left(c_{\rm s}{\bf L}+{\bf H}\right){\bf V}.$ (9)

Here, $c_{\rm s}(h)=\frac{a\epsilon_{\rm out}}{4}\ln(a/h)$ is the singular
capacitance for conductors derived above (assuming that all gap distances are
$h$). Since $c_{\rm s}$ becomes significant only for small gaps, the entries
in the array $\bf L$ are nonvanishing only for closely neighboring particles.
Specifically, if particle $i$ and $j$ form a close contact, we set the entries
$L_{ij}=L_{ji}=-1$. The diagonal entry $L_{ii}$ equals the number of close
neighbors of the particle $i$. All other entries of $\bf L$ are set to zero.
The singular capacitance $\bf L$ as defined here is the same as the adjacency
matrix representing the connectivity of clusters (see ref.5 for explicit
examples). In practice, for a given cluster configuration, we first
numerically solve the full capacitance array $\bf C$ at finite but small $h$
values, then subtract from $\bf C$ the singular term $c_{\rm s}(h){\bf L}$, to
get the regular capacitance ${\bf H}$. The $h$-dependence of this regular
capacitance is then fitted to a straight line when $h$ is small, allowing us
to obtain the full $h$-dependence for the capacitance array $\bf C$ and,
consequently, the full $h$-dependence of energy, down to $h=0$.

Generalizing eq. (9) to dielectrics requires two nontrivial modifications.
First, the decomposition in eq. (9) is valid because the potential of the
conductors can be used to evaluate the contact potential difference. But
dielectrics are not equipotential, and the surface potential at the contact
points $\tilde{V}_{i}$ generally differ from the average potential $V_{i}$ by
a numerical factor that depends on the dielectric permittivity and the cluster
configuration. Its variation with gap distance $h$ is weak, and approaches a
constant as $h\to 0$. So we generalize eq. (9) to

${\bf Q}=\left(\hbox{$\bf\tilde{\phantom{L}}$\kern-5.50003pt\hbox{$\bf
L$}}+{\bf H}\right){\bf V},\quad\tilde{L}_{ij}\equiv\frac{1-\epsilon_{{\rm
r},j}^{-1}}{V_{j}}\begin{cases}\tilde{V}_{j}^{(i)}c_{ij},&i\neq j;\\\\[5.0pt]
{\sum_{k=1}^{N}}\tilde{V}_{j}^{(k)}c_{jk},&i=j.\end{cases}$ (10)

Above, the superscript ‘$(i)$’ in $\tilde{V}_{j}^{(i)}$ indicates that the
contact potential is evaluated on the particle $j$ at the contact formed with
the particle $i$. The singular capacitance $c_{ij}$ is given from eq. (7) by
$c_{ij}=\frac{\epsilon_{\rm
out}a_{ij}}{4}\left(\ln\frac{a_{ij}}{h_{ij}}-P(\lambda_{ij})\right)$, which
depends on the mean radius of curvature $a$, gap distance $h$, and $\lambda$
value evaluated for the particle pair $i$ and $j$. From the definition, it is
clear that $\bf\tilde{\phantom{L}}$$\bf L$ is generally not symmetric, which
however reduces to the (symmetric) form identical to eq. (9) in the conductor
limit, since $\tilde{V}_{j}^{(i)}=V_{j}$ for all contacts on the particle $j$.

Second, it turns out that the contact potentials exhibit strong dependence on
the second order and more distant neighbors, because the dielectric screening
is comparatively weak (see Fig. 4 below). Therefore, to ensure rapid
convergence of the correct contact potentials, our construction of the
adjacency array $\bf\tilde{\phantom{L}}$$\bf L$ contains all pairs of
particles. Fortunately, as we shall see below, this construction requires no
extra computation other than evaluating the surface potentials at additional
contact points.

In order to identify the regular part ${\bf H}$ using eq. (10), we numerically
compute the full capacitance matrix ${\bf C}(h)$ and all the contact
potentials $\tilde{V}_{j}^{(i)}$, for a range of small but nonzero $h$ values.
We used the boundary element method (BEM) implemented in the package COPSS 18
to solve Poisson’s equation. The solution is expressed in terms of the induced
bound (polarization) charges $\sigma_{\rm b}$, which satisfies the following
boundary integral equation,

$\frac{1+\epsilon_{\rm r}}{2}\sigma_{\rm b}+(1-\epsilon_{\rm r})\epsilon_{\rm
out}\left({\bf E}_{\rm b}+{\bf E}_{\rm f}\right)\cdot\hat{\bf
n}=\frac{1-\epsilon_{\rm r}}{2}\sigma_{\rm f}.$ (11)

Above, the electric field strength at the surface due to induced bound charge
$\sigma_{\rm b}$ and free charges $\sigma_{\rm f}$ are expressed as surface
integrals ($\alpha={\rm b,f}$),

${\bf E}_{\alpha}=\frac{1}{4\pi\epsilon_{\rm out}}\int\\!{\rm
d}S^{\prime}\,\frac{{\bf r}-{\bf r}^{\prime}}{|{\bf r}-{\bf
r}^{\prime}|^{3}}\sigma_{\alpha}({\bf r}^{\prime}).$ (12)

By our convention, the surface normal $\hat{\bf n}$ points outward. For
clusters of multiple particles, it is understood that different particles may
have different values of permittivity $\epsilon_{\rm r}$. Equation (11) is
linear in $\sigma_{\rm b}$ and $\sigma_{\rm f}$. After discretization, the
bound charge is obtained by inverting the coefficient array. The apparent
divergent diagonal entries while discretizing the surface integral eq. (12)
can be re-parameterized to reproduce the correct self-energy. Further
numerical details can be found in ref.8. In all calculations by BEM presented
in this work (the dimer example used the exact series expansion.10), we meshed
each spherical surface into $17342$ triangular patches such that mesh size is
about $0.03\,a$.

In general, for a cluster of $n$ particles, $n$ separate numerical
calculations with $n$ independent charge vectors $\bf{Q}$ are needed to
determine the full capacitance matrix ${\bf C}$. The computational cost can be
reduced by imposing the symmetry of the cluster configuration. Subtracting
from $\bf C$ the singular contribution $\bf\tilde{\phantom{L}}$$\bf L$ is
expected to give the $h$-dependence of the regular part $\bf H$. However,
unlike the conductor case, where the coefficient to the logarithmic term is
exactly known, the singular capacitance $\bf\tilde{\phantom{L}}$$\bf L$ for
dielectrics depends on the calculated contact potentials
$\tilde{V}_{j}^{(i)}$, which is susceptible to the numerical precision and
what is meant by ‘contact point’ on a surface mesh. In practice, we vary the
gap distance $h$ and compute a few trial values of contact potentials, then
select the one that ensures the difference $H_{ij}(h)=C_{ij}-\tilde{L}_{ij}$
converges linearly with $h$ as $h\to 0$ for all the entries. In the following,
three examples are presented to demonstrate the nature of the contact
singularity, and to show that the regularization scheme allows us to obtain
the energy of cluster of dielectric particles with small $h$.

The first example is a dimer of identical dielectric spheres. It is the
simplest example for demonstrating the effects of interfacial polarization.6,
9, 10, 21 Except a study on the interaction between a conducting sphere and a
dielectric plane,20 very few past work tried to evaluate the energy at small
separation, presumably due to the difficulty of resolving the aforementioned
contact singularity. Therefore, we analyzed the polarization effect for
dielectric dimers in close contact, verified the singular behavior in eq. (7),
and showed that eq. (10) yields the full $h$-dependence of energy. We studied
both symmetric (${\bf Q}=(q,q)$) and asymmetric (${\bf Q}=(q,-q)$) cases, for
a range of permittivity values ($1\leq\epsilon_{\rm r}\leq 10^{6}$). The
energy at $h=0$ is presented as a function of relative permittivity, which
agrees with the exact result (see SI) found using the tangent-sphere
coordinate.

Figure 3: Electrostatic energy of a pair of identical spheres. (a) and (b):
asymmetric charge ${\bf Q}=(q,-q)$. (c) and (d): symmetric charge ${\bf
Q}=(q,q)$. In (b) and (d), dots are numerical results obtained using our
proposed method, and curves are exact prediction obtained using the tangent-
sphere coordinate (SI).

For the asymmetric case, interfacial polarization enhances the inter-particle
attraction. To illustrate this effect, we subtract the energy of two isolated
spheres from that of a dimer, then plot it against the gap distance in Fig. 3,
for three representative values of $\epsilon_{\rm r}$: $10$, $10^{2}$ and
$10^{4}$. In absence of the polarization effect, the energy curves would
exhibit no dependence on $\epsilon_{\rm r}$, and appear flat over this narrow
range of $h/a$ values. As $\epsilon_{\rm r}$ is increased from $1$, we confirm
the expected, stronger distance dependence, as $h\to 0$. In particular, for
the case of $\epsilon_{\rm r}=10^{4}$, an apparent logarithmic dependence on
$h$ is seen in Fig. 3a, which is precisely the singularity revealed in eq. (7)
and is stronger for higher permittivity values. On the other hand, when the
permittivity $\epsilon_{\rm r}$ is decreased, this logarithmic dependence is
cut off at an increasingly larger separation and becomes almost invisible when
$\epsilon_{\rm r}=10$.

The normalized energy difference, as shown in Fig. 3a, evaluated at $h=0$, is
show in Fig. 3b for different $\epsilon_{\rm r}$ values. This contact energy
represents the energy gain for bringing two particles from infinity to
contact. In terms of the normalizing energy unit, $q^{2}/(8\pi\epsilon_{\rm
out}a)$, its value is $-1$ without the polarization effect, and becomes more
negative as $\epsilon_{\rm r}$ increases, eventually reaching $-2$ at
$\epsilon_{\rm r}=\infty$.5 The contribution from the polarization effect for
$\epsilon_{\rm r}\gtrsim 100$ is comparable to the coulombic attraction alone.
Moreover, we notice a slow convergence of contact energy for permittivities
$\epsilon_{\rm r}\gtrsim 10^{3}$, owing to the weak $\ln\epsilon_{\rm r}^{2}$
dependence in the contact capacitance eq. (7) at $h=0$. Finally, we affirm our
results obtained from eq. (10), by comparing the contact energies to the exact
analytical results (SI).

For the symmetric case, the interfacial polarization weakens the inter-
particle repulsion, which is seen from the $h$-dependence of dimer energy.
However, by symmetry, the potential difference in the gap region vanishes and
thus no logarithmic behavior is observed in Fig. 3c. The energy scales
linearly with $h$ and a straightforward extrapolation gives the contact energy
for all permittivity values plotted in Fig. 3d. A much weaker polarization
effect is found. Even for the conductor case, only about $10\%$ contact energy
can be attributed to polarization effects. The contact energies converge for
$\epsilon_{\rm r}\gtrsim 100$, and reach $2(1/\ln 2-1)\simeq 0.88$, in
agreement with Maxwell’s result22 for the conducting dimers. As for the
asymmetric case, the full $\epsilon_{\rm r}$-dependence matches the analytical
calculation (SI). To conclude this example, we note that the symmetric case is
a peculiar example where the contact singularity is strictly absent. For all
other charge ratios, using the decomposition in eq. (10) to isolate the strong
$h$-dependence is necessary.

The second example concerns the energy of $8$ identical spheres placed at the
vertices of a cube, which illustrates the importance of the second order
contact for dielectric clusters. As discussed above, the singular capacitance
$\bf\tilde{\phantom{L}}$$\bf L$ in eq. (10) contains not only contributions
from the nearest neighbors, but also those from the second-order and higher
order neighbors. Therefore, a sphere at a vertex of a cube forms a secondary
contact with its $3$ plane diagonal vertices, and a third order contact with
the body diagonal vertex. Even though the higher order contributions are
minor, because at such large separation, these contacts cease to be singular,
keeping the second order contribution is essential for obtaining the correct
linear scaling of the regular capacitance with gap distance $h$ shown in Fig.
4a.

As in the dimer case, these regularized capacitance allow us to evaluate the
energy at arbitrarily small gap distance. Figure 4b shows the normalized
energies when only one particle is charged, for which all the energetic
contributions come from the interfacial polarization. Comparing the results
from BEM and our regularization schemes containing varying order of contact
contributions, it is clear that keeping the higher order singular contribution
is necessary for correctly evaluating the energy for $h/a\lesssim 0.05$. In
contrast, no such terms are needed for the conducting case.

Figure 4: (a) Variation of entry $H_{11}$ in the regular capacitance array
against separation. The results for the other entries are provided in Fig. S3.
(b) Variation of electrostatic energy against separation for dielectric
spheres with $\epsilon_{\rm r}=100$ placed on cubic vertices, where one sphere
is charged.

The third example is our main result on the energy of clusters, from which the
cohesive energy is obtained by subtracting the self-energy of particles at
infinite separation. We considered two limiting configurations: the most
extended one with all spheres arranged along a straight line (string), and the
most compact one with all spheres posited at the vertices of the platonic
solids (polyhedron). In our earlier work on the conducting spheres, the
charges are allowed to flow freely between contacting spheres. The energy of
polyhedron packing is found to be lower than that of the string packing at a
finite separation. However, as the gap distance decreases, the polarization
effects due to contact singularity become increasingly relevant, which
ultimately makes the string packing to be energetically more favorable than
the polyhedron packing. For the dielectric cases, we show that the weakened
contact singularity causes another crossover between the relative stability of
string-like and polyhedral packings.

Figure 5: Size dependence of electrostatic energy $\mathcal{E}$ for string-
like and polyhedral configurations, with fixed total charge. (a) Charge
transfer is permitted. Dashed curve: $\mathcal{E}=(2/n)^{1/3}/(2\ln 2)$, the
energy for polyhedral configurations of conducting spheres.5 Solid curve:
$\mathcal{E}=\ln(2na/r_{0})/[n\ln 2\,\ln(2a/r_{0})]$, the energy for a
conducting cylinder of the same volume with length $L=2na$ and radius
$r_{0}=0.816\,a$.23 (b) Charge transfer is prohibited. The total charge $q$
resides on the middle particle of the string (string-$2$) and an arbitrary
vertex of polyhedron (polyhedron).

As the previous example, to highlight the polarization effects, a unit charge
$q$ is placed on one sphere. For symmetric polyhedron packing, this charged
sphere is arbitrarily chosen. For the string packing, we considered two
extreme placings: at one end (string-1) and in the middle (string-2). For each
particle configuration and charge assignment, we considered two scenarios:
charge transfer between contacting particles is permitted (Fig. 5a) and
prohibited (Fig. 5b).

When inter-particle charge transfer is permitted, while maintaining uniform
charge distribution on each individual particle, the energy $\mathcal{E}$ can
be calculated using ${\bf Q}\cdot{\bf C}^{-1}\cdot{\bf Q}/2$. Minimizing this
energy subject to the constraint of constant total charge $q$, the optimal
charge assignment is found to be ${\bf Q}_{e}=q\,{\bf C}\cdot{\bf u}/C_{\rm
e}$, where ${\bf u}\equiv(1,1,\cdots,1)$ and $C_{\rm
e}\equiv\sum_{i,j=1}^{n}C_{ij}$. This implies that the optimal charge ${\bf
Q}_{e}$ produces identical average surface potentials for all the particles,
which generalizes the ‘equipotential’ concept for conductors.5
Correspondingly, the minimized energy is given by $q^{2}/(2C_{\rm e})$, where
the mean capacitance $C_{\rm e}$ weakly depends on the dielectric
permittivity. In this scenario, there is no need to differentiate string-$1$
and string-$2$ configurations. The energies of string and polyhedron
configurations at $h=0$ are plotted in Fig. 5a, for $\epsilon_{\rm r}=10$ and
$500$. For all configurations, the energy decreases with the cluster size $n$
monotonically since the redistributed surface charges are further apart. The
dependence on permittivity is surprisingly weak, and the size-dependence is
nearly indistinguishable from that of the conductor case.5 Two curves are
obtained by treating the cluster as a single conductor respectively, whose
capacitance scales with the length scale of the cluster, i.e., $n^{1/3}$ for
the polyhedron packing and $n/\ln(n)$ for the string packing.5 The extended
string packing has a lower cohesive energy because its effective capacitance
is higher than that of the compact polyhedron packing.

More interesting behavior is found (Fig. 5b) when charge transfer is
prohibited. The dependence on permittivity is evident for all three cases:
string-$1$, string-$2$, and polyhedron. To better visualize energies of
string-$2$ and polyhedron packing, energies of string-$1$ is presented in Fig.
S4) The energy is lower for the cluster with higher permittivity, because the
polarization effect is stronger. The energy of string-$2$ is lower than that
of string-$1$ for identical $n$, because the charged particles in the middle
of the string can polarize the particles in two half strings. Further, the
energies for both string-$1$ and string-$2$ configurations saturate as the
cluster size grows beyond $n=12$, because the polarization effect is short-
ranged. To assess the relative stability of compact or extended
configurations, we then only need to focus on the string-$2$ and the
polyhedron cases.

Our results indicate that the relative energetic stability depends on the
permittivity. For $\epsilon_{\rm r}=10$, the energy of the polyhedron packing
is lower than the string-2 packing, which is opposite to the characteristic
result for conductors shown in Fig. 5a. Therefore, we expect a crossover from
stable compact packing to stable open packing at intermediate permittivity
values. This is indeed the case shown for $\epsilon_{\rm r}=500$, whereby the
energies of compact and open configurations closely trace each other, and the
relative energetic stability changes at $n=4$ and $n=8$ respectively. As
$\epsilon_{\rm r}$ is further increased, the energy curves of corresponding
configuration will converge to those in Fig. 5a, reversing the relative
stability of open and compact packings.

In summary, we generalized our previous work on conducting particles,5 and
developed a scheme to resolve the singular contact charges between touching
dielectric spheres, which regularizes the full capacitance array by isolating
the singular contributions, i.e., eq. (10). Using this scheme, we obtain the
cohesive energies for dielectric clusters at zero separation containing up to
$n=20$ particles, which is difficult to resolve with brutal force numerical
calculations. Our results show that the shape of stable clusters formed from
dielectric particles depends on the permittivity ratio $\epsilon_{\rm r}$:
open clusters is more stable for large $\epsilon_{\rm r}$, and compact
clusters is more stable for small $\epsilon_{\rm r}$. Our scheme is applicable
to systems with arbitrary packing geometry, charge distribution, and
containing asymmetric interfaces that have different permittivities and radii
of curvature on the contacting particles. At the heart of our scheme is the
contact singularity, which resembles those encountered in the study of thermal
conduction 16 and momentum transport across a narrow gap.24

## References

  * 1 J. Kolehmainen, A. Ozel, Y. Gu, T. Shinbrot, and S. Sundaresan.  Phys. Rev. Lett., 121:124503, 2018.
  * 2 V. Lee, S. R. Waitukaitis, M. Z. Miskin, and H. M. Jaeger.  Nature Phys., 11:733, 2015.
  * 3 E. V. Shevchenko, D. V. Talapin, N. A. Kotov, S. O’Brien, and C. B. Murray.  Nature, 439:55, 2006.
  * 4 J. W. Merrill, S. K. Sainis, and E. R. Dufresne.  Phys. Rev. Lett., 103:138301, 2009.
  * 5 J. Qin, N. W. Krapf, and T. A. Witten.  Phys. Rev. E, 93:022603, 2016.
  * 6 A. Russell.  P. R. Soc. Lond. A-Conta., 82:524, 1909.
  * 7 J. Lekner.  P. Roy. Soc. A-Math. Phy., 468:2829, 2012.
  * 8 K. Barros, D. Sinkovits, and E. Luijten.  J. Chem. Phys., 140:064903, 2014.
  * 9 E. B. Lindgren, A. J. Stace, E. Polack, Y. Maday, B. Stamm, and E. Besley.  J. Comput. Phys., 371:712, 2018.
  * 10 H. Lian and J. Qin.  Mol. Syst. Des. Eng., 3:197, 2018.
  * 11 J. Qin, J. J de Pablo, and K. F. Freed.  J. Chem. Phys., 145:124903, 2016.
  * 12 M. Shen, H. Li, and M. O. De La Cruz.  Phys. Rev. Lett., 119:138002, 2017.
  * 13 K. Barros and E. Luijten.  Phys. Rev. Lett., 113:017801, 2014.
  * 14 Z. M. Sherman, D. Ghosh, and J. W Swan.  Langmuir, 34:7117, 2018.
  * 15 L. G. Leal.  Advanced Transport Phenomena: Fluid Mechanics and Convective transport Processes.  Cambridge University Press, 2007.
  * 16 G. K. Batchelor and R. W. O’Brien.  P. Roy. Soc. A-Math. Phy., 355:313, 1977.
  * 17 O. D. Kellogg.  Foundations of Potential Theory.  Courier Corporation, 1953.
  * 18 X. Jiang, J. Li, X. Zhao, J. Qin, D. Karpeev, J. Hernandez-Ortiz, J. J. de Pablo, and O. Heinonen.  J. Chem. Phys., 145:064307, 2016.
  * 19 L. Poladian.  Q. J. Mech. Appl. Math., 41:395, 1988.
  * 20 S. V. Kalinin, E. Karapetian, and M. Kachanov.  Phys. Rev. B, 70:184101, 2004.
  * 21 J. D. Love.  Q. J. Mech. Appl. Math., 28:449, 1975.
  * 22 J. C. Maxwell.  A Treatise on Electricity and Magnetism.  Oxford: Clarendon Press, 1873.
  * 23 J. D. Jackson.  Classical Electrodynamics.  John Wiley & Sons, 1999.
  * 24 R. H. Davis, J. A. Schonberg, and J. M. Rallison.  Phys. Fluids. A-Fluid., 1:77, 1989.