††thanks: email<EMAIL_ADDRESS>
# Cycle equivalence classes, orthogonal Weingarten calculus,
and the mean field theory of memristive systems

F. Caravelli  Theoretical Division (T4), Los Alamos National Laboratory, Los
Alamos, New Mexico 87545, USA

###### Abstract

It has been recently noted that for a class of dynamical systems with explicit
conservation laws represented via projector operators, the dynamics can be
understood in terms of lower dimensional equations. This is the case, for
instance, of memristive circuits. Memristive systems are important classes of
devices with wide-ranging applications in electronic circuits, artificial
neural networks, and memory storage. We show that such mean-field theories can
emerge from averages over the group of orthogonal matrices, interpreted as
cycle-preserving transformations applied to the projector operator describing
Kirchhoff’s laws. Our results provide insights into the fundamental principles
underlying the behavior of resistive and memristive circuits and highlight the
importance of conservation laws for their mean-field theories. In addition, we
argue that our results shed light on the nature of the critical avalanches
observed in quasi-two dimensional nanowires as boundary phenomena.

††preprint: LA-URed

##  I Introduction

With the increasing interest in developing neuromorphic computers, it is
crucial to understand physical devices that exhibit similar structures and
functionalities to biological neural networks. The nonlinear interactions and
complex connectivity of biological neuronal networks are well-known
characteristics [1]. Yet, it is still a mystery why low-dimensional
representations of certain rather complex dynamical systems exist [2, 3] even
if in certain regimes. This study sheds light on the existence of such
representations in the context of memristive circuits, which we use as a toy
model for more complex neuromorphic systems [4].

Circuits composed of nanodevices with memory are at the forefront of
neuromorphic computing research, as their behavior often mimics synaptic
plasticity observed in biological neuronal circuits. However, a comprehensive
theory that effectively describes the behavior of these circuits is currently
lacking. Memristive devices are a promising area of research for the
development of next-generation computing systems [5, 6]. These devices are
passive 1-port resistive components that have the ability to remember past
voltages and currents and can change their resistance based on the history of
the input signals [7, 8]. The experimental existence of switching in physical
systems dates back to the late 60’s [9], but the connection to memristive
behavior has been made just over a decade ago [10, 11, 12]. One peculiar
feature of the memristive devices is that their resistance changes between two
limiting values $R_{on}\leqslant R_{off}$ (or analogously conductance value).
The development of circuits of memristive devices has become an important area
of research, as it enables the creation of neuromorphic devices that can
support the existing von Neumann architecture [13, 6, 5] in a variety of tasks
more prone to an analog computing approach. From an experimental perspective,
nanowire networks have emerged as a promising material for the fabrication of
disordered memristive networks. They exhibit reversible resistance switching
behavior when subjected to an external electric field, making them ideal for
use in memristive devices. Additionally, silver nanowires are low-cost, have
high aspect ratios, and can be synthesized using a variety of techniques,
making them highly versatile. The avoided crossings between the nanowires act
as tunneling junctions [14, 15], and for coated silver nanowires the
phenomenon of filament formation is the main driver between the memristive
effects that emerge [16, 17, 14, 18]. Their behavior is particularly similar
to the behavior of neuronal circuits, first and foremost, their connectivity
strongly resembles the one of a neuronal connectome. In addition, the
formation of filaments and their effective memristive behavior strongly
parallels the plasticity of neuronal circuits. In this respect, then,
memristive networks are an area of research in Physics that parallels the
study of neuronal networks in biology. For instance, memristive networks have
Lyapunov functions [19] similar to recurrent neural networks [20]. Studying
memristive circuits can provide important insights into the non-trivial
dynamics of large assemblies of neuronal networks.

As an example of this, it became clear only recently that certain disordered
circuits with memory exhibit a certain regularity [7, 10, 21, 22, 23] in their
dynamical behavior as the applied voltage is constant in time. This should be
somewhat surprising, given that these are rather nonlinearly interacting
systems, characterized by conservation laws that lead to effective nonlocal
interactions, although exponentially bounded [24]. Such mean field techniques
have also been tested in experimentally viable systems such as nanowire
connectomes [25]. Although it is very well-known that circuits composed of
purely memristive devices must also be memristive in 2-probe experiments [8],
it is less obvious why the collective behavior should be so much similar to a
single device, in particular when the network of devices is disordered. This
observation is not only numerical or model-based, but it has been recently
shown to be true also in the experimental setup of silver nanowires [25],
where a mean field theoretical description of a disordered network of
nanowires could fit the potentiation and depression of the conductance.
Following these results, the author proposed that there is a class of
dynamical systems that, thanks to the presence of projector operators in their
dynamics, have a well-defined mean field theory [26]. However, rigorous
results could be obtained only for simple projector operators.

On a side note, understanding the phenomenology of memristive networks has
direct applications for the employment of these physical systems as reservoir
computing devices at the edge of chaos, an idea proposed in the early 90’s
[27, 28], which has seen a direct application in memristive circuits [29, 30,
31, 32]. In fact, memristive systems seem to perform better near critical
transitions [33, 34]. In the case of the brain, the idea that the brain is in
a constant critical state has been proposed in the mid 90’s [35, 36, 37, 38,
39, 40].

However, in the case of memristive circuits, one can have both first and
second-order dynamical transitions in the conductance state. Avalanches of
memristive switching are a phenomenology also discussed in the literature [41,
42, 43] both theoretically and experimentally. However, it is worth pointing
out that these switching phenomena are strongly due to what we wish to call
here boundary phenomena: they depend on how the memristive device approaches
dynamically the boundary of their resistive or conductance values. We call
instead a bulk phenomenon the switching is purely due to the presence of
hysteresis, e.g. a pinched hysteresis loop. Boundary phenomena are based on a
rapid switching of memristive devices between two resistive states. Bulk
phenomena are instead due only to the presence of memory. This is for instance
the case of the recent symmetry-breaking transition of [44]. There, a first
order transition in the conductance state is purely due to the presence of
hysteresis, rather than the phenomenology of the boundary of the device. In
fact, the presence or not of this transition depends only on the ratio
$r=R_{off}/R_{on}$, which at present is the best indicator of a bulk-induced
transition. Previously, a mean-field theory for this type of transition was
developed in a series of works [44, 25, 26]. Other mean field theories have
appeared in the literature [45, 19, 46], describing the Ising-like behavior.
The present manuscript is an attempt of overcoming some of the limitations of
previous works, incorporating the circuit properties into mean field theory.
Bulk in deriving a mean-field theory for memristive network bulk transitions,
by averaging only in a subclass with the same number of cycles. This is
important for experimental reasons. We know for instance that in quasi-two
dimensional materials that these second order transitions are only a cross
over [25]. One reason is that the mean field theoretical results imply
implicitly the homogeneity of the network. We thus wish to go beyond this
approximation, applying techniques that keep track of the density of cycles in
the network. This result is important, even if for a toy model, to develop
techniques to explain the phenomenology of more realistic systems.

The paper is organized as follows. First, we establish a connection between
cycle-preserving transformations and the group of orthogonal transformations,
introducing an equivalence relation on the cycle space. We then utilize the
Weingarten calculus to average over orthogonal transformations, leading to a
mean-field theory that extends previous results and renormalizes the ratio
$r$. We show that the cycle density directly affects the transition from low
to high conductance. Conclusions follow.

##  II Cycle equivalence classes and memristive circuits

###  II.1 Circuits, Kirchhoff’s laws, and projector operators

In circuits, Kirchhoff’s laws are manifestations of the conservation of
physical quantities such as charge and energy. Mathematically, these can be
expressed via the introduction of projection operators [47, 21], i.e. matrices
$\bm{\mathbf{\Omega}}$ satisfying the constraint
$\bm{\mathbf{\Omega}}^{2}=\bm{\mathbf{\Omega}}$, and directly connected to
circuit topology. For instance, for a resistive circuit made of identical
resistances of value $r$ in series with voltage generators in series, Ohm’s
law for the network can be expressed as

$\displaystyle\vec{i}=\frac{1}{r}\bm{\mathbf{\Omega}}\vec{s},$ (1)

where $\vec{s}$ is the collection of voltage generators connected in series to
each resistance, while $\vec{i}$ contains the branch currents. The value of
the current is connected to the effective graph resistance [48]. For a circuit
of identical resistors, if one picks two nodes $a,b$ and applies a unital
voltage, the effective resistance $R(a,b)$ is equivalent to the inverse of the
effective current that flows through the generator. The effective graph
resistance is $R_{eff}=\sum_{(a,b)}R_{ab}.$ The graph resistance can be
evaluated via the graph Laplacian, defined as the matrix $Q_{ab}=1$ if $a=b$,
$-1$ if $(a,b)$ is an edge of the graph and 0 otherwise. Then, if
$\vec{e}_{a}$ is the vector with all zeros and $1$ in position $a$, we have
that

$\displaystyle
R(a,b)=(\vec{e}_{a}-\vec{e}_{b})^{t}Q^{-1}(\vec{e}_{a}-\vec{e}_{b})$ (2)

where $Q^{-1}$ is intended as the pseudo-inverse of $Q$, e.g. the matrix
obtained by diagonalizing the matrix, and inverting the diagonal matrix of
eigenvalues for only the non-zero eigenvalues. It follows that
$R_{eff}=V\sum_{i=2}^{V}\mu_{i}^{-1}$, $\mu_{i}$ are the eigenvalues of $Q$.
It is known using the spectrum of the Laplacian that for a complete graph
$R_{eff}=V-1$.

Figure 1: Cycles in graphs. Given an oriented graph, we assign an orientation
to the edges and each cycle. Given a spanning tree $\mathcal{T}$, the number
of fundamental cycles are associated to the element not in $\mathcal{T}$.

.

###  II.2 Graph theoretic formalism

The underlying assumption of (1) is that the voltage generators $s_{i}$’s are
in series to the resistances $i$’s, while the circuit can be represented as a
graph with $E$ edges. Given the branch currents and a certain orientation of
the (fundamental) graph cycles $1,\dots,L$, we can obtain the so-called cycle
matrix of the circuit $A$, of size $L\times E$, such that
$\bm{\mathbf{\Omega}}=\bm{A}^{t}(\bm{A}\bm{A}^{t})^{-1}\bm{A}$, where t
denotes the transpose. The fundamental cycles $L$ can be obtained by picking
any spanning tree $\mathcal{T}$, and associating the fundamental cycles in the
complement of the tree $\bar{\mathcal{T}}$.111Details about this construction
can be found in [21].

Although the details of the derivation of $\bm{\mathbf{\Omega}}$ from the
circuit topology are beyond the scope of this paper, where
$\bm{\mathbf{\Omega}}$ will be kept generic and unrelated to any underlying
graph, it is worth stressing its relation to the graph structure of the
circuit. We begin by introducing some concepts from graph theory before
defining the projector operator $\mathbf{\Omega}$. A _walk_ on a directed
graph $G$ is a sequence of vertices, denoted as $v_{1},v_{2},\dots,v_{n}$,
such that for every pair of consecutive vertices in the sequence, there exists
an edge that connects them. Formally, for every $i$ such that $1<i\leqslant
n$, either $(v_{i-1},v_{i})\in E(G)$ or $(v_{i},v_{i-1})\in E(G)$, where
$E(G)$ is the edge set of the graph $G$.

A _cycle_ is a walk denoted as $W=v_{1}v_{2}\dots v_{n}$, such that
$n\geqslant 3$, $v_{0}=v_{n}$, and the vertices $v_{i}$, $0<i<n$, are distinct
from each other and from $v_{0}$. An example of a cycle is shown in Figure 1.
The space spanned by the edges of the graph has a vector space structure. The
_cycle space_ is the subset of the edge space that is spanned by all the
cycles of the graph.

A _cycle matrix_ $A$ is a matrix whose columns form the basis of the cycle
space. For example, if ${\vec{c}_{1},\dots,\vec{c}_{n}}$ is a set of column
vectors that form a basis of the cycle space, then the cycle matrix is
$A=(\vec{c}_{1},\dots,\vec{c}_{n})$. Finally, the _projector operator
$\mathbf{\Omega}$ on the cycle space of the graph_ is defined as
$\mathbf{\Omega}=A(A^{T}A)^{-1}A^{T}$. The issue with this formalization is
that calculating the cycle matrix is cumbersome, as this is a nonlocal
quantity. Evaluating $\mathbf{\Omega}$ is easier instead if one uses instead a
local quantity, such as the local node connectivity. Let $G$ be a directed
graph. One way of representing $G$ is by specifying where each of its edges
starts and where it ends. It is convenient to do this using a matrix. We call
such a matrix an _incidence matrix_. Each column of its incidence matrix
represents an edge: the first edge starts at vertex $1$ and ends at vertex
$2$, so the first column of the matrix has entry $1$ in the first row and
entry $-1$ in the second row. All the other entries in the first column are
$0$ because none of the other vertices are a part of that edge. By continuing
this process for every edge, we get the incidence matrix $B$ of the given
graph:

$B=\begin{pmatrix}1&-1&-1&0&0&0&0&0\\\ -1&0&0&1&1&0&0&0\\\
0&1&0&-1&0&-1&-1&0\\\ 0&0&0&0&0&1&0&-1\\\ 0&0&1&0&-1&0&1&1\\\ \end{pmatrix}$
(3)

More formally, if a graph $G$ has $v$ vertices and $e$ edges, then the
incidence matrix $B$ of $G$ is a $v\times e$ matrix (i.e. a matrix with $v$
rows and $e$ columns), whose entry $(i,j)$ is defined as

$B_{ij}\coloneqq\begin{cases}1&\text{if $v_{i}$ is the initial vertex of the
edge $e_{j}$,}\\\ -1&\text{if $v_{i}$ is the terminal vertex of the edge
$e_{j}$,}\\\ 0&\text{otherwise.}\end{cases}$ (4)

If $B$ is an incidence matrix, we can define the projector operator
$\mathbf{\Omega}_{B^{T}}$:

$\mathbf{\Omega}_{B^{T}}=B^{T}\left(BB^{T}\right)^{-1}B,\ \ \
\mathbf{\Omega}_{B^{T}}^{2}=\mathbf{\Omega}_{B^{T}}.$ (5)

If we try to compute $\mathbf{\Omega}_{B^{T}}$ from the definition, we will
find that the inverse $\left(BB^{T}\right)^{-1}$ does not exist in general.
This can be solved by either considering the reduced incidence matrix
(obtained by removing a row from the original incidence matrix) or by taking
the pseudoinverse of the expression instead of the “regular” inverse.

The following useful identity connects the projector operator
$\mathbf{\Omega}$ (based on the matrix $A$) to the projector operator
$\mathbf{\Omega}_{B^{T}}$ (based on the incidence matrix $B$):
$\mathbf{\Omega}=\mathbb{I}-\mathbf{\Omega}_{B^{T}}$, where $\mathbb{I}$ is
the identity matrix. Such duality between the incidence and cycle adjacency
matrix is instead a manifestation of the cycle and nodal analysis in circuit
theory. Such a relationship is based on the fact that $\mathbf{\Omega}$ and
$\mathbf{\Omega}_{B^{T}}$ are orthogonal matrices and the fact that
$AB^{t}=0$, an identity related to the conservation of energy called
Tellegen’s theorem [24].

Figure 2: Example of an isospectral transformation for the projector operator
$\mathbf{\mathbf{\Omega}}$. The number of fundamental cycles is preserved, and
thus $|Span(A)|$ and the number of edges is preserved.

###  II.3 Orthogonal group and cycle preserving transformations

In essence, the projector $\mathbf{\Omega}$ captures the part of a vector that
can be represented as a combination of cycles in the graph, and discards any
components that do not lie within the cycle space. Let us then introduce a set
of transformations that preserve its spectrum.

######  Definition 1 (Cycle preserving transformation).

. Let $G=(\mathcal{E},\mathcal{V})$ be a graph. Let $\theta:G\rightarrow G$ be
a map, such that for $\mathcal{G}\in G$, we have
$|\mathcal{L}(\mathcal{G})|=|\mathcal{L}\big{(}\theta(\mathcal{G})\big{)}|$
and
$|\mathcal{E}(\mathcal{G})|=|\mathcal{E}\big{(}\theta(\mathcal{G})\big{)}|.$
Then, $\theta$ is cycle preserving.

