# On the Degeneracy of Spin Ice Graphs, and Its Estimate via the Bethe
Permanent

Francesco Caravelli  Theoretical Division (T4), Los Alamos National
Laboratory, Los Alamos, New Mexico 87545, USA  Michael Saccone  Center for
Nonlinear Studies, Los Alamos National Laboratory, Los Alamos, New Mexico
87545, USA  Theoretical Division (T4), Los Alamos National Laboratory, Los
Alamos, New Mexico 87545, USA  Cristiano Nisoli  Theoretical Division (T4),
Los Alamos National Laboratory, Los Alamos, New Mexico 87545, USA

###### Abstract

The concept of spin ice can be extended to a general graph. We study the
degeneracy of spin ice graph on arbitrary interaction structures via graph
theory. Via the mapping of spin ices to the Ising model, we clarify whether
the inverse mapping is possible via a modified Krausz construction. From the
gauge freedom of frustrated Ising systems, we derive exact, general results
about frustration and degeneracy. We demonstrate for the first time that every
spin ice graph, with the exception of the 1D Ising model, is degenerate. We
then study how degeneracy scales in size, using the mapping between Eulerian
trails and spin ice manifolds, and a permanental identity for the number of
Eulerian orientations. We show that the Bethe permanent technique provides
both an estimate and a lower bound to the frustration of spin ices on
arbitrary graphs of even degree. While such technique can be used also to
obtain an upper bound, we find that in all the examples we studied but one,
another upper bound based on Schrijver inequality is tighter.

##  I Introduction

Ever since the discovery of degeneracy of ground states in constrained,
disordered systems obeying the so-called ice rule bernal ; pauling ; lieb  and
their subsequent experimental implementation in magnetic systems called Spin
Ices ramirez , there has been an active interest in frustrated materials.
Recently, the idea has been extended both theoretically and experimentally to
artificial realizations called artificial spin ices (ASI) colloq  which have
allowed for design of frustration to generate exotic behaviors in their
collective physics. ASIs are arrays of interacting, shape-anisotropic nano-
islands, each of which can be modeled as a binary Ising spin. They can be
characterized directly in real time and real space via a variety of techniques
Wang1 ; Sandra ; Bader .

A spin ice can be described abstractly as a set of binary spins arranged on
the edge of a lattice, such that its low energy configuration obeys the ice
rule. This rule dictates that for each vertex the absolute difference between
spins pointing toward the vertex and spins pointing out of the vertex is zero
if the vertex has even coordination, or one if the vertex has odd
coordination.

Recently, there has been an intense investigation both in the physical
underpinnings and control of artificial spin ices and their emergent
interactions reddim ; Heyderman ; Canals1 ; Nisoli1 ; Morgan ; Budrikis ;
Branford ; Ryzhkin ; Moeller ; Chern2 ; Le ; Chern3 ; Gliga ; Nisoli4 ;
Gilbert2 ; Bhat ; CaravelliMASI , with a broad interest in applications,
ranging from topological order Castelnovo1 ; topor , memory in materials
Lammert2 ; GilbertMem , disordered systems and slow relaxation Cugliandolo2 ,
novel resistive switching memristors ; caravelli , and embedding logic
circuits in the magnetic substrate logic2 ; logic3 ; logic4 ; logic5 . These
new materials are highly controllable gartside ; WangYL2 ; WangYL ; colloq ;
Nisoli4 ; vavassori0 ; vavassori  and can be used to realize novel models via
engineering geometric frustration Nisoli4 ; Schanilec ; Nisoli8 ; mol .

The collective behavior of these artificial structures typically depends on
the geometry, which is open to design Morrison ; Nisoli4  via lithographic
printing. As novel lithographic techniques are discovered, the control of the
dimensionality of these materials require new theoretical tools to understand
the frustration of non-planar artificial spin ices Ladak3d ; Ladak3d2 .
Moreover, it is now possible to embed general spin ice graphs into quantum
annealers QuantumASI . For these reasons, this paper is focused on a more
theoretical and extensible approach for calculating lower and upper bounds to
the degeneracy of the ground state of generic spin ices.

Previous work Field  set up the concept of spin ice on a general graph. Spin
ice concepts are often translatable into the language of graph theory, and
vice versa. For instance in the mathematical literature a balanced graph (with
even degree nodes) is a directed graph whose indegrees is equal to the
outdegrees Graph1 . In spin ice language, this corresponds to a configuration
of the ice manifold Field . Thus, the problem of finding the Pauling entropy
pauling  of a spin ice graph is equivalent to the problem of counting the
number of balanced digraphs. Moreover, it is one of the many celebrated
results by Euler that only graphs which can be balanced via an orientation
support an Eulerian trail Euler1 . In this sense, many results from graph
theory can be borrowed to study frustration in spin ice.

In this paper we use and generalize some of the known results from graph
theory. We demonstrate for the first time that a general spin ice graph is
always degenerate, with the exception of the trivial case: the one-dimensional
Ising model. There, because the scaling of degeneracy with size is fundamental
to the notion of Pauling entropy in spin ices, we compute lower and upper
bound the ground state degeneracy of several spin ices.

In the first part of the paper, we use the Line Graph dual representation
graph , and derive properties for the effective Ising model on arbitrary
graphs. In particular, we use the Krausz coarse-graining procedure to show
that from the effective Ising model there is a well-known technique to obtain
the original spin ice interaction graph.

In the second part of the paper, we focus on estimation techniques. For even-
degree graphs there exists a permanental identity which in principle, but not
in practice, allows for the evaluation of the degeneracy of the spin ice.
However, using the Bethe Permanent and the Schrijver bound, it is possible to
obtain lower and upper bounds to geometric frustration for the cases in which
the degree of the graph is even. We apply these bounds to the square lattice,
the triangular tiling, the cubic lattice, and tournaments.

##  II Graph theory for Spin Ice

General graph theoretic approaches are common tools in Statistical Physics
Baxter ; Fisher ; CaravelliMarkopoulou . Previously, one of us has discussed
the notion of spin ice on a general undirected graph $\mathcal{G}$, where a
spin configuration may be thought of as a directed graph, considered for its
Coulomb phase properties, and shown that charge correlations are computable
from graph spectral analysis Field . We use the same approach here, but unlike
the previous approximated study, we derive exact results for the degeneracy of
the ice manifold. Later in the paper we will also introduce a way to obtain
also estimates (which are lower bounds), but here we prefer to keep a general
discussion.

Consider a set of spins $s_{j}$ lying on the bonds of a graph, and the
Hamiltonian

$H=J\sum_{v}(\sum_{j\rightarrow v}\pm s_{j})^{2},$ (1)

where $v$ are the vertices of the graph $\mathcal{G}$, and $s_{j}$ have a
certain orientation, such that the minimum of the energy corresponds to
minimal absolute value of charge, defined as the difference between vertices
pointing in or out.

