# Positive-Energy Dirac Particles and Dark Matter

Eugene Bogomolny  Université Paris-Saclay, CNRS, LPTMS, 91405 Orsay, France

###### Abstract

The relativistic positive-energy wave equation proposed by P. Dirac in 1971 is
an old but largely forgotten subject. The purpose of this note is to speculate
that particles described by this equation (called here Dirac particles) are
natural candidates for the dark matter. The reasoning is based on a fact that
the internal structure of such particles simply prohibits their interaction
with electromagnetic fields (at least with the minimal coupling) which is
exactly what is required for dark matter. Dirac particles have quite unusual
properties. In particular, they are transformed by an infinite-dimensional
representation of the homogeneous Lorentz group, which clearly distinguishes
them from all known elementary particles described by finite-dimensional
representations and hints to a physics beyond the Standard Model. To clarify
the topic, a brief review of the main features of the above-mentioned Dirac
equation is given.

##  I Introduction

A tremendous number of various astrophysical observations collected during the
last century are explained by postulating the existence of invisible (i.e.,
not interacting with electromagnetic fields) matter. Such a hypothetical
substance goes under the name of dark matter, and according to different
estimates, it constitutes $\sim 85\%$ of the total matter content of the
universe (see, e.g., reviews supersymmetricDM ; particleDM ; nature  and
references therein). Though dark matter is a predominant ingredient of the
universe, its composition remains unknown. A great number of different
candidates have been proposed supersymmetricDM ; particleDM ; nature , but
none have of course been confirmed so far.

The purpose of this short note is to suggest a (seemingly unexplored)
possibility that dark matter is made up of particles obeying the so-called new
Dirac equation proposed by P. Dirac Dirac_I  40 years after his famous (old)
equation Dirac . Both equations are relativistic wave equations. The usual
(old) Dirac equation describes particles with positive energy as well as
antiparticles with negative energy and is one of the best known equations in
physics. The principal feature of the new equation is that all its solutions
have only positive energies (thus, antiparticles are absent or, at least, not
described by the same equation). For relativistic equations, such property is
unusual. The relativistic invariance implies that $E^{2}=\vec{p}^{\,2}+m^{2}$,
and in general, there exist solutions with positive as well as negative
energies: $E=\pm\sqrt{\vec{p}^{\,2}+m^{2}}$. Though today negative energy
solutions simply indicate the existence of antiparticles, the search for
relativistic equations with only positive energies continues to persist (see,
e.g., cirilo  and references therein).

It was E. Majorana majorana ; fradkin  who first constructed a relativistic
wave equation with only positive energy solutions. It appears that it becomes
feasible provided one considers particles transformed by an infinite-
dimensional representation of the homogeneous (i.e., without translations)
Lorenz group. Usual relativistic particles (scalars, spinors, vectors, etc.)
realise finite-dimensional representations of the Lorentz group, and their
equations necessarily have positive and negative energy solutions. The
positive-energy Majorana equation formally has form of the usual Dirac
equation:

$\big{(}\mathbf{\Gamma}^{\mu}\partial_{\mu}-m\big{)}\mathbf{\Psi}(x)=0$ (1)

but with the difference that its eigenfunctions $\mathbf{\Psi}(x)$ have an
infinite number of components (labeled by two non-negative integers) and
matrix $\mathbf{\Gamma}^{\mu}$ is an infinite-dimensional matrix derived from
the requirement of the relativistic invariance (whose explicit form could be
obtained from expressions in Section IV). This equation explicitly describes
an infinite number of particles with all integer and half-integer spins and
masses decreasing with spin: $m(j)=m/(j+1/2)$ majorana ; fradkin . Equation
(1) should not be confused with another much better-known Majorana equation
for self-conjugate particles with spin-$\tfrac{1}{2}$ usual_majorana .

The distinction between finite-dimensional and infinite-dimensional
representations of a group is that for the former, the action of group
generators mixes a finite number of wave components (e.g., four for the usual
spinors) but for the latter, such action requires an infinite number of
components. Mathematically, it is well established that the homogenous Lorentz
group has both finite and infinite-dimensional representations (see, e.g.,
harish ; gelfand  and references therein). All known elementary particles are
classified according to finite-dimensional representations of this group, yet
particles associated with its infinite-dimensional representations are not
usually taken into account.

In the new equation, Dirac found an elegant way (see below) to select one
state from an infinite tower of Majorana solutions. Particles described by the
new Dirac equation are also transformed by an infinite-dimensional
representation of the Lorenz group, which markedly differentiates them from
all known elementary particles. Such particles can be called positive-energy
Dirac particles, but for the sake of brevity, in this note, they are referred
to as Dirac particles.

Another characteristic feature of these Dirac particles is a nontrivial fact
that they cannot interact with electromagnetic fields (at least with the
minimal coupling). Such a no-interaction theorem is related with the internal
structure of the new Dirac equation and is not imposed by simply postulating
that the charge of Dirac particles equals zero. This propriety emerges as a
byproduct of Dirac’s construction, and in the past has been considered as a
serious defect of the new Dirac equation, but it is exactly what makes Dirac
particles a natural candidate for dark matter, which by definition is a
substance not interacting with electromagnetic fields.

In the literature, there are articles that relate Dirac particles and dark
matter (see, e.g., dirac_dark_matter , besides others), but in these papers,
the term ’Dirac particles’ indicates solely spin-$\tfrac{1}{2}$ particles
obeying the usual Dirac equation Dirac . Dirac particles discussed here are of
a different nature as they satisfy the new Dirac equation Dirac_I  and are
transformed by an infinite-dimensional representation of the Lorentz group. To
the author’s knowledge, such particles were not considered previously as
candidates for dark matter.

The purpose of the paper is to attract wider attention to this potentially
important subject and to serve as a brief introduction to Dirac particle
properties. As such issues are largely undeveloped, the discussion below is
restricted only to principal features of Dirac particles leaving aside many
important questions.

The plan of the paper is the following. In Section II, the main properties of
the new Dirac equation are discussed. Section III is devoted to the
investigation of the relativistic invariance of this equation and to the
meaning of the infinite-dimensional representation of its solutions. By
construction, the new Dirac equation is an overdetermined system of equations,
and consequently, one can derive additional equations that are indicated in
Section IV. Section V deals with the no-electromagnetic-interaction property
of the new Dirac equation, which is of primary importance for the conjectured
relation of Dirac particles with dark matter. In Section VI, the interaction
of Dirac particles with gravitational fields is shortly treated. Section VII
gives an overview of the discussed topics. For clarity, certain useful
formulas are presented in Appendix A, and for completeness, the tetrad
formalism is recalled in Appendix B.

##  II New Dirac Equation

The usual Dirac equation Dirac  describes relativistic fermions with spin
$\tfrac{1}{2}$ and mass $m$ and is one of the most important and the best-
known equations of mathematical physics. Much less attention has been
attracted to another Dirac work Dirac_I , written about 40 years after the
famous paper Dirac , in which a different relativistic equation (called the
new Dirac equation) was proposed and analysed.

The new equation has the functional form similar to the usual Dirac equation:

$\left(\gamma^{\mu}\partial_{\mu}-m\right)\Psi=0$ (2)

but with an important difference in that column-vector $\Psi$ depends on two
sets of variables: four space-time coordinates $x^{\mu}$ with $\mu=0,1,2,3$
and two ’internal’ (auxiliary) variables $q_{1}$ and $q_{2}$, arranged in the
following manner:

$\Psi(x,q)=Q\,\Phi(x,q)\,,\qquad Q=\left(\begin{array}[]{c}q_{1}\\\ q_{2}\\\
q_{3}\\\ q_{4}\end{array}\right)$ (3)

where components $q_{3}$ and $q_{4}$ are the momenta in the ’internal’ space

$q_{3}=-i\frac{\partial}{\partial q_{1}}\,,\qquad
q_{4}=-i\frac{\partial}{\partial q_{2}}\,.$ (4)

The choice of such column-vector $Q$ corresponds to the following commutation
relations:

$[q_{j},q_{k}]=i\beta_{jk}\,,\qquad\beta=\left(\begin{array}[]{cc}{0}&{1}\\\
-{1}&{0}\end{array}\right).$ (5)

As for the usual Dirac equation, $\gamma^{\mu}$ are $4\times 4$ matrices
obeying $\gamma^{\mu}\gamma^{\nu}+\gamma^{\nu}\gamma^{\mu}=2\eta^{\mu\nu}$
with the metric tensor of the Minkowski space
$\eta^{\mu\nu}=\mathrm{diag}(-1,1,1,1)$ (cf. (48)).

Function $\Phi(x,q)$ in (3) depends on both sets of coordinates (i.e. on six
variables), and it is the only one unknown function. Thus, the new Dirac
equation corresponds to an overdetermined system of four equations in one
unknown function:

$\hat{P}_{j}\,\Phi(x,q)=0\,,\qquad\hat{P}_{j}=(\gamma^{\mu}\partial_{\mu}-m)_{jk}q_{k}\,.$
(6)

To ensure the consistency of these equations, it is necessary that the
commutator of any two $\hat{P}_{j}$ equals zero. Using (5) one finds

$[\hat{P}_{j},\,\hat{P}_{k}]=\big{(}(\gamma^{\mu}\partial_{\mu}-m)i\beta(\gamma^{\mu\,T}\partial_{\mu}-m)\big{)}_{jk}$
(7)

where $\gamma^{\mu\,T}$ are the transposed gamma-matrices. In the new Dirac
equation, these matrices have to be chosen to fulfil the requirement

$\beta\,\gamma^{\mu\,T}=-\gamma^{\mu}\,\beta\,.$ (8)

By construction, matrix $\beta$ is an antisymmetric matrix (cf. (5)), and this
requirement means that $\beta\gamma^{\mu}$ are symmetric matrices.

Under these conditions, the above commutator becomes

$[\hat{P}_{j}\,\hat{P}_{k}]=i\beta_{jk}\big{(}-\eta^{\mu\nu}\partial_{\mu}\partial_{\nu}+m^{2}\big{)}$
(9)

which implies that all four equations (6) are consistent iff $\Phi(x,q)$ as a
function of $x$ obeys the Klein–Gordon equation:

$\big{(}-\eta^{\mu\nu}\partial_{\mu}\partial_{\nu}+m^{2}\big{)}\Phi(x,q)=0\,.$
(10)

It is this additional relation that differs the new Dirac equation from the
Majorana equation (1) majorana  and permits us to select one state from an
infinite number of Majorana solutions.

It is plain that (10) is a necessary condition for the existence of a solution
of the new Dirac equation. It could be obtained by, e.g., the multiplication
of (2) with (3) by $Q^{T}\beta\,\big{(}\gamma^{\nu}\partial_{\nu}+m\big{)}$
and using the fact that $Q^{T}\beta\,Q=2i$ (see (51)). The sufficiently of
this equation is more subtle. Usually, one refers to the Frobebius theorem
(see, e.g., foliations  and references therein), which (in the simplest
setting) states that in order for a system of linear differential equations to
have a solution, it is necessary and sufficient that their commutators (more
general, Lie brackets) lie in the their span. The important ingredient of this
theorem is the fact that commutators cancel second derivatives, which is not
valid for the considered problem. Therefore, the statement that (10) is the
sufficient condition for the existence of solutions of the new Dirac equation
requires additional investigations.

It has been proved in Dirac_I  that these equations determine massive
relativistic particles (the Dirac particles) with only positive energy and
zero spin in the rest frame. In the next publication Dirac_II , Dirac showed
that when the noncommutativity of different variables is ignored, Dirac
particles could be considered as vibrating spherical shells similar to
Zitterbewegung of particles obeying the usual Dirac equation.

The simplest solution of the new Dirac equation is the plane wave of the
following form Dirac_I :

$\Phi(x,q)=C\exp\left(-\frac{1}{2}aq_{1}^{2}-\frac{1}{2}bq_{2}^{2}+cq_{1}q_{2}+ip_{\mu}x^{\mu}\right)$
(11)

where $p_{0}=\sqrt{\mathbf{p}^{2}+m^{2}}>0$, $C$ is a normalisation constant,
and coefficients $a,b,c$ are

$a=\frac{m+ip_{1}}{p_{0}+p_{3}}\,,\quad b=\frac{m-ip_{1}}{p_{0}+p_{3}}\,,\quad
c=\frac{ip_{2}}{p_{0}+p_{3}}\,.$ (12)

The positivity of Dirac particles’ energy is attested by the fact that
negative energy solutions with $p_{0}=-\sqrt{m^{2}+\mathbf{p}^{2}}$ lead to
non-normalised functions and have to be excluded.

##  III Infinite-Dimensional Representation

