# The Classical Equations of Motion of Quantized Gauge Theories, Part 2:
Electromagnetism

David E. Kaplan  Department of Physics & Astronomy, The Johns Hopkins
University, Baltimore, MD 21218, USA Tom Melia  Kavli IPMU (WPI), UTIAS, The
University of Tokyo, Kashiwa, Chiba 277-8583, Japan  Surjeet Rajendran
Department of Physics & Astronomy, The Johns Hopkins University, Baltimore, MD
21218, USA

###### Abstract

In this and companion papers, we show that quantum field theories with gauge
symmetries permit a broader class of classical dynamics than typically
assumed. In this article, we show that the quantization of electromagnetism
permits the existence of classical electric field states that do not obey
Gauss’s law. These states are gauge invariant and their time evolution can be
consistently described using the Schrödinger equation. The time evolution of
these states is such that at the classical level, the full set of Maxwell’s
equations would appear to hold, with the physical effects of these states
being attributable to an auxiliary, static “shadow” charge density with no
internal degrees of freedom. This density could affect the dynamics of charged
particles in our universe and it may thus be of observational interest.

###### Contents

  1. 1 Introduction
  2. 2 Executive Summary
  3. 3 Weyl Gauge
    1. 3.1 Canonical Formulation
    2. 3.2 Path Integral Formulation
    3. 3.3 Gauss’s Law and Its Violation
  4. 4 Other Gauges
  5. 5 Gravitation and Cosmology
  6. 6 Discussion

##  1 Introduction

Quantum electrodynamics, with its better than part-per-billion agreement with
experiment, is an extraordinarily successful theory. The central claim of this
paper is that despite this success, the conventional quantization procedure
used to describe the quantum mechanics of electromagnetism has overlooked a
theoretically consistent possibility that likely has interesting observational
and theoretical implications. This possibility is the existence of quantum
states that violate Gauss’s law, but can nevertheless be consistently time
evolved using the Schrödinger equation in a gauge-invariant manner. The
existence of these states is tied to quantum mechanics. In classical
electrodynamics, since Gauss’s law is one of Maxwell’s equations, such a state
cannot be described. However, unlike the classical equations of Maxwell, time
evolution in quantum mechanics is described by a single equation, namely, the
Schrödinger equation.

When a quantum state evolves, the expectation values of various quantum
mechanical operators in that state will automatically obey corresponding
classical equations. Since the Schrödinger equation is a dynamical equation,
it can force the expectation values of quantum states to obey the dynamical
equations of the corresponding classical theory, i.e., Ampere’s law in the
case of electromagnetism. But Gauss’s law is not a dynamical equation - it is
a constraint on the initial states of the theory and is thus not a consequence
of the Schrödinger equation. The initial state of a quantum theory can be any
state that can be described and time evolved in a gauge invariant manner. The
main point of this paper is to show that there are such initial states that
violate Gauss’s law but still lead to gauge invariant physics. We summarize
these arguments in Section 2.

In a companion paper [1], we pointed out that similar states also exist in
General Relativity (also see [2]). In that paper, since the quantum theory of
General Relativity has been poorly understood [3, 4], we described this point
in detail. In this paper, our interest is electromagnetism - a theory whose
quantization is well understood (as it appears in standard textbooks, e.g.,
[5, 6, 7]). We will thus focus on the specific arguments in the quantization
of electromagnetism to pinpoint how states that violate Gauss’s law can also
be consistently described in this theory. The existence of these states is
most easily understood in the Weyl gauge - we will thus begin our discussion
in this gauge in Section 3. Following this, we show how such states can also
be described in other gauges (such as the Coulomb gauge) in Section 4. We then
describe the gravitational and cosmological implications of such states in
Section 5 and conclude with a discussion about the phenomenological relevance
of these states and potential experimental ways to constrain or discover them.

Throughout we use Greek letters for four-vector indices, Latin letters for
three-vector (spatial) indicies, and overdots for coordinate time derivatives.
Coordinates ${\bf x}$ and vectors in bold refer to spatial coordinates and
vectors respectively while coordinates $x$ in regular fonts represent space-
time coordinates. Throughout this paper, when we describe Hamiltonians, we
will work in the Schrödinger picture where the operators are time independent
and the states evolve in time.

##  2 Executive Summary

The physics of electromagnetism is contained in the following Lagrangian:

$\mathcal{L_{EM}}=-\frac{1}{4}F_{\mu\nu}F^{\mu\nu}+A_{\mu}J^{\mu}+\mathcal{L}_{J}$
(1)

which describes the interactions of the electromagnetic potential $A_{\mu}$
with a current $J^{\mu}$.111We have absorbed the electromagnetic coupling
$e_{c}$ in the definition of the current $J^{\mu}$. The other interactions of
this current (such as its kinetic terms) are described by $\mathcal{L}_{J}$.
When we perform a gauge transformation $A_{\mu}\rightarrow
A_{\mu}+\partial_{\mu}\alpha$, the kinetic term of this current transforms
covariantly so that the Lagrangian (1) is invariant. At the classical level,
one obtains the equations of motion by varying the action:

$\mathcal{S_{EM}}=\int d^{4}x\,\mathcal{L_{EM}}$ (2)

with respect to the gauge field $A_{\mu}$ and setting the variations
$\frac{\partial\mathcal{S_{EM}}}{\partial A_{\mu}}$ to zero. Accordingly, we
get the following classical equations for the electric (${\bf E}$) and
magnetic (${\bf B}$) fields:

${\bf\nabla}\cdot{\bf E}=J^{0}$ (3)

otherwise known as Gauss’s law (obtained by varying $A_{0}$), and

${\bf\nabla}\times{\bf B}-\frac{d{\bf E}}{dt}={\bf J}$ (4)

otherwise known as Ampere’s law (obtained by varying $A_{i}$).

Even at the classical level, observe the following fact. Gauss’s law (3) is a
constraint on the initial values of the electric field ${\bf E}$ whereas
Ampere’s law (4) is a dynamical equation. Suppose we ignore (3) and solve (4)
with some initial values for ${\bf E}$ and ${\bf B}$ and the currents
$J^{\mu}$ (whose dynamical evolution is governed by its own equation of
motion, such as the Lorentz force law).222While these are second order
equations for $A_{\mu}$, they are first order equations for ${\bf E}$ and
${\bf B}$ \- thus given initial conditions for ${\bf E}$ and ${\bf B}$ there
is no mathematical obstacle to solving these dynamical equations. What does
this evolution look like? By taking a divergence of (4) and using the fact
that the dynamical equations of the current $J^{\mu}$ will enforce current
conservation $\partial_{\mu}J^{\mu}=0$, we see that the following identity
holds:

$\frac{d\left({\bf\nabla}\cdot{\bf E}-J^{0}\right)}{dt}=0$ (5)

Thus, if the initial value of ${\bf E}$ were such that it obeyed Gauss’s law
$({\bf\nabla}\cdot{\bf E}-J^{0})=0$, the time evolution generated by (4) will
automatically continue to enforce it. But, suppose we took an initial state of
${\bf E}$ so that Gauss’s law was violated by a function of only space,
$({\bf\nabla}\cdot{\bf E}-J^{0})=J_{s}^{0}\neq 0$. Then, the evolution would
be such that:

$\frac{d\left({\bf\nabla}\cdot{\bf
E}-J^{0}\right)}{dt}=\frac{dJ_{s}^{0}}{dt}=0$ (6)

That is, the state looks like the time evolution of a system where in addition
to the known current $J^{\mu}$, there is an additional “shadow” charge density
$J_{s}^{0}$ that is somehow unmovable. The interesting fact about a state that
violates Gauss’s law is that even at the classical level it is described in
terms of the gauge invariant observable ${\bf E}$ \- there is thus no logical
issue involved in time evolving such a state.

At the classical level, we would reject such states simply because we believe
in Gauss’s law and would thus require $({\bf\nabla}\cdot{\bf E}-J^{0})=0$.
But, classical mechanics is not the correct description of nature - the
underlying theory is quantum mechanics and classical physics is a limit of
quantum mechanics. The key question that we need to ask is if Gauss’s law
follows from quantum mechanics. We argue that it does not. Instead, we show
that quantum mechanics allows for the existence of gauge invariant states of
electromagnetism that violate Gauss’s law. At the classical level, the time
evolution of these states would be identical to that of states containing an
unmovable shadow charge $J_{s}^{0}$ described in (6) - but there is no new
degrees of freedom associated with $J_{s}^{0}$. It is simply a state of
electromagnetism. Further, these states can also be consistently coupled to
gravity.

How could Gauss’s law not be true in quantum mechanics but appear to be true
in classical electromagnetism? As a first peek at this issue, observe the
following. In classical physics, we obtained Maxwell’s equations by varying
the action $\mathcal{S}_{EM}$ along four independent variations of the
potential $A_{\mu}$. But, due to gauge redundancy, there are not four
independent variations of $A_{\mu}$. By Helmholtz’s theorem, any potential
$A_{\mu}$ can be decomposed as $A_{\mu}=K_{\mu}+\partial_{\mu}\alpha$ where
$K_{\mu}$ is divergence-less, $\partial^{\mu}K_{\mu}=0$, and thus only
contains three degrees of freedom. Write the classical action
$\mathcal{S}_{EM}$ in terms of $K_{\mu}$ and $\alpha$ instead of $A_{\mu}$.
Due to the gauge invariance of the action $\mathcal{S}_{EM}$ under the gauge
transformations $A_{\mu}\rightarrow A_{\mu}+\partial_{\mu}\alpha$ and the
associated covariant transformations on $\mathcal{L}_{J}$, the action:

$\mathcal{S_{EM}}=\int
d^{4}x\,\mathcal{L_{EM}}\left(A_{\mu},\partial_{\nu}A_{\mu}\right)=\int
d^{4}x\,\mathcal{L_{EM}}\left(K_{\mu},\partial_{\nu}K_{\mu}\right)$ (7)

is only a function of $K_{\mu}$. But since the divergence-less four vector
$K_{\mu}$ only has three degrees of freedom, we do not have four independent
variations to obtain four independent equations.

The reader will observe that in writing (7), we have effectively picked the
Lorentz gauge where we set $\partial^{\mu}A_{\mu}=0$ and as a consequence, we
naively lost an equation of motion. The definition of the quantum theory,
either at the level of the Hamiltonian or the path integral, requires us to
pick a gauge in order to define operators and states. It is this picking of
operators and states which affords additional freedom in the derived classical
equations. This results in loosening restrictions on the allowed quantum
states of the theory, as we’ll see, permitting states analogous to (6) that
violate Gauss’s law.