Cycle preserving transformations can be thought as transformations in which an
edge is disconnected from the graph, and connected to two other nodes where
such edge was not present, as we consider for simplicity simple graphs. An
example of such transformation is shown in Fig. 2. As we can see, the
transformation preserves the total number of edges.

Let us now introduce a very simple result, but key to the rest of this
manuscript. The significance of the following lemma lies in its ability to
establish a crucial connection between edge-preserving graph transformations
and the orthogonal group. Despite its simplicity, this lemma plays a pivotal
role in the following, by providing a fundamental insight that forms the basis
for our further analyses and conclusions.

######  Proposition 1.

Let $\mathcal{G}$ be a graph $\mathcal{G}=(\mathcal{E},\mathcal{V})$. Let
$\theta$ be cycle preserving. If $\mathbf{\mathbf{\Omega}}_{A(\mathcal{G})}$
is the projector operator on $\mathcal{L}=\text{Span}(A)$, we have

$\mathbf{\mathbf{\Omega}}^{\prime}\equiv\mathbf{\mathbf{\Omega}}_{A\big{(}\theta(\mathcal{G})\big{)}}=\mathscr{O}^{t}\mathbf{\mathbf{\Omega}}_{A(\mathcal{G})}\mathscr{O}$
(6)

for some orthogonal matrix $\mathscr{O}$.

###### Proof.

. Since $\theta$ preserves the cycle space, we have
$L=|\mathcal{L}|=|\mathcal{L}^{\prime}|$ and
$|\mathcal{E}|=|\mathcal{E}^{\prime}|$. Then, we know that
$\text{Span}(A)=\text{Span}(A^{\prime})$, since
$\mathbf{\mathbf{\Omega}}^{2}=\mathbf{\mathbf{\Omega}}$, we have that the
spectrum is given by $\Lambda(\mathbf{\mathbf{\Omega}})=\\{0^{E-L},1^{L}\\}$
is preserved. As such, $\mathbf{\mathbf{\Omega}}$ and
$\mathbf{\mathbf{\Omega}}^{\prime}$ are isospectral, and thus since both
$\mathbf{\mathbf{\Omega}}$ and $\mathbf{\mathbf{\Omega}}^{\prime}$ are real,
$\mathbf{\mathbf{\Omega}}^{\prime}=\mathscr{O}^{t}\mathbf{\mathbf{\Omega}}\mathscr{O}$
for some orthogonal matrix $\mathscr{O}$. ∎

We say that two graphs $\mathcal{G}$ and $\mathcal{G}^{\prime}$ are conjugate
if $\mathbf{\Omega}$ and $\mathbf{\Omega}^{\prime}$ are related by an
orthogonal transformation. It is easy to see that, in fact, such
transformations form an equivalence class. Such a simple result will go a long
way in this paper, as it allows us to perform averages over the orthogonal
transformations and derive a mean field result in each conjugacy class, e.g.
averaging on graphs with the same number of cycles.

To see a direct application of this result, consider again the solution for a
resistive circuit, eqn. (1).

Figure 3: Examples of partitioning in equivalence classes for $V=5$ for graphs
representing circuits. As we can see, the elements in each class are related
by an edge-preserving transformation. We characterize each equivalence class
as $\mathrm{C}(E,L)$, or alternatively $\mathrm{C}(E,V)$.

If two graphs representing the resistive circuits are conjugate, then because
of Proposition 1 we know that

$\displaystyle\vec{i}$ $\displaystyle=$ $\displaystyle\mathbf{\Omega}\vec{s}$
(7) $\displaystyle\vec{i}^{\prime}$ $\displaystyle=$
$\displaystyle\mathscr{O}^{t}\mathbf{\Omega}\mathscr{O}\vec{s}$ (8)

Which implies immediately that

$\displaystyle\mathscr{O}\vec{i}^{\prime}=\mathbf{\Omega}\mathscr{O}\vec{s},$
(9)

from which we get that up to a rotation of currents and voltages, the two are
identical. This implies that there is a class of circuits whose currents and
applied voltages can be mapped to each other via an orthogonal transformation.
This is why we can introduce an equivalence class induced by such
transformations. We can then classify the number of such equivalence classes.

######  Definition 2.

We say that $\mathcal{G}\sim\mathcal{G}^{\prime}$ if
$\mathbf{\Omega}(\mathcal{G})=\mathscr{O}^{t}\mathbf{\Omega}(\mathcal{G}^{\prime})\mathscr{O}$.

Using the equivalence relation, we have the following lemma for circuits

######  Proposition 2.

Let $\mathcal{G}$ be a graph on $V$ nodes representing a circuit, e.g. there
are no nodes of degree 1 and $\mathcal{G}$ is connected. Then there are at
most $\frac{1}{2}(V^{2}-3V+2)$ and minimum $2V-1-\big{(}V\ \text{mod}\ 2)/2$
cycles respectively. This classifies the equivalence classes with fixed nodes.

###### Proof.

The number of equivalence classes is the number of possible spectral
configurations of $\mathbf{\Omega}$. This is given in the most general case,
$\Lambda(\mathbf{\Omega})=\\{1^{L},0^{E-L}\\}$. Circuits must satisfy certain
relationships between the number of cycles, the number of fundamental faces
and the number of vertices. First, note that the ordering of spectrum is not
important, as we can permute the eigenvalues via a permutation matrix $P$, and
then introduce the orthogonal matrix $\mathscr{O}\rightarrow
P\mathscr{O}P^{t}$ to reorder them. Since $P\mathscr{O}P^{t}$ is still an
orthogonal matrix, permutations of spectra belog to the same class. Thus, the
eigenvalue ordering is not important, but only the number of 1’s and 0’s, the
spectral fingerprint. The number of fundamental cycles is equal to the number
$L=E-|\mathcal{T}|$ where $\mathcal{T}$ is the cardinality of a spanning tree.
Because the graph of a circuit is connected, a spanning tree always exists
with the number of edges $|\mathcal{T}|=V-1$ (e.g. we do not consider
forests). It is easy to see that the number of leaves is even if $V\
\text{mod}\ 2=0$ and the number of leaves is odd if $V\ \text{mod}\ 2=1$. Let
us assume $V=2k$ even or $V=2k+1$ odd with $k\geqslant 1$. Then, we need to
add at least $k$ edges between the leaves to form $k$ cycles. Thus, the
minimum number of fundamental cycles is $k$. We can then add one edge at a
time until we obtain a complete graph. Then, for a graph with $V=2k(+1)$
nodes, we have at most $L=V(V-1)/2-V+1=\frac{1}{2}(V^{2}-3V+2)$ (the graph
$\mathcal{K}_{V}$), and minimum $L=V-1+V-\big{(}V\ \text{mod}\
2)/2\big{)}=2V-1-\big{(}V\ \text{mod}\ 2)/2\big{)}$ cycles (the graph in which
we added the minimum number of edges to remove the leaves from the spanning
tree). ∎

######  Corollary 1.

All graphs with fixed number of edges and vertices
$\mathcal{G}=(\mathcal{E},\mathcal{V})$ representing a circuit belong to the
same equivalence class under $\sim$. We call this equivalence class
$\mathrm{C(E,V)}$.

###### Proof.

This is a consequence of Proposition 1 and Proposition 2. If $V$ is fixed, and
$E$ is fixed, then $L=E-V+1$. Then, if $E$ and $V$ are fixed, then
$\forall\mathcal{G}$ such that $|\mathcal{E}|=E$ and $|\mathcal{V}|=V$ the
number of fundamental cycles is fixed. Then all graphs in the same equivalence
class are related to each other via local edge moves as in Fig. 2 which
preserve the number of edges and nodes. ∎

The previous results establish that in order to remain into the same
equivalence class it is sufficient to perform moves that preserve the number
of edges, given the number of vertices. Then, it follows that
$\mathbf{\Omega}\sim\mathbf{\Omega}^{\prime}$, then we know the number of
edges and the number of vertices must be identical. An example of such
equivalence class partitioning is shown in Fig. 3.

Before we continue, let us make some comments on the relation between number
of cycles, edges and nodes. In particular, it will be important to calculate
the ration $L/E$. For a connected graph, the cardinality of a spanning tree is
equal to $V-1$, and thus

$\displaystyle\frac{L}{E}=\frac{E-V+1}{E}=1-\frac{V}{E}+\frac{1}{E}.$ (10)

Because of the handshaking lemma, we know that $V\bar{d}=2E$, from which we
obtain that $\frac{V}{E}=\frac{2}{\bar{d}}$, where $\bar{d}$ is the average
degree. We then have

$\displaystyle\frac{L}{E}=1-\frac{2}{\bar{d}}+\frac{1}{E}.$ (11)

We see from the equation above that for complete graphs this ratio goes to
one.

###  II.4 Averaging over the equivalence class

One question we might ask is then: what is the average current vector on the
equivalence class determined by cycle isomorphism? What we can do is average
over all orthogonal matrices which span an equivalence class
$\mathrm{C}(E,L)$, and write

$\displaystyle\langle\vec{i}\rangle_{\mathscr{O}}=\frac{1}{r}\langle\mathscr{O}^{t}\mathbf{\Omega}\mathscr{O}\rangle_{\mathscr{O}}\vec{s}.$
(12)

We can identify one element in the class $\mathrm{C}(E,V)$, e.g. the graph
$\mathcal{G}_{0}$ and associate it to the identity element
$\mathscr{O}_{0}=\mathbb{I}\in\mathbb{O}(E)$, the group of orthogonal matrices
of size $E$. Then, given $\mathbf{\Omega}_{0}$ associated to $G_{0}$, $\forall
G_{i}\in\mathrm{C}(E,V)$, there exist an orthogonal matrix $\mathscr{O}_{i}$
s.t.
$\mathbf{\Omega}_{i}=\mathscr{O}_{i}^{t}\mathbf{\Omega}_{0}\mathscr{O}_{i}$.
For generic graphs, such average is over a finite number of elements. As an
example, consider Fig. 3. For $V=5,E=5$, there is a unique representative for
this equivalence class. For $V=5,E=6$, there are instead 5 representatives.
This means that there is, up to permutations of nodes and edges, only a finite
number of orthogonal matrices to average over, and thus the average above is
over a finite sum. Let $S(E,L)$ be the number of elements in
$\mathrm{C}(E,V)$. Then, our average of eqn. (12) can be written as

$\displaystyle\langle\cdot\rangle_{\mathscr{O}}=\frac{1}{S(E,L)}\sum_{i=0}^{S(E,L)-1}\mathscr{O}_{i}^{t}\cdot\mathscr{O}_{i}$
(13)

and thus we obtain that

$\displaystyle\langle\vec{i}\rangle_{\mathscr{O}}=\frac{1}{r}\frac{1}{S(E,L)}\sum_{i=0}^{S(E,L)-1}\mathscr{O}_{i}^{t}\mathbf{\Omega}\mathscr{O}_{i}\vec{s}$
(14)

For finite $V,E$, we are not able to perform such an average in this paper.
However, a possibility might be to perform the average when the number of
edges and loops goes to infinity. In this case, the number of orthogonal
matrices also goes to infinity. A mathematical question might then be what is
the measure in the limit of $E,L$ going to infinity? Our working assumption
going forward will be that such average converges to the continuous group of
orthogonal matrices, but this is a rather non-trivial statement. In the
following, this will be our working assumption and in the worst-case scenario
our approximation. The average over continuous orthogonal matrices can be
performed using the Haar measure, e.g. we will perform the replacement

$\displaystyle\lim_{E,L\rightarrow\infty}\frac{1}{S(E,L)}\sum_{i=0}^{S(E,L)-1}\mathscr{O}_{i}^{t}\cdot\mathscr{O}_{i}\rightarrow\int_{\mathbb{O}(E)}\mathscr{O}_{i}^{t}\cdot\mathscr{O}_{i}d\mathscr{O},$

where $d\mathscr{O}$ is the continuous Haar measure over the orthogonal group
[50]. Thus the group $\mathbb{O}(E)$ is equipped with a Haar probability
measure, we can define the average

$\displaystyle\langle\cdot\rangle_{\mathscr{O}}\equiv\int_{\mathbb{O}(E)}\cdot
dO.$ (15)

Then, performing this average is equivalent to averaging over random
orthogonal matrices with a constant measure over $\mathbb{O}(E)$.

To motivate this approach, let us consider a very simple example, based on the
exact current solution of resistive network, eqn. (8). We have an expression
of the form

$\displaystyle\langle\mathscr{O}^{t}\mathbf{\Omega}\mathscr{O}\rangle_{\mathscr{O}}$
(16)

where the average is performed via the Haar measure over the orthogonal group.
Let us anticipate one result here: it is known that the operation above is the
first order isospectral twirling of the matrix $\mathbf{\Omega}$, and if one
is acquainted with these methods, it is immediate to see that
$\langle\mathscr{O}^{t}\mathbf{\Omega}\mathscr{O}\rangle_{\mathscr{O}}\propto\mathbb{I}$.
The reason why this is the case will be clear in the next section, where we
describe in detail the techniques used to perform these averages in detail. We
will just state for the sake of clarity here that such average gives

$\displaystyle\langle\vec{i}\rangle_{\mathscr{O}}=\frac{{\rm
Tr}(\mathbf{\Omega})}{rE}\vec{s}.$ (17)

Above, $E$ is the size of the matrix $\mathscr{O}$, which in our case is the
number of edges of the graph, and thus resistors. The quantity ${\rm
Tr}(\mathbf{\Omega})$ represents the matrix trace, and is simply the sum over
the eigenvalues of $\mathbf{\Omega}$. We know from the introduction that,
because the $\mathbf{\Omega}$ is a projector over the number of fundamental
cycles, then the number sum of 1’s and 0’s in the spectrum is simply $L$. It
follows that

$\displaystyle\langle\vec{i}\rangle_{\mathscr{O}}=\frac{L}{rE}\vec{s}.$ (18)

We see that the average over rotations of the matrix $\mathbf{\Omega}$ leads
to a trace, which contains the information over the network topological
properties, such as the number of cycles. The amount above is an average over
all possible real projectors with the same spectrum, and it provides
information on the expected size of the current as a function of graph
topological quantities, such as the number of fundamental cycles and the
number of edges/resistors. Thus, this result suggests that our system of
currents is equivalent to a set of uncoupled resistances of effective value
$r^{\prime}=\frac{rE}{L}$. To see whether we can make sense of this, consider
$E$ parallel resistance. We know from earlier discussions that for complete
graphs, such a ratio goes to one. This is the same result obtained from the
effective resistance analysis obtained in [48] using eqn. (2). We thus expect
our Haar average to be informative for very dense circuits.

The second question one might ask is whether this result is typical. This
means asking what is the probability that, given two isospectral graphs, a
certain measure between the two vectors is far from a certain value. For
instance, let us pick

$\displaystyle\text{Var}(\vec{i})=\|\vec{i}-\langle\vec{i}\rangle_{\mathscr{O}}\|^{2}.$
(19)

We can expand the expression above, and write it as

$\displaystyle\text{Var}(\vec{i})$ $\displaystyle=$
$\displaystyle\vec{i}\cdot\vec{i}+\langle\vec{i}\rangle_{\mathscr{O}}\cdot\langle\vec{i}\rangle_{\mathscr{O}}-2\langle\vec{i}\rangle_{\mathscr{O}}\cdot\vec{i}$
(20) $\displaystyle=$ $\displaystyle\vec{s}\
^{t}\cdot\mathbf{\Omega}\vec{s}+\frac{L^{2}}{E^{2}}\vec{s}\
^{t}\cdot\vec{s}-2\frac{L}{E}\vec{s}\ ^{t}\cdot\mathbf{\Omega}\vec{s}$ (21)

where we used $\mathbf{\Omega}^{2}=\mathbf{\Omega}$. We can average this
expression again, and obtain

$\displaystyle\langle\text{Var}(\vec{i})\rangle_{\mathscr{O}}=\frac{L}{E}\|\vec{s}\|^{2}\big{(}1-\frac{L}{E}\big{)}.$
(22)

There are now two possible cases. First, consider $\|\vec{s}\|^{2}=O(1)$ e.g.
not scaling with $E$. In the limit in which $E\gg L$, or $E=L$ such a result
is typical because of concentration inequalities such as Chebyshev
inequalities. This is the case when for instance we have a single cycle and a
large number of resistors. The second is the case when for instance the graph
is extremely dense, for instance in the case of a complete graph, in which
$L\propto E$. Using this technique, we can thus make estimates about the
typical behavior of a certain system for very dense graphs.