The relativistic (co)variance of the new Dirac equation means that a linear
Lorentz transformation $x^{\prime\mu}=\Lambda^{\mu}_{\nu}x^{\nu}$
($\eta^{\mu\nu}=\eta^{\sigma\rho}\Lambda^{\mu}_{\sigma}\Lambda^{\nu}_{\rho}$)
of a solution

$\Phi(\Lambda x,q)=U(\Lambda)\Phi(x,q)$ (13)

is also a solution of the same equations. Here, $U(\Lambda)$ is an operator
acting on $q$ variables.

The inspection of Equation (3) shows that it will be the case if the following
two conditions are fulfilled:

$U^{-1}\,Q\,U=V\,Q\,,\qquad V\,\gamma\,V^{-1}=\gamma\,\Lambda$ (14)

where matrix $V\equiv V_{jk}(\Lambda)$ acts on spinor indices. Indeed, under
such changes, Equation (3) transforms linearly; thus,

$\hat{P}\,\Phi(\Lambda x,q)=U\,V\,\hat{P}\,\Phi(x,q)=0\,.$ (15)

The explicit expressions of $V$ and $U$ can easily be obtained from the
infinitesimal Lorentz transformations

$\Lambda^{\mu}_{\nu}=\delta^{\mu}_{\nu}+a^{\mu}_{\nu}=\delta^{\mu}_{\nu}+\frac{1}{2}a_{cd}\big{(}I^{cd}\big{)}^{\mu}_{\
\nu}$ (16)

where $4\times 4$ antisymmetric matrix $a_{cd}$ describes six parameters of
the Lorentz transformation and matrices $I^{cd}$ are generators of Lorentz
transformation

$\big{(}I^{cd}\big{)}^{\mu}_{\
\nu}=\eta^{c\mu}\delta^{d}_{\nu}-\eta^{d\mu}\delta^{c}_{\nu}\,.$ (17)

These generators obey the commutator relations of the Lie algebra of the
Lorenz group

$\big{[}I^{cd},I^{kr}\big{]}=\eta^{dk}I^{cr}+\eta^{cr}I^{dk}-\eta^{dr}I^{ck}-\eta^{ck}I^{dr}\,.$
(18)

It is plain that $V$ and $U$ should be of the following form:

$V=\exp\Big{(}\frac{1}{2}a_{cd}\Sigma^{cd}\Big{)},\qquad
U=\exp\Big{(}\frac{i}{2}a_{cd}S^{cd}\Big{)}.$ (19)

In the first order of $a$, unknown $\Sigma$ and $S$ are determined from the
following commutators:

$[\Sigma,\gamma]=\gamma I,\qquad i[Q,S]=\Sigma Q$ (20)

whose solutions are

$\Sigma^{cd}=\frac{1}{4}\Big{(}\gamma^{c}\gamma^{d}-\gamma^{d}\gamma^{c}\Big{)},\qquad
S^{cd}=\frac{1}{2}\big{(}\,\overline{Q}\,\Sigma^{cd}\,Q\,\big{)}$ (21)

where $\overline{Q}=Q^{T}\beta$. Because $\beta\,\Sigma^{T}=-\Sigma\,\beta$
one has

$\beta V^{T}=V^{-1}\beta$ (22)

which implies that $Q^{\prime}=VQ$ also obeys the commutation relations (5)
and $U^{-1}\,\overline{Q}\,U=\overline{Q}\,V^{-1}$.

It is well known that $\Sigma^{cd}$ obey Equation (18), thus realising the
four-dimensional (spinor) representation of the Lie algebra of the Lorenz
group. As matrices $\beta\Sigma^{cd}$ are symmetric (cf. (5)), it follows from
(52) that spin operators $iS^{cd}$ from the above equation obey the same
commutator relations (18), i.e., they also realise a representation of the Lie
algebra of the Lorentz group.

The explicit form of these operators is presented in (• ‣ A). These operators
are a distinct feature of the new Dirac equation. When acting on an $L_{2}$
function $\Phi(x,q)$ depended on auxiliary variables $q_{1}$ and $q_{2}$
(i.e., $\int_{-\infty}^{\infty}|\Phi(x,q)|^{2}dq<\infty$ for arbitrary $x$),
they lead to an infinite-dimensional representation of the homogeneous Lorentz
group. A simple way to see it explicitly is to write the function as an
infinite series in suitable functions of $q_{1}$ and $q_{2}$. A convenient
choice is eigenfunctions of the harmonic oscillators
$\hat{H}_{1}=\frac{1}{2}(q_{3}^{2}+q_{1}^{2})$ and
$\hat{H}_{2}=\frac{1}{2}(q_{4}^{2}+q_{2}^{2})$ (the Hermite functions):
$\hat{H}_{i}\phi_{m}(q_{i})=(m+\frac{1}{2})\phi_{m}(q_{i})$ (cf. (66))

$\Phi(x,q)=\sum_{n,m=0}^{\infty}A_{mn}(x)\phi_{m}(q_{1})\phi_{n}(q_{2})\,.$
(23)

(Due to a symmetry of the new Dirac equation, integers $m$ and $n$ are of the
same parity.)

Coefficients $A_{mn}(x)$ form an infinite number of wave function components
of Dirac particles. The action of operators $S^{cd}$ on functions
$\phi_{m}(q_{1})$ and $\phi_{n}(q_{2})$ induces a transformation of $A_{mn}$
corresponded to an infinite (and unitary) representation of the homogeneous
Lorenz group (see, e.g., majorana ; fradkin , harish ; gelfand , and Dirac_III
). Its explicit form can be obtained from Equation (68). In a sense, an
infinite ladder of coefficients $A_{mn}$ is hidden in the dependence of an
eigenfunction on variables $q_{1}$ and $q_{2}$, and function $\Phi(x,q)$ is a
concise generating function of these coefficients.

Spin matrices $\Sigma^{cd}$ equal the commutator of two $\gamma$-matrices:
$\Sigma^{cd}=\tfrac{1}{4}[\gamma^{c},\gamma^{d}]$. Similarly, spin operators
$S^{cd}$ also can be written as the commutator of two operators (cf. (52),
(53)):

$S^{cd}=i\big{[}\Gamma^{c}\,,\Gamma^{d}\big{]},\qquad\Gamma^{c}=-\frac{1}{4}\big{(}\overline{Q}\gamma^{c}Q\big{)}.$
(24)

Operators $S^{cd}$ together with $\Gamma^{c}$ form a representation of the
group $SO(3,2)$ Dirac_IV ; lectures  (see also Appendix A).

From (13), it follows that

$\Phi(x,q)=U(\Lambda)\,\Phi(\Lambda^{-1}x,q)\,.$ (25)

This invariance leads to the conservation of the total momenta Dirac_I
$J^{cd}=M^{cd}-iS^{cd}$. $M^{cd}=x^{c}\partial^{d}-x^{d}\partial^{c}$ are the
usual momenta and $S^{cd}$ represent the spin momenta.

##  IV Additional Equations

The new Dirac equation (6) is an overdetermined system of equations and any
linear combination of these equations is also a valid equation; thus,
$M\hat{P}_{j}\,\Phi=0$ with any $4\times 4$ matrix $M$. Multiplying this
equality by $\overline{Q}$ and summing over spin indices, one obtains $16$
equations quadratic in $Q$:
$\overline{Q}M(\gamma^{\mu}\partial_{\mu}-m)Q\Phi=0$. When $M=\gamma^{5}$, one
obtains identical zero and the remaining $15$ equations (corresponded to
$M=1,\;\gamma^{\mu},\;\Sigma^{\mu\nu},\;\gamma^{5}\gamma^{\mu}$, respectively)
take the following form Dirac_I ; staunton :

$\displaystyle\big{(}\Gamma^{\mu}\partial_{\mu}+\frac{i}{2}m\big{)}\Phi=0\,,$
(26)
$\displaystyle\big{(}\frac{i}{2}\partial^{\mu}+S^{\mu\nu}\partial_{\nu}+m\Gamma^{\mu}\big{)}\Phi=0\,,$
(27)
$\displaystyle\big{(}\Gamma^{\mu}\partial^{\nu}-\Gamma^{\nu}\partial^{\mu}+mS^{\mu\nu}\big{)}\Phi=0\,,$
(28)
$\displaystyle\big{(}e_{\mu\nu\alpha\beta}S^{\nu\alpha}\partial^{\beta}\big{)}\Phi=0\,.$
(29)

The first of these equations is (up to a notation) the positive-energy
Majorana equation (1) majorana , and the last one coincides (up to a factor)
with the Pauli–Lubanski (pseudo)vector, indicating that field $\Phi$ has zero
spin.

The Majorana equation (26) is a characteristic feature of positive-energy
equations. All solutions of such equations, in a way or another, obey this
equation. Though the form of the Majorana equation is similar to the usual
Dirac equation (2), they lead to different conclusions. The main point is
that, contrary to the usual $\gamma^{0}$ matrix, which has positive and
negative eigenvalues (see (48)), operator $\Gamma^{0}$ has the form (cf. (55))

$\Gamma^{0}=\frac{1}{4}\Big{(}(-\partial_{q_{1}}^{2}+q_{1}^{2})+(-\partial_{q_{2}}^{2}+q_{2}^{2})\Big{)}.$
(30)

When acting on $L_{2}$ functions of $q_{1},q_{2}$ invariant over inversion, as
in (23), it has only positive eigenvalues equal to $\tfrac{1}{4}(2j+1)$ with
integer $j\geq 0$. It is this property that implies that particles described
by the Majorana equation have masses decreasing with the spin:
$m(j)=m/(j+1/2)$.

In Dirac’s approach, the above equations are just consequences of Equation
(6). It was argued in staunton  that Equation (27) may be chosen as the
fundamental one, leading to the same plane wave solution as in (11). A more
general equation has been considered in that reference:

$\Big{(}i\kappa\partial^{\mu}+S^{\mu\nu}\partial_{\nu}+m\Gamma^{\mu}\Big{)}\Phi=0$
(31)

with a parameter $\kappa$. It has been shown that this system of equations is
consistent only for two values $\kappa=\tfrac{1}{2}$ and $\kappa=1$. The first
one corresponds to the Dirac equation (27), and the second one gives a new
equation describing particles transforming by an infinite-dimensional
representation of the Lorentz group but with spin $\tfrac{1}{2}$ in the rest
frame. Contrary to the Dirac particles described in this note, such
spin-$\tfrac{1}{2}$ particles permit the interaction with electromagnetic
fields.

##  V No-Interaction Theorem

The principal goal of the new Dirac equation was the construction of a
relativistic equation with only positive energy solutions. But for the
relation with dark matter advocated here, the most important is another
unusual fact: that Dirac particles cannot interact with electromagnetic
fields.

Usually, the interaction of a particle with an electromagnetic field $A_{\mu}$
is introduced within the so-called minimal coupling by the substitution
$\partial_{\mu}\to D_{\mu}=\partial_{\mu}+ieA_{\mu}(x)$ into the wave
equation. The remarkable propriety of the new Dirac equation is that this
standard procedure leads to the inconsistency of resulting equations, which
means that that Dirac particles cannot interact with electromagnetic fields
(at least with minimal coupling). This fact had been briefly mentioned in the
very end of the Dirac paper Dirac_I  and has been discussed in detail in
biedenharn  (see also mukunda ). For completeness, that reasoning is presented
below.

Let us consider the new Dirac equation with minimal coupling:

$\big{(}\gamma^{\mu}D_{\mu}-m\big{)}Q\Phi(x,q)=0\,,\qquad
D_{\mu}=\partial_{\mu}+ieA_{\mu}(x)\,.$ (32)

Applying the operator $\big{(}\gamma^{\mu}D_{\mu}+m\big{)}$ to this equation,
one finds

$\Big{(}\big{(}\eta^{\mu\nu}D_{\mu}D_{\nu}-m^{2}\big{)}+ie\Sigma^{\mu\nu}F_{\mu\nu}\Big{)}Q\Phi(x,q)=0$
(33)

where $F_{\mu\nu}=\partial_{\mu}A_{\mu}-\partial_{\nu}A_{\mu}$ and
$\Sigma^{\mu\nu}=\frac{1}{4}(\gamma^{\mu}\gamma^{\nu}-\gamma^{\nu}\gamma^{\mu})$.

To remove the first term in (33), (which corresponds to the Klein–Gordon
equation with the minimal coupling), let us multiply (33) by an antisymmetric
matrix $G$ such that $\mathrm{Tr}(G\beta)=0$ and convolute the result with
$Q^{T}$. According to (51), it will cancel the Klein–Gordon term. Besides $16$
independent products of the above $\gamma$ matrices, there are five matrices
with such proprieties: $\beta\gamma^{5}\gamma^{\mu}$ and $\beta\gamma^{5}$.
Let $G_{j}$ be one of these five matrices. After the indicated
transformations, Equation (33) takes the form