##  3 Weyl Gauge

To quantize electromagnetism in Weyl (or temporal) gauge, the following
procedure is adopted [6] to specify the operators, Hamiltonian and physical
states. First, we set $A_{0}=0$ in the classical Lagrangian. The spatial
components $A_{j}$ of the vector potential and its conjugate momentum
$\Pi_{j}=\frac{\partial\mathcal{L_{EM}}}{\partial A_{j}}=-E_{j}$ get promoted
to operators $\hat{A}_{j},\hat{\Pi}_{j}$ with canonical equal time commutation
relations:

$[\hat{A}_{j}\left(\bf{x}\right),\hat{\Pi}_{j^{\prime}}\left(\bf{x^{\prime}}\right)]=i\,\delta\left({\bf
x-x^{\prime}}\right)\delta_{jj^{\prime}}$ (8)

Using $\Pi_{j}=-E_{j}$, we have

$[\hat{A_{j}}\left(\bf{x}\right),\hat{E_{j^{\prime}}}\left(\bf{x^{\prime}}\right)]=-i\,\delta\left({\bf
x-x^{\prime}}\right)\delta_{jj^{\prime}}$ (9)

The Hamiltonian constructed from these operators is:

$\hat{H}_{W}=\int d^{3}{\bf
x}\,\left(\frac{1}{2}\left(\hat{\bf{E}}\cdot\hat{\bf{E}}+\hat{\bf{B}}\cdot\hat{\bf{B}}\right)+\hat{\bf{J}}\cdot\hat{\bf{A}}+\hat{\cal
H}_{J}\right)$ (10)

where $\hat{\bf{B}}\equiv{\bf\nabla}\times\hat{\bf{A}}$ and $\hat{\cal H}_{J}$
is the remaining Hamiltonian density of the degrees of freedom in the current
${\bf J}$. If a quantum state $|\Psi\rangle$ obeys the Schrödinger equation

$i\frac{\partial|\Psi\rangle}{\partial t}=\hat{H}_{W}|\Psi\rangle$ (11)

what equations of motion are automatically obeyed by the expectation values of
various operators such as $\langle\Psi|{\bf E}|\Psi\rangle$? In the next two
subsections, using the canonical and path integral formulations, we will show
that Ampere’s law follows from the quantum dynamics while Gauss’s law does
not. Following this discussion, we will show how Gauss’s law is obtained in
the Weyl gauge and argue that this prescription permits a broader class of
quantum states than previously considered.

###  3.1 Canonical Formulation

In the canonical language (see, for example, [5]), the classical equations
should be true in expectation value (Ehrenfest’s theorem). This can be seen by
taking a time derivative of expectation values of fields and using the
Schrödinger equation:

$\displaystyle\partial_{t}\langle\hat{E}^{j}({\bf x})\rangle$ $\displaystyle=$
$\displaystyle i\langle\left[\hat{H}_{W},\hat{E}^{j}({\bf x})\right]\rangle$
(12) $\displaystyle=$ $\displaystyle\int d^{3}{\bf
x}^{\prime}\,\langle\frac{1}{2}\left[({\bf\nabla}\times\hat{\bf{A}}({\bf
x}^{\prime}))^{2},{\hat{E}}^{j}({\bf x})\right]+\left[\hat{\bf{A}}({\bf
x}^{\prime}),\hat{{E}}^{j}({\bf x})\right]\cdot\hat{\bf{J}}({\bf
x}^{\prime})\,\rangle$ $\displaystyle=$
$\displaystyle\langle({\bf\nabla}\times\hat{\bf{B}})^{j}({\bf
x})\rangle-{\langle\hat{{J}}^{j}({\bf x})\,\rangle}$

where we integrated by parts to generate the last step. We see here that
Ampere’s law is reproduced in expectation value. This, and
$\partial_{t}\langle\hat{\bf{A}}\rangle=-\langle\hat{\bf{E}}\rangle$, are the
only equations of motion for the electromagnetic field predicted by this
method.

But what about Gauss’s law? The absence of Gauss’s law reflects the missing
conjugate momentum of $A_{0}$, and thus the missing equation of motion. As in
the classical case, a divergence of (12) reproduces the time-derivative of the
expectation value of Gauss’s law (using
$\langle\partial_{\mu}J^{\mu}\rangle=0$):

$\partial_{t}\langle{\bf\nabla}\cdot\hat{\bf
E}\rangle=-\langle{\bf\nabla}\cdot\hat{\bf
J}\rangle=\partial_{t}\langle\hat{J}^{0}\rangle$ (13)

At this point, the case is usually made to impose Gauss’s law by fixing the
remaining spatial gauge invariance. Classically, in Weyl gauge, the spatial
transformations ${\bf A}({\bf x})\rightarrow{\bf A}({\bf
x})+{\bf\nabla}\alpha({\bf x})$ are still an invariance of the action. Quantum
mechanically, the generator of this transformation is the operator
$\hat{G}\equiv{\bf\nabla}\cdot\hat{\bf E}-\hat{J}^{0}$. Not surprisingly, it
can be checked that this operator commutes with the Hamiltonian:

$\left[\hat{G},\hat{H}_{W}\right]=0$ (14)

which means this operator and the Hamiltonian can be simultaneously
diagonalized. One can now impose Gauss’s law by requiring physical state
vectors to be invariant under spatial gauge transformations, namely:

$\displaystyle\hat{G}\,|\Psi_{EM}\rangle=0$ (15) $\displaystyle e^{-i\int
d^{3}{\bf x}\,\alpha\,\hat{G}}\,|\Psi_{EM}\rangle=|\Psi_{EM}\rangle$ (16)

and thus all physical states are constrained to obey Gauss’s law. Thus we see,
to obtain Gauss’s law, a constraint equation, from the quantum field theory,
one has to impose it by hand on the physical states. We will see this again
below in the path-integral formulation. Then we will see that this constraint
on states is not a requirement for a consistent quantum theory.

###  3.2 Path Integral Formulation

It is useful to reproduce this in the path integral language (see, for
example, [6]). Let us naively construct the path integral that solves (11) in
the field basis:

$\mathcal{T}=\langle{\bf A}_{f}|T\left(t_{2};t_{1}\right)|{\bf
A}_{i}\rangle=\int_{{\bf A}\left(t_{1}\right)={\bf A}_{i}}^{{\bf
A}\left(t_{2}\right)={\bf
A}_{f}}\,DA_{\mu}\,D\lambda\,e^{i\int_{t_{1}}^{t_{2}}d^{4}x\,\left(\mathcal{L_{EM}}-\lambda
A_{0}\right)}$ (17)

where the Lagrange multiplier $\lambda$ enforces the Weyl gauge $A_{0}=0$.
This path integral yields the transition matrix element for the field basis
state $|{\bf A}_{i}\rangle$ at time $t_{1}$ to evolve to $|{\bf A}_{f}\rangle$
at time $t_{2}$. This transition matrix element (17) should be invariant when
we compute this path integral with a variable redefinition $A_{\mu}\rightarrow
A_{\mu}+\delta A_{\mu}$ with $\delta A_{\mu}$ vanishing at the boundaries.
This yields a set of Schwinger-Dyson equations which show how the classical
field equations arise as identities automatically obeyed by the expectation
values of various field operators when the quantum state evolves as per (17).
One can check that when this procedure is applied to the spatial variations
$\delta A_{i}$, one obtains the result that Ampere’s law is obeyed by the
expectation values of the quantum operators. But, for the variation $\delta
A_{0}$, this yields the equation:

$\langle\Psi|{\bf\nabla}\cdot{\bf E}-J^{0}+\lambda|\Psi\rangle=0$ (18)

This is not an equation of motion or a constraint on the physical state
$|\Psi\rangle$ \- instead, it describes how the unphysical Lagrange multiplier
$\lambda$ evolves in the path integral to maintain the gauge $A_{0}=0$. Thus
Gauss’s law does not immediately follow from the quantum Hamiltonian
$\hat{H}_{W}$ .

###  3.3 Gauss’s Law and Its Violation

The origin of Gauss’s law in the Weyl gauge is tied to the elimination of
residual spatial gauge transformations in the theory. The gauge choice
$A_{0}=0$ does not eliminate all the gauge freedom in the theory - in
principle we would like to identify quantum states $|{\bf A}\rangle$ that are
related to each other by purely spatial gauge transformations: $|{\bf
A}\rangle\equiv|{\bf A}+{\bf\nabla}\alpha\rangle$ for any $\alpha\left({\bf
x}\right)$. For this equivalence to hold, we need the following to be true:

$\langle{\bf A}_{f}|T\left(t_{2};t_{1}\right)|{\bf A}_{i}\rangle=\langle{\bf
A}_{f}+{\bf\nabla}\alpha_{f}|T\left(t_{2};t_{1}\right)|{\bf
A}_{i}+{\bf\nabla}\alpha_{i}\rangle$ (19)

with $\alpha_{i}\left({\bf x}\right)\neq\alpha_{f}\left({\bf x}\right)$. It
can be checked that when $\alpha_{i}\left({\bf
x}\right)\neq\alpha_{f}\left({\bf x}\right)$, the path integral (17) does not
maintain this equality, violating the equivalence $|{\bf A}\rangle\equiv|{\bf
A}+{\bf\nabla}\alpha\rangle$. To maintain this equivalence, the physical
Hilbert space of the theory is restricted to a smaller space of states that
are invariant under these spatial gauge transformations and the time evolution
operator (17) is projected onto this restricted space, leading to physics that
is also invariant under such general spatial gauge transformations.

Let us see how this works. We want the physics of the states $|{\bf A}\rangle$
to be identical to that of the states $|{\bf A}+{\bf\nabla}\alpha\rangle$. As
we’ve seen, these spatial gauge transformations are generated by the operator
$\hat{G}={\bf\nabla}\cdot\hat{{\bf E}}-\hat{J}^{0}$. Thus, under a spatial
gauge transformation, the eigenstates of this operator will transform with an
overall, physically irrelevant, phase. Now, consider the physics of a subspace
of the Hilbert space where eigenstates of $\hat{G}$ with the same eigenvalue.
Since the operator $\hat{G}$ commutes with $\hat{H}_{W}$, the time evolution
of an initial state in this subspace will remain in the same subspace.
Together, these facts imply that the physics of this subspace is invariant
under spatial gauge transformations.