While resistive circuits might not be considered as interesting as other
physical systems, the point of our paper is that similar techniques can also
be used in the case of dynamical resistive circuits, such as memristive
circuits. Unfortunately, these are typically more complicated, and as we show
in this paper, even for the case of the simplest model of memristive dynamics,
exact results can be obtained only in the case of dense graphs and in the
asymptotic limit. The question is then what happens then in the case of
memristive circuits?

###  II.5 Toy model of memristive circuit dynamics

Let us now discuss the explicit equation of memristive dynamics for linear and
current controlled devices, as a toy model for the conductance transitions.
The equation of motion for a circuit of memristors has been derived in a
previous study [24, 44]. It models a flow network that adheres to current and
energy conservation laws, and the dynamics of its edges are bounded within an
interval. A memristor with memory can be described by an effective resistance
that depends on an internal parameter $x$. For example, $TiO_{2}$ memristors
can be approximated by the functional form $R(x)=R_{on}(1-x)+xR_{off}$, where
$R_{on}<R_{off}$ are the limiting resistances, and $x\in[0,1]$ physically
represents the size of the oxygen-deficient conducting layer [10]. The
internal memory parameter $x$ evolves, to the lowest order of description,
according to a simple equation:

$\frac{d}{dt}x=\frac{R_{off}}{\beta}I-\alpha
x=\frac{R_{off}}{\beta}\frac{V}{R(x)}-\alpha x,$ (23)

with hard boundaries. Here, $\alpha$ and $\beta$ are the decay constant and
the effective activation voltage per unit of time, respectively, and they
determine the timescales of the dynamical system. Many extensions of this
basic model have been considered in the literature. For example, diffusive
effects near the boundaries can be approximated by removing the hard
boundaries and multiplying by a window function [51, 52, 53]. Nonlinear
conductive effects can also be included by replacing $I$ with a function
$f(x,I)$ or introducing new parameter dependencies, such as temperature in the
case of thermistors [54]. Comparisons between these models [55, 56, 57, 58]
show that many of them are more faithful to the precise IV curves of physical
devices, but most share the basic pinched hysteresis phenomenology of the
linear model. In analytical work, it is often assumed that the dynamics are
linear in the currents in order to study the behavior in a wide context.

For a single memristor under an applied voltage $S$, Ohm’s law $S=RI$ can be
used to obtain an equation for $x(t)$ in adimensional units , given by:

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

where $\chi=\frac{R_{off}-R_{on}}{R_{off}}$ and $s=\frac{S}{\alpha\beta}$,
with $0\leqslant\chi\leqslant 1$ in physically relevant cases, and $V(x,s)$
represents an effective potential.

The dynamics of the one-dimensional dynamical system is described by an ODE,
and is fully characterized by gradient descent in the potential:

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

The equation above, in the case of a circuit of purely memristive devices with
a voltage generator in series, is generalized to a set of coupled ODE written
compactly as [44]

$\displaystyle\frac{d\vec{x}}{dt}=\frac{1}{\beta}(\mathbb{I}-\chi\mathbf{\Omega}\mathbb{X})^{-1}\mathbf{\Omega}\vec{s}-\alpha\vec{x}.$
(26)

We can see that the (26) is written for arbitrary network topologies, as these
enter only via $\mathbf{\Omega}$.

To see why $\mathbf{\Omega}$ is connected to Kirchhoff’s laws for memristive
devices, consider again eqn. (26). It is easy to see that the transformation
$\vec{s}\rightarrow\vec{s}+(I-\mathbf{\Omega})\vec{k}$ for arbitrary vectors
$\vec{k}$ leaves the dynamics invariant, since
$\mathbf{\Omega}(\mathbb{I}-\mathbf{\Omega})=0$. This is a manifestation of
Kirchhoff’s laws [24], which are however also present in eqn. (1). However, in
memristive devices described by eqn. (26) this is a subset of a larger
invariance set, which are more generally of the form

$\displaystyle(\vec{x},\vec{s})\rightarrow(\vec{x}\ ^{\prime},\vec{s}\
^{\prime})=(\vec{x}+\delta\vec{x},\vec{s}+\delta\vec{s})$ (27)

which preserves the dynamics, e.g. such that the time derivative is invariant,
e.g. $\frac{d\vec{x}}{dt}=\frac{d\vec{x}\ ^{\prime}}{dt}$. To see this, let us
equate

$\displaystyle(\mathbb{I}-\chi\mathbf{\Omega}\mathbb{X})^{-1}\mathbf{\Omega}\frac{\vec{s}}{\beta}$
$\displaystyle-$ $\displaystyle\alpha\vec{x}$ (28) $\displaystyle=$
$\displaystyle(\mathbb{I}-\chi\mathbf{\Omega}(\mathbb{X}+\delta\mathbb{X}))^{-1}\mathbf{\Omega}\frac{(\vec{s}+\delta\vec{s})}{\beta}$
$\displaystyle-\alpha(\vec{x}+\delta\vec{x})$

After a bit of algebra, we get the following generalized relationship:

$\displaystyle(\mathbb{I}-\mathbf{\Omega})\delta\vec{x}$ $\displaystyle=$
$\displaystyle 0$ $\displaystyle\delta\vec{s}$ $\displaystyle=$
$\displaystyle\alpha\beta\big{(}\mathbb{I}-\chi(\mathbb{X}+\delta\mathbb{X})\big{)}\delta\vec{x}-\chi\delta\mathbb{X}(\mathbb{I}-\chi\mathbf{\Omega}\mathbb{X})^{-1}\mathbf{\Omega}\vec{s}$
(29) $\displaystyle\ +(\mathbb{I}-\mathbf{\Omega})\vec{k}$

for any vector $\vec{k}$. To see the origin of this generalized invariance,
consider for instance resistances $R_{1},\cdots,R_{n}$ in series with $k$
generators $s_{i}$. We have

$\displaystyle\frac{\sum_{j=1}^{k}s_{j}}{\sum_{j=1}^{n}R_{j}}=i.$ (30)

To preserve the currents, we can either change $R_{j}$ and $s_{j}$ in such a
way that the ratio is preserved. For networks of memristors, eqn. (29)
generalizes this simple invariance.

###  II.6 Existing results on the mean field theories and limitations

The property that we are interested in in this paper is the mean-field
behavior eqn. (26). Interestingly, it has been shown that in the laminar
regime, the dynamics of the whole system can be described in terms of a
dynamical scalar order parameter $\bar{x}$.

The first result of this type was based on a random matrix theory
approximation of the resolvent. Recent studies have provided evidence for
statistical regularities in the resolvent of large matrices, suggesting that
it can be approximated by an effectively one-dimensional matrix that is
universal at the zeroth-order in the limit of weak correlations [59]. This
result has broad applicability across various domains of complexity science
[60]. By defining the vector $\vec{f}$ as $\Omega\vec{x}$, and introducing the
matrix $\tilde{A}$ given by

$\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},$
(31)

we can establish an approximate relation given by

$\displaystyle(I-\chi\mathbf{\Omega}\mathbb{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),$

where $\chi$ is a parameter, $X$ is the matrix of variables $\vec{x}$, and $N$
is the size of the matrix. This approximation allows us to express the
dynamics in terms of the effective one-dimensional variable $x_{cg}$, defined
as the coarse-grained average of $\vec{f}(\vec{x})$, and an operator mean
$\langle\vec{\square}\rangle=N^{-1}\sum_{i=1}^{N}\square_{i}$. The resulting
effective dynamics can be written as

$\displaystyle\frac{d}{d\tau}\bar{x}$ $\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})$ (32) $\displaystyle=$
$\displaystyle-\partial_{x_{cg}}V(x_{cg},\chi)+\mathcal{L}(\vec{x})$

where $\alpha$, $\beta$, and $\vec{S}$ are parameters, and
$\mathcal{L}(\vec{x})$ represents an effective force arising from imperfect
coarse-graining. Notably, this dynamics resembles that of a single memristor,
with the parameter $\frac{s}{\alpha\beta}$ replaced by the mean-field value
$s=\frac{\langle\Omega\vec{S}\rangle}{\alpha\beta}$ at the zeroth-order
approximation. The effective potential above develops two minima when the
parameter $\bar{s}$, defined as the mean of $\vec{s}\beta$ crosses a certain
threshold.

Instead, in [25] the experimental potentiation-depression fit of the two-probe
conductance was based on a different approximation, which is the mean-field
approximation of the matrix inverse which minimizes the Frobenius norm, e.g.
the quantity $\bar{x}$ such that
$\|(I+\chi\mathbf{\Omega}\mathbb{X})^{-1}-(I+\chi\mathbf{\Omega}\bar{x})^{-1}\|^{2}$
is minimized, which is a more brute but effective approximation. The circuit
connectivity was then shown to reabsorbed into effective parameters that can
then be fit a posteriori.

More recently, in [26] it was shown that memristive circuits described by (26)
belong to a larger class of dynamical systems called PEDS (Projective
embedding of dynamical systems) where fixed points of larger systems can be
mapped to fixed points of a lower dynamical systems, using the relationship
$\mathbf{\Omega}^{2}=\mathbf{\Omega}$. However, also there, exact results
could be obtained only for very specific types of projector operators which
are essentially a slight generalization of a mean-field coarse graining
operator.

We see then that in both cases the topological properties of the network
become less important. This should feel to be dissatisfactory, as we do expect
that the topological properties of the network should matter. Thus, in what
follows, we will average while preserving the global properties of the
network, that in our case means preserving the spectral properties of the
matrix $\mathbf{\Omega}$.

##  III Haar measure and averages

###  III.1 von Neumann series and its tensor representation

The key difference between the purely resistive and the purely memristive case
is that the dynamics of the purely memristive case contains an infinite sum
and terms of the form
$\mathscr{O}^{t}\mathbf{\Omega}\mathscr{O}^{t}\mathbb{X}$. To simplify the
calculations that follow, we note that the memristive equation can be written
in terms of the adimensional time variable $\tau=\alpha t$ and $\mathscr{O}$,
and we perform the following change of variable: $\vec{x}=e^{-\tau}\vec{g}$.
This simply removes the term $-\alpha\vec{x}$ in the equation, and allows us
to work temporarily only with the infinite sum. We now write the von Neumann
series for the matrix inverse, in powers of $\chi$, as

$\displaystyle\frac{d\vec{g}}{d\tau}$ $\displaystyle=$
$\displaystyle\sum_{k=0}^{\infty}\chi^{k}e^{k\tau}(\mathscr{O}^{t}\mathbf{\mathbf{\Omega}}\mathscr{O}\mathbb{G})^{k}\vec{s}.$
(33)

It is thus immediate to see that performing the average over the orthogonal
group as we set to do at the beginning of this manuscript is slightly more
complicated than the case of the resistive average.

Let us now perform the following matrix manipulations, using the properties of
tensor products, which we introduce in A. These techniques, used commonly in
quantum information, are borrowed from the theory of linear algebra when
applied to tensor products of linear operators. Let us assume here that
$A_{i}$’s are linear operators on $\mathbb{R}^{E}$. Since we are dealing with
matrices of the form
$(\mathscr{O}^{t}\mathbf{\Omega}\mathscr{O}^{t}\mathbb{X})^{k}$, let us
rewrite this matrix power in a way in which we can perform the average using
well known results.

Let us begin with the simplest non-trivial case, referring to App. A for the
introduction to the tensor product. We have

$A^{2}={\rm Tr}_{2}((A\otimes A)\mathbb{S})$ (34)

Above, the trace is partial, and over the second tensor space. The matrix
$\mathbb{S}$ is the swap operator.

It follows that we can write

$\displaystyle A^{k}$ $\displaystyle=$ $\displaystyle{\rm
Tr}_{2}\big{(}(A^{k-1}\otimes A)\mathbb{S}_{12}\big{)}$ (35) $\displaystyle=$
$\displaystyle{\rm Tr}_{23}((A^{k-2}\otimes A)\otimes
A)\mathbb{S}_{12}\mathbb{S}_{23})$ $\displaystyle\vdots$ $\displaystyle=$
$\displaystyle{\rm Tr}_{2\cdots k}\big{(}(\underbrace{A\otimes\cdots\otimes
A}_{k\ times})\mathbb{S}_{12}\mathbb{S}_{23}\cdots\mathbb{S}_{k-1,k}\big{)}$

where above, with an abuse of notation, we wrote at each step
$\mathbb{S}_{ii+1}$ intending that it swaps only the ith and i-$th+1$ line,
and leaves intact the other lines. Let us call $\mathbb{S}_{2\cdots
k}=\prod_{i=1}^{k-1}\mathbb{S}_{i,i+1}$. Graphically, the identity above can
be seen in Fig. 4, with the swap operators acting on each tensor index.

Figure 4: Graphical representation of the identity of eqn. (35) for $k=4$.

We now note that, then

$\displaystyle(\mathscr{O}^{t}\mathbf{\mathbf{\Omega}}\mathscr{O}\mathbb{G})^{k}$
$\displaystyle=$ $\displaystyle{\rm Tr}_{2\cdots
k}\big{(}(\underbrace{\mathscr{O}^{t}\mathbf{\mathbf{\Omega}}\mathscr{O}\mathbb{G}\otimes\cdots\otimes\mathscr{O}^{t}\mathbf{\mathbf{\Omega}}\mathscr{O}\mathbb{G}}_{k\
times})\mathbb{S}_{2\cdots k}\big{)}.$

Let us call $(\mathscr{O}^{t})^{\otimes
k}=\mathscr{O}^{t}\otimes\cdots\otimes\mathscr{O}^{t}$ and similarly for
$\mathbf{\mathbf{\Omega}}^{\otimes k}$, $\mathscr{O}^{\otimes k}$ and
$\mathbb{G}^{\otimes k}$. We have then that

$\displaystyle(\mathscr{O}^{t}\mathbf{\mathbf{\Omega}}\mathscr{O}\mathbb{G})^{k}$
$\displaystyle=$ $\displaystyle{\rm Tr}_{2\cdots
k}\big{(}(\mathscr{O}^{t})^{\otimes k}\mathbf{\mathbf{\Omega}}^{\otimes
k}\mathscr{O}^{\otimes k}\mathbb{G}^{\otimes k}\mathbb{S}_{2\cdots k}\big{)}.$
(36)

Using this formalism, we then have

$\frac{d\vec{g}}{d\tau}=\sum_{k=0}^{\infty}\chi^{k}e^{k\tau}{\rm Tr}_{2\cdots
k}\big{(}(\mathscr{O}^{t})^{\otimes k}\mathbf{\mathbf{\Omega}}^{\otimes
k}\mathscr{O}^{\otimes k}\mathbb{G}^{\otimes k}\mathbb{S}_{2\cdots
k}\big{)}\vec{s}.$ (37)

The equation above is written in the form of a twirl. A twirl is an expression
of the form [61]

$\displaystyle\langle(\mathscr{O}^{t})^{\otimes
k}\mathbf{\mathbf{\Omega}}^{\otimes k}\mathscr{O}^{\otimes
k}\rangle_{\mathscr{O}}=\mathbb{M}_{k}(\mathbf{\Omega})$ (38)

where $\mathbb{M}_{k}$ is a certain linear operator on
$(\mathbb{R}^{E})^{\otimes k}$ following the average. In tensorial graphical
form, the expression is shown in Fig. 5. The internal index contractions will
be shown as blue lines, and external indices to the twirl shown in red.

We now ask for the following question, e.g. what is the average behavior of
the circuit when we take the isospectral twirling average with respect to
orthogonal operator $\mathscr{O}$? The techniques used to perform these
averages are well established [62, 63, 64, 65].

We have then that

$\big{\langle}\frac{d\vec{g}}{d\tau}\big{\rangle}_{\mathscr{O}}=\sum_{k=0}^{\infty}\chi^{k}e^{k\tau}{\rm
Tr}_{2\cdots k}\big{(}\mathbb{M}_{k}\mathbb{G}^{\otimes k}\mathbb{S}_{2\cdots
k}\big{)}\vec{s}.$ (39)

The goal will be now to calculate the operators $\mathbb{M}_{k}$.

Figure 5: Example of a twirl average in tensorial graphic representation. The
operator $\mathbf{\Omega}^{\otimes k}$ is squeezed in between
$(\mathscr{O}^{t})^{\otimes k}$ and $\mathscr{O}^{\otimes k}$ which is being
averaged over.

###  III.2 Haar measure, Weingarten calculus and mean field theory