$F_{\mu\nu}h^{\mu\nu}_{j}\Phi=0\,,\qquad
h^{\mu\nu}_{j}=Q^{T}G_{j}\Sigma^{\mu\nu}Q$ (34)

where $h^{\mu\nu}_{j}$ are certain operators constructed from ten possible
symmetric bilinear combinations of $q_{j}$. Using properties of
$\gamma$-matrices (cf. (56) and (57)), one obtains

$\big{(}\overline{Q}\gamma^{5}\gamma^{\lambda}\Sigma^{\mu\nu}Q\big{)}=\frac{1}{2}e^{\lambda\mu\nu\rho}\Gamma_{\rho}\,,\qquad\big{(}\overline{Q}\gamma^{5}\Sigma^{\mu\nu}Q\big{)}=\frac{1}{8}e^{\mu\nu\lambda\rho}S_{\lambda\rho}$
(35)

where $S^{\mu\nu}$ are antisymmetric spin operators introduced in (21)), and
$\Gamma^{\mu}$ is defined by (24).

Therefore, the above equations take the form

$\tilde{F}_{\mu\nu}\Gamma^{\nu}\Phi=0\,,\qquad\tilde{F}_{\mu\nu}S^{\mu\nu}\Phi=0\,.$
(36)

Here, $\tilde{F}_{\mu\nu}=\frac{1}{2}e_{\mu\nu}^{\hskip
6.02777pt\rho\delta}F_{\rho\delta}$ is the dual electromagnetic tensor.

In biedenharn , it was argued that from (36), it follows that if
$F_{\mu\nu}\neq 0$, all operators $\Gamma^{\mu}$ and
$S^{\mu\nu}=i\big{[}\Gamma^{\mu},\,\Gamma^{\nu}\big{]}$ annihilate
$\Phi(x,q)$, which means that $\Phi(x,q)$ is independent on $q$, and
consequently, the new Dirac equation with the minimal coupling of
electromagnetic fields has no nontrivial solution. In other words, Dirac
particles cannot interact with nonzero electromagnetic fields.

Practically all works on this subject have stressed this ’drawback’, and many
papers have been devoted to generalisations of the new Dirac equation, which
would permit the electromagnetic interaction. It appears that it can be
carried out only by considerable complications of the equation, e.g., by
considering higher-order momentum terms staunton ; biedenharn , or by
increasing the number of ’intrinsic’ variables chiang  or by introducing
parabosonic constituents lectures ; sudarshan .

The no-electromagnetic-interaction theorem has been considered in the past as
a serious disadvantage of the new Dirac equation, and to a great extent, it
was responsible for the loss of interest to this subject. But it is exactly
what is required for dark matter, which, by definition, is a substance whose
characteristic property is that it does not interact with electromagnetic
fields. Dirac particles cannot be connected to electromagnetic fields not by
simply postulating that they have zero charge but because their relativistic
equation is inconsistent in the presence of the fields. It is this property
that makes Dirac particles an ideal and natural candidate for the illusive
dark matter. In a sense, merely the existence of particles transforming by an
infinite-dimensional representation of the Lorenz group implies the existence
of dark matter.

##  VI Gravity and Dirac Particles

The interaction of elementary particles with gravitational fields has
attracted a certain attention from the very beginning of wave mechanics fock ;
schrodinger . Due to the smallness of the gravitational interaction, only very
special experiments are possible at the present time (see, e.g., antigravity
), and such a subject is mostly of a pure theoretical interest. But for Dirac
particles, this problem is more fundamental, as they seem to interact directly
only with the gravity (as some other dark matter candidates).

A usual way to describe spinors in gravitational fields (i.e., in a curved
space) is the tetrad formalism (see, e.g., fock ; weinberg ). In a nutshell,
it consists to introduce tetrads (vierbeins) $t_{\mu}^{a}(x)$ instead of
(symmetric) metric tensor $g_{\mu\nu}(x)$, such that

$g_{\mu\nu}(x)=t_{\mu}^{\hskip 3.01389pta}(x)t_{\nu}^{\hskip
3.01389ptb}(x)\eta_{ab}\,,\qquad\eta=\mathrm{diag}(-1,1,1,1)\,.$ (37)

For completeness, the main formulas of the tetrad formalism are presented in
Appendix B. Using the indicated rules, any given relativistic equation in a
flat space (without gravity) can be transformed in a covariant manner to a
curved space (with gravitational fields).

However, the new Dirac equation is not just one equation but an overdetermined
system of equations, and it is not clear which of the many (almost) equivalent
equations in a flat space (cf., (2), (3) and (26)–(29)) have to be converted
to a curved space. Different equations may and will lead to different answers.
For example, it is claimed in ahner  that the use of Equation (31) with
$\kappa=1$ (i.e., for spin-$\tfrac{1}{2}$ particles) leads to a consistent
theory in gravitational fields.

A natural way of introducing gravitational fields into the new Dirac Equations
(2) and (3) consists of transforming these equations themselves to a curved
space. Using formulas from Appendix B, it is easy to check that the
transformed equations can be written in two equivalent forms:

$\big{(}\hat{\gamma}^{\mu}QD_{\mu}-m\,Q\big{)}\Phi=0,\qquad\big{(}\hat{\gamma}^{\mu}D_{\mu}^{\prime}-m\big{)}Q\Phi=0$
(38)

where

$D_{\mu}=\partial_{\mu}+\frac{i}{2}\omega_{\mu ab}S^{ab}\,,\qquad
D_{\mu}^{\prime}=\partial_{\mu}+\frac{1}{2}\omega_{\mu
ab}\big{(}\Sigma^{ab}+iS^{ab}\big{)}\,.$ (39)

As $\Sigma Q=i[Q,S]$ (cf. (20)), these two expressions are consistent. (For
clarity, covariant derivatives acted on different quantities are denoted here
by different symbols.)

As above, one can derive many consequences of that equation. In particular, by
acting on this equation by
$\big{(}\hat{\gamma}^{\nu}D_{\nu}^{\prime}+m\big{)}$, one concludes that

$\Big{(}Q\big{(}g^{\mu\nu}D_{\mu}D_{\nu}-m^{2}\big{)}+\Sigma^{\nu\mu}Q\big{[}D_{\nu},D_{\mu}\big{]}\Big{)}\Phi=0\,.$
(40)

Using (85), it can be rewritten as follows:

$\Big{(}Q\big{(}g^{\mu\nu}D_{\mu}D_{\nu}-m^{2}\big{)}+\frac{i}{2}\Sigma^{\mu\nu}\,Q\,R_{\nu\mu
ab}S^{ab}\Big{)}\Phi=0\,.$ (41)

where $R_{\lambda\rho ab}$ is the Riemann curvature tensor.

This equation has a form similar to (33), and the same arguments can be
applied to cancel the first (scalar) term. Therefore, the necessary conditions
for the existence of solutions require that the following two expressions are
zero (cf. with Equation (36)):

$e^{\mu\nu\lambda\rho}\Gamma^{\nu}R_{\lambda\rho ab}S^{ab}\Phi=0\,,\qquad
e^{\mu\nu\lambda\rho}S^{\mu\nu}R_{\lambda\rho ab}S^{ab}\Phi=0\,.$ (42)

Contrary to the interaction with electromagnetic fields, these equations may
have nontrivial solutions. The point is that $R_{\lambda\rho ab}S^{ab}\Phi$ is
an antisymmetric tensor of the second degree. When the metric is without
torsion, there is only one such tensor, $S_{\lambda\rho}$. It means that

$R_{\lambda\rho ab}S^{ab}\Phi=S_{\lambda\rho}\,\hat{W}\Phi$ (43)

where $\hat{W}$ is a scalar operator. By convoluting the above equation with
$S^{\lambda\rho}$ and using (61), one obtains

$\hat{W}\Phi=-\frac{2}{3}S^{\lambda\rho}R_{\lambda\rho ab}S^{ab}\Phi\,.$ (44)

Substituting (43) into (42) and using identities (62) shows that all Equations
(42) become automatically zero.

Multiplying (41) by $\overline{Q}$, one finds that

$\Big{(}g^{\mu\nu}D_{\mu}D_{\nu}-m^{2}+\frac{1}{2}S^{\mu\nu}R_{\nu\mu
ab}S^{ab}\Big{)}\Phi=0\,.$ (45)

It is plain that the same equation can be obtained directly by the
transformation Equation (27) to the curved space and applying the covariant
derivative to the result.

The corresponding equation for the usual Dirac equation has been derived in
schrodinger , and it has the similar form:

$\Big{(}g^{\mu\nu}D_{\mu}D_{\nu}-m^{2}+\frac{1}{2}\Sigma^{\mu\nu}R_{\nu\mu
ab}\Sigma^{ab}\Big{)}\Psi=0$ (46)

with the difference being that $\Psi$ is a usual four component spinor, and
due to (71), (82), and (83), the term with the Riemann tensor can be reduced
to a scalar curvature $R=g^{\nu\lambda}g^{\mu\rho}R_{\nu\mu\lambda\rho}$

$\frac{1}{2}\Sigma^{\mu\nu}R_{\nu\mu ab}\Sigma^{ab}=-\frac{1}{4}R\,.$ (47)

The above arguments demonstrate how the main obstacle to the interaction of
Dirac particles with electromagnetic fields can be overcome in gravitational
fields, but unfortunately, they do not prove the existence of its solutions.
The author is not aware of mathematical results, which could be applied for
such overdetermined equations, and further investigations are desirable.

##  VII Conclusions

The relativistic invariance implies that elementary particles are transformed
by irreducible representations of the Lorenz group. All known elementary
particles correspond to finite-dimensional representations. But infinite-
dimensional representations of the homogeneous Lorenz group do exist, and it
is natural to conjecture that particles related with such representations
exist as well.

The new Dirac equation describes the simplest of such particles with zero spin
in the rest frame. Originally, this equation has been developed as an example
of a relativistic wave equation whose all solutions have only positive
energies. Similar to all known positive-energy equations, solutions of the new
Dirac equation obey also the infinite-dimensional Majorana equation. But the
latter has too large a spectrum of possible solutions. To select a single
solution, the new Dirac equation is constructed to be not an isolated equation
but an overdetermined system of equations such that one of its consequences is
the Klein–Gordon equation with a fixed mass. This construction necessarily
requires that all these different equations are consistent (i.e., have the
same solutions), which is a nontrivial statement.

A Dirac particle is described by one function depended on four usual space-
time coordinates and two additional variables. The later are only dummy
variables introduced in such a way that the expansion over them generates an
infinite number of particle wave components corresponding to the infinite-
dimensional representation of the Lorenz group.

It was proved that the resulting equations have no common solution when
electromagnetic fields are introduced into the new Dirac equation within the
minimal coupling. It means that Dirac particles cannot, in principle, interact
with electromagnetic fields (with the minimal coupling) but it comes as a
consequence of their internal structure and has not been imposed ad hoc.
Consequently, if Dirac particles exist they will automatically have dark
matter properties.

The distinctive property of Dirac particles is that they have only positive
energy which infers the absence of Dirac antiparticles. Thus, Dirac particles
break the CPT theorem but without violating any fundamental law. It signifies
that such particles might be a new source of observed matter–antimatter
asymmetry.

Dirac particles (as well as other particles transformed by infinite-
dimensional representations of the homogeneous Lorentz group) are unusual,
interesting, and potentially important objects, which may be a clue to a new
physics. Their characteristic properties are peculiar, and to a large extent,
counterintuitive. To summarise:

  * •

They have spin zero (in the rest frame) but are not scalars.

  * •

They obey covariant relativistic equations but have only positive energies,
and consequently do not have antiparticles.

  * •

Their equations are linear in the momenta but do not permit the minimal
coupling to electromagnetic fields.

Unfortunately, Dirac particle investigations are undeveloped. No usual
textbooks have mentioned their existence, no standard models incorporate them,
and many important questions remain unanswered. At the present stage of the
study, it seems that Dirac particles are sterile and can interact directly
only with the gravity. Detailed examinations of the interaction of Dirac
particles with usual matter through gravitational loops and their production
during the Big Bang are essential to give them the status of established
elementary particles capable of explaining dark matter.

Such inquiries are not straightforward and require the development of new
methods, but the possibility that Dirac particles could dominate the universe
clearly stresses the importance of the topic and the necessity to investigate
it seriously.

###### Acknowledgements.

The author is grateful to M. Cirelli for useful comments.