In the traditional quantization procedure (e.g., [6]), the physical states
$|\Psi_{EM}\rangle$ are taken to be eigenstates of $\hat{G}$ with zero
eigenvalue. In this subspace, Gauss’s law is preserved. But, this is a choice.
In a subspace of eigenstates of $\hat{G}$ with a non-zero eigenvalue
(function) $J_{s}^{0}\left({\bf x}\right)$,

$\hat{G}({\bf x})|\Psi_{EM}\rangle=J_{s}^{0}\left({\bf
x}\right)|\Psi_{EM}\rangle$ (20)

the states all transform with the same phase:

$e^{-i\int d^{3}{\bf x}\,\alpha\,\hat{G}}\,|\Psi_{EM}\rangle=e^{-i\int
d^{3}{\bf x}\,\alpha\,J_{s}^{0}}|\Psi_{EM}\rangle$ (21)

which thus also leads to gauge invariant physics. For these states, Gauss’s
law is not obeyed. Instead, the time evolution is such that

$\frac{d\langle\Psi_{EM}|\,\hat{G}\,|\Psi_{EM}\rangle}{dt}=0$ (22)

which is exactly the form of (6).

We thus see that there are gauge invariant quantum states that violate Gauss’s
law and there are no difficulties in time evolving these states. The initial
quantum state of the universe could have been a state where Gauss’s law was
preserved, i.e., a state that was an eigenstate of $\hat{G}$ with eigenvalue
zero. But it could just as easily have been an eigenstate of $\hat{G}$ with a
non-zero eigenvalue $J_{s}^{0}\left({\bf x}\right)$. In this case, Gauss’s law
would be violated, and it is a matter for experiment to decide which of these
scenarios is realized in our universe. We are simply choosing a Hilbert space
in which the states have a background, static, longitudinal electric field.
While such a field chooses a rest frame (and generally breaks Poincare
invariance), the current state of dynamical degrees of freedom (e.g., photons
and electrons) also seem to choose a rest frame. The only difference is that
we have turned on something additional that evolves trivially.

It is possible that it is not necessary to fix the remaining spatial gauge
symmetry. It would remain a time-independent symmetry of the Hamiltonian and
thus should not generate, for example, a photon mass. Without requiring this
fix, it should be possible to write physical states which are not eigenstates
of $\hat{G}$. While it is not clear this is important physically, it may play
an important role in the non-Abelian version of our story, which we will
explore in the next paper.

##  4 Other Gauges

In the Weyl gauge, we have shown the existence of gauge invariant quantum
states that violate Gauss’s law. In this section, we show how they can be
described in other gauges, such as the Coulomb gauge.

To describe this construction, we begin by reviewing the conventional methods
used to translate the physics of electromagnetism in the Weyl gauge to other
gauge choices. Since this is a review, we simply sketch the major steps and
refer the reader to [6] for the details. The first step in this procedure is
to construct the Lagrangian that describes the physics of the Weyl gauge from
the Hamiltonian $\hat{H}_{W}$ in (10). As described in Section 3, we restrict
the Hilbert space - to states that are eigenstates of
$\hat{G}={\bf\nabla}\cdot\hat{{\bf E}}-\hat{J}^{0}$ with the same eigenvalue
making the theory trivial under spatial gauge transformations. We can
construct the path integral (17) while restricting to such gauge-invariant
states. In the conventional quantization procedure, this is enforced by
inserting a projection operator $\hat{P}$ that projects the basis states
$|\hat{{\bf A}}\rangle$ onto the eigenspace of zero eigenvalue of $\hat{G}$.
Thus, the generating functional of the theory is:

$\mathcal{Z}=\text{tr}\left(Te^{-i\int dt\hat{H}_{W}}\hat{P}\right)$ (23)

where the projection operator

$\hat{P}=\Pi_{t,{\bf x}}\,\delta\left(\hat{G}\left({\bf x}\right)\right)$ (24)

This projection operator acts on each point in time and it can be implemented
in the path integral using an integral representation of the delta function
(24)

$\delta\left({\bf\nabla}\cdot{\bf E}\left({\bf x},t\right)-{J}^{0}\left({\bf
x},t\right)\right)=\int DA_{0}\,e^{i\delta t\int d^{3}{\bf x}\,A_{0}\left({\bf
x},t\right)\left({\bf\nabla}\cdot{\bf E}\left({\bf x},t\right)-J^{0}\left({\bf
x},t\right)\right)}$ (25)

We thus see the role of $A_{0}$ \- it enforces the constraint that the Hilbert
space of the theory is restricted to a specific space of states that are
annihilated by the operator $\hat{G}$. Constructing a path integral via the
standard procedure of inserting complete sets of states, and including the
above delta function, one obtains the action:

$S\left(A_{\mu},E\right)=\int d^{4}x\left(-{\bf E}\cdot\partial_{t}{\bf
A}-\frac{1}{2}\left({\bf E}^{2}+{\bf B}^{2}\right)-{\bf A}\cdot{\bf
J}+A_{0}\left({\bf\nabla}\cdot{\bf E}-J^{0}\right)+\cdots\right)$ (26)

where the ellipsis represents the remaining matter terms. From this, one
performs the Gaussian path integral over $E$ and other conjugate variables and
obtains the conventional Lagrangian (1).

Let us see how this procedure would change if we had picked a subspace of
states where $\hat{G}$ has a non-zero eigenvalue $J_{s}^{0}\left({\bf
x}\right)$. In this case, the projection operator would be:

$\hat{P}=\Pi_{t,{\bf x}}\,\delta\left(\hat{G}\left({\bf
x}\right)-J_{s}^{0}\left({\bf x}\right)\right)$ (27)

The subsequent mathematical procedure (following the steps in [6]) results in
the effective Lagrangian:

$\mathcal{\tilde{L}_{EM}}=-\frac{1}{4}F_{\mu\nu}F^{\mu\nu}+A_{\mu}J^{\mu}+A_{\mu}J_{s}^{\mu}+\mathcal{L}_{J}$
(28)

where the field $J_{s}^{\mu}$ is the background classical field
$J_{s}^{\mu}=\left(J_{s}^{0}\left({\bf x}\right),0,0,0\right)$. Thus, the
theory is identical to that of electromagnetism coupled to a classical
background charge density $J_{s}^{0}\left({\bf x}\right)$ \- this charge
density picks a rest frame and it is thus the theory of electromagnetism in
such a Lorentz breaking background. While it is well recognized that such a
classical background can be added to quantum electrodynamics, the key point of
our paper is to point out that there is no additional microphysics associated
with this background. The field theory that describes $J_{s}^{0}\left({\bf
x}\right)$ is simply a quantum state of electromagnetism.

The quantum theory of states with $J_{s}^{0}\left({\bf x}\right)\neq 0$ in any
gauge can now be described by applying the conventional quantization procedure
specific to that gauge starting with the effective Lagrangian (28). For
example, to obtain the quantum theory in Coulomb gauge, construct the
canonical Hamiltonian corresponding to the effective Lagrangian (28). Then,
impose the operator requirement ${\bf\nabla}\cdot\hat{{\bf A}}=0$ and solve
for $\hat{A}_{0}$ in terms of the effective charge density
$\hat{J}^{0}+J_{s}^{0}$:

$\hat{A}_{0}\left({\bf x}\right)=\int d^{3}{\bf
x^{\prime}}\frac{\left(\hat{J}^{0}\left({\bf
x^{\prime}}\right)+J_{s}^{0}\left({\bf
x^{\prime}}\right)\right)}{4\pi{\bf|x-x^{\prime}|}}$ (29)

The resulting Hamiltonian describes the quantum theory of states with
$J_{s}^{0}\left({\bf x}\right)\neq 0$.

##  5 Gravitation and Cosmology

In this section, we describe how states that violate Gauss’s law couple to
gravity. Our treatment of the gravitational interactions parallels the flat
space treatment - we begin by describing these states in the Weyl gauge where
the existence of these states is most easily understood. From the Weyl gauge,
we construct the effective Lagrangian that describes these states and the
gravitational dynamics can be readily read off from this Lagrangian.

Our principal interest here is to understand the cosmological implications of
such states. We will thus specialize to the case of a FRW cosmology - but, the
methods we describe can be extended to any space-time. Accordingly, we take
the metric of the space-time from the interval:

$ds^{2}=-N(t)^{2}dt^{2}+a(t)^{2}\left(dx^{2}+dy^{2}+dz^{2}\right)$ (30)

where $a(t)$ is the scale factor. $N(t)$ is a gravitational gauge degree of
freedom that needs to be fixed in order to define a Hamiltonian [1, 2]. While
we will eventually fix $N(t)$, we retain it for now for the sake of clarity.
Since we know how regular charged particles behave in the presence of these
states, we will focus our attention solely on the cosmological evolution of
these states, neglecting the physics of charged matter that may be coupled to
these states. The Lagrangian that describes this system is:

$\mathcal{L_{T}}=\sqrt{-g}\left(-\frac{1}{4}F_{\mu\nu}F^{\mu\nu}+\frac{M_{pl}^{2}}{2}R\right)$
(31)

With this Lagrangian, we obtain the conjugate momenta:

$\Pi^{j}=\frac{\partial\mathcal{L_{T}}}{\partial\dot{A_{j}}}$ (32)

for the electromagnetic degrees of freedom and

$\Pi_{a}=\frac{\partial\mathcal{L_{T}}}{\partial\dot{a}}$ (33)

for the scale factor. We pick the Weyl gauge by setting $A_{0}=0$. We then set
equal time commutation relations:

$[\hat{a},\hat{\Pi}_{a}]=i$

$[\hat{A}_{j}\left({\bf x}\right),\hat{\Pi}^{j^{\prime}}\left({\bf
x^{\prime}}\right)]=i\delta_{j}^{j^{\prime}}\delta\left({\bf
x-x^{\prime}}\right)$ (34)

With these definitions, the electric field operator is:

$\hat{E}_{i}=\hat{F_{0i}}=N\left(t\right)\frac{\hat{\Pi}_{i}}{\hat{a}}$ (35)

and the Hamiltonian for the full system (gravity and electromagnetism) is333A
discussion of the subtleties of defining operators with inverse scale factors
are discussed in [1]:

$\hat{H}_{T}=N\left(t\right)\left(\int d^{3}{\bf
x}\left(\frac{\hat{\bf{\Pi}}^{2}+\hat{\bf{B}}^{2}}{2\,\hat{a}}\right)-\frac{1}{12\,M_{pl}^{2}}\frac{\hat{\Pi}_{a}^{2}}{\hat{a}}\right)$
(36)

where we define ${\bf B}$ using the Maxwell tensor $F_{\mu\nu}$. The time
evolution of physical states $|\Psi\rangle$ is given by the equation:

$i\,\frac{\partial|\Psi\rangle}{\partial t}=\hat{H}_{T}|\Psi\rangle$ (37)

Similar to the situation in flat space, notice that
${\bf\nabla}\cdot{\hat{\bf\Pi}}$ commutes with the Hamiltonian $\hat{H}_{T}$.
Further, due to the commutation relations (34), ${\bf\nabla}\cdot\hat{\bf\Pi}$
is also the generator of spatial gauge transformations, where ${\bf\nabla}$ is
the gradient operator on the comoving coordinates ${\bf x}$. Following the
arguments in Section 3, to obtain gauge invariant physics, we demand that the
physical states $|\Psi\rangle$ are eigenstates of the operator
${\bf\nabla}\cdot\hat{\bf{\Pi}}$, i.e.:

${\bf{\nabla}}\cdot\hat{\bf{\Pi}}|\Psi\rangle=J_{s}^{0}\left({\bf
x}\right)|\Psi\rangle$ (38)

The above equations fully describe how these states interact with gravitation
for any choice of time parameterization $N(t)$. For simplicity, we will choose
$N(t)=1$.

Now we can examine what the states that violate Gauss’s law look like in an
expanding background. To understand this physics, we want to parallel the
discussion in Section 4 and obtain the effective Lagrangian that describes
these states in a covariant manner. For simplicity, we assume that the energy
density in these states is small and that the dynamical evolution of the scale
factor $a(t)$ is governed by other, more dominant energy densities. We are
thus interested in understanding how these states respond to a pre-determined
dynamical evolution of $a(t)$, i.e., we want to understand how these states
“redshift” in the semi-classical limit. With these assumptions, the physical
state $|\Psi\rangle=|a_{cl}\rangle\otimes|\Psi_{EM}\rangle$ where
$|a_{cl}\rangle$ describes the state of the metric and $|\Psi_{EM}\rangle$ is
the state of the electromagnetic fields. $|a_{cl}\rangle$ is a coherent state
of the gravitational field and in the semi-classical limit that we are in, we
can simply replace the operators $\hat{a}$ and $\hat{\Pi}_{a}$ in (36) with
their corresponding classical field values $a(t)$ and $\Pi_{a}(t)$.

We can now easily adapt our discussion to the discussion in Section 4.
Accordingly, for the generating functional:

$\mathcal{Z}=\text{tr}\left(Te^{-i\int dt\hat{H}_{T}}\hat{P}\right)$ (39)

with the projection operator:

$\hat{P}=\Pi_{t,{\bf x}}\,\delta\left({\bf\nabla}\cdot\hat{\bf{\Pi}}({\bf
x})-J_{s}^{0}({\bf x})\right)$ (40)

Enforcing this projection operator via the integral representation of the
delta function and proceeding as in Section 4, we get the effective
Lagrangian:

$\mathcal{\tilde{L}_{EM}}=\sqrt{-g}\left(-\frac{1}{4}F_{\mu\nu}F^{\mu\nu}+A_{\mu}{\mathbb{J}}^{\mu}\right)$
(41)

where the current ${\mathbb{J}}^{\mu}$ is of the form:

${\mathbb{J}}^{\mu}=\left(\frac{J_{s}^{0}\left({\bf
x}\right)}{a(t)^{3}},0,0,0\right)$ (42)

Notice that the current density ${\mathbb{J}}^{\mu}({\bf x})$ obeys the
conservation equation:

$\nabla_{\mu}{\mathbb{J}}^{\mu}=0$ (43)

where $\nabla_{\mu}$ is the covariant derivative with respect to the
background metric.

The physical effects of these states can be computed using the effective
Lagrangian (41) in concert with the conservation equation (43) - these effects
are identical to that of a conserved, static background charge density that
redshifts with the expansion of the universe. For example, breaking
$\hat{\bf\Pi}$ into transverse and longitudinal parts,
$\hat{\bf\Pi}=\hat{\bf\Pi}_{\perp}+{\bf\nabla}\hat{\Phi}$, the expectation
value of the Hamiltonian contains a contribution from the longitudinal part:

$\langle\Psi|\hat{H}_{T}|\Psi\rangle\supset-\frac{1}{2a(t)}\int d^{3}{\bf
x}\,d^{3}{\bf x^{\prime}}\,\frac{J_{s}^{0}({\bf x^{\prime}})J_{s}^{0}({\bf
x})}{4\pi|{\bf x}-{\bf x^{\prime}}|}\equiv\frac{V_{coulomb}^{s}}{a(t)}$ (44)

and we see that the effective Coulomb potential energy of the shadow charge
density redshifts as $a^{-1}$, as it would for any two charges a fixed
coordinate distance apart.

One can also read off the redshift from the energy-momentum tensor. In the
Weyl gauge, the time-time component of the energy-momentum tensor is:

$\displaystyle T_{0}^{0}$ $\displaystyle=$ $\displaystyle F_{0\mu}F^{\mu
0}-\frac{1}{4}\delta_{0}^{0}F_{\alpha\beta}F^{\alpha\beta}$ $\displaystyle=$
$\displaystyle-\frac{1}{2}\frac{F_{0i}^{2}}{a(t)^{2}}=-\frac{{\bf\Pi}^{2}}{2a(t)^{4}}$

where in the second line we used the quantum relation between $\hat{\Pi^{i}}$
and $\hat{F}_{0i}$ and inserted the classical expectation values. For the
Gauss’s law violation ${\bf\nabla}\cdot{\bf\Pi}=J_{s}^{0}$, the $\bf\Pi$ is
time-independent and thus the energy density redshifts as radiation.

##  6 Discussion

In this paper, we have shown that the theory of quantum electrodynamics
permits a broader class of quantum states that can be time evolved in a gauge
invariant manner than traditionally considered. These states violate Gauss’s
law, leading to electric fields in the universe that act as though they were
sourced by an immovable, conserved background classical charge density. But
there is no additional micro-physics associated with this charge density -
these are simply allowed states of electromagnetism. These states pick a rest
frame and thus break Lorentz invariance.

Since there is nothing logically wrong with such states, it is a matter of
experiment and observation to see if our universe is in a quantum state where
Gauss’s law is preserved or violated. It is likely that there would be rich
phenomenology associated with these states since the effective “shadow” charge
density that is associated with these states does not represent new dynamics
or degrees of freedom. The latter is typically constrained by a variety of
stringent astrophysical and collider limits, limiting its ability to
significantly impact standard model particles. But, these are simply states of
electromagnetism and thus they can have a significant impact on standard model
particles without being subject to such constraints. It would be interesting
to develop the phenomenology of these states and identify their cosmological,
astrophysical and laboratory signatures. This is especially important since
the discovery of such states would rule out a period of cosmic inflation in
the past history of our universe.

For simplicity, in this paper, we also assumed that the rest frame picked by
these states was the same as the cosmic rest frame. It would be interesting to
see if novel phenomenology could arise if this is not the case. It might be
that the relative motion between the cosmological background and these states
would effectively source long range, coherent magnetic fields and these might
be of cosmological significance. In addition to cosmology, it would also be
interesting to develop the phenomenology of these states in the vicinity of
black hole horizons. Since these states act as a fixed background, it is
plausible that they source new physical divergences around such horizons. This
may provide new opportunities to discover such states and potentially lead to
additional insight into divergences caused by quantum mechanics [8, 9] in the
vicinity of horizons.

While we have focused on the existence of these broader class of quantum
states in electromagnetism (in this paper) and gravitation (in our companion
paper [1]), it is likely that such states also exist for non-Abelian gauge
theories. These theories also possess non-dynamical degrees of freedom and
their quantum description requires suitably fixing various gauge degrees of
freedom. It is thus possible these theories also allow for the existence of
such “shadow” non-Abelian charge densities in a cosmological rest frame. In
future work, we intend to develop the theoretical framework to analyze these
effects and extract the associated phenomenology.

## Acknowledgements

We thank Michael Peskin and Raman Sundrum for fruitful discussions. This work
was supported by the U.S. Department of Energy (DOE), Office of Science,
National Quantum Information Science Research Centers, Superconducting Quantum
Materials and Systems Center (SQMS) under Contract No. DE-AC02-07CH11359.
D.E.K. and S.R. are supported in part by the U.S. National Science Foundation
(NSF) under Grant No. PHY-1818899. S.R. is also supported by the DOE under a
QuantISED grant for MAGIS. The work of S.R. was also supported by the Simons
Investigator Award No. 827042. T.M. is supported by the World Premier
International Research Center Initiative (WPI) MEXT, Japan, and by JSPS
KAKENHI grants JP19H05810, JP20H01896, JP20H00153, and JP22K18712.

## References

  * [1] David E. Kaplan, Tom Melia, and Surjeet Rajendran.  The Classical Equations of Motion of Quantized Gauge Theories, Part I: General Relativity.  5 2023.
  * [2] Anne-Katherine Burns, David E. Kaplan, Tom Melia, and Surjeet Rajendran.  Time Evolution in Quantum Cosmology.  4 2022.
  * [3] Paul Adrien Maurice Dirac.  The theory of gravitation in hamiltonian form.  Proceedings of the Royal Society A, 246, Aug 1958.
  * [4] Bryce S. DeWitt.  Quantum theory of gravity. i. the canonical theory.  Phys. Rev., 160:1113–1148, Aug 1967.
  * [5] Roman Jackiw.  Topological investigations in quantum gauge theories.  In R. Jackiw S. B. Treiman, E. Witten and B. Zumino, editors, Current Algebra And Anomalies, page 537p, Singapore, Singapore, 1985. World Scientific.
  * [6] Eduardo Fradkin.  Quantum Field Theory: An Integrated Approach.  Princeton University Press, August 2019.
  * [7] Steven Weinberg.  Quantum Theory of Fields, volume 1.  Cambridge University Press, Cambridge, UK, 1995.
  * [8] David E. Kaplan and Surjeet Rajendran.  Firewalls in General Relativity.  Phys. Rev. D, 99(4):044033, 2019.
  * [9] Stefan Hollands, Robert M. Wald, and Jochen Zahn.  Quantum instability of the Cauchy horizon in Reissner–Nordström–deSitter spacetime.  Class. Quant. Grav., 37(11):115009, 2020.