In the present manuscript, we are concerned with performing averages of the
form

$\displaystyle\langle\cdot\rangle_{\mathscr{O}}$ (40)

over orthogonal matrices. Let $\mathscr{O}$ be a $E\times E$ orthogonal
matrix, e.g. such that $\mathscr{O}^{t}=\mathscr{O}^{-1}$. The set of
orthogonal matrices forms a group, which we call $\mathbb{O}(E)$. For any
orthogonal matrix $\mathscr{O}$, there exists an inverse
$\mathscr{O}^{-1}=\mathscr{O}^{t}$ such that
$\mathscr{O}\mathscr{O}^{t}=\mathscr{O}^{t}\mathscr{O}=\mathbb{I}$. In fact,
if $\mathscr{O}_{1},\mathscr{O}_{2}$ are orthogonal, it is easy to see that
$\mathscr{O}_{1}\mathscr{O}_{2}$ is also orthogonal, and the identity matrix
is associated to the identity in the group. This is a subspace of
$\mathcal{M}_{N}(\mathbb{C})$, the algebra of $E\times E$ complex matrices
with operation given by the standard matrix multiplication, and which is a
compact topological space. Haar’s theorem [50] states that for any locally
compact Hausdorff topological group, there exists a left-translation-invariant
measure and a right-translation-invariant measure that is unique up to a
positive multiplicative constant. In the case of a compact group, these
measures coincide and are known as the Haar measure. Let $G$ be a group. A
representation $R$ of $\mathscr{O}$ on $\mathcal{H}$ is defined as a group
homomorphism $R:\mathscr{O}\to GL(\mathcal{H})$ [66]. Given $R$, we can
construct the $k$-fold tensor representation of $G$ on $\mathcal{H}^{\otimes
k}$ by acting with $R(g)^{\otimes k}$ for any $g\in G$. It can be observed
that if $R$ is a valid representation, then its $k$-fold tensor product also
forms a representation. The Haar average over the k-fold tensor product can be
performed using the procedure of Weingarten calculus [67, 68, 62], which we
will briefly review here. Consider a bounded operator $A$ on the Hilbert space
$\mathcal{H}$, and let $\mathscr{O}$ be an orthogonal operator chosen
uniformly at random from the orthogonal group $\mathbb{O}(E)$. We are
interested in particular in the isospectral twirling. This is the following
operator. The $2k$-Isospectral twirling of $X$ is defined as [61]:

$\displaystyle\hat{\mathcal{R}}_{\mathbb{O}}^{(2k)}(X):=\int\,d\mathscr{O}\,\mathscr{O}^{{\dagger}\otimes
2k}\left(X^{\otimes 2k}\right)\mathscr{O}^{\otimes 2k}$ (41)

Such an operator arises from the linearization of the expectation values of
polynomials of $\mathscr{O}^{t}X\mathscr{O}$. We see immediately that the
operator $\mathbb{M}_{k}(\mathbf{\Omega})$ in eqn. (38) is in the form of an
isospectral twirling. The Haar average will be computed by the Weingarten
functions method for the orthogonal group. The Haar measure has the properties
$\int d\mathscr{O}=1$ and
$d\mathscr{O}=d(\mathscr{O}\mathscr{O}^{\prime})=d(\mathscr{O}^{\prime}\mathscr{O})$
for any orthogonal matrix $\mathscr{O}^{\prime}\in{\mathbb{O}}(E)$, which are
referred to as the left(right)-invariance of the Haar measure. The general
formula to compute the Haar average is given by [62], and for the orthogonal
group reads:

$\langle\mathscr{O}^{t,\otimes k}A^{\otimes k}\mathscr{O}^{\otimes
k}\rangle_{\mathscr{O}}=\sum_{\sigma,\tau\in\mathcal{P}_{2k}}W^{O}_{g}(\sigma\tau^{-1},E){\rm
Tr}(A^{\otimes k}\mathbb{B}_{\sigma})\mathbb{B}_{\pi}$ (42)

where $W^{O}_{g}(\sigma\tau^{-1},E)$ denotes the Weingarten function of the
group $\mathbb{O}(E)$. For the case of the orthogonal group
$\mathbb{B}_{\sigma}$ and $\mathbb{B}_{\tau}$ are operators associated with
elements of the Brauer algebra $\sigma$ and $\tau$, respectively, and we will
describe them in a moment. Eqn. (42) will be motivated in a moment by showing
explicitly the contractions, and the explicit expression for the Weingarten
functions in terms of the Brauer elements further below.

The Brauer algebra is a mathematical structure that was introduced by Richard
Brauer in the 1930s, and it has connections to both algebra and representation
theory, and is denoted here by $\mathcal{P}_{2k}$. One interesting feature of
the Brauer algebra is its connection to pairings in set theory. In particular,
the Brauer algebra can be used to study pairings of elements from two sets in
a meaningful way. This connection arises from the fact that the Brauer algebra
can be used to represent certain operations on partitions, which are
combinatorial objects that encode the ways in which a set can be divided into
smaller subsets. Specifically, the Brauer algebra can be used to describe
pairings of elements from two sets when the sets have a certain kind of
symmetry. In set theory, a pairing is a function or a relation that maps pairs
of elements from two sets to a third set. The Brauer algebra provides a way to
describe and analyze these pairings in a systematic manner. Examples of Bauer
pairings between $2k$ elements are shown in Fig. 6.

Figure 6: Pairings over 8 elements. Given a set of $2k$ elements, the set of
pairings $\mathcal{P}_{2k}$ is a set of the form
$\tau=\\{\\{i_{1},i_{2}\\},\cdots,\\{i_{2k-1},i_{2k}\\}\\}$, with $i_{a}\neq
i_{b}$. The set can be represented graphically as in the three pairings above.
It is common to divide the set in two set of indices,
$A_{k}=\\{i_{1},\cdots,i_{k}\\},B_{k}=\\{i_{k+1},\cdots,i_{2k}\\}$.

In the following, we will introduce a coloring scheme to understand better
this average. As we can see in Fig. 5, the twirling involves external and
internal indices, the latter being contrated with between $\mathscr{O}^{t}$
$\mathbf{\Omega}$ and $\mathscr{O}$. We will denote internal indices with a
blue color, and external indices with a red color. Another way of expressing
eqn. (42) is by making the indices explicit. The average over the orthogonal
group is given by the formula

$\displaystyle\langle\prod_{r=1}^{2k}\mathscr{O}_{{\color[rgb]{0,0,1}i_{r}}{\color[rgb]{1,0,0}j_{r}}}\rangle_{\mathscr{O}}=\sum_{{\color[rgb]{0,0,1}\sigma},{\color[rgb]{1,0,0}\tau}\in\mathcal{P}_{2k}}\Delta_{{\color[rgb]{0,0,1}\sigma}}({\color[rgb]{0,0,1}\vec{i}})\Delta_{{\color[rgb]{1,0,0}\tau}}({\color[rgb]{1,0,0}\vec{j}})W_{g}^{O}({\color[rgb]{0,0,1}\sigma}{\color[rgb]{1,0,0}\tau^{-1}},E)$

Above, we have the following definition:

$\displaystyle\Delta_{\sigma}(\vec{i})=\prod_{\\{a,b\\}\in\sigma}\delta_{i_{a},i_{b}}.$
(43)

The function above is an index explicit definition of the operators
$\mathbb{B}_{\sigma}^{\prime}s$ involved in eqn. (42). Using the formula
above, the isospectral twirling is written in the form

$\displaystyle(\hat{\mathcal{R}}_{\mathbb{O}}^{(2k)})_{{\color[rgb]{1,0,0}j_{1}\cdots
j_{k}}}^{{\color[rgb]{1,0,0}j_{k+1}\cdots j_{2k}}}(\mathbf{\Omega})$
$\displaystyle=$
$\displaystyle\sum_{{\color[rgb]{0,0,1}i_{1},\cdots,i_{2k}}}\langle\mathscr{O}_{{\color[rgb]{0,0,1}i_{1}}{\color[rgb]{1,0,0}j_{1}}}\cdots\mathscr{O}_{{\color[rgb]{0,0,1}i_{k}}{\color[rgb]{1,0,0}j_{k}}}\mathbf{\Omega}_{{\color[rgb]{0,0,1}i_{1}}{\color[rgb]{0,0,1}i_{k+1}}}\cdots\mathbf{\Omega}_{{\color[rgb]{0,0,1}i_{k}}{\color[rgb]{0,0,1}i_{2k}}}\mathscr{O}_{{\color[rgb]{0,0,1}i_{1}}{\color[rgb]{1,0,0}j_{k+1}}}\cdots\mathscr{O}_{{\color[rgb]{0,0,1}i_{k}}{\color[rgb]{1,0,0}j_{2k}}}\rangle_{\mathscr{O}}$
(44) $\displaystyle=$
$\displaystyle\sum_{{\color[rgb]{0,0,1}i_{1},\cdots,i_{2k}}}\langle\mathscr{O}_{{\color[rgb]{0,0,1}i_{1}}{\color[rgb]{1,0,0}j_{1}}}\cdots\mathscr{O}_{{\color[rgb]{0,0,1}i_{k}}{\color[rgb]{1,0,0}j_{k}}}\mathscr{O}_{{\color[rgb]{1,0,0}i_{1}}{\color[rgb]{0,0,1}j_{k+1}}}\cdots\mathscr{O}_{{\color[rgb]{1,0,0}i_{k}}{\color[rgb]{0,0,1}j_{2k}}}\rangle_{\mathscr{O}}\mathbf{\Omega}_{{\color[rgb]{1,0,0}i_{1}}{\color[rgb]{1,0,0}i_{k+1}}}\cdots\mathbf{\Omega}_{{\color[rgb]{1,0,0}i_{k}}{\color[rgb]{1,0,0}i_{2k}}}$
$\displaystyle=$
$\displaystyle\sum_{{\color[rgb]{0,0,1}i_{1},\cdots,i_{2k}}}\sum_{{\color[rgb]{0,0,1}\sigma},{\color[rgb]{1,0,0}\tau}\in\mathcal{P}_{2k}}\Delta_{{\color[rgb]{0,0,1}\sigma}}({\color[rgb]{0,0,1}\vec{i}})\Delta_{{\color[rgb]{1,0,0}\tau}}({\color[rgb]{1,0,0}\vec{j}})W_{g}^{O}({\color[rgb]{0,0,1}\sigma}{\color[rgb]{1,0,0}\tau^{-1}},E)\mathbf{\Omega}_{{\color[rgb]{0,0,1}i_{1}}{\color[rgb]{0,0,1}i_{k+1}}}\cdots\mathbf{\Omega}_{{\color[rgb]{0,0,1}i_{k}}{\color[rgb]{0,0,1}i_{2k}}}$
$\displaystyle=$
$\displaystyle\sum_{{\color[rgb]{1,0,0}\tau}\in\mathcal{P}_{2k}}\Delta_{{\color[rgb]{1,0,0}\tau}}({\color[rgb]{1,0,0}\vec{j}})\sum_{{\color[rgb]{0,0,1}i_{1},\cdots,i_{2k}}}\sum_{{\color[rgb]{0,0,1}\sigma}\in\mathcal{P}_{2k}}\Delta_{{\color[rgb]{0,0,1}\sigma}}({\color[rgb]{0,0,1}\vec{i}})W_{g}^{O}({\color[rgb]{0,0,1}\sigma}{\color[rgb]{1,0,0}\tau^{-1}},E)\mathbf{\Omega}_{{\color[rgb]{0,0,1}i_{1}}{\color[rgb]{0,0,1}i_{k+1}}}\cdots\mathbf{\Omega}_{{\color[rgb]{0,0,1}i_{k}}{\color[rgb]{0,0,1}i_{2k}}}$
$\displaystyle=$
$\displaystyle\sum_{{\color[rgb]{1,0,0}\tau}\in\mathcal{P}_{2k}}\Delta_{{\color[rgb]{1,0,0}\tau}}({\color[rgb]{1,0,0}\vec{j}})R_{E}({\color[rgb]{1,0,0}\tau},\mathbf{\Omega})$

where we defined

$\displaystyle R_{E}({\color[rgb]{1,0,0}\tau},\mathbf{\Omega})$
$\displaystyle=$
$\displaystyle\sum_{{\color[rgb]{0,0,1}i_{1},\cdots,i_{2k}}}\sum_{{\color[rgb]{0,0,1}\sigma}\in\mathcal{P}_{2k}}\Delta_{{\color[rgb]{0,0,1}\sigma}}({\color[rgb]{0,0,1}\vec{i}})W_{g}^{O}({\color[rgb]{0,0,1}\sigma}{\color[rgb]{1,0,0}\tau^{-1}},E)\cdot$
(45) $\displaystyle\hskip
56.9055pt\cdot\mathbf{\Omega}_{{\color[rgb]{0,0,1}i_{1}}{\color[rgb]{0,0,1}i_{k+1}}}\cdots\mathbf{\Omega}_{{\color[rgb]{0,0,1}i_{k}}{\color[rgb]{0,0,1}i_{2k}}}.$

We see then that one of the summation over the Brauer’s algebra elements is
associated with external indices (red, $\tau$), while the other is associated
with internally contracted indices (blue, $\sigma$). The Weingarten function
depends instead only on the element $\alpha=\sigma\tau^{-1}$. The element
$\alpha$ is essentially the superposition of $\sigma$ and $\tau$, for
instance, in Fig. 6, we simply overlay two Brauer elements on top of each
other. To obtain the Weingarten matrix elements, we first construct the Gram
matrix $Gm(\sigma,\tau)$ as follows [63, 62]. At the order $2k$, let $m$ be
the number of Brauer’ pairings. The Gram matrix is $m\times m$, and in each
corresponding matrix elements associated to $\sigma,\tau$ we have
$E^{c(\sigma,\tau)}$, where $c(\sigma,\tau)$ is the number of closed cycles
obtained by overlaying $\sigma$ and $\tau$ on the same diagram. Then, the
Weingarten matrix $W_{g}^{O}(\sigma\tau^{-1},E)$ is the corresponding element
of the pseudo-inverse of the matrix $Gm(\sigma,\tau)$.

We now have all the elements to perform the average. We need to see the effect
of the two summations and the delta function on the indices of the product of
$\mathbf{\Omega}$’s. Let us first consider the case $k=1$. In this case, there
are only two elements, and thus the only available pairing is
${\color[rgb]{0,0,1}\sigma},{\color[rgb]{1,0,0}\tau}=\\{1,2\\}$.

Since there is only one component, we have $W_{g}^{O}(\sigma,E)=\frac{1}{E}$,
and thus

$\displaystyle R_{E}({\color[rgb]{1,0,0}\tau},\mathbf{\Omega})$
$\displaystyle=\sum_{{\color[rgb]{0,0,1}i_{1},i_{2}}}\frac{\delta_{{\color[rgb]{0,0,1}i_{1},i_{2}}}}{E}\mathbf{\Omega}_{{\color[rgb]{0,0,1}i_{1}i_{2}}}=\frac{1}{E}\sum_{{\color[rgb]{0,0,1}i_{1}}}\mathbf{\Omega}_{{\color[rgb]{0,0,1}i_{1}i_{1}}}=\frac{{\rm
Tr}(\mathbf{\Omega})}{E}$

Thus

$\displaystyle(\hat{\mathcal{R}}_{\mathbb{O}}^{(2)})_{{\color[rgb]{1,0,0}j_{1}}}^{{\color[rgb]{1,0,0}j_{2}}}(\mathbf{\Omega})=\delta_{{\color[rgb]{1,0,0}j_{1}j_{2}}}\frac{{\rm
Tr}(\mathbf{\Omega})}{E}.$ (46)

This is shown in Fig. 7 (top).

Let us now focus on the case $k=2$. This is already not as simple as the case
of $k=1$. The pairing on $4$ elements are now 3:
$\sigma_{1}=\\{\\{1,3\\},\\{2,4\\}\\}$,
$\sigma_{2}=\\{\\{1,4\\},\\{2,3\\}\\}$,
$\sigma_{3}=\\{\\{1,2\\},\\{3,4\\}\\}$. Thus, as we see in Fig. 7 (bottom),
the average is projected over three operators, and we have