##  Appendix A Useful Relations

  * •

The choice in Dirac_I  is equivalent to the following gamma matrices obeying
(8): $\gamma^{0}=\beta$,

$\gamma^{1}=\Big{(}\begin{array}[]{cc}{-\sigma_{3}}&{0}\\\
{0}&{\sigma_{3}}\end{array}\Big{)}\,,\quad\gamma^{2}=\Big{(}\begin{array}[]{cc}{\sigma_{1}}&{0}\\\
{0}&{-\sigma_{1}}\end{array}\Big{)}\,,\quad\gamma^{3}=\Big{(}\begin{array}[]{cc}{0}&{-1}\\\
{-1}&{0}\end{array}\Big{)}\,,$ (48)

and

$\gamma^{5}=\gamma^{0}\gamma^{1}\gamma^{2}\gamma^{3}=\Big{(}\begin{array}[]{cc}{i\sigma_{2}}&{0}\\\
{0}&{-i\sigma_{2}}\end{array}\Big{)}$ (49)

where $\sigma_{i}$ are the usual $2\times 2$ Pauli matrices

$\sigma_{1}=\Big{(}\begin{array}[]{cc}{0}&{1}\\\
{1}&{0}\end{array}\Big{)}\,,\quad\sigma_{2}=\Big{(}\begin{array}[]{cc}{0}&{-i}\\\
{i}&{0}\end{array}\Big{)}\,,\quad\sigma_{3}=\Big{(}\begin{array}[]{cc}{1}&{0}\\\
{0}&{-1}\end{array}\Big{)}\,.$ (50)

  * •

When working with quadratic operators in $Q$, the following two formulas are
useful.

    * (a)

For an antisymmetric matrix $G_{jk}=-G_{kj}$ there is an operator identity

$Q^{T}GQ\equiv q_{j}G_{jk}q_{k}=-\frac{i}{2}\mathrm{Tr}(G\beta)\,.$ (51)

    * (b)

Let us consider the following transformation of a $4\times 4$ matrix $A$ (see,
e.g., lectures ):

$\mathcal{D}:A\to\mathcal{D}(A)=\frac{i}{2}\overline{Q}AQ\,,\qquad\overline{Q}=Q^{T}\beta\,.$
(52)

If $\beta A$ and $\beta B$ are symmetric matrices, then the direct
calculations prove that this transformation preserves commutator relations

$[\mathcal{D}(A),\mathcal{D}(B)]=\mathcal{D}([A,B])\ .$ (53)

  * •

The spin operators, which realise the infinite-dimensional representation of
the Lorentz group, are defined in (21):
$S^{cd}=\frac{1}{2}\overline{Q}\Sigma^{cd}Q$. Rewriting them by components
gives

$\displaystyle S^{01}$ $\displaystyle=$
$\displaystyle\frac{1}{4}(q_{1}^{2}-q_{3}^{2}-q_{2}^{2}+q_{4}^{2})\,,\qquad
S^{02}=\frac{1}{2}(q_{3}q_{4}-q_{1}q_{2})\,,$ $\displaystyle S^{03}$
$\displaystyle=$ $\displaystyle\frac{1}{2}(q_{1}q_{3}+q_{4}q_{2})\,,\qquad
S^{12}=\frac{1}{2}(q_{2}q_{3}-q_{1}q_{4})\,,$ (54) $\displaystyle S^{13}$
$\displaystyle=$
$\displaystyle\frac{1}{4}(q_{2}^{2}+q_{4}^{2}-q_{1}^{2}-q_{3}^{2}),\qquad
S^{23}=\frac{1}{2}(q_{1}q_{2}+q_{3}q_{4})\,.$

  * •

The components of vector operator
$\Gamma^{c}=-\frac{1}{4}\overline{Q}\gamma^{c}Q$ are

$\displaystyle\Gamma^{0}$ $\displaystyle=$
$\displaystyle\frac{1}{4}(q_{1}^{2}+q_{2}^{2}+q_{3}^{2}+q_{4}^{2})\,,\quad\Gamma^{1}=\frac{1}{2}(q_{2}q_{4}-q_{1}q_{3})\,,$
(55) $\displaystyle\Gamma^{2}$ $\displaystyle=$
$\displaystyle\frac{1}{2}(q_{1}q_{4}+q_{2}q_{3})\,,\qquad\Gamma^{3}=\frac{1}{4}(q_{1}^{2}+q_{2}^{2}-q_{3}^{2}-q_{4}^{2})\,.$

  * •

In the calculations, the following identities are also helpful (see, e.g.,
lectures ):

$\gamma^{5}\Big{(}\gamma^{\mu}\gamma^{\nu}-\gamma^{\nu}\gamma^{\mu}\Big{)}=e^{\mu\nu}_{\hskip
6.02777pt\rho\sigma}\gamma^{\rho}\gamma^{\sigma}$ (56)

and

$\frac{1}{2}\gamma^{\lambda}\Big{(}\gamma^{\mu}\gamma^{\nu}-\gamma^{\nu}\gamma^{\mu}\Big{)}=\eta^{\lambda\mu}\gamma^{\nu}-\eta^{\lambda\nu}\gamma^{\mu}+e^{\lambda\mu\nu}_{\hskip
9.04166pt\rho}\gamma^{5}\gamma^{\rho}\,.$ (57)

Here, $e^{\mu\nu\rho\sigma}$ is the completely antisymmetric tensor with
$e^{0123}=1$.

  * •

For convenience, the commutators of spin operators $S$ and $\Gamma$ are listed
below:

$\displaystyle\big{[}S^{cd},\,S^{kr}\big{]}$ $\displaystyle=$ $\displaystyle
i\big{(}\eta^{ck}S^{dr}+\eta^{dr}S^{ck}-\eta^{dk}S^{cr}-\eta^{cr}S^{dk}\big{)}\,,$
(58) $\displaystyle\big{[}\Gamma^{c},\,\Gamma^{d}\big{]}$ $\displaystyle=$
$\displaystyle-iS^{cd}\,,$ (59)
$\displaystyle\big{[}S^{cd},\,\Gamma^{b}\big{]}$ $\displaystyle=$
$\displaystyle i\big{(}\eta^{cb}\Gamma^{d}-\eta^{db}\Gamma^{c}\big{)}\,.$ (60)

These commutators are the periphrases of the corresponding relations for the
usual $\gamma$-matrices.

  * •

The above operators obey a number of identities (see e.g., staunton ):

$\displaystyle
S_{ab}\,S^{ab}=-\frac{3}{2},\qquad\Gamma_{a}\,\Gamma^{a}=\frac{1}{2}\,,$ (61)
$\displaystyle e^{abcd}S_{ab}\,S_{cd}=0,\qquad
e^{abcd}S_{ab}\,\Gamma_{c}=0\,,$ (62)
$\displaystyle\Gamma_{a}\,S^{ab}=-\frac{3i}{2}\Gamma^{b},\qquad
S^{ab}\,\Gamma_{b}=-\frac{3i}{2}\Gamma^{a},$ (63) $\displaystyle S^{a}_{\hskip
3.01389ptb}\,S^{bc}=\frac{1}{2}\big{(}\eta^{ac}-2\Gamma^{a}\,\Gamma^{c}-3iS^{ac}\big{)}\,.$
(64)

  * •

Let us introduce a $5\times 5$ metric tensor
$\eta^{AB}=\mathrm{diag}(-1,1,1,1-1)$ with $A,B=0,\ldots 4$ and denote
$S^{4\,a}=\Gamma^{a}$. Ten antisymmetric operators $S^{ab}$ and $S^{4\,a}$
with $a,b=0,\ldots 3$ form a representation of the group O($3,2$) as the above
commutation relations can be rewritten in the form (cf. with Equation (18))

$\big{[}S^{CD},\,S^{KR}\big{]}=i\big{(}\eta^{CK}S^{DR}+\eta^{DR}S^{CK}-\eta^{DK}S^{CR}-\eta^{CR}S^{DK}\big{)}$
(65)

where majuscule letters are from $0$ to $4$.

  * •

The Hermite functions used in the expansion (23) have the form (see, e.g.,
bateman )

$\phi_{n}(q)=\Big{(}2^{n}n!\sqrt{\pi}\Big{)}^{-1/2}H_{n}(q)e^{-q^{2}/2}$ (66)

where $H_{n}(q)$ are the Hermite polynomials. These functions are orthogonal:

$\int_{-\infty}^{\infty}\phi_{m}(q)\phi_{n}(q)dq=\delta_{mn}$ (67)

and obey the following relations:

$\displaystyle q\phi_{n}(q)$ $\displaystyle=$
$\displaystyle\sqrt{\frac{n}{2}}\phi_{n-1}(q)+\sqrt{\frac{n+1}{2}}\phi_{n+1}(q)\,,$
(68) $\displaystyle\phi_{n}^{\prime}(q)$ $\displaystyle=$
$\displaystyle\sqrt{\frac{n}{2}}\phi_{n-1}(q)-\sqrt{\frac{n+1}{2}}\phi_{n+1}(q)\,.$

The use of these identities permits one to rewrite any equation in auxiliary
variables $q_{j}$ in the form of infinite-dimensional matrices without these
variables.

##  Appendix B Tetrad Formalism

A tetrad (vierbein) is, by definition, a set of (usually orthogonal) vectors
$t^{\mu}_{a}(x)$ defined in all space points $x$ which transform the true
symmetric metric tensor of a curved space $g^{\mu\nu}(x)$ locally into the
flat Minkowski metric tensor $\eta=\mathrm{diag}(-1,1,1,1)$ (see, e.g.,
weinberg )

$\eta_{ab}=t^{\mu}_{a}(x)t^{\nu}_{b}(x)g_{\mu\nu}(x)\,.$ (69)

Here Latin, (resp., Greek) letters indicate flat (resp., curved) coordinates.

The curved-space gamma matrices are determined as follows:

$\tilde{\gamma}^{\mu}(x)=t^{\mu}_{a}(x)\gamma^{a}\,.$ (70)

These $4\times 4$ matrices obey the relation

$\tilde{\gamma}^{\mu}(x)\tilde{\gamma}^{\nu}(x)+\tilde{\gamma}^{\nu}(x)\tilde{\gamma}^{\mu}(x)=2g^{\mu\nu}(x)\,.$
(71)

For clarity, all quantities multiplied by tetrads will be indicated by the
same letter as in the flat space but with a tilde:

$\tilde{\Gamma}_{\mu}(x)=t^{\mu}_{a}(x)\Gamma^{a}\,,\qquad\tilde{\Sigma}^{\mu\nu}(x)=t^{\mu}_{a}(x)t^{\nu}_{b}(x)\Sigma^{ab}\,,\qquad\tilde{S}^{\mu\nu}(x)=t^{\mu}_{a}(x)t^{\nu}_{b}(x)S^{ab}\,.$
(72)

At a fixed point, a (symmetric) metric tensor is determined by ten components,
but there exist sixteen tetrads. The remaining six free parameters correspond
to six-parameter local Lorentz transformations, which preserve the Minkowski
metric tensor.

It is well known that to conserve a local invariance, it is necessary to
change the usual derivatives to covariant ones:

$\partial_{\mu}\Psi(x)\longrightarrow
D_{\mu}\Psi(x)=\big{(}\partial_{\mu}+C_{\mu})\Psi(x)$ (73)

where $C_{\mu}$ is a connection whose explicit form depends on a quantity
$\Psi(x)$ on which the derivative acts. In general, if under the infinitesimal
local Lorentz transformations $x^{\mu}=x^{\mu}+a^{\mu}_{\nu}(x)x^{\nu}$ (cf.
(16)) a field $\Psi(x)$ transforms in the following way
$\Psi(x)\to\Psi(x)+a_{ab}(x)\hat{s}^{ab}\Psi(x)$, then

$C_{\mu}=\omega_{\mu ab}\,\hat{s}^{ab}$ (74)

with certain coefficients $\omega_{\mu ab}$ depending on quantities
considered. It is common to denote a covariant derivative by $D_{\mu}$, but a
care should be taken when its explicit expression is used in calculations. For
vector (and tensor) fields, $\omega_{\mu ab}=\Gamma_{\mu ab}$, and it is
called the Christoffel symbol. Different cases are indicated in Table 1.

