# Ratio of Symmetries Relating Graphs Differing by an Edge Set

Justus Isaiah Hibshman<EMAIL_ADDRESS>University of Notre Dame

###### Abstract

Given an undirected graph $G=(V,E)$ and a subset of its edges
$E^{\prime}\subseteq E$, one can consider the subgraph
$G^{\prime}=(V,E\setminus E^{\prime})$. We relate the total number of
automorphisms of $G$ and of $G^{\prime}$ via the automorphism orbit of
$E^{\prime}$. That is, we consider the automorphism orbit of $E^{\prime}$ in
$G$ (i.e. $\text{AO}_{G}(E^{\prime})$) as an orbit of a set of edges and the
automorphism orbit of $E^{\prime}$ in $G^{\prime}$ (i.e.
$\text{AO}_{G^{\prime}}(E^{\prime})$) as an orbit of a set of non-edges. We
then prove that the following ratio holds for any $G$ and any $E^{\prime}$:
$\frac{|\text{Aut}(G)|}{|\text{AO}_{G}(E^{\prime})|}=\frac{|\text{Aut}(G^{\prime})|}{|\text{AO}_{G^{\prime}}(E^{\prime})|}$.
The proof is relatively brief and makes use of the Erdős-Rényi graph
generation model.

##  I Introduction

The automorphism groups of a graph and of the same graph with some edges
removed can be arbitrarily different. For instance, Hartke, Kolb, Nishikawa,
and Stolee show this when they proved that given any two finite groups $A$ and
$B$ (which can be arbitrarily different), they construct a graph $G$ such that
$G$’s automorphism group is isomorphic to $A$ and if a particular edge $e$ in
$G$ is removed then the resulting graph’s automorphism group is isomorphic to
$B$ hartke2009automorphism . Since this holds for _all_ pairs of groups $A$
and $B$, it means that merely knowing a graph’s automorphism group (the “$A$”)
tells you nothing about what the automorphism group will be if an edge is
deleted (the “$B$”).

Until now the relationship between the automorphism group of an arbitrary
graph and the automorphism group of the same graph with some edges removed has
been relatively unknown. We provide such a link in this paper: The size of the
automorphism orbit of the deleted edge set. Our result applies to all finite
undirected graphs and the removal of any set of edges.

In addition to the fact that we find this result fascinating in its own right,
we believe it may have implications for the longstanding open problem of graph
reconstruction (and the simpler but still open problem of edge reconstruction)
bondy1977graph ; greenwell1969reconstructing ; bollobas1990almost ;
kocay1987family ; mcmullen2006graph ; andersen1996set . The graph
reconstruction conjecture was originally formulated in 1941 by Kelly and Unam,
and postulates that for any graph $G=(V,E)$, $G$ is uniquely identified by the
multiset of $|V|$ subgraphs obtained by deleting a single vertex from $G$
bondy1977graph . The multiset of subgraphs is called the “deck” (as in a deck
of cards). The conjecture has been described as a “disease” due to the way it
infects mathematicians with interest in the problem harary1969four .

In the reconstruction conjecture, removing a vertex involves removing both the
vertex itself _and_ the edges that were incident on the vertex. However, if
one constructs an augmented deck where one removes the edges but not the
vertex itself, one effectively has the same problem, and our theoretical
result directly applies. We spell out this connection a bit more fully in our
conclusion (Section V).

Beyond the graph reconstruction conjecture itself, there has been a fairly
large amount of work on vertex-deleted subgraphs lauri1992vertex ;
kirkland2010algebraic ; schmeichelm1979connectivity ; bowler2011recognizing ;
farrugia2021graphs , and less on edge-deleted (or edge-set-deleted) subgraphs,
though there is some andersen1992edge . We expect this paper’s result will
have relevance to much of this work.

We provide preliminary material in Section II, our main result in Section III,
and our proof in Section IV. Wrap-up discussion is in Section V.

##  II Preliminaries

###  II.1 Shorthand Notation

Given a graph $G$ with a set of vertices $V$, it is common to refer to the set
of edges in $G$ as being a subset $E\subseteq V\times V$. This is not quite
correct when the graph is undirected, but we will use the notation anyways, in
which case for undirected graphs a pair of nodes $(a,b)\in V\times V$ is
considered to be equivalent to its reverse (i.e. $(a,b)\equiv(b,a)$).

Given a set $S$ and a bijection $f:S\rightarrow S$, and given two elements
$a,b\in S$, we use $f((a,b))$ to denote $(f(a),f(b))$. Likewise, given a set
$S^{\prime}\subseteq S$ we use $f(S^{\prime})$ to denote $\\{f(x)\ |\ x\in
S\\}$. Lastly, given a set of pairs $X\subseteq S\times S$ we use $f(X)$ to
denote $\\{f((c,d))\ |\ (c,d)\in X\\}$.

Give an graph $G=(V,E)$ and an edge $e\in E$, we use “$G-e$” to denote
$(V,E\setminus\\{e\\})$. Likewise, given $E^{\prime}\subseteq E$, we use
“$G-E^{\prime}$” to denote $(V,E\setminus E^{\prime})$.

###  II.2 Iso- and Auto-morphisms

Given an undirected graph $G=(V,E)$ and a graph
$G^{\prime}=(V^{\prime},E^{\prime})$, an isomorphism between $G$ and
$G^{\prime}$ is a bijection $f:V\rightarrow V^{\prime}$ such that $(a,b)\in
E\leftrightarrow f((a,b))\in E^{\prime}$.

An automorphism of $G$ is simply an isomorphism of $G$ with itself.

We use Aut($G$) to denote the set of all automorphisms of $G$.

We use $[G]$ to denote the set of all graphs with vertex set $V$ which are
isomorphic to $G=(V,E)$. Formally,

$[G]=\\{(V,g(E))\ |\ g\in\text{Aut}(G)\\}$

We can indicate that graphs $G$ and $G^{\prime}$ are isomorphic by writing
“$G\cong G^{\prime}$” or by writing $G\in[G^{\prime}]$. We use “$\cong$” in
our introduction and conclusion due to it being standard notation, and we use
“$\in[\cdot]$” in our proof because it invokes the idea of a collection of
different edge sets all of which are isomorphic to each other - an idea that
our proof makes heavy use of.

###  II.3 Automorphism Orbits

Given a graph $G=(V,E)$ and some “entity” $x$ where $x$ could be a node ($x\in
V$), a node pair ($x\in V\times V$), or a _set_ of node pairs ($x\subseteq
V\times V$), the automorphism orbit of $x$ is the set of all entities
structurally equivalent to $x$:

$\text{AO}_{G}(x)=\\{f(x)\ |\ f\in\text{Aut}(G)\\}$

Example Graph $G$Automorphism Orbit ofEdge Set $\\{(a,e),(e,f)\\}$ in
$G$AO${}_{G}\left(\\{(a,e),(e,f)\\}\right)=\\{{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\textbf{\\{(a,
e), (e, f)\\}}},\ \
{\color[rgb]{0,0.6,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.6,0}\textbf{\\{(f,
g), (g,
c)\\}}},$${\color[rgb]{0.95,0.1,0.05}\definecolor[named]{pgfstrokecolor}{rgb}{0.95,0.1,0.05}\textbf{\\{(b,
e), (e, f)\\}}},\ \
{\color[rgb]{0.4,0.4,0.4}\definecolor[named]{pgfstrokecolor}{rgb}{0.4,0.4,0.4}\pgfsys@color@gray@stroke{0.4}\pgfsys@color@gray@fill{0.4}\textbf{\\{(f,
g), (g, d)\\}}}\\}$fgeabcdfgeabcd Figure 1: Example of an Automorphism Orbit
of a Set of Edges: On the right, each distinct set of edges in the
automorphism orbit is shown in a unique color. The multi-edges on the right
are shown solely to indicate that the single edges on the left participate in
multiple edge sets in the orbit. Note that sets such as $\\{(a,e),(f,g)\\}$
are not part of the example orbit.

Note that if $x=(a,b)$, $x$ could be an edge in $E$ or $x$ could just as
easily be a non-edge. Similarly, if $x$ is a set of node pairs, $x$ could be a
set of edges, a set of non-edges, or a mixture. Lastly, note that the
automorphism orbit of an entity $x$ contains the same kind of entities that
$x$ is, per our notation in Section II.1; the orbit of a node is a set of
nodes, the orbit of an edge is a set of edges, etc. To see an example of an
automorphism orbit of a _set_ of edges, see Figure 1.

###  II.4 Erdős-Rényi Graph Generation

The Erdős-Rényi graph generation model erdHos1960evolution  takes two
parameters $n$ and $m$ where $n\geq 1$ and $0\leq m\leq\binom{n}{2}$. It
generates a graph with $n$ nodes (labeled 1 through $n$) and $m$ vertices by
randomly sampling $m$ edges, each with equal probability. That is, we sample
$m$ pairs uniformly without replacement from $(V\times V)\setminus\\{(a,a)\ |\
a\in V\\}$ where $V=\\{1,2,...,n\\}$.

Given these generation parameters $n$ and $m$, the probability of generating a
graph $H$ isomorphic to $G=(V,E)$ where $|V|=n$ and $|E|=m$ is found via the
equation:

$\mathbb{P}(H\in[G])=\frac{1}{\binom{\binom{n}{2}}{m}}\cdot\frac{n!}{|\text{Aut}(G)|}$
(1)

The first fraction is the probability of generating a specific set of $m$
edges on $n$ nodes. That is because there are $\binom{n}{2}$ possible pairs of
nodes to sample as edges and $m$ that are actually selected. So there are
$\binom{\binom{n}{2}}{m}$ distinct sets of edges and thus the probability of
getting a particular one of those is just 1 divided by this quantity.

The second fraction is the number of distinct sets of edges that form a graph
isomorphic to $G$ harary1967number ; harary2014graphical . The intuition
behind this number is relatively straightforward: There are a total of $n!$
different ways to permute the $n$ vertices of graph $G$. Each permutation
relabels the nodes and thereby implicitly relabels the edges. Thus at first it
may seem that there are $n!$ different edge sets which form a graph $G$.
However, permutations which are _automorphisms_ re-arrange the nodes in such a
way that the overall labeled set of edges remains the same; that is precisely
the definition of an automorphism. Thus, for every distinct _edge set_ which
forms a graph isomorphic to $G$, there are $|\text{Aut}(G)|$ _permutations_ of
$G$’s nodes that yeild that _particular_ edge set. In sum, the full formula
above can be thought of as, “(the probability of a specific edge set) $\times$
(the number of distinct permutations of the nodes) $\times$ (the ratio of
distinct edge sets to distinct permutations [when the edge set in question is
isomorphic to $G$])”

abccbabaccabacbbca$\left.\vbox{\hrule
height=6.86601pt,depth=6.86601pt,width=0.0pt}\right\\}$a b c a b c 0 1 0 1 0 1
0 1 0 $\left[\vbox{\hrule
height=23.19867pt,depth=23.19867pt,width=0.0pt}\right.$$\left.\vbox{\hrule
height=23.19867pt,depth=23.19867pt,width=0.0pt}\right]$$\left.\vbox{\hrule
height=6.86601pt,depth=6.86601pt,width=0.0pt}\right\\}$a b c a b c 0 1 1 1 0 0
1 0 0 $\left[\vbox{\hrule
height=23.19867pt,depth=23.19867pt,width=0.0pt}\right.$$\left.\vbox{\hrule
height=23.19867pt,depth=23.19867pt,width=0.0pt}\right]$$\left.\vbox{\hrule
height=6.86601pt,depth=6.86601pt,width=0.0pt}\right\\}$a b c a b c 0 0 1 0 0 1
1 1 0 $\left[\vbox{\hrule
height=23.19867pt,depth=23.19867pt,width=0.0pt}\right.$$\left.\vbox{\hrule
height=23.19867pt,depth=23.19867pt,width=0.0pt}\right]$All Distinct Node Label
PermutationsAll Distinct Edge Sets Figure 2: Example of Counting the Number of
Distinct Edge Sets Forming a Structure \- This figure illustrates part of the
meaning and intuition behind Equation 1. The wedge graph has $n=3$ nodes, so
there are $n!=3!=6$ permutations of the three node labels $a,b,c$. The wedge
graph has two automorphisms ($|\text{Aut}(G)|=2$) Thus there are
$\frac{6}{2}=3$ distinct edge sets (i.e. distinct adjacency matrices) that
form a wedge graph.

To see this in an example, it may be helpful to think of edge sets in terms of
adjacency matrices. Each edge set $E$ can be thought of as an adjacency matrix
$A$ where $A_{ij}=\mathbbm{1}\left((i,j)\in E\right)$. A permutation $f$ on
the node labels corresponds to relabeling the columns _and_ rows of the
adjacency matrix; that is, $f(A)_{ij}=A_{f^{-1}(i)f^{-1}(j)}$. The permutation
$f$ is an automorphism of $A$ if $f(A)=A$. See Figure 2 for a visual
interpretation.

##  III Main Result

For any undirected graph $G=(V,E)$ and any subset of edges
$E^{\prime}\subseteq E$ where $|E^{\prime}|>0$, we prove that:

$\frac{|\text{Aut}(G)|}{|\text{AO}_{G}(E^{\prime})|}=\frac{|\text{Aut}(G-E^{\prime})|}{|\text{AO}_{G-E^{\prime}}(E^{\prime})|}$

Remember that $E^{\prime}$ is basically a set of node pairs and thus that
$E^{\prime}$ can both be thought of as a set of edges in graph $G$ and as a
set of non-edges in graph $G-E^{\prime}$. The above formula states that the
ratio of automorphisms of $G$ to the number of edge sets that $E^{\prime}$ is
equivalent to in $G$ is equal to the ratio of automorphisms of $G-E^{\prime}$
to the number of non-edge sets that $E^{\prime}$ is equivalent to in
$G-E^{\prime}$.

When $E^{\prime}$ contains a single edge $e$, we can simply write:

$\frac{|\text{Aut}(G)|}{|\text{AO}_{G}(e)|}=\frac{|\text{Aut}(G-e)|}{|\text{AO}_{G-e}(e)|}$

##  IV Proof

###  IV.1 Erdős-Reyni Generation and Edge Sampling Order

Our proof revolves around analyzing the probability that a graph $H$ generated
by the Erdős-Rényi model is isomorphic to some graph $G$ and the probability
that when $H$ was generated “the edges in $H$ corresponding to the edges
$E^{\prime}$ in $G$” were the last edges sampled. However, given that there is
not necessarily just one isomorphism between $G$ and $H$, there is not
necessarily a single set of edges in $H$ which uniquely corresponds to an edge
set $E^{\prime}$ in $G$.