To see the connection between frustrated spins in the system consider
$e=|E(\mathcal{G})|$ the total number of edges (ordered pairs of adjacent
vertices) of the graph graph , while $n=|V(\mathcal{G})|$ is the total number
of nodes. Let us first introduce a few graph theoretical constructs in order
to fix the notation. For a generic and undirected graph $\mathcal{G}$ consider
the (undirected) incidence matrix $B$ of size $n\times e$ with entries
$B_{i\beta}$, where $i$ is an integer between $1$ and $n$ on the set of
vertices and $\beta$ is an integer between $1$ and $e$ on the set of edges,
such that:

$B_{i\beta}=\left\\{\begin{array}[]{rl}1&\text{if the edge $\beta$ contains
the vertex $i$},\\\ 0&\text{otherwise }.\end{array}\right.$ (2)

Let us now consider instead the directed incidence matrix $B_{i\beta}$, of
size $n\times p$ constructed as follows. First, assign an orientation to each
edge of the graph. This can be thought as a possible spin configuration Field
. Given such orientation $\mathcal{O}$, we assign the matrix elements of
$B^{\mathcal{O}}_{v\beta}$ as

$B_{v\beta}=\begin{cases}0&\text{if the edge }\beta\text{ is not incident to
the vertex }v\\\ 1&\text{if the edge, given the orientation $\mathcal{O}$,
enter $v$}\\\ -1&\text{if the edge, given the orientation $\mathcal{O}$,
leaves $v$}\end{cases}$ (3)

(we use latin indices for vertices and greek for edges or spins).

Crucially, we can rewrite the Hamiltonian for a generic spin ice as

$H=J\sum_{v=1}^{n}(\sum_{\beta=1}^{p}B_{v,\beta}s_{\beta})^{2}=J\sum_{v=1}^{n}\sum_{\beta,\beta^{\prime}=1}^{e}pB_{v,\beta}B_{v,\beta^{\prime}}s_{\beta}s_{\beta^{\prime}}$
(4)

Swapping the vertex and spin summation, we write

$H=J\sum_{\beta,\beta^{\prime}=1}^{p}Q_{\beta,\beta^{\prime}}s_{\beta}s_{\beta^{\prime}},$
(5)

where
$Q_{\beta,\beta^{\prime}}=\sum_{v}B_{v,\beta}B_{v,\beta^{\prime}}\equiv(B^{t}B)_{\beta,\beta^{\prime}}$
is symmetric.

We now note the following. Given the directed incidence matrix, we have

$Q_{\beta,\beta^{\prime}}=\begin{cases}2&:\text{if }\beta=\beta^{\prime}\\\
0&:\text{if the spins $\beta$ and $\beta^{\prime}$}\\\ &\text{have no vertex
in common}\\\ -1&:\text{if both spins $\beta$, $\beta^{\prime}$ }\\\
&\text{leave or enter a common vertex $v$}\\\ 1&:\text{if one spin $\beta$
leaves a common vertex}\\\ &\text{$v$ and $\beta^{\prime}$ enters it, and
viceversa}.\end{cases}$ (6)

Given the matrix $Q$, we can write $Q=2I-A$ and define $A$, which is a
directed incidence matrix with support on the line graph
$\mathcal{L}(\mathcal{G})$. In order to gain some intuition about the matrix
$A$, let us discuss the non-directed case first.

###  II.1 Undirected Line Graphs

We start by defining line graphs, which are graph constructed from an
undirected graph $\mathcal{G}$ and such that the edges of the graph
$\mathcal{G}$ become the vertices of the graph $\mathcal{L}(\mathcal{G})$. The
edges (or edges) of the line graph are constructed based on the connectivity
of the original graph, as follows Whitney ; Krausz ; Harary ; Beineke .111It
is interesting to note that here there is a mismatch between the original
literature in graph theory, starting with the original work of Harary Harary
(1965). The original Line Graph of a digraph did not have any negative values,
e.g. if the two edges do not have zero sum, then we assign a value zero.

Let $\mathcal{G}=(V,E)$ denote a graph with vertex set
$V=\\{v_{1},v_{2},...v_{n}\\}$ and edge set $E=\\{e_{1},e_{2},...,e_{p}\\}$

Figure 1: The Line Graph construction. Black vertices and dashed lines
correspond the original graph $\mathcal{G}$, while grey vertices and the solid
lines correspond to $\mathcal{L}(\mathcal{G})$.

Each vertex $\tilde{v}\in{\widetilde{V}}({\mathcal{L}}(\mathcal{G}))$
corresponds to an edge $e\in E(\mathcal{G})$. Two vertices $\tilde{v}_{1}$ and
$\tilde{v}_{2}$ in $\widetilde{V}(\mathcal{L}(\mathcal{G}))$ are adjacent if
and only if the edges in $\mathcal{G}$ (corresponding to $\tilde{e}_{1}$ and
$\tilde{e}_{2}$) share a vertex. The correspondence between $\mathcal{G}$ and
$\mathcal{L}(\mathcal{G})$ is injective but not surjective. From a given graph
$\mathcal{G}$ we can construct only one $\mathcal{L}(\mathcal{G})$, and an
example is provided in Fig. 1. Yet, in general it is not true that any graph
can be thought as the line graph or another graph. In fact, according to the
Beineke classification, there are 9 non-minimal graphs that are not line
graphs of another graph, and each graph containing them is thus not a line
graph descendent of any other graph. We will discuss this later in detail
Beineke .

Given a graph $\mathcal{G}$, we can construct its line graph using the
following procedure:

  1. 1.

Enumerate the vertices of $\mathcal{G}$.

  2. 2.

Enumerate the edges of $\mathcal{G}$ with a fixed prescription (see example
below).

  3. 3.

If two edges share a vertex, draw a bold line between them.

  4. 4.

Remove $\mathcal{G}$ and its enumeration.

What is left is the line graph of $\mathcal{G}$, $\mathcal{L}(\mathcal{G})$.

Consider now the Kirchhoff matrix, obtained from the (undirected) incidence
matrix of the graph $\mathcal{L}(\mathcal{G}).$

The Kirchhoff matrix $K$ is the $p\times p$ matrix built from $P$, such that:

$K=B^{t}B,$ (7)

$B^{t}$ being the transpose of $B$. A well-known theorem now gives the
relationship between the incidence matrix and the adjacency matrix of the line
graph $\mathcal{L}(\mathcal{G})$:

Let $\mathcal{G}$ be a graph with $p$ edges and $n$ vertices and let
$\mathcal{L}(\mathcal{G})$ be its line graph. Then we have

$K=A-2\ I,$ (8)

where $I$ is the $p\times p$ identity matrix, and $A$ is the adjacency matrix
of $\mathcal{L}(\mathcal{G})$.

###  II.2 Directed Line Graphs

We see immediately that the definition of $Q$ and $K$ are very similar, with
an important difference. The matrix $Q$ can be written as

$Q=2I-A,$ (9)

where $A$ is called the weighted adjacency matrix, has the same support as the
undirected line graph $\mathcal{L}(\mathcal{G})$, but can take both positive
and negative values on edges of the line graph, depending on the orientation
$\mathcal{O}$ (see Fig. 1 and Fig. 2). In fact, $A$ takes a positive value
$+1$ if the two edges have zero sum on the vertex according to the
orientation, e.g. if one leaves and one enters, while $+1$ if they both enter
and leave.

A fundamental result follows: in general, we can write any spin ice model
(written in charge formulation), as

$H=-J\sum_{\beta\beta^{\prime}}A_{\beta\beta^{\prime}}s_{\beta}s_{\beta^{\prime}}$
(10)

where $A$ is the weighted adjacency matrix according to the rule above. Thus,
in the case of directed graphs we can have both ferromagnetic
$(JA_{\beta\beta^{\prime}}>0)$ and antiferromagnetic values
$(JA_{\beta\beta^{\prime}}<0)$, and is thus a weighted adjacency matrix. If
the element of $A$ is positive, the interaction on the line graph is
ferromagnetic (e.g. the two spins are aligned in the ground state), while if
it is negative the interaction is antiferromagnetic (anti-aligned in the
ground state).

Note that spin ices are generally frustrated, but frustration cannot be
reabsorbed by a spin redefinition, as it invariant under the Ising model gauge
freedom Hey  which in our graph-theoretical language corresponds to
$s_{\beta}\rightarrow\xi_{\beta}s_{\beta}$, $B_{v\beta}\to
B_{v\beta}\xi_{\beta}$ for $\xi_{\beta}=\pm 1$. As such, the couplings that
one obtains in the procedure depend on the gauge transformation. What does not
change is the frustration, which cannot be removed.

As an example, consider the Hexagonal spin ice of Fig. 2. We see that with the
orientation we have used, the only frustrated cycles are those associated to
the vertices of the Hexagonal model. This implies naturally that the
degeneracy of the ground state in the model must scale with the size of the
vertices.

Figure 2: The weighted line graph $A$ superimposed to the hexagonal lattice,
for a given orientation, where blue are antiferromagnetic couplings (negative)
and red are ferromagnetic (positive). Figure 3: Equivalent interaction model
on the directed line graph, with red lines equal to antiferromagnetic
interactions. Every single fundamental cycle (the triangles) are frustrated.

###  II.3 Spin Ice, Frustration, and Degeneracy: General Facts

We are now in a position to state some general facts for a spin ice on a
graph.

Remark 1 Not all frustrated Ising models are spin ices. Because the directed
line graph dual has support on the undirected line graph, we can borrow the
results from the undirected case Beineke . If the line graph contains any of
the subgraphs contained in the Beineke classification, then we know that the
graph is not the line graph of any root graph.

Remark 2 Vertices are mapped to complete graph interactions. This is a well
known fact that we restate graph theoretically. Vertices in a spin ice are
mapped to a complete graph with a number of vertices equal to the degree of
the vertex. This implies immediately that if a spin ice is composed by a
sequence of vertices of degree $d=\\{d_{1},\cdots,d_{n}\\}$, if for any $i$ we
have $d_{i}>3$, the line graph dual will not be planar.

Remark 3 The only purely ferromagnetic spin ice graph is the one dimensional
Ising model. It is also the only spin ice whose ground state is non-
degenerate. Consider the following one-dimensional spin ice:

$H=J\sum_{i=1}^{n}(s_{i}-s_{i+1})^{2}.$ (11)

The ground state has a $Z_{2}$ symmetry: these are all right or all left
spins, which is equivalent to a 1-dimensional ferromagnetic Ising model. This
implies that at least one spin ice is non-degenerate. Interestingly, it is the
only one.

In order to see this, consider a vertex with $d$ edges or spins. Then, the
total number of interaction terms are $d(d-1)/2$. Assume an orientation in
which $d_{1}$ spins go in and $d_{2}$ go out, with $d=d_{1}+d_{2}$. Then, the
number of ferromagnetic and antiferromagnetic interactions are

antiferromagnetic $\displaystyle:$
$\displaystyle\frac{d_{1}(d_{1}-1)}{2}+\frac{d_{2}(d_{2}-1)}{2},$
ferromagnetic $\displaystyle:$ $\displaystyle d_{1}d_{2}.$ (12)

The only case in which we have no antiferromagnetic interaction is
$d_{1}=d_{2}=1$, which is a vertex of degree $2$. It follows that the only
graphs that can be formed with vertices of the degree $2$ are either line or
circle graphs, and thus one dimensional spin ices are the only ferromagnetic
models.

This does not mean that all models are necessarily extensively degenerate,
that is possess a nonzero Pauling entropy bernal ; pauling . This is a more
complicated notion which depends on the product of signs of interaction on a
loop. We discuss this next.

Remark 4: All spin ices with $d>2$ have frustration at the vertex level.

Since frustration is gauge invariant, we can pick any orientation of the
vertex configuration and calculate frustration along a certain loop (a closed
sequence of edges, or a loop) in that particular configuration. Let us choose
$d_{2}=0$, thus all spins going into the vertex. Now, all the fundamental
cycles at the vertex interactions are of length 3, as the effective
interaction is a complete graph $K_{d}$. There are $m=\frac{d(d-1)(d-2)}{3!}$
fundamental circuits of length $3$. Since all interactions are anti-
ferromagnetic, we have that the product of the signs in every cycle is $-1$,
and thus frustrated.

Remark 5: For planar spin ices, frustration is only at the vertex level.

This is a byproduct of the following fact. Consider a cycle in a spin ice (in
the original lattice). We can always choose an orientation of the lattice such
that, for a given cycle, the arrows are chosen head to tail. Thus, when we
construct the directed line graphs, the interactions bordering two spins are
going to be ferromagnetic. Thus, if we go around a cycle in the line graph
dual, we only have ferromagnetic interactions, and thus the cycle is not
frustrated. In order to see that this is true always, note that for a planar
graph we can always choose orientations of the spin such that such
configuration is consistent. Thus, frustration is due only to vertices
(cliques) in the directed line graph dual.

Note that this implies that so-called vertex-frustration Morrison , that is
the inability to arrange collectively all vertices in a lowest energy
configuration, cannot exist in a graph spin ice whose Hamiltonian depends only
on the vertex charge. Indeed, all the vertex-frustrated systems Morrison ,
many of which have been realized tetris ; shakti  and depend upon a lifting of
degeneracy within vertices of the same charge. They are therefore not pure
spin ice graphs.

###  II.4 Spin ice reconstruction via Krausz clique partitions

One of the most interesting byproducts of the direct construction is that
there is an inverse procedure, known as Krausz partitioning. We know that if
the original spin ice interaction is planar, then vertices are mapped to fully
frustrated cliques (condition 1). Also, any cycle subgraph which is not a
clique must not be frustrated (condition 2). If these conditions are
satisfied, and if none of the Beineke graphs are present, then we can
reconstruct the original spin ice interaction via the Krausz decomposition,
which goes as follows (condition 3).

Figure 4: The coarse graining procedure according to the Krausz partitioning.

Given a graph $\mathcal{Q}$, if condition $1$ and $2$ are satisfied, consider
the unweighted graph $|\mathcal{Q}|$, and

  1. 1.

Enumerate all complete subgraph of the graph $|\mathcal{Q}|$, and defined as
partitions $\mathcal{K}$;

  2. 2.

If all partitions $\mathcal{K}$ have only one vertex in common, contract the
cliques into a vertex, and connect the partitions by an edge

  3. 3.

The resulting graph is the spin ice interaction matrix: assign spins to the
edges of the resulting graph and add a term $J(\sum_{i}s_{i})^{2}$ to the
corresponding interaction.

  4. 4.

Because of gauge invariance, the directionality of the original spin is
irrelevant.

It is important to note that a graph is a line graph of a root graph
$\mathcal{G}$ if and only if there is a Krausz partitioning of the graph
$\mathcal{Q}$. An example of such procedure is in Fig. 4. Each complete
subgraph is identified and the coarse graining procedure is performed.

##  III Lower and upper bounds to degeneracy

Extending Pauling’s estimate pauling  can become extremely challenging on
arbitrary graphs. In order to estimate the degeneracy of the ice manifold
beyond the case of planar graphs or non-bipartite graphs, we will use a
general graph-theoretic approach, and upper bound the entropy associated to
the ice manifold. First, let us note that the maximum degeneracy that a spin
ice can have is given by $2^{N_{spins}}$. If the spin ice is degree regular
and has $N_{v}$ vertices, then the maximum entropy of the ice manifold is
naturally given by $\epsilon_{max}=N_{v}\frac{d_{v}}{2}\ln 2.$ Below we
provide a procedure to systematically calculate upper bounds based on the
theory of Eulerian tours on graphs. Given a certain spin ice graph
$\mathcal{G}$, we are interested in calculating the number of configurations
in the ice manifold, $\epsilon(\mathcal{G})$, and its entropy
$S(\mathcal{G})=\ln\epsilon(\mathcal{G})$ pauling .

###  III.1 Case of even degrees

We first focus on the case in which all vertices have even degree,
independently from planarity. Consider a graph $\mathcal{G}$. It is very well
known that Euler got interested in the problem of walks on graphs with the
following property: starting from a certain vertex $v$, perform a walk on the
graph $\mathcal{G}$ such that you never use the same edge twice. A Euler cycle
starts at a vertex $v$ and ends in the vertex $v$. Euler proved the following
theorem Euler1 :

Theorem (Euler) Let $\mathcal{G}$ be a connected graph. Then, $\mathcal{G}$
has an Euler cycle if and only if every vertex $v$ has even degree $d_{v}$.

One direction of this theorem is rather obvious, as if a Euler cycle exists,
necessary $d_{v}$ must be even. The theorem is powerful because it ensures
that if all $d_{v}$ are even, the converse also applies. While an enumeration
of the number of eulerian tours for undirected graph is an open problem, a
formula for the number of eulerian orientations exists. Let $\mathcal{G}$ be a
graph with degrees $d_{v}$. Then, the total number of Eulerian orientations
$\epsilon(\mathcal{G})$ is given by schrijver :

$\displaystyle\epsilon(\mathcal{G})=\frac{\text{perm}(A)}{\prod_{v\in
V}(\frac{d_{v}}{2})!},$ (13)

where the matrix $A$ is constructed as follows. Let us consider the incidence
matrix $B$ of the undirected graph $\mathcal{G}$. For each vertex $v$, $A$
contains a $\frac{d_{v}}{2}$ number of identical rows of $B$. This implies
that $A$ is square and of size equal to the number of edges of $\mathcal{G}$.
The formula above is exact but hard to use, given the fact that the permanent
is rather hard to calculate, being $\\#P$-complete valiant . However, one can
upper bound the permanent of $(0,1)$ matrices using the Bregman-Minc result
for the permanent schrijver , which is

$\displaystyle\text{perm}(A)\leq\prod_{i=1}^{m}(r_{i}!)^{\frac{1}{r_{i}}},$
(14)

where $r_{i}$ is the row sum of the i-th row of $A$. This implies that, if
$d_{v}$ is the degree of the graph, one has Schrijver’s inequality schrijver

$\displaystyle\epsilon(\mathcal{G})\leq\prod_{v\in
V}\sqrt{\binom{d_{v}}{\frac{d_{v}}{2}}}.$ (15)

We note that the upper bound above is base on the fact that the graph has an
eulerian orientation, and thus $d_{v}$ must be even. For degree regular
graphs, we have

$\displaystyle\ln\epsilon(\mathcal{G})\leq\frac{N}{2}\ln\binom{d_{v}}{\frac{d_{v}}{2}}.$
(16)

which is a bound of fairly general nature and depends on the graph properties,
such as the vertex degree.

###  III.2 Approximating the number of Eulerian configurations from the
permanent

For completeness, we discuss first a technique which proved unfruitful for us,
but which deserves to be mentioned. One way to calculate such quantity
exploits the Godsil-Gutman theorem GGe ; Lovasz . Let $A$ be a non-negative
matrix. Then, if we define $B_{ij}=\epsilon_{ij}\sqrt{a_{ij}}$, where
$\epsilon_{ij}$’s are uncorrelated random variables distributed according to
$P(\epsilon_{ij})$, with mean $0$ and variance $1$, we have

$\displaystyle\text{perm}(A)=\langle\text{det}B^{t}B\rangle_{P(\epsilon)}.$
(17)

Since we are interested in the logarithm, we have

$\displaystyle\ln\text{perm}(A)=\ln\langle\text{det}(B)^{2}\rangle_{P(\epsilon)}=\ln\langle\text{det}(B^{t}B)\rangle_{P(\epsilon)}.$
(18)

If $\epsilon_{ij}=\\{1,-1\\}$ the estimator above is called Godsil-Gutman, but
if $\epsilon_{ij}\in\mathcal{N}(0,1)$ is is called Barvinok estimator. We have
tested both the Gutman-Godsil and Barvinok estimators for the case of the
triangular, square and cubic lattice degeneracies, but we have found that the
variance of the estimates does not fall fast enough with the number of Monte
Carlo samples.

Figure 5: Hypergraph construction for the spin ice degeneracy.

The second approximation method for the permanent of a non-negative matrix is
based on Belief Propagation, which is the one we present here huang . Consider
the square matrix $A$ obtained via the Schrijver augmentation. The permanent
of the matrix $A$ is defined via the sum over all possible permutations, as

$\text{perm}(A)=\sum_{\pi\in S_{n}}f(\pi;A),$ (19)

where $f(\pi;W)=\prod_{i=1}^{n}A_{i\pi(i)}$. Another way of thinking of these
permutations is in terms of perfect matching between two sets $A$ and $B$, and
in particular “double dimerizations”, as follows.

Given a certain permutation $\pi$, a matching between the set A (the first
index) and the set B (the second index), can be represented as
$\Big{(}i,\sigma(i)\Big{)}$. Similarly, a particular valid configuration of
the permanent is a set of $n$ non-overlapping dimers
$A_{1,\sigma(1)}....A_{n,\sigma(n)}$. So one constructs a bipartite graph in
which one places a dimer between $(i,j)$ if $A_{ij}>0$. The permanent is a sum
over all possible dimerizations.

Based on this idea, Huang and Jebara mapped the permanent to a set of double-
dimerizations which correspond to a valid choice of the permanent as follows
huang . A dimer between the set $A$ and $B$ is an assignment between the
variables $X=\\{x_{1},\cdots,x_{n}\\}$ and the variables
$Y=\\{y_{1},\cdots,y_{n}\\}$. We now introduce the potentials
$\phi(x_{i})=\sqrt{A_{ix_{i}}}$ and $\phi(y_{j})=\sqrt{A_{y_{j}j}}$ and
introduce the function

$\displaystyle
f(X,Y)=\prod_{ij}\psi(x_{i},y_{j})\prod_{k}\phi_{x_{k}}\phi(y_{k})$ (20)

which is a function of the assignment. If the function $\psi(x_{i},y_{j})$
enforces that, given two sets of assignments (two possible dimer
configurations) between $A$ and $B$, the two assignments are identical, then
ones obtains

$\displaystyle\prod_{k}\phi_{x_{k}}\phi(y_{k})=A_{1,\sigma(1)}....A_{n,\sigma(n)},$
(21)

and

$\displaystyle Z\equiv\text{\text{Bperm}}(A)=\sum_{\sigma,\pi\in
S_{n}}f(X,Y).$ (22)

where $\sigma=(i,x_{i})$ and $\pi=(y_{j},j)$. In terms of the dimer
representation, a valid configuration is such that two dimers either overlap
completely or do not overlap at all. Then, a logic function which ensures such
condition is the negation of the XOR function $I(\neg(j=x_{i}\oplus
i=y_{j}))$, where the function $I(\cdot)$ is zero if the condition is false
and one otherwise.

It is interesting at this point to note that eqn.(22) can be interpreted as
the partition function for a particular factor graph with pairwise
interactions, and can be analyzed in terms of the belief propagation. Let us
define the Bethe free energy $F_{\text{Bethe}}$ as

$\displaystyle F_{\text{Bethe}}$ $\displaystyle=$
$\displaystyle-\sum_{ij}\sum_{x_{i},y_{k}}b(x_{i},y_{j})\ln\Big{(}\psi(x_{i},y_{j})\phi(x_{i})\phi(y_{j})\Big{)}$
(23) $\displaystyle+$
$\displaystyle\sum_{ij}\sum_{x_{i},y_{j}}b(x_{i},y_{j})\ln b(x_{i},y_{j})$
$\displaystyle-$ $\displaystyle(n-1)\sum_{i}\sum_{x_{i}}b(x_{i})\ln b(x_{i})$
$\displaystyle-$ $\displaystyle(n-1)\sum_{j}\sum_{y_{j}}b(y_{j})\ln b(y_{j})$

Using Belief propagation, the partition function can be written as a
minimization of the Bethe free energy

$\displaystyle Z=e^{-\text{min}_{b}F_{\text{Bethe}}(b)}.$ (24)

The minima of the belief can be obtained via a message passing algorithm. The
beliefs must satisfy $b(x_{i})=\sum_{y_{j}}b(x_{i},y_{j})$ and
$b(y_{j})=\sum_{x_{i}}b(x_{i},y_{j})$, and
$\sum_{x_{i},y_{j}}b(x_{i},y_{j})=1$, and these functions can be obtained
iteratively as

$b(x_{i},y_{j})\propto\psi(x_{i},y_{j})\phi(x_{i})\phi(y_{j})\prod_{k\neq
j}m_{y_{k}}(x_{i})\prod_{l\neq i}m_{x_{l}}(y_{j})$ (25)

and,

$\displaystyle b(x_{i})$ $\displaystyle\propto$
$\displaystyle\phi(x_{i})\prod_{l\neq i}m_{y_{l}}(x_{i})$ $\displaystyle
b(y_{j})$ $\displaystyle\propto$ $\displaystyle\phi(y_{j})\prod_{l\neq
i}m_{x_{l}}(y_{j}).$ (26)

The messages can be obtained iteratively, starting from a random initial
state, one has

$\displaystyle
m_{x_{i}}^{t+1}(y_{j})=\sum_{x_{i}}\Big{(}\phi(x_{i})\psi(x_{i},y_{j})\prod_{k\neq
j}m_{y_{k}}^{t}(x_{i})\Big{)}$ (27)

An interesting byproduct is that it is possible to prove that the Bethe
permanent can serve as a lower and upper bound vontobel ; gurvits :

Theorem (Vontobel-Gurvits)

$\displaystyle\text{\text{Bperm}}(A)\leq\text{perm}(A)\leq\sqrt{2}^{n}\text{\text{Bperm}}(A).$
(28)

where $n$ is the size of $A$ and $\text{Bperm}(A)$ is its Bethe permanent.

The theorem above implies that we can obtain, via the Bethe permanent Bp$(A)$,
a level of confidence on the value of the permanent and in particular a
certificate for the lower bound scaling. Note that in practice, we have found
that Schrijver’s upper bound is typically lower than the one obtained via the
Bethe Permanent.

Given a certain lattice described by the graph $\mathcal{G}$ we call
$B\epsilon(\mathcal{G})$ the lower bound obtained via the permanent. We thus
have

$\displaystyle
B\epsilon(\mathcal{G})\leq\epsilon(\mathcal{G})\leq\prod_{v}\sqrt{\binom{d_{v}}{\frac{d_{v}}{2}}},$
(29)

where $d_{v}$’s are the vertex degrees.

The numerical results are shown in Fig. 6 for the planar cases of the square
lattice and the triangular tiling, which are two perfect Archimedean lattices
harrison . For the case of the square ice we have Lieb’s exact result lieb
$\epsilon^{Exact}=(\frac{4}{3})^{3L^{2}/2}$. For the case of the triangular
tiling there is no exact solution. Pauling’s argument pauling  does not apply,
as it relies on the bipartiteness of the lattice. For the case of the
triangular tiling, we find that $\frac{\ln\epsilon}{L^{2}}\geq 2.33[..]$,
using the fact that the Bethe permanent is a lower bound.

Figure 6: Numerical scaling of the degeneracy of the ground state via the
Bethe Permanent for the square spin ice and the triangular tiling (on the
torus), of linear size $L$ (equal to the number of nodes). The lower bound
(gold solid) corresponds to the Bethe Permanent $\text{Bperm}(A)$, while the
upper bound (purple dashed) corresponds to $\sqrt{2^{n}}\text{Bperm}(A)$. A
numerical fit shows that the Bethe permanent scales as
$\epsilon_{SI}\approx\text{Bperm}(A_{L})/(\prod(d_{v}/2)!)\approx(1.419[..])^{L^{2}}$,
while Lieb’s exact result is
$\epsilon_{SI}^{Exact}=(\frac{8}{3\sqrt{3}})^{L^{2}}\approx(1.53[..])^{L^{2}}$.
For the case of the triangular tiling, the scaling of the Bethe Permanent can
be fit as $\epsilon_{TT}\approx(2.3396[..])^{L^{2}}$. In both figures, the
shaded area is the bound on the (logarithm) of the degeneracy of the balanced
configuration according to the Bethe Permanent bounds.

For the cubic lattice (with toroidal boundary conditions), Pauling’s
calculation suggests that the entropy of the spin ice scales with the number
of vertices $L^{3}$. In fact, given a certain node, we have $2^{6}=64$
possible configurations, but only $20$ of them satisfy the ice rule. It
follows that, from Pauling’s argument, the spin ice ground state degeneracy
should be

$\displaystyle\epsilon_{\text{Pauling}}=2^{3L^{3}}(\frac{20}{64})^{L^{3}}=(5/2)^{L^{3}}.$
(30)

or $\ln\epsilon_{\text{Pauling}}=L^{3}\ln\frac{5}{2}\approx 0.916(3)\cdot
L^{3}$. Using the Bethe permanent (see Fig. 7), we observe that Pauling’s
estimate is not too far from the Bethe permanent lower bound, which gives
$B\epsilon\approx 2.41^{L^{3}}$, but provides a certificate for a lower bound.

Figure 7: Numerical scaling of the degeneracy of the ground state via the
Bethe Permanent for the cubic spin ice with a number of nodes $L^{3}$. The
lower bound (gold solid) corresponds to the Bethe Permanent $\text{Bperm}(A)$,
while the upper bound (purple dashed) corresponds to
$\sqrt{2^{n}}\text{Bperm}(A)$. A numerical fit shows that the Bethe permanent
scales as
$\epsilon_{SI}\approx\text{Bperm}(A_{L})/(\prod(d_{v}/2)!)\approx(2.41(0))^{L^{3}}$,
while Pauling’s estimate is $\epsilon_{\text{Pauling}}=(2.5)^{L^{3}}$. The
lower bound obtained via the Bethe permanent is thus within $3\%$ of Pauling’s
estimate.

As a last application, we consider regular tournaments, which is the total
number of Eulerian configurations for a complete graph with $L$ nodes, as in
Fig. 8, or equivalently the degeneracy of the spin ice configurations on a
complete graph. This number was calculated by McKay mckay  and is given by

$\displaystyle\epsilon=\left(\frac{2^{L+1}}{\pi
L}\right)^{\frac{L-1}{2}}\frac{\sqrt{L}}{\sqrt{e}}\left[1+O(\frac{1}{\sqrt{L}})\right].$
(31)

It can be easily seen that also in this case the Bethe permanent provides a
good estimate for the number of tournaments.

Figure 8: Number of Eulerian orientations for the complete graph $K_{L}$, also
called regular tournaments. We provide a comparison between the upper and
lower bounds given by the Bethe Permanent and the McKay exact calculation
given in eqn. (31). We note that for this particular case, Schrijver upper
bound is looser than the Bethe permanent upper bound.

A summary of the results is provided in the Table 1.

##  IV Conclusions

We have discussed some theoretical results for spin ice on arbitrary graphs.
In the first part of the paper we have provided a graph theoretical mapping
between spin ices and Ising models based on the directed incidence matrix of
the Line graph. Specifically, we have shown something that was already known:
while all spin ices can be mapped to Ising models, not all Ising models can be
mapped to spin ices; but here we show a general a procedure to extract the
original spin ice. We have also shown that all spin ices are degenerate
(except the trivial case: the 1D Ising model). We have proved this result
using the known gauge transformation for the Ising model. Another fact which
we managed to prove in the first part of this paper, and that we found
surprising, is that for planar spin ices, degeneracy is only at the vertex
level. This implies that the effective Ising model, cycles between different
vertices of the original spin ice are unfrustrated, as it can be seen via a
gauge transformation to align these vertices. The same is not true for the
interactions at the vertices, implying that the frustration scales with the
number of vertices of the original spin ice. In the second part of the paper
we have focused on techniques to bound the degeneracy of the ice manifold. The
method is based on the number of Eulerian tours, as for every balanced graph
(ice manifold state) there exists at least one Eulerian trail. The advantage
of using the Eulerian trails is that there are spectral techniques to count
the number of Eulerian trails of a directed graph. We applied these techniques
for even-degree and connected graphs, showing that these bounds can be
extended to the case of odd-degree reducible graphs. To conclude, we have used
an exact formulation for the degeneracy of spin ices configurations using a
permanent identity. Given the fact that the permanent is $\\#P$-complete
quantity to compute, we have employed numerical methods based on the Bethe
free energy. Such method provides numerical lower bounds to the permanent, and
we have thus obtained lower bounds to the spin ice degeneracy for various spin
ice regular lattices. The advantage of such procedure is that these certified
lower bound estimates can be obtained relatively quickly and without using
loop Monte Carlo techniques barkema ; loop ; giawei  which, however, would
give a more precise estimate.

In this respect, it is worth mentioning that there are also other
implementations using fractional belief propagation chertkov , with an extra
parameter which is known (given the matrix) to provide the exact result of the
permanent, and for values above and below give upper and lower bounds
respectively. These will be considered in the future.

Graph | Exact | Pauling | Bethe | Maximum
---|---|---|---|---
Square | 1.53[..] | $\frac{3}{2}$ | 1.41 | $4$
Triangular |  |  | 2.33 | 8
Cubic |  | $\frac{5}{2}$ | 2.41 | 8
Table 1: Normalized ground state entropy $s=\frac{\ln\epsilon}{N_{v}}$. The
Bethe permanent typically provides a certified lower bound given the matrix,
and the scaling is extracted numerically. For the complete graph, tournaments
do not scale simply exponentially with the number of vertices (see eqn. (31)).
The maximum is obtained via the relationship for regular graphs, $s=\ln
2^{\frac{d}{2}}$, with $N_{v}$ the total number of vertices.

###### Acknowledgements.

We thank Prof. A. Schrijver for some clarifications regarding the permanent
formula for the eulerian orientations of eqn. (13), and Prof. M. Chertkov for
some clarifications regarding the Bethe Permanent. We also acknowledge the
support of NNSA for the U.S. DoE at LANL under Contract No. DE-AC52-06NA25396.
This work was carried out under the auspices of the U.S. DoE through the Los
Alamos National Laboratory, operated by Triad National Security, LLC (Contract
No. 892333218NCA000001). FC was also financed via DOE-LDRD grants PRD20190195.
M. Saccone acknowledges the Center for Nonlinear Studies for a fellowship.

## References

  * (1) J. Bernal and R. Fowler, A theory of water and ionic solution, with particular reference to hydrogen and hydroxyl ions,” The Journal of Chemical Physics 1, 515-548 (1933).
  * (2) L. Pauling, The structure and entropy of ice and of other crystals with some randomness of atomic arrangement,” Journal of the American Chemical Society 57, 2680-2684 (1935).
  * (3) Lieb E.H. (2004) Residual Entropy of Square Ice. In: Nachtergaele B., Solovej J.P., Yngvason J. (eds) Condensed Matter Physics and Exactly Soluble Models. Springer, Berlin, Heidelberg.
  * (4) A. P. Ramirez, A. Hayashi, R. J. Cava, R. Siddharthan, and B. S. Shastry, Zero-point entropy in ‘spin ice’,” Nature 399, 333-335 (1999).
  * (5) C. Nisoli et al.,“Colloquium: Artificial spin ice: Designing and imaging magnetic frustration”, Rev. Mod. Phys. 85(1473), (2013)
  * (6) R. F. Wang et al., Artificial spin ice sin a geometrically frustrated lattice of nanoscale ferromagnetic islands, Nature 439(7074):303-6, (2006).
  * (7) S. H. Skjærvø, C. H. Marrows, RL Stamps, LJ Heyderman “Advances in artificial spin ice” Nature Reviews Physics 1-16 (2019).
  * (8) S.D. Bader, Colloquium: Opportunities in nanomagnetism, Rev. Mod. Phys., 78(1):1, (2006).
  * (9) I. Gilbert et al., “Emergent reduced dimensionality by vertex frustration in artificial spin ice”, Nature Phys. 12, 162-165 (2016)
  * (10) L. J. Heyderman, R. L. Stamps, Artificial ferroic systems: novel functionality from structure, interactions and dynamics, J. of Phys.: Condensed Matter, 25(36):363201 (2013)
  * (11) B. Canals et al., Fragmentation of magnetism in artificial kagome dipolar spin ice, Nat. Comm. 7 (2016)
  * (12) C. Nisoli et al., Ground State Lost but Degeneracy Found: The Effective Thermodynamics of Artificial Spin Ice, Phys. Rev. Lett., 98(21):217203 (2007)
  * (13) J. P. Morgan et al., Thermal ground state ordering and elementary excitations in artificial magnetic square ice, Nat. Phys. 7(1):75-70 (2010)
  * (14) Z. Budrikis et al., Disorder strength and field-driven ground state domain formation in artificial spin ice: experiment, simulation and theory, Phys. Rev. Lett 109(30):037203 (2012)
  * (15) W. R. Branford et al., Emerging Chirality in Artificial Spin Ice, Science, 335(6076):1597-1600 (2012)
  * (16) I.A. Ryzhkin. Zhurnal Ehksperimentalnoj i Teoreticheskoj Fiziki, 128(3):559-566 (2005)
  * (17) G. Moeller, R. Moessner, Magnetic multipole analysis of kagome and artificial spin-ice dipolar arrays, Phys. Rev. B, 80(14):140409 (2009)
  * (18) G.-W. Chern, P. Mellado, Magnetic monopole polarons in artificial spin ices, EPL 114 (3): 37004 (2016)
  * (19) B. L. Le et al., Understanding magnetotransport signatures in networks of connected permalloy nanowires. Phys. Rev. B, 95:060405 (2017)
  * (20) G.-W. Chern, Magnetotransport in Artificial Kagome Spin Ice, Phys. Rev. App. 8(6) : 064006 (2017)
  * (21) S. Gliga, et al., Spectral analysis of topological defects in an artificial spin-ice lattice, Phys. Rev. Lett, 110(11):117205 (2013).
  * (22) I. Gilbert et al., Emergent ice rule and magnetic charge screening from vertex frustration in artificial spin ice, Nat Phys. 10(9):670-675 (2014)
  * (23) V. S. Bhat et al., Controlled magnetic reversal in permalloy films patterned into artificial quasicrystals, Phys. Rev. Lett. 111(7):077201 (2013)
  * (24) F. Caravelli, A model for the Mediated Artificial Square Ice phenomenology, EPL 120(4),2020
  * (25) C. Nisoli, V. Kapaklis, P. Schiffer, Deliberate exotic magnetism via frustration and topology, Nature Phys.13(3):200-203 (2017)
  * (26) M J Morrison, et al., Unhappy vertices in artificial spin ice: new degeneracies from vertex frustration. New Journal of Physics, 15(4):045009 (2013)
  * (27) C. Castelnovo et al., Spin ice, fractionalization, and topological order, Ann. Rev. Condens. Matter Phys., 3(1): 35-55 (2012)
  * (28) Y. Lao et al., “Classical topological order in the kinetics of artificial spin ice”, Nature Phys. 14, 723-727 (2018)
  * (29) I. Gilbert et al., Direct visualization of memory effects in artificial spin ice. Phys. Rev. B, 92(10):104417 (2015)
  * (30) P. E. Lammert et al., Direct entropy determination and application to artificial spin ice. Nat. Phys., 6(10):786-789 (2010)
  * (31) D. Levis et al., Thermal phase transitions in artificial spin ice, Phys. Rev. Lett., 110(20):207206 (2013)
  * (32) F. Caravelli, G.-W. Chern, C. Nisoli, Artificial Spin Ice Memory Resistors, arXiv:1908.08073
  * (33) F. Caravelli, J. Carbajal, Memristors for the curious outsiders, Technologies 6(4):118 (2019)
  * (34) H. Arava et al., ”Computational logic with square rings of nanomagnets.” Nanotechnology 29, no. 26 265205 (2018)
  * (35) J. H. Hensen, E. Folven, G. Tufte, Computation in artificial spin ice, Proc. of ALIFE 2018, pp. 15-22, MIT Press, 10.1162/isal-a-00011 (2018)
  * (36) M. T. Niemier, et al.,‘Nanomagnet logic: progress toward system-level integration”, J. of Phys.: Condensed Matter, 23(49):493202 (2011)
  * (37) F. Caravelli, C. Nisoli, Logical gates embedding in Artificial Spin Ice, arXiv:1810.09190 (to appear in the NJP)
  * (38) J. C. Gartside et al., “Magnetic topological lithography: Gateway to the artificial spin ice manifold”, Nature Nano., 13(1):53-58 (2018)
  * (39) Y.-L. Wang et al., Rewritable artificial magnetic charge ice, Science 352, 6288: 962-966 (2016)
  * (40) Y.-L. Wang et al., Switchable geometric frustration in an artificial-spin-ice-superconductor heterosystem, Nature Nano. 13(7): 560 (2018)
  * (41) P. Vavassori, ”Towards plasmon-assisted thermal excitations in artificial spin ices”, invited talk at DPG Fruehjahrstagung, 19-24/03/2017, Dresden, Germany.
  * (42) Z. Li et al., “Simultaneous Local Heating/Thermometry Based on Plasmonic Magnetochromic Nanoheaters”, Small 14, 1800868 (2018)
  * (43) V. Schanilec, Y. Perrin, S. Le Denmat, B. Canals, N. Rougemaille, Artificial vertex systems by design, arxiv:1902.00452
  * (44) C. Nisoli, Write it as you like it, Nature Nano. 13(1): 5. (2018)
  * (45) L. A. S. Mol, A. R. Pereira, W. A. Moura-Melo, Extending spin ice concepts to another geometry: The artificial triangular spin ice, Phys. Rev. B 85, 184410 (2012)
  * (46) A. May, M. Hunt, A. Van Den Berg, A. Hejazi, S. Ladak , Realisation of a frustrated 3D magnetic nanowire lattice, Comm. Phys. (13) (2019)
  * (47) A. May, M. Saccone, A. van den Berg, J. Askey, M. Hunt, S. Ladak, Magnetic Charge Propagation upon a 3D Artificial Spin-ice. arXiv:2007.07618
  * (48) AD King, C Nisoli, ED Dahl, G Poulin-Lamarre, A Lopez-Bezanilla arXiv:2007.10555 (2020)
  * (49) C. Nisoli AIP Advances 10 (11), 115102 (2020).
  * (50) J. Bang-Jensen, G. Gutin, Digraphs:Theory, Algorithms and Applications, Springer-Verlag Berlin (2007)
  * (51) L. Euler, ”Commentarii academiae scientiarum petropolitanae,” Solutio problematis ad geometriam situs pertinentis 8, 128-140 (1736).
  * (52) Godsil, C. and Royle, G. “Algebraic Graph Theory”, New York: Springer-Verlag, 2001.
  * (53) R. Baxter, Exactly Solved Models in Statistical Mechanics, Academic Press (London), 1989
  * (54) M. Fisher, Transformations of Ising Models, Phys. Rev. 113:4, pp 969-981 (1959)
  * (55) F. Caravelli, F. Markopoulou, “Properties of Quantum Graphity at Low Temperature”, Phys. Rev. D 84(2), 2010
  * (56) H. Whitney, ”Congruent graphs and the connectivity of graphs”, American Journal of Mathematics, 54 (1): 150-168 (1932);
  * (57) J. Krausz, ”Démonstration nouvelle d’un théorème de Whitney sur les réseaux”, Mat. Fiz. Lapok, 50: 75–85 (1943); Rendiconti del Circolo Matematico di Palermo, 9 (2): 161–169 (1960);
  * (58) F. Harary, R. Z Norman, ”Some properties of line digraphs”,
  * (59) L. W. Beineke, ”Derived graphs of digraphs”, in Sachs, H.; Voss, H.-J.; Walter, H.-J. (eds.), Beiträge zur Graphentheorie, Leipzig: Teubner, pp. 17-33 (1968);
  * (60) A. J G Hey, Lattice gauge theory - an introductory review, Surveys in High Energy Physics, 5:4, 287-327 (1987)
  * (61) I. Gilbert, Y. Lao, I. Carrasquillo, L. O’Brien, J. D. Watts, M. Manno, C. Leighton, A. Scholl, C. Nisoli, P. Schiffer, Emergent reduced dimensionality by vertex frustration in artificial spin ice, Nature Physics 12, p; 162–165(2016)
  * (62) I. Gilbert, G.-W. Chern, S. Zhang, L. O’Brien, B. Fore, C. Nisoli, P. Schiffer, Emergent ice rule and magnetic charge screening from vertex frustration in artificial spin ice, Nature Physics 10, p 670-675 (2014)
  * (63) K. G. Valiant (1979), The Complexity of Computing the Permanent, Theor. Comp. Sci. 8 (2): 189–201.
  * (64) A. Schrijver, Bounds on the number of Eulerian Orientations, Combinatorica 3(3-4) 375-380, 1983.
  * (65) C. Godsil, I. Gutman, On the matching polynomial of a graph, in Algebraic Methods in Graph Theory, Vol. I,II (Szeged, 1978), North-Holland, Amsterdam–New York, 1981, pp. 241-249.
  * (66) L. Lovasz and M. D. Plummer, Matching Theory, North Holland, Amsterdam (1986).
  * (67) B. Huang, T. Jebara, Approximating the Permanent with Belief Propagation, arXiv:0908.1769
  * (68) P.O. Vontobel, The Bethe permanent of a non-negative matrix, in Proc. of Communication, Control, and Computing (Allerton), 2010 48th Annual Allerton Conference on, 29 2010-Oct. 1 2010; available at Pascal Vontobel home page.
  * (69) L. Gurvits. Unleashing the power of Schrijver’s permanental inequality with the help of the Bethe approximation, arXiv:1106.2844 (2011)
  * (70) A.Harrison, First catch your hare: the design and synthesis of frustrated magnets, Journal of Physics: Condensed Matter, 16(11), S553–S572 (2004)
  * (71) B. D. McKay, Combinatorica 10, pages 367–377 (1990)
  * (72) G. T. Barkema, M. E. J. Newman, Monte Carlo simulation of ice models, Phys. Rev. E 57, 1155–1166 (1998)
  * (73) R. G. Melko, B. den Hertog, M. J. P. Gingras, Long-Range Order at Low Temperatures in Dipolar Spin Ice, Phys. Rev. Let.. 87(6):067203 (2001)
  * (74) G.-W. Chern, O. Tchernyshyov, Magnetic charge and ordering in kagome spin ice, Phil. Trans. R. Soc. A (2012) 370, 5718–5737
  * (75) A. May et al.,Realisation of a frustrated 3D magnetic nanowire lattice., Comm. Physics , 2:13 (2019)
  * (76) M. Chertkov and A. B. Yedidia, “Approximating the permanent with fractional belief propagation”, In: The Journal of Machine Learning Research 14.1 (2013),