$\displaystyle(\hat{\mathcal{R}}_{\mathbb{O}}^{(4)})_{{\color[rgb]{1,0,0}j_{1}j_{2}}}^{{\color[rgb]{1,0,0}j_{3}j_{4}}}(\mathbf{\Omega})$
$\displaystyle=$ $\displaystyle
R_{E}({\color[rgb]{1,0,0}\tau_{1}},\mathbf{\Omega})\delta_{{\color[rgb]{1,0,0}j_{1},j_{3}}}\delta_{{\color[rgb]{1,0,0}j_{2},j_{4}}}$
(47)
$\displaystyle+R_{E}({\color[rgb]{1,0,0}\tau_{2}},\mathbf{\Omega})\delta_{{\color[rgb]{1,0,0}j_{1},j_{4}}}\delta_{{\color[rgb]{1,0,0}j_{2},j_{3}}}$
$\displaystyle+R_{E}({\color[rgb]{1,0,0}\tau_{3}},\mathbf{\Omega})\delta_{{\color[rgb]{1,0,0}j_{1},j_{2}}}\delta_{{\color[rgb]{1,0,0}j_{3},j_{4}}}$

where we now need to calculate

$\displaystyle
R_{E}({\color[rgb]{1,0,0}\tau_{j}},\mathbf{\Omega})=\sum_{{\color[rgb]{0,0,1}i_{1},\cdots,i_{4}}}\sum_{{\color[rgb]{0,0,1}\sigma}\in\mathcal{P}_{2k}}\Delta_{{\color[rgb]{0,0,1}\sigma}}({\color[rgb]{0,0,1}\vec{i}})W_{g}^{O}({\color[rgb]{0,0,1}\sigma}{\color[rgb]{1,0,0}\tau_{j}^{-1}},E)\mathbf{\Omega}_{{\color[rgb]{0,0,1}i_{1}}{\color[rgb]{0,0,1}i_{3}}}\mathbf{\Omega}_{{\color[rgb]{0,0,1}i_{2}}{\color[rgb]{0,0,1}i_{4}}}$
(48)

As we have seen in the previous example, the quantity
$\sum_{{\color[rgb]{0,0,1}i_{1},\cdots,i_{4}}}\Delta_{{\color[rgb]{0,0,1}\sigma}}({\color[rgb]{0,0,1}\vec{i}})$
results in contractions over the indices
${\color[rgb]{0,0,1}i_{j}}^{\prime}s$, and the rest are summed over. These
thus become immediate that these are traces and products of traces of
$\mathbf{\Omega}$. The contractions are shown in Fig. 8. These corresponds to
the elements ${\color[rgb]{0,0,1}\sigma_{1}}\rightarrow{\rm
Tr}(\mathbf{\Omega})^{2}$, ${\color[rgb]{0,0,1}\sigma_{2}}\rightarrow{\rm
Tr}(\mathbf{\Omega}^{2})$, ${\color[rgb]{0,0,1}\sigma_{3}}\rightarrow{\rm
Tr}(\mathbf{\Omega}^{2})$, where we used the fact that $\mathbf{\Omega}$ is a
symmetric matrix. However, these must be multiplied by the orthogonal
Weingarten functions, which are now non-trivial. The table of the Weingarten
coefficients can be obtained by composing ${\color[rgb]{1,0,0}\sigma}$ and
${\color[rgb]{0,0,1}\tau}$. This is simplified by the observation that
$\tau^{-1}=\tau$ since these are pairings. The orthogonal Weingarten function
is the pseudo-inverse of the Gram matrix, whose elements
$Gm({\color[rgb]{1,0,0}\sigma},{\color[rgb]{0,0,1}\tau})=E^{c({\color[rgb]{1,0,0}\sigma},{\color[rgb]{0,0,1}\tau})}$,
where $c({\color[rgb]{1,0,0}\sigma},{\color[rgb]{0,0,1}\tau})$ is the number
of connected components of the graph resulting from the composition of
${\color[rgb]{1,0,0}\sigma},{\color[rgb]{0,0,1}\tau}$. The number of connected
components is given, inspecting Figure 9, by the Gram matrix below

$\displaystyle Gm(\sigma,\tau)=\left(\begin{array}[]{ccc}E^{2}&E&E\\\
E&E^{2}&E\\\ E&E&E^{2}\\\ \end{array}\right)$ (52)

We thus find that, ordering the columns and rows according to
$\tau_{1},\tau_{2},\tau_{3}$ described above, the inverse of the Gram matrix
yields the following Weingarten matrix for $k=2$.

$\displaystyle
W_{g}^{O}({\color[rgb]{0,0,1}\sigma}{\color[rgb]{1,0,0}\tau}^{-1},E)$
$\displaystyle\ \ \ \
=\left(\begin{array}[]{ccc}\frac{E+1}{(E-1)E(E+2)}&-\frac{1}{(E-1)E(E+2)}&-\frac{1}{(E-1)E(E+2)}\\\
-\frac{1}{(E-1)E(E+2)}&\frac{E+1}{(E-1)E(E+2)}&-\frac{1}{(E-1)E(E+2)}\\\
-\frac{1}{(E-1)E(E+2)}&-\frac{1}{(E-1)E(E+2)}&\frac{E+1}{(E-1)E(E+2)}\\\
\end{array}\right)$ (56)

The Weingarten matrix is a symmetric matrix with identical elements on the
diagonal, and for $k=2$ identical elements on the off-diagonal. Let us call
$C_{1}=\frac{E+1}{(E-1)E(E+2)}$ the element on the diagonal, and
$C_{2}=-\frac{1}{E(E-1)(E+2)}$ the elements off-diagonal.

$\displaystyle W_{g}^{O}=\begin{pmatrix}C_{1}&C_{2}&C_{2}\\\
C_{2}&C_{1}&C_{2}\\\ C_{2}&C_{2}&C_{1}\end{pmatrix}$ (57)

Then, we have

$\displaystyle\ \ \ \ \ {\color[rgb]{0,0,1}\sigma_{1}}\hskip
36.98866pt{\color[rgb]{0,0,1}\sigma_{2}}\hskip
42.67912pt{\color[rgb]{0,0,1}\sigma_{3}}$ $\displaystyle R_{1}\equiv
R_{E}({\color[rgb]{1,0,0}\tau_{1}},\mathbf{\Omega})$ $\displaystyle=$
$\displaystyle C_{1}{\rm Tr}(\mathbf{\Omega})^{2}+C_{2}{\rm
Tr}(\mathbf{\Omega}^{2})+C_{2}{\rm Tr}(\mathbf{\Omega}^{2})$ $\displaystyle
R_{2}\equiv R_{E}({\color[rgb]{1,0,0}\tau_{2}},\mathbf{\Omega})$
$\displaystyle=$ $\displaystyle C_{2}{\rm Tr}(\mathbf{\Omega})^{2}+C_{1}{\rm
Tr}(\mathbf{\Omega}^{2})+C_{2}{\rm Tr}(\mathbf{\Omega}^{2})$ $\displaystyle
R_{3}\equiv R_{E}({\color[rgb]{1,0,0}\tau_{3}},\mathbf{\Omega})$
$\displaystyle=$ $\displaystyle C_{2}{\rm Tr}(\mathbf{\Omega})^{2}+C_{2}{\rm
Tr}(\mathbf{\Omega}^{2})+C_{1}{\rm Tr}(\mathbf{\Omega}^{2})$

from which we

$\displaystyle(\hat{\mathcal{R}}_{\mathbb{O}}^{(4)})_{{\color[rgb]{1,0,0}j_{1}j_{2}}}^{{\color[rgb]{1,0,0}j_{3}j_{4}}}(\mathbf{\Omega})$
$\displaystyle=$ $\displaystyle(C_{1}{\rm Tr}(\mathbf{\Omega})^{2}+2C_{2}{\rm
Tr}(\mathbf{\Omega}^{2}))\delta_{{\color[rgb]{1,0,0}j_{1},j_{3}}}\delta_{{\color[rgb]{1,0,0}j_{2},j_{4}}}$
$\displaystyle+\big{(}(C_{1}+C_{2}){\rm Tr}(\mathbf{\Omega}^{2})+C_{2}{\rm
Tr}(\mathbf{\Omega})^{2}\big{)}\delta_{{\color[rgb]{1,0,0}j_{1},j_{4}}}\delta_{{\color[rgb]{1,0,0}j_{2},j_{3}}}$
$\displaystyle+\big{(}(C_{1}+C_{2}){\rm Tr}(\mathbf{\Omega}^{2})+C_{2}{\rm
Tr}(\mathbf{\Omega})^{2}\big{)}\delta_{{\color[rgb]{1,0,0}j_{1},j_{2}}}\delta_{{\color[rgb]{1,0,0}j_{3},j_{4}}}$

Note that we used here the notation in which input indices are at the bottom,
while output indices are at the top. We note that
$\delta_{{\color[rgb]{1,0,0}j_{1},j_{3}}}\delta_{{\color[rgb]{1,0,0}j_{2},j_{4}}}\equiv\mathbb{I}^{\otimes
2}$, $\delta_{{\color[rgb]{1,0,0}j_{1},j_{4}}}\equiv\mathbb{S}$ while
$\delta_{{\color[rgb]{1,0,0}j_{1},j_{2}}}\delta_{{\color[rgb]{1,0,0}j_{3},j_{4}}}\equiv\Pi$,
and information theory $\Pi$ is proportional to the projector onto the
maximally entangled Bell state between the two copies of a Hilbert space
$\mathcal{H}$. In the bra-ket notation, if one defines $|\Phi^{+}\rangle$ as
the Bell state between the two Hilbert spaces,

$\displaystyle|\Phi^{+}\rangle=\frac{1}{\sqrt{E}}\sum_{j_{1}=0,j_{2}=0}^{E-1}|j_{1},j_{1}\rangle\langle
j_{2},j_{2}|,$ (59)

then one has $\Pi=E|\Phi^{+}\rangle\langle\Phi^{+}|$.

Figure 7: Graphical representation of the orthogonal twirling at the first and
second order. At the first, the only contribution is proportional to the the
identity operator. Figure 8: Contractions of the internal lines at the second
order of the twirling. We see that given the pairing $\sigma_{1}$ contributes
a term ${\rm Tr}(\mathbf{\Omega})^{2}$, while $\sigma_{2}$ and $\sigma_{3}$
both contribute ${\rm Tr}(\mathbf{\Omega}^{2})$. Figure 9: The Gram matrix
$Gm$ from which the Weingarten coefficients are derived via a pseudo-inverse.
The element $Gm(\sigma,\tau)$ is the number of loops produced by superposing
$\sigma,\tau$ [63].

###  III.3 Perturbative average

Obtaining exact expressions up to an arbitrary order is rather cumbersome. For
this reason, let us focus on expressions up to order $k=2$. In this case, we
have

$\displaystyle\langle\frac{d\vec{g}}{d\tau}\rangle_{\mathscr{O}}$
$\displaystyle=$ $\displaystyle\vec{s}$ $\displaystyle+$ $\displaystyle\chi
e^{\tau}\langle\mathscr{O}^{t}\mathbf{\mathbf{\Omega}}\mathscr{O}\rangle_{\mathscr{O}}\mathbb{G}\vec{s}$
$\displaystyle+$ $\displaystyle\chi^{2}e^{2\tau}{\rm
Tr}_{2}\big{(}\langle(\mathscr{O}^{t})^{\otimes
2}\mathbf{\mathbf{\Omega}}^{\otimes 2}\mathscr{O}^{\otimes
2}\rangle_{\mathscr{O}}\mathbb{G}^{\otimes 2}\mathbb{S}\big{)}\vec{s}$
$\displaystyle+O(\chi^{3})$

We see in the expression above that in order to perform this calculation, we
need to calculate explicitly

$\displaystyle\langle\mathscr{O}^{t}\mathbf{\mathbf{\Omega}}\mathscr{O}\rangle_{\mathscr{O}}$
$\displaystyle=$ $\displaystyle\mathbb{M}_{1}$
$\displaystyle\langle(\mathscr{O}^{t})^{\otimes
2}\mathbf{\mathbf{\Omega}}^{\otimes 2}\mathscr{O}^{\otimes
2}\rangle_{\mathscr{O}}$ $\displaystyle=$ $\displaystyle\mathbb{M}_{2}$ (60)

Using the methodology described in the previous section, it is easy to see
that $\mathbb{M}_{1}=\frac{{\rm Tr}(\mathbf{\Omega})|}{E}\mathbb{I}$.

We consider the average up to the second order in $\chi$, eqn. (III.3). Let us
note that since $\mathbf{\Omega}^{2}=\mathbf{\Omega}$, we have ${\rm
Tr}(\mathbf{\Omega}^{2})={\rm Tr}(\mathbf{\Omega})$.

From eqn. (LABEL:eqapp:secondorder), we know instead that

$\displaystyle\mathbb{M}_{2}=R_{1}\mathbb{I}^{\otimes
2}+R_{2}\mathbb{S}+R_{3}\Pi.$ (61)

The result of these contractions can be seen in a graphical representation in
Fig. 11. Inserting these expressions into eqn.(LABEL:eq:secondocc2), we obtain

$\displaystyle\langle\frac{d\vec{g}}{d\tau}\rangle_{\mathscr{O}}$
$\displaystyle=$ $\displaystyle(\mathbb{I}+\chi e^{\tau}\frac{{\rm
Tr}(\mathbf{\Omega})}{E}\mathbb{G})\vec{s}+\chi^{2}e^{2\tau}{\rm
Tr}_{2}\big{(}(R_{1}\mathbb{I}^{\otimes
2}+R_{2}\mathbb{S}+R_{3}\Pi)\mathbb{G}^{\otimes
2}\mathbb{S}_{12}\big{)}\vec{s}+O(\chi^{3})$ $\displaystyle=$
$\displaystyle(\mathbb{I}+\chi e^{\tau}\frac{{\rm
Tr}(\mathbf{\Omega})}{E}\mathbb{G})\vec{s}+\chi^{2}e^{2\tau}\Big{(}(R_{1}+R_{3})\mathbb{G}^{2}+R_{2}{\rm
Tr}(\mathbb{G})\mathbb{G}\Big{)}\vec{s}+O(\chi^{3})$ $\displaystyle\approx$
$\displaystyle(\mathbb{I}+\chi e^{\tau}\frac{{\rm
Tr}(\mathbf{\Omega})}{E}\mathbb{G})\vec{s}$ $\displaystyle+$
$\displaystyle\chi^{2}e^{2\tau}\Big{(}((C_{1}+C_{2}){\rm
Tr}(\mathbf{\Omega})^{2}+(3C_{2}+C_{1}){\rm
Tr}(\mathbf{\Omega}))\mathbb{G}^{2}+(C_{2}{\rm
Tr}(\mathbf{\Omega})^{2}+(C_{1}+C_{2}){\rm Tr}(\mathbf{\Omega})){\rm
Tr}(\mathbb{G})\mathbb{G}\Big{)}\vec{s}$ (63)

Using the assumption about the size of the system, e.g. that $E\gg 1$, and the
fact that the graph is loop dense, e.g. ${\rm Tr}(\mathbf{\Omega})=L\gg 1$, we
have $C_{1}\approx\frac{1}{E^{2}}$ and $C_{2}\approx-\frac{1}{E^{3}}$, From
which we get

$\displaystyle\lim_{E\gg 1,L\gg
1}\langle\frac{d\vec{g}}{d\tau}\rangle_{\mathscr{O}}$ $\displaystyle\approx$
$\displaystyle(\mathbb{I}+\chi e^{\tau}\frac{L}{E}\mathbb{G})\vec{s}$ (65)
$\displaystyle+\chi^{2}e^{2\tau}\Big{(}(\frac{L^{2}}{E^{2}}+\frac{L}{E^{2}})\mathbb{G}^{2}$
$\displaystyle\ \ \ +(-\frac{L^{2}}{E^{3}}+\frac{L}{E^{2}}){\rm
Tr}(\mathbb{G})\mathbb{G}\Big{)}\vec{s}$

We now define $\bar{x}(t)=\frac{1}{E}\sum_{j=1}^{E}x_{i}(t)$. It is easy to
see then that $e^{\tau}{\rm Tr}(\mathbb{G})/E=\bar{x}$. Summing the left and
right equation, and performing the mean field replacement
$\vec{s}\rightarrow\bar{s}\vec{1}$, we then obtain the mean-field
approximation, using the fact that $L/E^{2}\rightarrow 0$, we obtain