Consequently, we consider the set of all possible correspondences. If the
generated graph $H$ is isomorphic to $G$, we can represent this as $H$ being
drawn uniformly from $[G]$, and then we consider a random correspondence
(isomorphism) $f:V\rightarrow V$ sampled uniformly from the isomorphisms
between $G$ and $H$.

We can consider the edges sampled during the Erdős-Rényi generation process as
having been sampled (without replacement) from the set of all possible edges
_one at a time_. Given any particular subset $s$ of the $m$ sampled edges, we
can observe that the probability it was the last edge sampled is
$\frac{1}{\binom{m}{|s|}}$. This is independent of the structure of the
generated graph $H$, since the _identity_ of the other sampled edges has
nothing to do with the conceptual _order_ in which they were sampled.

Given a set of edges $s$ in graph $H$, we will use the expression $l(s)$ to
denote the event that the edges in $s$ were the last $|s|$ edges sampled.

Similarly, if we are given a _set_ of edge sets $S$, we can talk about the
event that _one_ of the edge sets in $S$ was the collection of the last edges
sampled. We denote this as $l(S)$, where $l(S)=\bigvee_{s\in S}l(s)$.

###  IV.2 The Derivation

Note that _if_ $|E|=|E^{\prime}|$, then the result follows trivially because
$Aut(G-E^{\prime})$ would be the symmetric group on $n$ with size $n!$,
$|\text{AO}_{G}(E^{\prime})|$ would be 1, and
$|\text{AO}_{G-E^{\prime}}(E^{\prime})|$ would be the number of distinct edge
sets that produce the graph $G$, which as discussed in Section II.4 is
$\frac{n!}{|\text{Aut}(G)|}$. From here on out we assume that
$|E|>|E^{\prime}|$.

We can now ask the question: “Given a set of edges $E^{\prime}\subseteq E$,
what is the probability that $H$ is isomorphic to $G$ _and_ that under a
random isomorphism $f$ between $G$ and $H$, the edge set corresponding to
$E^{\prime}$ contains the last edges added?” This event can be expressed as
the event that $H\in[G]$ and the event that $l(f(E^{\prime}))$ given that
$H\in[G]$. Let $k=|E^{\prime}|$. Because the edge set that $f$ maps
$E^{\prime}$ to had a probability of $\frac{1}{\binom{m}{k}}$ of being sampled
last, this quantity is:

$\mathbb{P}(l(f(E^{\prime}))\ |\
H\in[G])\cdot\mathbb{P}(H\in[G])=\frac{1}{\binom{m}{k}}\mathbb{P}(H\in[G])$
(2)

However, we can also break the quantity down in a different way. Let
$H^{\prime}$ denote the graph generated by including all but the last $k$
sampled edges. Note that the following holds:

$H\in[G]\land l(f(E^{\prime}))\rightarrow H^{\prime}\in[G-E^{\prime}]$

That is, if $H$ is isomorphic to $G$ and $E^{\prime}$ were the last $k$ edges
added, then it follows by the definition of $H^{\prime}$ that $H^{\prime}$ was
isomorphic to $G-E^{\prime}$. In the same vein, note that the following
conditions also hold:

$l(f(E^{\prime}))\rightarrow l(f(\text{AO}_{G}(E^{\prime})))$

and

$H\in[G]\land l(f(\text{AO}_{G}(E^{\prime})))\rightarrow
H^{\prime}\in[G-E^{\prime}]$

These observations allow us to break down the event “$H\in[G]\land
l(f(E^{\prime}))$” into the equivalent event “$H\in[G]\land
l(f(E^{\prime}))\land l(f(\text{AO}_{G}(E^{\prime})))\land
H^{\prime}\in[G-E^{\prime}]$” so as to make use of standard probability
conditioning:

$\displaystyle\mathbb{P}(l(f(E^{\prime}))$ $\displaystyle|\
H\in[G])\cdot\mathbb{P}(H\in[G])$ (3)
$\displaystyle=\mathbb{P}(l(f(E^{\prime}))\ |\ H\in[G]\land
l(f(\text{AO}_{G}(E^{\prime})))\land H^{\prime}\in[G-E^{\prime}])$
$\displaystyle\ \ \ \ \cdot\mathbb{P}(H\in[G]\land
l(f(\text{AO}_{G}(E^{\prime})))\ |\ H^{\prime}\in[G-E^{\prime}])$
$\displaystyle\ \ \ \ \cdot\mathbb{P}(H^{\prime}\in[G-E^{\prime}])$
$\displaystyle=\frac{1}{|\text{AO}_{G}(E^{\prime})|}\cdot\frac{|\text{AO}_{G-E^{\prime}}(E^{\prime})|}{\binom{\binom{n}{2}-(m-k)}{k}}\cdot\mathbb{P}(H^{\prime}\in[G-E^{\prime}])$

The three final quantities in the last line of the equation correspond to the
three probabilities in the middle. To explain the first of these three
quantities, observe that:

$\displaystyle\mathbb{P}\left(l(f(E^{\prime}))\ |\ H\in[G]\land
l(f(\text{AO}_{G}(E^{\prime})))\land H^{\prime}\in[G-E^{\prime}]\right)$
$\displaystyle\ \ \ =\mathbb{P}\left(l(f(E^{\prime}))\ |\ H\in[G]\land
l(f(\text{AO}_{G}(E^{\prime})))\right)$ $\displaystyle\ \ \
=\frac{1}{|\text{AO}_{G}(E^{\prime})|}$

In other words, this is the probability that _if_ the final graph is
isomorphic to $G$ and _if_ the last $k$ edges added are structurally
equivalent to $E^{\prime}$ in $G$ (i.e. if the last set of the $k$ edges are
in $\text{AO}_{H}(f(E^{\prime}))$), that those $k$ edges in $H$ correspond
(via $f$) to $E^{\prime}$ in particular. Note that
$\text{AO}_{H}(f(E^{\prime}))=f(\text{AO}_{G}(E^{\prime}))$. The reason we can
let the “$H^{\prime}\in[G-E^{\prime}]$” drop out from the equation is that as
we discussed earlier, $H\in[G]\land l(f(\text{AO}_{G}(E^{\prime})))\rightarrow
H^{\prime}\in[G-E^{\prime}]$.

Further,
$\frac{|\text{AO}_{G-E^{\prime}}(E^{\prime})|}{\binom{\binom{n}{2}-(m-k)}{k}}=\mathbb{P}(H\in[G]\land
l(f(\text{AO}_{G}(E^{\prime})))\ |\ H^{\prime}\in[G-E^{\prime}])$. The
numerator $|\text{AO}_{G-E}(E^{\prime})|$ is the number of non-edge sets in
$G-E^{\prime}$ which are equivalent to the missing edge set $E^{\prime}$ (i.e.
the number of potential edge sets which, if added to a graph isomorphic to
$G-E^{\prime}$, will yield a graph $G^{\prime}$ isomorphic to $G$ _where_ that
added edge set plays the same role in $G^{\prime}$ that $E^{\prime}$ plays in
$G$). Because $H^{\prime}$ has $m-k$ edges, then
$\binom{\binom{n}{2}-(m-k)}{k}$ is the total number of $k$-size non-edge sets
in $H^{\prime}$ from which a $k$-size edge set will be uniformly sampled by
the ER generation process. Thus the ratio of $|\text{AO}_{G-E}(E^{\prime})|$
to $\binom{\binom{n}{2}-(m-k)}{k}$ is the probability that if $H^{\prime}$ is
isomorphic to $G-E^{\prime}$ and if $k$ edges are added to $H^{\prime}$, then
the resulting graph will be isomorphic to $G$ _and_ the $k$ added edges will
have the same role in $H$ as $E^{\prime}$ has in $G$.

Putting equations 2 and 3 together, we get that:

$\frac{1}{\binom{m}{k}}\cdot\mathbb{P}(H\in[G])=\frac{1}{|\text{AO}_{G}(E^{\prime})|}\cdot\frac{|\text{AO}_{G-E^{\prime}}(E^{\prime})|}{\binom{\binom{n}{2}-(m-k)}{k}}\cdot\mathbb{P}(H^{\prime}\in[G-E^{\prime}])$
(4)

Multiplying both sides by $\binom{m}{k}$ gives us:

$\mathbb{P}(H\in[G])=\frac{\binom{m}{k}}{|\text{AO}_{G}(E^{\prime})|}\cdot\frac{|\text{AO}_{G-E^{\prime}}(E^{\prime})|}{\binom{\binom{n}{2}-(m-k)}{k}}\cdot\mathbb{P}(H^{\prime}\in[G-E^{\prime}])$
(5)

If we substitute in the formula for the probability of generating a type of
graph according to the Erdős-Rényi model, this equation becomes:

$\frac{1}{\binom{\binom{n}{2}}{m}}\cdot\frac{n!}{|\text{Aut}(G)|}=\frac{\binom{m}{k}}{|\text{AO}_{G}(E^{\prime})|}\cdot\frac{|\text{AO}_{G-E^{\prime}}(E^{\prime})|}{\binom{\binom{n}{2}-(m-k)}{k}}\cdot\frac{1}{\binom{\binom{n}{2}}{m-k}}\cdot\frac{n!}{|\text{Aut}(G-E^{\prime})|}$
(6)

Re-arranging and expanding gives us:

$\displaystyle\frac{|\text{AO}_{G}(E^{\prime})|}{|\text{Aut}(G)|}$
$\displaystyle=\frac{\binom{\binom{n}{2}}{m}}{\binom{\binom{n}{2}}{m-k}}\cdot\frac{\binom{m}{k}}{\binom{\binom{n}{2}-(m-k)}{k}}\cdot\frac{|\text{AO}_{G-E^{\prime}}(E^{\prime})|}{|\text{Aut}(G-E^{\prime})|}$
(7)
$\displaystyle=\frac{\frac{\binom{n}{2}\cdot\left(\binom{n}{2}-1\right)\cdot...\cdot\left(\binom{n}{2}+1-m\right)}{m!}}{\frac{\binom{n}{2}\cdot\left(\binom{n}{2}-1\right)\cdot...\cdot\left(\binom{n}{2}+1-(m-k)\right)}{(m-k)!}}$
$\displaystyle\ \ \ \
\cdot\frac{\frac{m\cdot(m-1)\cdot...\cdot(m+1-k)}{k!}}{\frac{\left(\binom{n}{2}-(m-k)\right)\cdot\left(\binom{n}{2}-((m-k)+1)\right)\cdot...\cdot\left(\binom{n}{2}+1-m\right)}{k!}}\cdot\frac{|\text{AO}_{G-E^{\prime}}(E^{\prime})|}{|\text{Aut}(G-E^{\prime})|}$

To simplify this expression, we define the following shorthand:

  * •

$a=\binom{n}{2}\cdot\left(\binom{n}{2}-1\right)\cdot...\cdot\left(\binom{n}{2}+1-m\right)$

  * •

$b=m!$

  * •

$c=\binom{n}{2}\cdot\left(\binom{n}{2}-1\right)\cdot...\cdot\left(\binom{n}{2}+1-(m-k)\right)$

  * •

$d=(m-k)!$

  * •

$e=m\cdot(m-1)\cdot...\cdot(m+1-k)$

  * •

$f=k!$

  * •

$g=\left(\binom{n}{2}-(m-k)\right)\cdot\left(\binom{n}{2}-((m-k)+1)\right)\cdot...\cdot\left(\binom{n}{2}+1-m\right)$

  * •

$h=k!$

Using this shorthand, Equation 7 becomes:

$\frac{|\text{AO}_{G}(E^{\prime})|}{|\text{Aut}(G)|}=\frac{\frac{a}{b}}{\frac{c}{d}}\cdot\frac{\frac{e}{f}}{\frac{g}{h}}\cdot\frac{|\text{AO}_{G-E^{\prime}}(E^{\prime})|}{|\text{Aut}(G-E^{\prime})|}=\frac{adeh}{bcfg}\cdot\frac{|\text{AO}_{G-E^{\prime}}(E^{\prime})|}{|\text{Aut}(G-E^{\prime})|}$
(8)

Observe that $a$, $b$, $c$, $d$, $e$, and $g$ all take the form of
$x\cdot(x-1)\cdot...\cdot(x-y)$ for different values of $x$ and $y$. Further,
observe that $c$ and $g$ combine to form $a$: That is, $c$ and $a$ both begin
with $\binom{n}{2}$, but $c$ stops before $a$ does; then $g$ picks up where
$c$ left off, and $g$ and $a$ end with the same value. Thus $a=cg$. In the
same manner, $b=ed$. Given that $f=h$, these observations give us that
$\frac{adeh}{bcfg}=1$, and thus our result has been derived.

##  V Conclusion

###  V.1 Extensions

We expect that this result can be extended to other kinds of graphs, including
directed graphs and graphs with node colors. The formula for the probability
of generating a directed graph with the Erdős-Rényi model would be very
similar to the one for undirected graphs shown in Section II.4 because the
same logic about symmetry and edge sets applies.

The main limitation of our result is that it only relates the symmetries of
graphs where one graph’s edges is a superset of the other graph’s edges; it
would be a significant addition if this kind of result could be obtained for
any two graphs. Such a modification might simultaneously consider a set of
edge additions and a set of edge deletions.

###  V.2 Application to the Graph Reconstruction Conjecture

As we began to discuss in the introduction, we believe this result may have a
bearing on the graph reconstruction conjecture. We begin to explore that
connection here.

Given a graph $G=(V,E)$ and a vertex $v\in V$, define the set of edges
incident on $v$ to be $e_{G}(v)=\\{(a,b)\ |\ (a,b)\in E\land(a=v\lor b=v)\\}$
and define $G-v$ to be the graph $(V\setminus\\{v\\},E\setminus e_{G}(v))$.

In the graph reconstruction problem, the deck $\mathcal{D}$ is defined to be
the multiset $\\{G-v\ |\ v\in V\\}$. We define the augmented deck
$\mathcal{D}_{A}$ to be the multiset $\\{G~{}-~{}e_{G}(v)\ |\ v\in V\\}$. Note
that the two multisets are equivalent in the sense that if you are given
$\mathcal{D}$ you can obtain $\mathcal{D}_{A}$ by adding a singleton vertex to
each graph, and given $\mathcal{D}_{A}$ you can obtain $\mathcal{D}$ by
removing a singleton vertex from each graph. Thus the graph reconstruction
conjecture can be framed as the question, “Given the augmented deck
$\mathcal{D}_{A}$, can you uniquely determine $G$?”