Table 1: Covariant derivatives of different quantities. Contravariant Lorentz vector | $V^{\mu}\to V^{\mu}+a^{\mu}_{\lambda}(x)V^{\lambda}$ | $D_{\nu}V^{\mu}=\partial_{\nu}V^{\mu}+\Gamma^{\mu}_{\nu\lambda}V^{\lambda}$
---|---|---
Contravariant tetrad vector | $V^{b}\to V^{b}+a^{b}_{c}(x)V^{c}$ | $D_{\nu}V^{b}=\partial_{\nu}V^{b}+\omega_{\nu c}^{\hskip 6.02777ptb}V^{c}$
Usual Lorentz spinor | $\Psi\to\Psi+\frac{1}{2}a_{ab}(x)\Sigma^{ab}\Psi$ | $D_{\nu}\Psi=\big{(}\partial_{\nu}+\frac{1}{2}\omega_{\nu ab}\Sigma^{ab}\big{)}\Psi$
Dirac particle | $\Phi\to\Phi+\frac{i}{2}a_{ab}(x)S^{ab}\Phi$ | $D_{\nu}\Phi=\big{(}\partial_{\nu}+\frac{i}{2}\omega_{\nu ab}S^{ab}\big{)}\Phi$

For covariant quantities, the connections change the sign. When tensors are
considered, it is necessary to add the corresponding connections for each
index. For matrices (and operators) such transformation corresponds to the
commutator with the connections.

Covariant derivatives transform in the same way as the corresponding
quantities in the flat space. In particular, they have to commute with the
rising and lowering indices with the metric tensor. Therefore,

$D_{\nu}g_{\mu\lambda}=0\,,\qquad\partial_{\nu}g_{\mu\lambda}-\Gamma_{\nu\mu}^{\beta}g_{\beta\lambda}-\Gamma_{\nu\lambda}^{\beta}g_{\mu\beta}=0\,.$
(75)

This relation implies that

$\Gamma^{\mu}_{\nu\lambda}=\frac{1}{2}g^{\mu\rho}\left(\partial_{\nu}g_{\rho\lambda}+\partial_{\lambda}g_{\rho\nu}-\partial_{\rho}g_{\mu\nu}\right),\qquad\Gamma^{\mu}_{\nu\lambda}=\Gamma^{\mu}_{\lambda\nu}.$
(76)

Similarly, it is convenient to require that

$D_{\nu}t_{\mu}^{\hskip 3.01389pta}=0\,,\qquad\partial_{\nu}t_{\mu}^{\hskip
3.01389pta}+\Gamma_{\nu\mu}^{\lambda}+\omega_{\nu c}t_{\mu}^{\hskip
3.01389ptc}=0$ (77)

which leads to

$\omega_{\nu c}^{\hskip
6.02777ptb}=t^{\mu}_{c}\big{(}-\partial_{\nu}t_{\mu}^{b}+\Gamma^{\lambda}_{\nu\mu}t^{b}_{\lambda}\big{)}\,,\qquad\omega_{\mu
cb}=-\omega_{\mu bc}\,.$ (78)

The same answer can be obtained from

$D_{\nu}\hat{\gamma}^{\mu}(x)=0\,,\qquad\partial_{\nu}\hat{\gamma}^{\lambda}(x)+\Gamma_{\nu\lambda}^{\mu}\hat{\gamma}^{\lambda}(x)+\big{[}g_{\nu},\hat{\gamma}^{\mu}(x)\big{]}=0\,,\qquad
g_{\mu}=\frac{1}{2}\omega_{\mu ab}\Sigma^{ab}\,.$ (79)

In the same way,

$D_{\nu}\hat{\Gamma}^{\mu}(x)=0\,,\qquad\partial_{\nu}\hat{\Gamma}^{\lambda}(x)+\Gamma_{\nu\lambda}^{\mu}\hat{\Gamma}^{\lambda}(x)+\big{[}G_{\nu},\hat{\Gamma}^{\mu}(x)\big{]}=0\,,\qquad
G_{\mu}=\frac{i}{2}\omega_{\mu ab}S^{ab}\,.$ (80)

By definition, the commutator of two covariant derivatives acting at a
covariant vector is determined by the Riemann curvature tensor
$R_{\nu\mu\rho}^{\hskip 9.04166pt\lambda}$

$\big{[}D_{\mu},\,D_{\nu}\big{]}A_{\rho}=R_{\nu\mu\rho}^{\hskip
9.04166pt\lambda}A_{\lambda}\,,\qquad R_{\nu\mu\rho}^{\hskip
9.04166pt\lambda}=\partial_{\nu}\Gamma_{\mu\rho}^{\lambda}-\partial_{\mu}\Gamma_{\nu\rho}^{\lambda}+\Gamma_{\mu\rho}^{\sigma}\Gamma_{\nu\sigma}^{\lambda}-\Gamma_{\nu\rho}^{\sigma}\Gamma_{\mu\sigma}^{\lambda}\,.$
(81)

It is known that the Riemann tensor obeys the symmetry relations (see, e.g.,
weinberg )

$R_{\nu\mu\lambda\rho}=-R_{\mu\nu\lambda\rho}=-R_{\nu\mu\rho\lambda}=R_{\lambda\rho\nu\mu}$
(82)

and the Bianchi identity

$R_{\nu\mu\lambda\rho}+R_{\nu\rho\mu\lambda}+R_{\nu\lambda\rho\mu}=0\,.$ (83)

As $D_{\nu}t_{\rho}^{\hskip 3.01389pta}=0$ it follows that
$\big{[}D_{\mu},\,D_{\nu}\big{]}t_{\rho}^{\hskip 3.01389pta}=0$, from which
one proves the validity of the following identity:

$\partial_{\mu}\omega_{\nu b}^{a}-\partial_{\nu}\omega_{\mu b}^{a}+\omega_{\mu
c}^{a}\omega_{\nu b}^{c}-\omega_{\nu c}^{a}\omega_{\mu
b}^{c}=-R_{\nu\mu\rho}^{\hskip 9.04166pt\lambda}\,t_{\lambda}^{\hskip
3.01389pta}\,t^{\rho}_{\hskip 3.01389ptb}\,.$ (84)

Using this expression, one obtains

$\big{[}D_{\mu},\,D_{\nu}\big{]}\Phi=\frac{i}{2}R_{\nu\mu ab}S^{ab}\Phi\,.$
(85)

## References

  * (1) Jungman, G.; Kamionkowski, M.; Griest, K. Supersymmetric dark matter. _Phys. Rep._ 1996, _267_ , 195.
  * (2) Bertone, G.; Hooper, D.; Silk, J. Particle dark matter: Evidence, candidates and constraints. _Phys. Rep._ 2005, _405_ , 279.
  * (3) Bertone, G.; Tait, T.M.P. A new era in the quest for dark matter. _Nature_ 2018, _562_ , 51.
  * (4) Dirac, P.A.M. A positive energy relativistic wave equation. _Proc. R. Soc. Lond. A_ 1971, _322_ , 435.
  * (5) Dirac, P.A.M. The quantum theory of the electron. _Proc. R. Soc. Lond. A_ 1928, _117_ , 610.
  * (6) Cirilo-Lombardo, D.J.; Sanchez, N.G. Quantum-spacetime symmetries: A principle of minimum group representation. _Universe_ 2024, _10_ , 22.
  * (7) Majorana, E. Teoria relativistica di particelle con momento intriciseco arbitrario. _Nuovo Cimento_ 1932, _9_ , 335.
  * (8) Fradkin, D.M. Comments on a paper by Majorana concerning elementary particles. _EJTP_ 2006, _3_ , 305.
  * (9) Majorana, E. Teoria simmetrica dell’ elettrone e del positrone. _Nuovo Cimento_ 1937, _14_ , 171.
  * (10) Harish-Chandra, H. Infinite irreducible representations of the Lorentz group. _Proc. R. Soc. A_ 1947, _189_ , 372.
  * (11) Gelfand, I.M.; Naimark, M.A. Unitary representations of the Lorentz group. _Izvestiya Akad. Nauk SSSR Ser. Mat._ 1947, _11_ , 411.
  * (12) Athron, P.; Kozar, N.A.; Balázs, C.; Beniwal, A.; Bloor, S.; Bringmann, T.; Brod, J.; Chang, C.; Cornell, J.M.; Farmer, B.; et al. Thermal WIMPs and the scale of new physics: Global fits of Dirac dark matter effective field theories. _Eur. Phys. J. C_ 2021, _81_ , 992.
  * (13) Lawson, H.B. Foliations. _Bull. Amer. Math, Soc._ 1974, _80_ , 369.
  * (14) Dirac, P.A.M. A positive-energy relativistic wave equation. II. _Proc. R. Soc. Lond. A_ 1972, _328_ , 1.
  * (15) Dirac, P.A.M. Unitary representations of the Lorentz group. _Proc. R. Soc. Lond. A_ 1945, _183_ , 284.
  * (16) Dirac, P.A.M. A remarkable representation of the 3 + 2 de Sitter group. _J. Math. Phys._ 1963, _4_ , 901.
  * (17) Bierderman, L.C.; van Dam, H.; Mukunda, N.; Böhm, A.; Dolland, J.D. Dirac’s new relativistic equation and its generalization, Chapter II. In _Relativistic Models of Extended Hadrons Obeying a Mass-Spin Trajectory Constraint_ ; Lecture Notes in Physics; Springer: Berlin/Heidelberg, Germany, 1982; Volume 165.
  * (18) Staunton, L.P. A spin-$\frac{1}{2}$ positive-energy relativistic equation. _Phys. Rev. D_ 1974, _10_ , 1760.
  * (19) Biedenharn, L.C.; Han, M.Y.; van Dam, H. Generalization and interpretation of Dirac’s positive-energy relativistic wave equation. _Phys. Rev. D_ 1973, _8_ , 1735.
  * (20) Mukunda, N. The Dirac equation—Old and new. _Curr. Sci._ 1982, _51_ , 78.
  * (21) Mukunda, N.; Chiang, C.C.; Sudarshan, E.C.G. Dirac’s new relativistic wave equation in interaction with an electromagnetic field. _Proc. R. Soc. Lond. A_ 1982, _379_ , 103.
  * (22) Sudarshan, E.C.G.; Mukunda, N.; Chiang, C.C. Dirac positive-energy wave equation with para-Bose internal variables. _Phys. Rev. D_ 1982, _25_ , 3237.
  * (23) Fock, V. Geometrization of the Dirac theory of electrons. _Zeit. Phys._ 1929, _57_ , 261.
  * (24) Schrödinger, E. Diracsches Elektron im Schwerefeld I. Sitzungsberichte der preussischen Akademie der Wissenschaften, physikalisch-mathematische Klasse. In _Quantum Gravity in the First Half of the Twentieth Century: A Sourcebook_ ; Blum, A.S., Rickles, D., Eds.; Max-Planck-Gesellschaft zur Förderung der Wissenschaften: Berlin, Germany, 2018. Available online: https://edition-open-sources.org/sources/10/
  * (25) Anderson, E.K.; Baker, C.J.; Bertsche, W.; Bhatt, N.M.; Bonomi, G.; Capra, A.; Carli, I.; Cesar, C.L.; Charlton, M.; Christensen, A.; et al. Observation of the effect of gravity on the motion of antimatter. _Nature_ 2023, _621_ , 716.
  * (26) Weinberg, S. _Gravitation and Cosmology. Principles and Applications of the Theory of General Relativity_ ; Wiley: New York, NY, USA, 1972.
  * (27) Ahner, H.F. Gravitation and positive-energy equation. _Phys. Rev. D_ 1975, 11, 3384.
  * (28) Erdelyi, A. _Higher Transcendental Functions_ ; McGraw-Hill : New York, 1953; Volume 2.


# An accurate and efficient Ehrenfest dynamics approach for calculating linear
and nonlinear electronic spectra

Austin O. Atsango  Department of Chemistry, Stanford University, Stanford,
California, 94305, USA  Andrés Montoya-Castillo
<EMAIL_ADDRESS>Department of Chemistry, University of
Colorado Boulder, Boulder, Colorado, 80309, USA  Thomas E. Markland
<EMAIL_ADDRESS>Department of Chemistry, Stanford University, Stanford,
California, 94305, USA

###### Abstract

Linear and nonlinear electronic spectra provide an important tool to probe the
absorption and transfer of electronic energy. Here we introduce a pure state
Ehrenfest approach to obtain accurate linear and nonlinear spectra that is
applicable to systems with large numbers of excited states and complex
chemical environments. We achieve this by representing the initial conditions
as sums of pure states and unfolding multi-time correlation functions into the
Schrödinger picture. By doing this we show that one can obtain significant
improvements in accuracy over the previously used projected Ehrenfest approach
and that these benefits are particularly pronounced in cases where the initial
condition is a coherence between excited states. While such initial conditions
do not arise when calculating linear electronic spectra, they play a vital
role in capturing multidimensional spectroscopies. We demonstrate the
performance of our method by showing that it is able to quantitatively capture
the exact linear, 2DES, and pump-probe spectra for a Frenkel exciton model in
slow bath regimes and is even able to reproduce the main spectral features in
fast bath regimes.