$\displaystyle\lim_{E\gg 1,L\gg
1}\langle\frac{d\vec{x}}{d\tau}\rangle_{\mathscr{O}}$ $\displaystyle\approx$
$\displaystyle(\mathbb{I}+\chi\frac{L}{E}\mathbb{X})\vec{s}$ (66)
$\displaystyle+$
$\displaystyle\chi^{2}\frac{L^{2}}{E^{2}}\Big{(}\mathbb{X}^{2}-\frac{{\rm
Tr}(\mathbb{X})}{E}\mathbb{X}\Big{)}\vec{s}$

From summing over each index on the left, and dividing by $E$, we obtain

$\displaystyle\frac{d\bar{x}}{dt}\approx\big{(}1+\chi\frac{L}{E}\bar{x}+\chi^{2}\frac{L^{2}}{E^{2}}(\overline{x^{2}}-{\bar{x}}^{2})\big{)}\bar{s}$
(67)

where $\overline{x^{2}}=\frac{1}{E}\sum_{i}x_{i}^{2}$. We introduce
$\text{var}(x)=\overline{x^{2}}-{\bar{x}}^{2}$. Using this definition, we then
have, for $E\gg 1$, that, introducing $\rho=\frac{L}{E}$, that

$\displaystyle\frac{d\bar{x}}{dt}\approx\Big{(}1+\rho\chi\bar{x}+\chi^{2}\rho^{2}\text{var}(x)\Big{)}\bar{s}$
(68)

which is the first equation used in the main text. We see however that this
equation is not closed, as it involves implicitly $\overline{x^{2}}$. We then
note that

$\displaystyle\frac{d}{dt}(\vec{x}^{2})=\frac{d}{dt}\mathbb{X}\vec{x}$
$\displaystyle=$ $\displaystyle 2\mathbb{X}\frac{d}{dt}\vec{x}\approx
2(\mathbb{X}+\chi\frac{L}{E}\mathbb{X}^{2})\vec{s}$ (69) $\displaystyle+$
$\displaystyle 2\chi^{2}\frac{L^{2}}{E^{2}}\Big{(}\mathbb{X}^{3}-\frac{{\rm
Tr}(\mathbb{X})}{E}\mathbb{X}^{2}\Big{)}\vec{s}$ (70)

From which we get, using the mean field approximation again for $\vec{s}$, and
summing cleverly on the left, we get

$\displaystyle\frac{d}{dt}\overline{x^{2}}\approx
2(\bar{x}+\chi\rho\overline{x^{2}})\bar{s}+2\chi^{2}\rho^{2}(\overline{x^{3}}-\overline{x^{2}}\overline{x})\bar{s}$
(71)

which now depends on higher moments, implying a tower of coupled equations. We
can impose closure by assuming that

$\displaystyle\frac{d}{dt}\overline{x^{k}}=0.$ (72)

We will imposing closure at $k=3$, imposing $\overline{x^{3}}=r_{3}$ constant.
Note also that

$\displaystyle\frac{d}{dt}\bar{x}^{2}$ $\displaystyle=$ $\displaystyle
2\bar{x}\frac{d}{dt}\bar{x}\approx
2\bar{x}\Big{(}1+\rho\chi\bar{x}+\chi^{2}\rho^{2}\text{var}(x)\Big{)}\bar{s}$
(73) $\displaystyle=$ $\displaystyle
2\Big{(}\bar{x}+\rho\chi\bar{x}^{2}+\bar{x}\chi^{2}\rho^{2}\text{var}(x)\Big{)}\bar{s}$

and thus, using $\text{var}(x)=\overline{x^{2}}-\bar{x}^{2}$, we get

$\displaystyle\frac{d}{dt}\text{var}(x)$ $\displaystyle=$ $\displaystyle
2(\bar{x}+\chi\rho(\text{var}(x)+\bar{x}^{2}))\bar{s}$ $\displaystyle+$
$\displaystyle
2\chi^{2}\rho^{2}(r_{3}-(\text{var}(x)+\bar{x}^{2})\overline{x})\bar{s}$
$\displaystyle-$ $\displaystyle
2\Big{(}\bar{x}+\rho\chi\bar{x}^{2}+\bar{x}\chi^{2}\rho^{2}\text{var}(x)\Big{)}\bar{s}$
$\displaystyle=$ $\displaystyle
2\rho\chi\bar{s}\Big{(}\text{var}(x)(1-3\bar{x}\rho\chi)+2\rho\chi(r_{3}-\bar{x}^{3})\Big{)}$

Thus, we obtain the set of coupled differential equations

$\displaystyle\frac{d\bar{x}}{dt}$ $\displaystyle\approx$
$\displaystyle\Big{(}1+\rho\chi\bar{x}+\chi^{2}\rho^{2}v\Big{)}\bar{s}$ (75)
$\displaystyle\frac{dv}{dt}$ $\displaystyle\approx$ $\displaystyle
2\rho\chi\bar{s}\Big{(}(1-3\bar{x}\rho\chi)v+2\rho\chi(r_{3}-\bar{x}^{3})\Big{)}$
(76)

The question is whether now these equations have physical fixed points.
However, it is not hard to see that there is no attractive fixed point in the
dynamical system at finite $\bar{x}$ and $v$. This is shown in Fig. 10, where
we plot the phase portrait of the dynamics. Note that the dynamics should be
constrained in the box $\bar{x}\in[0,1]$. Moreover, from the Bhatia-Davis
inequality [69], since $0\leqslant g_{i}\leqslant 1$, we have from the
identity $\text{var}(g)=\frac{1}{E}({\rm Tr}(\mathbb{G}^{2})-{\rm
Tr}(\mathbb{G})^{2})$

$\displaystyle\text{var}(g)\leqslant(1-\bar{g})\bar{g}\leqslant 1.$ (77)

In particular

$\displaystyle\frac{1}{E}\sum_{i}g_{i}^{2}\leqslant\bar{g},$ (78)

or

$\displaystyle{\rm Tr}(\mathbb{G}^{2})\leqslant{\rm Tr}(\mathbb{G}).$ (79)

However, the trajectories are not constrained to this box, and thus this must
imposed in the variance. More importantly, we can see that the only fixed
point is a saddle, which is unphysical as we know that there must be an
attractive fixed point. Modifying the parameters of the equation only modifies
the location of such a fixed point, but not the spectrum of the Jacobian. We
interpret this failure as the necessity to obtain non-perturbative results,
which however we can only obtain using the asymptotics of the Weingarten
calculus.

Figure 10: Phase portrait for the dynamical system of eqns. (75-76), for
$r_{3}=0.4$, $\rho=\chi=0.8$, $s_{0}=0.1$. As we can see there is a saddle
fixed point at an unphysical value of the variance.

###  III.4 Asymptotic regime

As we have seen, the asymptotic results obtained in perturbation theory
expanding in the power of $\chi$ lead to interesting but unsatisfactory
results. In order to obtain non-perturbative results, we are then forced to
consider the asymptotic results in the Weingarten calculus. In particular, we
will use the following result [70]:

$\displaystyle\lim_{E\rightarrow\infty}W_{g}^{O}({\color[rgb]{0,0,1}\sigma}{\color[rgb]{1,0,0}\tau^{-1}},E)=\frac{1}{E^{k}}\delta_{{\color[rgb]{0,0,1}\sigma},{\color[rgb]{1,0,0}\tau}}+O\big{(}\frac{1}{E^{k+1}}\big{)}.$
(80)

This can be intuitively obtained from the definition of the Gram matrix. The
dominant elements of the Gram matrix are in fact on the diagonal, as these
always contribute $E^{k}$, while the off-diagonal elements are of the form
$E,E^{2},\cdots,E^{k-1}$.

Inserting this expression into eqn. (45), we then obtain that [65]

$\displaystyle R^{as}_{E}({\color[rgb]{1,0,0}\tau},\mathbf{\Omega})$
$\displaystyle=$
$\displaystyle\sum_{{\color[rgb]{0,0,1}i_{1},\cdots,i_{2k}}}\sum_{{\color[rgb]{0,0,1}\sigma}\in\mathcal{P}_{2k}}\Delta_{{\color[rgb]{0,0,1}\sigma}}({\color[rgb]{0,0,1}\vec{i}})\frac{1}{E^{k}}\delta_{{\color[rgb]{0,0,1}\sigma},{\color[rgb]{1,0,0}\tau}}\mathbf{\Omega}_{{\color[rgb]{0,0,1}i_{1}}{\color[rgb]{0,0,1}i_{k+1}}}\cdots\mathbf{\Omega}_{{\color[rgb]{0,0,1}i_{k}}{\color[rgb]{0,0,1}i_{2k}}}$
(81) $\displaystyle=$
$\displaystyle\sum_{{\color[rgb]{0,0,1}i_{1},\cdots,i_{2k}}}\Delta_{{\color[rgb]{1,0,0}\tau}}({\color[rgb]{0,0,1}\vec{i}})\frac{1}{E^{k}}\mathbf{\Omega}_{{\color[rgb]{0,0,1}i_{1}}{\color[rgb]{0,0,1}i_{k+1}}}\cdots\mathbf{\Omega}_{{\color[rgb]{0,0,1}i_{k}}{\color[rgb]{0,0,1}i_{2k}}}$
$\displaystyle(\hat{\mathcal{R}}_{\mathbb{O}}^{(2k)})_{{\color[rgb]{1,0,0}j_{1}\cdots
j_{k}}}^{{\color[rgb]{1,0,0}j_{k+1}\cdots j_{2k}}}(\mathbf{\Omega})$
$\displaystyle\hskip
28.45274pt=\sum_{{\color[rgb]{1,0,0}\tau}\in\mathcal{P}_{2k}}\Delta_{{\color[rgb]{1,0,0}\tau}}({\color[rgb]{1,0,0}\vec{j}})R^{as}_{E}({\color[rgb]{1,0,0}\tau},\mathbf{\Omega})$
$\displaystyle\hskip
28.45274pt=\frac{1}{E^{k}}\sum_{{\color[rgb]{1,0,0}\tau}\in\mathcal{P}_{2k}}\Delta_{{\color[rgb]{1,0,0}\tau}}({\color[rgb]{1,0,0}\vec{j}})\sum_{{\color[rgb]{0,0,1}i_{1},\cdots,i_{2k}}}\Delta_{{\color[rgb]{1,0,0}\tau}}({\color[rgb]{0,0,1}\vec{i}})\mathbf{\Omega}_{{\color[rgb]{0,0,1}i_{1}}{\color[rgb]{0,0,1}i_{k+1}}}\cdots\mathbf{\Omega}_{{\color[rgb]{0,0,1}i_{k}}{\color[rgb]{0,0,1}i_{2k}}}.$
(82)

Let us now look at the expressions

$\displaystyle{\rm Tr}_{2\cdots k}\Big{(}\langle\mathscr{O}^{t,\otimes
k}\mathbf{\Omega}^{\otimes k}\mathscr{O}^{\otimes
k}\rangle_{\mathscr{O}}\mathbb{G}^{\otimes k}\mathbb{S}_{2\cdots k}\Big{)}$
$\displaystyle\hskip 56.9055pt={\rm Tr}_{2\cdots
k}\Big{(}\mathbb{M}_{k}\mathbb{G}^{\otimes k}\mathbb{S}_{2\cdots k}\Big{)}$
(83)

Let us now call $\mathbb{B}_{\tau}$ the operator corresponding to the element
$\tau$ of the Brauer algebra associated to the contractions
$\Delta_{{\color[rgb]{1,0,0}\tau}}({\color[rgb]{1,0,0}\vec{j}})$. Using eqn.
(82), we obtain the asymptotic expression

$\displaystyle{\rm Tr}_{2\cdots k}\Big{(}\mathbb{M}_{k}\mathbb{G}^{\otimes
k}\mathbb{S}_{2\cdots k}\Big{)}$ $\displaystyle\hskip
28.45274pt=\frac{1}{E^{k}}\sum_{{\color[rgb]{1,0,0}\tau}\in\mathcal{P}_{2k}}{\rm
Tr}_{2\cdots k}\Big{(}\mathbb{B}_{{\color[rgb]{1,0,0}\tau}}\mathbb{G}^{\otimes
k}\mathbb{S}_{2\cdots k}\Big{)}{\rm
Tr}\Big{(}B_{{\color[rgb]{1,0,0}\tau}}\mathbf{\Omega}^{\otimes k}\Big{)}$

Now, note that $\mathbb{S}_{2\cdots k}$ is associated to an element of the
Brauer algebra, with pairing
$\\{\\{1,2k\\},\\{2,k+1\\},\\{3,k+2\\},\cdots,\\{k,2k-1\\}\\}$. Let us call
this element ${\color[rgb]{1,0,0}\rho}$. Then, we can write the expression
above as

$\displaystyle{\rm Tr}_{2\cdots k}\Big{(}\mathbb{M}_{k}\mathbb{G}^{\otimes
k}\mathbb{S}_{2\cdots k}\Big{)}$ $\displaystyle\hskip
14.22636pt=\frac{1}{E^{k}}\sum_{{\color[rgb]{1,0,0}\tau}\in\mathcal{P}_{2k}}{\rm
Tr}_{2\cdots
k}\Big{(}\mathbb{B}_{{\color[rgb]{1,0,0}\rho\tau}}\mathbb{G}^{\otimes
k}\Big{)}{\rm Tr}\Big{(}B_{{\color[rgb]{1,0,0}\tau}}\mathbf{\Omega}^{\otimes
k}\Big{)}$ (84)

Now, note that we can write

$\displaystyle{\rm
Tr}\Big{(}B_{{\color[rgb]{1,0,0}\tau}}\mathbf{\Omega}^{\otimes
k}\Big{)}=\prod_{s=1}^{c({\color[rgb]{1,0,0}\tau})}{\rm
Tr}(\mathbf{\Omega}^{j_{s}({\color[rgb]{1,0,0}\tau})}),\ \ \
\sum_{s=1}^{c(\tau)}{j_{s}({\color[rgb]{1,0,0}\tau})}=k,$ (85)

where $c({{\color[rgb]{1,0,0}\tau}})$ is the number of connected components of
$\tau$. Since we have $\mathbf{\Omega}^{k}=\mathbf{\Omega}$, the expression
above is maximized when $c({\color[rgb]{1,0,0}\tau})=k$. This is true when
${\color[rgb]{1,0,0}\tau}$ is the identity over the Brauer algebra. As a
result, we have proven the following

######  Proposition 3.

If ${\rm Tr}(\mathbf{\Omega})=L>1$, then $\text{max}_{\tau}{\rm
Tr}\Big{(}B_{\tau}\mathbf{\Omega}^{\otimes k}\Big{)}={\rm
Tr}\Big{(}\mathbf{\Omega}\Big{)}^{k}$, obtained for $\tau=e$, the identitity
in the Brauer algebra.

We now want to argue that if $L(E)$ grows as a function of $E$, then we have

$\displaystyle\lim_{E\rightarrow\infty}{\rm Tr}_{2\cdots
k}\Big{(}\mathbb{M}_{k}\mathbb{G}^{\otimes k}\mathbb{S}_{2\cdots k}\Big{)}$
$\displaystyle=$ $\displaystyle\frac{L^{k}}{E^{k}}{\rm Tr}_{2\cdots
k}\Big{(}\mathbb{B}_{{\color[rgb]{1,0,0}\rho}}\mathbb{G}^{\otimes k}\Big{)}$
(86) $\displaystyle=$
$\displaystyle\frac{L^{k}}{E^{k}}\mathbb{G}^{k}+O(\frac{L^{k-1}}{E^{k}}).$

which is the expression we need in order to derive the large $E$ mean field
theory.

Figure 11: Tensorial representation of the partial trace occurring after the
twirl orthogonal average at the order $\chi^{2}$. As we can see, the first and
third terms give $\mathbb{G}^{2}$, while the middle term contributes ${\rm
Tr}(\mathbb{G})\mathbb{G}$.

We now want to show why this is the dominant contribution for each term at the
order $\chi^{k}$. Let us now discuss quantities of the form ${\rm
Tr}(\mathbb{G}^{j})$. First, note that we have that $|x_{i}|\leqslant 1$, and
thus, since $\alpha\geqslant 0$, we have $|g_{i}|=|e^{-\alpha
t}x_{i}|\leqslant 1$ for $t\geqslant 0$. Note that ${\rm
Tr}(\mathbb{G})=E\bar{g}$. It then means that

$\displaystyle 0\leqslant{\rm Tr}(\mathbb{G})\leqslant E,$ (87)