The reconstruction conjecture only considers graphs with 3 or more nodes, and
we know that a graph can be reconstructed if it is disconnected bondy1977graph
. Thus the conjecture is only an open question for graphs which are connected
and have 3 or more nodes, which in turn entails that each vertex either has a
degree $>1$ or the _neighbor_ of each vertex of degree of 1 has a degree $>1$.
From here on out we assume that $G$ meets these basic conditions. These
conditions give us that no two vertices have the same edge set. Formally,
$u\neq v\rightarrow e_{G}(u)\neq e_{G}(v)$.

We also note that for any vertex $v\in V$, its edge set $e_{G}(v)$ can only be
mapped by automorphism to another edge set $E^{\prime}$ in $G$ if $v$ can be
mapped by automorphism to a vertex whose incident edge set is $E^{\prime}$.
Formally, this can be expressed as $f\in\text{Aut}(G)\rightarrow
f(e_{G}(v))=e_{G}(f(v))$.

These observations yeild that under automorphisms, $e_{G}(v)$ maps to all (and
only) other $e_{G}(u)$ where $u\in\text{AO}(v)$, and as we noted above, $u\neq
v\rightarrow e_{G}(u)\neq e_{G}(v)$. Consequently we get that
$|\text{AO}_{G}(e_{G}(v))|=|\text{AO}_{G}(v)|$. Further, for any graph
$G^{\prime}$ in the deck corresponding to deleting some vertex $v$, someone
with the full deck _knows_ the quantity $|\text{AO}_{G}(v)|$ because it is the
number of times the graph appears in the deck. Let $M(G^{\prime})$ denote this
quantity (the multiplicity of $G^{\prime}$ in $\mathcal{D}_{A}$).

Given a graph $G^{\prime}$ in the augmented deck, we let $v(G^{\prime})$
denote a (arbitrary) vertex for which $G-e_{G}(v)\cong G^{\prime}$. Applying
our result from this paper gives us that for any graph
$G^{\prime}\in\mathcal{D}_{A}$:

$\frac{|\text{Aut}(G)|}{|\text{AO}_{G}(e_{G}(v(G^{\prime})))|}=\frac{|\text{Aut}(G^{\prime})|}{|\text{AO}_{G^{\prime}}(e_{G}(v(G^{\prime})))|}$

Re arranging gives us:

$|\text{Aut}(G)|=\frac{|\text{Aut}(G^{\prime})|}{|\text{AO}_{G^{\prime}}(e_{G}(v(G^{\prime})))|}\cdot|\text{AO}_{G}(e_{G}(v(G^{\prime})))|$

As we said above, for relevant graphs
$|\text{AO}_{G}(e_{G}(v(G^{\prime})))|=|\text{AO}_{G}(v(G^{\prime}))|=M(G^{\prime})$,
which can be determined from the (augmented) deck. Thus, our equation becomes:

$|\text{Aut}(G)|=\frac{|\text{Aut}(G^{\prime})|}{|\text{AO}_{G^{\prime}}(e_{G}(v(G^{\prime})))|}\cdot
M(G^{\prime})$

This applies to the reconstruction conjecture as follows: We know that the
number of missing edges is reconstructible bondy1977graph . Given a graph
$G^{\prime}$ in the (augmented) deck, let $\mathcal{E}(G^{\prime})$ denote the
set of all edge sets that represent adding a vertex with the appropriate
degree. _If_ for a subset of the augmented deck
$\\{G_{1},G_{2},...,G_{k}\\}\subseteq\mathcal{D}_{A}$ we get that there are
respective (automorphically) unique sets of edges
$E_{1}\in\mathcal{E}(G_{1}),E_{2}\in\mathcal{E}(G_{2}),...,E_{k}\in\mathcal{E}(G_{k})$
where:

$\forall\ 1\leq i,j\leq k.\
\frac{M(G_{i})\cdot|\text{Aut}(G_{i})|}{|\text{AO}_{G_{i}}(E_{i})|}=\frac{M(G_{j})\cdot|\text{Aut}(G_{j})|}{|\text{AO}_{G_{j}}(E_{j})|}$

…Then we know there is only one way to add edges to the deck graphs which
yeilds the same amount of symmetry in all the resulting graphs, which means
$G$ can be uniquely reconstructed from the deck. As a reminder,
$|\text{AO}_{G_{i}}(E_{i})|$ represents the number of distinct sets of edges
that can be added to $G_{i}$ that will result in a graph isomorphic to $G_{i}$
with $E_{i}$ added.

###  V.3 Wrap-Up

We have proven an elegant combinatorial result on graphs relating the
symmetries of any two graphs differing by an edge set. Our proof utilizes the
Erdős-Rényi graph generation model, but the combinatorial result itself makes
no explicit reference to the model. We expect that due to its fundamental
nature and simplicity, our equation will have much to say about various
mathematical questions, perhaps including the graph reconstruction conjecture.

## Acknowledgements