# Reachability Analysis of Neural Feedback Loops

Michael Everett, Golnaz Habibi, Chuangchuang Sun, and Jonathan P. How Authors
are with the Aerospace Controls Laboratory, Massachusetts Institute of
Technology, Cambridge, MA, 02139 USA e-mail<EMAIL_ADDRESS>version of
paper appeared at IEEE International Conference on Robotics and Automation
(ICRA) on June 1, 2021 [1]. Code:https://github.com/mit-
acl/nn˙robustness˙analysis

###### Abstract

Neural Networks (NNs) can provide major empirical performance improvements for
closed-loop systems, but they also introduce challenges in formally analyzing
those systems’ safety properties. In particular, this work focuses on
estimating the forward reachable set of neural feedback loops (closed-loop
systems with NN controllers). Recent work provides bounds on these reachable
sets, but the computationally tractable approaches yield overly conservative
bounds (thus cannot be used to verify useful properties), and the methods that
yield tighter bounds are too intensive for online computation. This work
bridges the gap by formulating a convex optimization problem for the
reachability analysis of closed-loop systems with NN controllers. While the
solutions are less tight than previous (semidefinite program-based) methods,
they are substantially faster to compute, and some of those computational time
savings can be used to refine the bounds through new input set partitioning
techniques, which is shown to dramatically reduce the tightness gap. The new
framework is developed for systems with uncertainty (e.g., measurement and
process noise) and nonlinearities (e.g., polynomial dynamics), and thus is
shown to be applicable to real-world systems. To inform the design of an
initial state set when only the target state set is known/specified, a novel
algorithm for backward reachability analysis is also provided, which computes
the set of states that are guaranteed to lead to the target set. The numerical
experiments show that our approach (based on linear relaxations and
partitioning) gives a $5\times$ reduction in conservatism in $150\times$ less
computation time compared to the state-of-the-art. Furthermore, experiments on
quadrotor, 270-state, and polynomial systems demonstrate the method’s ability
to handle uncertainty sources, high dimensionality, and nonlinear dynamics,
respectively.

###### Index Terms:

Reachability Analysis, Neural Networks, Deep Learning, Safety, Verification,
Safe Learning

##  I Introduction

Neural Networks (NNs) are pervasive in many fields because of their ability to
express highly general input-output relationships, such as for perception,
planning, and control tasks in robotics. However, before deploying NNs on
safety-critical systems, there must be techniques to guarantee that the
closed-loop behavior of systems with NNs, which we call neural feedback loops
(NFLs), will meet desired specifications. The goal of this paper is to develop
a framework for guaranteeing that systems with NN controllers will reach their
goal states while avoiding undesirable regions of the state space, as in Fig.
1.

Despite the importance of analyzing closed-loop behavior, much of the recent
work on formal NN analysis has focused on NNs in isolation (e.g., for image
classification) [2, 3, 4, 5, 6, 7], with an emphasis on efficiently relaxing
NN nonlinearities [8, 9, 10, 11, 12, 13, 14]. On the other hand, closed-loop
system reachability has been studied for decades, but traditional approaches,
such as Hamilton-Jacobi methods [15, 16], do not consider NNs in the loop.

Several recent works [17, 18, 19, 20, 21, 22] propose methods that compute
forward reachable sets of NFLs. A key challenge is in maintaining
computational efficiency while still providing tight bounds on the reachable
sets. In addition, large initial state sets present challenges as the
underlying NN relaxations become loose, leading to loose reachable set
estimates. The literature also typically assumes perfect knowledge of system
dynamics, with no stochasticity. Furthermore, the literature focuses on
forward reachability analysis (what states could the system end in, starting
from a given $\mathcal{X}_{0}$?), but backward reachability analysis (what
states can the system start from, to end in a given $\mathcal{X}_{T}$?) is an
equally important problem in safety analysis. While forward and backward
reachability differ only by a change of variables in classical systems [16,
23], both the nonlinearities and matrix dimensions of NN controllers raise
fundamental challenges that prevent propagating sets backward through a NFL.

Figure 1: Forward Reachability Analysis. The objective is to compute the blue
sets $\mathcal{R}_{t}(\mathcal{X}_{0})$, to ensure a system starting from
$\mathcal{X}_{0}$ (yellow) ends in $\mathcal{G}$ (green) and avoids
$\mathcal{A}_{0},\mathcal{A}_{1}$ (red). This is especially challenging for
systems with NN control policies.

To address these challenges, this work’s contributions include:

  * •

A linear programming-based formulation of forward reachability analysis for
NFLs, which provides a computationally efficient method for verifying safety
properties,

  * •

The use of input set partitioning techniques to provide tight bounds on the
reachable sets despite large initial state sets,

  * •

The consideration of measurement noise, process noise, and nonlinear dynamics,
which improves the applicability to real systems with uncertainty, and

  * •

A framework for backward reachability analysis of NFLs, which enables
estimating which starting states will lead to a target set despite non-
invertible NN weight matrices and activations.

Numerical experiments show that the new method provides $5\times$ better
accuracy in $150\times$ less computation time compared to [22] and that the
method can handle uncertainty sources, high dimensionality, and nonlinear
dynamics via applications on quadrotor, 270-state, and polynomial systems.

This article extends our conference paper [1] by providing an improved problem
formulation, notation, and algorithmic details. It also includes:

  * •

A closed-form solution to provide an order-of-magnitude computational speedup
for the same bounds in Section IV-E,

  * •

An extension of state-of-the-art partitioning techniques to the NFL setting in
Section V,

  * •

Formulations that handle nonlinearities (control limits and polynomial terms)
in the dynamics in Section VI,

  * •

A framework for backward reachability analysis in Section VII.

##  II Related Work

Related work on reachability analysis can be categorized into works on NNs in
isolation, closed-loop systems without NNs, and closed-loop systems with NNs.
For instance, machine learning literature includes many methods to verify
properties of NNs, often motivated by defending against adversarial examples
[24]. These methods broadly range from exact [25] to tight [14] to efficient
[11] to fast [8]. Although these tools are not designed for closed-loop
systems, the NN relaxations from [11] provide a key foundation here.

For closed-loop systems, reachability analysis is a standard component of
safety verification. Modern methods include Hamilton-Jacobi Reachability
methods [15, 16], SpaceEx [26], Flow* [27], CORA [28], and C2E2 [29, 30], but
these do not account for NN control policies. Orthogonal approaches that do
not explicitly estimate the system’s forward reachable set, but provide other
notions of safety, include Lyapunov function search [31] and control barrier
functions (CBFs) [32].

Recent reachability analysis approaches that do account for NN control
policies face a tradeoff between computation time and conservatism. [17, 18,
19] use polynomial approximations of NNs to make the analysis tractable. Most
works consider NNs with ReLU approximations, whereas [20] considers sigmoidal
activations. [21, 33] introduce conservatism by assuming the NN controller
could output its extreme values at every state. Most recently, [22] formulated
the problem as a SDP, called Reach-SDP. This work builds on both [21, 22] and
makes the latter more scalable by re-formulating the SDP as a linear program,
introduces sources of uncertainty in the closed-loop dynamics, and shows
further improvements by partitioning the input set.

##  III Preliminaries

###  III-A Closed-Loop System Dynamics

Figure 2: Neural Feedback Loop. Many real world deployments of NNs are part of
complex, closed-loop systems (e.g., robots with learned policies take actions
that modify the world state and influence future observations).

Consider a discrete-time linear time-varying system,

$\displaystyle\begin{split}\mathbf{x}_{t+1}&=\mathbf{A}_{t}\mathbf{x}_{t}+\mathbf{B}_{t}\mathbf{u}_{t}+\mathbf{c}_{t}+\bm{\omega}_{t}\\\
\mathbf{y}_{t}&=\mathbf{C}_{t}^{T}\mathbf{x}_{t}+\bm{\nu}_{t},\end{split}$ (1)

where
$\mathbf{x}_{t}\in\mathds{R}^{n_{x}},\,\mathbf{u}_{t}\in\mathds{R}^{n_{u}},\,\mathbf{y}_{t}\in\mathds{R}^{n_{y}}$
are state, control, and output vectors,
$\mathbf{A}_{t},\,\mathbf{B}_{t},\,\mathbf{C}_{t}$ are known system matrices,
$\mathbf{c}_{t}\in\mathds{R}^{n_{x}}$ is a known exogenous input, and
$\bm{\omega}_{t}\sim\Omega$ and $\bm{\nu}_{t}\sim N$ are process and
measurement noises sampled at each timestep from unknown distributions with
known, finite support (i.e.,
$\bm{\omega}_{t}\in[\underaccent{\bar}{\bm{\omega}}_{t},\,\bar{\bm{\omega}}_{t}],\bm{\nu}_{t}\in[\underaccent{\bar}{\bm{\nu}}_{t},\,\bar{\bm{\nu}}_{t}]$
element-wise).

We assume an output-feedback controller $\mathbf{u}_{t}=\pi(\mathbf{y}_{t})$
parameterized by an $m$-layer feed-forward NN, optionally subject to control
constraints, $\mathbf{u}_{t}\in\mathcal{U}_{t}$. We denote the closed-loop
system with dynamics Eq. 1 and control policy $\pi$ as

$\mathbf{x}_{t+1}=f(\mathbf{x}_{t};\pi).$ (2)

This NFL is visualized in Fig. 2. At each timestep, the state $\mathbf{x}$
enters the perception block and is perturbed by sensor noise $\bm{\nu}$ to
create observation $\mathbf{y}$. The observation is fed into the learned
policy, which is a NN that selects the control input $\mathbf{u}$. The control
input enters the plant, which is defined by $\mathbf{A},\mathbf{B},\mathbf{c}$
matrices and subject to process noise $\bm{\omega}$. The output of the plant
is the state vector at the next timestep.

###  III-B Reachable Sets

For the closed-loop system Eq. 2, we denote $\mathcal{R}_{t}(\mathcal{X}_{0})$
the forward reachable set at time $t$ from a given set of initial conditions
$\mathcal{X}_{0}\subseteq\mathds{R}^{n_{x}}$, which is defined by the
recursion