and thus $0\leqslant\bar{g}\leqslant 1.$ Note that in general, ${\rm
Tr}(\mathbb{G}^{j})=\sum_{i}g_{i}^{j}$.

######  Proposition 4.

Let $\chi<1$, $0\leqslant|\mathbb{X}_{ii}|\leqslant 1$. Then

$\displaystyle\lim_{E\rightarrow\infty}\langle(I-\chi\mathscr{O}^{t}\mathbf{\Omega}\mathscr{O}\mathbb{X})^{-1}\rangle_{\mathscr{O}}=(I-\chi\langle\mathscr{O}^{t}\mathbf{\Omega}\mathscr{O}\rangle_{\mathscr{O}}\mathbb{X})^{-1}.$

###### Proof.

Since both $x_{i}$ and $g_{i}$ are constrained to the interval $[0,1]$, such
result holds for both $\mathbb{G}$ and $\mathbb{X}$. Note that we have, since
$g_{i}\in[0,1]$, $g_{i}^{r}\leqslant g_{i}$, $\forall r\geqslant 1$. It
follows that [71, 72],

$\displaystyle\frac{{\rm Tr}(\mathbb{G}^{r})}{E^{r}}\leqslant\frac{{\rm
Tr}(\mathbb{G})}{E^{r}}\leqslant\frac{\bar{g}}{E^{r-1}}\leqslant\frac{1}{E^{r-1}},$
(88)

since ${\rm Tr}(\mathbb{G})\leqslant E$.

Let us now analyze terms of the form

$\displaystyle{\rm Tr}_{2\cdots
k}\Big{(}\mathbb{B}_{{\color[rgb]{1,0,0}\rho\tau}}\mathbb{G}^{\otimes
k}\Big{)}$ (89)

for $\tau\neq e$. This trace can be written in the form

$\displaystyle\prod_{j=1}^{s}{\rm Tr}(\mathbb{G}^{m_{j}})\mathbb{G}^{m_{0}}$
(90)

where $m_{0}+\sum_{j=1}^{s}m_{j}=k$, where $s$ is the number of closed loops
in the Brauer diagram resulting in the partial trace, and $m_{j}$ is the
number of $\mathbb{G}$ inserted in the loop $j$. Using eqn. (88), we have

$\displaystyle\prod_{j=1}^{s}\frac{{\rm Tr}(\mathbb{G}^{m_{j}})}{E^{m_{j}}}$
$\displaystyle\leqslant$ $\displaystyle\prod_{j=1}^{s}\frac{{\rm
Tr}(\mathbb{G})}{E^{m_{j}}}=\frac{\bar{g}^{s}}{E^{\sum_{j=1}^{s}m_{j}-s}}=\frac{\bar{g}^{s}}{E^{k-m_{0}-s}}$

Now note that since $m_{j}\geqslant 1$, we have $k-m_{0}-s\geqslant 0$. Note
that we assume $\bar{g}<1$, which is the interesting case, as $\bar{g}=1$
means that all memristors are at the fixed point, and the system is not
evolving. As a result, we have shown that every loop contribution, for
$\bar{g}<1$ will be sub-dominant with respect to the identity element in the
Brauer algebra. Since identity is the only term not containing loops, the
result follows in the limit $E\rightarrow\infty$. It then follows that

$\displaystyle\lim_{E\rightarrow\infty}{\rm Tr}_{2\cdots
k}\Big{(}\mathbb{M}_{k}\mathbb{G}^{\otimes k}\mathbb{S}_{2\cdots
k}\Big{)}=\frac{L^{k}}{E^{k}}\mathbb{G}^{k}+O(\frac{c_{k}}{E^{k+1}}).$ (91)

where $c_{k}$ is a constant.

Now, the question is whether we can swap the series with such limit. Note that
the (scalar) series

$\displaystyle\frac{1}{1-x}=\sum_{j=1}^{\infty}x^{j}$ (92)

converges uniformly on $x\in[0,1)$. Similarly, the von Neumann series
$(I-A)^{-1}$ converges uniformly if $\forall\lambda\in\Lambda(A)$,
$\lambda<1$. Note then that since $|x_{i}|\leqslant 1$ and $\chi<1$, we have
$\Lambda(\chi\mathscr{O}^{t}\mathbf{\Omega}\mathscr{O}\mathbb{X})\in
R(0,\chi)$, where $R(0,a)$ is the area in the complex plane of radius $a$ and
centered in zero. Using Proposition 3, the contribution due to the
$\mathbf{\Omega}$ traces is obtained. Then, at order $k$, the dominant
operator contribution is of the form $\mathbb{G}^{k}{\rm
Tr}(\mathbf{\Omega})^{k}$. Going to the $\mathbb{X}$ variables, we have then
obtained the final result

$\displaystyle\lim_{E\rightarrow\infty}\langle(I-\chi\mathscr{O}^{t}\mathbf{\Omega}\mathscr{O}\mathbb{X})^{-1}\rangle_{\mathscr{O}}$
$\displaystyle=$
$\displaystyle\sum_{k=0}^{\infty}\lim_{E\rightarrow\infty}\langle(\chi\mathscr{O}^{t}\mathbf{\Omega}\mathscr{O}\mathbb{X})^{k}\rangle_{\mathscr{O}}$
$\displaystyle=$
$\displaystyle(I-\chi\langle\mathscr{O}^{t}\mathbf{\Omega}\mathscr{O}\rangle_{\mathscr{O}}\mathbb{X})^{-1}\rangle_{\mathscr{O}}$

which proves the proposition. ∎

Figure 12: Effective mean field potential as a function of $x$ and $\rho/\chi$
for the values of (a) $s=0.1$ (b) $s=0.23$ and (c) $s=0.3$. As we can see, for
values of $\rho/\chi$ above $\approx 0.8$, the potential develops two
competing minima as we increase $s$, the effective voltage. For lower values,
instead, the single minimum moves smoothly as a function of $s$.

We use the proposition above to obtain the average dynamics over each cycle
class $C(E,L)$. We have

$\displaystyle\lim_{E\rightarrow\infty}\big{\langle}\frac{d\vec{g}}{d\tau}\big{\rangle}$
$\displaystyle\approx$
$\displaystyle\sum_{k=0}^{\infty}\chi^{k}e^{k\tau}\frac{{\rm
Tr}(\mathbf{\mathbf{\Omega}})^{k}}{E^{k}}{\rm Tr}_{2\cdots
k}\big{(}\mathbb{G}^{\otimes k}\mathbb{S}_{2\cdots k}\big{)}\vec{s}$ (93)
$\displaystyle=$ $\displaystyle\sum_{k=0}^{\infty}\chi^{k}e^{k\tau}\frac{{\rm
Tr}(\mathbf{\mathbf{\Omega}})^{k}}{E^{k}}\mathbb{G}^{k}\vec{s}$
$\displaystyle=$ $\displaystyle\frac{1}{\mathbb{I}-e^{\tau}\frac{{\rm
Tr}(\mathbf{\mathbf{\Omega}})\chi}{E}\mathbb{G}}\vec{s}$

A result which refines those of [44, 59] although for particular type of
random matrices and only in the large $E$ limit.

Because of this, this implies that each memristive device is decoupled from
another one on average, as $\mathbb{G}$ is a diagonal matrix. We can then
analyze the behavior of the system by simply looking independently at each
equation, and we get for the $\langle x_{i}\rangle$ variables that

$\displaystyle\frac{d\langle x_{i}\rangle}{d\tau}$ $\displaystyle=$
$\displaystyle-\langle
x_{i}\rangle+\frac{s_{i}}{\alpha\beta}\frac{1}{1-\chi\frac{L}{E}\langle
x_{i}\rangle}$ (94) $\displaystyle=$ $\displaystyle-\partial_{\langle
x_{i}\rangle}V(\langle x_{i}\rangle)$ (95)

where $L={\rm Tr}(\mathbf{\mathbf{\Omega}})<E$ is the number of cycles in the
graph. We have that $L=E$ if the number memristors is equal to the number of
cycles, and thus if all memristors are completely decoupled. In this case,
$\mathbf{\mathbf{\Omega}}=I$ and in fact equation (95) is exact. Note at this
point that we have the effective mean field potential

$\displaystyle V(\langle x_{i}\rangle)=\frac{1}{2}\langle
x_{i}\rangle^{2}+\frac{s_{i}}{\alpha\beta}\frac{E}{\chi
L}\log(1-\frac{L\chi}{E}\langle x_{i}\rangle)$ (96)

We see the resemblance between eqn. (96) and eqn. (32). The equation is
identical as long as we replace $s_{i}\rightarrow\bar{s}$ and
$\chi\rightarrow\chi\frac{L}{E}$, and thus we can simply analyze the behavior
of the effective potential as a function of voltage as we did in previous
papers.

Let us now make some comments on the effective value of $\chi$, using eqn.
(11) introduced earlier. In the limit $E$ for a connected graph, we obtain
that

$\displaystyle\chi_{eff}=\chi(1-\frac{2}{\bar{d}})$ (97)

where $\bar{d}$ is the average degree. This result allows us to connect the
geometrical properties of the graph to the transition properties. For a
connected graph, $\bar{d}\geqslant 2$. This implies that
$1-\frac{2}{\bar{d}}\leqslant 1$. As a result, we have that if $\bar{d}\propto
V$ as in the case of a complete graph, then in the limit $E\rightarrow\infty$
we have $\chi_{eff}\rightarrow\chi$. For planar graphs instead, we have
$\bar{d}\leqslant 6$ because of Kuratowski’s lemma. It follows that

$\displaystyle 0\leqslant\frac{\chi{eff}}{\chi}\leqslant\frac{2}{3}$ (98)

The effective mean field potential as a function of $\rho/\chi$ and three
values of $s$ are shown in Fig. 12. As we can see, for low value of $s$ the
minimum is located at $x=0$ for all values of $\chi$. However, for larger
values of $s$ the potential has a minimum at intermediate values of $x$. At
large values of $s$, the potential has two competing minima only for values of
$\rho/\chi$ above $0.8$. However, for quasi-planar graphs we know that these
values are not attainable. This provides an explanation for the absence of
first order bulk induced [44] transitions in the effective 2-probe conductance
of silver nanowire experiments [16, 25, 42].

##  IV Conclusions

In conclusion, this manuscript utilized techniques from the orthogonal
Weingarten calculus to derive the mean field theory of memristive systems,
with a specific application in mind of nanowires of low dimensionality. The
results obtained shed light on the behavior of quasi two-dimensional systems,
revealing that they exhibit only a crossover and not a first-order switching
transition in conductance as a bulk phenomenon. This finding provides
important insights into the fundamental physics governing memristive behavior
in nanowires, and suggests that first order transitions are possibly a
boundary phenomenon, e.g. driven by avalanches in switching near the boundary
for each memristive device. Of course, although our results apply to the case
of a simple toy model, the techniques developed here can be used in more
realistic case.

Furthermore, the developed technique has broader applicability in the theory
of neuromorphic devices, as it draws connections with the physics of brain-
like materials. This implies that the derived mean field theory can be
extended to understand and potentially design other neuromorphic devices
beyond nanowires, opening up new possibilities for the development of advanced
electronic devices with brain-inspired functionalities.

The findings presented in this paper contribute to the understanding of
memristive systems and their behavior in quasi two-dimensional systems, while
also highlighting the broader applicability of the developed mean field theory
to the field of neuromorphic devices. These results have the potential to
impact the design and development of future electronic devices with
applications in areas such as artificial intelligence, cognitive computing,
and brain-computer interfaces. Further research and experimental validation of
the proposed theory in various material systems are warranted to fully
comprehend the potential of this approach in advancing the field of memristive
devices and their applications. Additionally, the connection between memory
and ergodicity in brain-inspired devices is worth noting in the context of the
derived mean field theory of memristive systems using cycle classes. The
understanding of how memory is encoded and processed in neuromorphic devices
is crucial for the development of advanced brain-inspired computing systems
[42, 33]. The insights obtained from the developed mean field theory provide
valuable information on the role of boundary-induced transitions such as
critical avalanches and ergodicity.

The findings suggest that the mean field theory can offer a deeper
understanding of the interplay between memory, ergodicity, and the conductance
switching behavior in quasi two-dimensional memristive systems such as
nanowire connectomes. This knowledge can be leveraged to design more efficient
and reliable neuromorphic devices that mimic the memory processing mechanisms
of brain-like devices.

###### Acknowledgements.

The work of FC was carried out under the auspices of the NNSA of the U.S. DoE
at LANL under Contract No. DE-AC52-06NA25396, and in particular support from
LDRD via 20230338ER and 20230627ER. We thank Marco Cerezo, Salvatore Olivero,
Lorenzo Leone and Alioscia Hamma for many useful discussions on the Weingarten
calculus.

