# Global minimization via classical tunneling assisted by collective force
field formation

$\ {}^{\dagger}$F. Caravelli  $\ {}^{\dagger,*}$F. C. Sheldon  †Theoretical
Division (T4) and ∗Center for Nonlinear Studies, Los Alamos National
Laboratory, Los Alamos, New Mexico 87545, USA  F. L. Traversa  MemComputing,
Inc. San Diego, California, 92121

###### Abstract

Simple dynamical models can produce intricate behaviors in large networks.
These behaviors can often be observed in a wide variety of physical systems
captured by the network of interactions. Here we describe a phenomenon where
the increase of dimensions self-consistently generates a force field due to
dynamical instabilities. This can be understood as an unstable (“rumbling")
tunneling mechanism between minima in an effective potential. We dub this
collective and nonperturbative effect a “Lyapunov force" which steers the
system towards the global minimum of the potential function, even if the full
system has a constellation of equilibrium points growing exponentially with
the system size. The system we study has a simple mapping to a flow network,
equivalent to current-driven memristors. The mechanism is appealing for its
physical relevance in nanoscale physics, and to possible applications in
optimization, novel Monte Carlo schemes and machine learning.

Introduction. The dynamics of simple elements interacting in networks can give
rise to _emergent_ behaviors, which are not exhibited in the dynamics of a
single element. These behaviors are inherently multi-element in nature and yet
may admit low-dimensional analytical descriptions in terms of coarse-grained
variables. Examples, such as synchronization and phase transitions are often
widely applicable, as many different systems may reduce to the same effective
model. In this work we demonstrate such a behavior in a model of input-output
networks inspired by memristors. When weakly driven, the system is described
by movement in an effective potential, but at higher driving, instabilities
will cause escapes from local minima of the potential which can be captured as
a tunneling phenomenon.

Tunneling and escape phenomena are a remarkable feature of quantum mechanical
and thermal systems [1, 2]. Thermal and quantum tunneling can occur when a
barrier separates two possible minima, and a particle can escape (tunnel) from
one minimum to the other aided by either thermal fluctuations or spreading of
the wavefunction. If the barrier is not infinitely high, both thermal and
quantum systems have a non-zero probability to tunnel between minima, with a
probability exponentially dependent on the height or width of the barrier. The
discovery of quantum tunneling was a radical paradigmatic shift from classical
mechanics and underlies many important physical phenomena and technologies
(_e.g._ alpha decay, scanning-tunneling microscopy). A similar behavior occurs
in non-equilibrium stochastic dynamical systems. Particles at a certain
temperature can in fact escape from a metastable state after a typical time
proportional to the inverse of the Kramers’ rate [2]. Intuitively, thermal
fluctuations provide sufficient energy for the particle to escape the local
minimum (e.g. a “kick"). An analogy between these behaviors is not surprising
if one considers that the Schrödinger equation can be mapped to a non-
Markovian stochastic system [3], or to a classical particle in a non-local
force field [4].

Figure 1: Representation of the tunneling phenomenon. On the left (a,c), the
memristor potential $V(x,s)$ for an high (a) and a low (c) value of the
voltage $s$ is sketched and two attractors are highlighted through the
directions that internal memristor state take depending on the initial
conditions. On the right (b,d), the memristor potential $V(x,s)$ now
represents either each one of memristors in the network in absence of the
coupling, or equivalently the mean field approximation of the system. In this
case, the emergent Lyapunv force (blue arrows) pushes the system towards the
global minimum creating an effective tunneling through the barrier. In the
middle inset (e) the equivalent schematic of the memristive system is reported
highlighting the role of the potential $V(x,s)$.

Tunneling phenomena are of current interest for their application to
computation, as in quantum annealing [5]. Metastability and tunneling are a
leitmotif for optimization schemes both in a classical and a quantum framework
[6, 7, 8]. In addition to quantum annealing, interest in alternative
approaches to computation and optimization is mounting [9, 10, 11, 12, 13, 14,
15, 16]. For instance, there are proposals to use oscillators or more general
encodings to process information in the frequency domain [15, 16, 10, 17].
Others leverage memory, ranging from near-memory [18] to in-memory computing,
[19, 20, 21] and finishing in memcomputing [11, 22]. These alternatives are
trying to more efficiently solve hard problems across optimization [23, 24,
25, 26] and machine learning [27].

Figure 2: Upper panel: Local ($x_{lm}$, dashed dotted red line) and global
($x_{gm}$, solid green line) minimum location as function of voltage $s$.
Lower panel: Energy Step (blue solid curve) and Barrier $\Delta E$ (dashed
yellow curve) as a function of voltage $s$. Figure 3:  Dynamics of the system
for $s=0.15$ and $s=0.25$, and the mean field basins of the attractions of the
effective potential are represented by the background blue and red coloring.
The initial condition for the elements of $x$ is a uniform distribution in
$[0,1]$. Left panel: trajectories of the dynamical system in the case with
$s=0.15$ and $\chi=0.9$ and in line with the mean field regime (equivalent to
a single memristor).  Right panel: rumbling trajectories of the dynamical
system for the case of the non-convex potential, for $s=0.25$, and in which
the mean field equations fail to capture the dynamics of the network.

Among these alternatives, interest in active or chaotic systems has
experienced a revival, also because of their applicability to optimization
problems [28]. Chaotic dynamics have been posed as an obstacle to reaching
asymptotic fixed points in dynamical system-based computing [29] but there has
also been some indirect evidence that unstable dynamics can improve the
efficiency of some optimization schemes. In particular it has been noted that
in crossbar arrays made of memristive elements, unstable dynamics can lead to
escaping local minima [30] in memristor-implemented simulated annealing.

While escape phenomena are familiar in quantum and thermal settings, there are
no examples, of which we are aware, of classical dynamical systems which
exhibit barrier escapes when the system is a-thermal and passive. In this
article we present an example of barrier tunneling as an emergent, multi-
particle effect in an a-thermal, passive [31] system moving in an effective
potential. Both thermal and quantum tunneling involve a distribution which is
either intrinsic (quantum) or induced by uncertainty (thermal). Here this
distribution is replaced by a coarse-grained representation of the system,
which transitions from a laminar phase (always negative Lyapunov exponents) to
an unstable or transiently unstable regime (local Lyapunov exponents positive
for short transients).

The model. The example considered is inspired by the dynamics of memristor
networks [32, 33, 34, 35, 36, 37], in which a low dimensional effective
potential representation shows barrier escapes. This property is represented
in Fig. 1 and is a novel finding of this paper which we characterize both
analytically and numerically. This effect comes from the intrinsic coupling
among elements, which in the memristor network interpretation is induced by
Kirchhoff’s laws. Due to the non-linear nature of the coupling, we use a
recently obtained exact large-$N$ formula for the matrix inverse, and show
analytically that the effective coarse grained dynamics can be approximated by
the one of a single memristor (e.g. a mean field approximation). When the
effective potential lacks convexity, the emergent Lyapunov force pushes the
system into the absolute minimum of the potential, thus effectively tunneling
through the barriers. Interestingly, this Lyapunov force is directly related
to instabilities in the dynamics.

The dynamical equation for the circuit has been derived in [38], but also
applies to physarum molds [39] or as a dynamic generalization of supply-chain
or input-output models. In particular, it models a flow network which obeys
current and energy conservation and whose edge dynamics are linear in the
currents. A resistor with memory can be described by an effective dynamical
resistance depending on an internal parameter $x$. For instance, $TiO_{2}$
memristors are approximately described by the functional form
$R(x)=R_{\text{off}}(1-x)+xR_{\text{on}}$, where
$R_{\text{on}}<R_{\text{off}}$ are the limiting resistances, and $x\in[0,1]$.
The internal memory parameter evolves according to a simple equation of the
form $\frac{d}{dt}x=\frac{R_{off}}{\beta}I-\alpha
x=\frac{R_{off}}{\beta}\frac{V}{R(x)}-\alpha x$. The parameters $\alpha$ and
$\beta$ are the decay constant and the effective activation voltage per unit
of time respectively, and determine the timescales of the system. For a single
memristor under an applied voltage $S$ we use Ohm’s law $S=RI$ to obtain an
equation for $x(t)$ in adimensional units ($\tau=\alpha t$) given by

$\displaystyle\frac{d}{d\tau}x=\frac{S}{\alpha\beta}\frac{1}{1-\chi
x}-x=-\partial_{x}V(x,s).$ (1)

Here we have defined
$\chi=\frac{R_{\text{off}}-R_{\text{on}}}{R_{\text{off}}}$, with
$0\leq\chi\leq 1$ physically, and $V(x,s)$ as an effective potential.

The dynamics of the one-dimensional system of Eqn. (1) are fully characterized
by gradient descent in the potential

$V(x,s)=\frac{1}{2}x^{2}+\frac{s}{\chi}\log(1-\chi x),$ (2)

with $s=\frac{S}{\alpha\beta}$. The potential can have two local minima and
the energy step between the two is shown in Fig. 2 for various values of $s$
and for $\chi=0.9$. We consider $s$ on the restricted interval for which a
barrier exists as shown in Fig. 2. In this range, and with $\chi$ near 1, the
local minimum can move inside the domain $[0,1]$, and an unstable fixed point
(i.e. the location of the peak of the barrier) emerges, leading to two basins
of attraction (Fig. 2). For $\alpha=\beta=1$ and $\chi=0.9$ this range is
$1/10<s<5/18$ . The requirement that $\chi$, which characterizes the
nonlinearity of the system, be near 1 implies that the phenomenon is
nonperturbative. The asymptotic behavior of this single-element dynamical
system is fully characterized by the simple basins of attraction of the
potential, and presents no exotic features 111 Typically, since one must have
$0\leq x\leq 1$, the equations of motion of the single variables are supplied
with (non-absorbing) cutoff functions, e.g. $\frac{d}{d\tau}x=-W(x)f(x)$, with
$W(x)=1$ for $0\leq x\leq 1$ and zero otherwise. . An analytical solution of
such differential equation has been recently obtained in implicit form [41].

When, instead of a single memory device, we have a circuit composed of many
memristors, each with resistance $R(x_{i})$ and with voltage generators
$S_{i}$ in series, (taken to be constant $S_{i}=\tilde{S}$), the differential
equation (1) generalizes to a system of coupled and nonlinear ordinary
differential equations. The network dynamics equation for the memory elements
$x_{i}(t)$ is [42, 43]:

$\displaystyle\frac{d}{dt}\vec{x}=\frac{1}{\beta}(I-\chi\Omega
X)^{-1}\Omega\vec{S}-\alpha\vec{x},$ (3)

with $\chi=\frac{R_{\text{off}}-R_{\text{on}}}{R_{\text{off}}}<1$, and
$X_{ij}(t)=x_{i}(t)\delta_{ij}$. The matrix $\Omega$ is the projection
operator ($\Omega^{2}=\Omega$) on the vector space of the cycles of
$\mathcal{G}$, the graph representing the circuit [42]. In practice, $\Omega$
can be determined via the directed incidence matrix $B$ of $\mathcal{G}$ as
$\Omega=I-B(B^{t}B)^{-1}B^{t}$ [44, 38, 45, 46, 47]. That $\Omega$ is a
projection operator is a mathematical representation of Kirchhoff’s circuit
laws. In this paper we choose $B$ to be a random matrix in order to abstract
the dynamical system from a particular circuit topology.

The system has a Lyapunov function when $\vec{S}$ is constant (see Supp. Mat.
A), and memristors are passive elements and thus cannot possess positive
Lyapunov exponents (e.g. the system cannot be unstable for long times [31]).
Passive components subject to DC voltage generators will approach an
equilibrium and thus only transient forms of instability are possible. Exact
information about the behavior of Eqn. (3) has been difficult to obtain as the
matrix inverse $(I-\chi\Omega X)^{-1}$ contains the variables $x_{i}(t)$ and
it is therefore difficult to solve analytically. Additionally, Bézout’s
theorem for quadrics of order $2$ with $N$ variables suggests the number of
equilibrium points (stable, unstable or saddle) can be exponential in the
number of memristors [41] (this can be justified intuitively, as each
memristor has 2 fixed points). Despite this complexity, the higher-dimensional
dynamics of networks of memristors can be well-approximated as an effective
one-dimensional system in a coarse-grained representation as we show in the
following section.

Interestingly, in the regime where the potential $V(x,s)$ has multiple minima
and the initial condition of the system lies outside of the attraction basin
of the global minimum, the system exhibits instability. This can be seen in
Fig. 3 (red curves), where for $\chi=0.9$ the system shows qualitatively
different dynamics at different values of $s$. Here we have assigned the
initial conditions of the elements of $\vec{x}$ according to a uniform
distribution on $[0,1]$ and considered two cases for $s$: $s=0.15$ for which
the global minimum is located at $x_{i}\approx 0.2$, and $s=0.25$ for which
the global minimum at $x_{i}=1$ (see Fig. 2). The former case shows laminar
behaviour with the trajectories converging smoothly to $x_{i}\approx 0.2$. On
the other hand, for $s=0.25$, initial conditions for the elements of $\vec{x}$
are almost all within the “local" minimum’s basin of attraction since now the
global minimum is at $x_{i}=1$; now the system shows unstable dynamics,
pushing trajectories to the other side of the barrier and converging to the
global minimum $x_{i}=1$ (see Fig. 3).

We quantify this behavior by examining the distance between copies of the
system with similar initial conditions $\vec{x}^{(1)}(\tau)$ and
$\vec{x}^{(2)}(\tau)$. As a distance, we use the absolute value, which is
essentially an $L_{1}$ measure of the Lyapunov exponents for the dynamics. For
these experiments, we initialize the two systems in the same state
$\vec{x}_{0}$, but with a small deviation $\vec{\epsilon}$, with
$||\vec{\epsilon}||_{1}=0.01$, and where $||\cdot||_{1}$ represents the
1-norm. It is natural to define the quantity
$\mathcal{M}_{\tau}=\frac{1}{N}\sum_{i=1}^{N}|x_{i}^{(1)}(\tau)-x_{i}^{(2)}(\tau)|$,
e.g. a normalized $1$-norm which quantifies the deviation of the trajectories
element by element, and which is convenient because it is naturally normalized
between $0$ and $1$. The results are shown in Fig. 4, where we see that for
the case in which $s=0.15$ (laminar), $\mathcal{M}_{\tau}$ decays to zero
rapidly. In the case $s=0.25$ (unstable), $\mathcal{M}_{\tau}$ grows for a
transient of time of approximately 7-fold longer before asymptotically
reaching zero. For this reason, we dub this phenomenon as a rumbling
transition due to the chaoticity of the individual trajectories, resulting
from the nonlinearity of the interaction.

Figure 4: Growth and decay of
$\mathcal{M}_{\tau}=\frac{1}{N}\sum_{i=1}^{N}|x_{i}(t)-x^{\epsilon}_{i}(t)|$
as a function of time for $\chi=0.9$ and $s=0.15$ and $s=0.25$. The
realization of the noise $\epsilon$ on the initial condition is such that
initially $\mathcal{M}_{0}=0.01$. In the laminar regime $(s=0.15)$,
$\mathcal{M}_{\tau}$ quickly decays to zero, while it grows for a transient of
time in the unstable regime. In order to measure the growth and decay of
$M_{\tau}$, we look at
$\lambda=\frac{1}{\bar{\tau}}\log\frac{\mathcal{M}_{\bar{\tau}}}{\mathcal{M}_{0}}$.
For $s=0.15$, we obtain $\lambda=-0.79(1)$ measures on the whole interval,
while for $s=0.25$ we obtain $\lambda=0.2(1)$, and thus positive. Inset:
Evolution of the Lyapunov force as a function of time and for
$s\in[0.15,0.27]$. Every curve is averaged over $100$ initial conditions and
the shadow area represents the standard deviation on the mean curves.

Coarse-grained dynamics. We now make the intuitive picture above more precise
by introducing an effective mean field potential derived from an approximation
for the matrix inverse. It has been recently proven that the resolvent of
large matrices exhibits statistical regularities, for which it can be
approximated by an effectively one-dimensional matrix [48] which is universal
at the zeroth-order in the limit of weak correlations. The result is rather
general and has applications to various domains of complexity science [49]. If
we define $\vec{f}=\Omega\vec{x}$ and the matrix

$\displaystyle\tilde{A}=\begin{pmatrix}f_{1}(\vec{x})&\cdots&f_{1}(\vec{x})\\\
f_{2}(\vec{x})&\cdots&f_{2}(\vec{x})\\\ \vdots&\ddots&\vdots\\\
f_{N}(\vec{x})&\cdots&f_{N}(\vec{x})\end{pmatrix}=\vec{f}\otimes\vec{1}^{t},$
(4)

we then have the approximate relation (details in Supp. Mat. B)

$\displaystyle(I-\chi\Omega
X)^{-1}=I+\frac{1}{N}\frac{\chi}{1-\frac{1}{N}\chi\sum_{i=1}^{N}f_{i}(\vec{x})}\tilde{A}+O\left(\frac{1}{N}\right).$
(5)

Thus, if we define a coarse-grained variable
$x_{cg}=\frac{1}{N}\sum_{i}f_{i}$, and define
$\langle\Omega\vec{S}\rangle=\frac{1}{N}\sum_{i=1}^{N}(\Omega\vec{S})_{i}$, we
can write an effective one dimensional dynamics:

$\displaystyle\frac{d}{d\tau}x_{cg}$ $\displaystyle=$
$\displaystyle\frac{1}{\alpha\beta}\Big{(}\langle\Omega\vec{S}\rangle+\frac{\chi\langle\Omega\vec{S}\rangle}{1-\chi
x_{cg}}x_{cg}\Big{)}-x_{cg}+\mathcal{L}(\vec{x})$ (6) $\displaystyle=$
$\displaystyle-\partial_{x_{cg}}V(x_{cg},\chi)+\mathcal{L}(\vec{x})$

where $\mathcal{L}(\vec{x})$ is an effective force due to the fact that the
coarse-graining is not exact. We thus see that at the zeroth order the mean
field variable obeys dynamics similar to those of a single memristor, where
now the parameter $\frac{s}{\alpha\beta}$ is replaced by the mean field value
$s=\frac{\langle\Omega\vec{S}\rangle}{\alpha\beta}$. Thus, the equation above
again represents gradient descent dynamics for a one-dimensional system, but
where an effective external force $\mathcal{L}(\vec{x})$ emerges from the
interaction between the large set of (hidden) variables $x_{i}$, and which we
can evaluate numerically on the dynamics. In a sense, while the $N$ particles
somehow feel a similar potential, these can interact as well, leading to a
non-trivial collective dynamics. As a result we find that the variable
$x_{cg}$ is a natural dynamical order parameter to study for our system.

Laminar versus Transient unstable dynamics. We numerically examined the
validity of the coarse-grained dynamics for various ranges of the parameters
$s$ and $\chi$.222The numerical results are obtained using a simple Euler
integration with a step $dt=0.1$, but the results are robust to smaller
values, and agree with independent simulations using Runge-Kutta 4 method and
similar step sizes.

The picture which emerges is represented in Fig. 3. As discussed previously,
we consider a uniform distribution in [0,1] for the initial condition of the
elements of $\vec{x}$ and the laminar case $s=0.15$ for which the global
minimum is located at around $x_{i}\approx 0.2$ and the unstable case $s=0.25$
with the global at $x_{i}=1$ (see Fig. 2). In the former case the majority of
the initial conditions, $x_{i}(0)$ lie in the basin of the global minimum. On
the other hand, for $s=0.25$, almost all initial conditions for the elements
of $\vec{x}$ are on the basin of the “local" minimum as now the global minimum
is at $x_{i}=1$. In the laminar regime the mean-field dynamics accurately fit
the dynamics of the system. This is always the case for $\chi\ll 1$, in which
very few trajectories cross the barrier and the dynamics are therefore smooth.
In this regime, the effective force $\mathcal{L}(x)$ (which we call Lyapunov
force for reasons which become clear soon) is typically small after a short
transient, as we can see the inset in Fig. 4.

This implies that the equilibrium points of the dynamics, which are those
satisfying $\frac{d}{d\tau}x_{cg}=0$, are well approximated by the mean field
dynamics (e.g. the particles follow the same basin of attraction). This gives
the equilibrium value

$\displaystyle
x_{cg}^{*}=\frac{1}{N}\sum_{ij}\Omega_{ij}x_{j}^{*}=\frac{\alpha\beta-\sqrt{\alpha^{2}\beta^{2}-4\alpha\beta\chi\langle\Omega\vec{S}\rangle}}{2\alpha\beta\chi},$
(7)

which matches simulations. The equation above is equivalent to finding values
of $x_{i}\in[0,1]$ such that $\vec{x}\cdot\vec{c}=a$, and $\vec{x}$ is such
that $\vec{x}=\Omega\vec{x}$ (See Supp. Mat. C). This equation has many
multiple possible solutions in high dimensions, explaining the large number of
asymptotic states of the memristive dynamics commonly seen in numerical
experiments [51]. Yet, this also implies that the dynamics can be succinctly
described by a scalar order parameter for the system, and that the mean field
and the effective potential do play a role.

Figure 5: Average asymptotic state
$E\left[\langle\vec{x}(t=\infty)\rangle\right]$ (red curve) of the memristor
coarse grained variable obtained from 400 Monte Carlo instances for each point
$(s,\vec{x})$, on a grid of $50$ points $s\in[0.15,0.27]$ and 30 points
$x_{i}\in[0,1]$ and for $\chi=0.9$, obtained with an Euler integration method
with $dt=0.1$, and tested against a Runge-Kutta 4 integration method for
stability. The transparent plane is for visual aid to show that the average
climbs the barrier.

\begin{overpic}[height=270.30118pt]{probability_vert.png} \put(0.0,98.0){(a)}
\put(0.0,65.0){(b)} \end{overpic}

\begin{overpic}[height=270.30118pt]{probability_sigma.png} \put(0.0,98.0){(c)}
\put(0.0,65.0){(d)} \end{overpic}

Figure 6: Tunneling probability $P_{t}$ from the local to the global minimum
of $V(x,s)$. Insets (a-b) are evaluated for the standard deviation of initial
condition distribution around the local minimum $\sigma=0.5$. (b) probability
$P_{t}$ as a function of the energy step between the two minima and the size
of the system (or equivalently the number of memristors). (a) section of the
plot in (b) for System Size = 100. Blu dots are the numerical evaluation of
$P_{t}$ vs the energy step while the solid red curve is the fitting using the
$\tanh$ functions. The dashed yellow curve is the barrier $\Delta E$ as a
function of the energy step. Insets (c-d) are evaluated for System Size = 100.
(c) probability $P_{t}$ as a function of the energy step between the two
minima and the standard deviation $\sigma$ of initial condition distribution
around the local minimum. (c) section of the plot in (d) for $\sigma=0$ and
0.5. Blu dots and squares are the numerical evaluation of $P_{t}$ vs the
energy step while the solid red curves are the fittings using the $\tanh$
functions. The dashed yellow curve is the barrier $\Delta E$ as a function of
the energy step.

The ability of the coarse-grained dynamics to capture the memristor network
evolution changes abruptly when most $x_{i}$ have initial conditions outside
of the attraction basin of the global minimum, as can be seen in the right
panel of Fig. 3. Fixing $\chi=0.9$ and $s=0.25$ such that the global minimum
is now at $x_{i}=1$ the one-dimensional representation which we obtained at
the leading order fails to captures the dynamics of the system. Individual
trajectories $x_{i}(t)$ now deviates from the coarse-grained value $x_{cg}(t)$
substantially, as we can see in Fig. 3 (right), with clear instabilities which
lead them to the absorbing boundary $x_{i}(t)=1$. While this might seem
irksome at first, we argue that this phenomenon is worthy of attention. In
fact, the effective Lyapunov force $\mathcal{L}(x)$ is now consistently large
along the dynamics, pushing the system towards the absolute minimum. Since the
resolvent expansion of Eqn. (5) is valid in the limit of weak correlations of
the matrix elements, the fact that the Lyapunov force is non-zero is naturally
associated to strong correlations between the dynamical variables.

As described, the system displays this behaviour even when randomly
initialized, and was investigated systematically. To highlight the role of the
basin of attraction of the global minimum we employed Monte Carlo simulations,
randomizing over the initial conditions on $[0,1]$. The results are shown in
Fig. 5. The surface shows the potential $V(x,s)$ as a function of $x_{cg}$ and
$s$. The superimposed red curve represents the asymptotic average position
$E\left[\langle\vec{x}(t=\infty)\rangle\right]$ where $E\left[\cdot\right]$ is
the average over the Monte Carlo samplings and
$\langle\vec{x}\rangle=\frac{1}{N}\sum_{i}x_{i}$ is the average over the
components. As can be seen, before the barrier disappears the dynamics reaches
the right minimum of the potential $V(x,s)$ with probability one. A glimpse of
the structure of the basins of attraction can be obtained by fixing the
initial conditions homogeneously for all but two variables $x_{i}$, and look
at the resulting asymptotic state in these two variables (see the Supp. Mat. E
for details) near the transition point.

Tunneling. To capture this as a tunneling effect, we initialized the system
around the local minimum using the initial condition
$\vec{x}(0)=x_{lm}(s)+\vec{\sigma}$ where $x_{lm}(s)$ is the location of the
coarse grained local minimum as depicted in figure 2, and $\vec{\sigma}$ being
a random vector drawn from the Gaußian distribution $\mathcal{N}(0,\sigma)$.
Given this, we performed Monte Carlo simulations and obtained the probability
$P_{t}$ that the system hops from the local minimum to the global minimum, or,
in other words, the probability $P_{t}$ that the initial condition
$\vec{x}(0)=x_{lm}(s)+\vec{\sigma}$ leads to $\vec{x}(t=\infty)=x_{gm}(s)$. We
also note that $x_{lm}(s)$ shows an abrupt change for $s=s_{crit}\approx
0.206$ as depicted int Fig. 2 as at that point the local and global minima
switch locations. This leads to two different effective tunneling directions
for $s>s_{crit}$ and $s<s_{crit}$ (see the sketch in Fig. 1).

We first investigated how the tunneling emerges as a collective behaviour of
the system. To this end we evaluated the dynamics of systems as the number of
elements increased. Fig. 6-b reports the tunneling probability for
$\sigma=0.5$ and for system sizes from 1 to 100 components and for the entire
range of $s$. We can observe that at small system size the tunneling
probability is proportional to the number of initial conditions that fall into
the global minimum attraction basin to end with an almost perfect sigmoid
($\tanh$), as a function of the energy barrier $\Delta E$ (see also Fig. 6-b).
This transition is smooth and qualitatively equivalent for both tunneling
directions. Therefore the system shows a size-depended transition from
gradient dynamics to a collective tunneling towards the global optimum, led by
the emergent Lyapunov force.

We also investigated how the emergent force depends on the spread of the
initial conditions around $x_{lm}$, i.e., how it depends on $\sigma$. Figure
6-d reports the tunneling probability $P_{t}$ as a function of $\sigma$ for a
system of 50 components. An interesting feature manifests itself: for
$\sigma=0$ there is a non-null height of the barrier $\Delta E$ for each
tunneling direction below which the system is still able to tunnel towards the
global minimum with probability 1 (this is highlighted in figure 6-c). This
shows that the Lyapunov force is an intrinsic collective large-system-size
feature present even in the absence of any randomness in the initial
conditions.

Based on the latter finding, we investigated whether escape from the minimum
might be due to a transient instability of non-normal dynamics [52]. One
possible mechanism would in fact be that via this (local) transient
instability perturbations are amplified for a short time and the system
reaches the basin of attraction of another minimum. To examine this we studied
the local stability of the system near $x_{lm}$ via the Jacobian matrix. While
the details are provided in App. D, we analyzed a large number of local minima
via Monte Carlo, finding that the non-normality of the Jacobian does not lead
to any amplification phenomenon. This also leaves us with the only option that
the instability of Fig. 4 is in fact due to a non-perturbative and cooperative
phenomenon between the dynamical variables. From the point of view of the
effective potential $V(x,s)$, this analysis implies that the global minimum
lying behind a barrier is reached (via the instabilities above mentioned) via
a dynamics assisted by the effective Lyapunov force which pushes the system
towards it.

Conclusions. Traditional escape phenomena in thermal and quantum systems are a
consequence of either thermal fluctuations or the spreading of a wavefunction.
We have presented an additional setting in which barrier escapes emerge in the
effective description of a multiparticle system. We presented a class of
dynamical system, derived from networks of memristors, that can be mapped to
an effective one-dimensional system. The resulting dynamics are characterized
by a potential which depends on the (mean-field) external applied voltage, and
an effective force which occurs when the system becomes Lyapunov unstable. The
result of this instability is that the effective representation may jump
between basins of attraction, even through a barrier. The explanation we
provide is simple but difficult to capture analytically and is due to the
large number of directions in which the system can travel, via a sequence of
saddle points in the dynamics, and from which the local instability emerges.
This is compatible with previous investigations of memristive circuit dynamics
applied to logic gates [53, 54], and the observation of instanton-like
behavior in the dynamics [55, 56].

Our analysis of the effective “Lyapunov" force provides an explanation of the
observed tunneling as an epiphenomenon of the interaction between the memory
elements. In principle this might fit into a coarse graining argument of the
hidden variables, due to a Zwanzig expansion [57], which will be considered in
the future. However, we wish to stress that the system we studied is a-thermal
(_i.e._ no stochastic forces were introduced in our analysis) and that the
effective Lyapunov force is far from random.

The main message of this paper is that the introduction of hidden variables in
a dynamical system can lead to transitions between local and global minima of
the effective description via instabilities in the full system. In this sense,
this type of behavior is novel and worth further exploration. In quantum and
statistical field theory, tunneling phenomena are described by instantons,
e.g. solutions in which the field tunnels in a finite time between local
minima. Instantons, originally introduced in cosmology [58], are nonlinear
solutions to the field equation which are non-local in time, and represent
tunneling between two local minima. Typically, instantonic behavior is
associated to some form of stochasticity, whether thermal or quantum in
nature. In this paper we have shown that instantonic behavior can also be an
emergent phenomenon [59] thanks to the nonlinear interaction between many
components. We hope that our paper can spark further interest in the study of
tunneling in dynamical systems via “hidden variables", with their interesting
properties, and with important applications to optimization [30, 36, 55, 56],
and machine learning [sompolinsky, ganguli].

Acknowledgements. We acknowledge the support of NNSA for the U.S. DoE at LANL
under Contract No. DE-AC52-06NA25396. FC was also financed via DOE-LDRD grants
PRD20170660 and PRD20190195, and would like to thank S. Bartolucci, F.
Caccioli and P. Vivo for collaborations on the resolvent formula which enabled
us to write this paper, and C. Coffrin and Y.-T. Lin for comments on the
draft. FCS was also supported by a Center for Nonlinear Studies fellowship.

## References

  * [1] J. J. Sakurai and J. Napolitano, Modern Quantum Mechanics.  Cambridge University Press, sep 2017.
  * [2] P. Hänggi, P. Talkner, and M. Borkovec, “Reaction-rate theory: fifty years after kramers,” Reviews of Modern Physics, vol. 62, pp. 251–341, apr 1990\.
  * [3] E. Nelson, Quantum Fluctuations.  Princeton University Press, sep 2020.
  * [4] D. Bohm, “A suggested interpretation of the quantum theory in terms of "hidden" variables. II,” Physical Review, vol. 85, pp. 180–193, jan 1952\.
  * [5] W. D. Oliver, “Quantum computing takes flight,” Nature, vol. 574, pp. 487–488, oct 2019.
  * [6] G. E. Santoro, “Theory of quantum annealing of an ising spin glass,” Science, vol. 295, pp. 2427–2430, mar 2002.
  * [7] R. Hamerly and et al., “Experimental investigation of performance differences between coherent ising machines and a quantum annealer,” Science Advances, vol. 5, p. eaau0823, may 2019.
  * [8] C. Baldassi and R. Zecchina, “Efficiency of quantum vs. classical annealing in nonconvex learning problems,” Proceedings of the National Academy of Sciences, vol. 115, pp. 1457–1462, jan 2018.
  * [9] J. L. Hennessy and D. A. Patterson, “A new golden age for computer architecture,” Communications of the ACM, vol. 62, pp. 48–60, jan 2019\.
  * [10] S. K. Vadlamani, T. P. Xiao, and E. Yablonovitch, “Physics successfully implements lagrange multiplier optimization,” Proceedings of the National Academy of Sciences, vol. 117, pp. 26639–26650, oct 2020.
  * [11] F. L. Traversa and M. Di Ventra, “Universal memcomputing machines,” IEEE Transactions on Neural Networks and Learning Systems, vol. 26, pp. 2702–2715, nov 2015.
  * [12] B. Sutton, K. Y. Camsari, B. Behin-Aein, and S. Datta, “Intrinsic optimization using stochastic nanomagnets,” Scientific Reports, vol. 7, mar 2017.
  * [13] J. Torrejon, M. Riou, F. A. Araujo, S. Tsunegi, G. Khalsa, D. Querlioz, P. Bortolotti, V. Cros, K. Yakushiji, A. Fukushima, H. Kubota, S. Yuasa, M. D. Stiles, and J. Grollier, “Neuromorphic computing with nanoscale spintronic oscillators,” Nature, vol. 547, pp. 428–431, jul 2017.
  * [14] B. Molnár, F. Molnár, M. Varga, Z. Toroczkai, and M. Ercsey-Ravasz, “A continuous-time MaxSAT solver with high analog performance,” Nature Communications, vol. 9, nov 2018.
  * [15] F. Böhm, G. Verschaffelt, and G. V. der Sande, “A poor man’s coherent ising machine based on opto-electronic feedback systems for solving optimization problems,” Nature Communications, vol. 10, aug 2019.
  * [16] D. Pierangeli, G. Marcucci, and C. Conti, “Large-scale photonic ising machine by spatial light modulation,” Physical Review Letters, vol. 122, p. 213902, may 2019.
  * [17] G. Csaba and W. Porod, “Coupled oscillators for computing: A review and perspective,” Applied Physics Reviews, vol. 7, p. 011302, mar 2020.
  * [18] G. Singh, L. Chelini, S. Corda, A. J. Awan, S. Stuijk, R. Jordans, H. Corporaal, and A.-J. Boonstra, “Near-memory computing: Past, present, and future,” Microprocessors and Microsystems, vol. 71, p. 102868, nov 2019\.
  * [19] D. Ielmini and H.-S. P. Wong, “In-memory computing with resistive switching devices,” Nature Electronics, vol. 1, pp. 333–343, jun 2018.
  * [20] F. L. Traversa, F. Bonani, Y. V. Pershin, and M. Di Ventra, “Dynamic computing random access memory,” Nanotechnology, vol. 25, p. 285201, 2014\.
  * [21] A. Sebastian, M. L. Gallo, R. Khaddam-Aljameh, and E. Eleftheriou, “Memory devices and applications for in-memory computing,” Nature Nanotechnology, vol. 15, pp. 529–544, mar 2020.
  * [22] M. Di Ventra and F. L. Traversa, “Perspective: Memcomputing: Leveraging memory and physics to compute efficiently,” Journal of Applied Physics, vol. 123, p. 180901, may 2018.
  * [23] S. Kirkpatrick, C. D. Gelatt, and M. P. Vecchi, “Optimization by simulated annealing,” Science, vol. 220, pp. 671–680, may 1983.
  * [24] F. Glover and G. A. Kochenberger, eds., Handbook of Metaheuristics.  Springer US, 2003.
  * [25] M. Dorigo and T. Stützle, Ant Colony Optimization.  The MIT Press, 2004.
  * [26] X.-S. Yang, Engineering Optimization.  John Wiley & Sons, Inc., jun 2010.
  * [27] Y. LeCun, Y. Bengio, and G. Hinton, “Deep learning,” Nature, vol. 521, pp. 436–444, may 2015.
  * [28] M. Ercsey-Ravasz and Z. Toroczkai, “The chaos within sudoku,” Scientific Reports, vol. 2, oct 2012.
  * [29] T. Tél, “The joy of transient chaos,” Chaos: An Interdisciplinary Journal of Nonlinear Science, vol. 25, p. 097619, sep 2015.
  * [30] J. Yang, Q. Duan, Y. Wang, T. Zhang, T. Yang, and R. Huang, “Transiently chaotic simulated annealing based on intrinsic nonlinearity of memristors for efficient solution of optimization problems,” Science Advances, vol. 6, p. eaba9901, aug 2020.
  * [31] S. Grivet-Talocia and B. Gustavsen, Passive macromodeling: Theory and applications.  John Wiley & Sons, 2015.
  * [32] L. Chua, “Memristor-the missing circuit element,” IEEE Transactions on Circuit Theory, vol. 18, no. 5, pp. 507–519, 1971.
  * [33] L. Chua and S. M. Kang, “Memristive devices and systems,” Proceedings of the IEEE, vol. 64, no. 2, pp. 209–223, 1976.
  * [34] D. B. Strukov, G. S. Snider, D. R. Stewart, and R. S. Williams, “The missing memristor found,” Nature, vol. 453, pp. 80–83, may 2008.
  * [35] F. Caravelli and J. P. Carbajal, “Memristors for the curious outsiders,” Technologies, vol. 6, p. 118, dec 2018.
  * [36] M. Di Ventra and Y. V. Pershin, “The parallel approach,” Nature Physics, vol. 9, pp. 200–202, apr 2013.
  * [37] C. Du, F. Cai, M. A. Zidan, W. Ma, S. H. Lee, and W. D. Lu, “Reservoir computing using dynamic memristors for temporal information processing,” Nature Communications, vol. 8, dec 2017.
  * [38] F. Caravelli, “Locality of interactions for planar memristive circuits,” Physical Review E, vol. 96, p. 052206, nov 2017.
  * [39] A. Johannson and J. Zou, “A slime mold solver for linear programming problems,” in Conference on Computability in Europe, pp. 344–354, Springer, 2012.
  * [40] Typically, since one must have $0\leq x\leq 1$, the equations of motion of the single variables are supplied with (non-absorbing) cutoff functions, e.g. $\frac{d}{d\tau}x=-W(x)f(x)$, with $W(x)=1$ for $0\leq x\leq 1$ and zero otherwise.
  * [41] F. C. Sheldon, F. Caravelli, and C. Coffrin, “Fully analog memristive circuits for optimization tasks: a comparison,” Sept. 2020.
  * [42] F. Caravelli, F. L. Traversa, and M. Di Ventra, “Complex dynamics of memristive circuits: Analytical results and universal slow relaxation,” Physical Review E, vol. 95, p. 022140, feb 2017.
  * [43] F. Caravelli, “Asymptotic behavior of memristive circuits,” Entropy, vol. 21, p. 789, aug 2019.
  * [44] F. Caravelli, “The mise en scéne of memristive networks: effective memory, dynamics and learning,” International Journal of Parallel, Emergent and Distributed Systems, vol. 33, pp. 350–366, may 2017.
  * [45] A. Zegarac and F. Caravelli, “Memristive networks: From graph theory to statistical physics,” EPL Perspectives, vol. 125, p. 10001, jan 2019\.
  * [46] J. W. Nilsson and S. Riedel, Electric Circuits (9th ed),.  Pearson Education, Saddle River NJ, 2011.
  * [47] B. Bollobás, Modern Graph Theory.  Springer New York, 1998.
  * [48] S. Bartolucci, F. Caccioli, F. Caravelli, and P. Vivo, “Inversion-free leontief inverse: statistical regularities in input-output analysis from partial information,” Sept. 2020.
  * [49] S. Bartolucci, F. Caccioli, F. Caravelli, and P. Vivo, “Universal rankings in complex input-output organizations,” Sept. 2020.
  * [50] The numerical results are obtained using a simple Euler integration with a step $dt=0.1$, but the results are robust to smaller values, and agree with independent simulations using Runge-Kutta 4 method and similar step sizes.
  * [51] F. Caravelli and F. C. Sheldon, “Phases of memristive circuits via an interacting disorder approach,” Aug. 2020.
  * [52] L. N. Trefethen and M. Embree, Spectra and pseudospectra: the behavior of nonnormal matrices and operators.  Princeton University Press, 2005.
  * [53] F. L. Traversa and M. Di Ventra, “Polynomial-time solution of prime factorization and NP-complete problems with digital memcomputing machines,” Chaos: An Interdisciplinary Journal of Nonlinear Science, vol. 27, p. 023107, 2017.
  * [54] M. Di Ventra and F. L. Traversa, “Absence of periodic orbits in digital memcomputing machines with solutions,” Chaos: An Interdisciplinary Journal of Nonlinear Science, vol. 27, p. 101101, oct 2017.
  * [55] M. Di Ventra, F. L. Traversa, and I. V. Ovchinnikov, “Topological field theory and computing with instantons,” Annalen der Physik, vol. 529, p. 1700123, aug 2017.
  * [56] S. R. Bearden, H. Manukian, F. L. Traversa, and M. Di Ventra, “Instantons in self-organizing logic gates,” Physical Review Applied, vol. 9, p. 034029, mar 2018.
  * [57] R. Zwanzig, Nonequilibrium statistical mechanics.  Oxford New York: Oxford University Press, 2001.
  * [58] S. Coleman, “Fate of the false vacuum: Semiclassical theory,” Physical Review D, vol. 15, pp. 2929–2936, may 1977.
  * [59] P. W. Anderson, “More is different,” Science, vol. 177, pp. 393–396, aug 1972.
  * [60] E. Ott, Chaos in Dynamical Systems.  Cambridge University Press, aug 2002.

Supplementary Material

## A: Proof of Lyapunov function for sharp boundaries

In this section we provide a proof of the fact that the system is dissipative
when each state variable is far from the boundary $x_{i}=1$ [41, 51]. We
consider the case of infinitely sharp boundaries, e.g.
$\lim_{p\rightarrow\infty}W_{p}(x)=\theta(x)\theta(1-x)$. In order to prove
it, we attempt to write a Lyapunov function of the form for the differential
equation (when boundaries are not considered)

$\displaystyle\frac{d}{dt}\vec{x}=\frac{1}{\beta}(I-\chi\Omega
X)^{-1}\Omega\vec{S}-\alpha\vec{x},$ (8)

which we can rewrite as

$\displaystyle\frac{1}{\alpha}(I-\chi\Omega
X)\frac{d}{dt}\vec{x}=\frac{1}{\alpha\beta}\vec{s}-(I-\chi\Omega X)\vec{x},$
(9)

where we wrote $\Omega\vec{S}=\vec{s}$.

$L=\frac{1}{3}\vec{x}^{T}X\vec{x}-\frac{\chi}{4}\vec{x}^{T}X\Omega
X\vec{x}-\frac{1}{2\alpha\beta}\vec{x}^{T}X\vec{s}.$ (10)

After taking a time derivative of the Lyapunov function above, we get

$\displaystyle\frac{dL}{dt}$
$\displaystyle=\dot{\vec{x}}^{T}\left(X\vec{x}-\chi X\Omega
X\vec{x}-\frac{1}{\alpha\beta}X\vec{s}\right)=\dot{\vec{x}}^{T}X\left(\vec{x}-\chi\Omega
X\vec{x}-\frac{1}{\alpha\beta}\vec{s}\right)$
$\displaystyle=-\dot{\vec{x}}^{T}\
X\left(\frac{1}{\alpha\beta}\vec{s}-(I-\chi\Omega X)\vec{x}\right)\text{ which
is Eqn. (\ref{eq:rewritten})}$
$\displaystyle=-\frac{1}{\alpha}\dot{\vec{x}}^{T}X(I-\chi\Omega
X)\dot{\vec{x}}=-\frac{1}{\alpha}\dot{\vec{x}}^{T}(X-\chi X\Omega
X)\dot{\vec{x}}$
$\displaystyle=-\frac{1}{\alpha}\dot{\vec{x}}^{T}\sqrt{X}(I-\chi\sqrt{X}\Omega\sqrt{X})\sqrt{X}\dot{\vec{x}}=-\frac{1}{\alpha}||\sqrt{X}\dot{\vec{x}}||^{2}_{(I-\chi\sqrt{X}\Omega\sqrt{X})}$
(11)

We obtain that if $I-\chi\sqrt{X}\Omega\sqrt{X}\succ 0$ (and $X\neq 0$), the
Lyapunov function has negative derivative always. This occurs if
$x_{i}\in[0,1]$, which is a key factor in the dynamics, it is not hard to see
that $\sqrt{X}\Omega\sqrt{X}\prec 1$ (if $x_{i}\in[0,1]$), and since $\chi<1$,
the Lyapunov property applies also in this case. The existance of a Lyapunov
function also implies that there is a global notion of “energy" beyond the
mean-field potential provided in the main text.

## B: Large-N resolvent and mean field equations

In this section we discuss the theorem used in the main text to derive the
mean field potential. Such theorem has been derived with the intention of
studying resolvent matrices for generic random variables whose correlations
are weak (in a sense defined below).

Theorem. Let $A=(a_{ij})$ be a random $N\times N$ matrix, characterized by the
joint probability density function of the entries
$P_{A}(a_{11},\ldots,a_{NN})$. Let $a_{ij}\geq 0$, $\rho(A)<1$, $\beta>1$ and

$\displaystyle\langle a_{ij}\rangle_{A}$ $\displaystyle=z_{i}/N>0$ (12)
$\displaystyle\langle a_{ij}a_{k\ell}\rangle_{A}-\langle
a_{ij}\rangle_{A}\langle a_{k\ell}\rangle_{A}$
$\displaystyle\sim\frac{C_{ij}}{N^{\beta}}\delta_{ik}\delta_{j\ell}\qquad\text{as
}N\to\infty$ (13)

for all $i,j$, with $\langle(\cdot)\rangle_{A}$ defined as

$\langle(\cdot)\rangle_{A}=\int\mathrm{d}a_{11}\cdots\mathrm{d}a_{NN}P_{A}(a_{11},\ldots,a_{NN})(\cdot)\
,$ (14)

and $C_{ij}$ be constants. Then we have:

$\displaystyle(I-A)^{-1}=I+\frac{1}{N}\frac{1}{1-\bar{z}}\vec{z}\otimes\vec{1}^{t}+O(\frac{1}{N^{\gamma}}).$
(15)

where $\bar{z}=(1/N)\sum_{i=1}^{N}z_{i}$, and $z_{i}=\sum_{j=1}^{N}A_{ij}$,
with $\gamma>0$. The results above hold irrespective of the precise form of
$P_{A}$. For a proof, see [48], and an application of the formula above to
complex datasets [49]. While the theorem is stated for a particular form of
the correlation between the elements, as a matter of fact it is more general,
and it applies to correlation functions which are sufficiently sparse.
Nonetheless, for finite $N$ the expansion has some corrections which are
important, and that we estimate numerically in the main text.

In our case, the ensemble $A$ is the one of the random projector operators. We
will make hereby the connection between the theorem and our setup. The formula
for the resolvent is exact when the span of the operator $\Omega$ is one
dimensional, as it can be promptly seen using the Sherman-Morrison formula.

In order to derive a mean-field (or coarse-grained) dynamics for the network
of memristors, we now use a recently obtained formula for the resolvent. Let
us now consider the equation:

$\frac{d}{dt}\vec{x}=\frac{1}{\beta}(I-\chi\Omega
X)^{-1}\Omega\vec{S}-\alpha\vec{x}.$ (16)

The technical condition $\lim_{N\rightarrow\infty}N\langle
a_{ij}\rangle_{A}=c_{ij}$ applies in the case of memristive dynamics if we
consider the fact that $\Omega$ a projector operator. As shown in [43],
typically, $\Omega_{ij}\approx\frac{c}{N}$ because of the condition
$\Omega^{2}=\Omega$. As such, if we identify $A=\chi\Omega X$, since we impose
$0<x_{i}<1$ dynamically, and $\chi<1$, then we have that
$a_{ij}<\frac{\chi|c|}{N}$, which satisfies the technical condition. The
numerical results are obtained for a number of memristors equal to $N=200$,
and the of the span of $|Span(\Omega)|=N-M$ with $M=50$. For random matrices,
we generate $\Omega=A^{t}(AA^{t})^{-1}A$, with $A$ a random matrix of size
$N\times M$, and $a_{ij}$ randomly distributed in $[0,1]$.

We now discuss how to apply the result above to the dynamics of memristors.
Let us define $f_{i}(\vec{x})=\sum_{j}\Omega_{ij}x_{j}$. It is not hard to see
from the definition above that $z_{i}=\chi f_{i}(\vec{x})$. Let us define
$\vec{f}=\\{f_{i}(\vec{x})\\}$. We then have

$\displaystyle\tilde{A}=\begin{pmatrix}f_{1}(\vec{x})&\cdots&f_{1}(\vec{x})\\\
f_{2}(\vec{x})&\cdots&f_{2}(\vec{x})\\\ \vdots&\ddots&\vdots\\\
f_{N}(\vec{x})&\cdots&f_{N}(\vec{x})\end{pmatrix}=\vec{f}\otimes\vec{1}^{t}.$
(17)

The expansion above can be rewritten as

$\displaystyle(I-\chi\Omega
X)^{-1}=I+\frac{1}{N}\frac{\chi}{1-\chi\frac{1}{N}\sum_{i=1}^{N}f_{i}(\vec{x})}\vec{f}\otimes\vec{1}^{t}+O(\frac{1}{N})$
(18)

from which we obtain, if we call $\Omega\vec{S}=\Omega\vec{S}$

$\frac{d}{dt}\vec{x}=\frac{1}{\beta}\Big{(}\Omega\vec{S}+\frac{\chi\frac{1}{N}\Omega\vec{S}\cdot\vec{1}}{1-\chi\frac{1}{N}\vec{f}\cdot\vec{1}}\vec{f}\Big{)}-\alpha\vec{x}+O(\frac{1}{N})$
(19)

Now that we have the identity $\Omega\vec{f}=\vec{f}$ and $\Omega\vec{S}=s$.
Thus, if we project the equation using $\Omega$, we have

$\frac{d}{dt}\vec{f}=\frac{1}{\beta}\Big{(}\Omega\vec{S}+\frac{\chi\frac{1}{N}\Omega\vec{S}\cdot\vec{1}}{1-\chi\frac{1}{N}\vec{f}\cdot\vec{1}}\vec{f}\Big{)}-\alpha\vec{f}+O(\frac{1}{N}).$
(20)

The asymptotic states are thus determined by the equation

$\displaystyle\frac{1}{\beta}\Big{(}\Omega\vec{S}+\frac{\chi\frac{1}{N}\Omega\vec{S}\cdot\vec{1}}{1-\chi\frac{1}{N}\vec{f}\cdot\vec{1}}\vec{f}\Big{)}-\alpha\vec{f}=0$
(21)

or

$\displaystyle\vec{f}=\Omega\vec{x}=\frac{1}{1-\frac{\chi}{\alpha\beta}\frac{\langle\Omega\vec{S}\rangle}{1-\chi
x_{cg}}}\frac{\Omega\vec{S}}{\alpha\beta},$ (22)

where we defined $\langle\vec{G}\rangle=\frac{1}{N}\vec{G}\cdot\vec{1}$. The
asymptotic states can be obtained up to a gauge transformation
$\vec{x}^{*}=\vec{x}^{*}+(I-\Omega)\vec{r}$ for an arbitrary $\vec{r}$. The
equation above however cannot be solved without first finding a value for
$x_{cg}.$ An equation for $x_{cg}=\frac{1}{N}\vec{f}\cdot\vec{1}$ can be
obtained by summing over all the indices on both sides of equation (20) and
dividing by $N$, obtaining

$\displaystyle\partial_{t}x_{cg}=\frac{1}{\beta}\Big{(}\langle\Omega\vec{S}\rangle+\frac{\chi\langle\Omega\vec{S}\rangle}{1-\chi
x_{cg}}x_{cg}\Big{)}-\alpha x_{cg}+\mathcal{L}(\vec{x})$ (23)

where $\mathcal{L}(\vec{x})$ is an effective force we will discuss in a
moment.

In order for the equilibrium points to be compatible , we thus have to solve
also for the equation

$\displaystyle\frac{1}{\beta}\Big{(}\langle\Omega\vec{S}\rangle+\frac{\chi\langle\Omega\vec{S}\rangle}{1-\chi
x_{cg}}x_{cg}\Big{)}-\alpha x_{cg}=0$ (24)

which is equivalent to a mean-field equation of the form:

$\displaystyle
x_{cg}=\frac{\langle\Omega\vec{S}\rangle}{\alpha\beta}(1+\frac{\chi
x_{cg}}{1-\chi x_{cg}})$ (25)

whose solution is given by

$\displaystyle
x_{cg}^{*}=\frac{1}{N}\sum_{ij}\Omega_{ij}x_{j}^{*}=\frac{\alpha\beta-\sqrt{\alpha^{2}\beta^{2}-4\alpha\beta\chi\langle\Omega\vec{S}\rangle}}{2\alpha\beta\chi}$
(26)

By setting $\frac{\langle\Omega\vec{S}\rangle}{\alpha\beta}=s$ we obtain a
mapping to the one dimensional case. Let us now discuss the emergence of the
lack of convexity of the potential $V(x,s)$ as a function of $s$;
specifically, let us analyze the situation in which the potential has a
minimum or a maximum. The coarse grained variable is given, as mentioned
earlier, by $x_{cg}=\frac{1}{N}\sum_{ij}\Omega_{ij}x_{j}$. First, we focus on
the potential on a compact support $D=[x_{min},x_{max}]$, determined via the
condition $x_{min}=min_{\vec{x}\in[0,1]^{N}}(x_{cg})$ and
$x_{max}=max_{\vec{x}\in[0,1]^{N}}(x_{cg})$. For $s<0$ the potential has only
a local minimum (which is also an absolute minimum) at
$x_{cg}=\text{max}(0,x_{min})$. For $s>0$, the potential of Eqn. (2) can have
at most two points in which its derivative in $x$ is zero, given by the values
$x^{\pm}=\frac{1\pm\sqrt{1-4s\chi}}{2\chi}$. While $x^{-}$ is always a local
minimum, $x^{+}$ is always a local maximum of the potential. However, there is
an intermediate set of values, given by $s=[0,\frac{1}{4\chi}]$, in which the
potential has only a single minimum in the domain $[x_{min},x_{max}]$, as
$x^{+}$ is outside the domain. In this case, the local minimum $x^{-}$ is also
an absolute minimum in the domain $[0,1]$. However, if $\chi$ is close enough
to $1$ and $s$ large enough, we have that $x^{+}$ moves inside the domain, and
thus the potential has two local minima in $[0,1]$, one at $x^{-}$ and one at
$x_{max}$. If $s$ increases further, we can have that $V(x_{max})<V(x^{-})$,
implying that the potential is not non convex, and a barrier of height $\Delta
E(s)=V_{s}(x^{+},\chi)-V_{s}(x^{-},\chi)$ emerges between the the local
minimum $x^{-}$ and the absolute minimum at $x_{xmax}$.

## C: Projected Dynamics and equilibrium points

One of the advantages of using an analytical framework to study memristor
dynamics is that many subtleties of nonlinear circuits become explicit in this
setup. In this regard, it is interesting to recall the fact that
$\Omega^{2}=\Omega$. As mentioned earlier, the fact that this matrix is a
projector operator is essentially a consequence of the Kirchhoff’ laws. Such
freedom implies a redundancy in the way in which we can effectively represent
the dynamics, as investigated in earlier papers [44, 42, 38, 43]. In fact, in
order to obtain the matrix $\Omega$ from the circuit, one has the liberty of
choosing a particular reference spanning tree, whose details are clarified in
[45]. Such freedom is not hidden, but is evident from the fact that if we
perform the transformation $\vec{S}^{\prime}=\vec{S}+(I-\Omega)\vec{k}$ for an
arbitrary vector $\vec{k}$, the dynamics is unchanged (this is reminiscent of
the “toy" gauge freedom which is typically present in circuits). We now see
that such freedom is also present in the representation of the fixed points.
We can in fact write $\vec{x}_{j}^{\prime}=\vec{x}_{j}+(I-\Omega)\vec{k}$,
without effectively changing the coarse grained variable fixed point state.

Let us consider the two formulations of the dynamics, based on the fact that
$\vec{x}=\Omega\vec{x}+(I-\Omega)\vec{x}=\vec{x}_{c}+\vec{x}_{e}$, in which
$\vec{x}_{c}\cdot\vec{x}_{e}=0$ because $\Omega$ is an orthogonal projector.

Given the dynamics introduced in the bulk of the paper, we have

$\displaystyle\frac{d\vec{x}}{dt}=\frac{1}{\beta}(I-\chi\Omega
X)^{-1}\Omega\vec{S}-\alpha\vec{x},$ (27)

from which, projecting via $\Omega$ and $\Omega_{B}=I-\Omega$, we have, using
the series expansion for the matrix expansion, that

$\displaystyle\frac{d\vec{x}_{c}}{dt}$ $\displaystyle=$
$\displaystyle\frac{1}{\beta}(I-\chi\Omega(X_{e}+X_{c}))^{-1}\Omega\vec{S}-\alpha\vec{x}_{c},$
$\displaystyle\frac{d\vec{x}_{e}}{dt}$ $\displaystyle=$
$\displaystyle-\alpha\vec{x}_{e}.$ (28)

It follows that since $\vec{x}_{e}$ falls exponentially to zero, the fixed
points are those such that $\vec{x}_{e}=0$. Given $\vec{x}_{c}=\vec{f}$, the
basin of attraction for the fixed points $\vec{f}^{*}$ is up to an arbitrary
vector $\vec{x}_{e}=(I-\Omega)\vec{k}$ added to the initial conditions. For
the asymptotic fixed points emerging from the mean field dynamics, this
implies that, since $x^{*}_{cg}$ must be such that
$\frac{1}{N}\sum_{ij}\Omega_{ij}x_{j}^{*}=a$, then necessarily
$\sum_{j}\Omega_{ij}x_{j}^{*}=x_{i}^{*}$. Thus, if $\Omega_{ij}$ has a span of
$N-M$ vectors, $\vec{x}^{*}=\sum_{k=1}^{M-N}a_{k}\vec{n}_{k}$, where
$\vec{n}_{k}$’s are the basis vectors of $\text{Span}(\Omega).$ It follows
that the solutions of the problem are such that, if we call
$\tilde{n}_{k}=\sum_{i}(\vec{n}_{k})_{i}$, the solutions are of the form
$\sum_{k}a_{k}\tilde{n}_{k}=a$.

## D: Local stability analysis and pseudospectral triviality

For the single memristor dynamical system, the stability of the equilibrium is
determined by the condition $-\partial_{x}^{2}V(x^{*})>0$, with
$\partial_{x}V(x)=0$. This is equivalent to the condition

$\displaystyle\frac{2\alpha\beta\left(-\alpha\beta+\sqrt{\alpha\beta(\alpha\beta-4s\chi)}+4s\chi\right)}{\left(\sqrt{\alpha\beta(\alpha\beta-4s\chi)}-\alpha\beta\right)^{2}}<0.$
(29)

For reference, in the main text we consider the values $\alpha=\beta=1$ and
$\chi=0.9$, and the dynamical instability occurs around $s\approx 0.27777(8)$.

The dynamics of our system can be described by the differential equations:

$\frac{d{\vec{x}}}{dt}={\vec{f}}({\vec{x}})$ (30)

where ${\vec{f}}=[f_{1},f_{2}\cdots,f_{S}]^{T}$ is a set of non-linear
functions [60].

Since we are interested in our system escaping from the equilibrium point, we
study numerically the equilibria ${\vec{x}*}$ which must satisfy

${\vec{x}}({\vec{x}^{*}})=0.$ (31)

Since the system can have an exponential number of fixed points, such search
cannot be exaustive. However, finding these points is relatively easy and can
be done numerically by initializing the dynamics randomly and letting the
system evolve into one of these minima.

The local dynamics near each equilibrium point is described by the Jacobian
matrix evaluated at an equilibrium point. An equilibrium point is stable if
under any infinitesimally small perturbation, $\Delta{\vec{x}}(0)$, eventually
decays to zero, i.e., $\lim_{t\rightarrow\infty}\Delta{\vec{x}}(t)=0$. In the
vicinity of an equilibrium point, the time evolution of a perturbation can be
written as

$\Delta{\vec{x}}(t)=e^{{DJ}t}\Delta{\vec{x}}(0).$ (32)

Therefore, the spectrum of $DJ$ is relevant for local stability analysis. If
$\Lambda($ M$)$ is the set of eigenvalues of ${\bf M}$, then the equilibrium
point is stable if all eigenvalues have negative real part. For stability near
the boundaries, the dynamical system is given by

$\displaystyle\frac{d}{dt}x_{j}=W_{j}(x_{j})f_{j}(\vec{x})$ (33)

where $W_{p}(\vec{x})$ is a window function enforcing $0\leq x\leq 1$ and
$f_{j}(\vec{x})=(\frac{1}{\beta}(I-\chi\Omega
X)^{-1}\Omega\vec{S}-\alpha\vec{x})_{j}$. In the following we consider the
non-absorbing window function
$W(x_{i})=\theta(-f_{i})\theta(x_{i})+\theta(f_{i})\theta(1-x_{i})$, where we
relaxed to $\theta(a)$ to $\theta_{p}(a)=\frac{1}{2}+\frac{1}{2}\tanh(pa)$,
for numerical stability[11]. The numerical results shown in the paper are
obtained for $p=100$, which is a good approximation of the Heaviside theta
function $\theta(x)$. The full Jacobian of the dynamics is given by

$\displaystyle
DJ_{ij}^{f}=\partial_{x_{i}}(W_{j}(x_{j})f_{j}(\vec{x}))=W_{j}(x_{j})\partial_{x_{i}}(f_{j}(\vec{x}))+f_{j}(\vec{x})\partial_{x_{i}}W_{j}(x_{j}).$
(34)

The Jacobian for the dynamical system without window functions can be
evaluated exactly, and is given by (using
$\partial_{\alpha}A^{-1}=-A^{-1}(\partial_{\alpha}A)A^{-1}$) by
($\delta_{ij})$ is the Kronecker delta):

$\displaystyle DJ_{ji}$ $\displaystyle=$
$\displaystyle\partial_{x_{i}}f_{j}(\vec{x})$ (35) $\displaystyle=$
$\displaystyle\frac{\chi}{\beta\alpha}\sum_{krpm}(I-\chi\Omega
X)^{-1}_{jk}\Omega_{kr}\delta_{ri}(I-\chi\Omega
X)^{-1}_{rm}(\Omega\vec{S})_{m}-\delta_{ji}$ $\displaystyle=$
$\displaystyle\frac{\chi}{\alpha\beta}P_{ji}V_{i}-\delta_{ji},$

where $P=(I-\chi\Omega X)^{-1}\Omega$ and $\vec{V}=P\vec{S}$, and which we
could use to evaluate the Local Lyapunov Exponents or study local stability.
The local flow divergence is obtained via the relation

$\displaystyle\delta\vec{x}(t_{n+1})=\big{(}I+DJ^{f}(t_{n})dt\big{)}\delta\vec{x}(t_{n}).$
(36)

The spectrum of the Jacobian at the fixed point has been studied first in
[44], showing that the spectrum is real even if the matrix is not symmetric.
Despite this, non-normal Jacobian matrices can in principle exhibit
amplification of perturbations on a stable equilibrium and it is worth of
careful scrutiny in this paper. Small perturbations of a stable equilibria
typically decay in exponential time. In the case of non-normal matrices
however, even if the eigenvalues have all real part negative, can exhibit a
transient instability. How non-normal a matrix $M$ is quantified by (any) norm
of the matrix $N=(M^{t}M-MM^{t})/2$. The Jacobian matrix of our system can be
written in the form $DJ(X)=P^{2}\mathcal{S}$, where $P=(I-\chi\Omega
X)^{-1}\Omega$ and $\mathcal{S}_{ij}=S_{i}\delta_{ij}$. It follows that
$N=\frac{P^{2}\mathcal{S}-\mathcal{S}(P^{t})^{2}}{2}$. This is not
inconsistent with the notion of stability. As perturbations increase in
magnitude via a transient phase, the effect on eigenvalues is more pronounced.
For long times, if the system is linear, the long term dynamics is still
dominated by the negativity of the real part of the spectrum. However, in the
case of non-normal matrices, the transient dynamics can be no longer
consistent with the traditional picture of asymptotic stability.

For this reason we go beyond a simple spectral analysis. Generalizations of
the notion of eigenvalues have been studied in the literature. In general, the
eigenvalues of ${DJ}$ can be defined:
$\Lambda({DJ})=\\{z\in\mathbb{C}:\text{det}(zI-DJ)=0\\}$, meaning that if $z$
is an eigenvalue of ${DJ}$ then by convention the norm of $(zI-DJ)^{-1}$ is
defined to be infinity [52]. But if $||(zI-DJ)^{-1}||$ is finite and very
large, as is often the case with perturbed non-normal matrices, then the
pseudospectrum of ${DJ}$ must be considered. The ‘$\epsilon$-pseudospectrum’
is a generalization of the notion of eigenvalues which depend on a real
parameter $\epsilon$, and can be defined in various equivalent ways [52]. We
use the following definition:

$\Lambda_{\epsilon}({DJ})=\\{z\in\mathbb{C}:||(zI-{DJ})^{-1}||\geq\epsilon^{-1}\\}.$
(37)

If a matrix is normal then its $\epsilon$-pseudospectrum (from now one we only
call it ‘pseudospectrum’) is somewhat trivial: it consists of closed balls of
radius $\epsilon$ surrounding the original eigenvalues of $DJ$. In the case of
non-normal matrices however, pseudospectra can be much larger and much more
intricate, which is the case we are interested in here.

Local asymptotic stability is determined in the same way for normal and non-
normal matrices. The ‘spectral abscissa’ of ${DJ}$ is defined as
$\alpha({DJ})=\sup_{z\in\Lambda({DJ})}\text{Re}(z)$, where the supremum is the
largest real part of $\Lambda({DJ})$; clearly, stability is guaranteed for
$\alpha({DJ})<0$. If ${DJ}$ is normal, then $||e^{{DJ}t}||=e^{\alpha({DJ})t}$
and dynamics is determined by $\alpha(DJ)$. The relationship between the
spectral abscissa and the system dynamics is evident from the bouns [52]:
$e^{\alpha({DJ})t}\leq||e^{{DJ}t}||\leq\kappa({\vec{V}})e^{\alpha({DJ})t}$
where $\kappa({\bf V})=||{\bf V}||\cdot||{\bf V}^{-1}||$ is called the
conditioning of the matrix ${\bf V}$, which is built from the eigenvectors of
${DJ}$, and is a measure of invertibility of ${\bf V}$. We see that for a
normal dynamics, the conditioning provides an upper bound to the maximum
amplification of the perturbations. Clearly, if the spectral abscissa is
positive, perturbations are unstable but within the bounds defined above.
These bounds can however generalized by introducing the notion of
‘$\epsilon$-pseudospectral abscissa’ of ${DJ}$. This is defined as
$\alpha_{\epsilon}({DJ})=\sup_{z\in\Lambda_{\epsilon}({DJ})}\text{Re}(z)$,
e.g. the largest real part of the spectrum of $DJ$ (for a given $\epsilon$).
The relevance of the $\epsilon$-pseudospectral abscissa is given by the fact
that it provides a lower bound to the maximum amplification of the
perturbation[52]:

$\sup_{\epsilon\geq 0}\frac{\alpha_{\epsilon}({DJ})}{\epsilon}\leq\sup_{t\geq
0}||e^{{DJ}t}||.$ (38)

The quantity

$f_{DJ}(\epsilon)=\frac{\alpha_{\epsilon}({DJ})}{\epsilon}$ (39)

is called the Kreiss constant of the matrix $DJ$, and is thus relevant in
order to understand the transient dynamics for short times. Given the bound of
Eqn. (38), it is important to look for values $\epsilon^{*}$, at which
$f_{DJ}(\epsilon^{*})=1$. Via a visual analysis of the contours of
pseudospectrum, a necessary condition of a critical value of $\epsilon^{*}$ is
that $\epsilon^{*}$-contour crosses the imaginary axis. At this point,
perturbations from the equilibrium population vector can be in principle
amplified, depending on the value of the Kreiss constant.

In order to study the non-normality, we looked at the Jacobian matrix
evaluated at the equilibrium points for the dynamics obtained numerically via
Monte Carlo. In Fig. 7 (left) we can observe the bulk of
$\epsilon$-pseudospectrum of the Jacobian matrix $DJ$, fixed $\Omega$ as in
the main text, and with $\alpha=\beta=1$ and $\chi=0.9$. However, such non-
normality is not enough to explain the transient instability we observe. A
careful numerical analysis of the Kreiss constant for varying $s$ gave us a
maximum value of the Kreiss matrix of $K=sup_{\epsilon}F_{DK}(\epsilon)=0.7$,
thus less than 1, the critical value for the amplification. We thus directly
studied the norm of $e^{DJt}$ in Fig. 7 (right), not observing any
amplification phenomenon up to the instability threshold for the matrix. This
implies that $\sup_{\epsilon}\frac{\alpha_{\epsilon}(DJ)}{\epsilon}\leq 1$ and
that transient instability is not the cause of the escape from the minimum.

Figure 7: Left: $\epsilon$-Pseudospectrum of the Jacobian matrix for
$\epsilon=[10^{-2},10^{-1}]$. Right: behavior of the norm $||e^{DJ\
t}||_{\infty}$ as a function of time for $s\in[0.24,0.27]$.

## E: Basins of attraction

First of all, let us first mention, as shown in the main text, that the
tunneling probability depends on the number of dynamical variables. In Fig. 6
we see that the larger the height of the barrier ($\Delta E$ decreases with
$s$), the larger the number of variables one needs to tunnel through. In order
to obtain a glimpse of the structure of the basins of attractions and the role
of the dimensionality in reaching the two asymptotic states we perform the
following numerical experiment. Given a random instance of $\Omega$ obtained
as $\Omega=A^{t}(AA^{t})^{-1}A$, for $A$ a random matrix (with flat
probability on $[0,1]$) of size $N_{c}\times N$ ($N=200$, $N_{c}=100$) we fix
$\chi$ and $s=\Omega\vec{S}$ in a region in which we know there is a mixing
between the local minimum and global minimum of the mean field potential
potential $V(x,\chi)$ in $x\in[0,1]$. Specifically, we study the basins of
attraction of the system near the mean field situation (e.g. a single
memristor), e.g. when many of the variables $x_{i}$ are initialized in the
same value. We consider this analysis in the regime reported in the main text,
for the values $\chi=0.9$ and $s\approx 0.22$ which is close to the boundary
between the two minima. We then perform the following two experiments. In the
first, we choose $x_{i}(0)$ randomly, and then follow the trajectory and mark
it as local or global minimum asymptotic state. Since $\Omega_{ij}$ is random
and there is no preference between the variables, we choose to plot initial
conditions of the variables $x_{1}$ and $x_{2}$ and mark them depending on
which asymptotic state they reached. In the second experiment instead we
choose $x_{3}\cdots x_{N}=\tilde{x}$, and $x_{1}$ and $x_{2}$ still at random,
following the procedure of the first experiment and following the trajectory.
The value of $\tilde{x}$ is an extra parameter and has to be obtained such
that system ends in both minima. A value for $\tilde{x}$ which we found to be
working well for random $\Omega$ is $\tilde{x}\approx.5135$. The results are
shown in Fig. 9, in which we see a complete mixing in the first experiment
between the asymptotic states, and a complete separation via a separatrix in
the second experiment. The separatrix moves with $\tilde{x}$, and for larger
values the system it tends to move towards the global minimum; for smaller
values towards the local minimum. Intuitively, this implies that for random
initial conditions, for every possible asymptotic state there is a initial
condition nearby to the global minimum.

Figure 8: Distributions $P(x^{1}_{cg})$ and $P(x^{2}_{cg})$ as a function of
$\Delta E$ of the potential. Inset: Fit of the Fermi distribution with a
chemical potential.

A more careful analysis has been done for random initial conditions. We
studied via Monte Carlo sampling the probability distribution, given a random
initial condition, of the system being in one or the other local minimum,
showing that it can be fit by a Fermi-like distribution depending on the
energy barrier $\Delta E$. Let $P(x;s,\chi)$ be the probability distribution
as a function of the height of the barrier $\Delta E(s)=E_{max}-E_{local\
min}$ from the local minimum. We obtained such function numerically, averaging
over 100 Monte Carlo simulations for each value of $s$, and shown in Fig. 8.
We fit $P(\Delta E)$ with a function of the form $P(x;s,\chi)=1/(1-e^{-{\Delta
E(s,\chi)-\mu}/{T_{eff}}})$, where we obtained numerically the value
$T_{eff}=0.0048$, and $\mu=0.036$. These results are shown in Fig. 8. We see
that the transition probabilities of Fig. 8 are consistent with the sigmoidal
probabilities described in the main text.

Figure 9: Projected basin of attraction for the coarse grained asymptotic
states according to the two controlled numerical experiments, as described in
the text, in order to visualize the basins of attractions. On the left, we see
the completely randomize initial condition, in which we see a complete mixing
in the initial conditions for $x_{1}$ and $x_{2}$ as a function of the
asymptotic states. On the right, we see that there is a neat separation
between the two asymptotic states via a separatrix bordering two basins of
attractions.


$|\tau_{{\mathrm{2}}}|\,=\,|\tau^{\prime}|<|\tau^{\prime}\operatorname{\mathbf{ref}}^{1}|$,
$\mathit{a}$ cannot be reachable from $\mathit{R}(\mathit{x})$ in $\mathit{H}$
(otherwise we would have $\mathit{a}$ reaches an integer along multiple heap
paths of differing lengths, a clear contradiction). Then the value of
$\mathit{a}$ in $\mathit{H}$ is irrelevant to the derivation of
$\mathit{H}\vdash\mathit{R}(\mathit{x})\Downarrow|\tau_{{\mathrm{2}}}|$,
whereby
$\mathit{H}\\{\mathit{a}\hookleftarrow\mathit{R}(\mathit{x})\\}\vdash\mathit{R}(\mathit{x})\Downarrow|\tau_{{\mathrm{2}}}|$
must hold.

Then, it’s to show
$\forall\,\mathit{a^{\prime}}\in\mathit{dom}(\mathit{H^{\prime}}).\mathbf{Own}\,(\mathit{H^{\prime}},\mathit{R},\Gamma^{\prime\prime})(\mathit{a^{\prime}})=(\mathit{O^{\prime}_{{\mathrm{1}}}}+\mathit{O^{\prime\prime}_{{\mathrm{1}}}})(\mathit{a^{\prime}})=\mathit{O^{\prime}_{{\mathrm{1}}}}(\mathit{a^{\prime}})+\mathit{O^{\prime\prime}_{{\mathrm{1}}}}(\mathit{a^{\prime}})\,\leq\,1$.
For every
$\mathit{z}\in\mathit{dom}(\Gamma)\setminus\\{\mathit{y},\mathit{x}\\}$ (and
similarly for $\Gamma^{\prime\prime}$), we have
$\Gamma(\mathit{z})=\Gamma^{\prime\prime}(\mathit{z})$. Further, from
$\mathit{O^{\prime}_{{\mathrm{0}}}}(\mathit{a})=0$ above, we must have
$\mathbf{own}\,(\mathit{H},\mathit{R}(\mathit{z}),\Gamma(\mathit{z}))(\mathit{a})=0$
for all such $\mathit{z}$. As
$\mathit{H}\approx_{\mathit{a}}\mathit{H^{\prime}}$, by Lemma 19, we have that
$\mathit{O^{\prime}_{{\mathrm{0}}}}=\mathit{O^{\prime}_{{\mathrm{1}}}}$. Then,
from
$\forall\,\mathit{a^{\prime}}\in\mathit{dom}(\mathit{H}).\mathbf{Own}\,(\mathit{H},\mathit{R},\Gamma)(\mathit{a^{\prime}})=\mathit{O^{\prime}_{{\mathrm{0}}}}(\mathit{a^{\prime}})+\mathit{O^{\prime\prime}_{{\mathrm{0}}}}(\mathit{a^{\prime}})\,\leq\,1$,
it suffices to show that
$\forall\,\mathit{a^{\prime}}\in\mathit{dom}(\mathit{H}).\mathit{O^{\prime\prime}_{{\mathrm{1}}}}(\mathit{a^{\prime}})\leq\mathit{O^{\prime\prime}_{{\mathrm{0}}}}(\mathit{a^{\prime}})$.

We first consider the case for $\mathit{a}$:

$\displaystyle\mathit{O^{\prime\prime}_{{\mathrm{1}}}}(\mathit{a})$
$\displaystyle=\mathbf{own}\,(\mathit{H^{\prime}},\mathit{H^{\prime}}(\mathit{R}(\mathit{y})),\tau_{{\mathrm{2}}})(\mathit{a})+\mathbf{own}\,(\mathit{H^{\prime}},\mathit{R}(\mathit{x}),\tau_{{\mathrm{1}}})(\mathit{a})+1$
$\displaystyle\mathit{O^{\prime\prime}_{{\mathrm{0}}}}(\mathit{a})$
$\displaystyle=\mathbf{own}\,(\mathit{H},\mathit{R}(\mathit{x}),\tau_{{\mathrm{1}}}+\tau_{{\mathrm{2}}})(\mathit{a})+\mathbf{own}\,(\mathit{H},\mathit{H}(\mathit{R}(\mathit{y})),\tau^{\prime})(\mathit{a})+1$

From above, we have
$\mathbf{own}\,(\mathit{H},\mathit{R}(\mathit{x}),\tau_{{\mathrm{2}}}+\tau_{{\mathrm{1}}})(\mathit{a})=\mathbf{own}\,(\mathit{H},\mathit{H}(\mathit{R}(\mathit{y})),\tau^{\prime})(\mathit{a})=0$.
By Lemma 19 and $\mathit{H}\approx_{\mathit{a}}\mathit{H^{\prime}}$, we have
$\mathbf{own}\,(\mathit{H},\mathit{R}(\mathit{x}),\tau_{{\mathrm{2}}}+\tau_{{\mathrm{1}}})=\mathbf{own}\,(\mathit{H^{\prime}},\mathit{R}(\mathit{x}),\tau_{{\mathrm{2}}}+\tau_{{\mathrm{1}}})$.
Also by Lemma 15, we have
$\mathbf{own}\,(\mathit{H^{\prime}},\mathit{R}(\mathit{x}),\tau_{{\mathrm{2}}}+\tau_{{\mathrm{1}}})=\mathbf{own}\,(\mathit{H^{\prime}},\mathit{R}(\mathit{x}),\tau_{{\mathrm{1}}})+\mathbf{own}\,(\mathit{H^{\prime}},\mathit{R}(\mathit{x}),\tau_{{\mathrm{2}}})$.
From $\mathit{H^{\prime}}(\mathit{R}(\mathit{y}))\,=\,\mathit{R}(\mathit{x})$,
we therefore have
$\mathbf{own}\,(\mathit{H^{\prime}},\mathit{H^{\prime}}(\mathit{R}(\mathit{y})),\tau_{{\mathrm{2}}})(\mathit{a})=\mathbf{own}\,(\mathit{H^{\prime}},\mathit{R}(\mathit{x}),\tau_{{\mathrm{1}}})(\mathit{a})=0$,
and thus:

$\mathit{O^{\prime\prime}_{{\mathrm{1}}}}(\mathit{a})=\mathbf{own}\,(\mathit{H^{\prime}},\mathit{H^{\prime}}(\mathit{R}(\mathit{y})),\tau_{{\mathrm{2}}})(\mathit{a})+\mathbf{own}\,(\mathit{H^{\prime}},\mathit{R}(\mathit{x}),\tau_{{\mathrm{1}}})(\mathit{a})+1=1=\mathit{O^{\prime\prime}_{{\mathrm{0}}}}(\mathit{a})$

Next, consider some $\mathit{a}\,\neq\,\mathit{a^{\prime}}$;

$\displaystyle\mathit{O^{\prime\prime}_{{\mathrm{1}}}}(\mathit{a^{\prime}})$
$\displaystyle=\mathbf{own}\,(\mathit{H^{\prime}},\mathit{H^{\prime}}(\mathit{R}(\mathit{y})),\tau_{{\mathrm{2}}})(\mathit{a^{\prime}})+\mathbf{own}\,(\mathit{H^{\prime}},\mathit{R}(\mathit{x}),\tau_{{\mathrm{1}}})(\mathit{a^{\prime}})$
$\displaystyle\mathit{O^{\prime\prime}_{{\mathrm{0}}}}(\mathit{a^{\prime}})$
$\displaystyle=\mathbf{own}\,(\mathit{H},\mathit{R}(\mathit{x}),\tau_{{\mathrm{2}}}+\tau_{{\mathrm{1}}})(\mathit{a^{\prime}})+\mathbf{own}\,(\mathit{H},\mathit{H}(\mathit{R}(\mathit{y})),\tau^{\prime})(\mathit{a^{\prime}})$

By reasoning similar to the case for $a=a^{\prime}$, we have
$\mathit{O^{\prime\prime}_{{\mathrm{1}}}}(\mathit{a^{\prime}})\leq\mathbf{own}\,(\mathit{H},\mathit{R}(\mathit{x}),\tau_{{\mathrm{2}}}+\tau_{{\mathrm{1}}})(\mathit{a^{\prime}})\leq\mathit{O^{\prime\prime}_{{\mathrm{0}}}}(\mathit{a^{\prime}})$.
We therefore conclude that
$\forall\,\mathit{a^{\prime}}\in\mathit{dom}(\mathit{H^{\prime}}).\mathbf{Own}\,(\mathit{H^{\prime}},\mathit{R},\Gamma^{\prime\prime})(\mathit{a^{\prime}})\leq
1$.

Case R-Alias:
$\quad\begin{array}[t]{@{}l@{\qquad}l@{\qquad}l@{\qquad}l}\vdash_{\mathit{conf}}^{D}\left\langle{\mathit{H},\mathit{R},\vec{F},\mathit{E}[\mathbf{alias}(\mathit{x}=\mathit{y})\operatorname{;}\mathit{e}]}\right\rangle\\\
\left\langle{\mathit{H},\mathit{R},\vec{F},\mathit{E}[\mathbf{alias}(\mathit{x}=\mathit{y})\operatorname{;}\mathit{e}]}\right\rangle\longrightarrow_{\mathit{D}}\left\langle{\mathit{H},\mathit{R},\vec{F},\mathit{E}[\mathit{e}]}\right\rangle\\\
\mathit{R}(\mathit{x})\,=\,\mathit{R}(\mathit{y})\end{array}$
By inversion (see R-Deref) we have for some $\Gamma$ that:

$\displaystyle\Theta\mid\vec{\ell}\mid\Gamma[\mathit{x}:\tau_{{\mathrm{1}}}\operatorname{\mathbf{ref}}^{r_{{\mathrm{1}}}}][\mathit{y}:\tau_{{\mathrm{2}}}\operatorname{\mathbf{ref}}^{r_{{\mathrm{2}}}}]\vdash\mathbf{alias}(\mathit{x}=\mathit{y})\operatorname{;}\mathit{e}:\tau\operatorname{\Rightarrow}\Gamma^{\prime}$
$\displaystyle\Theta\mid\vec{\ell}\mid\Gamma[\mathit{x}\hookleftarrow\tau^{\prime}_{{\mathrm{1}}}\operatorname{\mathbf{ref}}^{r^{\prime}_{{\mathrm{1}}}}][\mathit{y}\hookleftarrow\tau^{\prime}_{{\mathrm{2}}}\operatorname{\mathbf{ref}}^{r^{\prime}_{{\mathrm{2}}}}]\vdash\mathit{e}:\tau\operatorname{\Rightarrow}\Gamma^{\prime}$
$\displaystyle\tau_{{\mathrm{1}}}\operatorname{\mathbf{ref}}^{r_{{\mathrm{1}}}}+\tau_{{\mathrm{2}}}\operatorname{\mathbf{ref}}^{r_{{\mathrm{2}}}}\approx\tau^{\prime}_{{\mathrm{1}}}\operatorname{\mathbf{ref}}^{r^{\prime}_{{\mathrm{1}}}}+\tau^{\prime}_{{\mathrm{2}}}\operatorname{\mathbf{ref}}^{r^{\prime}_{{\mathrm{2}}}}$
$\displaystyle\mathbf{Cons}\,(\mathit{H},\mathit{R},\Gamma)$

We give
$\Gamma^{\prime\prime}=\Gamma[\mathit{x}\hookleftarrow\tau^{\prime}_{{\mathrm{1}}}\operatorname{\mathbf{ref}}^{r^{\prime}_{{\mathrm{1}}}}][\mathit{y}\hookleftarrow\tau^{\prime}_{{\mathrm{2}}}\operatorname{\mathbf{ref}}^{r^{\prime}_{{\mathrm{2}}}}]$,
and must show
$\Theta\mid\vec{\ell}\mid\Gamma^{\prime\prime}\vdash\mathit{e}:\tau\operatorname{\Rightarrow}\Gamma^{\prime}$
and $\mathbf{Cons}\,(\mathit{H},\mathit{R},\Gamma^{\prime\prime})$. The first
is immediate.

To show $\mathbf{Cons}\,(\mathit{H},\mathit{R},\Gamma^{\prime\prime})$ we
first define:

$\displaystyle\tau_{\mathit{p}\,{\mathrm{1}}}=$
$\displaystyle\tau_{{\mathrm{1}}}\operatorname{\mathbf{ref}}^{r_{{\mathrm{1}}}}$
$\displaystyle\tau_{\mathit{p}\,{\mathrm{2}}}=$
$\displaystyle\tau_{{\mathrm{2}}}\operatorname{\mathbf{ref}}^{r_{{\mathrm{2}}}}$
$\displaystyle\tau_{\mathit{q}\,{\mathrm{1}}}=$
$\displaystyle\tau^{\prime}_{{\mathrm{1}}}\operatorname{\mathbf{ref}}^{r^{\prime}_{{\mathrm{1}}}}$
$\displaystyle\tau_{\mathit{q}\,{\mathrm{2}}}=$
$\displaystyle\tau^{\prime}_{{\mathrm{2}}}\operatorname{\mathbf{ref}}^{r^{\prime}_{{\mathrm{2}}}}$
$\displaystyle\tau_{\mathit{q}}=$
$\displaystyle\tau_{\mathit{q}\,{\mathrm{1}}}+\tau_{\mathit{q}\,{\mathrm{2}}}$
$\displaystyle\tau_{\mathit{p}}=$
$\displaystyle\tau_{\mathit{p}\,{\mathrm{1}}}+\tau_{\mathit{p}\,{\mathrm{2}}}$

We thus have $\tau_{\mathit{q}}\approx\tau_{\mathit{p}}$.

We know that $\mathbf{Cons}\,(\mathit{H},\mathit{R},\Gamma)$, from which we
have $\mathbf{SAT}\,(\mathit{H},\mathit{R},\Gamma)$, in particular
$\mathbf{SATv}(\mathit{H},\mathit{R},\mathit{R}(\mathit{y}),\Gamma(\mathit{y})=\tau_{{\mathrm{2}}}\operatorname{\mathbf{ref}}^{r_{{\mathrm{2}}}})$
and
$\mathbf{SATv}(\mathit{H},\mathit{R},\mathit{R}(\mathit{x}),\Gamma(\mathit{x})=\tau_{{\mathrm{1}}}\operatorname{\mathbf{ref}}^{r_{{\mathrm{1}}}})$.
From
$\tau_{\mathit{p}\,{\mathrm{1}}}+\tau_{\mathit{p}\,{\mathrm{2}}}=\tau_{\mathit{p}}$
and Lemma 16, we have
$\mathbf{SATv}(\mathit{H},\mathit{R},\mathit{v},\tau_{\mathit{p}\,{\mathrm{1}}})$
and
$\mathbf{SATv}(\mathit{H},\mathit{R},\mathit{v},\tau_{\mathit{p}\,{\mathrm{2}}})$
imply $\mathbf{SATv}(\mathit{H},\mathit{R},\mathit{v},\tau_{\mathit{p}})$,
where
$\mathit{v}=\mathit{H}(\mathit{R}(\mathit{y}))=\mathit{H}(\mathit{R}(\mathit{x}))$.
From $\tau_{\mathit{q}}\approx\tau_{\mathit{p}}$ and Lemma 13, we have that
$\mathbf{SATv}(\mathit{H},\mathit{R},\mathit{v},\tau_{\mathit{p}})$ implies
$\mathbf{SATv}(\mathit{H},\mathit{R},\mathit{v},\tau_{\mathit{q}})$. From
Lemma 16 we also have that
$\mathbf{SATv}(\mathit{H},\mathit{R},\mathit{v},\tau_{\mathit{q}})$ implies
$\mathbf{SATv}(\mathit{H},\mathit{R},\mathit{v},\tau_{\mathit{q}\,{\mathrm{1}}})$
and
$\mathbf{SATv}(\mathit{H},\mathit{R},\mathit{v},\tau_{\mathit{q}\,{\mathrm{2}}})$,
where again
$\mathit{v}=\mathit{H}(\mathit{R}(\mathit{y}))=\mathit{H}(\mathit{R}(\mathit{x}))$.

Then from the reasoning above, the refinements of
$\tau_{\mathit{q}\,{\mathrm{1}}}$ and $\tau_{\mathit{q}\,{\mathrm{2}}}$ are
valid and $\mathbf{Cons}\,(\mathit{H},\mathit{R},\Gamma^{\prime\prime})$
holds.

Then, it’s to show
$\forall\,\mathit{a}\in\mathit{dom}(\mathit{H}).\mathbf{Own}\,(\mathit{H},\mathit{R},\Gamma^{\prime\prime})(\mathit{a})\leq
1$. To prove that
$\mathbf{Own}\,(\mathit{H},\mathit{R},\Gamma)=\mathbf{Own}\,(\mathit{H},\mathit{R},\Gamma^{\prime\prime})$
follows from:

$\displaystyle\mathbf{own}\,(\mathit{H},\mathit{R}(\mathit{x}),\tau_{{\mathrm{1}}}\operatorname{\mathbf{ref}}^{r_{{\mathrm{1}}}})+\mathbf{own}\,(\mathit{H},\mathit{R}(\mathit{y}),\tau_{{\mathrm{2}}}\operatorname{\mathbf{ref}}^{r_{{\mathrm{2}}}})=$
$\displaystyle\,\,\,\,\mathbf{own}\,(\mathit{H},\mathit{R}(\mathit{x}),\tau^{\prime}_{{\mathrm{1}}}\operatorname{\mathbf{ref}}^{r^{\prime}_{{\mathrm{1}}}})+\mathbf{own}\,(\mathit{H},\mathit{R}(\mathit{y}),\tau^{\prime}_{{\mathrm{2}}}\operatorname{\mathbf{ref}}^{r^{\prime}_{{\mathrm{2}}}})$

which follows immediately from the conditions
$\tau_{{\mathrm{1}}}\operatorname{\mathbf{ref}}^{r_{{\mathrm{1}}}}+\tau_{{\mathrm{2}}}\operatorname{\mathbf{ref}}^{r_{{\mathrm{2}}}}\approx\tau^{\prime}_{{\mathrm{1}}}\operatorname{\mathbf{ref}}^{r^{\prime}_{{\mathrm{1}}}}+\tau^{\prime}_{{\mathrm{2}}}\operatorname{\mathbf{ref}}^{r^{\prime}_{{\mathrm{2}}}}$,
$\mathit{R}(\mathit{x})\,=\,\mathit{R}(\mathit{y})$, and Lemmas 15 and 14.

Case R-AliasPtr:
By reasoning similar to the R-Alias case.

Case R-AliasFail,R-AliasPtrFail:
The result configuration $\mathbf{AliasFail}$ is trivially well-typed.

Case R-Assert:
$\quad\begin{array}[t]{@{}l@{\qquad}l@{\qquad}l@{\qquad}l}\vdash_{\mathit{conf}}^{D}\left\langle{\mathit{H},\mathit{R},\vec{F},\mathit{E}[\mathbf{assert}(\varphi)\operatorname{;}\mathit{e}]}\right\rangle\quad\quad\Gamma\models[\mathit{R}]\,\varphi\\\
\left\langle{\mathit{H},\mathit{R},\vec{F},\mathit{E}[\mathbf{assert}(\varphi)\operatorname{;}\mathit{e}]}\right\rangle\longrightarrow_{\mathit{D}}\left\langle{\mathit{H},\mathit{R},\vec{F},\mathit{E}[\mathit{e}]}\right\rangle\\\
\end{array}$
By inversion (see R-Deref) we can obtain
$\Theta\mid\vec{\ell}\mid\Gamma\vdash\mathbf{assert}(\varphi)\operatorname{;}\mathit{e}:\tau\operatorname{\Rightarrow}\Gamma^{\prime}$
and
$\Theta\mid\vec{\ell}\mid\Gamma\vdash\mathit{e}:\tau\operatorname{\Rightarrow}\Gamma^{\prime}$,
and the result follows immediately by taking $\Gamma^{\prime\prime}=\Gamma$.

Case R-AssertFail:
$\quad\begin{array}[t]{@{}l@{\qquad}l@{\qquad}l@{\qquad}l}\vdash_{\mathit{conf}}^{D}\left\langle{\mathit{H},\mathit{R},\vec{F},\mathit{E}[\mathbf{assert}(\varphi)\operatorname{;}\mathit{e}]}\right\rangle\\\
\left\langle{\mathit{H},\mathit{R},\vec{F},\mathit{E}[\mathbf{assert}(\varphi)\operatorname{;}\mathit{e}]}\right\rangle\longrightarrow_{\mathit{D}}\mathbf{AssertFail}\\\
\Theta\mid\vec{\ell}\mid\Gamma\vdash\mathbf{assert}(\varphi)\operatorname{;}\mathit{e}:\tau\operatorname{\Rightarrow}\Gamma^{\prime}\end{array}$
By inversion (see the R-Deref case) we have that $\Gamma\models\varphi$, i.e.,
$\models\left\llbracket{\Gamma}\right\rrbracket\implies\varphi$, for some
$\Gamma$ such that $\mathbf{Cons}\,(\mathit{H},\mathit{R},\Gamma)$. From Lemma
11 we therefore have
$\models[\mathit{R}]\,\left\llbracket{\Gamma}\right\rrbracket$. From the
precondition of R-AssertFail we have that $\not\models[\mathit{R}]\,\varphi$.
But from $\models\left\llbracket{\Gamma}\right\rrbracket\implies\varphi$ and
$\models[\mathit{R}]\,\left\llbracket{\Gamma}\right\rrbracket$ we can conclude
that $\models[\mathit{R}]\,\varphi$, yielding a contradiction. We therefore
conclude that this case is impossible.

Case R-Call:
$\quad\begin{array}[t]{@{}l@{\qquad}l@{\qquad}l@{\qquad}l}\vdash_{\mathit{conf}}^{D}\left\langle{\mathit{H},\mathit{R},\vec{F},\mathit{E}[\operatorname{\mathbf{let}}\mathit{x}=\mathit{f}^{\ell}(\mathit{y_{{\mathrm{1}}}},\ldots,\mathit{y_{\mathit{n}}})\operatorname{\mathbf{in}}\mathit{e^{\prime}}]}\right\rangle\\\
\mathit{f}\mapsto(\mathit{x_{{\mathrm{1}}}},\,..\,,\mathit{x_{\mathit{n}}})\mathit{e}\in\mathit{D}\\\
\begin{array}[]{l}\left\langle{\mathit{H},\mathit{R},\vec{F},\mathit{E}[\operatorname{\mathbf{let}}\mathit{x}=\mathit{f}^{\ell}(\mathit{y_{{\mathrm{1}}}},\ldots,\mathit{y_{\mathit{n}}})\operatorname{\mathbf{in}}\mathit{e^{\prime}}]}\right\rangle\\\
\quad\longrightarrow_{\mathit{D}}\left\langle{\mathit{H},\mathit{R},\mathit{E}[\operatorname{\mathbf{let}}\mathit{x}=[]^{\ell}\operatorname{\mathbf{in}}\mathit{e^{\prime}}]:\vec{F},[\mathit{y_{{\mathrm{1}}}}/\mathit{x_{{\mathrm{1}}}}]\cdots[\mathit{y_{\mathit{n}}}/\mathit{x_{\mathit{n}}}]\mathit{e}}\right\rangle\end{array}\\\
\end{array}$
We must show that
$\vdash_{\mathit{conf}}^{D}\left\langle{\mathit{H},\mathit{R},\mathit{E}[\operatorname{\mathbf{let}}\mathit{x}=[]^{\ell}\operatorname{\mathbf{in}}\mathit{e^{\prime}}]:\vec{F},[\mathit{y_{{\mathrm{1}}}}/\mathit{x_{{\mathrm{1}}}}]\cdots[\mathit{y_{\mathit{n}}}/\mathit{x_{\mathit{n}}}]\mathit{e}}\right\rangle$
for some $\Gamma^{\prime\prime}$.

By inversion on the configuration typing, we have that, for some $\Gamma$:

$\Theta\mid\vec{\ell}\mid\Gamma\vdash\mathit{E}[\operatorname{\mathbf{let}}\mathit{x}=\mathit{f}^{\ell}(\mathit{y_{{\mathrm{1}}}},\ldots,\mathit{y_{\mathit{n}}})\operatorname{\mathbf{in}}\mathit{e^{\prime}}]:\tau_{\mathit{n}}\operatorname{\Rightarrow}\Gamma_{\mathit{n}}.$

By Lemma 7, we then have for some $\tau$, and $\Gamma^{\prime}$ that:

$\displaystyle\Theta\mid\vec{\ell}\mid\Gamma\vdash\operatorname{\mathbf{let}}\mathit{x}=\mathit{f}^{\ell}(\mathit{y_{{\mathrm{1}}}},\ldots,\mathit{y_{\mathit{n}}})\operatorname{\mathbf{in}}\mathit{e^{\prime}}:\tau\operatorname{\Rightarrow}\Gamma^{\prime}$
$\displaystyle\Theta\mid[]:\tau\operatorname{\Rightarrow}\Gamma^{\prime}\mid\vec{\ell}\vdash_{\mathit{ectx}}\mathit{E}:\tau_{\mathit{n}}\operatorname{\Rightarrow}\Gamma_{\mathit{n}}$

Taking
$\tau_{{\mathrm{1}}}=\tau,\Gamma_{{\mathrm{1}}}=\Gamma^{\prime},\Gamma_{{\mathrm{0}}}=\Gamma,\Gamma_{{\mathrm{2}}}=\Gamma_{\mathit{n}},\tau_{{\mathrm{2}}}=\tau_{\mathit{n}}$,
by Lemma 28 we have, for some
$\tau^{\prime\prime\prime},\Gamma^{\prime\prime\prime}$:

$\displaystyle\Theta\mid\ell:\vec{\ell}\mid\Gamma\vdash\sigma_{x}\mathit{e}:\tau^{\prime\prime\prime}\operatorname{\Rightarrow}\Gamma^{\prime\prime\prime}$
$\displaystyle\Theta\mid[]:\tau^{\prime\prime\prime}\operatorname{\Rightarrow}\Gamma^{\prime\prime\prime}\mid\vec{\ell}\vdash_{\mathit{ectx}}\mathit{E}[\operatorname{\mathbf{let}}\mathit{x}=[]^{\ell}\operatorname{\mathbf{in}}\mathit{e^{\prime}}]:\tau_{\mathit{n}}\operatorname{\Rightarrow}\Gamma_{\mathit{n}}$

where:

$\displaystyle\sigma_{x}$
$\displaystyle=[\mathit{y_{{\mathrm{1}}}}/\mathit{x_{{\mathrm{1}}}}]\cdots[\mathit{y_{\mathit{n}}}/\mathit{x_{\mathit{n}}}]$
$\displaystyle\Theta(\mathit{f})$
$\displaystyle=\forall\lambda.\left\langle{\mathit{x_{{\mathrm{1}}}}\operatorname{:}\tau_{\mathit{i}},\dots,\mathit{x_{\mathit{n}}}\operatorname{:}\tau_{\mathit{n}}}\right\rangle\operatorname{\rightarrow}\left\langle{\mathit{x_{{\mathrm{1}}}}\operatorname{:}\tau^{\prime}_{{\mathrm{1}}},\dots,\mathit{x_{\mathit{n}}}\operatorname{:}\tau^{\prime}_{\mathit{n}}\mid\tau_{\mathit{p}}}\right\rangle$

We therefore take $\Gamma^{\prime\prime}=\Gamma$.

We must also prove that $\forall
i\in\left\\{{1..n+1}\right\\}.\Theta\mid[]:\tau_{\mathit{i}}\operatorname{\Rightarrow}\Gamma_{\mathit{i}}\mid\vec{\ell}_{{\mathit{i}-1}}\vdash_{\mathit{ectx}}\mathit{E^{\prime}_{{\mathit{i}-1}}}:\tau_{{\mathit{i}-1}}\operatorname{\Rightarrow}\Gamma_{{\mathit{i}-1}}$
where
$\mathit{E^{\prime}_{\mathit{n}}}=\mathit{E}[\operatorname{\mathbf{let}}\mathit{x}=[]^{\ell}\operatorname{\mathbf{in}}\mathit{e^{\prime}}]$
and $\mathit{E^{\prime}_{\mathit{i}}}=\mathit{E_{\mathit{i}}}(0\leq i<n)$,
which can be divided into proving $\forall
i\in\left\\{{1..n}\right\\}.\Theta\mid[]:\tau_{\mathit{i}}\operatorname{\Rightarrow}\Gamma_{\mathit{i}}\mid\vec{\ell}_{{\mathit{i}-1}}\vdash_{\mathit{ectx}}\mathit{E^{\prime}_{{\mathit{i}-1}}}:\tau_{{\mathit{i}-1}}\operatorname{\Rightarrow}\Gamma_{{\mathit{i}-1}}$
and
$\Theta\mid[]:\tau_{n+1}\operatorname{\Rightarrow}\Gamma_{n+1}\mid\vec{\ell}_{\mathit{n}}\vdash_{\mathit{ectx}}\mathit{E^{\prime}_{\mathit{n}}}:\tau_{\mathit{n}}\operatorname{\Rightarrow}\Gamma_{\mathit{n}}$.
The first follows by inversion on
$\vdash_{\mathit{conf}}^{D}\left\langle{\mathit{H},\mathit{R},\vec{F},\mathit{E}[\operatorname{\mathbf{let}}\mathit{x}=\mathit{f}^{\ell}(\mathit{y_{{\mathrm{1}}}},\ldots,\mathit{y_{\mathit{n}}})\operatorname{\mathbf{in}}\mathit{e^{\prime}}]}\right\rangle$.
To show the latter, we define $\Gamma_{n+1}=\Gamma^{\prime\prime\prime}$ and
$\tau_{n+1}=\tau^{\prime\prime\prime}$, whereby the well-typing holds from the
result of applying Lemma 28 above.

Finally, $\mathbf{Cons}\,(\mathit{H},\mathit{R},\Gamma^{\prime\prime})$
follows immediately from $\mathbf{Cons}\,(\mathit{H},\mathit{R},\Gamma)$ and
$\Gamma^{\prime\prime}=\Gamma$.

##  Appendix 0.D Proof of Progress

We first state the standard decomposition lemma.

######  Lemma 29 (Decomposition)

For any term $\mathit{e}$, either $\mathit{e}=\mathit{x}$ for some
$\mathit{x}$ or there exists some $\mathit{E}$ and $\mathit{e^{\prime}}$ where
$\mathit{E}[\mathit{e^{\prime}}]=\mathit{e}$ and one of the following cases
hold:

  1. 1.

$\mathit{e^{\prime}}=\operatorname{\mathbf{let}}\mathit{x}=\operatorname{\mathbf{mkref}}\mathit{y}\operatorname{\mathbf{in}}\mathit{e^{\prime\prime}}$

  2. 2.

$\mathit{e^{\prime}}=\operatorname{\mathbf{let}}\mathit{x}=\mathit{y}\operatorname{\mathbf{in}}\mathit{e^{\prime\prime}}$

  3. 3.

$\mathit{e^{\prime}}=\operatorname{\mathbf{let}}\mathit{x}=n\operatorname{\mathbf{in}}\mathit{e^{\prime\prime}}$

  4. 4.

$\mathit{e^{\prime}}=\operatorname{\mathbf{let}}\mathit{x}=*\mathit{y}\operatorname{\mathbf{in}}\mathit{e^{\prime\prime}}$

  5. 5.

$\mathit{e^{\prime}}=\operatorname{\mathbf{let}}\mathit{x}=\mathit{f}^{\ell}(\mathit{y_{{\mathrm{1}}}},\ldots,\mathit{y_{\mathit{n}}})\operatorname{\mathbf{in}}\mathit{e^{\prime\prime}}$

  6. 6.

$\mathit{e^{\prime}}=\mathit{x}\operatorname{;}\mathit{e^{\prime\prime}}$

  7. 7.

$\mathit{e^{\prime}}=\mathbf{alias}(\mathit{x}=\mathit{y})\operatorname{;}\mathit{e^{\prime\prime}}$

  8. 8.

$\mathit{e^{\prime}}=\mathbf{alias}(\mathit{x}=*\mathit{y})\operatorname{;}\mathit{e^{\prime\prime}}$

  9. 9.

$\mathit{e^{\prime}}=\operatorname{\mathbf{ifz}}\mathit{x}\operatorname{\mathbf{then}}\mathit{e_{{\mathrm{1}}}}\operatorname{\mathbf{else}}\mathit{e_{{\mathrm{2}}}}$

  10. 10.

$\mathit{e^{\prime}}=\mathbf{assert}(\varphi)\operatorname{;}\mathit{e^{\prime\prime}}$

  11. 11.

$\mathit{e^{\prime}}=\mathit{x}\operatorname{:=}\mathit{y}\operatorname{;}\mathit{e}$

###### Proof

Straightforward induction on $\mathit{e}$.

######  Proof (Progress; Lemma 4)

By inversion on $\vdash_{\mathit{conf}}^{D}\mathbf{C}$, either
$\mathbf{C}=\mathbf{AliasFail}$ or
$\mathbf{C}=\left\langle{\mathit{H},\mathit{R},\vec{F},\mathit{e}}\right\rangle$.
In the former case the result is immediate. In the latter case we have that
$\Theta\mid\vec{\ell}\mid\Gamma\vdash\mathit{e}:\tau\operatorname{\Rightarrow}\Gamma^{\prime}$
for some $\tau,\Gamma$ and $\Gamma^{\prime}$, and further from Lemma 29, we
have that either $\mathit{e}=\mathit{x}$ for some $\mathit{x}$ or there exists
some $\mathit{E}$ or $\mathit{e^{\prime}}$ where
$\mathit{e}=\mathit{E}[\mathit{e^{\prime}}]$ and $\mathit{e^{\prime}}$ meets
one of the cases in Lemma 29.

In the case $\mathit{e}=\mathit{x}$, we further make case analysis on the form
of $\vec{F}$. The case where $\vec{F}=\epsilon$ is immediate; In the other
case where $\vec{F}=F:\vec{F}^{\prime}$, the configuration can step to
$\left\langle{\mathit{H},\mathit{R},\vec{F},F[\mathit{x}]}\right\rangle$
according to R-Var.

For the remaining cases where $\mathit{e}=\mathit{E}[\mathit{e^{\prime}}]$, by
the well-typing of $\mathit{e}$ with respect to $\Gamma$ and Lemma 7, we have
that
$\Theta\mid\mathcal{L}\mid\Gamma\vdash\mathit{e^{\prime}}:\tau_{{\mathrm{0}}}\operatorname{\Rightarrow}\Gamma_{{\mathrm{0}}}$
some $\tau_{{\mathrm{0}}}$ and $\Gamma_{{\mathrm{0}}}$.

We now treat the remaining forms of $\mathit{e^{\prime}}$

Case:
$\quad\begin{array}[t]{@{}l@{\qquad}l@{\qquad}l@{\qquad}l}\mathit{e^{\prime}}=\operatorname{\mathbf{let}}\mathit{x}=*\mathit{y}\operatorname{\mathbf{in}}\mathit{e^{\prime\prime}}\end{array}$
By inversion (Lemma 6) and Lemma 12 we must have that for some
$\Gamma_{\mathit{p}}$ where
$\mathbf{Cons}\,(\mathit{H},\mathit{R},\Gamma_{\mathit{p}})$ that
$\mathit{y}\in\mathit{dom}(\Gamma_{\mathit{p}})$ and
$\Gamma_{\mathit{p}}(\mathit{y})=\tau^{\prime}\operatorname{\mathbf{ref}}^{r}$.
From $\mathbf{Cons}\,(\mathit{H},\mathit{R},\Gamma_{\mathit{p}})$ we must have
$\mathit{y}\in\mathit{dom}(\mathit{R})$ and further
$\mathbf{SATv}(\mathit{H},\mathit{R},\mathit{R}(\mathit{y}),\tau^{\prime}\operatorname{\mathbf{ref}}^{r^{\prime}})$,
from which we must have $\mathit{R}(\mathit{y})\,=\,\mathit{a}$ and
$\mathit{a}\in\mathit{dom}(\mathit{H})$. Then $\mathbf{C}$ can step according
to R-Deref.

Case:
$\quad\begin{array}[t]{@{}l@{\qquad}l@{\qquad}l@{\qquad}l}\mathit{e^{\prime}}=\operatorname{\mathbf{let}}\mathit{x}=\mathit{y}\operatorname{\mathbf{in}}\mathit{e^{\prime\prime}}\end{array}$
Again, by Lemmas 6 and 12 and the definition of $\mathbf{Cons}$, we must have
that $\mathit{y}\in\mathit{dom}(\mathit{R})$, and the system can step
according to R-LetVar.

Case:
$\quad\begin{array}[t]{@{}l@{\qquad}l@{\qquad}l@{\qquad}l}\mathit{e^{\prime}}=\operatorname{\mathbf{let}}\mathit{x}=\operatorname{\mathbf{mkref}}\mathit{y}\operatorname{\mathbf{in}}\mathit{e^{\prime\prime}}\end{array}$
Similar to the R-LetVar case above.

Case:
$\quad\begin{array}[t]{@{}l@{\qquad}l@{\qquad}l@{\qquad}l}\mathit{e^{\prime}}=\operatorname{\mathbf{let}}\mathit{x}=n\operatorname{\mathbf{in}}\mathit{e^{\prime\prime}}\\\
\mathit{e^{\prime}}=\mathit{x}\operatorname{;}\mathit{e^{\prime\prime}}\\\
\mathit{e^{\prime}}=\mathbf{assert}(\varphi)\operatorname{;}\mathit{e^{\prime\prime}}\end{array}$
The first two can trivially step according to R-LetInt and R-Seq respectively.
the last can step according to R-Assert or R-AssertFalse (although by Lemmas 3
and 2 the latter is impossible).

Case:
$\quad\begin{array}[t]{@{}l@{\qquad}l@{\qquad}l@{\qquad}l}\mathit{e^{\prime}}=\mathbf{alias}(\mathit{x}=\mathit{y})\operatorname{;}\mathit{e^{\prime\prime}}\end{array}$
Again by Lemmas 6 and 12 and that
$\mathbf{Cons}\,(\mathit{H},\mathit{R},\Gamma_{\mathit{p}})$ implies
$\mathit{x}$ and $\mathit{y}$ are bound to addresses in the register file, we
have that the configuration can step according to R-Alias or R-AliasFail.

Case:
$\quad\begin{array}[t]{@{}l@{\qquad}l@{\qquad}l@{\qquad}l}\mathit{e^{\prime}}=\mathbf{alias}(\mathit{x}=*\mathit{y})\operatorname{;}\mathit{e^{\prime\prime}}\end{array}$
Similar to the case above, we must have that $\mathit{x}$ is bound to an
address in the register file, and that $\mathit{y}$ is bound to an address
that is itself mapped to an address in the heap $\mathit{H}$. Then the
configuration may step according to R-AliasPtr or R-AliasPtrFail

Case:
$\quad\begin{array}[t]{@{}l@{\qquad}l@{\qquad}l@{\qquad}l}\mathit{e}=\operatorname{\mathbf{ifz}}\mathit{x}\operatorname{\mathbf{then}}\mathit{e_{{\mathrm{1}}}}\operatorname{\mathbf{else}}\mathit{e_{{\mathrm{2}}}}\end{array}$
As above, from the well-typing we must have that $\mathit{x}$ is bound in
$\mathit{R}$ to some integer $n$. Then the configuration may step according to
R-IfTrue or R-IfFalse depending on whether $n\,=\,0$ or $n\,\neq\,0$.

Case:
$\quad\begin{array}[t]{@{}l@{\qquad}l@{\qquad}l@{\qquad}l}\mathit{e^{\prime}}=\mathit{x}\operatorname{:=}\mathit{y}\operatorname{;}\mathit{e^{\prime\prime}}\end{array}$
From the well-typing of $\mathit{e^{\prime}}$, Lemmas 6 and 12 and the
definition of $\mathbf{Cons}$, we must have that
$\mathit{y}\in\mathit{dom}(\mathit{R})$,
$\mathit{x}\in\mathit{dom}(\mathit{R})$,
$\mathit{R}(\mathit{x})\,=\,\mathit{a}$, and
$\mathit{a}\in\mathit{dom}(\mathit{H})$. Then the configuration can step
according R-Assign.

Case:
$\quad\begin{array}[t]{@{}l@{\qquad}l@{\qquad}l@{\qquad}l}\mathit{e^{\prime}}=\operatorname{\mathbf{let}}\mathit{x}=\mathit{f}^{\ell}(\mathit{y_{{\mathrm{1}}}},\ldots,\mathit{y_{\mathit{n}}})\operatorname{\mathbf{in}}\mathit{e^{\prime\prime}}\end{array}$
From the well-typing of the function call we must have that
$\mathit{f}\in\mathit{dom}(\Theta)$. From $\Theta\vdash\mathit{D}$ in the
precondition of $\vdash_{\mathit{conf}}^{D}\mathbf{C}$, we must have that
$\mathit{f}\mapsto(\mathit{x_{{\mathrm{1}}}},\,..\,,\mathit{x_{\mathit{j}}})\mathit{e^{\prime\prime\prime}}\in\mathit{D}$.
Then from T-FunDef we must have that $j=n$ whereby the configuration can step
according to R-Call.