$\mathcal{R}_{t+1}(\mathcal{X}_{0})=f(\mathcal{R}_{t}(\mathcal{X}_{0});\pi),\qquad\mathcal{R}_{0}(\mathcal{X}_{0})=\mathcal{X}_{0}.$
(3)

###  III-C Finite-Time Reach-Avoid Verification Problem

The finite-time reach-avoid property verification problem is defined as
follows: Given a goal set $\mathcal{G}\subseteq\mathds{R}^{n_{x}}$, a sequence
of avoid sets $\mathcal{A}_{t}\subseteq\mathds{R}^{n_{x}}$, and a sequence of
reachable set estimates $\mathcal{R}_{t}\subseteq\mathds{R}^{n_{x}}$,
determining that every state in the final estimated reachable set will be in
the goal set and any state in the estimated reachable sets will not enter an
avoid set requires computing set intersections,
$\texttt{VERIFIED}(\mathcal{G},\mathcal{A}_{0:N},\mathcal{R}_{0:N})\equiv{\mathcal{R}_{N}\subseteq\mathcal{G}}\
\mathrm{\&}\ {\mathcal{R}_{t}\cap\mathcal{A}_{t}=\emptyset},\forall
t\in\\{0,\ldots,N\\}$.

In the case of our nonlinear closed-loop system Eq. 2, where computing the
reachable sets exactly is computationally intractable, we can instead compute
outer-approximations of the reachable sets,
$\bar{\mathcal{R}}_{t}(\mathcal{X}_{0})\supseteq\mathcal{R}_{t}(\mathcal{X}_{0})$.
This is useful if the finite-time reach-avoid properties of the system as
described by outer-approximations of the reachable sets are verified, because
that implies the finite-time reach-avoid properties of the exact closed loop
system are verified as well. Tight outer-approximations of the reachable sets
are desirable, as they enable verification of tight goal and avoid set
specifications, and they reduce the chances of verification being unsuccessful
even if the exact system meets the specifications.

###  III-D Control Policy Neural Network Structure

Using notation from [11], for the $m$-layer neural network used in the control
policy, the number of neurons in each layer is $n_{l},\ \forall l\in[m]$,
where $[a]$ denotes the set $\\{1,2,\ldots,a\\}$. Let the $l$-th layer weight
matrix be $\mathbf{W}^{(l)}\in\mathds{R}^{n_{l}\times n_{l-1}}$ and bias
vector be $\mathbf{b}^{(l)}\in\mathds{R}^{n_{l}}$, and let
$\mathbf{F}_{l}:\mathds{R}^{n_{x}}\to\mathds{R}^{n_{l}}$ be the operator
mapping from network input (measured output vector $\mathbf{y}_{t}$) to layer
$l$. We have
$\mathbf{F}_{l}(\mathbf{y}_{t})=\sigma(\mathbf{W}^{(l)}\mathbf{F}_{l-1}(\mathbf{y}_{t})+\mathbf{b}^{(l)}),\forall
l\in[m-1]$, where $\sigma(\cdot)$ is the coordinate-wise activation function.
The framework applies to general activations, including ReLU,
$\sigma(\mathbf{z})=\mathrm{max}(0,\mathbf{z})$. The network input
$\mathbf{F}_{0}(\mathbf{y}_{t})=\mathbf{y}_{t}$ produces the (unclipped)
control input,

$\mathbf{u}_{t}=\pi(\mathbf{y}_{t})=\mathbf{F}_{m}(\mathbf{y}_{t})=\mathbf{W}^{(m)}\mathbf{F}_{m-1}(\mathbf{y}_{t})+\mathbf{b}^{(m)}.$
(4)

###  III-E Neural Network Robustness Verification

A key step in quickly computing reachable sets of the NFL Eq. 2 is to relax
the nonlinear constraints induced by the NN’s nonlinear activation functions.
The relaxation converts each nonlinearity into a linear upper and lower bound,
where each bound holds within the known range of inputs to the activation.

To represent the range of inputs to the activation functions, we first define
the $\epsilon$-ball. Denote the $\epsilon$-ball (also called the
$\ell_{p}$-ball) under the $\ell_{p}$ norm, centered at
$\mathring{\mathbf{x}}$, with scalar radius, $\epsilon$,

$\displaystyle\mathcal{B}_{p}(\mathring{\mathbf{x}},\epsilon)$
$\displaystyle=\\{\mathbf{x}\ \lvert\
\lvert\lvert\mathbf{x}-\mathring{\mathbf{x}}\rvert\rvert_{p}\leq\epsilon\\}.$
(5)

The $\epsilon$-ball is extended to the case of vector
$\bm{\epsilon}\in\mathds{R}^{n}_{\geq 0}$ (i.e., $\bm{\epsilon}$-ball),
defined as

$\displaystyle\mathcal{B}_{p}(\mathring{\mathbf{x}},\bm{\epsilon)}$
$\displaystyle=\\{\mathbf{x}\ \lvert\
\lim_{\bm{\epsilon}^{\prime}\to\bm{\epsilon}^{+}}\lvert\lvert(\mathbf{x}-\mathring{\mathbf{x}})\oslash\bm{\epsilon}^{\prime}\rvert\rvert_{p}\leq
1\\},$ (6)

where $\oslash$ denotes element-wise division.

######  Theorem III.1 (From [11], Convex Relaxation of NN).

Given an $m$-layer neural network control policy
$\pi:\mathds{R}^{n_{y}}\to\mathds{R}^{n_{u}}$, there exist two explicit
functions $\pi_{j}^{L}:\mathds{R}^{n_{y}}\to\mathds{R}^{n_{u}}$ and
$\pi_{j}^{U}:\mathds{R}^{n_{y}}\to\mathds{R}^{n_{u}}$ such that $\forall
j\in[n_{m}],\forall\mathbf{y}\in\mathcal{B}_{p}(\mathring{\mathbf{y}},\epsilon)$,
the inequality
$\pi_{j}^{L}(\mathbf{y})\leq\pi_{j}(\mathbf{y})\leq\pi_{j}^{U}(\mathbf{y})$
holds true, where

$\displaystyle\pi_{j}^{U}(\mathbf{y})$
$\displaystyle=\mathbf{\Psi}_{j,:}\mathbf{y}+\bm{\alpha}_{j}$ (7)
$\displaystyle\pi_{j}^{L}(\mathbf{y})$
$\displaystyle=\mathbf{\Phi}_{j,:}\mathbf{y}+\bm{\beta}_{j},$ (8)

where $\mathbf{\Psi},\mathbf{\Phi}\in\mathds{R}^{n_{u}\times n_{y}}$ and
$\bm{\alpha},\bm{\beta}\in\mathds{R}^{n_{u}}$ are defined recursively using NN
weights, biases, and activations (e.g., ReLU, sigmoid, tanh), as detailed in
[11].

In a closed-loop system, Theorem III.1 bounds the control output for a
particular measurement $\mathbf{y}$. Moreover, if all that is known is
$\mathbf{y}\in\mathcal{B}_{p}(\mathring{\mathbf{y}},\epsilon)$, Theorem III.1
provides affine relationships between $\mathbf{y}$ and $\mathbf{u}$ (i.e.,
bounds valid within the known set of possible $\mathbf{y}$). These
relationships enable efficient calculation of NN output bounds, using
Corollary 3.3 of [11].

We could leverage [11] to compute reachable sets by first bounding the
possible controls, then bounding the next state set by applying the extreme
controls from each state. This is roughly the approach in [21, 33], for
example. However, this introduces excessive conservatism, because both
extremes of control would not be applied at every state (barring pathological
examples). To produce tight bounds on the reachable sets, we leverage the
relationship between measured output and control in Section IV.

##  IV Forward Reachability Analysis

Figure 3: Approach Overview for simple 1D system. Theorem III.1 relaxes the NN
to give affine relationships between observation $\mathbf{y}_{t}$ and control:
$\pi^{U},\pi^{L}$. Lemma IV.1 uses the system dynamics to associate
$\pi^{U},\pi^{L}$ with the next state set. Lemma IV.2 optimizes the closed-
loop dynamics over all states $\mathbf{x}_{t}\in\mathcal{X}_{t}$ to compute
bounds on the next state, $\bm{\gamma}^{U}_{t+1},\bm{\gamma}^{L}_{t+1}$.

Recall that our goal is to find the set of all possible next states,
$\mathbf{x}_{t+1}\in\mathcal{X}_{t+1}$, given that the current state lies
within a known set, $\mathbf{x}_{t}\in\mathcal{X}_{t}$. This will allow us to
compute reachable sets recursively starting from an initial set
$\mathcal{X}_{t}=\mathcal{X}_{0}$.

The approach follows the architecture in Fig. 3. After first relaxing the NN
controller using Theorem III.1, we then associate linearized extreme
controllers with extreme next states in Section IV-B. Then, using the
linearized extreme controller, we optimize over all states in the input set to
find extreme next states in Section IV-C. We extend the formulation to handle
control limits in Section VI-A, then describe how to convert the solutions of
the optimization problems into reachable set descriptions in Section IV-D.

###  IV-A Assumptions & Problem Statement

This work assumes that $\mathcal{X}_{t}$ is described by either:

  * •


# Inverse design of dissipative quantum steady-states with implicit
differentiation

Rodrigo A. Vargas-Hernández 22footnotemark: 2 33footnotemark: 3
<EMAIL_ADDRESS>
&Ricky T. Q. Chen 11footnotemark: 1 33footnotemark: 3
<EMAIL_ADDRESS>Kenneth A. Jung 11footnotemark: 1 22footnotemark: 2
<EMAIL_ADDRESS>&Paul Brumer 11footnotemark: 1 22footnotemark: 2
<EMAIL_ADDRESS>University of Toronto. †Chemical Physics Theory Group,
Department of Chemistry. ‡Vector Institute.

###### Abstract

Inverse design of a property that depends on the steady-state of an open
quantum system is commonly done by grid-search type of methods. In this paper
we present a new methodology that allows us to compute the gradient of the
steady-state of an open quantum system with respect to any parameter of the
Hamiltonian using the implicit differentiation theorem. As an example, we
present a simulation of a spin-boson model where the steady-state solution is
obtained using Redfield theory.

##  1 Introduction

The field of open quantum systems (OQS) is focused on studying the interaction
of a _system_ (S) with its surroundings, typically referred to as the _bath_
(B) (Breuer and Petruccione, 2002). The Hamiltonian for the complete problem
is defined as the sum of the isolated system Hamiltonian, $H_{S}$, the bath
that surrounds the system, $H_{B}$, and the interaction between the system and
the bath, $H_{SB}$,