##  I Introduction

Nonlinear optical spectroscopy is a powerful tool for probing the structure
and dynamics of chemical systems Ginsberg, Cheng, and Fleming (2009); Cho
(2008); Khalil, Demirdöven, and Tokmakoff (2003), and it typically provides
significantly more information than its linear counterparts. For example, 2D
electronic spectroscopy (2DES) Hybl  _et al._ (1998); Mukamel (1995) is
frequently used to elucidate inter-chromophoric couplings, energy transfer and
relaxation processes, and environmental effects in systems such as
photosynthetic molecular aggregates with femtosecond time resolution Brixner
_et al._ (2005); Engel  _et al._ (2007); Dostál  _et al._ (2012); Abramavicius
and Mukamel (2010). However, interpreting these spectra in terms of the
specific molecular structures and motions that give rise to the electronic
states and their relaxation pathways often requires the assistance of
simulations to disentangle spectral features such as short-time coherent
oscillations and long-time population decay pathways. Accurately and
efficiently simulating 2DES signals from atomistic simulations and uncovering
how they arise from the underlying quantum dynamics of the nuclear and
electronic states remains a significant challenge.

Simulating nonlinear optical spectra requires obtaining higher-order response
functions and hence the number of methods that have successfully been applied
to generate them is significantly more limited than for the less
computationally demanding linear spectra. Exact methods such as the
Hierarchical Equations of Motion (HEOM) Tanimura and Kubo (1989); Ishizaki and
Tanimura (2007); Kreisbeck  _et al._ (2011); Chen  _et al._ (2010) and
multiconfiguration time-dependent Hartree (MCTDH) Meyer (2011); Schulze  _et
al._ (2017) provide important insights and benchmarks for approximate methods.
However, the desire to treat large condensed-phase systems containing many
electronic states with a fully atomistic treatment of the nuclear motions and
even on-the-fly evaluation of the electronic surfaces has spurred the recent
development of approximate dynamics methods. When approximate methods are
used, 2DES has been shown to provide a much stricter test of the quantum
dynamics method employed than linear electronic spectroscopy. For example,
when using Redfield theoryREDFIELD (1965); Bloch (1957) it has been shown that
in parameter regimes where the linear electronic spectrum can be
quantitatively captured, the 2DES spectrum significantly deviates from the
exact results Fetherolf and Berkelbach (2017), although this can be somewhat
alleviated by freezing the low-frequency bath modes.Montoya-Castillo,
Berkelbach, and Reichman (2015); Fetherolf and Berkelbach (2017)

Trajectory-based methods provide a particularly appealing approach to simulate
2DES since they are typically compatible with a fully atomistic treatment of
the nuclear motions and even on-the-fly evaluation of the electronic surfaces
on which the nuclei evolve. Many recent applications of trajectory-based
methods to 2DES have been based on semiclassical treatments of the Meyer-
Miller-Stock-Thoss Meyer and Miller (1979); Stock and Thoss (1997) mapping of
the electronic degrees of freedom. In particular, the partially linearized
density matrix (PLDM) approach Huo and Coker (2011) and its more recent spin-
PLDM variant Mannouch and Richardson (2020a, b) have been shown to produce
accurate 2DES spectra in parameter regimes where perturbative methods break
down Provazza  _et al._ (2018); Mannouch and Richardson (2022). Further work
has shown how mapping-based semiclassical methods can be used to simulate 2DES
beyond the perturbative treatment of the field-matter interaction Gao and Geva
(2020). The optimized mean trajectory approach has also been recently coupled
with trajectories based on the Meyer-Miller Hamiltonian to compute two-
dimensional vibrational-electronic spectra Polley and Loring (2021); Loring
(2022). Other trajectory-based approaches include the numerical integration of
the Schrodinger equation (NISE) Jansen and Knoester (2006); Liang and Jansen
(2012); Torii (2006) and the stochastic Liouville Equation methods Jansen,
Zhuang, and Mukamel (2004); Jansen  _et al._ (2005), which involve an explicit
treatment of the bath degrees of freedom but do not account for the back-
reaction of the electronic degrees of freedom on the bath, leading to
inaccuracies in linear and nonlinear spectra. To correct this, NISE has been
combined with the fewest switches surface hopping approach, allowing it to
incorporate the quantum back-reaction on the bath and thus improve its
accuracy.Tempelaar  _et al._ (2013)

Here, we present an accurate method to simulate 2DES using Ehrenfest dynamics
and contrast it with a previously reported method detailed in Ref. van der
Vegte  _et al._ , 2013. By applying our method, we show that significantly
more accurate 2DES spectra can be obtained from Ehrenfest dynamics if one uses
an appropriate initialization of coherence and mixed states. In particular, we
express all initial conditions as sums of pure states, which have important
implications for the accuracy of the resulting dynamics and allow the
simulation to be done unambiguously in either the wavefunction or density
matrix formulations of Ehrenfest theory. We show that when formulated in this
way, Ehrenfest theory can closely reproduce the HEOM result, and that this
result is not achieved when initialization is not done from pure states.

##  II Pure State and Projected Ehrenfest Methods

It has previously been shown that in the density matrix formulation of the
Ehrenfest method Kapral and Ciccotti (1999), a carefully constructed
initialization is required to obtain unambiguous results because density
matrices initialized from non-pure states do not yield unique dynamics
Montoya-Castillo and Reichman (2016). Conversely, the wavefunction formulation
of Ehrenfest does not suffer from this problem because all initial conditions
are automatically pure states. Our approach to the density matrix formulation
of Ehrenfest avoids this previously reported ambiguity by expanding any given
initial density matrix $\rho_{0}$ as a sum of pure states,

$\rho_{0}=\sum_{\rm{i}}a_{\rm{i}}\ket{\psi_{\rm{i}}}\bra{\psi_{\rm{i}}}.$ (1)

For an arbitrary initial density matrix $\rho_{0}=\ket{a}\bra{b}$, the pure
state decomposition can be accomplished as follows:

$\ket{a}\bra{b}=\frac{1}{2}\left(\rho^{+}_{\rm{ab}}+i\rho^{-}_{\rm{ab}}-\left[1+i\right]\left(\rho_{\rm{aa}}+\rho_{\rm{bb}}\right)\right)$
(2)

where

$\begin{split}\rho_{\rm{aa}}&=\ket{a}\bra{a}\\\
\rho_{\rm{bb}}&=\ket{b}\bra{b}\\\
\rho^{+}_{\rm{ab}}&=\frac{1}{2}\left(\ket{a}+\ket{b}\right)\left(\bra{a}+\bra{b}\right)\\\
\rho^{-}_{\rm{ab}}&=\frac{1}{2}\left(\ket{a}+i\ket{b}\right)\left(\bra{a}-i\bra{b}\right)\end{split}$
(3)

are pure states. Here we will refer to this approach as pure state Ehrenfest
and contrast it with the unmodified alternative, which we refer to as
projected Ehrefest, where all initial conditions $\ket{a}\bra{b}$ are
propagated directly without having been decomposed into pure states.

##  III Simulation Details: Frenkel exciton model

We consider a Frenkel exciton model of coupled chromophores where the full
matter Hamiltonian is

$\hat{H}_{\rm{mat}}=\hat{H}_{\rm{s}}+\hat{H}_{\rm{b}}+\hat{H}_{\rm{sb}}.$ (4)

Here, $\hat{H}_{\rm{s}}$ is the system Hamiltonian, $\hat{H}_{\rm{b}}$ is the
bath Hamiltonian, and $\hat{H}_{\rm{sb}}$ is the system-bath coupling. The
system Hamiltonian, $\hat{H}_{\rm{s}}$, consists of individual chromophores,
each with a site energy $\epsilon_{\rm{m}}$, that couple electronically via
the transfer integrals $J_{\rm{mn}}$:

$\hat{H}_{\rm{s}}=\sum_{\rm{m}}\epsilon_{\rm{m}}\ket{m}\bra{m}+\sum_{\rm{m\neq
n}}J_{\rm{mn}}\ket{m}\bra{n}.$ (5)

The bath Hamiltonian, $\hat{H}_{\rm{b}}$, consists of independent sets of
phonon modes coupled to each chromophore. It is expressed in terms of the
phonon mode frequencies $\omega$ and their mass-weighted momenta and
coordinates $\hat{P}$ and $\hat{Q}$ as:

$\hat{H}_{b}=\sum_{\rm{m}}\sum_{\rm{j=1}}^{\rm{N^{m}_{b}}}\left(\frac{1}{2}\hat{P}^{2}_{\rm{mj}}+\frac{1}{2}\omega^{2}_{\rm{mj}}\hat{Q}^{2}_{\rm{mj}}\right).$
(6)

Each chromophore site $\ket{m}\bra{m}$ is linearly coupled to its local bath
of phonon modes, such that $\hat{H}_{\rm{sb}}$ takes the form

$\hat{H}_{sb}=\sum_{\rm{m}}\sum_{\rm{j=1}}^{\rm{N^{m}_{b}}}c_{\rm{mj}}\hat{Q}_{\rm{mj}}\ket{m}\bra{m},$
(7)

where $c_{\rm{mj}}$ is the coupling strength of the $j^{th}$ phonon mode
attached to chromophore $m$. The characteristics of the baths and their effect
on the chromophores are specified via the spectral density,

$J_{\rm{m}}(\omega)=\frac{\pi}{2}\sum_{\rm{j}}\frac{c^{2}_{\rm{mj}}}{\omega_{\rm{j}}}\delta(\omega-\omega_{\rm{j}}).$
(8)

All chromophore sites are assumed to have identical spectral densities
($J_{\rm{m}}(\omega)=J(\omega)$). Here, we use an Ohmic spectral density with
a Lorentzian cutoff (Debye),

$J(\omega)=\frac{2\lambda\omega_{\rm{c}}\omega}{\omega_{\rm{c}}^{2}+\omega^{2}},$
(9)

where $\omega_{\rm{c}}$ is the bath cut-off frequency and
$\lambda=(\hbar\pi)^{-1}\int_{0}^{\infty}d\omega J(\omega)/\omega$ is the
reorganization energy.

To obtain linear and non-linear spectra, we treat the light-matter interaction
perturbatively Mukamel (1995), i.e.,

$\hat{H}_{\rm{spec}}=\hat{H}_{\rm{mat}}-\bm{\mu}\cdot\bm{E(t)},$ (10)

where $\bm{E(t)}$ is the classical electric field and $\bm{\mu}$ is the total
dipole operator,

$\bm{\mu}=\sum_{\rm{m}}\mu_{\rm{m}}(\ket{m}\bra{0}+\ket{0}\bra{m}).$ (11)

We consider a system of two coupled chromophores with $\epsilon_{1}=50$ cm-1,
$\epsilon_{2}=-50$ cm-1, and $J_{12}=$ 100 cm-1, consistent with that used in
previous studies Ishizaki and Fleming (2009); Fetherolf and Berkelbach (2017).
In this system, there are four states, all of which are accessible: the ground
state, two singly excited states, and one doubly excited state. We report
results for two sets of bath parameters: a slow bath with a relaxation time
$\omega_{c}^{-1}=300$ fs and a fast bath with a relaxation time
$\omega_{c}^{-1}$ = 17.7 fs, both at $k_{B}T=208$ cm-1. Each bath was
discretized using 300 phonon modes via the method employed in Ref. Berkelbach,
Reichman, and Markland, 2012 that is designed to yield the exact
reorganization energy regardless of the number of modes used. The
reorganization energy $\lambda$ was set to 50 cm-1. All spectra were
calculated with a transition dipole matrix where $\mu_{1}/\mu_{2}=-5$. The
slow bath parameter regime was previously investigated via HEOM and Redfield
theory and its frozen mode variant in Ref. Fetherolf and Berkelbach, 2017.

Our HEOM and projected Ehrenfest calculations were conducted using the python
package pyrho Berkelbach  _et al._ (2020). For the HEOM calculations,
converged results for the regimes studied here were obtained by employing the
high-temperature approximation, i.e. using zero Matsubara frequencies ($K=0$),
and truncating the auxiliary density matrices at $L=15$. Both the HEOM and
projected Ehrenfest equations of motion were propagated using a fourth-order
Runge-Kutta integrator, while the pure state Ehrenfest method was propagated
via the split Liouvillian method, allowing for the use of larger time steps.
In the slow bath parameter regime, we used time steps of 10 fs for HEOM and 2
fs for both versions of Ehrenfest, while in the fast bath parameter regime, we
used a time step of 2.5 fs for HEOM and 10 fs for pure state Ehrenfest. The
linear spectra were generated using 20000 trajectories of length 400 fs, while
the nonlinear spectra used 20000 trajectories of length 300 fs for the slow
bath regime and 200 fs for the fast bath regime.