This research was supported by a grant from the US National Science Foundation
(#1652492).

Thanks to Daniel Gonzalez Cedre and Tim Weninger for their useful feedback in
clarifying the narrative of this document.

## References

  * [1] Lars Døvling Andersen, Songkang Ding, Gert Sabidussi, and Preben Dahl Vestergaard.  Edge orbits and edge-deleted subgraphs.  Graphs and Combinatorics, 8(1):31–44, 1992.
  * [2] Lars Døvling Andersen, Songkang Ding, and Preben Dahl Vestergaard.  On the set edge-reconstruction conjecture.  JCMCC-Journal of Combinatorial Mathematics and Combinatorial Computing, 20:3–9, 1996.
  * [3] Béla Bollobás.  Almost every graph has reconstruction number three.  Journal of Graph Theory, 14(1):1–4, 1990.
  * [4] J Adrian Bondy and Robert L Hemminger.  Graph reconstruction—a survey.  Journal of Graph Theory, 1(3):227–268, 1977.
  * [5] Andrew Bowler, Paul Brown, Trevor Fenner, and Wendy Myrvold.  Recognizing connectedness from vertex-deleted subgraphs.  Journal of Graph Theory, 67(4):285–299, 2011.
  * [6] Paul Erdős, Alfréd Rényi, et al.  On the evolution of random graphs.  Publ. Math. Inst. Hung. Acad. Sci, 5(1):17–60, 1960.
  * [7] Alexander Farrugia.  Graphs having most of their eigenvalues shared by a vertex deleted subgraph.  Symmetry, 13(9):1663, 2021.
  * [8] DL Greenwell and RL Hemminger.  Reconstructing graphs.  In The many facets of graph theory, pages 91–114. Springer, 1969\.
  * [9] Frank Harary.  The four color conjecture and other graphical diseases.  Proof Techniques in Graph Theory, pages 1–9, 1969.
  * [10] Frank Harary and Edgar M Palmer.  Graphical enumeration.  Elsevier, 2014.
  * [11] Frank Harary, Edgar M Palmer, and Ronald C Read.  The number of ways to label a structure.  Psychometrika, 32(2):155–156, 1967.
  * [12] Stephen G Hartke, Hannah Kolb, Jared Nishikawa, and Derrick Stolee.  Automorphism groups of a graph and a vertex-deleted subgraph.  arXiv preprint arXiv:0909.3252, 2009.
  * [13] Steve Kirkland.  Algebraic connectivity for vertex-deleted subgraphs, and a notion of vertex centrality.  Discrete Mathematics, 310(4):911–921, 2010.
  * [14] William L Kocay.  A family of nonreconstructible hypergraphs.  Journal of Combinatorial Theory, Series B, 42(1):46–63, 1987.
  * [15] J Lauri.  Vertex-deleted and edge-deleted subgraphs.  A collection of papers by members of the University of Malta on the occasion of its quartercentenary celebrations (editors: R. Ellul-Micallef and S. Fiorini), Malta, 1992.
  * [16] Brian McMullen.  Graph reconstruction numbers.  2006\.
  * [17] Edward F Schmeichelm and Gary S Bloom.  Connectivity, genus, and the number of components in vertex-deleted subgraphs.  Journal of Combinatorial Theory, Series B, 27(2):198–201, 1979\.


published online,SuperConga team (2022a) containing numerous pedagogical
examples, tutorials and guides. The framework has been developed and tested
for Unix-based environments, and generally runs on most modern laptops,
desktops, as well as in cluster environments. The framework relies on the
high-performance capabilities offered by GPU acceleration and CUDANVIDIA
Corporation (2021a); Nickolls  _et al._ (2008), and therefore requires a CUDA-
capable device (i.e. NVIDIA GPU)NVIDIA Corporation (2021b) to run. To the best
of our knowledge, SuperConga is the only open-source code that uses
quasiclassical theory to describe mesoscopic superconducting grains and
ballistic devices. It therefore adds the missing “mesoscopic link” between
more phenomenological methods, such as London-MaxwellHorn and Moler (2022) and
Ginzburg-Landau,Sadovskyy  _et al._ (2015); Rybakov and Babaev (2022, ) and
fully microscopic methods, such as density-functional theoryKresse and
Furthmüller (1996, 1996); Giannozzi  _et al._ (2009, 2017); Nakhaee, Ketabi,
and Peeters (2020) and tight-binding approaches.Groth  _et al._ ; Björnson
(2019); Löthman  _et al._ (2021); Nagai (2020)

This paper outlines the functionality of version 1.0 of SuperConga. We foresee
that the project website will evolve continuously, and that a community of
developers can be formed to enable further improvements that increase the
functionality and scope of the framework.

###### Acknowledgements.

N.W-W. and P.H. have contributed equally to this work. We thank Anton B.
Vorontsov and Kevin Marc Seja for valuable discussions. We thank the Swedish
Research Council (VR) for financial support. The computations were enabled by
resources provided by the Swedish National Infrastructure for Computing (SNIC)
at the computing centers C3SE, HPC2N and NSC, partially funded by the Swedish
Research Council through grant agreement no. 2018-05973.

##  Appendix A Dependencies

This appendix lists external dependencies and packages that SuperConga
offloads functionality to. For a more detailed description, please refer to
the installation guide in the main README of the repository.SuperConga team
(2022b) SuperConga uses git as a version control system, with the central
repository being hosted on GitLab.SuperConga team (2022b) The SuperConga
backend is written in C++, with most of the heavy lifting being done by
CUDANVIDIA Corporation (2021a); Nickolls  _et al._ (2008) kernels and Thrust
Bell and Hoberock (2012); Kaczmarski and Rzkażewski (2013); Wynters (2013)
algorithms running on the GPU. This functionality is provided by the CUDA
toolkit, which also provides the NVCC compiler for compiling the CUDA code.
Alternatively, Clang can be used for compilation. To run the compiled code, a
CUDA-enabled device (an NVIDIA GPU) is required. Note that we use
CMakeCedilnik, E. and Hoffman, B. and King, B. and Martin, K. and Neundorf, A.
(2021) as the build-system generator, and NinjaMartin, E. (2020) as the build
system, but that it is possible to replace either of these with another
choice. Linear algebra on the CPU uses ArmadilloSanderson and Curtin (2016,
2018), and all tests are implemented in the doctestKirilov (2021) framework.
The real-time plotting is done using OpenGLKhronos Group (2021); Woo  _et al._
(1999) and ArrayFire-ForgeYalamanchili  _et al._ (2015). The package
JsonCppDunn (2020) is used to handle JSONPezoa  _et al._ (2016) objects on the
backend side, and HighFiveBlueBrain (2021) for HDF5The HDF Group (NNNN)
support. The SuperConga frontend requires a Python3 environment, in particular
with the modules NumPyHarris  _et al._ (2020), MatplotlibHunter (2007),
h5pyCollette (2013) and pandasMcKinney  _et al._ (2010).

##  Appendix B The simulation-parameter file

The simulation configuration is specified in a simulation_config.json file. It
is a JSONPezoa  _et al._ (2016) file, which consists of a number of key-value
pairs on the form "key": value. The value of a key can be more key-value
pairs, numbers, strings, etc. Below is the JSON file with the parameters used
in the demonstration example in Section IV.

Listing 10: examples/swave_disc_vortex/simulation_config.json.

⬇

{

"physics": {

"temperature": 0.5,

"external_flux_quanta": 1.5,

"penetration_depth": 5.0,

"crystal_axes_rotation": 0.0,

"gauge": "symmetric",

"charge_sign": -1,

"order_parameter": {

"s": {

"critical_temperature": 1.0,

"initial_phase_shift": 0.0,

"initial_noise_stddev": 0.0,

"vortices": [

{

"center_x": 0.0,

"center_y": 0.0,

"winding_number": -1.0

}

]

}

}

},

"geometry": [

{

"disc": {

"add": true,

"center_x": 0.0,

"center_y": 0.0,

"radius": 15.0

}

}

],

"numerics": {

"convergence_criterion": 1e-05,

"energy_cutoff": 16.0,

"norm": "l2",

"num_fermi_momenta": 32,

"num_iterations_burnin": 0,

"num_iterations_max": 10000,

"num_iterations_min": 0,

"points_per_coherence_length": 20.0,

"vector_potential_error": 1e-06

},

"accelerator": {

"name": "congacc"

},

"misc": {

"data_format": "h5",

"load_path": "",

"save_frequency": -1,

"save_path": "data/examples/swave_disc_vortex",

"verbose": true,

"visualize": true

}

}

Under physics are all parameters directly related to the physics of the
simulation, except the geometry which is specified separately. The temperature
is given in units of the highest critical temperature of the order-parameter
components, $T_{\mathrm{c}}\equiv\max\\{T_{\mathrm{c},\Gamma}\\}$. The
external_flux_quanta is given in units of $\Phi_{0}$. The penetration_depth is
given in units of $\xi_{0}$. If the penetration depth is zero or negative it
is considered to be infinite and the induced vector potential will be zero.
The crystal_axes_rotation, i.e. the rotation of the crystal $ab$-axes relative
to the simulation coordinate system, is given in units of a full turn, i.e.
$2\pi$. The choices for the gauge is either symmetric, landau or solenoid. The
sign of the carrier charge is controlled by charge_sign.

The field order_parameter is a list of components. The supported order-
parameter components are $d_{x^{2}-y^{2}}$, $d_{xy}$, $g_{xy(x^{2}-y^{2})}$,
and $s$. The critical_temperature is only meaningful with multiple components
as it is the _relative_ $T_{\mathrm{c},\Gamma}$ that is relevant. Each
component is initialized as

$\displaystyle\Delta_{\Gamma}({\bm{R}})$ $\displaystyle=$
$\displaystyle\Delta^{(\mathrm{bulk})}_{\Gamma}e^{i\chi_{\Gamma}({\bm{R}})}+Z\;,$
(76) $\displaystyle\chi_{\Gamma}({\bm{R}})$ $\displaystyle=$ $\displaystyle
2\pi\delta_{\Gamma}+\sum_{v}w^{(v)}_{\Gamma}\phi\left({\bm{R}}-{\bm{C}}^{(v)}_{\Gamma}\right)\;,$
(77) $\displaystyle Z$ $\displaystyle\sim$
$\displaystyle\mathcal{N}(0,\sigma^{2}_{\Gamma})\;,$ (78)

where $\delta_{\Gamma}$ is controlled by initial_phase_shift, and $Z$ is
complex additive Gaussian noise with a standard deviation $\sigma_{\Gamma}$
controlled by initial_noise_stddev in units of $2\pi
k_{\mathrm{B}}T_{\mathrm{c}}$. The field vortices is list of vortices, with
each vortex being specified by its center position,
${\bm{C}}_{\Gamma}=(\mathtt{center\\_x},\mathtt{center\\_y})$, and phase
winding $w_{\Gamma}=\mathtt{phase\\_winding}$. Furthermore,
$\phi({\bm{R}})=\mathrm{atan2}(y,x)$ is the polar angle of the coordinate
${\bm{R}}$. We may add additional components by adding entries to the
"order_parameter" field in the above JSON file. For example, we may add a
subdominant $d_{xy}$-wave component initialized with a phase shift of $\pi/2$,
some Gaussian noise, and no vortices:

Listing 11: Example of a multi-component order parameter.

⬇

"order_parameter": {

"s": {

"critical_temperature": 1.0,

"initial_phase_shift": 0.0,

"initial_noise_stddev": 0.0,

"vortices": [

{

"center_x": 0.0,

"center_y": 0.0,

"winding_number": -1.0

}

]

},

"dxy": {

"critical_temperature": 0.1,

"initial_phase_shift": 0.25,

"initial_noise_stddev": 0.01,

"vortices": [

{

"center_x": 0.0,

"center_y": 0.0,

"winding_number": -1.0

}

]

}

}

The geometry is built of three different kinds of components; disc,
regular_polygon, and polygon. A disc is specified by its center which is set
by center_x and center_y, and its radius. All parameters are in units of
$\xi_{0}$. A regular_polygon is specified by its center coordinates center_x
and center_y, the number of edges num_edges, its rotation, and the
side_length. The side length and center is given in units of $\xi_{0}$, and
the rotation in units of a full turn, i.e. $2\pi$. A polygon is specified by
the $x$\- and $y$-coordinates of its vertices, vertices_x and vertices_y
respectively. The vertices should be given in counter-clockwise order. All
components has a parameter add, which determines if the area of the component
should be added (true), or subtracted (false), from the geometry. For more
general geometries it is possible to add and subtract several objects by
adding to the geometry field. For example, we may put a square hole of side
length $5\xi_{0}$ in the disk at $10\xi_{0}$ to the right of the center:

Listing 12: Example of a compound geometry.

⬇

"geometry": [

{

"disc": {

"add": true,

"center_x": 0.0,

"center_y": 0.0,

"radius": 15.0

}

},

{

"regular_polygon": {

"add": false,

"center_x": 10.0,

"center_y": 0.0,

"num_edges": 4,

"rotation": 0,

"side_length": 5.0

}

}

Note that the order of the components is important if any component is
subtracted.

Under numerics are parameters controlling the numerics of the simulation. The
convergence_criterion determines when the simulation is considered converged.
When all computed quantities; each order-parameter component, the charge
current density, the induced vector potential, and the free energy, has a
residual (as given by Eq. (27)) smaller than convergence_criterion the
simulation should stop. The energy_cutoff is the largest energy, in units of
$2\pi k_{\mathrm{B}}T_{\mathrm{c}}$, to include in the energy sums. The
$p$-norm to use when computing the residuals is controlled by norm, with
$p\in\\{1,2,\infty\\}$. The num_fermi_momenta is the number of discrete
momenta to use when approximating the Fermi-surface averages
$\langle\ldots\rangle_{{\bm{p}}_{\mathrm{\mbox{$\scriptscriptstyle{\rm{F}}$}}}}$.
Note that the terms in the Fermi-surface averages are computed sequentially in
contrast to the terms in the energy sum which are computed in parallel on the
GPU.

The field num_iterations_burnin controls the number of iterations during which
the order-parameter and vector potential are kept constant, while the
coherence functions are converged and updated self-consistently, in particular
on the boundary, i.e. $\gamma_{\partial\mathcal{D}}$ and
$\tilde{\gamma}_{\partial\mathcal{D}}$. Once the burn-in period is over, all
quantities are updated self-consistently again. This is particularly useful
when resuming a previous simulation from file, since the order parameter and
vector potentials are stored, but the coherence functions are not (due to
occupying too much disk space). If the number of burn-in iterations is
negative, it will run until convergence. The maximum and minimum of iterations
to run is controlled by num_iterations_max and num_iterations_min,
respectively. When num_iterations_max iterations has been run the simulation
stops regardless of whether convergence has been achieved. Similarily, the
simulation does _not_ stop until num_iterations_min iterations has been run
even if the simulation is considered converged. This is useful for avoiding
endless simulations, or prematurely stopping a simulation that is moving
slowly on a plateau, respectively. The spatial resolution of the simulation is
set by points_per_coherence_length. The vector_potential_error is explained in
App. D.

All parameters of the accelerator are explained in App. E.

Lastly, under misc are miscellaneous parameters not affecting the simulation
_per se_. data_format specifies which data format to use when saving files.
The options are h5 and csv yielding compressed HDF5The HDF Group (NNNN) files
or plain-text CSV files, respectively. Both data formats are standard and have
excellent support in e.g. MATLABMATLAB (2021) and Python (via h5pyCollette
(2013) and pandasMcKinney  _et al._ (2010)). By setting load_path SuperConga
will read the files located there and use them as the initial values for the
order parameter and the vector potential. How often to save the results are
controlled by save_frequency. If save_frequency is negative the results will
only be saved at the end of the simulation, and if it is zero the results will
not be saved at all, otherwise it will save every $n$th iteration. The results
are saved to save_path. If verbose is true the progress of the simulation will
be written to the terminal, and if visualize is true the simulation will be
visualized live.

##  Appendix C Solving the Riccati equations

Using the trajectory coordinate, $s$, Eq. (25) in standard form is

$\displaystyle\frac{\partial\gamma}{\partial s}$ $\displaystyle=$
$\displaystyle f(s,\gamma(s))\;,$ (79) $\displaystyle f(s,\gamma(s))$
$\displaystyle=$
$\displaystyle\frac{i}{\mathnormal{v}_{\mathrm{\mbox{$\scriptscriptstyle{\rm{F}}$}}}}\left[\Delta^{*}(s)\gamma^{2}(s)+2z(s)\gamma(s)+\Delta(s)\right]\;,$
(80)

with the dependence on
${\bm{p}}_{\mathrm{\mbox{$\scriptscriptstyle{\rm{F}}$}}}$ and $\varepsilon$
dropped, and

$z(s)\equiv\varepsilon+\tfrac{e}{c}{\bm{v}}_{\mathrm{\mbox{$\scriptscriptstyle{\rm{F}}$}}}\cdot{\bm{A}}(s)\;,$
(81)

for brevity. Instead of using the stepping method,Eschrig (2009); Grein,
Löfwander, and Eschrig (2013) i.e. using an analytical formula for $\gamma(s)$
but approximating $\Delta(s)$ and $z(s)$ as piecewise constant, we solve Eq.
(79) along a trajectory using the implicit mid-point method, which we find to
be both faster and more accurate. The implicit mid-point method has the
following update rule,

$\displaystyle\gamma_{j+1}$ $\displaystyle=$
$\displaystyle\gamma_{j}+hf\left(s_{j+\frac{1}{2}},\frac{1}{2}\left(\gamma_{j}+\gamma_{j+1}\right)\right)\;,$
(82)

with the discretized trajectory coordinate $s_{j}=jh$, where $j$ is an integer
and $h$ is the length between points. This yields a quadratic equation in
$\gamma_{j+1}$ to solve for each step along the trajectory;

$\displaystyle 0$ $\displaystyle=$ $\displaystyle
c_{2}\gamma^{2}_{j+1}+c_{1}\gamma_{j+1}+c_{0}\;,$ (83) $\displaystyle c_{2}$
$\displaystyle=$
$\displaystyle\frac{ih}{4\mathnormal{v}_{\mathrm{\mbox{$\scriptscriptstyle{\rm{F}}$}}}}\Delta^{*}_{j+\frac{1}{2}}\;,$
(84) $\displaystyle c_{1}$ $\displaystyle=$ $\displaystyle
2c_{2}\gamma_{j}+\frac{ih}{\mathnormal{v}_{\mathrm{\mbox{$\scriptscriptstyle{\rm{F}}$}}}}z_{j+\frac{1}{2}}-1\;,$
(85) $\displaystyle c_{0}$ $\displaystyle=$
$\displaystyle\left(c_{1}-c_{2}\gamma_{j}+2\right)\gamma_{j}-4c^{*}_{2}\;,$
(86)

with the solution

$\displaystyle\gamma_{j+1}$ $\displaystyle=$
$\displaystyle\frac{2c_{0}}{-c_{1}+\sqrt{c^{2}_{1}-4c_{2}c_{0}}}\;.$ (87)

Half-step values, $\Delta_{j+\frac{1}{2}}$ and $z_{j+\frac{1}{2}}$, are
calculated by linear interpolation. Nearest-neighbour extrapolation is used at
the boundary. If we have $\gamma(s_{j})$ at one point we obtain
$\gamma(s_{j+1})$ via Eq. (87). The solution along the whole trajectory is
obtained by starting with the initial boundary value
$\gamma(s_{0}={\bm{R}}_{\mathrm{min}})$ and stepping along the trajectory to
obtain $\gamma(s_{j})$ with $j>0$. The stepping proceeds until a second
boundary is reached (denoted ${\bm{R}}_{\mathrm{max}}$, see Fig 1). We compute
$\tilde{\gamma}(s)$ by solving Eq. (26) in the same manner, but stepping in
the opposite direction from ${\bm{R}}_{\mathrm{max}}$ to
${\bm{R}}_{\mathrm{min}}$.

##  Appendix D Computing the vector potential

The induced vector potential is obtained by solving Eq. (7), in the Coloumb
gauge, using the Green’s function of the 2D Laplacian operator. The current is
computed on a $N\times N$ grid, with every grid cell having the area $h^{2}$.
We treat the current as being piecewise constant within each grid cell. With
our choice of units, this yields

$\displaystyle{\bm{A}}^{\mathrm{ind}}_{m,n}$ $\displaystyle=$
$\displaystyle-\frac{1}{\kappa^{2}}\sum^{N-1}_{m^{\prime}=0}\sum^{N-1}_{n^{\prime}=0}G_{|m^{\prime}-m|,|n^{\prime}-n|}{\bm{j}}_{m^{\prime}n^{\prime}}\;,$
(88) $\displaystyle G_{m,n}$ $\displaystyle=$
$\displaystyle\frac{1}{4\pi}\int^{+\frac{1}{2}h}_{-\frac{1}{2}h}\int^{+\frac{1}{2}h}_{-\frac{1}{2}h}\ln\left(\frac{\left(x^{\prime}+hm\right)^{2}+\left(y^{\prime}+hn\right)^{2}}{\xi^{2}_{0}}\right)\mathrm{d}x^{\prime}\mathrm{d}y^{\prime}$
(89) $\displaystyle=$
$\displaystyle\frac{1}{4\pi}\left(\frac{h}{\xi_{0}}\right)^{2}\left(2\ln\left(\frac{h}{\xi_{0}}\right)-\ln\left(2\right)-3+\frac{1}{4}\sum_{a=\pm}\sum_{b=\pm}\Gamma(am,bn)\right)\;,$
(90) $\displaystyle\Gamma(m,n)$ $\displaystyle=$
$\displaystyle\alpha(m,n)+\beta(m,n)+\beta(n,m)\;,$ (91)
$\displaystyle\alpha(m,n)$ $\displaystyle=$
$\displaystyle\left(1+2m\right)\left(1+2n\right)\ln\left(\frac{1}{2}\left[1+2m\right]^{2}+\frac{1}{2}\left[1+2n\right]^{2}\right)\;,$
(92) $\displaystyle\beta(m,n)$ $\displaystyle=$
$\displaystyle\left(-2m^{2}+2n^{2}+4n+1\right)\tan^{-1}\left(\frac{1+2m}{1+2n}\right)\;.$
(93)

$G$ is a normal matrix, so it can be eigendecomposed as $G=V\Lambda V^{\top}$,
where $V$ is a unitary matrix with (normalized) eigenvectors as its columns,
and $\Lambda$ is a diagonal matrix with the eigenvalues as diagonal elements.
Equivalently,

$\displaystyle G_{m,n}$ $\displaystyle=$
$\displaystyle\sum^{N-1}_{k=0}\lambda_{k}V^{(k)}_{m}V^{(k)}_{n}\;,$ (94)

where $V^{(k)}_{m}$ is the eigenvector corresponding to the eigenvalue
$\lambda_{k}$. Thus

$\displaystyle{\bm{A}}^{\mathrm{ind}}_{m,n}$ $\displaystyle=$
$\displaystyle-\frac{1}{\kappa^{2}}\sum^{N-1}_{k=0}\sum^{N-1}_{m^{\prime}=0}\sum^{N-1}_{n^{\prime}=0}\lambda_{k}V^{(k)}_{|m^{\prime}-m|}V^{(k)}_{|n^{\prime}-n|}{\bm{j}}_{m^{\prime}n^{\prime}}\;,$
(95)

and the cross-correlation with the eigenvectors can be done separately. The
eigendecomposition is done with ArmadilloSanderson and Curtin (2016, 2018). In
order to increase performance only $N^{*}$ eigenvalues are kept. It is
computed as the smallest integer obeying

$\displaystyle\sqrt{1-\frac{\sum^{N^{*}-1}_{n=0}\lambda^{2}_{n}}{\sum^{N-1}_{n=0}\lambda^{2}_{n}}}$
$\displaystyle<$ $\displaystyle\epsilon\;,$ (96)

with the eigenvalues sorted
$|\lambda_{0}|\geq|\lambda_{1}|\geq...\geq|\lambda_{N-1}|$, and $\epsilon>0$
being a tolerance chosen by the user via the parameter vector_potential_error
in the configuration file. If $\epsilon=0$ no eigendecomposition is done, and
the current density is cross-correlated with the full 2D Green’s function.

Lastly, the mean of the charge-current density is in general not zero due to
numerical errors or bad initialization. This has the undesired effect of
introducing a constant term in the induced vector potential, which gives rise
to an overall phase gradient through the grain and an associated current. For
small penetration depths, this can prevent the simulation from converging. In
order to remedy this, the mean of the charge-current density is subtracted
before solving Eq. (7).

##  Appendix E Convergence accelerators

SuperConga solves Eqs. (20) and (7) self-consistently. With our choice of
units, the solution to Eq. (7) is given by

$\displaystyle\kappa^{2}{\bm{A}}_{\mathrm{ind}}({\bm{R}})$ $\displaystyle=$
$\displaystyle-\int_{\mathcal{A}}\mathrm{d}{\bm{R}}^{\prime}\,G({\bm{R}},{\bm{R}}^{\prime}){\bm{j}}({\bm{R}}^{\prime})\;,$
(97)

in the Coloumb gauge, where $G({\bm{R}},{\bm{R}}^{\prime})$ is the Green’s
function of the 2D Laplacian. In order to simplify the notation we concatenate
the left- and right-hand sides of Eqs. (20) and (97), yielding

$\displaystyle{\bm{x}}$ $\displaystyle=$ $\displaystyle{\bm{g}}({\bm{x}})\;,$
(98) $\displaystyle{\bm{x}}$ $\displaystyle\equiv$
$\displaystyle\left(\mathrm{Re}\left[\Delta_{\Gamma_{1}}\right],\mathrm{Im}\left[\Delta_{\Gamma_{1}}\right],\ldots,\kappa^{2}A^{\mathrm{ind}}_{x},\kappa^{2}A^{\mathrm{ind}}_{y}\right)^{\top}\;.$
(99)

Note that the vector potential is scaled by $\kappa^{2}$ in order to have
similar magnitudes of all elements in ${\bm{x}}$. SuperConga provides several
different methods of solving Eq. (98). Namely, basic Picard iterations,
Polyak’sPolyak (1964) "small heavy sphere", a variant of the Barzilai-Borwein
methodBarzilai and Borrwein (1988), and a custom method. The user controls
which method to use by setting the name in the accelerator in Listing. 10. By
only setting the name, all internal parameters of the accelerator will be set
to default values. The user can change the internal parameters by specifying
them in the configuration file, or via the CLI.

In the following discussion of the different accelerators, it is convenient to
define the difference between the right- and left-hand sides of Eq. (98),

${\bm{d}}\equiv{\bm{g}}({\bm{x}})-{\bm{x}}\;,$ (100)

which enters in all accelerator methods.

###  E.1 Picard

The simplest way of solving Eq. (98) is by using Picard iterations,

$\displaystyle{\bm{x}}_{i+1}$ $\displaystyle=$
$\displaystyle{\bm{x}}_{i}+\alpha{\bm{d}}_{i}\;,$ (101)

where $\alpha>0$ is the the step size. Picard iterations can converge very
slowly, or not at all. Hence, the need for a method accelerating the
iterations, and preferably stabilizing them. Calling the Picard method an
accelerator is a misnomer as it is the baseline. The other methods should
improve on this method. How to use this accelerator is shown in Listing 13.

Listing 13: Picard accelerator with default values.

⬇

"accelerator": {

"name": "picard",

"step_size": 1.0

}

###  E.2 Polyak

In Polyak’s method, ${\bm{x}}$ behaves similarly to the position of a particle
moving through a viscous fluid, with a low Reynold’s number, and a potential
field (Eq. (100) taking the role of the negative potential gradient),

$\displaystyle{\bm{v}}_{i+1}$ $\displaystyle=$
$\displaystyle(1-\beta){\bm{v}}_{i}+\alpha{\bm{d}}_{i}\;,$ (102)
$\displaystyle{\bm{x}}_{i+1}$ $\displaystyle=$
$\displaystyle{\bm{x}}_{i}+{\bm{v}}_{i+1}\;,$ (103)

where $\beta\in(0,1)$ is the drag, and $\alpha>0$ is the step size. The
velocity is initialized to zero, ${\bm{v}}_{0}=\mathbf{0}$. How to use this
accelerator is shown in Listing 14.

Listing 14: Polyak accelerator with default values.

⬇

"accelerator": {

"name": "polyak",

"step_size": 2.0,

"drag": 0.5

}

###  E.3 Barzilai-Borwein

The Barzilai-Borwein (BB) methodBarzilai and Borrwein (1988) is similar to the
Picard method but with an adaptive step size,

$\displaystyle s^{\mathrm{BB1}}_{i}$ $\displaystyle=$
$\displaystyle\frac{\|{\bm{d}}_{i-1}\|^{2}_{2}}{\left|{\bm{d}}_{i-1}\cdot\left({\bm{d}}_{i}-{\bm{d}}_{i-1}\right)\right|}\;,$
(104) $\displaystyle s^{\mathrm{BB2}}_{i}$ $\displaystyle=$
$\displaystyle\frac{\left|{\bm{d}}_{i-1}\cdot\left({\bm{d}}_{i}-{\bm{d}}_{i-1}\right)\right|}{\|{\bm{d}}_{i}-{\bm{d}}_{i-1}\|^{2}_{2}}\;,$
(105) $\displaystyle\alpha_{i+1}$ $\displaystyle=$
$\displaystyle\alpha_{i}\cdot\left\\{\begin{matrix}s^{\mathrm{BB1}}_{i}&\mathrm{if}\;i\;\mathrm{odd}\\\
s^{\mathrm{BB2}}_{i}&\mathrm{if}\;i\;\mathrm{even}\end{matrix}\right.\;,$
(106) $\displaystyle{\bm{x}}_{i+1}$ $\displaystyle=$
$\displaystyle{\bm{x}}_{i}+\alpha_{i+1}{\bm{d}}_{i}\;,$ (107)

where $s^{\mathrm{BB1}}$ and $s^{\mathrm{BB2}}$ are two different variants of
how to scale the step size. Note that we force the step size to be positive.
Alternating between the variants gives better performance on the examples
provided by SuperConga than using either one of them on its own. What step
size to use during the first iteration is a free parameter and can be set by
the user. How to use this accelerator is shown in Listing 15.

Listing 15: Barzilai-Borwein accelerator with default values.

⬇

"accelerator": {

"name": "barzilai-borwein",

"step_size": 1.0,

"step_size_max": 100.0,

"step_size_min": 0.001

}

###  E.4 CongAcc

CongAcc, short for (Super)Conga Accelerator, is our custom made accelerator.
Introducing the _component_ index, $c$, where a component is either the real
or imaginary part of a order-parameter component, or the $x$\- or
$y$-component of the induced vector potential, in Eq. (99), the update rule is
as follows,

$\displaystyle s_{i,c}$ $\displaystyle=$ $\displaystyle
S_{C}({\bm{m}}_{i,c},{\bm{d}}_{i,c})\;,$ (108) $\displaystyle\alpha_{i+1,c}$
$\displaystyle=$ $\displaystyle\alpha_{i,c}k^{s_{i,c}}\;,$ (109)
$\displaystyle w_{i+1,c}$ $\displaystyle=$
$\displaystyle\left\\{\begin{matrix}w_{i,c}&\mathrm{if}\;s_{i,c}\geq
s_{\mathrm{min}}\\\ w^{*}&\mathrm{otherwise}\;\end{matrix}\right.\;,$ (110)
$\displaystyle{\bm{m}}_{i+1,c}$ $\displaystyle=$ $\displaystyle
w_{i+1,c}{\bm{m}}_{i,c}+(1-w_{i+1,c}){\bm{d}}_{i,c}\;,$ (111)
$\displaystyle{\bm{x}}_{i+1,c}$ $\displaystyle=$
$\displaystyle{\bm{x}}_{i,c}+\alpha_{i+1,c}{\bm{m}}_{i+1,c}\;,$ (112)

where $S_{C}$ is the cosine similarity, ${\bm{m}}_{i,c}$ is an exponential
moving average of the component difference ${\bm{d}}_{i,c}$, and $k$ is a
constant determining how much the step size should maximally increase
(decrease) which can be set by the user. If $s_{i,c}<s_{\mathrm{min}}$, where
$s_{\mathrm{min}}$ is the minimum tolerated similarity, the weight $w^{*}$ is
obtained by solving

$\displaystyle
S_{C}(w^{*}{\bm{m}}_{i,c}+(1-w^{*}){\bm{d}}_{i,c},{\bm{d}}_{i,c})$
$\displaystyle=$ $\displaystyle s_{\mathrm{min}}\;$ (113)

ensuring that the step will be similar to ${\bm{d}}_{i,c}$. Equation (113) is
solved approximately using the bisection method. Note that we set $w_{0,c}=0$
during the first iteration.

The reasoning is that if ${\bm{m}}_{i,c}$ and ${\bm{d}}_{i,c}$ roughly point
in the same direction, then we should have taken a larger step the previous
iteration, and thus the step size is increased. Analogously, if they roughly
point in opposite directions, the step size should decrease.

How to use this accelerator is shown in Listing 16.

Listing 16: CongAcc accelerator with default values.

⬇

"accelerator": {

"name": "congacc",

"step_size": 0.5,

"step_size_factor": 1.234,

"cos_similarity_min": 0.7,

"step_size_max": 100.0,

"step_size_min": 0.001

}

This accelerator should be regarded as experimental as no thorough analysis of
it has been done.

###  E.5 Comparing the accelerators

SuperConga examples and accelerator comparison
---
| Name | $T/T_{\mathrm{\mbox{$\scriptscriptstyle{\rm{c}}$}}}$ | $\Phi/\Phi_{0}$ | $\kappa$ | Symmetry | Initialization | Geometry | Picard | Polyak | BB | CongAcc
(a) | dwave_chiral | $0.5$ | $0$ | $\infty$ | $d_{xy}+id_{x^{2}-y^{2}}$ | Bulk | Disc: $\mathcal{R}=12.5\xi_{0}$ | $80$ | $34$ | $25$ | $37$
(b) | dwave_octagon | $0.5$ | $0$ | $\infty$ | $d_{x^{2}-y^{2}}$ | Bulk | Octagon: $\mathcal{S}=10\xi_{0}$ | $71$ | $33$ | $25$ | $32$
(c) | dwave_plus_swave | $0.5$ | $0$ | $\infty$ | $d_{xy}+s$ | Bulk | Disc: $\mathcal{R}=12.5\xi_{0}$ | $69$ | $36$ | $30$ | $33$
(d) | dwave_phase_crystal |  $0.1$111The temperature is low, $T<T^{*}$, so the solution is a phase crystal.Håkansson, Löfwander, and Fogelström (2015); Holmvall  _et al._ (2020) | $0$ | $100$ | $d_{x^{2}-y^{2}}$ | Bulk + vortex222The vortex has a phase winding $n=-1$ and is positioned outside of the grain, yielding a soft phase-gradient through the grain. | Irregular polygon333The shape of the grain is a square, $\mathcal{S}=25\xi_{0}$, with one corner removed, yielding four $[100]$-interfaces, and one pair-breaking $[110]$-interface.Wennerdal  _et al._ (2020) | $3027$ | $809$ | $1054$ | $439$
(e) | swave_abrikosov_lattice | $0.2$ | $20$ | $10$ | $s$ | Bulk + giant vortex444The giant vortex in the middle has a phase winding $n=-13$. It is unstable and decays into a lattice of singly-quantized Abrikosov vortices. | Square: $\mathcal{S}=25\xi_{0}$ | $1001$ | $313$ | $211$ | $133$
(f) | swave_disc_meissner | $0.5$ | $0.5$ | $5$ | $s$ | Bulk | Disc: $\mathcal{R}=15\xi_{0}$ | F | F | $20$ | $43$
(g) | swave_disc_vortex | $0.5$ | $1.5$ | $5$ | $s$ | Bulk + vortex555The vortex in the middle has a phase winding $n=-1$ | Disc: $\mathcal{R}=15\xi_{0}$ | F | F | $45$ | $42$
Table 1: How many iterations each accelerator, with default parameters,
require in order to reach convergence on the examples provided with
SuperConga. An F means that it did not converge during the $10^{4}$ iterations
the simulation was run. The convergence criterion is set to $10^{-5}$ for all
examples. In the geometry column, a radius is denoted $\mathcal{R}$, and a
polygon side length $\mathcal{S}$.

Table 1 shows the number of iterations needed to reach convergence in all the
examples, (a)–(g), for each of the different accelerators available in
SuperConga (using default parameters). The examples (a)–(c) are fairly easy as
even Picard iterations converge quickly. It is examples (d)–(g) that really
show the strength of the more advanced methods. In (d) Polyak, BB and CongAcc
provide roughly a speedup of $3\sim 7$ compared to Picard. Both Picard and
Polyak fail to converge in example (f) and (g). This is due to the small
penetration depth in those examples, making the default step sizes too large.
They can be made to converge by manually changing their parameters, however.
The adaptive methods, BB and CongAcc, have no problem with bad initial step
sizes; they will quickly adjust the step size to something appropriate.

Listing 17 shows how to run one of the available examples with a specific
accelerator, where <example> and <accelerator> should be replaced with the
example and accelerator names, respectively.

Listing 17: Running an examples with one of the accelerators.

⬇

python superconga.py simulate -C examples/<example> -A <accelerator>

##  Appendix F Units and basis functions

We summarize in Table 2 our choice of units and the resulting natural scales
for observables and other quantities. The available order parameter basis
functions are listed in Table 3

Input to quasiclassical theory
---
$\mathnormal{v}_{\mathrm{\mbox{$\scriptscriptstyle{\rm{F}}$}}}$ |  | Fermi velocity
$N_{\mathrm{F}}$ |  | normal-state density of states per spin at the Fermi level
$T_{\mathrm{c}}$ |  | superconducting transition temperature (critical temperature)
Geometry
$\mathcal{A}$ |  | total superconducting computational area666We use calligraphic capital letters for quantities related to geometry, such as area ${\cal A}$ or disc radius ${\cal R}$.
$\mathcal{V}$ | $=N_{l}d\mathcal{A}$ | volume of a stack of $N_{l}$ layers with interlayer distance $d$
Natural units
$T_{\mathrm{c}}$ |  |  temperature in units of critical temperature
$2\pi k_{\mathrm{B}}T_{\mathrm{c}}$ |  |  energy in units of critical temperature
$\xi_{0}$ | $=\frac{\hbar\mathnormal{v}_{\mathrm{\mbox{$\scriptscriptstyle{\rm{F}}$}}}}{2\pi k_{\mathrm{B}}T_{\mathrm{c}}}$ |  length in units of superconducting coherence length
$\Phi_{0}$ | $=\frac{hc}{2|e|}$ |  magnetic flux in units of superconducting flux quantum
Derived units and observables
$\lambda_{0}$ | $\frac{1}{\lambda_{0}^{2}}=\frac{4\pi e^{2}}{c^{2}}N_{\mathrm{F}}\mathnormal{v}_{\mathrm{\mbox{$\scriptscriptstyle{\rm{F}}$}}}^{2}$ | penetration depth. Here we use Gaussian units. In SI units $4\pi/c^{2}\rightarrow\mu_{0}$ with $\mu_{0}\approx 4\pi\times 10^{-7}J/A^{2}m$
$\kappa$ | $=\frac{\lambda_{0}}{\xi_{0}}$ | penetration depth in units of coherence length
$\alpha$ | $=\frac{B_{\mathrm{ext}}{\cal A}}{\Phi_{0}}$ | symmetric gauge: external magnetic flux penetrating the superconducting computational area in units of the flux quantum
|  | solenoid gauge: external magnetic flux penetrating a circular hole centered at origo in units of the flux quantum
$A_{0}$ | $=\frac{\Phi_{0}}{\pi\xi_{0}}$ | vector potential
$B_{0}$ | $=\frac{\Phi_{0}}{\pi\xi_{0}^{2}}$ | magnetic flux density
$j_{0}$ | $=2\pi k_{\mathrm{B}}T_{\mathrm{c}}|e|\mathnormal{v}_{\mathrm{\mbox{$\scriptscriptstyle{\rm{F}}$}}}N_{\mathrm{F}}$ | charge-current density
$m_{0}$ | $=j_{0}\xi_{0}\mathcal{V}$ | magnetic moment
$\Omega_{\mathcal{A}}$ | $=(2\pi k_{\mathrm{B}}T_{\mathrm{c}})^{2}N_{\mathrm{F}}{\cal A}$ | free energy
Table 2: Natural units and scales for observables used in this paper and in
the framework SuperConga. We use the convention $e=-|e|$ for the charge of the
electron, but allow for simulations with either positive or negative charge
carriers. Often, we let $\hbar=k_{\mathrm{B}}=c=1$. Order-parameter basis
functions
---
Symmetry class for $D_{4h}$ | Order parameter $\Delta({\bf k})$ | Basis function $\eta(\varphi)$ for circular Fermi surface
$A_{1g}$ ($s$) | $1$ | $1$
$A_{2g}$ ($g_{xy(x^{2}-y^{2})}$) | $k_{x}k_{y}(k_{x}^{2}-k_{y}^{2})$ | $\sqrt{2}\sin(4\varphi)$
$B_{1g}$ ($d_{x^{2}-y^{2}}$) | $k_{x}^{2}-k_{y}^{2}$ | $\sqrt{2}\cos(2\varphi)$
$B_{2g}$ ($d_{xy}$) | $k_{x}k_{y}$ | $\sqrt{2}\sin(2\varphi)$
Table 3: The available symmetries of the order parameter that can be chosen in
SuperConga. Note that any linear combination forming a multi-component order
parameter is also allowed.

## References

  * Tsuei and Kirtley (2000) C. C. Tsuei and J. R. Kirtley, “Pairing symmetry in cuprate superconductors,” Rev. Mod. Phys. 72, 969–1016 (2000).
  * Kirtley  _et al._ (2005) J. R. Kirtley, C. Tsuei, C. Ariando, S. Harkema, and H. Hilgenkamp, “Angle-resolved phase-sensitive determination of the in-plane gap symmetry in yba2cu3o7-$\delta$,” Nature Physics 2, 190 – 194 (2005).
  * Kokubo  _et al._ (2010) N. Kokubo, S. Okayasu, A. Kanda, and B. Shinozaki, “Scanning squid microscope study of vortex polygons and shells in weak-pinning disks of an amorphous superconducting film,” Phys. Rev. B 82, 014501 (2010).
  * Bert  _et al._ (2011) J. A. Bert, B. Kalisky, C. Bell, M. Kim, Y. Hikita, H. Y. Hwang, and K. A. Moler, “Direct imaging of the coexistence of ferromagnetism and superconductivity at the LaAlO3/SrTiO3 interface,” Nature Physics 7, 1 – 5 (2011).
  * Vasyukov  _et al._ (2013) D. Vasyukov, Y. Anahory, L. Embon, D. Halbertal, J. Cuppens, L. Neeman, A. Finkler, Y. Segev, Y. Myasoedov, M. L. Rappaport, M. E. Huber, and E. Zeldov, “A scanning superconducting quantum interference device with single electron spin sensitivity,” Nature Nanotechnology 8, 1 – 6 (2013).
  * Hess  _et al._ (1989) H. Hess, R. Robinson, R. Dynes, J. M. J. Valles, and J. V. Waszczak, “Scanning-Tunneling-Microscope observation of the Abrikosov flux lattice and the density of states near and inside a fluxoid.”  (1989).
  * Renner  _et al._ (1991) C. Renner, A. D. Kent, P. Niedermann, Ø. Fischer, and F. Lévy, “Scanning tunneling spectroscopy of a vortex core from the clean to the dirty limit,” Physical Review Letters 67, 1650–1652 (1991).
  * Maggio-Aprile  _et al._ (1995) I. Maggio-Aprile, C. Renner, A. Erb, E. Walker, and O. Fischer, “Direct Vortex Lattice Imaging and Tunneling Spectroscopy of Flux Lines on YBa2Cu3O7-δ,” Phys. Rev. Lett. 75, 2754–2757 (1995).
  * Yazdani  _et al._ (1997) A. Yazdani, B. Jones, C. Lutz, M. F. Crommie, and D. M. Eigler, “Probing the Local Effects of Magnetic Impurities on Superconductivity,” Science 275, 1767 – 1770 (1997).
  * Pan  _et al._ (2000a) S. Pan, E. Hudson, A. Gupta, K. Ng, H. Eisaki, S. Uchida, and J. Davis, “STM Studies of the Electronic Structure of Vortex Cores in Bi2Sr2CaCu2O8+δ,” Phys. Rev. Lett. 85, 1536 – 1539 (2000a).
  * Pan  _et al._ (2000b) S. Pan, E. Hudson, K. Lang, H. Eisaki, S. Uchida, and J. Davis, “Imaging the effects of individual zinc impurity atoms on superconductivity in Bi2Sr2CaCu2O8+δ,” Nature 403, 746 – 750 (2000b).
  * Hoffman  _et al._ (2002) J. E. Hoffman, K. McElroy, D.-H. Lee, K. M. Lang, H. Eisaki, S. Uchida, and J. C. Davis, “Imaging Quasiparticle Interference in Bi2Sr2CaCu2O8+δ,” Science 297, 1148–1151 (2002), cond-mat/0209276 .
  * Guillamón  _et al._ (2008) I. Guillamón, H. Suderow, S. Vieira, L. Cario, P. Diener, and P. Rodière, “Superconducting Density of States and Vortex Cores of 2H-NbS2,” Physical Review Letters 101, 166407 (2008), 0807.1809 .
  * Roditchev  _et al._ (2015) D. Roditchev, C. Brun, L. Serrier-Garcia, J. C. Cuevas, V. H. L. Bessa, M. V. Milošević, F. Debontridder, V. Stolyarov, and T. Cren, “Direct observation of Josephson vortex cores,” Nature Physics 11, 332 – 337 (2015).
  * Berthod  _et al._ (2017) C. Berthod, I. Maggio-Aprile, J. Bruér, A. Erb, and C. Renner, “Observation of Caroli–de Gennes–Matricon Vortex States in ${\mathrm{YBa}}_{2}{\mathrm{Cu}}_{3}{\mathrm{O}}_{7-\delta}$,” Phys. Rev. Lett. 119, 237001 (2017).
  * Bolle  _et al._ (1999) C. A. Bolle, V. Aksyuk, F. Pardo, P. L. Gammel, E. Zeldov, E. Bucher, R. Boie, D. J. Bishop, and D. R. Nelson, “Observation of mesoscopic vortex physics using micromechanical oscillators,” Nature 399, 43–46 (1999).
  * Bleszynski-Jayich  _et al._ (2009) A. C. Bleszynski-Jayich, W. E. Shanks, B. Peaudecerf, E. Ginossar, F. v. Oppen, L. Glazman, and J. G. E. Harris, “Persistent Currents in Normal Metal Rings,” Science 326, 272–275 (2009).
  * Jang  _et al._ (2011) J. Jang, D. G. Ferguson, V. Vakaryuk, R. Budakian, S. B. Chung, P. M. Goldbart, and Y. Maeno, “Observation of Half-Height Magnetization Steps in Sr2RuO4,” Science 331, 186 – 188 (2011).
  * Geim  _et al._ (1997) A. K. Geim, I. V. Grigorieva, S. V. Dubonos, J. G. S. Lok, J. C. Maan, A. E. Filippov, and F. M. Peeters, “Phase transitions in individual sub-micrometre superconductors,” Nature 390, 259–262 (1997).
  * Morelle, Bekaert, and Moshchalkov (2004) M. Morelle, J. Bekaert, and V. V. Moshchalkov, “Influence of sample geometry on vortex matter in superconducting microstructures,” Phys. Rev. B 70, 094503 (2004).
  * Grigorieva  _et al._ (2006) I. V. Grigorieva, W. Escoffier, J. Richardson, L. Y. Vinnikov, S. Dubonos, and V. Oboznov, “Direct observation of vortex shells and magic numbers in mesoscopic superconducting disks,” Phys. Rev. Lett. 96, 077005 (2006).
  * Khotkevych, Milošević, and Bending (2008) V. V. Khotkevych, M. V. Milošević, and S. J. Bending, “A scanning hall probe microscope for high resolution magnetic imaging down to 300 mk,” Review of Scientific Instruments 79, 123708 (2008), https://doi.org/10.1063/1.3046285 .
  * Curran  _et al._ (2014) P. J. Curran, S. J. Bending, W. M. Desoky, A. S. Gibbs, S. L. Lee, and A. P. Mackenzie, “Search for spontaneous edge currents and vortex imaging in ${\mathrm{sr}}_{2}\mathrm{Ru}{\mathrm{o}}_{4}$ mesostructures,” Phys. Rev. B 89, 144504 (2014).
  * Ge  _et al._ (2017) J.-Y. Ge, V. N. Gladilin, J. Tempere, V. S. Zharinov, J. V. d. Vondel, J. T. Devreese, and V. V. Moshchalkov, “Direct visualization of vortex ice in a nanostructured superconductor,” Phys Rev B 96, 134515 – 7 (2017).
  * Tuominen  _et al._ (1992) M. T. Tuominen, J. M. Hergenrother, T. S. Tighe, and M. Tinkham, “Experimental evidence for parity-based 2e periodicity in a superconducting single-electron tunneling transistor,”  69, 1997 – 2000 (1992).
  * Lafarge  _et al._ (1993) P. Lafarge, P. Joyez, D. Esteve, C. Urbina, and M. H. Devoret, “Measurement of the even-odd free-energy difference of an isolated superconductor,” Phys. Rev. Lett. 70, 994–997 (1993).
  * Ralph, Black, and Tinkham (1995) D. C. Ralph, C. T. Black, and M. Tinkham, “Spectroscopic Measurements of Discrete Electronic States in Single Metal Particles,” Physical Review Letters 74, 3241–3244 (1995).
  * Kubatkin  _et al._ (1996) S. E. Kubatkin, A. Y. Tzalenchuk, Z. G. Ivanov, P. Delsing, R. I. Shekhter, and T. Claeson, “Coulomb blockade electrometer with a high-Tc island,” JETP Letters 63, 126–132 (1996).
  * Gustafsson  _et al._ (2013) D. Gustafsson, D. Golubev, M. Fogelström, T. Claeson, S. Kubatkin, T. Bauch, and F. Lombardi, “Fully gapped superconductivity in a nanometre-size YBa2Cu3O7-δ island enhanced by a magnetic field,” Nature Nanotechnology 8, 25 – 30 (2013).
  * (30) J. Bardeen, L. N. Cooper, and J. R. Schrieffer, “Theory of Superconductivity,” Phys. Rev. 108, 1175–1204.
  * Schrieffer (1999) J. R. Schrieffer,  _Theory of Superconductivity_ , Advanced Book Classics (Perseus Books, Reading, Massachusetts, USA, 1999).
  * de Gennes (1989) P. G. de Gennes,  _Superconductivity of Metals and Alloys_ , Advanced Book Classics (Addison-Wesley Publishing Company, Reading, Massachusetts, USA, 1989).
  * Tinkham (2004) M. Tinkham,  _Introduction to Superconductivity: Second Edition_, Dover Books on Physics (Dover Publications, 2004).
  * Mineev and Samokhin (1999) V. P. Mineev and K. V. Samokhin,  _Introduction to Unconventional Superconductivity_ (Gordon and Breach Science Publishers, Amsterdam, 1999).
  * Serene and Rainer (1983) J. Serene and D. Rainer, “The quasiclassical approach to superfluid 3He,” Physics Reports 101, 221 – 311 (1983).
  * Eilenberger (1968) G. Eilenberger, “Transformation of Gorkov’s equation for type II superconductors into transport-like equations,” Zeitschrift für Physik A Hadrons and Nuclei 214, 195–213 (1968).
  * Larkin and Ovchinnikov (1968) A. I. Larkin and Y. N. Ovchinnikov, “Quasiclassical method in the theory of superconductivity,” Zh. Eksp. Teor. Fiz. 55, 2262 (1968), [Sov. Phys. JETP 28, 1200 (1969)].
  * Eliashberg (1971) G. M. Eliashberg, “Inelastic electron collisions and nonequilibrium stationary states in superconductors,” Zh. Eksp. Teor. Fiz. 61, 1254 (1971), [Sov. Phys. JETP 34, 668 (1972)].
  * Abrikosov, Gorkov, and Dzyaloshinski (1975) A. A. Abrikosov, L. P. Gorkov, and I. E. Dzyaloshinski,  _Methods of quantum field theory in statistical physics_ (Dover Publications, Inc., New York, 1975).
  * (40) D. Rainer and J. A. Sauls, “Strong-coupling theory of superconductivity,” in  _Superconductivity: From Basic Physics to the Latest Developments_ (WORLD SCIENTIFIC) pp. 45–78.
  * Kopnin (2001) N. Kopnin,  _Theory of Nonequilibrium Superconductivity_ (Clarendon Press, Oxford, 2001).
  * Seja and Lofwander (2022) K. M. Seja and T. Lofwander, “A finite element method for the quasiclassical theory of superconductivity,”  (2022).
  * Andreev (1964) A. A. Andreev, “The thermal conductivity of the intermediate state in superconductors,” Zh. Eksp. Teor. Fiz. 46, 1823 (1964), [Engl. Transl. 1964 Sov. Phys. JETP 19, 1228].
  * (44) A. Shelankov and M. Ozana, “Quasiclassical theory of superconductivity: A multiple-interface geometry,” Phys. Rev. B 61, 7077–7100.
  * Usadel (1970) K. D. Usadel, “Generalized Diffusion Equation for Superconducting Alloys,” Phys. Rev. Lett. 25, 507–509 (1970).
  * Alexander  _et al._ (1985) J. A. X. Alexander, T. P. Orlando, D. Rainer, and P. M. Tedrow, “Theory of fermi-liquid effects in high-field tunneling,” Phys. Rev. B 31, 5811–5825 (1985).
  * Belzig  _et al._ (1999) W. Belzig, F. K. Wilhelm, C. Bruder, G. Schön, and A. D. Zaikin, “Quasiclassical green’s function approach to mesoscopic superconductivity,” Superlattices and Microstructures 25, 1251–1288 (1999).
  * Bergeret  _et al._ (2018) F. S. Bergeret, M. Silaev, P. Virtanen, and T. T. Heikkilä, “Colloquium: Nonequilibrium effects in superconductors with a spin-splitting field,” Rev. Mod. Phys. 90, 041001 (2018).
  * Sauls (2022) J. A. Sauls, “Theory of disordered superconductors with applications to nonlinear current response,” Progress of Theoretical and Experimental Physics 2022 (2022), 10.1093/ptep/ptac034, 033I03, https://academic.oup.com/ptep/article-pdf/2022/3/033I03/43152929/ptac034.pdf .
  * Nazarov (1999) Y. V. Nazarov, “Novel circuit theory of andreev reflection,” Superlattices and Microstructures 25, 1221–1231 (1999).
  * Amundsen and Linder (2016) M. Amundsen and J. Linder, “General solution of 2d and 3d superconducting quasiclassical systems: coalescing vortices and nanoisland geometries,” Scientific Reports 6, 22765 (2016).
  * Timmermans  _et al._ (2016) M. Timmermans, L. Serrier-Garcia, M. Perini, J. Van de Vondel, and V. V. Moshchalkov, “Direct observation of condensate and vortex confinement in nanostructured superconductors,” Phys. Rev. B 93, 054514 (2016).
  * Vorontsov and Sauls (2003) A. B. Vorontsov and J. Sauls, “Thermodynamic properties of thin films of superfluid 3He-A,” 68, 064508 (2003).
  * Håkansson, Löfwander, and Fogelström (2015) M. Håkansson, T. Löfwander, and M. Fogelström, “Spontaneously broken time-reversal symmetry in high-temperature superconductors,” Nature Physics 11, 755 (2015).
  * Holmvall  _et al._ (2020) P. Holmvall, M. Fogelström, T. Löfwander, and A. B. Vorontsov, “Phase crystals,” Phys. Rev. Research 2, 013104 (2020).
  * Yip and Garg (1993) S. Yip and A. Garg, “Superconducting states of reduced symmetry: General order parameters and physical implications,” Phys. Rev. B 48, 3304–3308 (1993).
  * Graf, Rainer, and Sauls (1993) M. Graf, D. Rainer, and J. Sauls, “Coupled two-dimensional Fermi liquids as a model for layered superconductors: Basic equations and elementary results.” Physical Review B (Condensed Matter) 47, 12089–12098 (1993).
  * Wennerdal  _et al._ (2020) N. W. Wennerdal, A. Ask, P. Holmvall, T. Löfwander, and M. Fogelström, “Breaking time-reversal and translational symmetry at edges of $d$-wave superconductors: Microscopic theory and comparison with quasiclassical theory,” Phys. Rev. Research 2, 043198 (2020).
  * Eschrig (2000) M. Eschrig, “Distribution functions in nonequilibrium theory of superconductivity and andreev spectroscopy in unconventional superconductors,” Phys. Rev. B 61, 9061–9076 (2000).
  * Luttinger and Ward (1960) J. M. Luttinger and J. C. Ward, “Ground-State Energy of a Many-Fermion System. II,” Physical Review 118, 1417–1427 (1960).
  * Thuneberg, Kurkijärvi, and Rainer (1984) E. V. Thuneberg, J. Kurkijärvi, and D. Rainer, “Elementary-flux-pinning potential in type-ii superconductors,” Phys. Rev. B 29, 3913–3923 (1984).
  * Zaitsev (1984) A. V. Zaitsev, “Quasiclassical equation of the theory of superconductivity for contiguous metals and the properties of constricted microstructures,” Zh. Eksp. Teor. Fiz. 86, 1742 (1984).
  * (63) G. Kieselmann, “Self-consistent calculations of the pair potential and the tunneling density of states in proximity contacts,” Phys. Rev. B 35, 6762–6770.
  * (64) A. Millis, D. Rainer, and J. A. Sauls, “Quasiclassical theory of superconductivity near magnetically active interfaces,” Phys. Rev. B 38, 4504–4515.
  * (65) M. Fogelström, “Josephson currents through spin-active interfaces,”  62, 11812–11819.
  * (66) E. Zhao, T. Löfwander, and J. A. Sauls, “Nonequilibrium superconductivity near spin-active interfaces,” Phys. Rev. B 70, 134510.
  * Eschrig (2009) M. Eschrig, “Scattering problem in nonequilibrium quasiclassical theory of metals and superconductors: General boundary conditions and applications,” Phys Rev B 80, 134511 (2009).
  * Nagato, Nagai, and Hara (1993) Y. Nagato, K. Nagai, and J. Hara, “Theory of the Andreev reflection and the density of states in proximity contact normal-superconducting infinite double-layer,” Journal of Low Temperature Physics 93, 33–56 (1993).
  * Schopohl and Maki (1995) N. Schopohl and K. Maki, “Quasiparticle spectrum around a vortex line in a d-wave superconductor,” Phys. Rev. B 52, 490–493 (1995).
  * (70) N. Schopohl, “Transformation of the eilenberger equations of superconductivity to a scalar riccati equation,”  cond-mat/9804064 .
  * Ozaki (2007) T. Ozaki, “Continued fraction representation of the fermi-dirac function for large-scale electronic structure calculations,” Phys. Rev. B 75, 035123 (2007).
  * NVIDIA Corporation (2021a) NVIDIA Corporation, “NVIDIA CUDA C programming guide,”  (2021a), version 11.2.0.
  * Nickolls  _et al._ (2008) J. Nickolls, I. Buck, M. Garland, and K. Skadron, “Scalable parallel programming with cuda,” Queue 6, 40–53 (2008).
  * SuperConga team (2022a) SuperConga team, “SuperConga Online Documentation and User Manual,” https://superconga.gitlab.io/superconga-doc/ (2022a).
  * Polyak (1964) B. Polyak, “Some methods of speeding up the convergence of iteration methods,” USSR Computational Mathematics and Mathematical Physics 4, 1–17 (1964).
  * Barzilai and Borrwein (1988) J. Barzilai and J. M. Borrwein, “Two-Point Step Size Gradient Methods,” IMA Journal of Numerical Analysis 8, 141–148 (1988), https://academic.oup.com/imajna/article-pdf/8/1/141/2402762/8-1-141.pdf .
  * SuperConga team (2022b) SuperConga team, “SuperConga GitLab repository,” https://gitlab.com/superconga/superconga (2022b).
  * Holmvall  _et al._ (2019) P. Holmvall, A. B. Vorontsov, M. Fogelström, and T. Löfwander, “Spontaneous symmetry breaking at surfaces of $d$-wave superconductors: Influence of geometry and surface ruggedness,” Phys. Rev. B 99, 184511 (2019).
  * Caroli, De Gennes, and Matricon (1964) C. Caroli, P. De Gennes, and J. Matricon, “Bound fermion states on a vortex line in a type ii superconductor,” Physics Letters 9, 307–309 (1964).
  * Kim  _et al._ (2021) H. Kim, Y. Nagai, L. Rózsa, D. Schreyer, and R. Wiesendanger, “Anisotropic non-split zero-energy vortex bound states in a conventional superconductor,” (2021), arXiv:2105.01354 [cond-mat.supr-con] .
  * Bean and Livingston (1964) C. P. Bean and J. D. Livingston, “Surface barrier in type-ii superconductors,” Phys. Rev. Lett. 12, 14–16 (1964).
  * Zhao  _et al._ (2008a) H. J. Zhao, V. R. Misko, F. M. Peeters, V. Oboznov, S. V. Dubonos, and I. V. Grigorieva, “Vortex states in mesoscopic superconducting squares: Formation of vortex shells,” Phys. Rev. B 78, 104517 (2008a).
  * Misko  _et al._ (2009) V. R. Misko, H. J. Zhao, F. M. Peeters, V. Oboznov, S. V. Dubonos, and I. V. Grigorieva, “Formation of vortex shells in mesoscopic superconducting squares,” Superconductor Science and Technology 22, 034001 (2009).
  * Zhang  _et al._ (2012) L.-F. Zhang, L. Covaci, M. V. Milošević, G. R. Berdiyorov, and F. M. Peeters, “Unconventional vortex states in nanoscale superconductors due to shape-induced resonances in the inhomogeneous cooper-pair condensate,” Phys. Rev. Lett. 109, 107001 (2012).
  * Zhang  _et al._ (2013) L.-F. Zhang, L. Covaci, M. V. Milošević, G. R. Berdiyorov, and F. M. Peeters, “Vortex states in nanoscale superconducting squares: The influence of quantum confinement,” Phys. Rev. B 88, 144501 (2013).
  * Zhao  _et al._ (2008b) H. J. Zhao, V. R. Misko, F. M. Peeters, S. Dubonos, V. Oboznov, and I. V. Grigorieva, “Vortex configurations in mesoscopic superconducting triangles: Finite-size and shape effects,” Europhys. Lett. 83, 17008 (2008b).
  * Cabral and Aguiar (2009) L. R. E. Cabral and J. A. Aguiar, “Vortex configurations in thin superconducting equilateral triangles,” Phys. Rev. B 80, 214533 (2009).
  * Kokubo  _et al._ (2016) N. Kokubo, H. Miyahara, S. Okayasu, and T. Nojima, “Vortex shells in mesoscopic triangles of amorphous superconducting thin films,” Physica C: Superconductivity and its Applications 530, 42–45 (2016), 28th International Symposium on Superconductivity.
  * Wu  _et al._ (2017) W. Wu, H. J. Zhao, B. Li, and Z. X. Shi, “Finite-size and shape effects on vortex configurations in superconducting isosceles right triangles,” Supercond. Sci. Technol. 30, 125004 (2017).
  * Huy, Kato, and Ishida (2013) H. T. Huy, M. Kato, and T. Ishida, “Vortex states in de factomesoscopic mo80ge20pentagon plates,” Superconductor Science and Technology 26, 065001 (2013).
  * Fink and Presson (1966) H. J. Fink and A. G. Presson, “Magnetic irreversible solution of the ginzburg-landau equations,” Phys. Rev. 151, 219–228 (1966).
  * Moshchalkov  _et al._ (1996) V. V. Moshchalkov, M. Baert, V. V. Metlushko, E. Rosseel, M. J. Van Bael, K. Temst, R. Jonckheere, and Y. Bruynseraede, “Magnetization of multiple-quanta vortex lattices,” Phys. Rev. B 54, 7385–7393 (1996).
  * Schweigert, Peeters, and Deo (1998) V. A. Schweigert, F. M. Peeters, and P. S. Deo, “Vortex phase diagram for mesoscopic superconducting disks,” Phys. Rev. Lett. 81, 2783–2786 (1998).
  * Kanda  _et al._ (2004) A. Kanda, B. J. Baelus, F. M. Peeters, K. Kadowaki, and Y. Ootuka, “Experimental evidence for giant vortex states in a mesoscopic superconducting disk,” Phys. Rev. Lett. 93, 257002 (2004).
  * Golubović  _et al._ (2005) D. S. Golubović, M. V. Milošević, F. M. Peeters, and V. V. Moshchalkov, “Magnetically induced splitting of a giant vortex state in a mesoscopic superconducting disk,” Phys. Rev. B 71, 180502 (2005).
  * Cren  _et al._ (2011) T. Cren, L. Serrier-Garcia, F. Debontridder, and D. Roditchev, “Vortex fusion and giant vortex states in confined superconducting condensates,” Phys. Rev. Lett. 107, 097202 (2011).
  * Chibotaru  _et al._ (2000) L. F. Chibotaru, A. Ceulemans, V. Bruyndoncx, and V. V. Moshchalkov, “Symmetry-induced formation of antivortices in mesoscopic superconductors,” Nature 408, 833–835 (2000).
  * Iavarone  _et al._ (2011) M. Iavarone, A. Scarfato, F. Bobba, M. Longobardi, G. Karapetrov, V. Novosad, V. Yefremenko, F. Giubileo, and A. M. Cucolo, “Imaging the spontaneous formation of vortex-antivortex pairs in planar superconductor/ferromagnet hybrid structures,” Phys. Rev. B 84, 024506 (2011).
  * Giorgio  _et al._ (2017) C. D. Giorgio, D. D’Agostino, A. M. Cucolo, M. Iavarone, A. Scarfato, G. Karapetrov, S. A. Moore, M. Polichetti, D. Mancusi, S. Pace, V. Novosad, V. Yefremenko, and F. Bobba, “Superconducting vortex-antivortex pairs: Nucleation and confinement in magnetically coupled superconductor ferromagnet hybrids,” in  _Vortex Dynamics and Optical Vortices_, edited by H. Perez-de Tejada (InTech, Rijeka, 2017) Chap. 03.
  * Simmendinger  _et al._ (2020) J. Simmendinger, M. Bihler, A. M. Ionescu, M. Weigand, G. Schütz, and J. Albrecht, “Bound and stable vortex–antivortex pairs in high-t c superconductors,” New Journal of Physics 22, 123035 (2020).
  * Vorontsov (2009) A. B. Vorontsov, “Broken translational and time-reversal symmetry in unconventional superconducting films,” Phys. Rev. Lett. 102, 177001 (2009).
  * Holmvall (2017) P. Holmvall,  _Modeling mesoscopic unconventional superconductors_ (Chalmers University of Technology, 2017) licentiate thesis.
  * Holmvall, Löfwander, and Fogelström (2018) P. Holmvall, T. Löfwander, and M. Fogelström, “Spontaneous generation of fractional vortex-antivortex pairs at single edges of high-Tc superconductors,” Journal of Physics: Conference Series 969, 012037 (2018).
  * Holmvall  _et al._ (2018) P. Holmvall, A. B. I. Vorontsov, M. Fogelström, and T. Löfwander, “Broken translational symmetry at edges of high-temperature superconductors,” Nature Communications 9, 2190 (2018).
  * Holmvall (2019) P. Holmvall,  _Crystallization of the superconducting phase in unconventional superconductors_ (Chalmers University of Technology, 2019) phD thesis.
  * Abrikosov (1957) A. A. Abrikosov, “On the magnetic properties of superconductors of the second group,” Soviet Phys JETP 5, 1174–1182 (1957).
  * Pearl (1964) J. Pearl, “Current distribution in superconducting films carrying quantized fluxoids,” Applied Physics Letters 5, 65–66 (1964), https://doi.org/10.1063/1.1754056 .
  * Xu, Yip, and Sauls (1995) D. Xu, S. K. Yip, and J. A. Sauls, “Nonlinear meissner effect in unconventional superconductors,” Phys. Rev. B 51, 16233–16253 (1995).
  * Sakurai and Napolitano (2011) J. Sakurai and J. Napolitano,  _Modern Quantum Mechanics_ (Addison-Wesley, 2011).
  * Aharonov and Bohm (1959) Y. Aharonov and D. Bohm, “Significance of electromagnetic potentials in the quantum theory,” Phys. Rev. 115, 485–491 (1959).
  * Rybakov and Babaev (2022) F. N. Rybakov and E. Babaev, “Flux qubit based on hybrid ferromagnetic-superconducting device,”  (2022).
  * Deaver and Fairbank (1961) B. S. Deaver and W. M. Fairbank, “Experimental evidence for quantized flux in superconducting cylinders,” Phys. Rev. Lett. 7, 43–46 (1961).
  * Doll and Näbauer (1961) R. Doll and M. Näbauer, “Experimental proof of magnetic flux quantization in a superconducting ring,” Phys. Rev. Lett. 7, 51–52 (1961).
  * Onsager (1961) L. Onsager, “Magnetic flux through a superconducting ring,” Phys. Rev. Lett. 7, 50–50 (1961).
  * Byers and Yang (1961) N. Byers and C. N. Yang, “Theoretical considerations concerning quantized magnetic flux in superconducting cylinders,” Phys. Rev. Lett. 7, 46–49 (1961).
  * Brenig (1961) W. Brenig, “Remark concerning quantized magnetic flux in superconductors,” Phys. Rev. Lett. 7, 337–339 (1961).
  * Little and Parks (1962) W. A. Little and R. D. Parks, “Observation of quantum periodicity in the transition temperature of a superconducting cylinder,” Phys. Rev. Lett. 9, 9–12 (1962).
  * Groff and Parks (1968) R. P. Groff and R. D. Parks, “Fluxoid quantization and field-induced depairing in a hollow superconducting microcylinder,” Phys. Rev. 176, 567–580 (1968).
  * Buisson  _et al._ (1990) O. Buisson, P. Gandit, R. Rammal, Y. Wang, and B. Pannetier, “Magnetization oscillations of a superconducting disk,” Physics Letters A 150, 36–42 (1990).
  * Hu (1994) C.-R. Hu, “Midgap surface states as a novel signature for ${\mathit{d}}_{\mathit{x}\mathit{a}}^{2}$-${\mathit{x}}_{\mathit{b}}^{2}$-wave superconductivity,” Phys. Rev. Lett. 72, 1526–1529 (1994).
  * Suzuki and Asano (2014) S.-I. Suzuki and Y. Asano, “Paramagnetic instability of small topological superconductors,” Phys. Rev. B 89, 184508 (2014).
  * Fetter (1980) A. L. Fetter, “Flux-Penetration in a Thin Superconducting Disk,” Physical Review B (Condensed Matter) 22, 1200–1213 (1980).
  * Kogan, Clem, and Mints (2004) V. G. Kogan, J. R. Clem, and R. G. Mints, “Properties of mesoscopic superconducting thin-film rings: London approach,” Physical Review B 69, 064516 (2004), cond-mat/0310393 .
  * Higashitani (1997) S. Higashitani, “Mechanism of paramagnetic meissner effect in high-temperature superconductors,” J. Soc. Phys. Jpn. 66, 2556 (1997).
  * Fogelström, Rainer, and Sauls (1997) M. Fogelström, D. Rainer, and J. A. Sauls, “Tunneling into current-carrying surface states of high- ${\mathit{t}}_{\mathit{c}}$ superconductors,” Phys. Rev. Lett. 79, 281–284 (1997).
  * Fogelström, Rainer, and Sauls (2004) M. Fogelström, D. Rainer, and J. Sauls, “Nonlinear magnetic field dependence of the conductance in d-wave normal-metal/insulator/superconductor tunnel junctions,” Phys Rev B 70, 012503 (2004).
  * Clem, Coffey, and Hao (1991) J. R. Clem, M. W. Coffey, and Z. Hao, “Lower critical field of a josephson-coupled layer model of high-${\mathit{t}}_{\mathit{c}}$ superconductors,” Phys. Rev. B 44, 2732–2738 (1991).
  * NVIDIA Corporation (2021b) NVIDIA Corporation, “NVIDIA CUDA enabled devices,”  (2021b).
  * Horn and Moler (2022) L. B.-V. Horn and K. A. Moler, “Superscreen: An open-source package for simulating the magnetic response of two-dimensional superconducting devices,”  (2022).
  * Sadovskyy  _et al._ (2015) I. Sadovskyy, A. Koshelev, C. Phillips, D. Karpeyev, and A. Glatz, “Stable large-scale solver for ginzburg–landau equations for superconductors,” Journal of Computational Physics 294, 639–654 (2015).
  * (131) F. N. Rybakov and E. Babaev, “Excalibur software,”
URL: http://quantumandclassical.com/excalibur/.

  * Kresse and Furthmüller (1996) G. Kresse and J. Furthmüller, “Efficiency of ab-initio total energy calculations for metals and semiconductors using a plane-wave basis set,” Computational Materials Science 6, 15–50 (1996).
  * Kresse and Furthmüller (1996) G. Kresse and J. Furthmüller, “Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set,” Phys. Rev. B 54, 11169–11186 (1996).
  * Giannozzi  _et al._ (2009) P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, G. L. Chiarotti, M. Cococcioni, I. Dabo, A. D. Corso, S. de Gironcoli, S. Fabris, G. Fratesi, R. Gebauer, U. Gerstmann, C. Gougoussis, A. Kokalj, M. Lazzeri, L. Martin-Samos, N. Marzari, F. Mauri, R. Mazzarello, S. Paolini, A. Pasquarello, L. Paulatto, C. Sbraccia, S. Scandolo, G. Sclauzero, A. P. Seitsonen, A. Smogunov, P. Umari, and R. M. Wentzcovitch, “QUANTUM ESPRESSO: a modular and open-source software project for quantum simulations of materials,” Journal of Physics: Condensed Matter 21, 395502 (2009).
  * Giannozzi  _et al._ (2017) P. Giannozzi, O. Andreussi, T. Brumme, O. Bunau, M. B. Nardelli, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, M. Cococcioni, N. Colonna, I. Carnimeo, A. D. Corso, S. de Gironcoli, P. Delugas, R. A. DiStasio, A. Ferretti, A. Floris, G. Fratesi, G. Fugallo, R. Gebauer, U. Gerstmann, F. Giustino, T. Gorni, J. Jia, M. Kawamura, H.-Y. Ko, A. Kokalj, E. Küçükbenli, M. Lazzeri, M. Marsili, N. Marzari, F. Mauri, N. L. Nguyen, H.-V. Nguyen, A. O. de-la Roza, L. Paulatto, S. Poncé, D. Rocca, R. Sabatini, B. Santra, M. Schlipf, A. P. Seitsonen, A. Smogunov, I. Timrov, T. Thonhauser, P. Umari, N. Vast, X. Wu, and S. Baroni, “Advanced capabilities for materials modelling with quantum ESPRESSO,” Journal of Physics: Condensed Matter 29, 465901 (2017).
  * Nakhaee, Ketabi, and Peeters (2020) M. Nakhaee, S. A. Ketabi, and F. M. Peeters, “Tight-binding studio: A technical software package to find the parameters of tight-binding hamiltonian,” Computer Physics Communications 254, 107379 (2020).
  * (137) C. W. Groth, M. Wimmer, A. R. Akhmerov, and X. Waintal, “Kwant: a software package for quantum transport,”  16, 063065, publisher: IOP Publishing.
  * Björnson (2019) K. Björnson, “Tbtk: A quantum mechanics software development kit,” SoftwareX 9, 205–210 (2019).
  * Löthman  _et al._ (2021) T. Löthman, C. Triola, J. Cayao, and A. M. Black-Schaffer, “Efficient numerical method for evaluating normal and anomalous time-domain equilibrium green’s functions in inhomogeneous systems,” Phys. Rev. B 104, 125405 (2021).
  * Nagai (2020) Y. Nagai, “N-independent localized krylov–bogoliubov-de gennes method: Ultra-fast numerical approach to large-scale inhomogeneous superconductors,” Journal of the Physical Society of Japan 89, 074703 (2020), https://doi.org/10.7566/JPSJ.89.074703 .
  * Bell and Hoberock (2012) N. Bell and J. Hoberock, “Chapter 26 - thrust: A productivity-oriented library for cuda,” in  _GPU Computing Gems Jade Edition_, Applications of GPU Computing Series, edited by W. mei W. Hwu (Morgan Kaufmann, Boston, 2012) pp. 359 – 371.
  * Kaczmarski and Rzkażewski (2013) K. Kaczmarski and P. Rzkażewski, “Thrust and cuda in data intensive algorithms,” in  _New Trends in Databases and Information Systems_ , edited by M. Pechenizkiy and M. Wojciechowski (Springer Berlin Heidelberg, Berlin, Heidelberg, 2013) pp. 37–46.
  * Wynters (2013) E. Wynters, “Examples from computational geometry that demonstrate the potential of using the thrust library to implement parallel processing on gpus,” J. Comput. Sci. Coll. 28, 148–155 (2013).
  * Cedilnik, E. and Hoffman, B. and King, B. and Martin, K. and Neundorf, A. (2021) Cedilnik, E. and Hoffman, B. and King, B. and Martin, K. and Neundorf, A., “CMake,” https://cmake.org/ (2021), version 3.21.3.
  * Martin, E. (2020) Martin, E., “Ninja,” https://ninja-build.org/ (2020), version 1.10.2.
  * Sanderson and Curtin (2016) C. Sanderson and R. Curtin, “Armadillo: a template-based C++ library for linear algebra,” The Journal of Open Source Software 1, 26 (2016).
  * Sanderson and Curtin (2018) C. Sanderson and R. Curtin, “A user-friendly hybrid sparse matrix class in c++,” in  _Mathematical Software – ICMS 2018_ , edited by J. H. Davenport, M. Kauers, G. Labahn, and J. Urban (Springer International Publishing, Cham, 2018) pp. 422–430.
  * Kirilov (2021) V. Kirilov, “doctest,” https://github.com/onqtam/doctest (2021).
  * Khronos Group (2021) Khronos Group, “Opengl overview,”  (2021), version 4.6.
  * Woo  _et al._ (1999) M. Woo, J. Neider, T. Davis, and D. Shreiner,  _OpenGL programming guide: the official guide to learning OpenGL, version 1.2_ (Addison-Wesley Longman Publishing Co., Inc., 1999).
  * Yalamanchili  _et al._ (2015) P. Yalamanchili, U. Arshad, Z. Mohammed, P. Garigipati, P. Entschev, B. Kloppenborg, J. Malcolm, and J. Melonakos, “ArrayFire - A high performance software library for parallel computing with an easy-to-use API,”  (2015).
  * Dunn (2020) C. Dunn, “Jsoncpp,” https://github.com/open-source-parsers/jsoncpp (2020).
  * Pezoa  _et al._ (2016) F. Pezoa, J. L. Reutter, F. Suarez, M. Ugarte, and D. Vrgoč, “Foundations of json schema,” in  _Proceedings of the 25th International Conference on World Wide Web_ (International World Wide Web Conferences Steering Committee, 2016) pp. 263–273.
  * BlueBrain (2021) BlueBrain, “HighFive - Header-only C++ HDF5 interface,”  (2021), version 2.3.1.
  * The HDF Group (NNNN) The HDF Group, “Hierarchical Data Format, version 5,” (1997-NNNN), https://www.hdfgroup.org/HDF5/.
  * Harris  _et al._ (2020) C. R. Harris, K. J. Millman, S. J. van der Walt, R. Gommers, P. Virtanen, D. Cournapeau, E. Wieser, J. Taylor, S. Berg, N. J. Smith, R. Kern, M. Picus, S. Hoyer, M. H. van Kerkwijk, M. Brett, A. Haldane, J. F. del Río, M. Wiebe, P. Peterson, P. Gérard-Marchant, K. Sheppard, T. Reddy, W. Weckesser, H. Abbasi, C. Gohlke, and T. E. Oliphant, “Array programming with NumPy,” Nature 585, 357–362 (2020).
  * Hunter (2007) J. D. Hunter, “Matplotlib: A 2d graphics environment,” Computing in Science & Engineering 9, 90–95 (2007).
  * Collette (2013) A. Collette,  _Python and HDF5_ (O’Reilly, 2013).
  * McKinney  _et al._ (2010) W. McKinney _et al._ , “Data structures for statistical computing in python,” in  _Proceedings of the 9th Python in Science Conference_ , Vol. 445 (Austin, TX, 2010) pp. 51–56.
  * MATLAB (2021) MATLAB,  _9.11.0.1769968 (R2021b)_ (The MathWorks Inc., Natick, Massachusetts, 2021).
  * Grein, Löfwander, and Eschrig (2013) R. Grein, T. Löfwander, and M. Eschrig, “Inverse proximity effect and influence of disorder on triplet supercurrents in strongly spin-polarized ferromagnets,” Phys. Rev. B 88, 054502 (2013).