$H=H_{S}+H_{B}+H_{SB}.$ (1)

For physical problems the bath is often represented as a continuum that
describes the vibrational or solvent degrees of freedom. Given the large
dimensionality of the composite Hilbert space of such systems, it is common to
construct a reduced representation by tracing out the bath degrees of freedom,
$\rho_{S}=\text{Tr}_{B}[\rho]$. Here, we refer to $\rho_{S}$ as the reduced
density matrix (RDM) and $\rho$ is the density matrix of the complete system.

In OQS, one is usually interested in the dissipative effects the bath induces
in the system over time. However, for many OQS the properties of interest are
related to the steady-state (SS) of the system ($\rho^{ss}$),
$\;d\rho(t)/dt=0$; for example, energy transfer efficiency in biological
systems excited by natural incoherent light (Tscherbul and Brumer, 2014, 2015;
Pachón et al., 2017) and the performance of quantum heat engines/refrigerators
(Kilgour and Segal, 2018; Linden et al., 2010). Here, we propose a novel
numerical algorithm to compute the gradient of the steady-state of an OQS with
respect to any parameter of the Hamiltonian ($\bm{\theta}_{i}$) using implicit
differentiation, $\partial\rho^{ss}/\partial\bm{\theta}_{i}$, allowing us to
tackle problems related to inverse design and sensitivity analysis.

In the following sections, we discuss how as an example, we differentiate the
steady-state of the spin-boson (SB) model with respect to any parameter using
an ordinary differential equation solver and implicit differentiation.

##  2 Differentiation of the steady-state

Conditioned on initial values $\rho(t_{0})\in\mathbb{C}^{d}$ and free
parameters $\bm{\theta}$, let $\rho(t)\in\mathbb{C}^{d}$ be the solution to a
homogeneous dynamical system parameterized by an ordinary differential
equation (ODE) $\frac{d\rho}{dt}=f(\rho,\bm{\theta})$ where $d$ is the
dimensionality of the RDM. In optimization, we want to compute the gradient of
the steady state with respect to the parameters $\bm{\theta}$. Let us denote
$\rho^{ss}$ as the steady-state, which satisfies $f(\rho^{ss},\bm{\theta})=0$.

We can solve for $\rho^{ss}$ by running an appropriate ODE solver for a
sufficiently long period of time. However, differentiating through the
internals of the ODE solver is computationally expensive as it requires
storing all intermediate quantities of the solver. The adjoint method for
computing gradients of ODE solutions requires either the trajectory $\rho(t)$
to be stored in memory or solving $\rho(t)$ in reverse time, which uses
constant memory (Chen et al., 2018). However, since we are interested in
$\rho^{ss}$, the reversing approach is not applicable as the steady state,
once reached, cannot be reversed. Instead, since $\rho^{ss}$ is the solution
of a fixed point problem, the Jacobian $\frac{d\rho^{ss}}{d\bm{\theta}}$ can
be expressed using the implicit function theorem (Krantz and Parks, 2012),

$\frac{d\rho^{ss}}{d\bm{\theta}}=-\left(\frac{df(\rho^{ss},\bm{\theta})}{d\rho}\right)^{-1}\left[\frac{df(\rho^{ss},\bm{\theta})}{d\bm{\theta}}\right].$
(2)

We show that it is possible to compute this gradient, i.e. Eq. (2), using
constant memory cost, and without the need to know how $\rho^{ss}$ is computed
as long as it satisfies the steady state criterion, i.e.
$f(\rho^{ss},\bm{\theta})=0$.

For this work, we would also have a scalar loss function ${\cal
L}(\bm{\theta},\rho^{ss})$ that we wish to minimize. Here, the gradient with
respect to parameters can be factored with the chain rule $\frac{d{\cal
L}}{d\rho^{ss}}\frac{d\rho^{ss}}{d\bm{\theta}}$. In general, we’re interested
in vector-Jacobian products of the form $v\frac{d\rho^{ss}}{d\bm{\theta}}$
where $v$ is any vector. In automatic differentiation libraries such as jax
(Bradbury et al., 2018), we can easily compute vector-Jacobian products for
large systems but recovering the full Jacobian is more costly.

Computing a vector-Jacobian product using, Eq. (2), requires solving $vJ^{-1}$
where $J=\left(\frac{df(\rho^{ss},\bm{\theta})}{d\rho}\right)^{-1}$. For large
systems where $J$ cannot be tractably computed, we compute $vJ^{-1}$ by
recognizing that it is the steady-state solution of the ODE
$\frac{dy}{dt}=yJ-v$. Notably, simulating this ODE only requires vector-
Jacobian products, which are inexpensive in automatic differentiation
libraries. In the following sections we illustrate how inverse design and
sensitivity analysis are possible for dissipative quantum systems by
differentiating through the steady state of a Spin-Boson model.

##  3 Example: Spin-Boson model

In the field of OQS, the SB Hamiltonian (Thoss et al., 2001) is a standard
model used to describe a wide variety of physical phenomena (Breuer and
Petruccione, 2002; Leggett et al., 1987) including electron transfer (Egorova
et al., 2003), heat transport (Xu and Cao, 2016), and energy transfer (Liu and
Segal, 2019; Jung and Brumer, 2020). The total Hamiltonian is,

$\displaystyle
H_{S}=\frac{\varepsilon}{2}\sigma_{z}+\frac{\Delta}{2}\sigma_{x},\;\;\;\;H_{B}=\sum_{k}\omega_{k}b_{k}^{\dagger}b_{k},\;\;\;\;H_{SB}=\sigma_{z}\sum_{k}\lambda_{k}(b_{k}^{\dagger}+b_{k}),$
(3)

where $b_{k}^{\dagger}$ ($b_{k}$) is the creation (annihilation) operator of
mode $k$ in the bath, $\sigma_{z}$ and $\sigma_{x}$ are Pauli matrices, and
$\\{\lambda_{k}\\}$ are the coupling strength parameters, and $\varepsilon$
and $\Delta$ are system parameters.

The exact equation of motion (Nakajima, 1958; Zwanzig, 1960) for the RDM is
not practical to solve since one requires knowledge of the interacting system-
bath dynamics. Redfield theory (RT) (Redfield, 1965; Palenberg et al., 2001)
is a practical alternative for regimes in which the system-bath coupling is
weak, i.e. $\\{\lambda_{k}\\}\ll 1$. The RT equations of motions written in
the system eigenbasis are,

$\displaystyle\frac{\partial\rho_{\mu,\nu}(t)}{\partial
t}=-i\omega_{\mu,\nu}\rho_{\mu,\nu}(t)+\sum_{\kappa,\lambda}R_{\mu,\nu,\kappa,\lambda}\;\rho_{\kappa,\lambda}(t),$
(4)

where $[\mu,\nu,\kappa,\lambda]$ are the index of the eigenstates of $H_{S}$,
i.e. they satisfy $H_{S}|\mu\ \rangle=\varepsilon_{\mu}|\mu\rangle$, and
$\omega_{\mu,\nu}$ is the difference between eigenvalues $\varepsilon_{\mu}$
and $\varepsilon_{\nu}$. $R_{\mu,\nu,\kappa,\lambda}$ are the Redfield tensors
which describe the interaction of the system and bath and are given by,

$\displaystyle
R_{\mu,\nu,\kappa,\lambda}=\Gamma^{+}_{\lambda,\nu,\mu,\kappa}+\Gamma^{-}_{\lambda,\nu,\mu,\kappa}-\delta_{\nu,\lambda}\sum_{\alpha}\Gamma^{+}_{\mu,\alpha,\alpha,\kappa}-\delta_{\mu,\kappa}\sum_{\alpha}\Gamma^{+}_{\lambda,\alpha,\alpha,\nu},$
(5)

which contain the transition rates,

$\displaystyle\Gamma^{+}_{\lambda,\nu,\mu,\kappa}$ $\displaystyle=$
$\displaystyle\langle\lambda|\sigma_{z}|\nu\rangle\langle\mu|\sigma_{z}|\kappa\rangle\int_{0}^{\infty}\mathrm{d}\tau
F(\tau)e^{-i\omega_{\mu,\kappa}\tau},$ (6)
$\displaystyle\Gamma^{-}_{\lambda,\nu,\mu,\kappa}$ $\displaystyle=$
$\displaystyle\langle\lambda|\sigma_{z}|\nu\rangle\langle\mu|\sigma_{z}|\kappa\rangle\int_{0}^{\infty}\mathrm{d}\tau
F^{*}(\tau)e^{-i\omega_{\lambda,\nu}\tau},$ (7)

that are in turn comprised of the bath correlation function,

$\displaystyle
F(\tau)=\int_{0}^{\infty}\mathrm{d}\omega\;g(\omega)\left[\coth(\beta\omega/2)\cos(\omega\tau)-i\sin(\omega\tau)\right].$
(8)

$g(\omega)$ is the spectral density function,
$g(\omega)=\eta\omega^{s}\omega_{c}^{1-s}e^{-\omega/\omega_{c}}$; for this
work we used the super Ohmic form of $s=3$. $\eta$ is the bath friction
parameter that is on the order of $\lambda_{k}^{2}$ and $\beta$ is the inverse
temperature.

##  4 Results

###  4.1 Sensitivity analysis

One of the most studied observables for the SB model is the population
difference at equilibrium,
$\langle\sigma_{z}\rangle=\text{Tr}[\sigma_{z}\rho^{ss}]$. Fig. (1)
illustrates the change of $\langle\sigma_{z}\rangle$ with respect to bath
parameters $\beta$, $\eta$ and the system parameters $\Delta$. For all
calculations we used $w_{c}=1$, $\varepsilon=0.1$ and the initial density
matrix was taken to be

$\displaystyle\rho_{S}(t_{0})=\begin{pmatrix}\frac{3}{4}&-i\frac{\sqrt{3}}{4}\\\
i\frac{\sqrt{3}}{4}&\frac{1}{4}\end{pmatrix}.$ (9)

We can observe that $\langle\sigma_{z}\rangle$ is independent of $\beta$ and
$\eta$, as the gradients are effectively zero. However,
$\frac{\partial\langle\sigma_{z}\rangle}{\partial\Delta}$ is nonzero which
confirms the effect $\Delta$ has in $\langle\sigma_{z}\rangle$.