##  IV Results

Here we compare the Ehrenfest results for both the linear and nonlinear
optical spectra of the Frenkel exciton model outlined in Sec. III to spectra
obtained using the numerically exact HEOM approach. We demonstrate that, while
in linear absorption spectra the differences between pure state and projected
Ehrenfest are subtle, they become much more pronounced for nonlinear spectra.

###  IV.1 Linear Optical Spectroscopy

Figure 1: Linear absorption spectra for the slow bath (a) and fast bath (b)
parameter regimes computed via HEOM, pure state Ehrenfest, and projected
Ehrenfest.

Linear absorption spectra can be obtained from the Fourier transform of the
first-order response function, $\chi(t)$ Mukamel (1995),

$\begin{split}\chi(t)=\mathrm{Tr}[\hat{\mu}e^{-i\hat{H}t}\hat{\mu}(0)\rho(0)e^{i\hat{H}t}]\\\
\sigma(\omega)=\int_{0}^{\infty}dte^{i\omega t}\chi(t).\end{split}$ (12)

Applying the dipole operator at $t=0$ results in an optical coherence between
the ground state and singly excited states,
$\tilde{\rho}(0)=\hat{\mu}(0)\rho(0)$. This initial condition can be
represented either as an unmodified projected state (i.e., $\ket{e}\bra{g}$)
or as a sum of pure states. Absorption spectra computed via pure state
Ehrenfest, projected Ehrenfest, and the numerically exact HEOM for both the
fast and slow bath parameter regimes are shown in Fig. 1. In the slow bath
regime (Fig. 1 (a)), where Ehrenfest theory is expected to be accurate, both
projected Ehrenfest and pure state Ehrenfest quantitatively reproduce the HEOM
result. In the fast bath parameter regime (Fig. 1 (b)), despite qualitative
agreement with the HEOM result, there are small inaccuracies for both
Ehrenfest methods, with the pure state Ehrenfest method giving a more accurate
result than the projected Ehrenfest method. These results can be explained by
examining the dynamics arising from the coherence initial condition,
$\rho_{S}(0)=\ket{0}\bra{1}$ shown in Figs. 2 (a) and (b). $\ket{0}\bra{1}$ is
one of the coherence initial conditions (the others being
$\rho_{S}(0)=\ket{0}\bra{2}$ and their complex conjugates) that is propagated
to give rise to the first-order response function, $\chi(t)$. In the slow bath
regime, there is near quantitative agreement between both versions of
Ehrenfest and the HEOM result, while in the fast bath regime, there are more
pronounced differences, with the pure state Ehrenfest dynamics matching the
HEOM result more closely and the projected Ehrenfest result being underdamped.
These differences mirror and point to the source of inaccuracies in the
absorption spectra for the fast bath parameter regime. The overall similarity
in coherence dynamics across different methods (Fig. 2) also reveals why the
linear spectra appear to be relatively insensitive to the initialization
scheme used for Ehrenfest dynamics: the dynamics starting from a density
matrix corresponding to a coherence with the ground state are not as sensitive
to the initialization scheme. The relative insensitivity of the linear
absorption spectra to methods that inaccurately treat coherence dynamics has
previously been observed for Redfield theory, which was shown to yield
accurate linear absorption spectra despite getting incorrect population
dynamics Fetherolf and Berkelbach (2017), albeit with a different dynamical
method rooted in perturbation theory.

Figure 2: Coherence dynamics for the slow bath (a) and fast bath (b) parameter
regimes computed via HEOM, pure state Ehrenfest, and projected Ehrenfest.

###  IV.2 Nonlinear Optical Spectroscopy

Two-dimensional electronic spectra offer a much stricter test than absorption
spectroscopy of a method’s accuracy because they require one to correctly
capture both the population and coherence dynamics for a wider range of
initial conditions.

Figure 3: Exact 2DES spectra for the slow bath parameter regime computed using
HEOM (top row). 2DES spectra computed using the pure state Ehrenfest and
projected Ehrenfest methods are shown in the middle and bottom rows
respectively. All spectra are normalized such that the maximum amplitude
equals 1.

2DES spectra can be obtained from a Fourier transform over the rephasing
($R_{rp}$) and non-rephasing ($R_{nr}$) terms of third-order response function
under the rotating wave approximation Mukamel (1995),

$\begin{split}S(\omega_{3},t_{2},\omega_{1})&=\mathrm{Re}\int_{0}^{\infty}dt_{1}\int_{0}^{\infty}dt_{3}[e^{i(\omega_{1}t_{1}+\omega_{3})}R_{\rm{nr}}(t_{3},t_{2},t_{1})\\\
&+e^{i(-\omega_{1}t_{1}+\omega_{3})}R_{\rm{rp}}(t_{3},t_{2},t_{1})],\end{split}$
(13)

where

$\begin{split}R_{\rm{rp}}&=\Phi_{1}(t_{3},t_{2},t_{1})+\Phi_{2}(t_{3},t_{2},t_{1})-\Phi_{3}(t_{3},t_{2},t_{1}),\\\
R_{\rm{nr}}&=\Phi_{4}(t_{3},t_{2},t_{1})+\Phi_{5}(t_{3},t_{2},t_{1})-\Phi_{6}(t_{3},t_{2},t_{1})\end{split}$
(14)

include the ground state bleaching, stimulated emission, and excited state
absorption contributions, and

$\begin{split}\Phi_{1}(t_{3},t_{2},t_{1})&=\langle\mu_{-}\mathcal{G}(t_{3})[\mathcal{G}(t_{2})[\mu_{+}\mathcal{G}(t_{1})(\rho_{0}\mu_{-})]\mu_{+}]\rangle\\\
\Phi_{2}(t_{3},t_{2},t_{1})&=\langle\mu_{-}\mathcal{G}(t_{3})[\mu_{+}\mathcal{G}(t_{2})[\mathcal{G}(t_{1})(\rho_{0}\mu_{-})\mu_{+}]]\rangle\\\
\Phi_{3}(t_{3},t_{2},t_{1})&=\langle\mu_{-}\mathcal{G}(t_{3})[\mu_{+}\mathcal{G}(t_{2})[\mu_{+}\mathcal{G}(t_{1})(\rho_{0}\mu_{-})]]\rangle\\\
\Phi_{4}(t_{3},t_{2},t_{1})&=\langle\mu_{-}\mathcal{G}(t_{3})[\mathcal{G}(t_{2})[\mathcal{G}(t_{1})(\mu_{+}\rho_{0})\mu_{-}]\mu_{+}]\rangle\\\
\Phi_{5}(t_{3},t_{2},t_{1})&=\langle\mu_{-}\mathcal{G}(t_{3})[\mu_{+}\mathcal{G}(t_{2})[\mu_{-}\mathcal{G}(t_{1})(\mu_{+}\rho_{0})]]\rangle\\\
\Phi_{6}(t_{3},t_{2},t_{1})&=\langle\mu_{-}\mathcal{G}(t_{3})[\mu_{+}\mathcal{G}(t_{2})[\mathcal{G}(t_{1})(\mu_{+}\rho_{0})]\mu_{-}]\rangle.\end{split}$
(15)

Here, $\mathcal{G}(t)\rho=e^{-iH_{\rm{mat}}t/\hbar}\rho
e^{iH_{\rm{mat}}t/\hbar}$ is the Liouville-space propagator and

$\begin{split}\mu_{-}&=\sum_{\rm{m}}\mu_{\rm{m}}\ket{0}\bra{m}\\\
\mu_{+}&=\sum_{\rm{m}}\mu_{\rm{m}}\ket{m}\bra{0}.\end{split}$ (16)

A procedure for obtaining the third-order response function from Ehrenfest
dynamics was previously outlined in Ref. van der Vegte  _et al._ , 2013, which
employs ensemble averages over the ground state wavefunction, $\ket{0}$. In
this scheme, the dipole operator is applied via a product with individual
dipole matrix elements, time propagation is done only over select population
states, and the effect of coherences is modeled via an average over two
wavefunctions formed from the ket and bra. However, the use of specific dipole
matrix elements as well as time-propagation over select population states make
it difficult to generalize that method to arbitrary systems with multiple
quantum states. Here, we present a generalizable method that exploits the
linearity of density matrices to compute third-order response functions as
follows:

  1. 1.

Apply the first dipole operator to $\rho_{0}$ at $t=0$ and split the resultant
density matrix into 4 pure states as described in Section II. Via the Condon
approximation, the bath degrees of freedom are untouched by the dipole
operator and are initially sampled from the Wigner transform of the Boltzmann
distribution on the ground state. Each resultant pure state inherits its own
copy of the bath. Propagate states and their corresponding baths independently
through $t_{1}$.

  2. 2.

Apply the second dipole operator at $t=t_{1}$, this time to all 4 propagated
states. This operation usually results in a new set of non-pure states, which
are subsequently split into 4 of their respective pure states, bringing the
total number of branches to 16. Each branch adopts a copy of bath coordinates
and momenta from its parent branch, ensuring continuity from $t=0$. Propagate
the pure states and baths independently through the waiting time, $t_{2}$.

  3. 3.

Apply the third dipole operator at $t=t_{1}+t_{2}$ to the 16 propagated
states. This splits them again into a total of 64 branches, each a pure state.
Propagate all branches through $t_{3}$.

  4. 4.

Consolidate all 64 branches at $t=t_{1}+t_{2}+t_{3}$ and apply the final
dipole operator to obtain the third-order response function.

Multiple instances of this procedure are run in order to average over an
ensemble of initial bath states. From the above, it is clear that the pure
state Ehrenfest method requires more sampling than projected Ehrenfest, which
does not involve any splitting of initial conditions into pure states. As we
will see below, this splitting into pure states is required to obtain accurate
2DES spectra. The procedure can also be trivially parallelized since all
trajectories in the scheme run independently. Additionally, as demonstrated in
Section II, we decompose each initial density matrix $\ket{a}\bra{b}$ as a sum
of only 4 pure states in the same Liouville space.

An alternative method for obtaining accurate Ehrenfest dynamics from non-pure
state density matrices is outlined in our previous work in the Appendices of
Ref. Pfalzgraff  _et al._ , 2019, and it involves stochastically sampling an
auxiliary wavefunction. For example, the initial condition $\ket{0}\bra{1}$
would be obtained from

$\ket{\psi}=\frac{1}{\sqrt{2}}\left(\ket{0}+e^{i\phi}\ket{1}\right),$ (17)

where $\phi$ is a parameter that is sampled over the range $0-2\pi$. In this
stochastic scheme, sampling is required both over values of $\phi$ and initial
positions and momenta of the bath. In the context of a third-order response
function, $\phi$ would require fresh sampling every time a measurement is made
(i.e., at times $t=0$, $t=t_{1}$, and $t=t_{1}+t_{2}$) because the auxiliary
wavefunction would otherwise return incorrect dynamics. In contrast, the pure
state Ehrenfest method detailed above only requires the usual sampling over
initial bath conditions, which can be done once at $t=0$ since the baths are
continuous through $t=t_{1}+t_{2}+t_{3}$.

Figure 4: Real parts of $\Phi_{2}$ for the slow bath parameter regime as
computed via HEOM, pure state Ehrenfest, and projected Ehrenfest. Figure 5:
Imaginary parts of $\Phi_{2}$ for the slow bath parameter regime as computed
via HEOM, pure state Ehrenfest, and projected Ehrenfest.

Fig. 3 shows the 2DES spectra obtained from both the pure state and projected
forms of Ehrenfest and how they compare to the HEOM spectra. The response
functions used to obtain the Ehrenfest spectra were truncated (windowed using
a step function) as outlined in the SI Sec. LABEL:sec:third_order_truncation.
From this, we observe that pure state Ehrenfest recovers almost quantitative
agreement with the HEOM result, with just a slight deterioration as the
waiting time $t_{2}$ increases. In contrast, projected Ehrenfest yields much
worse agreement with the HEOM spectrum. While most of the spectral features
from the HEOM spectrum are present at $t_{2}=0$, the projected Ehrenfest
spectrum becomes progressively worse as $t_{2}$ increases and bears little
resemblance to the HEOM result by the time $t_{2}=600$ fs.