## References

  * Gerstner  _et al._ [2014] W. Gerstner, W. M. Kistler, R. Naud, and L. Paninski,  _Neuronal Dynamics: From Single Neurons to Networks and Models of Cognition_ (Cambridge University Press, USA, 2014).
  * Cueva and et. al. [2020] C. J. Cueva and et. al., Low-dimensional dynamics for working memory and time encoding, Proc. Nat. Ac. Sc. 117, 23021 (2020).
  * MacDowell and Buschman [2020] C. J. MacDowell and T. J. Buschman, Low-dimensional spatiotemporal dynamics underlie cortex-wide neural activity, Current Biology 30, 2665 (2020).
  * Mead [1990] C. Mead, Neuromorphic electronic systems, Proceedings of the IEEE 78, 1629 (1990).
  * Yang  _et al._ [2012] J. J. Yang, D. B. Strukov, and D. R. Stewart, Memristive devices for computing, Nature Nanotechnology 8, 13 (2012).
  * Caravelli and Carbajal [2018] F. Caravelli and J. Carbajal, Memristors for the curious outsiders, Technologies 6, 118 (2018).
  * Chua [1971] L. Chua, Memristor-the missing circuit element, IEEE Transactions on Circuit Theory 18, 507 (1971).
  * Chua and Kang [1976] L. Chua and S. M. Kang, Memristive devices and systems, Proceedings of the IEEE 64, 209 (1976).
  * Argall [1968] F. Argall, Switching phenomena in titanium oxide thin films, Solid-State Electronics 11, 535 (1968).
  * Strukov  _et al._ [2008] D. B. Strukov, G. S. Snider, D. R. Stewart, and R. S. Williams, The missing memristor found, Nature 453, 80 (2008).
  * Nagashima  _et al._ [2011] K. Nagashima, T. Yanagida, K. Oka, M. Kanai, A. Klamchuen, J.-S. Kim, B. H. Park, and T. Kawai, Intrinsic mechanisms of memristive switching,  _Nano Letters_ , Nano Letters 11, 2114 (2011).
  * He  _et al._ [2011] L. He, Z.-M. Liao, H.-C. Wu, X.-X. Tian, D.-S. Xu, G. L. W. Cross, G. S. Duesberg, I. V. Shvets, and D.-P. Yu, Memory and threshold resistance switching in ni/nio core–shell nanowires,  _Nano Letters_ , Nano Letters 11, 4601 (2011).
  * Di Ventra and Pershin [2013] M. Di Ventra and Y. V. Pershin, The parallel approach, Nature Physics 9, 200 (2013).
  * Diaz-Alvarez  _et al._ [2019] A. Diaz-Alvarez, R. Higuchi, P. Sanz-Leon, I. Marcus, Y. Shingaya, A. Z. Stieg, J. K. Gimzewski, Z. Kuncic, and T. Nakayama, Emergent dynamics of neuromorphic nanowire networks, Scientific Reports 9, 14920 (2019).
  * Kuncic and Nakayama [2021] Z. Kuncic and T. Nakayama, Neuromorphic nanowire networks: principles, progress and future prospects for neuro-inspired information processing, Advances in Physics: X 6, 10.1080/23746149.2021.1894234 (2021).
  * Ohno  _et al._ [2011] T. Ohno, T. Hasegawa, A. Nayak, T. Tsuruoka, J. K. Gimzewski, and M. Aono, Sensory and short-term memory formations observed in a ag2s gap-type atomic switch, Applied Physics Letters 99, 203108 (2011).
  * Wang  _et al._ [2016] Z. Wang, S. Joshi, S. E. Savel’ev, H. Jiang, R. Midya, P. Lin, M. Hu, N. Ge, J. P. Strachan, Z. Li, Q. Wu, M. Barnell, G.-L. Li, H. L. Xin, R. S. Williams, Q. Xia, and J. J. Yang, Memristors with diffusive dynamics as synaptic emulators for neuromorphic computing, Nature Materials 16, 101 (2016).
  * Milano  _et al._ [2020] G. Milano, G. Pedretti, M. Fretto, L. Boarino, F. Benfenati, D. Ielmini, I. Valov, and C. Ricciardi, Brain-inspired structural plasticity through reweighting and rewiring in multi-terminal self-organizing memristive nanowire networks, Advanced Intelligent Systems 2, 2000096 (2020).
  * Caravelli [2019] F. Caravelli, Asymptotic behavior of memristive circuits, Entropy 21, 789 (2019).
  * Hopfield and Tank [1986] J. J. Hopfield and D. W. Tank, Neural computation of decisions in optimization problems, Biological Cybernetics 52 (1986).
  * Caravelli  _et al._ [2017] F. Caravelli, F. L. Traversa, and M. Di Ventra, The complex dynamics of memristive circuits: Analytical results and universal slow relaxation, Physical Review E 95, 022140 (2017).
  * Caravelli [2017a] F. Caravelli, The mise en scéne of memristive networks: effective memory, dynamics and learning, International Journal of Parallel, Emergent and Distributed Systems 33, 350 (2017a).
  * Zegarac and Caravelli [2019] A. Zegarac and F. Caravelli, Memristive networks: From graph theory to statistical physics, EPL (Europhysics Letters) 125, 10001 (2019).
  * Caravelli [2017b] F. Caravelli, Locality of interactions for planar memristive circuits, Physical Review E 96, 052206 (2017b).
  * Caravelli  _et al._ [2023] F. Caravelli, G. Milano, C. Ricciardi, and Z. Kuncic, Mean field theory of self-organizing memristive connectomes, to appear in Annalen der Physik (2023).
  * Caravelli  _et al._ [2022] F. Caravelli, F. L. Traversa, M. Bonnin, and F. Bonani, Projective embedding of dynamical systems: uniform mean field equations, to appear in Physica D, arXiv:2201.02355 (2022).
  * Packard [1988] N. H. Packard, Dynamic patterns in complex systems, in Dynamic Patterns in Complex Systems, edited by J. A. S. Kelso, A. J. Mandell, and M. F. Shlesinger (World Scientific, Singapore) , 293 (1988).
  * Langton [1990] C. Langton, Computation at the edge of chaos: Phase transitions and emergent computation, Physica D 42, 12–37 (1990).
  * Jaeger [2001] H. Jaeger, The echo state approach to analysing and training recurrent neural networks, German National Research Center for Information Technology (GMD) Technical Report 148, 34 (2001).
  * Natschlaeger  _et al._ [2002] T. Natschlaeger, W. Maass, and H. Markram, The liquid computer: a novel stratey for real-time copmuting on time-series, Special Issue on Foundations of Information Processing of TELEMATIK 8, 39–43 (2002).
  * Carroll [2020] T. L. Carroll, Do reservoir computers work best at the edge of chaos?, Chaos 30 (2020).
  * Sheldon  _et al._ [2022] F. C. Sheldon, A. Kolchinsky, and F. Caravelli, Computational capacity of $lrc$, memristive, and hybrid reservoirs, Phys. Rev. E 106, 045310 (2022).
  * Loeffler  _et al._ [2023] A. Loeffler, A. Diaz-Alvarez, R. Zhu, N. Ganesh, J. M. Shine, T. Nakayama, and Z. Kuncic, Neuromorphic learning, working memory, and metaplasticity in nanowire networks, Science Advances 9, eadg3289 (2023), https://www.science.org/doi/pdf/10.1126/sciadv.adg3289 .
  * Baccetti  _et al._ [2023] V. Baccetti, R. Zhu, F. Caravelli, and Z. Kuncic, to appear (2023).
  * Stassinopoulos and Bak [1995] D. Stassinopoulos and P. Bak, Democratic reinforcement: a principle for brain function, Phys Rev E 51, 10.1103/physreve.51.5033 (1995).
  * Chialvo and Bak [1999] D. Chialvo and P. Bak, Learning from mistakes, Neuroscience 90, 1137 (1999).
  * Bak and Chialvo [2001] P. Bak and D. Chialvo, Adaptive learning by extremal dynamics and negative feedback, Phys Rev E 63 (2001).
  * Chialvo [2010] D. Chialvo, Emergent complex neural dynamics, Nature Physics 6, 744 (2010).
  * Carbajal  _et al._ [2022] J. P. Carbajal, D. A. Martin, and D. R. Chialvo, Learning by mistakes in memristor networks, Phys. Rev. E 105, 054306 (2022).
  * Jensen [1998] H. J. Jensen,  _Self-Organized Criticality: Emergent Complex Behavior in Physical and Biological Systems_, Cambridge Lecture Notes in Physics (Cambridge University Press, 1998).
  * Avizienis  _et al._ [2012] A. Avizienis, H. Sillin, C. Martin-Olmos, H. Shieh, M. Aono, A. Stieg, and J. Gimzewski, Neuromorphic atomic switch networks, PLoS One 7 (2012).
  * Hochstetter  _et al._ [2021] J. Hochstetter, R. Zhu, A. Loeffler, A. Diaz-Alvarez, T. Nakayama, and Z. Kuncic, Avalanches and edge-of-chaos learning in neuromorphic nanowire networks, Nature Communications 12, 4008 (2021).
  * Mallinson  _et al._ [2019] J. B. Mallinson, S. Shirai, S. K. Acharya, S. K. Bose, E. Galli, and S. A. Brown, Avalanches and criticality in self-organized nanoscale networks, Science Advances 5, eaaw8438 (2019), https://www.science.org/doi/pdf/10.1126/sciadv.aaw8438 .
  * Caravelli  _et al._ [2021] F. Caravelli, F. C. Sheldon, and F. L. Traversa, Global minimization via classical tunneling assisted by collective force field formation, Science Advances 7, 10.1126/sciadv.abh1542 (2021).
  * Caravelli and Barucca [2018] F. Caravelli and P. Barucca, A mean-field model of memristive circuit interaction, EPL (Europhysics Letters) 122, 40008 (2018).
  * Dowling and Pershin [2022] V. J. Dowling and Y. V. Pershin, Memristive ising circuits, Phys. Rev. E 106, 054156 (2022).
  * Bollobás [1998] B. Bollobás,  _Modern Graph Theory_ (Springer New York, 1998).
  * Ellens  _et al._ [2011] W. Ellens, F. Spieksma, P. Van Mieghem, A. Jamakovic, and R. Kooij, Effective graph resistance, Linear Algebra and its Applications 435, 2491 (2011).
  * Note [1] Details about this construction can be found in [21].
  * Haar [1933] A. Haar, Der massbegriff in der theorie der kontinuierlichen gruppen, Annals of Mathematics 34(2) (1933).
  * Joglekar and Wolf [2009] Y. N. Joglekar and S. J. Wolf, The elusive memristor: properties of basic electrical circuits, Eur. J. of Phys. 30, 661 (2009).
  * Biolek and et. al. [2013] D. Biolek and et. al., Some fingerprints of ideal memristors, in  _2013 IEEE Int. Symp. on Circ. and Sys.)_ (IEEE, 2013).
  * Prodromakis and et. al. [2011] T. Prodromakis and et. al., A versatile memristor model with nonlinear dopant kinetics, IEEE Trans. on El. Dev. 58, 3099 (2011).
  * Ginoux and et. al. [2020] J.-M. Ginoux and et. al., A physical memristor based muthuswamy–chua–ginoux system, Sci. Rep. 10, 10.1038/s41598-020-76108-z (2020).
  * Ascoli and et. al. [2013] A. Ascoli and et. al., Memristor model comparison, IEEE Circ. and Sys. Mag. 13, 89 (2013).
  * Corinto  _et al._ [2012] F. Corinto, A. Ascoli, and M. Gilli, Analysis of current-voltage characteristics for memristive elements in pattern recognition systems, Int. J. of Circ. Th. and App. 40, 1277 (2012).
  * Corinto and Ascoli [2012] F. Corinto and A. Ascoli, A boundary condition-based approach to the modeling of memristor nanostructures, IEEE Trans. on Circ. and Sys. 59, 2713 (2012).
  * Ascoli  _et al._ [2015] A. Ascoli, F. Corinto, and R. Tetzlaff, Generalized boundary condition memristor model, Int. J. of Circ. Th. and App. 44, 60 (2015).
  * Bartolucci  _et al._ [2020a] S. Bartolucci, F. Caccioli, F. Caravelli, and P. Vivo, Inversion-free leontief inverse: statistical regularities in input-output analysis from partial information, arXiv:2009.06350 (2020a).
  * Bartolucci  _et al._ [2020b] S. Bartolucci, F. Caccioli, F. Caravelli, and P. Vivo, Universal rankings in complex input-output organizations, to appear in Phys. Rev. Research, arXiv.org:2009.06307 (2020b).
  * Oliviero  _et al._ [2021] S. Oliviero, L. Leone, F. Caravelli, and A. Hamma, Random matrix theory of the isospectral twirling, SciPost Phys. 10 (2021).
  * Collins and Śniady [2006] B. Collins and P. Śniady, Integration with respect to the Haar measure on unitary, orthogonal and symplectic group, Communications in Mathematical Physics 264, 773 (2006).
  * Zinn-Justin [2010] P. Zinn-Justin, Jucys–murphy elements and weingarten matrices, Letters in Mathematical Physics 91, 119 (2010).
  * Collins and Matsumoto [2017] B. Collins and S. Matsumoto, Weingarten calculus via orthogonality relations: new applications, arXiv preprint arXiv:1701.04493 (2017).
  * Banica [2010] T. Banica, The orthogonal weingarten formula in compact form, Letters in Mathematical Physics 91, 105 (2010).
  * Fulton and Harris [1991] J. Fulton and W. Harris,  _Representation Theory: A First Course_ (Springer, New York, 1991).
  * Weingarten [1978] D. Weingarten, Asymptotic behavior of group integrals in the limit of infinite rank, Journal of Mathematical Physics 19, 999 (1978), https://doi.org/10.1063/1.523807 .
  * Collins [2003] B. Collins, Moments and cumulants of polynomial random variables on unitary groups, the Itzykson-Zuber integral, and free probability, International Mathematics Research Notices 2003, 953 (2003).
  * Bhatia and Davis [2000] R. Bhatia and C. Davis, A better bound on the variance, The American Mathematical Monthly 107, 353 (2000), https://doi.org/10.1080/00029890.2000.12005203 .
  * Collins  _et al._ [2009] B. Collins, A. Guionnet, and E. Maurel-Segala, Asymptotics of unitary and orthogonal matrix integrals, Advances in Mathematics 222, 172 (2009).
  * Hardy  _et al._ [1934] G. H. Hardy, J. E. Littlewood, and G. Pólya,  _Inequalities_ (1934).
  * Bernstein [2018] D. S. Bernstein,  _Scalar, Vector, and Matrix Mathematics: Theory, Facts, and Formulas - Revised and Expanded Edition_ (PRINCETON UNIV PR, 2018).

##  Appendix A Brief introduction to tensor products

A matrix, or linear operator $M_{ij}$ will be denoted in the following as
$M_{j}^{i}$, to make clear which indices are input and which indices are
output. Matrix multiplication is then denoted by
$\sum_{j}M_{ij}N_{jk}\equiv\sum_{j}M^{i}_{j}N^{j}_{k}$, indicating that
contractions occur between bottom and upper indices respectively. Matrix
products between matrices, in general, can be written as

$(A_{1}A_{2})_{ij}=\sum_{k}(A_{1})^{i}_{k}(A_{2})^{k}_{j}=\sum_{kt}(A_{1})^{i}_{k}(A_{2})^{k}_{j}\delta_{kt}$
(99)

The reason for the introduction of this notation will be clear in a moment,
when we introduce tensor products.

We can write the identity above in a different way. We introduce the tensor
product of the operator $A_{1}$ and $A_{2}$, and define

$A_{1}\otimes
A_{2}:\mathbb{R}^{E}\otimes\mathbb{R}^{E}\rightarrow\mathbb{R}^{E}\otimes\mathbb{R}^{E}.$
(100)

Thus, instead of acting on a single copy of $\mathbb{R}^{E}$ it acts on two
copies, $\mathbb{R}^{E}\otimes\mathbb{R}^{E}\equiv\mathbb{R}^{d,\otimes 2}$
(we will assume $d$ is the dimensionality of the single linear space from now
on) e.g.

$(A_{1}\otimes
A_{2})(\vec{v}_{1}\otimes\vec{v}_{2})=(A_{1}\vec{v}_{1})\otimes(A_{2}\vec{v}_{2}).$
(101)

Any linear operator on $\mathcal{R}^{d,\otimes n}$ can be written as

$\displaystyle A^{i_{1}\cdots i_{n}}_{j_{1}\cdots j_{n}}$ $\displaystyle=$
$\displaystyle\sum_{j_{1}\cdots j_{n},k_{1}\cdots k_{n}=1}^{n}c^{j_{1}\cdots
j_{n}}_{k_{1}\cdots k_{n}}A_{j_{1}k_{1}}\otimes\cdots\otimes A_{j_{1}k_{n}}$
(102) $\displaystyle=$ $\displaystyle\sum_{j_{1}\cdots
j_{n}=1}^{n^{2}}\tilde{c}_{j_{1}\cdots
j_{n}}\tilde{A}_{j_{1}}\otimes\cdots\otimes\tilde{A}_{j_{n}}$

where $\tilde{A}_{j_{i}}$ are a basis for $GL(E)$. We used the notation in
which the lower index is contracted with the vector’s index. Thus

$\displaystyle
A(\vec{v}_{1}\otimes\cdots\otimes\vec{v}_{n})=\sum_{j_{1},\cdots,j_{n}}A^{i_{1}\cdots
i_{n}}_{j_{1}\cdots j_{n}}(v_{1}^{j_{1}}\otimes\cdots\otimes v_{n}^{j_{n}})$

As made explicit above, we use the notation for a tensor the bottom indices
are input and the top indices are output.

We now introduce the swap operator on $\mathbb{R}^{E}\otimes\mathbb{R}^{E}$,
$\mathbb{S}$, or $\mathbb{S}_{12}$ to introduce a notation that will be clear
in a moment. The operator $\mathbb{S}_{12}$ performs the following action

$\mathbb{S}(\vec{v}_{1}\otimes\vec{v}_{2})=\vec{v}_{2}\otimes\vec{v}_{1}$
(104)

These operations can be written in graphical terms. Tensor products can be
written in terms of lines and boxes, as in Fig. 13. The last operation we wish
to discuss is the partial trace. Consider a tensor $A^{i_{1}\cdots
i_{n}}_{j_{1}\cdots j_{n}}$ acting on $\mathbb{R}^{\otimes n}$. A partial
trace ${\rm Tr}_{k}(\cdot):GL(E)^{\otimes n}\rightarrow GL(E)^{\otimes n-1}$
is the operation ${\rm Tr}_{k}(A^{i_{1}\cdots i_{n}}_{j_{1}\cdots
j_{n}})=\sum_{n}A^{i_{1}\cdots i_{k-1}ni_{k+1}i_{n}}_{j_{1}\cdots
j_{k-1}nj_{k+1}\cdots j_{n}}$. This is a generalization of the trace of a
matrix, in which ${\rm Tr}(A_{ij})={\rm Tr}(A^{i}_{j})=\sum_{n}A_{nn}$. The
partial trace can be generalized to multiple index contractions, depending on
the subspace.

These operations are shown in Fig. 13 in a graphical manner. Tensor products
can be represented as lines and boxes, while swap operators simply line
inversions. The partial or full trace contracts lines from the left to the
equivalent line on the diagram. Contracted lines simply mean summing over the
indices of that particular line. For instance, ${\rm Tr}(A)$ is shown in Fig.
14.

In the main text, we use the combined swap operator

$\displaystyle S_{2\cdots k}=S_{1,2}S_{2,3}S_{3,4}\cdots S_{k-1,k}$ (105)

to construct

$\displaystyle A^{k}=\text{Tr}_{2\cdots k}(A^{\otimes k}S_{2\cdots k}).$ (106)

Such identity, for $k=4$, is shown graphically in Fig. 4.

Figure 13: Graphical representations of the identities involving partial
traces over tensor products of linear operators, and the swap operator. Figure
14: A linear operator $A$ and the graphical representation of its trace.