While finite differences with respect to $\beta$, $\eta$, and $\Delta$ can be
computed one at a time, this requires $\mathcal{O}(d)$ evaluations of the
steady state. The automatic differentiation process discussed in Section 2
allows us to compute analytical derivatives for all parameters simultaneously
with just one evaluation of the steady-state.

To construct the Redfield tensors we must know the eigenvalues and eigenstates
of $H_{S}$ beforehand. As $H_{S}$ is not diagonal for $\Delta>0$, the
eigenvalues and eigenstates do not have trivial analytical solutions. To
compute these, we used jax’s eigendecomposition. This ensures this method can
be applied to arbitrary systems.

Figure 1: $\langle\sigma_{z}\rangle=\text{Tr}[\sigma_{z}\rho^{ss}]$ as a
function of different parameters, $\beta$, $\eta$ and $\Delta$. The inset of
each figure compares the gradient computed with Eq. (2) (blue solid curve) and
finite differences (orange solid curve). For all calculations, except for the
parameter in play, we used $\beta=0.1$, $\eta=0.01$, $w_{c}=1$,
$\varepsilon=0.1$, and $\Delta=0$. The initial density matrix used was
$\rho_{S}(t_{0})=[3/4,-i\frac{\sqrt{3}}{4},i\frac{\sqrt{3}}{4},1/4]$ (Eq. 9).

###  4.2 Inverse design for system’s Hamiltonian

Gradient based methods have proven to be powerful numerical tools to find the
minimizer of loss functions. Given the possibility to compute quantities like
$\partial{\cal L}/\partial\theta_{i}$, we could reformulate the inverse design
problem into an optimization one. For example, what values of $\varepsilon$
and $\Delta$ reproduce a target observable. In order to do so we define an
error function, e.g., ${\cal
L}(\varepsilon,\Delta)=\left\|\langle\sigma_{z}\rangle-\hat{\langle\sigma_{z}\rangle}\right\|_{2}$.

Fig. 2 illustrates how gradient-based methods can be used to search for the
optimal values of $\varepsilon$ and $\Delta$ that reproduce a target
$\langle\sigma_{z}\rangle$; as a proof of principle calculation we used
$\varepsilon=0.1$, $\Delta=0$ where $\langle\sigma_{z}\rangle=0.04995847$. For
these simulations we fixed the values of $\beta$ and $\eta$. Furthermore, in
Table 1 we report the values of $\varepsilon$ and $\Delta$ for 5 different
optimizations. To minimize ${\cal L}$ we used the Adam algorithm, which is a
first-order gradient-based optimization algorithm (Kingma and Ba, 2017).

Table 1:  $\varepsilon$ | $\Delta$ | $\langle\sigma_{z}\rangle$ | ${\cal L}$
---|---|---|---
0.0835 | 0.0578 | 0.049947 | 1.2E-5
0.0867 | 0.0526 | 0.049942 | 1.6E-5
0.0739 | 0.0692 | 0.049913 | 4.5E-5
0.0490 | 0.0809 | 0.049921 | 3.8E-5
0.0822 | 0.0597 | 0.049940 | 1.9E-5

(a)

(b)

Figure 2: We optimize the parameters of $H_{S}$, $\varepsilon$ and $\Delta$,
with respect to a target $\hat{\langle\sigma_{z}\rangle}$ computed with
$\varepsilon=0.1$ and $\Delta=0$. We used the Adam optimization algorithm to
minimize ${\cal
L}(\varepsilon,\Delta)=\left\|\langle\sigma_{z}\rangle-\hat{\langle\sigma_{z}\rangle}\right\|_{2}$.
(a) The value of $\langle\sigma_{z}\rangle$ at each iteration of the
optimization routine for different random initialization. The black-dashed
curves is the target $\langle\sigma_{z}\rangle$. (b) The iterations of the
optimization procedure in the parameter space for different random
initialization. The white-$x$ symbol represents the target set of parameters
and the white-$\bullet$ is the set of parameters found at 100th iteration of
the Adam algorithm. $\frac{\partial{\cal L}}{\partial\varepsilon}$ and
$\frac{\partial{\cal L}}{\partial\Delta}$ used in the Adam optimizer were
computed with Eq. (2). For all calculations, we used $\beta=0.1$, $\eta=0.01$,
$w_{c}=1$ and $\rho(t_{0})$ was the same as in Fig. (1). The learning rate was
set to 0.1, and we parameterized $\varepsilon$ and $\Delta$ to be strictly
positive by using the softplus function.

##  5 Summary

We have presented a novel numerical methodology capable of computing the
gradient of quantum observables of the form
$\langle\hat{O}\rangle=\text{Tr}[\hat{O}\rho^{ss}]$, with respect to any
parameter of the Hamiltonian. As we stated, this procedure does not depend on
the numerical procedure used to solve for $\rho^{ss}$ and it uses constant
memory. By computing the gradient of $\langle\sigma_{z}\rangle$, for the
particular parameter set and model used we found an independence of
$\langle\sigma_{z}\rangle$ with respect to the bath parameters.
$\langle\sigma_{z}\rangle$ increases as a function of $\Delta$ and the
gradient computed with Eq. (2) matches the gradient computed with finite
differences.

Inverse design for physical systems is one of the most common applications for
gradient-based search algorithms. We have demonstrated that the Adam algorithm
can be used to find the optimal values of $\varepsilon$ and $\Delta$ that
reproduce a target $\langle\sigma_{z}\rangle$.

For all the results presented here, we considered a unique $\rho_{S}(t_{0})$;
however, one can apply this technique to search for the optimal
$\rho_{S}(t_{0})$ given a scalar observable that depends on the steady-state.
For example, in the SB model one can define $\bm{\theta}$ to include the
initial values of the RDM to study the dependence of the SS on initial
conditions.

##  6 Broader Impact

Open quantum systems appear in many areas in physics. Our work can potentially
be used for designing quantum heat engines or materials with specific
desirable properties. This work itself is unlikely to have immediate ethical
issues.

## Acknowledgement

Support by the US Air Force (AFOSR) under grant number FA9550-17-1-0310 is
gratefully acknowledged.

## References

  * Bradbury et al. [2018] J. Bradbury, R. Frostig, P. Hawkins, M. J. Johnson, C. Leary, D. Maclaurin, and S. Wanderman-Milne.  JAX: composable transformations of Python+NumPy programs, 2018\.
  * Breuer and Petruccione [2002] H. P. Breuer and F. Petruccione.  _The theory of open quantum systems_.  Oxford University Press, Great Clarendon Street, 2002.
  * Chen et al. [2018] R. T. Q. Chen, Y. Rubanova, J. Bettencourt, and D. K. Duvenaud.  Neural ordinary differential equations.  In _Advances in neural information processing systems_ , pages 6571–6583, 2018.
  * Egorova et al. [2003] D. Egorova, M. Thoss, W. Domcke, and H. Wang.  Modeling of ultrafast electron-transfer processes: Validity of multilevel redfield theory.  _J. Chem. Phys._ , 119(5):2761–2773, 2003.
  * Jung and Brumer [2020] K. A. Jung and P. Brumer.  Energy transfer under natural incoherent light: Effects of asymmetry on efficiency.  _J. Chem. Phys._ , 153(11):114102, 2020.
  * Kilgour and Segal [2018] M. Kilgour and D. Segal.  Coherence and decoherence in quantum absorption refrigerators.  _Phys. Rev. E_ , 98:012117, 2018.
  * Kingma and Ba [2017] D. P. Kingma and J. Ba.  Adam: A method for stochastic optimization.  _arXiv_ , 1412.6980, 2017.
  * Krantz and Parks [2012] S. G. Krantz and H. R. Parks.  _The implicit function theorem: history, theory, and applications_.  Springer Science & Business Media, 2012.
  * Leggett et al. [1987] A. J. Leggett, S. Chakravarty, A. T. Dorsey, M. P. A. Fisher, Anupam Garg, and W. Zwerger.  Dynamics of the dissipative two-state system.  _Rev. Mod. Phys._ , 59:1–85, 1987.
  * Linden et al. [2010] N. Linden, S. Popescu, and P. Skrzypczyk.  How small can thermal machines be? the smallest possible refrigerator.  _Phys. Rev. Lett._ , 105:130401, 2010.
  * Liu and Segal [2019] J. Liu and D. Segal.  Interplay of direct and indirect charge-transfer pathways in donor–bridge–acceptor systems.  _J. Phys. Chem. B_ , 123(28):6099–6110, 2019\.
  * Nakajima [1958] S. Nakajima.  On quantum theory of transport phenomena.  _Progress of Theoretical Physics_ , 20(6):948–959, 1958.
  * Pachón et al. [2017] L. A. Pachón, J. D. Botero, and P. Brumer.  Open system perspective on incoherent excitation of light-harvesting systems.  _J. Phys. B: At. Mol. Opt. Phys_ , 50(18):184003, 2017.
  * Palenberg et al. [2001] M. A. Palenberg, R. J. Silbey, C. Warns, and P. Reineker.  Local and nonlocal approximation for a simple quantum system.  _J. Chem. Phys._ , 114(10):4386, 2001.
  * Redfield [1965] A. G. Redfield.  The theory of relaxation processes.  In J. S. Waugh, editor, _Advances in Magnetic Resonance_ , volume 1 of _Advances in Magnetic and Optical Resonance_ , pages 1 – 32. Academic Press, 1965.
  * Thoss et al. [2001] M. Thoss, H. Wang, and W. H. Miller.  Self-consistent hybrid approach for complex systems: Application to the spin-boson model with debye spectral density.  _J. Chem. Phys._ , 115(7):2991–3005, 2001.
  * Tscherbul and Brumer [2014] T. V. Tscherbul and P. Brumer.  Long-lived quasistationary coherences in a $v$-type system driven by incoherent light.  _Phys. Rev. Lett._ , 113:113601, 2014.
  * Tscherbul and Brumer [2015] T. V. Tscherbul and P. Brumer.  Partial secular bloch-redfield master equation for incoherent excitation of multilevel quantum systems.  _J. Chem. Phys._ , 142(10):104107, 2015.
  * Xu and Cao [2016] D. Xu and J. Cao.  Non-canonical distribution and non-equilibrium transport beyond weak system-bath coupling regime: A polaron transformation approach.  _Frontiers of Physics_ , 11(4), 2016.
  * Zwanzig [1960] R. Zwanzig.  Ensemble method in the theory of irreversibility.  _J. Chem. Phys._ , 33(5):1338–1341, 1960.