The poor performance of the projected Ehrenfest method can be tracked down to
discrepancies in the contributions of the response function. Figs. 4 and 5
illustrate this using plots of $\Phi_{2}$, with the rest of the contributions
of the response function being shown in SI Figs. LABEL:fig:phi_1_slow\-
LABEL:fig:phi_6_fast. Here, we observe that, as expected, pure state Ehrenfest
accurately reproduces the HEOM result while projected Ehrenfest fails to do
so, and that the discrepancy is much more striking in the imaginary part of
the response function. These discrepancies can be further traced to the
single-time dynamics obtained from coherence initial conditions between
excited states. Since the Frenkel exciton model is in the global ground state
$\rho_{0}=\ket{0}\bra{0}$ at $t=0$, such coherence states can only be obtained
when the dipole operator is applied multiple times, as is the case when
computing the third-order response. As such, the dynamics of coherence initial
conditions between excited states (e.g. $\ket{1}\bra{2}$) are relevant to
nonlinear spectra but not linear spectra. Figure 6 shows both the population
and coherence dynamics for $\ket{1}\bra{2}$. In both the slow and fast bath
parameter regimes, we see from the first two rows that the population dynamics
obtained from pure state Ehrenfest more closely match the HEOM result,
especially at longer times. In contrast, the projected Ehrenfest result
matches the HEOM result until $\sim$200 fs, after which there are notable
deviations. This poor performance of projected Ehrenfest is also observed to a
lesser extent in the coherence dynamics of the $\ket{1}\bra{2}$ initial state
(bottom two rows of Fig. 6). We emphasize that the coherence dynamics
presented here are distinct from those involved in the absorption spectra
since they do not involve the ground state. The worsening performance of
projected Ehrenfest with time in reproducing the dynamics arising from
$\ket{1}\bra{2}$ also explains why the respective 2DES spectrum degrades
progressively as the waiting time $t_{2}$ increases (Fig. 3). As such, the
differences in the performance between projected and pure state Ehrenfest can
be attributed to inaccuracies in the projected Ehrenfest dynamics when
populations, and coherences to a lesser extent, arise from an initial
condition of coherences between two excited states.

Figure 6: Single-time dynamics for an initial condition of $\ket{1}\bra{2}$
for both the slow (left) and fast (right) bath parameter regimes. The first
two rows contain population dynamics (i.e. $\ket{1}\bra{1}$ and
$\ket{2}\bra{2}$), while the bottom two rows show coherence dynamics
($\ket{2}\bra{1}$ and $\ket{1}\bra{2}$).

Pump-probe spectra, generated by integrating over $\omega_{1}$, are shown in
Fig. 7. The results here mirror those of the 2DES spectrum, where the pure
state Ehrenfest method accurately reproduces the HEOM result, albeit with
slightly higher peaks at 200 cm-1 for $t_{2}\geq$ 200 fs. In contrast, while
the projected Ehrenfest method yields accurate results at $t_{2}=0$, the
accuracy of the pump-probe spectra degrades at $t_{2}\geq$ 200 fs, with the
peak around 200 cm-1 losing definition.

Figure 7: Pump-probe spectra for the slow bath parameter regime obtained from
integrating the 2DES spectra over $\omega_{1}$.

Next we examine how well pure state Ehrenfest performs in a more challenging
fast bath regime. Fig. 8 shows 2DES spectra for the fast bath parameters
computed via HEOM and pure state Ehrenfest. The range of values for $t_{2}$ (0
- 150 fs) is more limited here than in the slow bath regime because we
observed virtually no change in the 2DES spectrum beyond $t_{2}=$200 fs (see
SI Figure LABEL:fig:fast_bath_extended). As with the slow bath parameter
regime, the pure state Ehrenfest response functions were truncated before
being Fourier transformed to produce the 2DES spectra (see SI Sec.
LABEL:sec:third_order_truncation). In this more challenging regime, pure state
Ehrenfest does not perform as well as it did for the slow bath parameters.
This is not surprising, as it aligns with the well-known deficiencies of the
Ehrenfest method, which struggles to capture the correct dynamics in the limit
of fast nuclei. Nonetheless, it is able to qualitatively capture most of the
features present in the HEOM spectrum. This result is corroborated by the
pump-probe spectra shown in Fig. 9.

Figure 8: Exact two-dimensional spectra for the fast bath parameters computed
using HEOM (top row), compared to the pure state Ehrenfest method (bottom
row). All spectra are normalized such that the maximum amplitude equals 1.
Figure 9: Pump-probe spectra for the fast bath parameter regime obtained from
integrating the 2DES spectra over $\omega_{1}$.

##  V Conclusion

Here we have introduced the pure state Ehrenfest method to obtain linear and
nonlinear spectra. We have shown that exploiting the linearity of the density
matrix to express all initial conditions as sums of pure states is essential
to obtaining accurate coherence dynamics. In particular, while our approach
provides modest improvements over projected Ehrenfest in computing linear
spectra, it is able to dramatically improve the accuracy with which nonlinear
spectra (2DES and pump-probe) can be obtained, especially at longer waiting
times. The physical origin of this improvement arises from the ability of our
pure state Ehrenfest approach to much more accurately capture the dynamics
when the initial condition is a coherence between excited states. These
conditions do not occur when calculating the linear spectra (where all the
initial density matrices are coherences with the ground state) but play a
vital role in obtaining the higher-order response functions needed to capture
nonlinear spectroscopy. We have shown that pure state Ehrenfest gives
excellent agreement with the exact linear and 2DES spectra in the slow bath
regime where the Ehrenfest method is expected to be accurate and have also
demonstrated that it still provides good results in a faster bath regime and
preserves the qualitative features observed in the exact spectra. These
results suggest that the pure state Ehrenfest approach provides a tractable
and accurate approach to calculate nonlinear spectra for multichromophoric
systems in a wide range of chemical environments.

## References

## References

  * Ginsberg, Cheng, and Fleming (2009) N. S. Ginsberg, Y.-C. Cheng, and G. R. Fleming, Accounts of Chemical Research 42, 1352 (2009).
  * Cho (2008) M. Cho, Chemical Reviews 108, 1331 (2008).
  * Khalil, Demirdöven, and Tokmakoff (2003) M. Khalil, N. Demirdöven, and A. Tokmakoff, The Journal of Physical Chemistry A 107, 5258 (2003).
  * Hybl  _et al._ (1998) J. D. Hybl, A. W. Albrecht, S. M. G. Faeder, and D. M. Jonas, Chemical Physics Letters 297, 307 (1998).
  * Mukamel (1995) S. Mukamel,  _Principles of Nonlinear Optical Spectroscopy_, Oxford series in optical and imaging sciences (Oxford University Press, 1995).
  * Brixner  _et al._ (2005) T. Brixner, J. Stenger, H. M. Vaswani, M. Cho, R. E. Blankenship, and G. R. Fleming, Nature 434, 625 (2005).
  * Engel  _et al._ (2007) G. S. Engel, T. R. Calhoun, E. L. Read, T.-K. Ahn, T. Mančal, Y.-C. Cheng, R. E. Blankenship, and G. R. Fleming, Nature 446, 782 (2007).
  * Dostál  _et al._ (2012) J. Dostál, T. Mančal, R. nas Augulis, F. Vácha, J. Pšenčík, and D. Zigmantas, Journal of the American Chemical Society 134, 11611 (2012).
  * Abramavicius and Mukamel (2010) D. Abramavicius and S. Mukamel, The Journal of Chemical Physics 133, 184501 (2010).
  * Tanimura and Kubo (1989) Y. Tanimura and R. Kubo, Journal of the Physical Society of Japan 58, 101 (1989).
  * Ishizaki and Tanimura (2007) A. Ishizaki and Y. Tanimura, The Journal of Physical Chemistry A 111, 9269 (2007).
  * Kreisbeck  _et al._ (2011) C. Kreisbeck, T. Kramer, M. Rodríguez, and B. Hein, Journal of Chemical Theory and Computation 7, 2166 (2011).
  * Chen  _et al._ (2010) L. Chen, R. Zheng, Q. Shi, and Y. Yan, The Journal of Chemical Physics 132, 024505 (2010).
  * Meyer (2011) H.-D. Meyer, Wiley Interdisciplinary Reviews: Computational Molecular Science 2, 351 (2011).
  * Schulze  _et al._ (2017) J. Schulze, M. F. Shibl, M. J. Al-Marri, and O. Kühn, Chemical Physics 497, 10 (2017).
  * REDFIELD (1965) A. REDFIELD, in  _Advances in Magnetic Resonance_ (Elsevier, 1965) pp. 1–32.
  * Bloch (1957) F. Bloch, Physical Review 105, 1206 (1957).
  * Fetherolf and Berkelbach (2017) J. H. Fetherolf and T. C. Berkelbach, The Journal of Chemical Physics 147, 244109 (2017).
  * Montoya-Castillo, Berkelbach, and Reichman (2015) A. Montoya-Castillo, T. C. Berkelbach, and D. R. Reichman, The Journal of Chemical Physics 143, 194108 (2015).
  * Meyer and Miller (1979) H.-D. Meyer and W. H. Miller, The Journal of Chemical Physics 70, 3214 (1979).
  * Stock and Thoss (1997) G. Stock and M. Thoss, Physical Review Letters 78, 578 (1997).
  * Huo and Coker (2011) P. Huo and D. F. Coker, The Journal of Chemical Physics 135, 201101 (2011).
  * Mannouch and Richardson (2020a) J. R. Mannouch and J. O. Richardson, The Journal of Chemical Physics 153, 194110 (2020a).
  * Mannouch and Richardson (2020b) J. R. Mannouch and J. O. Richardson, The Journal of Chemical Physics 153, 194109 (2020b).
  * Provazza  _et al._ (2018) J. Provazza, F. Segatta, M. Garavelli, and D. F. Coker, Journal of Chemical Theory and Computation 14, 856 (2018).
  * Mannouch and Richardson (2022) J. R. Mannouch and J. O. Richardson, The Journal of Chemical Physics 156, 024108 (2022).
  * Gao and Geva (2020) X. Gao and E. Geva, Journal of Chemical Theory and Computation 16, 6491 (2020).
  * Polley and Loring (2021) K. Polley and R. F. Loring, The Journal of Chemical Physics 154, 194110 (2021).
  * Loring (2022) R. F. Loring, Annual Review of Physical Chemistry 73, 273 (2022).
  * Jansen and Knoester (2006) T. L. C. Jansen and J. Knoester, The Journal of Physical Chemistry B 110, 22910 (2006).
  * Liang and Jansen (2012) C. Liang and T. L. C. Jansen, Journal of Chemical Theory and Computation 8, 1706 (2012).
  * Torii (2006) H. Torii, The Journal of Physical Chemistry A 110, 4822 (2006).
  * Jansen, Zhuang, and Mukamel (2004) T. L. C. Jansen, W. Zhuang, and S. Mukamel, The Journal of Chemical Physics 121, 10577 (2004).
  * Jansen  _et al._ (2005) T. L. C. Jansen, T. Hayashi, W. Zhuang, and S. Mukamel, The Journal of Chemical Physics 123, 114504 (2005).
  * Tempelaar  _et al._ (2013) R. Tempelaar, C. P. van der Vegte, J. Knoester, and T. L. C. Jansen, The Journal of Chemical Physics 138, 164106 (2013).
  * van der Vegte  _et al._ (2013) C. P. van der Vegte, A. G. Dijkstra, J. Knoester, and T. L. C. Jansen, The Journal of Physical Chemistry A 117, 5970 (2013).
  * Kapral and Ciccotti (1999) R. Kapral and G. Ciccotti, The Journal of Chemical Physics 110, 8919 (1999).
  * Montoya-Castillo and Reichman (2016) A. Montoya-Castillo and D. R. Reichman, The Journal of Chemical Physics 144, 184104 (2016).
  * Ishizaki and Fleming (2009) A. Ishizaki and G. R. Fleming, The Journal of Chemical Physics 130, 234111 (2009).
  * Berkelbach, Reichman, and Markland (2012) T. C. Berkelbach, D. R. Reichman, and T. E. Markland, The Journal of Chemical Physics 136, 034113 (2012).
  * Berkelbach  _et al._ (2020) T. Berkelbach, J. Fetherolf, P. Shih, and Ian S Dunn, “berkelbach-group/pyrho v1.0,”  (2020).
  * Pfalzgraff  _et al._ (2019) W. C. Pfalzgraff, A. Montoya-Castillo, A. Kelly, and T. E. Markland, The Journal of Chemical Physics 150, 244109 (2019).

## Acknowledgements

This work was funded by the U.S. Department of Energy, Office of Science,
Office of Basic Energy Sciences (DE-SC0020203 to T.E.M.). A.O.A. acknowledges
support from the Edward Curtis Franklin Award and the Stanford Diversifying
Academia, Recruiting Excellence Fellowship. A.M.C. acknowledges the start-up
funds from the University of Colorado, Boulder.