results.tex

%results

\acresetall

\part{Results}
\label{sec:results}

\chapter[Golden Dual Fullerenes]{Golden Dual Fullerenes\footnote{This chapter is
    composed of sections previously published in the article
    \citetitle*{Trombach_HollowGoldCages_2016}\autocite{Trombach_HollowGoldCages_2016}
    and is reprinted by permission from the publisher \textcopyright 2016
    John Wiley and Sons. Some sections have been modified to fit
    the style of this thesis.}}
\label{sec:goldendualfullerenes}

\section{Introduction}
\label{sec:introGold}

With the discovery of the catalytic activity of gold
nano-clusters\autocite{Haruta1987,Haruta2003,Haruta2007,Haruta2007a}, research
interest in this field has resurged over the recent
years\autocite{Schwerdtfeger_Goldgoesnano_2003,Hakkinen2008,Maity2012,Zhang2012a,Gong2012,Miao2012,Kyoungweon-2013}.
Gold compounds can show rather interesting topologies, like barrel shaped
structures\autocite{Chen-2015} and planar
sheets,\autocite{Bravo-Perez-1999,Hakkinen2000,Landman2002} mainly because of
strong relativistic effects compared to its lighter congeners copper and
silver\autocite{Pyykko-1988,Schwerdtfeger-2002HA,Pyykko-2004,Pyykko-2007a,Huang-2008,Schwerdtfeger-Lein-2009,pyykko-2012relativistic}.
These effects are also responsible for an unusually high electronegativity,
allowing gold to act as an electron acceptor in mixed-metal
complexes.\autocite{Schwerdtfeger-2002HA} This property could be used for
electronic fine-tuning of physical and chemical properties in gold containing
nano-materials of a certain size.\autocite{Schwerdtfeger_Goldgoesnano_2003} The
growth behaviour of such clusters is, however, still debated
heavily\autocite{Zhao-2010,Barnard-2010,Tian-2011} and even the exact nature of
the transition from planar structures to three-dimensional compact geometries in
small gold clusters is not entirely
resolved.\autocite{Johansson_2D3Dtransitiongold_2008,Fa-Luong-2008,Assadollahzadeh_systematicsearchminimum_2009,Wang-Pal-2010,Wang-Wang-2011,Wang-Wang-2011,Barnard-2012,Gotz_performancedensityfunctional_2013,Kinaci_UnravelingPlanarGlobularTransition_2016}

In 2004 the first hollow gold cluster \ce{Au32} was proposed by
\citeauthor{Johansson_Au3224CaratGolden_2004}\autocite{Johansson_Au3224CaratGolden_2004}
adopting an $I_\text{h}$ symmetric structure that can be created via a dual
transformation of $I_\text{h}$-\ce{C60}, effectively replacing every face in the
carbon fullerene with a gold atom, resulting in a triangulated surface.
\citeauthor{Karttunen_IcosahedralAu72_2008}\autocite{Karttunen_IcosahedralAu72_2008}
have predicted another cage-like gold cluster $I$-\ce{Au72}, which they expect
to be spherically aromatic. For clusters of copper or silver such hollow
structures are not very
stable.\autocite{Johansson_Au3224CaratGolden_2004,FERNANDEZ_DENSITYFUNCTIONALSTUDIES_2006}
The discovery of these types of structures has sparked interest in this field
and many more hollow
cages\autocite{Gu-2004,Fernandez-2006,Fa-Dong-2006,Fa-Zhou-2006,Karttunen_IcosahedralAu72_2008,Fa-Luong-2008,Chen_Structuresneutralanionic_2010,Tian-2011,De-2012,Ning-2014,Joshi-2015}
and clusters enclosing a central metal
atom\autocite{Autschbach_PropertiesWAu12_2004,Zhai-2004,Gao-Bulusu-2005,Wang_Dopinggoldencage_2007,Wang_DopingGoldenBuckyballs_2007,Fa_StructuresMAu16Ag_2008,Munoz-2013,Manna-2013,Tang-2013}
have been found. Most importantly, $I_\text{h}$-\ce{Au32^-},
$T_\text{d}$-\ce{Au16^-}, $C_\text{2v}$-\ce{Au17^-} and
$C_\text{2v}$-\ce{Au18^-} were found to sufficiently explain gas phase
photoelectron spectra of small gold
clusters.\autocite{Ji-2005,Bulusu_Evidencehollowgolden_2006}

On the following pages the relationship between carbon and gold fullerene cages
in terms of their topology is investigated. The similarities arise from the
fact, that topological features known for carbon
fullerenes\autocite{Cataldo-Ori-2011,Schwerdtfeger_topologyfullerenes_2015,Fowler-atlas-2006},
like the Goldberg-Coxeter
transformation,\autocite{Goldberg_ClassMultiSymmetricPolyhedra_1937,Coxeter-1971}
can also be applied to golden dual fullerenes to construct larger structures. A
new class of gold clusters emerges naturally from a one-to-one mapping of the
isomer space of fullerenes to hollow gold clusters. In the following sections,
the stability of such clusters as well as their photoelectron
spectra are investigated. 


%Ever since Haruta discovered that gold nano-clusters are catalytically
%active,\autocite{Haruta1987,Haruta2003,Haruta2007,Haruta2007a} we have
%experienced a new ``gold rush'' in nano
%science\autocite{Schwerdtfeger_Goldgoesnano_2003,Hakkinen2008,Maity2012,Zhang2012a,%Gong2012,Miao2012,Kyoungweon-2013}
%with the discovery of many interesting and often unexpected gold
%nano-structures.\autocite{Chen-2015} Gold shows indeed very unusual properties
%compared to its lighter congeners copper and silver due to pronounced
%relativistic effects within the Group 11 series of
%elements.\autocite{Pyykko-1988,Schwerdtfeger-2002HA,Pyykko-2004,Pyykko-2007a,%Huang-2008,Schwerdtfeger-Lein-2009,pyykko-2012relativistic}
%Albeit these effects increase with the expected $\sim Z^2$ scaling down a group
%in the periodic table, the late transition metals such as gold or mercury have
%rather large relativistic enhancement factors originating from the filling of
%the underlying valence
%$d$-shell.\autocite{Autschbach-2002,Schwerdtfeger-Lein-2009} As a result of
%relativistic effects, smaller gold clusters prefer a planar
%arrangement,\autocite{Bravo-Perez-1999,Hakkinen2000,Landman2002} and mixed
%metal-gold clusters experience strong electron donation toward the gold atoms
%due to its relativistically increased
%electronegativity.\autocite{Schwerdtfeger-2002HA} This makes mixed gold-cluster
%systems ideal for electronically fine-tuning chemical and physical
%properties.\autocite{Schwerdtfeger_Goldgoesnano_2003} Here we mention that the
%transition of 2D gold triangulated networks to 3D compact gold structures
%towards the growth into the fcc bulk gold arrangement is the subject of much
%discussion and
%debate.\autocite{Johansson_2D3Dtransitiongold_2008,Fa-Luong-2008,%Assadollahzadeh_systematicsearchminimum_2009,Wang-Pal-2010,Wang-Wang-2011,%Wang-Wang-2011,Barnard-2012,Gotz_performancedensityfunctional_2013}
%In other words, it is currently challenging to understand the growth of
%metallic clusters toward the bulk by using quantum chemical
%methods.\autocite{Zhao-2010,Barnard-2010,Tian-2011}
%
%Gold clusters can show very unusual and unexpected structures such as the
%pyramidal Au$_{20}$
%cluster\autocite{Li-2003,Fielicke-2008,Assadollahzadeh_systematicsearchminimum_2009}
%or the ``golden fullerene'' $I_\mathrm{h}-$Au$_{32}$ postulated in 2004 by Johansson %et
%al.\autocite{Johansson_Au3224CaratGolden_2004} This unique $I_\mathrm{h}-$Au$_{32}$
%hollow cage can be constructed by replacing each face of the $I_\mathrm{h}-$C$_{60}$
%fullerene polyhedron by a gold atom resulting in a triangulated surface of
%icosahedral symmetry.\autocite{Johansson_Au3224CaratGolden_2004} More recently,
%Karttunen et al. predicted a chiral $I-$Au$_{72}$ cage which is spherically
%aromatic.\autocite{Karttunen_IcosahedralAu72_2008} For both copper and silver
%such a hollow cage becomes rather
%unstable.\autocite{Johansson_Au3224CaratGolden_2004,%FERNANDEZ_DENSITYFUNCTIONALSTUDIES_2006} A number
%of such golden fullerenes, i.e. $I_\mathrm{h}-$Au$_{32}^-$, $T_\mathrm{d}-$Au$_{16}%^-$,
%$C_\mathrm{2v}-$Au$_\mathrm{17}^-$ and $C_\mathrm{2v}-$Au$_{18}^-$, have been %identified by
%photoelectron spectroscopy by Lai-Sheng Wang and
%co-workers.\autocite{Ji-2005,Bulusu_Evidencehollowgolden_2006} Since the
%publication of Johansson et al's
%paper\autocite{Johansson_Au3224CaratGolden_2004}, a number of other studies on
%golden fullerenes appeared, either with a hollow
%cage,\autocite{Gu-2004,Fernandez-2006,Fa-Dong-2006,Fa-Zhou-2006,%Karttunen_IcosahedralAu72_2008,Fa-Luong-2008,Chen_Structuresneutralanionic_2010,%Tian-2011,De-2012,Ning-2014,Joshi-2015}
%or with a central metal
%enclosed\autocite{Autschbach_PropertiesWAu12_2004,Zhai-2004,Gao-Bulusu-2005,%Wang_Dopinggoldencage_2007,Wang_DopingGoldenBuckyballs_2007,Fa_StructuresMAu16Ag_2008,Munoz-2013,%Manna-2013,Tang-2013}
%extending on the original work of Pyykk\"o and Runeberg on
%W@Au$_{12}$.\autocite{Pyykko_IcosahedralWAu12Predicted_2002,%Li_Experimentalobservationconfirmation_2002}
%For a recent review see Wang and Wang.\autocite{Wang-Wang-2012}
%
%In this study we explore the relationship between carbon and golden fullerene
%cages in detail as many interesting topological features known for
%fullerenes,\autocite{Cataldo-Ori-2011,Schwerdtfeger_topologyfullerenes_2015}
%such as the Goldberg-Coxeter transformation to construct larger fullerene
%cages,\autocite{Goldberg_ClassMultiSymmetricPolyhedra_1937,Coxeter-1971} can
%also be applied to the golden fullerenes. We show that a new class of golden
%fullerene structures evolve from a one-to-one mapping into the isomer space of
%fullerene graphs. With this knowledge we re-analyse the experimental
%photoelectron spectrum of the negatively charged Au$_{16}^-$ cage structure. We
%also show that stability of such hollow cage structures is not always
%guaranteed and depends on the sphericity of such systems, but is related to the
%unusual stability of the (111) fcc sheet of gold. We also explore an
%interesting topological relationship between Mackay icosahedra and halma
%transformations recently investigated for
%fullerenes.\autocite{Schwerdtfeger_topologyfullerenes_2015}

%Topological Aspects
\section{\label{sec:TopAsp}Topological Aspects}

The construction of carbon fullerenes can be explained by starting from a
graphene sheet and wrapping it around a sphere\footnote{Or any surface with
genus 0.}, which requires 12 of the hexagonal faces to be replaced by pentagons. This is a requirement imposed by Euler's polyhedral formula.
%
\begin{align}
    |N|-|E|+|F|=\chi
\end{align}
%
Here, $|N|$ is the number of vertices (or atoms), $|E|$ the number of edges (or
bonds), $|F|$ the number of faces and $\chi=2-2g$ the Euler characteristic,
which is $2$ for genus $g=0$ surfaces as in convex polyhedra. As shown in
section~\ref{sec:PlanarGraphs} the number of faces and vertices can be exchanged
without changing the result of Euler's formula. This is also called a dual
transformation, and in the case of the graphene sheet this transformation results in
a (111) \ac{fcc} sheet of, for example, gold bulk. Because the symmetry is
preserved by this transformation both objects belong to the hexagonal 2D lattice
group \textit{p3m1}. To distinguish the two sheets, the graphene sheet will be
denoted \textit{p3m1}-G and the gold sheet \textit{p3m1}-T (figure~\ref{fig:graphenedual}).

\begin{figure}[htb]
    \begin{center}
        \subfloat[\label{subfig:graphene-sheet}]{\includegraphics[width=.4\textwidth,frame]{golddual/sheets/graphene.png}}\hspace{.05\textwidth}
        \subfloat[\label{subfig:gold-sheet}]{\includegraphics[width=.4\textwidth,frame]{golddual/sheets/gold.png}}
        \caption{\protect\subref{subfig:graphene-sheet} \textit{p3m1}-G graphene
        and \protect\subref{subfig:gold-sheet} its dual sheet \textit{p3m1}-T
        adopted in the (111) surface of fcc gold.}
    \label{fig:graphenedual}
    \end{center}
\end{figure}

Small cut-outs of the \textit{p3m1}-T sheet can be found as global minima for smaller
gold clusters, indicating that this represents a very stable structural motif
for gold compounds.\autocite{Assadollahzadeh_systematicsearchminimum_2009}

This concept can be extended to non-spherical structures like carbon nano-tubes
to construct gold nanowires, and there has been experimental evidence
supporting the existence of such structures\autocite{Kondo-2000}. Because they
are the duals of the carbon nano-tubes they can be constructed in the same
way.\autocite{Dresselhaus-1992} Two examples of cylindrically shaped carbon and gold structures are shown in figure~\ref{fig:nanotubedual}.

\begin{figure}[htb]
    \begin{center}
        \subfloat[\label{subfig:nanotube}]{\includegraphics[width=.8\textwidth]{golddual/C144.png}} \\
        \subfloat[\label{subfig:gold-nanotube}]{\includegraphics[width=.8\textwidth]{golddual/Au74.png}}
        \caption{\protect\subref{subfig:nanotube} $D_\mathrm{6d}-$C$_{144}$ zig-zag
        fullerene nanotube and \protect\subref{subfig:gold-nanotube} its dual
        $D_\mathrm{6d}-$Au$_{74}$.}
    \label{fig:nanotubedual}
    \end{center}
\end{figure}

As fullerenes need to have exactly 12 pentagons, a dual fullerene will have 12
vertices of degree five instead. All other vertices will have degree six and
there are exactly as many as there are hexagons in the corresponding carbon
fullerene. The smallest carbon fullerene \ce{C20} has $|F_h|=0$ hexagons, and
all larger ones at least $|F_h|>1$. The fullerene \ce{C22}, which would contain
exactly one hexagon, is
non-existent,\autocite{Grunbaum_numberhexagonssimplicity_1963} thus, the
hypothetical \ac{GDF} \ce{Au13}\footnote{The relation between the
number of vertices in a fullerene $|N_f|$ and the number of vertices in the
corresponding dual fullerene $|N_d|$ is
$|N_d|=|F_f|=|N_f|/2+2$.\autocite{Schwerdtfeger_topologyfullerenes_2015}} also
cannot exist. Fullerenes often have much more than one stable isomer
(non-isomorphic graphs)\autocite{Fowler-atlas-2006} and because of the dual
relationship there should be as many isomers for the \acp{GDF}.
Additionally, the growth of this isomer space for fullerenes should scale the same with respect
to the number of vertices, which was found to be
$\mathcal{O}\left({|N|}^9\right)$.\autocite{Thurston_Shapespolyhedratriangulations_1998}

Both \ce{C60} and its dual \ce{Au32}, as well as their graph representations are
depicted in figure~\ref{fig:C60dual}. This relationship was first noticed in
conjunction with the prediction of
\ce{Au32},\autocite{Johansson_Au3224CaratGolden_2004} and it allows the usage of
the same algorithms used to construct fullerenes to create \acp{GDF}. For
example, using the generalised face-spiral
algorithm\autocite{Fowler-atlas-2006,Schwerdtfeger_Programfullerenesoftware_2013,Wirz-2014,Schwerdtfeger_topologyfullerenes_2015}
followed up by an embedding of the graph on a genus 0 surface and a dual
transformation.
%
\begin{figure}[htb]
    \begin{center}
        \subfloat[\label{subfig:c60graph}]{\includegraphics[width=.216\textwidth]{golddual/C60graph.png}}\hspace{0.05\textwidth}
        \subfloat[\label{subfig:c60ih}]{\includegraphics[width=.25\textwidth]{golddual/C60Ih.png}}\hspace{0.034\textwidth}
        \subfloat[\label{subfig:au32ih}]{\includegraphics[width=.25\textwidth]{golddual/Au32Ih.png}}
        \caption{\protect\subref{subfig:c60graph} Schlegel diagram of C$_{60}$ (red
        vertices) and its dual (blue vertices and dashed edges),
        \protect\subref{subfig:c60ih} the C$_{60}$ structure, and
        \protect\subref{subfig:au32ih} its dual Au$_{32}$ structure.}
    \label{fig:C60dual}
    \end{center}
\end{figure}
%

A recent investigation of photoelectron spectra of gold clusters considered the
existence of a $T_\text{d}$-\ce{Au16^-} cluster to explain the experimental
findings.\autocite{Bulusu_Evidencehollowgolden_2006} This cluster would be dual
to \ce{C28} and has exactly two isomers: $T_\text{d}$-\ce{Au16} and
$D_2$-\ce{Au16}. In the above study, the $D_2$-symmetric isomer has not been
considered to explain the observed spectra, which naturally raises the question
whether this isomer's photoelectron spectra is similar or even capable of
explaining the observations better. 

The question of which structure is dominating the experimental spectrum is
closely related to the question of which structure is thermodynamically more stable. For regular
carbon fullerenes there exists an ``isolated pentagon rule'', that states that a
carbon fullerene is more stable when none of the pentagons are in direct contact
with each other.\autocite{Kroto_stabilityfullerenesCn_1987} It is hitherto
unknown if there is an equivalent ``isolated vertex of degree five rule'' for
dual fullerene structures.

As mentioned before, methods like the Goldberg-Coxeter transformation can be used
to construct larger dual fullerenes from smaller
ones.\autocite{Goldberg_ClassMultiSymmetricPolyhedra_1937,Coxeter-1971,Dutour_GoldbergCoxeterConstructionvalent_2004}
The original Goldberg-Coxeter transformation was carried out on the dodecahedron (\ce{C20}
fullerene)\autocite{Goldberg_ClassMultiSymmetricPolyhedra_1937,Coxeter-1971}, but
it can be shown that it can be applied to any fullerene
graph.\autocite{Schwerdtfeger_topologyfullerenes_2015} The transformation
$GC_{k,l}$ can be controlled by two integer parameters $k,l$ describing the
scaling and rotation of the mesh on which the transformation is carried out. The
symmetry of the original fullerene is preserved if $k=l$ (leapfrog
transformation) or $l=0$  (halma transformation). Some important transformations
are for example
$GC_{1,1}$[$I_\mathrm{h}-$C$_{20}$]=$I_\mathrm{h}-$C$_{60}$\autocite{Fowler-atlas-2006}
and $GC_{2,0}$[$I_\mathrm{h}-$C$_{20}$]=$I_\mathrm{h}-$C$_{80}$, both preserving
the initial point group symmetry. In case of the gold clusters the same
transformations results in the respective dual representations, i.e.
$GC_{1,1}$[$I_\mathrm{h}-$Au$_{12}$]=$I_\mathrm{h}-$Au$_{32}$ and
$GC_{2,0}$[$I_\mathrm{h}-$Au$_{12}$]= $I_\mathrm{h}-$Au$_{42}$. Both of these
structures have been proposed previously to be stable
hollow cages.\autocite{Johansson_Au3224CaratGolden_2004} The new vertex count of
$GC_{k,l}[\ce{Au_{|N_d|}}]$ is
%
\begin{align}
    |N_\text{d}'| = (k^2+kl+l^2)(|N_\text{d}|-2) + 2.\label{eqn:dualvertex}
\end{align}

An often encountered structural motif in gold clusters is the Mackay
icosahedron.\autocite{Nam2002,Wang-Wang-2011} Although this is not a hollow
structure, it is related to dual fullerenes as it is made up of multiple
icosahedral shells. Each individual shell $m$ consists of
%
\begin{align}
    |N_\text{shell}|=10m^2+2\label{eqn:MackayShells}
\end{align}
%
atoms, resulting, when summing up, in the magical cluster numbers 13, 55, 147, 309, and so
on.\autocite{Mackay-1962,Kuo_MackayAntiMackayDoubleMackay_2002}
Figure~\ref{fig:mackaylarge} shows one such icosahedral structure with $m=7$
shells and 1415 atoms.
%
\begin{figure}[htb]
    \begin{center}
    \includegraphics[width=4.9cm]{golddual/ico.jpg}
        \caption{Mackay icosahedron with 7 shells and 1415 atoms. The outer icosahedral shell is the dual of the halma transform $GC_{7,0}$[$I_\mathrm{h}-$C$_{20}$]=$I_\mathrm{h}-$C$_{980}$.}
    \label{fig:mackaylarge}
    \end{center}
\end{figure}
%
The number of shells can be deduced from the number of spheres on one edge of
the icosahedron,  including the spheres marked in red. There is exactly one
sphere more on the edges than there are shells, thus $m=|N_\text{edge}|-1$. The
halma pattern of a $GC_{k,0}$ transformation is clearly visible on the faces of
the icosahedron, and it turns out the icosahedral shells are in fact related to
the smallest fullerene \ce{C20} by such a transformation and a subsequent
dualisation. For this process equation~\eqref{eqn:dualvertex} becomes
%
\begin{align}
    |N_\text{d}'| = k^2(|N_\text{d}|-2) + 2,
\end{align}
%
and with $|N_\text{d}|=12$ (as \ce{Au12} is the dual of \ce{C20}) this is equal
to equation~\eqref{eqn:MackayShells}. The parameter $k$ of the transformation
$GC_{k,0}[\ce{C20}]$ therefore defines which shell of the icosahedron is
created by the Goldberg-Coxeter transformation with subsequent dualisation.

The relationship between carbon fullerenes and hollow gold clusters can be used
to name the latter in the same way as the carbon fullerenes. For example, this
can be achieved by using the canonical face spiral pentagon indices (FSPI) in
combination with the numbering scheme introduced by
Manolopoulos.\autocite{Fowler-atlas-2006} A complete and unique method for
naming polyhedra extending the original algorithm by Manolopoulos has been
developed recently.\autocite{Wirz_Namingpolyhedrageneral_2018} In the following
sections the golden dual fullerenes from \ce{Au12} to \ce{Au20} (excluding
\ce{Au13}) will be investigated by means of \ac{DFT} calculations.

%Fullerenes show rich and mathematically interesting topological
%features\autocite{Cataldo-Ori-2011,Schwerdtfeger_topologyfullerenes_2015},
%which have been described for example in the works of Fowler and
%Manolopoulos\autocite{Fowler-atlas-2006}, and most recently by our research
%group in Auckland.\autocite{Schwerdtfeger_topologyfullerenes_2015} They can be
%thought of by wrapping a graphene sheet around a sphere (or more generally a
%genus 0 surface), but introducing 12 pentagons (no more and no less) to fulfil
%Euler's polyhedral formula, 
%%
%\begin{equation}
%  \label{eq:euler} 
%  N + F -E = \chi 
%\end{equation}
%%
%where $N$ is the number of vertices (atoms), $F$ is the number of faces
%(rings), $E$ is the number of edges (bonds) and $\chi$ is called the Euler
%characteristic with $\chi=2$ for convex
%polyhedra.\autocite{Kotschick_TopologyCombinatoricsSoccer_2006} Euler's formula
%already shows the symmetry between the number of vertices $N$ and the number of
%faces $F$, as their role can be interchanged without violating Euler's theorem.
%Interchanging the roles of vertices and faces in a graphene sheet leads to a
%(111) sheet (surface) of an (for example) fcc structure adopted in bulk gold
%(both belonging to the hexagonal 2D lattice group \textit{p3m1}), where the
%dual vertex is in the center of the hexagon connected by edges to the
%neighboring dual vertices. Several smaller gold clusters found in the search
%for global minima are in fact cut-outs from this (111) fcc
%sheet,\autocite{Assadollahzadeh_systematicsearchminimum_2009} denoted as
%\textit{p3m1}-T in the following (see figure~\ref{fig:graphenedual}). 

%To view
%it in a different way, the hexagons in the graphene sheets are exactly the
%Voronoi cells in the \textit{p3m1}-T sheet. As an interesting side aspect we
%mention the helical multi-shell (chiral) gold nanowires found experimentally by
%Kondo and Takayanagi,\autocite{Kondo-2000} which are duals of multi-shell
%(chiral) carbon nanotubes.\autocite{Johansson_Au3224CaratGolden_2004} These
%gold nanowires can be constructed exactly in the same way as carbon nanotubes
%using the chiral vector $C_h(n,m)$ on a hexagonal sheet as described in detail
%for example by Dresselhaus and co-workers.\autocite{Dresselhaus-1992} As an
%example we show the $D_\mathrm{6d}$ fullerene nanotube and its dual structure in
%figure~\ref{fig:nanotubedual}.

%The requirement to have 12 pentagons in a fullerene graph with $F_h$ hexagons
%($F_h$=0 for C$_{20}$ and $F_h > 1$ for all other fullerenes) implies for a
%fullerene dual to have exactly 12 vertices of degree five and the remaining of
%degree six. In fact it is well known that C$_{22}$ cannot exist as a
%fullerene,\autocite{Grunbaum_numberhexagonssimplicity_1963} which implies that
%its hypothetical dual Au$_{13}$ does not exist either (the number of vertices
%$N_d$ in the dual is identical to the number of faces in a fullerene,
%$N_d=F_f=N_f/2+2$,\autocite{Schwerdtfeger_topologyfullerenes_2015} using
%symbols $f$ and $d$ for the fullerene and its dual respectively). Because there
%is a one-to-one correspondence between a fullerene and its dual graph, we have
%as many isomers (nonisomorphic graphs) for C$_{N_f}$ as we have for
%Au$_{N_f/2+2}$ and dualization preserves the point group symmetry. Here we
%mention that according to Thurston, the number of isomers increases
%polynomially in ninth leading order with the number of vertices, i.e.
%$\sim\mathcal{O}({N_f^9})$.\autocite{Thurston_Shapespolyhedratriangulations_1998} 

%Au$_{32}$ was the first of such golden fullerenes postulated by Johansson et
%al. to be a rather stable hollow cluster, and they were the first ones
%mentioning that these golden dual fullerenes (GDF) are obtained from fullerene
%graphs.\autocite{Johansson_Au3224CaratGolden_2004} Au$_{32}$ is shown in
%figure~\ref{fig:C60dual} together with its dual, C$_{60}$ and the corresponding
%graph representation (twice the dual transformation leads back to the original
%polyhedron or graph). Now that we established an isomorphism between a
%fullerene graph and its dual, we can easily construct isomers of golden
%fullerenes by using standard algorithms for the construction of fullerenes,
%such as the generalized face-spiral
%algorithm,\autocite{Fowler-atlas-2006,Schwerdtfeger_Programfullerenesoftware_2013,%Wirz-2014,Schwerdtfeger_topologyfullerenes_2015}
%embedding the graph on a genus 0 surface and finally transforming the cage to
%its dual.



%As an example we mention Au$_{16}$ as the dual of C$_{28}$. Checking the list
%of possible isomers\autocite{Brinkmann_HouseGraphsdatabase_2013} we see that
%there are two possible non-isomorphic structures, $D_2-$Au$_{16}$ and
%$T_\mathrm{d}-$Au$_{16}$. However, only the more symmetric $T_\mathrm{d}-$Au$_{16}%^-$ has been
%considered as a possible candidate in recent photoelectron spectroscopy
%experiments.\autocite{Bulusu_Evidencehollowgolden_2006} The question naturally
%arises if the other negatively charged $D_\mathrm{2}$ isomer has a similar %photoelectron
%spectrum and is more stable or not compared to the $T_\mathrm{d}$ isomer. In fact,
%Au$_{32}$ which is the dual of C$_{60}$ has 1812 different isomers with only
%one fulfilling the isolated pentagon rule (IPR) as proposed by
%Kroto.\autocite{Kroto_stabilityfullerenesCn_1987} For example, in a recent
%paper Fa and Dong reported on a hollow gold $D_\mathrm{6d}-$Au$_{26}$
%cluster.\autocite{Fa-Dong-2006,Fa-Luong-2008} Looking at the list of possible
%fullerenes we see that there are 199 possible isomers for dual fullerene
%structures, and in fact there are two possible isomers having $D_\mathrm{6d}$
%symmetry. This just highlights the rich topology of such dual fullerenes.

%For the dual structures we do not know if a similar rule applies, that is an
%``isolated vertex rule'' of degree five (IVR5). In fullerenes the pentagons are
%responsible for the curvature of the carbon cage and for the overall symmetry
%and structure, with connected hexagons building planar sub-structures on the
%polyhedron. Here we mention that the Mackay icosahedron (discussed below) shows
%exactly that feature.

%\textit{p3m1}-G sheets have been considered by Goldberg and Coxeter for the
%construction of larger
%fullerenes.\autocite{Goldberg_ClassMultiSymmetricPolyhedra_1937,Coxeter-1971,%Dutour_GoldbergCoxeterConstructionvalent_2004}
%The original Goldberg-Coxeter transformation superimposes a hexagonal mesh on
%the surface of the C$_{20}$ dodecahedron forming a new polyhedron with leaving
%the number of pentagons at exactly
%12.\autocite{Goldberg_ClassMultiSymmetricPolyhedra_1937,Coxeter-1971} This
%transformation can be applied to any fullerene
%isomer.\autocite{Schwerdtfeger_topologyfullerenes_2015} The Goldberg-Coxeter
%transformation $GC_{k,l}$ increases the number of vertices for a fullerene by a
%factor of $(k^2+kl+l^2)$, where $k$ and $l$ are integers describing the scale
%and orientation of the
%mesh.\autocite{Dutour_GoldbergCoxeterConstructionvalent_2004,%Schwerdtfeger_topologyfullerenes_2015}
%If $k=l$ or $l=0$ the point group symmetry is preserved. For example, we have
%$GC_{1,1}$[$I_\mathrm{h}-$C$_{20}$]=$I_\mathrm{h}-$C$_{60}$ (leapfrog
%transformation)\autocite{Fowler-atlas-2006} and
%$GC_{2,0}$[$I_\mathrm{h}-$C$_{20}$]=$I_\mathrm{h}-$C$_{80}$ (halma transformation), or in the
%dual case applying the same procedure to the \textit{p3m1}-T sheet for our gold
%fullerenes $GC_{1,1}$[$I_\mathrm{h}-$Au$_{12}$]=$I_\mathrm{h}-$Au$_{32}$ and
%$GC_{2,0}$[$I_\mathrm{h}-$Au$_{12}$]= $I_\mathrm{h}-$Au$_{42}$. Simple algebra shows that
%$GC_{k,l}$[Au$_{N_d}$] has a new vertex count of
%%
%\begin{equation}
%  \label{eq:dualvertex} 
%N_d'=(k^2+kl+l^2)(N_d-2)+2 
%\end{equation}
%
%Both $I_\mathrm{h}-$Au$_{32}$ and $I_\mathrm{h}-$Au$_{42}$ have been postulated as stable %golden
%fullerenes before,\autocite{Johansson_Au3224CaratGolden_2004} and more recently
%the chiral $I-$Au$_{72}$\autocite{Karttunen_IcosahedralAu72_2008} which is
%nothing else as the dual of $GC_{2,1}$[$I_\mathrm{h}-$C$_{20}$] = $I-$C$_{140}$, which
%is chiral as well as symmetry is conserved upon dualization.

%Now it almost seems trivial to relate Mackay icosahedra\autocite{Mackay-1962}
%well known for gold clusters\autocite{Nam2002,Wang-Wang-2011} to fullerenes. We
%might call them "multi-walled halma-transformed icosahedral" dual fullerenes in
%analogous way to the multi-walled gold nanowires. A Mackay icosahedron is a
%closed packed multi-shell structure each shell being an icosahedron with
%%
%\begin{equation}
%  \label{eq:mackay} 
%N_\text{shell}=10k^2+2
%\end{equation}
%%
%number of atoms in each shell with increasing $k$. This gives the well known
%magic cluster numbers of (including the central atom) 13, 55, 147, 309, 561,
%... derived from
%%
%\begin{equation}
%  \label{eq:mackaytotal} 
%N_\text{total}=1+2\sum_{k=1}^{m_\text{shell}}\left( 5k^2+1 \right)=(10m_\text{shell}^3%+15m_\text{shell}^2+11m_\text{shell}+3)/3
%\end{equation}
%
%with $m_\text{shell} \ge 1$. One such Mackay icosahedron with $m_\text{shell}$=7 and
%$N_\text{total}$=1415 is shown in figure~\ref{fig:mackaylarge}. The triangles
%clearly show the halma pattern of a Goldberg-Coxeter $GC_{k,0}$ transformation.
%Because each transformation brings a new vertex on the icosahedral edge, we can
%just deduct the number of shells by counting the number of atoms at one edge of
%a triangle $N_\text{edge}$, i.e. $m_\text{shell}=N_\text{edge}-1$, which gives 8 atoms and 7
%icosahedral shells for the Mackay icosahedron shown in figure~\ref{fig:mackaylarge}.
%%
%
%%
%Mackay pointed out that the packing density (or atomic packing factor) for $N
%\rightarrow \infty$ with $\rho$=0.68818 is not too different from a closed
%packed structure such as fcc with
%$\rho=\pi/\sqrt{18}=0.74048$,\autocite{Mackay-1962} one reason why icosahedral
%cluster growth is often seen. The number of atoms in a shell $N_\text{shell}$
%directly corresponds to the number of faces in a halma transformed C$_{20}$,
%i.e. for $l=0$ in the Goldberg-Coxeter transformation and starting with
%$N_d=12$ in equation~\eqref{eq:dualvertex} we have $N_d'=k^2(N_d-2)+2 = 10k^2+2$
%identical with the formula given by Mackay. 

%Finally we mention that we can name the different isomers of the golden dual
%fullerenes exactly in the same way as we do for the fullerenes by using the
%canonical face spiral pentagon indices (FSPI) and the numbering scheme
%introduced by Fowler and Manolopoulos,\autocite{Fowler-atlas-2006} keeping in
%mind that the face spiral for fullerenes now becomes a vertex spiral for the
%dual triangulated surface. We are now turning to a detailed analysis of all
%possible golden dual fullerenes from $I_\mathrm{h}-$Au$_{12}$, Au$_{14}$ to Au$_{20}$
%and $I_\mathrm{h}-$Au$_{32}$ by quantum chemical calculations.

\section{Computational Details}

Program \textsc{Fullerene}\autocite{Schwerdtfeger_Programfullerenesoftware_2013}
has been used to construct initial structures of all isomers of the golden dual
fullerenes from Au$_{12}$ to Au$_{20}$ using a recently developed force-field
for fullerenes\autocite{Wirz_smallfullerenesgraphene_2015} (excluding the
non-existing golden dual fullerene Au$_{13}$). The following isomers need to be
considered according to the isomer list for the fullerenes (number in
parenthesis gives the number of different isomers of same
symmetry):\autocite{Brinkmann_HouseGraphsdatabase_2013,Schwerdtfeger_Programfullerenesoftware_2013}
$I_\mathrm{h}-$Au$_{12}$, $D_\mathrm{6d}-$Au$_{14}$, $D_\mathrm{3h}-$Au$_{15}$,
$D_\mathrm{2}-$Au$_{16}$, $T_\mathrm{d}-$Au$_{16}$, $D_\mathrm{5h}-$Au$_{17}$,
$C_\mathrm{2v}-$Au$_{17}$(2), $D_\mathrm{3h}-$Au$_{18}$,
$D_\mathrm{3d}-$Au$_{18}$, $D_\mathrm{3}-$Au$_{18}$, $D_\mathrm{2}-$Au$_{18}$,
$C_\mathrm{3}-$Au$_{18}$(2), $C_\mathrm{3v}-$Au$_{19}$,
$C_\mathrm{2}-$Au$_{19}$(3), $C_\mathrm{s}-$Au$_{19}$(2),
$D_\mathrm{6h}-$Au$_{20}$, $D_\mathrm{3h}-$Au$_{20}$,
$D_\mathrm{2d}-$Au$_{20}$(2), $C_\mathrm{2v}-$Au$_{20}$,
$D_\mathrm{2}-$Au$_{20}$(2), $C_\mathrm{2}-$Au$_{20}$(3),
$C_\mathrm{2}-$Au$_{20}$(2), $C_\mathrm{1}-$Au$_{20}$(2) and
$I_\mathrm{h}-$Au$_{32}$. The initial force-field optimised structures scaled to
an approximate internuclear distance were then refined by using the
Predew-Burke-Ernzerhof \ac{GGA}
functional\autocite{Perdew_GeneralizedGradientApproximation_1996,Perdew_GeneralizedGradientApproximation_1997}
corrected for dispersion interactions using Grimme's method
(PBE-D3)\autocite{Grimme_consistentaccurateinitio_2010,Grimme_Effectdampingfunction_2011}
together with a Los-Alamos scalar relativistic effective core potential for gold
and the accompanying double-zeta basis sets.\autocite{Wadt1985} Note that the
PBE functional was recently considered to perform well for gold
clusters.\autocite{Mancera_alternativemethodologyassess_2015} For several
selected clusters the geometries obtained were checked for accuracy by carrying
out calculations using a small core scalar relativistic Stuttgart
pseudopotential\autocite{Figgen_Energyconsistentpseudopotentialsgroup_2005}
together with an augmented valence double-zeta basis set by Peterson and
Puzzarini.\autocite{Peterson-2005} For comparison, the compact global minimum
cluster structures recently published for the neutral
compounds\autocite{Assadollahzadeh_systematicsearchminimum_2009} and for the
negatively charged
species\autocite{Schooss_Determiningsizedependentstructure_2010,Lechtken_Structuredeterminationgold_2009}
were calculated.

The simulation of the photoelectron spectra has been carried out by artificial
broadening the spectrum of orbital energies with Gaussian functions. The
standard deviation $\sigma$ for these functions was chosen to be 0.035~eV in
qualitative agreement with the experimental spectra. The orbital energies were
calculated using the PBE density functional with the
def2-SVP\autocite{Weigend_Balancedbasissets_2005} double-zeta basis implemented
in \textsc{Turbomole} 7.0.\autocite{_TURBOMOLEV72015_} The core region was
described using an effective core potential including scalar relativistic
effects.  The calculated electron affinities were used as the onset value for
simulating the photoelectron spectra. 


For the calculation of the (111) fcc sheet and the fcc bulk structure of gold
the program package VASP5\autocite{Kresse_Efficiencyabinitiototal_1996} was
used, utilizing a plane-wave basis set (cut-off energy $E_c=350$~eV) and the
standard \ac{PAW} datasets for the elements to model the electron-ion
interaction\autocite{Blochl_Projectoraugmentedwavemethod_1994,Kresse_ultrasoftpseudopotentialsprojector_1999}.
The electron-electron interaction was modelled within the \ac{GGA} to the
exchange-correlation energy functional as described above and dispersive effects
were taken into account by employing Grimme's D3 dispersion correction with
Becke-Johnson
damping.\autocite{Grimme_consistentaccurateinitio_2010,Grimme_Effectdampingfunction_2011}
Brillouin zone integrations were carried out on $\Gamma$-centred Monkhorst-Pack
grids of $k$-points with a distance of 0.2~\AA$^{-1}$. The cohesive energy is
defined as the atomisation energy per atom keeping in mind that one gold atom is
negatively charged for the anionic clusters.

In order to discuss how much the gold cages deviate from sphericity compared to
the dual fullerene structure, the previously introduced definition of a \ac{MDS}
was used,\autocite{Schwerdtfeger_Programfullerenesoftware_2013}
%
\begin{equation} 
\min\limits_{c_\mathrm{MDS} \in \mathrm{CH}(S)} \frac{1}{N} \sum _{i} \left|R_\mathrm{MDS} -\| \mathbf{p}_{i}-\mathbf{c}_\mathrm{MDS} \| \right|  
\end{equation}
with the \ac{MDS} radius defined as
\begin{equation} 
	R_{\mathrm{MDS}} =\frac{1}{N} \sum _{i}\| \mathbf{p}_{i} -\mathbf{c}_{\mathrm{MDS}} \|. 
	\label{eq:RMDS}
\end{equation}
%
Here $S$ is the set of $n$ points $\mathbf{p}_i$ ($i=1,\ldots ,n$) in
$3$-dimensional space, $\mathrm{CH}(S)$ its convex hull, $\|\cdot\| $ the
Euclidean norm, and $\mathbf{c}_\mathrm{MDS}$ is the barycentre of the \ac{MDS}
with radius $R_\mathrm{MDS}$. In other words, the procedure tries to locate a
sphere that approximates the position of the vertices well. A measure for
distortion from spherical symmetry through the
\ac{MDS} is defined as\autocite{Schwerdtfeger_Programfullerenesoftware_2013}
%
\begin{equation}
  \label{eq:DMDS}
  D_{\mathrm{MDS}} = \frac{100}{N R_\mathrm{min}} \sum_{i=1}^N \left|R_{\mathrm{MDS}} - \|\mathbf{p}_i - \mathbf{c}_{\mathrm{MDS}}\| \right|,
\end{equation}
where $R_\mathrm{min}$ is the smallest bond distance found in the cluster. The pentagon index $N_p$ is defined as
\begin{equation}
  \label{pentindex}
  N_p = \frac{1}{2}\sum_{k=1}^{5} kp_k \quad \text{ with } \quad  \sum_{k=0}^{5} p_k = 12
\end{equation}
%
where the pentagon indices $(p_i | i=0, \dots , 5)$ define the number of
pentagons attached to another pentagon.\autocite{Fowler-atlas-2006}

\section{Structure and Stability}

The results for the neutral and negatively charged gold clusters are collected
in tables~\ref{tab:neutral} and \ref{tab:anion} respectively. The dual fullerene
structures are compared to the known global minimum structures in these tables,
and the different isomers are numbered according to their canonical degree 5
vertex spiral, identical to the canonical face spiral pentagon indices for
fullerenes.\autocite{Fowler-atlas-2006} Calculations for the most stable neutral
and anionic compact Au$_n$ clusters for comparison are also included and are
listed in table~\ref{tab:Aun}. The investigated structures for the negatively
charged gold clusters are depicted in figures~\ref{fig:Au1219-} and
\ref{fig:Au20-}, and the energy differences compared to the global minimum
structures are shown in figure~\ref{fig:AunMinus2}.
%
%
\begin{figure}[htbp]
	\begin{center}
	\subfloat[12:1 ]{\includegraphics[width=0.25\textwidth]{golddual/anions/Au12.png}}
	\subfloat[14:1 ]{\includegraphics[width=0.25\textwidth]{golddual/anions/Au14.png}}
	\subfloat[15:1 ]{\includegraphics[width=0.25\textwidth]{golddual/anions/Au15.png}}
	\subfloat[16:1 ]{\includegraphics[width=0.25\textwidth]{golddual/anions/Au16-D2.png}}\\
	\subfloat[16:2 ]{\includegraphics[width=0.25\textwidth]{golddual/anions/Au16-Td.png}}
	\subfloat[17:1 ]{\includegraphics[width=0.25\textwidth]{golddual/anions/Au17-D5h.png}}
	\subfloat[17:2 ]{\includegraphics[width=0.25\textwidth]{golddual/anions/Au17-C2v-1.png}}
	\subfloat[17:3 ]{\includegraphics[width=0.25\textwidth]{golddual/anions/Au17-C2v-2.png}}\\
	\subfloat[18:1 ]{\includegraphics[width=0.25\textwidth]{golddual/anions/Au18-C2-1.png}}
	\subfloat[18:2 ]{\includegraphics[width=0.25\textwidth]{golddual/anions/Au18-D2.png}}
	\subfloat[18:3 ]{\includegraphics[width=0.25\textwidth]{golddual/anions/Au18-D3d.png}}
	\subfloat[18:4 ]{\includegraphics[width=0.25\textwidth]{golddual/anions/Au18-C2-2.png}}\\
	\subfloat[18:5 ]{\includegraphics[width=0.25\textwidth]{golddual/anions/Au18-D3h.png}}
	\subfloat[18:6 ]{\includegraphics[width=0.25\textwidth]{golddual/anions/Au18-D3.png}}
	\subfloat[19:1 ]{\includegraphics[width=0.25\textwidth]{golddual/anions/Au19-C2-1.png}}
	\subfloat[19:2 ]{\includegraphics[width=0.25\textwidth]{golddual/anions/Au19-Cs-1.png}}\\
	\subfloat[19:3 ]{\includegraphics[width=0.25\textwidth]{golddual/anions/Au19-Cs-2.png}}
	\subfloat[19:4 ]{\includegraphics[width=0.25\textwidth]{golddual/anions/Au19-C2-2.png}}
	\subfloat[19:5 ]{\includegraphics[width=0.25\textwidth]{golddual/anions/Au19-C2-3.png}}
	\subfloat[19:6 ]{\includegraphics[width=0.25\textwidth]{golddual/anions/Au19-C3v.png}}
	\caption{Structures of anionic gold clusters (Au$_{12}^-$ to Au$_{19}^-$).}
	\label{fig:Au1219-}
\end{center}
\end{figure}
%
\begin{figure}[htbp]
	\begin{center}
	\subfloat[20:1 ]{\includegraphics[width=0.25\textwidth]{golddual/anions/Au20-C2-1.png}}
	\subfloat[20:2 ]{\includegraphics[width=0.25\textwidth]{golddual/anions/Au20-D2-1.png}}
	\subfloat[20:3 ]{\includegraphics[width=0.25\textwidth]{golddual/anions/Au20-C1-1.png}}
	\subfloat[20:4 ]{\includegraphics[width=0.25\textwidth]{golddual/anions/Au20-Cs-1.png}}\\
	\subfloat[20:5 ]{\includegraphics[width=0.25\textwidth]{golddual/anions/Au20-D2-2.png}}
	\subfloat[20:6 ]{\includegraphics[width=0.25\textwidth]{golddual/anions/Au20-D2d-1.png}}
	\subfloat[20:7 ]{\includegraphics[width=0.25\textwidth]{golddual/anions/Au20-C1-2.png}}
	\subfloat[20:8 ]{\includegraphics[width=0.25\textwidth]{golddual/anions/Au20-Cs-2.png}}\\
	\subfloat[20:9 ]{\includegraphics[width=0.25\textwidth]{golddual/anions/Au20-C2v.png}}
	\subfloat[20:10]{\includegraphics[width=0.25\textwidth]{golddual/anions/Au20-C2-2.png}}
	\subfloat[20:11]{\includegraphics[width=0.25\textwidth]{golddual/anions/Au20-C2-3.png}}
	\subfloat[20:12]{\includegraphics[width=0.25\textwidth]{golddual/anions/Au20-C2-4.png}}\\
	\subfloat[20:13]{\includegraphics[width=0.25\textwidth]{golddual/anions/Au20-D3h.png}}
	\subfloat[20:14]{\includegraphics[width=0.25\textwidth]{golddual/anions/Au20-D2d-2.png}}
	\subfloat[20:15]{\includegraphics[width=0.25\textwidth]{golddual/anions/Au20-D6h.png}}
	\caption{Structures of anionic gold clusters (Au$_{20}^-$).}
	\label{fig:Au20-}
\end{center}
\end{figure} 

The optimised gold clusters can be sorted according to whether they can be
derived from a dual fullerene structure, or more generally from a cubic
polyhedral graph, or not. In this case Euler's polyhedral formula can be
simplified, which upon dualisation gives a triangulation of a sphere obeying the
formula
%
\begin{equation}
\label{eq3valent} 
\Gamma=\sum_{n=3}(6-n)|N_n| = 12,
\end{equation}
%
where $|N_n|$ denotes the number of $n$-valent vertices. Any deviation from
$\Gamma=12$ implies that the polyhedron is not a triangulation of a sphere. As mentioned before, for
dual fullerenes only values of $N_5=12$ and $N_6=\{0,2,3,4,5,\dots\}$ are
allowed. Hence, a true dual fullerene structure is obtained in case of a
complete triangulation and 12 vertices of degree five.
%
\begin{table}[htbp]
	\centering
    \setlength{\tabcolsep}{1.5pt}
    \footnotesize{
    \caption{Topological parameters for the neutral gold clusters. Number of
    gold atoms and isomer numbers of the corresponding fullerene in canonical
    order of the pentagon spiral indices,\autocite{Fowler-atlas-2006} ideal and
    actual point group symmetry, energy differences $\Delta E_g$ to the most
    stable neutral cluster of same size and binding energy per atom $\Delta E_n
    = [E(\textrm{Au}_n)-nE(\textrm{Au})]/n$  (in eV), shortest and largest bond
    distance (in \AA), pentagon index (PI) $N_p$, and distortion parameter $D$
    (in \%) for the initial force-field optimised fullerene structure (F) and
    the \acs{GDF}.}
	\label{tab:neutral}
	\begin{tabular}{lllrrrrrrrrrrrr}
\toprule
\multicolumn{1}{c}{  } & \multicolumn{2}{c}{ symmetry  }  & \multicolumn{2}{c}{stability} & \multicolumn{4}{c}{ vertices } & & \multicolumn{2}{c}{ bondlengths } &  PI & \multicolumn{2}{c}{ $D$ } \\
isomer & ideal  & actual  & $\Delta E_n$ &$\Delta E_g$ & \multicolumn{1}{c}{$|N_4|$} & \multicolumn{1}{c}{$|N_5|$} & \multicolumn{1}{c}{$|N_6|$} & \multicolumn{1}{c}{$|N_7|$} & $\Gamma$ & shortest & largest  & $N_p$ & F & GDF\\\midrule
12:1    & $I_\mathrm{h}$  & $D_\mathrm{4h}$ & $-2.058$ & $0.485$  & $8$ & $0$  & $4$      & $0$ & 16  & $2.798$ & $2.895$ & 30  & 0    & 21.1  \\
14:1    & $D_\mathrm{6d}$ & $D_\mathrm{2d}$ & $-2.134$ & $1.173$  & $0$ & $12$ & $2$      & $0$ & 12  & $2.739$ & $3.048$ & 24  & 6.1  & 23.4  \\
15:1    & $D_\mathrm{3h}$ & $C_\mathrm{2v}$ & $-2.192$ & $-0.083$ & $0$ & $12$ & $3$      & $0$ & 12  & $2.786$ & $2.901$ & 21  & 5.1  & 29.2  \\
16:1    & $D_\mathrm{2}$  & $D_\mathrm{2 }$ & $-2.247$ & $0.223$  & $0$ & $12$ & $4$      & $0$ & 12  & $2.770$ & $2.917$ & 20  & 7.9  & 24.3  \\
16:2    & $T_\mathrm{d}$  & $D_\mathrm{2d}$ & $-2.233$ & $0.440$  & $0$ & $12$ & $4$      & $0$ & 12  & $2.716$ & $2.996$ & 18  & 1.3  & 28.5  \\
17:1    & $D_\mathrm{5h}$ & $C_\mathrm{s}$    & $-2.259$ & $0.177$  & $2$ & $8$  & $3$      & $3$ & 12  & $2.747$ & $3.026$ & 20  & 11.5 & 17.3  \\
17:2    & $C_\mathrm{2v}$ & $C_\mathrm{2v}$ & $-2.272$ & $-0.038$ & $0$ & $12$ & $5$      & $0$ & 12  & $2.769$ & $2.931$ & 18  & 7.6  & 19.1  \\
17:3    & $C_\mathrm{2v}$ & $C_\mathrm{2v}$ & $-2.277$ & $-0.128$ & $0$ & $12$ & $5$      & $0$ & 12  & $2.762$ & $3.139$ & 17  & 5.5  & 20.8  \\
18:1    & $C_\mathrm{2}$  & $C_\mathrm{2 }$ & $-2.307$ & $0.321$  & $0$ & $12$ & $6$      & $0$ & 12  & $2.736$ & $2.934$ & 17  & 9.2  & 16.9  \\
18:2    & $D_\mathrm{2}$  & $D_\mathrm{2 }$ & $-2.290$ & $0.627$  & $0$ & $12$ & $6$      & $0$ & 12  & $2.733$ & $2.935$ & 18  & 11.6 & 17.2  \\
18:3    & $D_\mathrm{3d}$ & $D_\mathrm{3d}$ & $-2.275$ & $0.896$  & $0$ & $12$ & $6$      & $0$ & 12  & $2.714$ & $2.894$ & 18  & 12.1 & 18.2  \\
18:4    & $C_\mathrm{2}$  & $C_\mathrm{2 }$ & $-2.321$ & $0.073$  & $0$ & $12$ & $6$      & $0$ & 12  & $2.749$ & $2.931$ & 16  & 7.2  & 18.7  \\
18:5    & $D_\mathrm{3h}$ & $D_\mathrm{3h}$ & $-2.303$ & $0.386$  & $0$ & $12$ & $6$      & $0$ & 12  & $2.763$ & $3.159$ & 8   & 15.1 & 27.3  \\
18:6    & $D_\mathrm{3}$  & $D_\mathrm{3 }$ & $-2.310$ & $0.270$  & $0$ & $12$ & $6$      & $0$ & 12  & $2.742$ & $2.945$ & 15  & 5.8  & 15.2  \\
19:1    & $C_\mathrm{2}$  & $C_\mathrm{2 }$ & $-2.298$ & $1.196$  & $0$ & $12$ & $7$      & $0$ & 12  & $2.745$ & $3.006$ & 17  & 14.9 & 26.0  \\
19:2    & $C_\mathrm{s}$  & $C_\mathrm{s }$ & $-2.307$ & $1.014$  & $0$ & $12$ & $7$      & $0$ & 12  & $2.747$ & $2.972$ & 15  & 7.5  & 20.0  \\
19:3    & $C_\mathrm{s}$  & $C_\mathrm{s }$ & $-2.304$ & $1.077$  & $0$ & $12$ & $7$      & $0$ & 12  & $2.737$ & $2.957$ & 15  & 11.9 & 28.3  \\
19:4    & $C_\mathrm{2}$  & $C_\mathrm{2 }$ & $-2.311$ & $0.935$  & $0$ & $12$ & $7$      & $0$ & 12  & $2.745$ & $2.905$ & 15  & 7.0  & 17.7  \\
19:5    & $C_\mathrm{2}$  & $C_\mathrm{2 }$ & $-2.313$ & $0.911$  & $0$ & $12$ & $7$      & $0$ & 12  & $2.734$ & $2.947$ & 14  & 6.6  & 18.7  \\
19:6    & $C_\mathrm{3v}$ & $C_\mathrm{3v}$ & $-2.316$ & $0.854$  & $0$ & $12$ & $7$      & $0$ & 12  & $2.765$ & $2.890$ & 15  & 12.7 & 30.6  \\
20:1    & $C_\mathrm{2}$  & $C_\mathrm{1 }$ & $-2.324$ & $1.684$  & $2$ & $8$  & $10$     & $0$ & 12  & $2.711$ & $2.984$ & 16  & 15.3 & 36.7  \\
20:2    & $D_\mathrm{2}$  & $D_\mathrm{2 }$ & $-2.295$ & $2.271$  & $0$ & $12$ & $8$      & $0$ & 12  & $2.699$ & $3.023$ & 18  & 20.4 & 22.0  \\
20:3    & $C_\mathrm{1}$  & $C_\mathrm{1 }$ & $-2.339$ & $1.395$  & $2$ & $8$  & $10$     & $0$ & 12  & $2.724$ & $2.954$ & 15  & 13.1 & 129.1 \\
20:4    & $C_\mathrm{s}$  & $C_\mathrm{s }$ & $-2.324$ & $1.695$  & $0$ & $12$ & $8$      & $0$ & 12  & $2.709$ & $3.023$ & 16  & 13.7 & 25.5  \\
20:5    & $D_\mathrm{2}$  & $D_\mathrm{2 }$ & $-2.332$ & $1.541$  & $0$ & $12$ & $8$      & $0$ & 12  & $2.749$ & $3.080$ & 16  & 18.5 & 17.3  \\
20:6    & $D_\mathrm{2d}$ & $C_\mathrm{2v}$ & $-2.337$ & $1.440$  & $2$ & $8$  & $10$     & $0$ & 12  & $2.752$ & $2.977$ & 14  & 9.8  & 26.6  \\
20:7    & $C_\mathrm{1}$  & $C_\mathrm{1 }$ & $-2.325$ & $1.663$  & $2$ & $9$  & $8$      & $1$ & 12  & $2.712$ & $3.019$ & 14  & 10.9 & 25.8  \\
20:8    & $C_\mathrm{s}$  & $C_\mathrm{s }$ & $-2.346$ & $1.256$  & $2$ & $8$  & $10$     & $0$ & 12  & $2.748$ & $3.057$ & 14  & 8.4  & 40.7  \\
20:9    & $C_\mathrm{2v}$ & $D_\mathrm{6h}$ & $-2.362$ & $0.938$  & $6$ & $0$  & $14$     & $0$ & 12  & $2.744$ & $2.971$ & 13  & 3.8  & 23.2  \\
20:10   & $C_\mathrm{2}$  & $C_\mathrm{2 }$ & $-2.344$ & $1.299$  & $2$ & $8$  & $10$     & $0$ & 12  & $2.726$ & $3.004$ & 14  & 12.6 & 23.9  \\
20:11   & $C_\mathrm{2}$  & $C_\mathrm{s }$ & $-2.346$ & $1.256$  & $2$ & $8$  & $10$     & $0$ & 12  & $2.747$ & $3.056$ & 13  & 8.1  & 35.8  \\
20:12   & $C_\mathrm{2}$  & $C_\mathrm{1}$    & $-2.366$ & $0.861$  & $3$ & $5$  & $4$      & $7$ & 4   & $2.719$ & $3.048$ & 13  & 5.4  & 21.1  \\
20:13   & $D_\mathrm{3h}$ & $D_\mathrm{6h}$ & $-2.362$ & $0.938$  & $6$ & $0$  & $14$     & $0$ & 12  & $2.744$ & $2.970$ & 15  & 6.5  & 27.9  \\
20:14   & $D_\mathrm{2d}$ & $D_\mathrm{2d}$ & $-2.311$ & $1.948$  & $0$ & $12$ & $8$      & $0$ & 12  & $2.779$ & $2.929$ & 12  & 3.7  & 22.0  \\
20:15   & $D_\mathrm{6h}$ & $D_\mathrm{6h}$ & $-2.362$ & $0.936$  & $6$ & $0$  & $14$     & $0$ & 12  & $2.744$ & $2.972$ & 12  & 4.5  & 25.6  \\
32:1082 & $I_\mathrm{h}$  & $I_\mathrm{h}$    & $-2.494$ & $1.537$  & $0$ & $12$ & $20$     & $0$ & 12  & $2.793$ & $2.835$ & 0   & 0    & 7.5   \\
(111)   & 2D              & sheet    & $-2.994$ & $-$      & $0$ & $0$  & $\infty$ & $0$ & $-$ & $2.722$ & $2.722$ & 0   & 0    & 0     \\
fcc     & 3D              & bulk     & $-3.677$ & $-$      & $-$ & $-$  & $-$      & $-$ & $-$ & $2.897$ & $2.897$ & $-$ & $-$  & $-$   \\
		\bottomrule
    \end{tabular}}
\end{table}
%
%
\begin{table}[htbp]
	\centering
    \setlength{\tabcolsep}{1.5pt}
    \footnotesize{
    \caption{Topological parameters for the anionic gold clusters. Number of
    gold atoms and isomer numbers of the fullerene in canonical order of the
    pentagon spiral indices,\autocite{Fowler-atlas-2006} ideal and actual point
    group symmetry, energy differences $\Delta E_g$ to the most stable anionic
    cluster of same size and binding energy per atom $\Delta E_n =
    [E(\textrm{Au}_n)-(n-1)E(\textrm{Au})-E(\textrm{Au}^-)]/n$  (in eV),
    shortest and largest bond distance (in \AA), and distortion parameter $D$
    (in \%) for the \acs{GDF}. }
	\label{tab:anion}
\begin{tabular}{lllrrrrrrrrrr}
\toprule
\multicolumn{1}{c}{  } & \multicolumn{2}{c}{ symmetry  }  & \multicolumn{2}{c}{stability} & \multicolumn{4}{c}{ vertices } & & \multicolumn{2}{c}{ bondlengths } &  \multicolumn{1}{c}{$D$} \\
isomer & ideal  & actual  & $\Delta E_n$ &$\Delta E_g$ & \multicolumn{1}{c}{$|N_4|$} & \multicolumn{1}{c}{$|N_5|$} & \multicolumn{1}{c}{$|N_6|$} & \multicolumn{1}{c}{$|N_7|$}& $\Gamma$ & shortest & largest   & GDF  \\\midrule
    12:1    & $I_\mathrm{h}$  & $D_\mathrm{2d}$ & $-2.137$ & $0.665$  & $8$ & $0$  & $4$  & $0$ & 16 & $2.780$ & $2.869$ & 23.0  \\
    14:1    & $D_\mathrm{6d}$ & $D_\mathrm{2d}$ & $-2.242$ & $-0.089$ & $0$ & $12$ & $2$  & $0$ & 12 & $2.758$ & $2.989$ & 20.3  \\
    15:1    & $D_\mathrm{3h}$ & $C_\mathrm{2v}$ & $-2.281$ & $0.473$  & $0$ & $12$ & $3$  & $0$ & 12 & $2.741$ & $3.029$ & 21.2  \\
    16:1    & $D_\mathrm{2}$  & $D_\mathrm{2 }$ & $-2.328$ & $0.020$  & $0$ & $12$ & $4$  & $0$ & 12 & $2.764$ & $2.905$ & 17.7  \\
    16:2    & $T_\mathrm{d}$  & $D_\mathrm{2d}$ & $-2.330$ & $0.000$  & $0$ & $12$ & $4$  & $0$ & 12 & $2.738$ & $2.907$ & 16.2  \\
    17:1    & $D_\mathrm{5h}$ & $D_\mathrm{5h}$ & $-2.353$ & $0.469$  & $0$ & $12$ & $5$  & $0$ & 12 & $2.757$ & $3.017$ & 13.2  \\
    17:2    & $C_\mathrm{2v}$ & $C_\mathrm{2v}$ & $-2.368$ & $0.215$  & $0$ & $12$ & $5$  & $0$ & 12 & $2.742$ & $2.994$ & 14.4  \\
    17:3    & $C_\mathrm{2v}$ & $C_\mathrm{2v}$ & $-2.376$ & $0.087$  & $0$ & $12$ & $5$  & $0$ & 12 & $2.731$ & $3.019$ & 14.2  \\
    18:1    & $C_\mathrm{2}$  & $C_\mathrm{2 }$ & $-2.360$ & $0.589$  & $0$ & $12$ & $6$  & $0$ & 12 & $2.734$ & $2.968$ & 16.8  \\
    18:2    & $D_\mathrm{2}$  & $C_\mathrm{2 }$ & $-2.346$ & $0.848$  & $0$ & $12$ & $6$  & $0$ & 12 & $2.733$ & $3.059$ & 16.8  \\
    18:3    & $D_\mathrm{3d}$ & $C_\mathrm{2 }$ & $-2.348$ & $0.817$  & $4$ & $4$  & $10$ & $0$ & 12 & $2.701$ & $3.038$ & 24.3  \\
    18:4    & $C_\mathrm{2}$  & $C_\mathrm{1 }$ & $-2.364$ & $0.529$  & $0$ & $12$ & $6$  & $0$ & 12 & $2.740$ & $3.048$ & 19.0  \\
    18:5    & $D_\mathrm{3h}$ & $D_\mathrm{3h}$ & $-2.364$ & $0.516$  & $0$ & $12$ & $6$  & $0$ & 12 & $2.710$ & $3.023$ & 27.5  \\
    18:6    & $D_\mathrm{3}$  & $D_\mathrm{3 }$ & $-2.357$ & $0.642$  & $0$ & $12$ & $6$  & $0$ & 12 & $2.734$ & $2.912$ & 14.7  \\
    19:1    & $C_\mathrm{2}$  & $C_\mathrm{2 }$ & $-2.384$ & $0.967$  & $4$ & $4$  & $11$ & $0$ & 12 & $2.732$ & $2.985$ & 28.5  \\
    19:2    & $C_\mathrm{s}$  & $C_\mathrm{s }$ & $-2.368$ & $1.268$  & $2$ & $9$  & $7$  & $1$ & 12 & $2.727$ & $2.989$ & 16.9  \\
    19:3    & $C_\mathrm{s}$  & $C_\mathrm{3v}$ & $-2.381$ & $1.022$  & $0$ & $12$ & $7$  & $0$ & 12 & $2.755$ & $3.046$ & 32.6  \\
    19:4    & $C_\mathrm{2}$  & $C_\mathrm{2 }$ & $-2.390$ & $0.853$  & $2$ & $8$  & $9$  & $0$ & 12 & $2.744$ & $3.003$ & 22.1  \\
    19:5    & $C_\mathrm{2}$  & $C_\mathrm{2 }$ & $-2.398$ & $0.698$  & $2$ & $8$  & $9$  & $0$ & 12 & $2.743$ & $2.963$ & 31.3  \\
    19:6    & $C_\mathrm{3v}$ & $C_\mathrm{3v}$ & $-2.381$ & $1.023$  & $0$ & $12$ & $7$  & $0$ & 12 & $2.756$ & $3.044$ & 32.4  \\
    20:1    & $C_\mathrm{2}$  & $C_\mathrm{1 }$ & $-2.386$ & $0.927$  & $2$ & $8$  & $10$ & $0$ & 12 & $2.748$ & $3.032$ & 25.0  \\
    20:2    & $D_\mathrm{2}$  & $D_\mathrm{2 }$ & $-2.365$ & $1.348$  & $0$ & $12$ & $8$  & $0$ & 12 & $2.731$ & $2.971$ & 43.8  \\
    20:3    & $C_\mathrm{1}$  & $C_\mathrm{1 }$ & $-2.396$ & $0.716$  & $2$ & $8$  & $10$ & $0$ & 12 & $2.745$ & $2.926$ & 23.4  \\
    20:4    & $C_\mathrm{s}$  & $C_\mathrm{s }$ & $-2.385$ & $0.950$  & $0$ & $12$ & $8$  & $0$ & 12 & $2.740$ & $2.939$ & 24.4  \\
    20:5    & $D_\mathrm{2}$  & $D_\mathrm{2 }$ & $-2.390$ & $0.850$  & $0$ & $12$ & $8$  & $0$ & 12 & $2.775$ & $2.907$ & 37.4  \\
    20:6    & $D_\mathrm{2d}$ & $C_\mathrm{s }$ & $-2.382$ & $1.001$  & $2$ & $8$  & $10$ & $0$ & 12 & $2.770$ & $2.948$ & 20.5  \\
    20:7    & $C_\mathrm{1}$  & $C_\mathrm{1 }$ & $-2.384$ & $0.974$  & $0$ & $12$ & $8$  & $0$ & 12 & $2.739$ & $3.079$ & 25.3  \\
    20:8    & $C_\mathrm{s}$  & $C_\mathrm{s }$ & $-2.388$ & $0.888$  & $2$ & $8$  & $10$ & $0$ & 12 & $2.761$ & $2.977$ & 22.0  \\
    20:9    & $C_\mathrm{2v}$ & $D_\mathrm{6h}$ & $-2.402$ & $0.610$  & $6$ & $0$  & $14$ & $0$ & 12 & $2.731$ & $2.971$ & 27.1  \\
    20:10   & $C_\mathrm{2}$  & $C_\mathrm{2 }$ & $-2.404$ & $0.568$  & $2$ & $8$  & $10$ & $0$ & 12 & $2.743$ & $2.984$ & 37.1  \\
    20:11   & $C_\mathrm{2}$  & $C_\mathrm{1}$  & $-2.378$ & $1.093$  & $3$ & $8$  & $7$  & $2$ & 12 & $2.732$ & $3.018$ & 26.5  \\
    20:12   & $C_\mathrm{2}$  & $C_\mathrm{s}$  & $-2.412$ & $0.407$  & $2$ & $8$  & $3$  & $6$ & 6  & $1.755$ & $2.996$ & 192.0 \\
    20:13   & $D_\mathrm{3h}$ & $D_\mathrm{6h}$ & $-2.402$ & $0.610$  & $6$ & $0$  & $14$ & $0$ & 12 & $2.731$ & $2.972$ & 27.0  \\
    20:14   & $D_\mathrm{2d}$ & $C_\mathrm{1}$  & $-2.405$ & $0.544$  & $2$ & $8$  & $1$  & $8$ & 12 & $2.710$ & $3.010$ & 199.2 \\
    20:15   & $D_\mathrm{6h}$ & $D_\mathrm{6h}$ & $-2.361$ & $1.430$  & $0$ & $12$ & $8$  & $0$ & 12 & $2.792$ & $2.933$ & 15.2  \\
    32:1082 & $I_\mathrm{h}$  & $D_\mathrm{2h}$ & $-2.524$ & $2.201$  & $0$ & $12$ & $20$ & $0$ & 12 & $2.766$ & $3.004$ & 10.4  \\
		\bottomrule
\end{tabular}}
\end{table}

Tables~\ref{tab:neutral} and \ref{tab:anion} show vertex counts as well as
results from equation~\eqref{eq3valent} for the neutral and anionic clusters,
respectively. Considering only the topological parameter $\Gamma$ it is clear
that most of the optimised structures can be derived from a dual planar cubic
graph and therefore only consist of triangles. The few notable exceptions are
the isomers 12:1 and 20:12 for both the anionic and neutral structure. The ideal
icosahedral structure for the Au$_{12}$ cluster is not stable under the present
level of theory, and the optimised structure does not correspond to a
triangulation of a sphere. However, it has already been shown that this cage can
be stabilised by inserting a transition metal (e.g. tungsten) atom into the
central position of the icosahedron such that the 18 valence electron rule is
fulfilled.\autocite{Pyykko_IcosahedralWAu12Predicted_2002,Autschbach_PropertiesWAu12_2004}
Additional stabilisation of such an endohedral gold cluster can be achieved by
attaching ligands to the surface of the cluster.\autocite{Laupp-1994} Structure
20:12 converges towards a more compact cluster with an 8-fold coordinated gold
atom in the centre for both the anionic and the neutral cluster. 
%
\begin{figure}[htbp]
    \centering
	\begin{tabular}{lp{3cm}|p{3cm}}\toprule
		isomer  & \multicolumn{1}{c}{neutral} & \multicolumn{1}{c}{anion} \\ \midrule
		12:1    & \cellcolor{myorange}        & \cellcolor{myorange}      \\
		14:1    & \cellcolor{mygreen}         & \cellcolor{mygreen}       \\
		15:1    & \cellcolor{mygreen}         & \cellcolor{mygreen}       \\
		16:1    & \cellcolor{mygreen}         & \cellcolor{mygreen}       \\
		16:2    & \cellcolor{mygreen}         & \cellcolor{mygreen}       \\
		17:1    & \cellcolor{myorange}        & \cellcolor{mygreen}       \\
		17:2    & \cellcolor{mygreen}         & \cellcolor{mygreen}       \\
		17:3    & \cellcolor{mygreen}         & \cellcolor{mygreen}       \\
		18:1    & \cellcolor{mygreen}         & \cellcolor{mygreen}       \\
		18:2    & \cellcolor{mygreen}         & \cellcolor{mygreen}       \\
		18:3    & \cellcolor{mygreen}         & \cellcolor{myorange}      \\
		18:4    & \cellcolor{mygreen}         & \cellcolor{mygreen}       \\
		18:5    & \cellcolor{mygreen}         & \cellcolor{mygreen}       \\
		18:6    & \cellcolor{mygreen}         & \cellcolor{mygreen}       \\
		19:1    & \cellcolor{mygreen}         & \cellcolor{myorange}      \\
		19:2    & \cellcolor{mygreen}         & \cellcolor{myorange}      \\
		19:3    & \cellcolor{mygreen}         & \cellcolor{mygreen}       \\
		19:4    & \cellcolor{mygreen}         & \cellcolor{myorange}      \\
		19:5    & \cellcolor{mygreen}         & \cellcolor{myorange}      \\
		19:6    & \cellcolor{mygreen}         & \cellcolor{mygreen}       \\
		20:1    & \cellcolor{myorange}        & \cellcolor{myorange}      \\
		20:2    & \cellcolor{mygreen}         & \cellcolor{mygreen}       \\
		20:3    & \cellcolor{myorange}        & \cellcolor{myorange}      \\
		20:4    & \cellcolor{mygreen}         & \cellcolor{mygreen}       \\
		20:5    & \cellcolor{mygreen}         & \cellcolor{mygreen}       \\
		20:6    & \cellcolor{myorange}        & \cellcolor{myorange}      \\
		20:7    & \cellcolor{myorange}        & \cellcolor{mygreen}       \\
		20:8    & \cellcolor{myorange}        & \cellcolor{myorange}      \\
		20:9    & \cellcolor{myorange}        & \cellcolor{myorange}      \\
		20:10   & \cellcolor{myorange}        & \cellcolor{myorange}      \\
		20:11   & \cellcolor{myorange}        & \cellcolor{myorange}      \\
		20:12   & \cellcolor{red}             & \cellcolor{red}           \\
		20:13   & \cellcolor{myorange}        & \cellcolor{myorange}      \\
		20:14   & \cellcolor{mygreen}         & \cellcolor{red}           \\
		20:15   & \cellcolor{myorange}        & \cellcolor{mygreen}       \\
		32:1082 & \cellcolor{mygreen}         & \cellcolor{mygreen}       \\ \bottomrule
	\end{tabular}
	\caption{Overview of PBE-D3 optimisation results for the dual fullerene
	structures. Green: dual fullerene structure, orange: hollow structure, red:
	non-hollow structure.}
	\label{fig:optOverview}
\end{figure}

Figure~\ref{fig:optOverview} gives an overview over all optimised structures.
A green field marks a dual fullerene structure with exactly 12 vertices of
degree five and the remaining vertices being of degree six. These are also the
structures used in figure~\ref{fig:cohesiveenergies2} and they are more
abundant for clusters of size 14 to 19 atoms. Structures with an orange mark do
not fulfil the requirement of being a dual fullerene as they contain vertices
of degree 4. However, they are still hollow gold cages and, as mentioned
before, show a value of $\Gamma=12$. These structures can be rather similar to
the initial dual fullerene structures obtained from a force-field optimization
of the corresponding carbon cage, and are usually a result of a flattening
towards a more oblate geometry.  Most of the clusters shown here preserve their
hollow cage structure with only few clusters optimising into more stable
compact structures. These are marked as red in figure~\ref{fig:optOverview}.

As illustrated by the distortion parameter $D$(F) in tables~\ref{tab:neutral}
and \ref{tab:anion}, carbon fullerenes try to adopt ``spherical'' shapes if
permitted by the distribution of pentagons. This is especially the case for
$I_\mathrm{h}$-C$_{20}$ and $I_\mathrm{h}$-C$_{60}$ with a distortion parameter
of exactly zero (i.e. all atoms lie on a sphere).  In contrast, the golden dual
fullerene structures have much larger distortion parameters $D$(GDF) than their
carbon equivalent and are therefore less spheroidal. The golden dual fullerenes
usually distort into less symmetric structures, for example into oblate
structures as mentioned above.

Figure~\ref{fig:AunMinus2} shows the relative energies $\Delta E_g$ per atom
compared to the most stable compact arrangement for all optimised hollow gold
clusters.
%
\begin{figure}[htbp]
    \begin{center}
    \includegraphics[width=.8\textwidth]{golddual/energies.pdf}
    \caption{Relative energies for the investigated dual fullerene clusters.
    Energy differences compared to the most stable compact cluster (per atom) are given in eV.}
    \label{fig:AunMinus2}
    \end{center}
\end{figure}
%
It is immediately apparent, that the most stable dual fullerene structures can
be found in the region of 14 to 18 atoms. Some clusters in this region even
exceed the stability of formerly proposed global minimum structures. For
example, for \ce{Au16^-} the global minimum has been proposed previously to be
the tetrahedral hollow
cluster,\autocite{Schooss_Determiningsizedependentstructure_2010,Lechtken_Structuredeterminationgold_2009}
which is the dual of the tetrahedral C$_{28}$ isomer as observed experimentally
in photoelectron spectra.\autocite{Bulusu_Evidencehollowgolden_2006} It should
be noted, that Chen et al. have found the tetrahedral structure to lie 0.22~eV
above a sheet-like structure.\autocite{Chen_Structuresneutralanionic_2010}
However, our results contradict these findings as the planar structure is
predicted to be 0.939~eV higher in energy. Another interesting result from the
investigation of the cohesive energies is that the $D_\mathrm{2}$ symmetric
isomer 16:1 lies only 0.02~eV above the tetrahedral structure. Therefore, it
should also be possible to observe this isomer by experimental methods.

Possible Au$_{32}$ structures have been investigated intensively by Jalbout et
al.\autocite{Jalbout_LowSymmetryStructuresAu_2008} Table~\ref{tab:Aun} shows
their results in comparison with results from this work.
%
\begin{table}[htbp]
	\centering
    \setlength{\tabcolsep}{3pt}
    \footnotesize{
    \caption{Binding energy per atom (in eV) for investigated neutral and
    anionic compact cluster compounds. For the definition of the binding energy
    see tables~\ref{tab:neutral} and \ref{tab:anion}, and for the definition of
    the isomers 1 and 10 for Au$_{32}$ see \citeauthor{Jalbout_LowSymmetryStructuresAu_2008}.\autocite{Jalbout_LowSymmetryStructuresAu_2008}}
	\label{tab:Aun}
	\begin{tabular}{llcllcllc}
		\toprule
		$N$  & sym.  & $\Delta E_n$(neutral) & $N$  & sym.  & $\Delta E_n$(neutral) & $N$  & sym.  & $\Delta E_n$(anion) \\
		\midrule
2  & $D_\mathrm{\infty h}$ & $-1.105$ & 13 & $C_\mathrm{2v}$ & $-2.087$ & 12 & $D_\mathrm{3h}$  & $-2.192$ \\
3  & $C_\mathrm{2v}$       & $-1.152$ & 14 & $C_\mathrm{2v}$ & $-2.218$ & 14 & $D_\mathrm{2h}$  & $-2.236$ \\
4  & $D_\mathrm{2h}$       & $-1.486$ & 15 & $C_\mathrm{s}$  & $-2.186$ & 15 & $C_\mathrm{1}$     & $-2.313$ \\
5  & $C_\mathrm{2v}$       & $-1.631$ & 16 & $C_\mathrm{s}$  & $-2.261$ & 16 & $D_\mathrm{2d}$  & $-2.330$ \\
6  & $D_\mathrm{3h}$       & $-1.875$ & 17 & $C_\mathrm{s}$  & $-2.270$ & 17 & $C_\mathrm{2v}$  & $-2.381$ \\
7  & $C_\mathrm{s}$        & $-1.833$ & 18 & $C_\mathrm{s}$  & $-2.325$ & 18 & $C_\mathrm{2v}$  & $-2.393$ \\
8  & $D_\mathrm{4h}$       & $-1.959$ & 19 & $C_\mathrm{3v}$ & $-2.361$ & 19 & $C_\mathrm{3v}$  & $-2.435$ \\
9  & $C_\mathrm{2v}$       & $-1.944$ & 20 & $T_\mathrm{d}$  & $-2.409$ & 20 & $T_\mathrm{d}$     & $-2.432$ \\
10 & $D_\mathrm{2h}$       & $-2.028$ & 32 & $C_\mathrm{3v}$ & $-2.491$ & 32 & $C_\mathrm{3v}$  & $-2.548$ \\
11 & $D_\mathrm{3h}$       & $-2.063$ & 32 & Isomer 1        & $-2.536$ & 32 & Isomer 1  & $-2.590$ \\
12 & $D_\mathrm{3h}$       & $-2.098$ & 32 & Isomer 10       & $-2.542$ & 32 & Isomer 10 & $-2.593$ \\
		\bottomrule
    \end{tabular}}
\end{table}
%
For both neutral and anionic clusters, isomer 10 in their work turns out to be
the most stable compact geometry and the icosahedral hollow structure 32:1082 is
less stable in both the neutral and the anionic cases. The
$C_\mathrm{3v}$-symmetric compact structure not investigated before is also
included in table~\ref{tab:Aun}. It is derived from the ideal Au$_{35}$
tetrahedron by removing three of the corner atoms of the tetrahedron and can
be viewed as a cut-out of the \ac{fcc} bulk structure. This cluster is also very
stable compared to the other structures proposed by Jalbout et al. As reflected
by the distortion parameter $D$ of the \ce{Au32} hollow cage
($D(\ce{Au32})=10.4$) it deviates slightly from an ideal icosahedral symmetry
and can be seen as pseudo-spherical.
  
\section{Convergence Towards the Infinite Structure}

The neutral gold clusters and their property convergence towards the bulk has
already been discussed in previous
papers.\autocite{Assadollahzadeh_systematicsearchminimum_2009} Increasing the
size of non-hollow compact clusters lowers the cohesive energy until the
clusters are large enough to be a valid representation of the bulk gold
structure. This can be seen in figure~\ref{fig:cohesiveenergies1}, where a
clear linear correlation between $N^{-1/3}$ and the cohesive energy is
depicted.
%
\begin{figure}[htb]\centering
    \subfloat[\label{fig:cohesiveenergies1}]{\includegraphics[width=.49\textwidth]{golddual/cohesive.pdf}}\hfill
    \subfloat[\label{fig:cohesiveenergies2}]{\includegraphics[width=.49\textwidth]{golddual/cohesive2.pdf}}
	\caption{Cohesive energies for \protect\subref{fig:cohesiveenergies1} the compact gold clusters with cluster size $N$ and convergence toward the bulk \acs{fcc} structure and \protect\subref{fig:cohesiveenergies2} for the hollow gold clusters with cluster size $N$ and convergence toward the (111) gold sheet.}
\end{figure}
%
Hollow gold clusters can be created by wrapping a cut-out from a (111) gold 2D
sheet around a sphere while introducing 12 vertices of degree 5 to satisfy
Euler's theorem. Therefore, an infinitely large 2D gold sheet represents a
golden dual fullerene cage with an infinite sphere radius. As the cohesive
energy of the compact structures converges towards the bulk cohesive energy, the
cohesive energy of the 2D triangulated gold sheet should represent the infinite limit for the
dual golden fullerene structures. This is indeed the case and is depicted in
figure~\ref{fig:cohesiveenergies2} using a $N^{-1}$ scaling law analogous to the
one used for fullerenes.\autocite{Wirz_smallfullerenesgraphene_2015}

An interesting result was the difference between the cohesive energy of the bulk
\ac{fcc} structure compared to the (111) 2D sheet. Creating the bulk structure
from stacking (111) sheets only accounts for $\sim$0.68~eV of the total cohesive
energy of the bulk which is 3.81 eV.\autocite{takeuchi_first-principles_1989}
This implies that most of the cohesive energy of bulk gold originates from the
(111) sheet, which is therefore exceptionally stable and can be seen as a reason
for the preferred planar arrangement of many small gold clusters. As pointed out
by Takeuchi et al, relativistic effects increase the cohesive energy of bulk
gold by 1.5 eV.\autocite{takeuchi_first-principles_1989} A similar large
relativistic effect is expected for the (111) sheet of gold.

\section{Simulation of Photoelectron Spectra}

Photoelectron spectra of several \acp{GDF} have been determined experimentally
and simulated with theoretical methods by Bulusu et
al.\autocite{Bulusu_Evidencehollowgolden_2006} Before the discussion of results
produced in this work can commence, the spin-orbit effects from substantial
5d-mixing into the 6s orbitals in gold need to be considered.
Figure~\ref{fig:photoSOAu17} shows a comparison of simulated photoelectron
spectra of the three golden dual fullerene isomers of \ce{Au17^-}.
%
\begin{figure}[htb]\centering
	\subfloat[]{\includegraphics[width=.49\textwidth]{golddual/photo/Au17/turbo/rel-nonrel/photo1.pdf}}\hfill
	\subfloat[]{\includegraphics[width=.49\textwidth]{golddual/photo/Au17/turbo/rel-nonrel/photo2.pdf}}\\
	\subfloat[]{\includegraphics[width=.49\textwidth]{golddual/photo/Au17/turbo/rel-nonrel/photo3.pdf}}
	\caption{Comparison of simulated photoelectron spectra of the three dual fullerene isomers of \ce{Au17^-} with (2c) and without spin-orbit coupling.} 
	\label{fig:photoSOAu17}
\end{figure}
%
The results clearly indicate that spin-orbit effects can be safely neglected in this
energy range.

Bulusu et al. considered only the $T_\text{d}$-\ce{Au16} structure. From the
simulations carried out in this section there is reason to believe that the
other possible isomer 16:1 is also present in the measured spectrum.
Figure~\ref{subfig:photo1} shows these simulation results for the isomers 16:1
and 16:2 and a simulation for a mixture of both compounds with a ratio of 1:1 as
the energy of both isomers is comparable. 
%
\begin{figure}[htb]
    \begin{center}
        \subfloat[\label{subfig:photo1}]{\includegraphics[width=.49\textwidth]{golddual/photo/Au16/turbo/photo1.pdf}}\hfill
        \subfloat[\label{subfig:photo2}]{\includegraphics[width=.49\textwidth]{golddual/photo/Au17/turbo/compare/compare.pdf}}\\
        \subfloat[\label{subfig:photo3}]{\includegraphics[width=.49\textwidth]{golddual/photo/Au18/turbo/nonrel/compare.pdf}}
    \caption{Simulated photoelectron spectra for the negatively charged hollow
        gold clusters (shifted to the experimental threshold energy).
        \protect\subref{subfig:photo1} The two possible dual fullerene isomers
        of \ce{Au16^-}. The green curve shows a combination of the
        $D_\mathrm{2}$ and $T_\mathrm{d}$ spectra with a ratio of 1:1.;
        \protect\subref{subfig:photo2} The three possible dual fullerene isomers
        of \ce{Au17^-}; \protect\subref{subfig:photo3} The six possible dual
        fullerene isomers of
        \ce{Au18^-} (shifted to the experimental threshold energy).}
      \label{fig:photo_Au16}
    \end{center}
    \end{figure}
%
Looking at the experimental data, a shoulder can be identified in the first
peak. This feature can be reproduced by shifting the spectra for 16:1 and 16:2
according to the corresponding vertical ionisation potential, superimposing both
spectra and shifting the result by 0.18~eV to better fit the experimental data
as pictured in figure~\ref{subfig:photo1}. This indicates that the second hollow
cage isomer has also been produced. Further evidence for this could be the
experimental peak at 5.51~eV. The simulated spectrum for the tetrahedral cluster
shows a dip at this energy, while the $D_\mathrm{2}$ structure has a clear
intensity maximum. 

Figure~\ref{subfig:photo2} shows the simulated spectra for the three possible
dual fullerene isomers for \ce{Au17^-}. The spectra have been shifted according
to the vertical ionisation potential of the negatively charged clusters. The
most stable structures are 17:2 (red) and 17:3 (green) of which 17:2 fits
reasonably well for the first 4 peaks. The peak at 4.73~eV could be accounted
to the 17:3 isomer identical to Bulusu et al.'s $C_\mathrm{2v}$
symmetric structure.

From the relative energies in figure~\ref{fig:AunMinus2} it is clear that
anionic dual fullerene structures start to become rather unstable for $N=18$.
Therefore, compact clusters might dominate the experimental spectrum.
Figure~\ref{subfig:photo3} shows our calculated spectra in comparison with the
experimental data. The calculated spectra have been shifted to the corresponding
vertical ionisation potential first and subsequently shifted by 0.25~eV to
better fit the experimental data. The most stable dual fullerene clusters are
18:1, 18:4 and 18:5. 18:1 and 18:5 could be responsible for the second peak in
the experimental data at 3.63~eV, while 18:4 agrees with the first peak. The
signal at 3.97~eV could be an indication that isomer 18:3 was produced as it is
the only structure that shows a peak in that area, however, it is the least
stable of the hollow structures.

Finally, for future experiments a simulated photoelectron spectrum for
Au$_{32}^-$ is shown in figure~\ref{fig:photo_Au32}.
%
\begin{figure}[htb]
\begin{center}
\includegraphics[width=.75\textwidth]{golddual/photo/Au32/nonrel/compare.pdf}
\caption{Simulated photoelectron spectra for isomer 32:1812 of Au$_{32}^-$.}
  \label{fig:photo_Au32}
\end{center}
\end{figure}

\section{Conclusion}

An interesting topological relationship between fullerenes and the cage-like
gold clusters resulting in a triangulation of a sphere with vertices of degrees
5 and 6 fulfilling Euler's polyhedral formula was found. Because of this
isomorphism between the two types of structures by dualisation, there are as
many golden fullerene isomers as there are fullerene isomers. Gold nano-tubes
and carbon nano-tubes and halma transforms of C$_{20}$ to the shells of a Mackay
icosahedron are related in the same way. The stability of these golden
fullerenes was investigated. While they perhaps may not compete in energy with
the more compact gold clusters at larger cluster size, the smaller cage
structures are stable as observed by photoelectron spectroscopy. The simulated
photoelectron spectra suggest that more than one golden fullerene isomer was
observed. 

A natural step in the next direction would be to stabilise such hollow gold
clusters by either endohedral enclosure of gold or other metal atoms, by
attaching appropriate ligands to the outside of the cage, or both. Results have
already shown that enclosing a transition metal with the right amount of valence
electrons can stabilise otherwise instable hollow structures, like the
icosahedral \ce{Au12},
drastically.\autocite{Pyykko_IcosahedralWAu12Predicted_2002}

\acresetall

\chapter[From Sticky-Hard-Sphere to Lennard-Jones-Type Clusters]{From
    Sticky-Hard-Sphere to Lennard-Jones-Type clusters\footnote{This chapter is
    composed of sections previously published in the article
    \citetitle*{Trombach_stickyhardsphereLennardJonestypeclusters_2018}\autocite{Trombach_stickyhardsphereLennardJonestypeclusters_2018}
    and is reprinted with permission from the publisher \textcopyright 2018
    American Physical Society. Some sections have been modified to fit the style
    of this thesis.}}
\label{sec:fromstickyhardspheretoLJtypeclusters}

%Introduction
\section{Introduction}

Nucleation is a phenomenon that is a part of many natural processes and is
present in many everyday phenomena. Naturally, there is a large research
interest in this field, especially with respect to the nucleation of atoms and
molecules to clusters, eventually leading to the solid
state.\autocite{Stillinger_Packingstructurestransitions_1984,
Martin-1996,Wales-1996, Vlieg_atomicscaleunderstandingcrystal_2007, Arkus-2010,
Woodley-2010, Karthika-2016, Holmes-Cerfon_StickySphereClusters_2017} In an
early Faraday Discussion taking place in Bristol in 1949 Rowland concluded the
meeting with the assessment ``that the gap between the theoretical and
experimental approaches has been too wide''.\autocite{Rowland-1949} Here, he was
referring to the subject of nucleation. In a Faraday discussion about half a
century later \citeauthor{Vlieg_atomicscaleunderstandingcrystal_2007} stated
that ``the gap between the quite detailed experimental information [...] and
theoretical models, though getting smaller, is still
large''.\autocite{Vlieg_atomicscaleunderstandingcrystal_2007}

One reason for this slow progress in the theoretical description of nucleation
processes is that it is related to global optimisation problems. Exploring the
multi-dimensional potential energy surface belongs to the computational
complexity class ``NP-complete'', as already mentioned in
chapter~\ref{sec:GlobalOptimisation}. The number of local minima is expected to
grow exponentially,\autocite{Stillinger_Packingstructurestransitions_1984,
Oganov-2006, Massen_Powerlawdistributionsareas_2007, wales10, Oganov-2011,
calvo12, Wales-2015} which is problematic, because interesting phase transitions
usually occur for larger cluster sizes $N$ that can not be treated by accurate
quantum mechanical methods as introduced in chapter~\ref{sec:basicsofQC}. One
such phase transition is the transformation of argon clusters from icosahedral
clusters into anti-Mackay clusters for sizes of $N>2000$ and finally into
\acf{fcc} or \acf{hcp} solid state structures for
$N>10^5$.\autocite{Krainyukova-2012} Similar results are predicted by the
\acf{LJ} potential.\autocite{Martin-1996,
Schwerdtfeger_ExtensionLennardJonespotential_2006, Krainyukova-2007}

Because of the exponentially growing potential energy landscapes and the large
cluster sizes required to model phase transitions, investigations of this type
often have to rely on approximate interaction potentials. In this part of the
thesis two interaction potentials introduced in
chapter~\ref{sec:energylandscapes} are used. The first and maybe simpler one is
the \acf{SHS} potential $V_\mathrm{SHS}$, originally introduced by
\citeauthor{baxter68}\autocite{baxter68}.
%
\begin{align}
    V_\mathrm{SHS}(r)=\begin{cases}
        \infty, & r < r_s\\
        -\varepsilon, & r = r_s\\
        0, & r > r_s
    \end{cases}
\label{eqn:KS}
\end{align}
%
$\varepsilon$ and the equilibrium distance $r_s$ can be set to unity without
changing the qualitative information contained in the \ac{PES}. 

One of the defining properties for clusters bound by the \ac{SHS} potential is
the contact number $N_c$, which is directly related to the energy of the cluster
$E=-N_c\varepsilon$. How the number of contact points a cluster possesses grows
with its size $N$ is still researched actively. From the Gregory-Newton
argument,\footnote{For a more thorough discussion of this argument please refer
to chapter~\ref{sec:thegregorynewtonclusters}.} which was proven in 1953, it
follows that each sphere can not be surrounded by more than 12 other spheres of
the same size. For small clusters, the maximum number is governed by the total
number of (unique) entries in the adjacency matrix, i.e. $\text{max}\leq
N(N-1)/2$. Combining these two facts leads to a loose upper bound
%
\begin{equation}
    N_c^\mathrm{max}(N) \le \text{min}\{N(N-1)/2,f(N)\}.
    \label{eqn:upperlimitNc}
\end{equation}
%
Using the Gregory-Newton argument results in $f(N)=6N$, however, a tighter upper
bound has been published more recently by
\citeauthor{Bezdek-2013}\autocite{Bezdek-2013}
%
\begin{equation}
    f(N)=6N-3(18)^{1/3}\pi^{-2/3}N^{2/3}.
    \label{eqn:upperlimitBR}
\end{equation}
%
Theoretical investigations of the cluster landscape by means of the exact
enumeration method revealed that the contact numbers for clusters of size $4
\leq N \leq 19$ are\autocite{Hoy_Structuredynamicsmodel_2015,Holmes-Cerfon_EnumeratingRigidSphere_2016}
%
\begin{align}
    N_c^\mathrm{max}(N) =\{6,9,12,15,18,21,25,29,33,36,40,44,48,52,56,60\}.
\end{align}
%
The exact solution for $N_c^\mathrm{max}(N)$ for arbitrary $N$ is called the
Erd\H{o}s unit distance problem, which remains unsolved.\autocite{Erdos-1946}

The number of non-isomorphic cluster structures
$|\mathcal{M}(N)|$\footnote{$\mathcal{M}(N)$ refers to the set of all
non-isomorphic cluster structures, while the size of the set is denoted
$|\mathcal{M}(N)|$.} is expected to grow
exponentially.\autocite{Stillinger_Exponentialmultiplicityinherent_1999,Oganov-2006,Forman_ModelingAggregationProcesses_2017}
The exact numbers for $|\mathcal{M}(N)|$ have been determined via exact
enumeration studies for clusters of size $N \leq
14$.\autocite{Hoy_Structuredynamicsmodel_2015,Holmes-Cerfon_EnumeratingRigidSphere_2016}
Studies of this type are difficult to carry out, because they are
computationally expensive.\autocite{Heiles_Globaloptimizationclusters_2013}

Another interaction potential often used in cluster science is the \acf{LJ}
potential, as introduced in section~\ref{sec:LennardJones}. It is most commonly
used in the $(6,12)$ form, but in this section it will be employed with
arbitrary integer exponents $(m,n)$.
%
\begin{equation}
    V_{m,n}^\mathrm{LJ}(r)=\frac{\varepsilon}{n-m}\left[m\left(\frac{r_e}{r}\right)^{n}-n\left(\frac{r_e}{r}\right)^{m}\right] \ \ \ \ \ \ \ \ \ \  (\mathrm{with}\ n > m)
\label{eqn:nmpot}
\end{equation}
%
The two parameters $\varepsilon$ and $r_e$ are the depth of the potential energy
well and the equilibrium distance, respectively. The values for $\varepsilon$
and $r_e$ will be set to unity in the following for the same reasons as for the
\ac{SHS} potential. As already mentioned in chapter~\ref{sec:energylandscapes},
the \ac{SHS} potential emerges from the \ac{LJ} potential as the limit for large
exponents $(m,n)$ (figure~\ref{fig:LJ}).
%
\begin{figure}[htb]\centering
    \includegraphics[width=0.8\columnwidth]{kslj/exampleLJ.pdf}
    \caption{Lennard-Jones potentials for different exponents $(m,n)$ with fixed
    $n=2m$. As the exponents grow larger, the well of attraction becomes
    narrower and its shape approaches the \acs{SHS} potential. The dashed line
    shows the extended Lennard-Jones potential for the xenon dimer
    \autocite{Jerabek_relativisticcoupledclusterinteraction_2017}.}
    \label{fig:LJ}
\end{figure}

At first, it seems surprising that the absolute number of structures for a
certain size $N$ differs substantially between the two potentials. For $N=13$
there are $|\mathcal{M}_\mathrm{SHS}|=97,221$ non-isomorphic \ac{SHS} clusters
\autocite{Hoy_Structuredynamicsmodel_2015,Holmes-Cerfon_EnumeratingRigidSphere_2016},
but only $|\mathcal{M}_\mathrm{LJ}|=1,510$ $(6,12)$-\ac{LJ}
clusters\autocite{Doye_Evolutionpotentialenergy_1999}. However, this is a known
behaviour of energy landscapes of long-range (\ac{LJ}) and short-range
(\ac{SHS}) potentials, with the latter generally supporting many more local
minima compared to the
former.\autocite{braier90,Wales_MicroscopicBasisGlobal_2001} Decreasing the
exponents $(m,n)$ increases the range of the potential, which leads to increased
second-nearest-neighbour interactions. Furthermore, fold
catastrophes\autocite{Wales_MicroscopicBasisGlobal_2001,Wales_Energylandscapes_2003}
lead to the collapse of several stable \ac{SHS} structure into a single \ac{LJ}
minimum, leading to a decrease in the overall size of $|\mathcal{M}(N)|$.

In the following sections, the evolution of \ac{LJ} clusters towards \ac{SHS}
clusters by gradually decreasing the range of the \ac{LJ} potential is explored.
Additionally, the results from optimising with a \ac{LJ} potential starting from
the \ac{SHS} cluster is compared to the traditional approach of global
optimisation.

%The nucleation of atoms and molecules in the gas phase, or liquid, to the solid
%state is still an active research
%field\autocite{Stillinger_Packingstructurestransitions_1984,
%Martin-1996,Wales-1996, Vlieg_atomicscaleunderstandingcrystal_2007, Arkus-2010,
%Woodley-2010, Karthika-2016, Holmes-Cerfon_StickySphereClusters_2017}. Rowland
%noted in 1949 that ``The gap between theory and the experimental approaches to
%nucleation has been too wide'' and ``the subject [nucleation] is still in the
%alchemical stage'' \autocite{Rowland-1949}. More than half a century later,
%despite all the advancements made in cluster physics, ``there is still a large
%gap between experiment and theory'' as Unwin noted \autocite{Unwin-2007}.

%The underlying reason for this rather slow progress is that cluster formation
%is a dynamic process, and fully characterizing the corresponding
%high-dimensional potential energy landscape is typically an NP-hard problem,
%since there are (presumably) exponentially many local minima at any given
%temperature and pressure\autocite{Stillinger_Packingstructurestransitions_1984,
%Oganov-2006, Massen_Powerlawdistributionsareas_2007, wales10, Oganov-2011,
%calvo12, Wales-2015}.  Moreover, phase transitions between different
%morphologies as a function of size $N$ usually occur where $N$ is too large for
%an accurate quantum-theoretical treatment \autocite{Waal89,
%Cleveland_energeticsstructurenickel_1991, vandewaal96a, Doye-1995,
%vandeWaalTd00}.  For example, Krainyukova experimentally studied the growth of
%argon clusters \autocite{Krainyukova-2012}, and found that small, initially
%icosahedral clusters transform into anti-Mackay clusters for $N>2000$, and
%finally into the \acf{fcc} or \acf{hcp} structures at $N>10^5$ atoms, in
%qualitative agreement with theoretical predictions using \ac{LJ} type
%potentials \autocite{Martin-1996,
%Schwerdtfeger_ExtensionLennardJonespotential_2006, Krainyukova-2007}.  The
%notorious \textit{rare gas problem} was solved only very recently by accurate
%relativistic quantum methods, correctly predicting a slight preference of the
%\ac{fcc} over the hcp phase due to phonon dispersion \autocite{Schwerdtfeger-2016}.

%Simple models often have to be used to simulate cluster growth and nucleation
%\autocite{Johnston-1999, Shibuta_Homogeneousnucleationmicrostructure_2015,
%Leitold-2016, Sweatman-2016}.  The simplest model potentials that can be
%applied to theoretical studies of atomic cluster formation are ``HCR-SRA''
%potentials with isotropic hard-core-like repulsive and short-range-attractive
%interactions \autocite{baxter68}.  The simplest HCR-SRA potential is the
%\ac{SHS} potential \autocite{Yuste_Stickyhardspheres_1993}

%where $r_s$ and $\varepsilon$ can be arbitrarily set to 1 (unit sphere and
%reduced units, respectively). Equation~\eqref{eqn:KS} can be used as a perturbative
%basis for finite-ranged HCR-SRA potentials \autocite{cochran06,
%Holmes-Cerfon_geometricalapproachcomputing_2013}.  Since sticky hard spheres
%are impenetrable and their energy $E = -N_c\varepsilon$ ~is a function only of
%the number of interparticle contacts $N_c$, \ac{SHS} cluster structure and
%energetics can be uniquely mapped to their adjacency matrices $\bar{A}$, where
%$N_c=\sum_{i<j}^N A_{ij}$.  This mapping allows them to be exactly
%characterized via complete enumeration
%\autocite{Arkus_Minimalenergyclusters_2009, Arkus_DerivingFiniteSphere_2011,
%Hoy_Structurefinitesphere_2012}; recent studies have identified all
%mechanically stable \ac{SHS} clusters for $N \leq 14$, and putatively complete sets
%for $N \leq 19$ \autocite{Arkus_Minimalenergyclusters_2009,
%Arkus_DerivingFiniteSphere_2011, Hoy_Structurefinitesphere_2012,
%Hoy_Structuredynamicsmodel_2015, Holmes-Cerfon_StickySphereClusters_2017,
%Kallus_Freeenergysingular_2017}.  Note, however, that different \ac{SHS} structures
%can have the same adjacency matrix for $N \geq 14$
%\autocite{Holmes-Cerfon_EnumeratingRigidSphere_2016}, and the mapping is
%therefore only surjective.

%From the Gregory-Newton kissing-number argument proved in 1953 by Sch\"utte and
%van der Waerden \autocite{Schutte_ProblemdreizehnKugeln_1952}, no sphere can be
%surrounded by more than 12 spheres of equal radius \autocite%{Conway_SpherePackingsLattices_1999}.  For
%small clusters, graph-theoretic arguments dictate $\mathrm{max}(N_c)\le
%N(N-1)/2$.  Thus a loose bound on the maximum contact number $N_c(N)$ is
%%
%%
%with $f(N)=6N$.  This upper bound has been tightened several times, most
%recently by Bezdek and Reid \autocite{Bezdek-2013} to
%

%
%In
%references~\cite{Hoy_Structuredynamicsmodel_2015,Holmes-Cerfon_EnumeratingRigidSphere_2016}
%it was shown that $N_c^\mathrm{max}(N) =
%\{6,9,12,15,18,21,25,29,33,36,40,44,48,52,56,60\}$ for $4 \leq N \leq 19$.
%While determining $N_c^\mathrm{max}(N)$ for arbitrary $N$ is equivalent to the
%still-unsolved Erd\"os unit distance problem \autocite{Erdos-1946}, it is clear
%that $N_c^\mathrm{max}(N) = 3N - 6 + m(N)$, where $m(N)$ grows slowly from zero
%to around $f(N) - (3N-6)$ with increasing $N$.

%While the maximum contact number increases (sub)linearly with $N$, the number
%of non-isomorphic cluster structures $|\mathcal{M}(N)|$ and transition states
%is assumed to increase exponentially
%\autocite{Stillinger_Exponentialmultiplicityinherent_1999,Oganov-2006,%Forman_ModelingAggregationProcesses_2017}
%(here we denote $\mathcal{M}(N)$ as the set of all non-isomorphic cluster
%structures of size $N$, and $|\mathcal{M}(N)|$ as the number of structures in
%$\mathcal{M}(N)$).  Stillinger showed that under certain conditions
%$\lim_{N\to\infty} |\mathcal{M}(N)| \propto \exp(\alpha N)$
%\autocite{Stillinger_Exponentialmultiplicityinherent_1999}.  For \ac{SHS}
%clusters, the complete set $\mathcal{M}_\mathrm{SHS}(N,N_c)$ has been exactly
%determined for $N \leq 14$ and $3N - 6 \leq N_c \leq N_c^\mathrm{max}(N)$ via
%exact enumeration studies employing geometric rejection rules
%\autocite{Hoy_Structuredynamicsmodel_2015,Holmes-Cerfon_EnumeratingRigidSphere_2016}.
%Unfortunately, such precise calculations are very difficult for finite-ranged
%potentials since exhaustive searches for energy minima are computationally
%intensive \autocite{Heiles_Globaloptimizationclusters_2013}.  Only a few such
%studies have been performed, e.g.\ recent studies of $N \leq 19$ clusters
%interacting via short-range Morse potentials
%\autocite{wales10,calvo12,C7CP03346J}.


%
%It remains unclear how the HCR-SRA models commonly
%used in cluster physics relate to more physically relevant, softer interaction potentials %such
%as the $(m,n)$-Lennard-Jones (LJ) form:
%
%Here $\varepsilon>0$ is the dissociation energy and $r_e$ the equilibrium
%two-body interparticle distance. To simplify the presentation,
%we (without loss of generality) adopt reduced units ($\varepsilon=1$, $r_e=1$) below.
%For $m,n\rightarrow \infty$,
%$V_{m,n}^\mathrm{LJ}(r) \rightarrow V_\mathrm{SHS}(r)$ (figure~\ref{fig:LJ}); the
%energy landscapes of the two potentials converge in this limit.  However, real
%systems are not in this limit.  For example, for $N = 13$, there are $|\mathcal{M}_\mathrm%{SHS}|=97,221$
%stable \ac{SHS} clusters \autocite{Hoy_Structuredynamicsmodel_2015,%Holmes-Cerfon_EnumeratingRigidSphere_2016},
%but only $|\mathcal{M}_\mathrm{LJ}|=1,510$ stable $(m,n) = (6,12)$ LJ clusters \autocite%{Doye_Evolutionpotentialenergy_1999}.  
%This difference is understood qualitatively -- energy landscapes are well known to support %more
%local minima
%as the range of the interaction potential decreases \autocite{braier90,%Wales_MicroscopicBasisGlobal_2001}.
%There are several effects that will cause the set of
%stable LJ clusters to increasingly deviate from the set of stable \ac{SHS} clusters
%as interactions become longer ranged.  As $n$ and $m$ decrease,
%second-nearest-neighbor attractions become increasingly important,
%producing stable structures with $r_{ij} \leq 1$.  Fold catastrophes
%\autocite{Wales_MicroscopicBasisGlobal_2001,Wales_Energylandscapes_2003} progressively %eliminate stable \ac{SHS} clusters, and several stable \ac{SHS}
%structures may collapse into a single stable LJ cluster.  However,
%detailed quantitative understanding of such effects remains rather limited.

%In this paper, we quantitatively examine how stable $N \leq 14$ LJ cluster
%structures evolve away from the \ac{SHS} limit as ($m, n$) decrease.  We focus on
%both the topography of the energy landscape (decreasing
%$|\mathcal{M}_{\mathrm LJ}(N)|$) and the evolving topologies of the stable cluster sets.
%We examine these changes in further detail for specific $N = 13-14$ clusters discussed by
%Gregory and Newton in the 1600s in the context of the kissing number problem
%\autocite{Schutte_ProblemdreizehnKugeln_1952}, and also for a more realistic two-body %potential that has
%been shown to accurately model rare-gas clusters \autocite{Schwerdtfeger_ExtensionLennardJonespotential_2006}.


\section{Computational Details}

The \textit{pele} program~\autocite{_pelePythonenergy_2017} was used to generate
putatively complete sets of local minima for $(m,n)$-Lennard-Jones potentials
$V_{mn}^{\mathrm LJ}(r)$ as defined in equation~\eqref{eqn:nmpot}. This program
applies a basin-hopping algorithm that divides the potential energy surface into
basins of attraction, effectively mapping each point in configuration space to a
local minimum structure
\autocite{Li_MonteCarlominimizationapproach_1987,waless99,Wales_GlobalOptimizationBasinHopping_1997}.
The results confirmed the number of local minima reported in previous work
\autocite{Doye_Saddlepointsdynamics_2002}. Finite computer time limited the
search to clusters of size $N \leq 13$.

Starting from the sticky hard sphere packings up to $N=14$, with Cartesian
coordinates given by the exact enumeration algorithm
\autocite{Hoy_Structurefinitesphere_2012} including rigid hypostatic clusters
($N_c<3N-6$) \autocite{Holmes-Cerfon_EnumeratingRigidSphere_2016}, geometry
optimisations with $(m,n)$-Lennard-Jones potentials using the multidimensional
function minimiser from the \Cpp library \textit{dlib}
\autocite{King_DlibmlMachineLearning_2009} were carried out with the previously
described program package \textsc{Spheres}
(Chapter~\ref{sec:theprogramspheres}). The optimisation scheme was either the
\acf{BFGS} or the conjugate gradient algorithm. The optimisations were terminated
when the change in energy (in reduced units) over the course of one
optimisation cycle was smaller than $10^{-15}$. 

Subsequently, the eigenvalues of the Hessian were checked for all stationary
points. If negative eigenvalues were found, the affected structures were
re-optimised following displacements in both directions along the corresponding
eigenvectors to locate true local minima. This procedure assures that the floppy
\ac{SHS} packings are successfully mapped into \ac{LJ} minima.

As the optimisations often result in many duplicates, especially for small
values of $n$ and $m$ where $|\mathcal{M}_{(m,n)\mathrm{-LJ}}| \ll
|\mathcal{M}_\mathrm{SHS}|$, the final structures were further analysed and
sorted. Non-isomorphic \ac{SHS} clusters can be distinguished (apart from
permutation of the particles) by their different adjacency matrices for $N \leq
13$ \autocite{Holmes-Cerfon_EnumeratingRigidSphere_2016}. This is not the case
for soft potentials like the \ac{LJ} potential since drawing edges (bonds)
between the vertices (atoms) becomes a matter of defining the distance cut-off
criterion for a bond to be drawn. Therefore, the clusters were compared based on
the \ac{EDM} (the matrix of inter-particle distances $\{r_{ij}\}$) as described
previously: two clusters are isomorphic (structurally identical) if they have
the same ordered set of inter-particle distances $\{r_{ij}\}$. While enantiomers
cannot be separated using this methodology, permutation-inversion isomers are
usually lumped together, since the number of distinct minima is analytically
related to the order of the corresponding point group
\autocite{Wales_Energylandscapes_2003}. To verify the number of distinct
structures a second ordering scheme using the energy and moment of inertia
tensor eigenvalues was introduced. 

Two sets of structures are obtained from the optimisation procedure: the first
set contains all possible \ac{LJ} minima $\mathcal{M}_\mathrm{LJ}$ from the
basin-hopping algorithm, while the second set $\mathcal{M}_\mathrm{SHS\to LJ}$
contains the \ac{LJ} minima obtained using only the $\mathcal{M}_\mathrm{SHS}$
sticky-hard-sphere cluster structures as starting points for the geometry
optimisation. To compare and identify corresponding structures between the two
sets, the $N(N-1)/2$ inter-particle distances $\{r_{ij}\}$ were again used as an
identifying fingerprint.

Two-body \acf{eLJ} potentials that accurately model two-body interactions in
rare-gas clusters can be written as expansions of inverse-power-law terms
\autocite{Schwerdtfeger_ExtensionLennardJonespotential_2006}:
%
\begin{equation} \label{eq:ELJ}
V_{\mathrm ELJ}(r)=\sum_{n} c_nr^{-n},
\end{equation}
%
where in reduced units the condition $\sum_{n} c_n=-1$ holds. For comparison
to the simple (6,12)-\ac{LJ} potential, the \ac{eLJ} potential derived from
relativistic coupled-cluster theory applied to the xenon dimer was used with the
following coefficients (in reduced units):
$c_6=-1.0760222355$; $c_8=-1.4078314494$; $c_9=-185.6149933139$;
$c_{10}=+1951.8264493941$; $c_{11}=-8734.2286559729$;
$c_{12}=+22273.3203327203$; $c_{13}=-35826.8689874832$;
$c_{14}=+37676.9744744424$; $c_{15}=-25859.2842295062$;
$c_{16}=+11157.4331408911$; $c_{17}=-2745.9740079192$; $c_{18}=+293.9003309498$
\autocite{Jerabek_relativisticcoupledclusterinteraction_2017}. The \ac{eLJ}
potential for xenon is shown in figure~\ref{fig:LJ} (dashed line).


%Exploring the limits of Lennard-Jones
\section{Exploring the Limits of Lennard-Jones}

\begin{table}[htbp]\centering
    \resizebox{\linewidth}{!}{%
        \begin{threeparttable}[para] \caption{Number of distinct local minima
        $|\mathcal{M}_\mathrm{SHS}|$ for cluster size $N$ (from
        references~\cite{Holmes-Cerfon_EnumeratingRigidSphere_2016,Hoy_Structurefinitesphere_2012,Hoy_Structuredynamicsmodel_2015})
        and contact number $N_c$ from the exact enumeration, compared to the
        number of different structures obtained from a geometry optimisation
        starting from the set $\mathcal{M}_{\mathrm{SHS}}(N,N_c)$ for a
        (6,12)-\acs{LJ} potential. The overall number of unique minima
        $|\mathcal{M}_\mathrm{SHS\to LJ}|  = \sum_{N_c}
        |\mathcal{M}_\mathrm{SHS\to LJ} (N_c)| - \text{(\# of duplicate
        structures)}$ is shown in the following column.  This result can be
        compared to the number of unique minima found using the basin-hopping
        method ($|\mathcal{M}_\mathrm{LJ}|$). The difference
        $\Delta\mathcal{M}=|\mathcal{M}_\mathrm{LJ}| -
        |\mathcal{M}_\mathrm{SHS\to LJ}|$ is also listed.}
        \label{tab:comp}
        \begin{tabular}{clccccc}\toprule
            $N$ & $N_c$ & $|\mathcal{M}_\mathrm{SHS} (N_c)|$ & $|\mathcal{M}_{\mathrm{SHS\to LJ}}(N_c)|$ & $|\mathcal{M}_\mathrm{SHS\to LJ}|$ & $|\mathcal{M}_\mathrm{LJ}|$     &  $\Delta\mathcal{M}$ \\ \midrule
    8  & 18  & 13     				& 8    & 8                     & 8                     & 0                   \\  %\midrule
    9  & 21  & 52     				& 20   & 20                    & 21                    & 1                   \\  \midrule
    10 & 23  & 1      				& 1    &                       &                       &                     \\
       & 24  & 259    				& 60   & 62                    & 64                    & 2                   \\
       & 25  & 3      				& 3    &                       &                       &                     \\  \midrule
    11 & 25  & 2      				& 2    & \multirow{5}{*}{165}  & \multirow{5}{*}{170}  & \multirow{5}{*}{5}  \\
       & 26  & 18     				& 6    &                       &                       &                     \\
            & 27  & 1620\tnote{a}	& 158  &                       &                       &                     \\
       & 28  & 20     				& 12   &                       &                       &                     \\
       & 29  & 1      				& 1    &                       &                       &                     \\  \midrule
    12 & 28  & 11     				& 6    & \multirow{6}{*}{504}  & \multirow{6}{*}{515}  & \multirow{6}{*}{11} \\
       & 29  & 148    				& 24   &                       &                       &                     \\
       & 30  & 11638  				& 483  &                       &                       &                     \\
       & 31  & 174    				& 69   &                       &                       &                     \\
       & 32  & 8      				& 6    &                       &                       &                     \\
       & 33  & 1      				& 1    &                       &                       &                     \\  \midrule
    13 & 31  & 87     				& 23   & \multirow{6}{*}{1476} & \multirow{6}{*}{1510} & \multirow{6}{*}{34} \\
       & 32  & 1221   				& 100  &                       &                       &                     \\
            & 33  & 95810\tnote{a}& 1418 &                       &                       &                     \\
            & 34  & 1318\tnote{a} & 293  &                       &                       &                     \\
       & 35  & 96     				& 49   &                       &                       &                     \\
       & 36  & 8      				& 6    &                       &                       &                     \\  \midrule
            14 & 33  & 1      				& 1    & \multirow{8}{*}{4093} & \multirow{8}{*}{(4187)\tnote{b}}    & \multirow{8}{*}{(94)\tnote{b}}  \\
       & 34  & 707    				& 101  &                       &                       &                     \\
       & 35  & 10537  				& 410  &                       &                       &                     \\
       & 36  & 872992 				& 3939 &                       &                       &                     \\
       & 37  & 10280  				& 1002 &                       &                       &                     \\
       & 38  & 878    				& 237  &                       &                       &                     \\
       & 39  & 79     				& 42   &                       &                       &                     \\
       & 40  & 4      				& 3    &                       &                       &                     \\\bottomrule
        \end{tabular}
            \begin{tablenotes}
            \item[a]{Largest value for $|\mathcal{M}_\mathrm{SHS}|$ taken from references~\cite{Holmes-Cerfon_EnumeratingRigidSphere_2016,Hoy_Structurefinitesphere_2012,Hoy_Structuredynamicsmodel_2015}.} \item[b]{Estimated.}
            \end{tablenotes}
        \end{threeparttable}}
    \end{table}%
%
To study the convergence behaviour of the number of distinct (non-isomorphic)
\ac{LJ} minima in the \ac{SHS} limit, geometry optimisations were carried out,
starting from all non-isomorphic \ac{SHS} structures. It will be shown later
that the number of unique minima obtained in this procedure
$|\mathcal{M}_\mathrm{SHS\to LJ}|$ only misses out on a small portion of
minima obtained from the more exhaustive basin-hopping approach,
i.e.~$|\mathcal{M}_\mathrm{SHS\to LJ}|\approx|\mathcal{M}_\mathrm{LJ}|$.
The results for a constant chosen ratio of \ac{LJ} exponents $n/m=2$ are shown
in figure~\ref{fig:expinfty} (top).
%
\begin{figure}[htbp]
    \centering
    \subfloat{\includegraphics[width=0.8\columnwidth]{kslj/expinfty.pdf}}\\
    \subfloat{\includegraphics[width=0.8\columnwidth]{kslj/repulsive.pdf}}
    \caption{Convergence of the number of distinct \acs{LJ} local minima
    $|\mathcal{M}_\mathrm{SHS\to LJ}|$ obtained through geometry optimisations
    starting from the non-isomorphic \acs{SHS} structures with increasing
    \acs{LJ} exponent $n$. Permutation-inversion isomers and enantiomers are not
    distinguished. The dashed line gives the exact \acs{SHS} limit
    $|\mathcal{M}_\mathrm{SHS}|$. Top panel: $m=n/2$. Bottom panel: fixed
    $m=6$.}
    \label{fig:expinfty}
\end{figure}
%
$|\mathcal{M}_\mathrm{SHS\to LJ}|$ smoothly converges towards the \ac{SHS} limit
(dashed line, values in table~\ref{tab:comp}) from below, thus demonstrating
that for \ac{LJ} systems the number of distinct minima does not grow faster than
exponentially. The (48,96)-\ac{LJ} potential has $\Delta\mathcal{M} \equiv
|\mathcal{M}_\mathrm{LJ}| - |\mathcal{M}_\mathrm{SHS\to LJ}| =
\{1,1,7,91,1019,14890,209938\}$ fewer stable minima than the \ac{SHS} potential.
The fractions of missing minima $\Delta\mathcal{M}/|\mathcal{M}_\mathrm{SHS}|$
for this potential grow with increasing $N$ and are, respectively,
$\{7.69,1.92,2.67,5.46,8.62,15.32,23.44\}\%$.  Note that for $N \geq 10$ most of
these missing minima correspond to high energy ($N_c < N_c^\mathrm{max}$)
structures.

If the exponent $n$ for the repulsive part of the \ac{LJ} potential is increased
with $m$ kept constant, the \ac{LJ} potential becomes equivalent to the \ac{SHS}
potential in the repulsive range but remains attractive at long range. This
limit is also called the Sutherland potential. Figure~\ref{fig:expinfty}
(bottom) shows the convergence of the number of unique structures with respect
to $n$ at set $m=6$ towards the \ac{SHS} limit. Here, the number of distinct
minima converges towards a number that is much smaller than the total number of
\ac{SHS} packings demonstrating that (as expected) the attractive part of the
potential contributes significantly to the decrease of the number of local
minima compared to the rigid \ac{SHS} model.

To see if the asymptotic increase in the number of distinct minima
$|\mathcal{M}(N)| \sim e^{\alpha N}$ is indeed exponential, an expression for
the asymptotic exponential rise rate parameter developed by Stillinger was
used:\autocite{Stillinger_Exponentialmultiplicityinherent_1999}
%
\begin{equation} \label{eq:Stil}
\alpha = \lim_{N\rightarrow \infty} \left( N^{-1} \mathrm{ln} |\mathcal{M}(N)| \right).
\end{equation}
%
Figure~\ref{fig:asympt} shows the number of distinct minima for \ac{SHS}
clusters obtained from the data shown in table~\ref{tab:comp}. 
%
\begin{figure}[htb]
    \centering
    \includegraphics[width=0.8\columnwidth]{kslj/growth.pdf}
    \caption{Growth behaviour of $|\mathcal{M}(N)|$ of \acs{SHS} and (6,12)-\acs{LJ}
    clusters and corresponding asymptotic exponential rise rate parameter
    $\alpha$ for $N \geq 12$ as defined in equation~\eqref{eq:Stil}.  The
    intercepts $\ln|\mathcal{M}(N=0)|$ are $-17.19$ and $-6.94$ for the \acs{SHS}
    and (6,12)-\acs{LJ} cases, respectively.}
    \label{fig:asympt}
\end{figure}
%
The $N \geq 12$ \ac{SHS} data gives $\alpha_\mathrm{SHS}\approx 2.21$.
Figure~\ref{fig:asympt} also shows the (6,12)-\ac{LJ} results obtained using
basin-hopping; these yield $\alpha_\mathrm{LJ}\approx 1.10$, which is close to
the $\alpha=0.8$ value estimated by Wallace \autocite{Wallace-1997} or to the
recently given value of $\alpha=1.04$ by Forman and Cameron
\autocite{Forman_ModelingAggregationProcesses_2017}. Note that the rapid
increase of $|\mathcal{M}_\mathrm{SHS}|/|\mathcal{M}_\mathrm{LJ}|$ with $N$ is
explained by the much larger values of $\alpha$ for the \ac{SHS} compared to the
\ac{LJ} clusters.

Using the results for $N \geq 13$ depicted in figure~\ref{fig:expinfty}, the
dependence of $\alpha$ on the \ac{LJ} range parameter $n$ can be calculated. As
shown in figure~\ref{fig:repulsive13-14}, a general function of the form
%
\begin{align}
\label{expgrowth}
    \alpha(n)=\alpha_\text{max}+\frac{a}{(n-n_0)^{p}}
\end{align}
%
fits the results nicely, allowing the prediction of growth behaviour for
different \ac{LJ} potentials. 
%
\begin{figure}[htb]
    \centering
    \includegraphics[width=0.8\columnwidth]{kslj/repulsive13-14.pdf}
    \caption{Convergence behaviour of the asymptotic exponential rise rate
    parameter $\alpha$ (equation~\eqref{eq:Stil}) towards the \acs{SHS} limit with
    respect to the \acs{LJ} exponent $n$. The inlet shows the ratio of the two
    quantities $\alpha(|\mathcal{M}_{\text{SHS}\to (n/2,n)-\text{LJ}}(N)|)
    / \alpha(|\mathcal{M}_{\text{SHS}\to (6,n)-\text{LJ}}(N)|)$.}
    \label{fig:repulsive13-14}
\end{figure}
%
For $|\mathcal{M}_{(n/2,n)-\text{LJ}}|$, $\alpha_\text{max}$ is equivalent to
$\alpha_\text{SHS}=2.207$. The other adjusted parameters are $a=-66.588$,
$n_0=-3.386$ and $p=1.473$ (figure~\ref{fig:repulsive13-14}). We also show the
ratio $\alpha(|\mathcal{M}_{\text{SHS}\to (n/2,n)-\text{LJ}}|) /
\alpha(|\mathcal{M}_{\text{SHS}\to (6,n)-\text{LJ}}|)$ between the two different
\ac{LJ} asymptotic exponential rise rate parameters, which shows that larger
cluster sizes need to be studied to correctly describe the asymptotic limit. 

The distribution of minima as a function of (free) energy was suggested to be
Gaussian \autocite{Sciortino-1999}. Figure~\ref{fig:N13-steps} shows the energy
distribution of minima for different $(n/2,n)$-\ac{LJ} potentials derived from
\ac{SHS} initial structures.
%
\begin{figure}[htb]
    \centering
    \includegraphics[width=0.8\columnwidth]{kslj/N13-steps.pdf}
    \caption{Histogram of the energies (bin size $\Delta E=0.1$) of
    minima $\mathcal{M}_{\text{SHS}\to (n/2,n)-\text{LJ}}(N)$ for $N=13$ and
    different exponents $n$ up to the \acs{SHS} limit. For better visibility,
    the height of the bars are set to $\Delta|\mathcal{M}|/|\mathcal{M}|$ in the
    interval $\Delta(E/\epsilon)$. The inlet shows the same data in logarithmic
    scale.}
    \label{fig:N13-steps}
\end{figure}
%
A Gaussian type distribution was not observed; this result does not change if
the free energy at finite temperatures is used instead. The results indicate a
``phase transition'' in the potential energy landscape away from low energy to
high energy minima as $n$ increases. The transition occurs at fairly small $n$.
Results for the $(9,18)$-\ac{LJ} potential indicate two \ac{SHS}-like maxima
that are not present for the $(6,12)$-\ac{LJ} potential; these are associated
with the $N_c = 34$ and $N_c = 35$ \ac{SHS} clusters, respectively. It is also
clear that (as expected) the distributions narrow with increasing $n$.

It is well known that the global minimum for rare gas clusters with 13 atoms is
the ideal Mackay icosahedron
\autocite{Hoare_Physicalclustermechanics_1975,Hoare_Statisticalmechanicsmorphology_1976,Hoare_StructureDynamicsSimple_2007}.
Simple geometric considerations imply that such a symmetric cluster is not
possible for sticky hard spheres; all vertices of a regular icosahedron with
unit edge length lie on a circumscribing sphere with radius $r_c\approx 0.951$,
making it impossible to insert a sphere of the same radius into the centre of
the polyhedron. Therefore, there must be well-defined \ac{LJ} exponents $(m,n)$
at which the icosahedral $N = 13$ \ac{LJ} cluster breaks symmetry to form a
rigid cluster. For the $n = 2m$ case considered above, this symmetry-breaking
occurs at $m \simeq 15$.

We also explored a more realistic \ac{eLJ} potential (equation~\eqref{eq:ELJ};
figure~\ref{fig:LJ}) for one of the rare gas dimers (xenon) in comparison with
other \ac{LJ} potentials. It can be seen that the repulsive part agrees nicely
with the conventional (6,12)-\ac{LJ} potential, while for $r > 1$ the extended
\ac{LJ} potential is slightly less attractive. This change should lead to an
increase in the number of local minima compared to the conventional
(6,12)-\ac{LJ} potential. This prediction could be confirmed, i.e.
$|\mathcal{M}_\mathrm{SHS\to ELJ}|=\{8,21,74,205,685,2179,6863\}$ for
$N=\{8,9,10,11,12,13,14\}$. For $N=13$ the number of distinct minima is 44\%
larger than for the simple (6,12)-\ac{LJ} potential, which shows that
$|\mathcal{M}(N)|$ is rather sensitive to the potential chosen. Hence, to
correctly describe the topology of real systems, one has to take care of the
correct form of the 2-body contribution (as well as higher $n$-body
contributions) \autocite{Schwerdtfeger-2016}.


%\subsection{(6,12)-Lennard-Jones cluster from basin-hopping}
\section{(6,12)-Lennard-Jones Clusters from Basin-Hopping} 
\label{sec:612LennardJonesClustersfromBasinHopping}

Table~\ref{tab:comp} shows the number of distinct minima found by the cluster
geometry optimisation procedure employed in this work using the (6,12)-\ac{LJ}
potential compared to results from exact enumeration for \ac{SHS}s and from
basin-hopping for the (6,12)-\ac{LJ} potential. As the \ac{SHS} clusters for a
specific $N$ value can be grouped by their contact number $N_c$, the geometry
optimisations were carried out separately for each group of
$\mathcal{M}_\mathrm{SHS}(N_c)$.
\citeauthor{Hoy_Structurefinitesphere_2012}
\autocite{Hoy_Structurefinitesphere_2012,Hoy_Structuredynamicsmodel_2015} and
\citeauthor{Holmes-Cerfon_EnumeratingRigidSphere_2016}\autocite{Holmes-Cerfon_EnumeratingRigidSphere_2016}
reported slightly different results for $N=11$ and $N=13$; however, upon
geometry optimisation, their datasets yield the same final clusters
$|\mathcal{M}_{\mathrm{SHS\to LJ}}(N_c)|$. As identical \ac{LJ} clusters appear
in multiple groups with different contact numbers, the duplicates were removed
to create the set $\mathcal{M}_\mathrm{SHS\to LJ}$ of distinct minima, which can
be directly compared to the set of \ac{LJ} minima $\mathcal{M}_\mathrm{LJ}$
obtained from the basin-hopping method. It should be noted that including the
hypostatic clusters and the different $|\mathcal{M}_\mathrm{SHS}|$ for $N=11$
and $N=13$ from
\citeauthor{Holmes-Cerfon_EnumeratingRigidSphere_2016}\autocite{Holmes-Cerfon_EnumeratingRigidSphere_2016}
did not change our results, implying that hypostatic clusters are not an
important feature for the \ac{LJ} energy landscape. 

Interestingly, the gradient-based minimisation procedure employed here does not
in general lead to a complete set of \ac{LJ} minima; the mapping from \ac{SHS}
minima to \ac{LJ} minima is non-injective and non-surjective. Clearly, some
structural motifs found in \ac{LJ} clusters are not found in \ac{SHS} clusters
and vice versa, and the topology of the hypersurface changes in a non-trivial
fashion from \ac{SHS} to \ac{LJ}. However, it is surprising that the fraction of
structures that are missed by this optimisation procedure is so small (see
table~\ref{tab:seeds}). 
%
\begin{table}[htb]\centering
    \begin{threeparttable}
    \caption{Number of missing structures after optimisation belonging to the
    same "seed" (figure~\ref{fig:seeds}). $N=8$ is excluded because all \acs{LJ}
    minima were found starting from the \acs{SHS} model.}
    \label{tab:seeds}
    \begin{tabular}{llllll}\toprule
        seed      & $N=9$   & $N=10$  & $N=11$  & $N=12$  & $N=13$  \\ \midrule
        a         & 1    & 1    & -    & 3    & 8    \\
        b         & -    & 1    & 3    & 4    & 12\tnote{a}   \\
        c         & -    & -    & 1    & 1\tnote{a}    & -    \\
        d         & -    & -    & 1    & 1    & 5    \\
        e         & -    & -    & -    & 1    & 6    \\
        f         & -    & -    & -    & 1    & 1    \\
        remaining & -    & -    & -    & -    & 2    \\ 
        total     & 1    & 2    & 5    & 11   & 34   \\
        \%        & 4.76 & 3.13 & 2.94 & 2.14 & 2.25 \\ \bottomrule
    \end{tabular}
        \begin{tablenotes}
        \item[a]{Some structures do not resemble a perfect capped
        cluster, but undergo a slight rearrangement. Specifically, two structures belonging to seed (b) and one structure belonging to seed (c) were found to deviate slightly from the perfect arrangement, but minor rearrangements of these structures lead to the desired geometry and they can be reasonably associated with these seeds.}
        \end{tablenotes}
    \end{threeparttable}
\end{table}%
%
To gain further insights, the energetics and structures of the unmatched
clusters were investigated in more detail.

Figure~\ref{fig:bondlength-variance} shows an analysis of the difference between
the longest to the shortest bond lengths $d_\Delta=d_{\mathrm max}-d_{\mathrm
min}$ obtained for the largest clusters in $\mathcal{M}_{\mathrm{LJ}}$ with
$N=\{11,12,13\}$.\footnote{We define spheres that have a equilibrium distance
between $0.9-1.1$ to be bound.}
%
\begin{figure}[htb]
    \centering
    \subfloat[$N=11$]{\includegraphics[width=.5\textwidth]{kslj/lj11var.pdf}}
    \subfloat[$N=12$]{\includegraphics[width=.5\textwidth]{kslj/lj12var.pdf}}\\
    \subfloat[$N=13$]{\includegraphics[width=.5\textwidth]{kslj/lj13var.pdf}}
    \caption{Histograms of the difference between the longest and shortest bond
    distances $d_\Delta=d_\text{max}-d_\text{min}$ for the complete set of
    distinct \acs{LJ} minima $\mathcal{M}_\text{LJ}(N)$ for $N=\{11,12,13\}$.
    Orange bars give the number of distinct structures not contained in
    $\mathcal{M}_\mathrm{LJ}$ as obtained from the basin-hopping algorithm.}
    \label{fig:bondlength-variance}
\end{figure}%
%
The histograms show that the clusters most commonly have a $d_\Delta$ of about
$0.03$.  In contrast, as shown by the orange bars, the unmatched structures have
significantly larger $d_\Delta$ values of at least $0.05$, with most of them
having $d_{\Delta} \simeq 0.06$. This is a first indication of why these
structures are not found by starting from \ac{SHS} packings. The latter only
form bonds of length one, and a large variation in bond length could imply that
a \ac{SHS} packing similar to the \ac{LJ} structure does not exist as the
\ac{SHS} boundary conditions are not satisfied. The data in
table~\ref{tab:energies} shows that the unmatched (UM) structures for a specific
$N$ value have much higher energies compared to the one of the global minimum
(which is set to zero, i.e. $E_0=0$). They are always positioned in the upper
half of the energy spectrum, making them energetically unfavourable. However, no
correlation between $d_\Delta$ and the energetic position of the \ac{LJ}
clusters was found.

\begin{table}[htb]\centering
    \caption{Range $[E_0=0,E_\text{max}]$ of the energy spectrum of all \acs{LJ}
    minima, position of the second lowest minimum structure $E_1$ and position
    of the first unmatched (UM) structure $E_0^\text{UM}$ relative to the
    respective global minimum (in reduced units).}
    \label{tab:energies}
        \begin{tabular}{lccc}\toprule
        $N$ & $E_\text{max}$ & $E_1$ & $E_0^\text{UM}$ \\\midrule
        8   & 1.04   & 0.06    & -           \\
        9   & 2.08   & 0.84    & 1.19        \\
        10  & 3.13   & 0.87    & 2.22        \\
        11  & 4.22   & 0.85    & 2.27        \\
        12  & 6.16   & 1.62    & 3.38        \\
        13  & 9.26   & 2.85    & 6.14        \\\bottomrule
        \end{tabular}
\end{table}%

Last, the geometries of the missing structures were investigated in more detail.
As it turns out, almost all of the missing stable \ac{LJ} clusters can be
created from a smaller set of missing clusters by capping some of their
triangular faces. Therefore, these groups of clusters can be referred to as
``seeds''.\autocite{Arkus_DerivingFiniteSphere_2011} The corresponding starting
structures of each seed are shown in figure~\ref{fig:seeds}.
%
\begin{figure}[htb]
    \centering
    \subfloat[]{\includegraphics[width=0.26\columnwidth]{kslj/A.png}}
    \subfloat[]{\includegraphics[width=0.26\columnwidth]{kslj/B.png}}
    \subfloat[]{\includegraphics[width=0.26\columnwidth]{kslj/C.png}}\\
    \subfloat[]{\includegraphics[width=0.26\columnwidth]{kslj/D.png}}
    \subfloat[]{\includegraphics[width=0.26\columnwidth]{kslj/E.png}}
    \subfloat[]{\includegraphics[width=0.26\columnwidth]{kslj/F.png}}
    \caption{Graphical representations of the structures that are starting new
    seeds, but are not contained in $\mathcal{M}_\mathrm{SHS\to LJ}$. See
    table~\ref{tab:seeds} and text for more details.}
    \label{fig:seeds}
\end{figure}%
%
None of these structures are stable \ac{SHS} packings. For example, structure
(d) can be described as three octahedra connected via triangular faces sharing
one edge. Geometric
considerations\autocite{Arkus_DerivingFiniteSphere_2011,Hoy_Structurefinitesphere_2012}
immediately show that this structure cannot be a stable \ac{SHS} packing; the
dihedral angle in an octahedron is approximately $109.5^\circ$, which means
three octahedra only fill $328.5^\circ$ of a full circle, leaving a gap between
two faces. 

Table~\ref{tab:seeds} shows the number of missing minima belonging to each seed.
Over 60~\% of the unmatched structures belong to seeds (a) and (b). From a graph
theoretical point of
view,\autocite{Arkus_Minimalenergyclusters_2009,Arkus_DerivingFiniteSphere_2011}
grouping structures into seeds means that all structures belonging to the same
seed contain the graph of the starting structures as a sub-graph in their
respective connectivity matrix. This approach simplifies the analysis to a great
extent, as the feature that prevents the structures from being found by geometry
optimisation is the same for each of the structures arising from a specific
seed. The smallest unmatched structures that cannot be associated with any of
seeds (a)--(f) have $N=13$; these could be the starting structures for two new
seeds.

Finally, it should be noted that the starting \ac{SHS} minima in the
optimisation procedure are not stationary points on the \ac{LJ} hypersurface,
and the optimisations therefore lead to most but not all local and available
\ac{LJ} minima. This observation explains why some high-energy structures were
not found by the optimisation procedure. For a smooth change in the topology of
the potential energy surface from \ac{SHS} to \ac{LJ} type clusters one has to
continuously vary the exponents $(m,n)$ in real space, which is computationally
too demanding.

\section{Conclusion}

The sets of $(m,n)$-\ac{LJ}-potential minima obtained using
complete sets of non-isomorphic \ac{SHS} packings with $8 \leq N \leq 14$
\autocite{Arkus_Minimalenergyclusters_2009,Arkus_DerivingFiniteSphere_2011,Hoy_Structurefinitesphere_2012,Hoy_Structuredynamicsmodel_2015,Holmes-Cerfon_EnumeratingRigidSphere_2016}
as initial states for energy minimisation have been characterised. The number of
distinct minima (i.e.~excluding permutation-inversion isomers) is far smaller
than the number of \ac{SHS} packings for the standard Lennard-Jones exponents
$(m,n) = (6,12)$, but approaches the \ac{SHS} limit from below as ($m,n$)
increase. How the number of distinct minima $\mathcal{M}(N)$ increases with
cluster size $N$ has been investigated by determining Stillinger's rise rate
parameter $\alpha$ (equation~\eqref{eq:Stil})
\autocite{Stillinger_Exponentialmultiplicityinherent_1999}. The increase of
$\alpha$ from $\approx 1.1$ for (6,12)-\ac{LJ} clusters to $\approx 2.2$ for
\ac{SHS} clusters is described by a simple functional form
(equation~\eqref{expgrowth}). All these results can be understood in terms of a
smooth progression of the $(m,n)$-\ac{LJ} energy landscape towards the \ac{SHS}
energy landscape as $(m,n)$ increase.

Using a more realistic \ce{eLJ} potential obtained from coupled cluster
calculations for the xenon dimer
\autocite{Schwerdtfeger_ExtensionLennardJonespotential_2006,Jerabek_relativisticcoupledclusterinteraction_2017}
leads to $\mathcal{M}$ values close to those obtained for the
(6,12)-\ac{LJ} potential, but the results indicate that the topology of
the energy hypersurface is very sensitive to the model potential applied.  

Finally, the optimisation results have been compared to the previously published
results for the (6,12)-\ac{LJ} potential. The mapping from
$\mathcal{M}_\text{SHS}$ to $\mathcal{M}_\mathrm{SHS\to LJ}$ is non-injective
and non-surjective, however, the number of structures missed by the optimisation
procedure is relatively small. The unmatched structures belong to the high
energy region of the potential energy hypersurface and possess rather large
variations in their bond lengths. An analysis of their geometries revealed that
most of the larger structures can be constructed from a smaller cluster by
capping some of the triangular faces. This procedure effectively sorts almost
all unmatched structures into six seeds for clusters up to $N=13$. 

Further investigations should focus on comparing the presented results with
other commonly used interaction potentials like the Morse potential. It would
also be interesting to see how three-body interactions would affect the results.
However, that would require a considerable amount of programming as the current
program was developed with only two-body forces in mind.

\acresetall

\chapter[The Gregory-Newton Clusters]{The Gregory-Newton Clusters\footnote{This
    chapter is composed of sections previously published in the articles
    \citetitle*{Trombach_stickyhardsphereLennardJonestypeclusters_2018}\autocite{Trombach_stickyhardsphereLennardJonestypeclusters_2018}
    and
    \citetitle*{Trombach_GregoryNewtonproblemkissing_2018}\autocite{Trombach_GregoryNewtonproblemkissing_2018}
    and is reprinted with permission from the publisher \textcopyright 2018
    American Physical Society. Some sections have been modified to fit the style
    of this thesis.}}
\label{sec:thegregorynewtonclusters}

\section{Introduction}
\label{sec:introductiongregorynewton}

In 1930, Tammes studied the distribution of pores on pollen grains, which
required him to find a solution to the problem of packing a number of circles
(or spheres) on the surface of a unit sphere, maximising their
distance.\autocite{Tammes_originnumberarrangement_1930} In graph theoretical
terms, in which the centers of all the circles correspond to the vertices of a
convex polyhedron, one tries to find the graph representing the polyhedron that
maximises the shortest edge lengths, while keeping the distance to the center of
the polyhedron fixed. Exact solutions to this problem are available for cluster
sizes of $3 \leq N \leq 14$ and
$N=24$.\autocite{Robinson_Arrangement24points_1961,Musin_TammesProblem14_2015}

A related problem that goes back to an argument between Newton and Gregory, one of his
apprentices, is about the \emph{maximum kissing number}
or \emph{Newton number} $N_k(d)$ of three-dimensional unit spheres ($d=3$) that
can simultaneously touch a central sphere of the same
size.\autocite{Pfender_KissingNumbersSphere_2004} While Newton believed
$N_k(3)=12$, Gregory thought a 13th sphere could be brought into contact with
the central sphere. It turns out that Newton was right, which was first proven
in 1953 by
\citeauthor{Schutte_ProblemdreizehnKugeln_1952}\autocite{Schutte_ProblemdreizehnKugeln_1952}
In the following a cluster with $N\geq 13$ and at least 12 spheres in contact
with a central sphere will be referred to as a \acf{GNC}. Figure~\ref{fig:GN}
shows the most symmetric icosahedral solution to the Gregory-Newton problem.
%
\begin{figure}[htb]
    \centering
    \includegraphics[width=.4\textwidth]{gregory-newton/Kissing-3d.png} \quad
    \includegraphics[width=.5\textwidth]{gregory-newton/ico.pdf}
    \caption{Left: Symmetric realization of $N_k(3)=12$ for unit hard spheres
    (icosahedral symmetry, $I_h$). The minimum distance between the outer
    spheres is $r=\sin^{-1}\left(\frac{2\pi}{5}\right)=1.05146222\dots$,
    hence they do not touch. Right: The corresponding icosahedral graph.
    Numbering refers to the respective node index.}
    \label{fig:GN}
\end{figure}
%
The Gregory-Newton problem has been solved for dimensions 1--9 and 24 in lattice
packings and 1--4, 8 and 25 for non-lattice
packings.\autocite{Musin_kissingnumberfour_2008,Conway_SpherePackingsLattices_1999,Musin_proof24cellconjecture_2018}
Lower and upper bounds have also been
published.\autocite{Mittelmann_Highaccuracysemidefiniteprogramming_2010,Conway_SpherePackingsLattices_1999}
However, for the problem of cluster nucleation, the three-dimensional problem,
as posed by Gregory and Newton, is more relevant and has recently been reviewed
by
\citeauthor{Kusner_ConfigurationSpacesEqual_2018}\autocite{Kusner_ConfigurationSpacesEqual_2018}

There are other problems similar to the Tammes problem like the Thomson problem,
that tries to find the optimal solution for charged particles (e.g. electrons)
on the surface of a
sphere.\autocite{Wales_Structuredynamicsspherical_2006,Wales_Defectmotifsspherical_2009}
The solutions to this class of problems all show icosahedral symmetry for a
size of 12 spheres. When thinking of these problems in terms of nucleation
phenomena there is usually not a repulsive force between the surrounding
spheres, but an attractive one. For example, such a system could be modelled in
a gravitational potential $V_G=Gm_im_jr_{ij}^{-2}$, where $G$ is the gravitational
constant, $m_i$ the mass of sphere $i$ and $r_{ij}$ the distance between spheres
$i$ and $j$. For spheres $i$ and $j$ with radii $R_i$ and $R_j$ a geometry
optimisation can be carried out under the constraint $r_{ij}\geq (R_i + R_j)$,
which will remove the gaps between the spheres in the icosahedral arrangement.
Such problems are often studied for crystallisation and sedimentation
phenomena.\autocite{Levin_Crystallizationhardspheres_2000,Pusey_Hardspherescrystallization_2009}

Two other potentials that enforce theses types of rigidity constraints on the
system are the \acf{LJ} and \acf{SHS} potentials already introduced in
chapters~\ref{sec:energylandscapes} and
\ref{sec:fromstickyhardspheretoLJtypeclusters}. Unlike in
chapter~\ref{sec:fromstickyhardspheretoLJtypeclusters} the exponents $(a,b)$ are
now any real positive number instead of integers.
%
\begin{equation}
    V_{a,b}^\mathrm{LJ}(r)=\frac{ar^{-b}-br^{-a}}{b-a} \quad (\mathrm{with} \ \ r,a,b \in \mathbb{R}^+ \ \ \mathrm{and} \ \ b>a).
\label{eqn:abpot}
\end{equation}
%
The \ac{SHS} potential is employed unchanged.
%
\begin{align}
    \lim_{a,b\rightarrow \infty} V_{a,b}^\mathrm{LJ}(r) \rightarrow V_\mathrm{SHS}(r)=
    \begin{cases}
        \infty   & r < 1\\
        -1  \quad \mathrm{for} & r = 1\\
        0       & r > 1
    \end{cases}
\label{eqn:KS1}
\end{align}
%
Note that for both potentials both the depth of the energy well and the
equilibrium distance have been arbitrarily set to 1.

In the following sections these two potentials will be used to investigate the
\acp{GNC}. Again, the investigation starts from the \ac{SHS} clusters derived
from the exact enumeration
method.\autocite{Hoy_Structuredynamicsmodel_2015,Holmes-Cerfon_EnumeratingRigidSphere_2016,Holmes-Cerfon_StickySphereClusters_2017}
This set of structures is searched for clusters fulfilling the requirements to
be \acp{GNC}. \ac{LJ} clusters of this size have an icosahedral global minimum,
however, this structure is impossible to realise in the \ac{SHS} potential. This
is due to the fact that in a perfect icosahedron the distance between the
vertices is always larger than the distance of all the vertices to the centre of
mass of the cluster. Therefore, as increasing \ac{LJ} exponents make the
potential more \ac{SHS}-like, there must be a point of symmetry breaking at which the icosahedral
structure cannot be supported any more by the \ac{LJ} potential energy surface.

%The arrangement of $N$ points on the surface of a sphere corresponding to the
%placement of $N$ identical spheres around a central sphere is called a
%spherical packing. To achieve optimal packings for spheres is known as the
%Tammes problem, originally posed in 1930 to study the distribution of pores on
%pollen grains \autocite{Tammes_originnumberarrangement_1930}: It is to
%determine the largest diameter and distribution that $N$ equal spheres may have
%when packed onto the surface of a sphere of radius 1 (unit sphere), without
%allowing for any overlap of these spheres. Alternatively, if the centre of each
%sphere is considered as the vertex of a polyhedron, the graph theoretical
%problem is to find the polyhedron that maximizes the shortest edge lengths with
%fixed distance to the central vertex.

%Newton and Gregory argued about the maximum possible number $N_k(d)$ (the
%\textit{maximum kissing number} or \textit{Newton number}) of three-dimensional
%unit spheres ($d=3$) that could be brought into contact with a central sphere
%\autocite{Pfender_KissingNumbersSphere_2004}. Sch\"utte and van der Waerden
%provided the first proof in 1953 that max$\{N_k(3)\}=12$
%\autocite{Schutte_ProblemdreizehnKugeln_1952}. Such a cluster of 12 unit
%spheres kissing a central unit sphere is called a \ac{GNC}, and is shown in its
%most symmetric icosahedral form in figure~\ref{fig:GN}. Exact Newton numbers
%for unit spheres in lattice packings are known for dimensions $d=1$ to 9 and
%$d=24$, and for non-lattice packings for $d=1-3$, 8 and 24
%\autocite{Conway_SpherePackingsLattices_1999,Musin_proof24cellconjecture_2018}.  Lower and %upper
%bounds for max$\{N_k(d)\}$ are also available
%\autocite{Mittelmann_Highaccuracysemidefiniteprogramming_2010,%Conway_SpherePackingsLattices_1999}.

%The Tammes, Thomson  or related models employ repulsive forces between points
%or spheres \autocite{Wales_Structuredynamicsspherical_2006,Wales_Defectmotifsspherical_2009} %and, for the
%three-dimensional problem with 12 kissing spheres, lead to ideal icosahedral
%symmetry (figure~\ref{fig:GN}). We may however pose the question of what
%happens if we let the outer kissing spheres of a \ac{GNC} touch each other to
%enforce rigidity? We could try to find the global and all local minima for the
%12 kissing hard spheres interacting through an attractive instead of a
%repulsive potential. For example, we can place the central hard sphere in a
%gravitational field of strength $F_G=Gm_im_jr_{ij}^{-2}$ and relax all
%positions $r_{ij}\ge (R_i+R_j)$ between the kissing hard spheres $i$ and $j$,
%in the most general case having sphere radii $R_i$ and masses $m_i$. It is
%clear that such a procedure leads to a less flexible and more rigid sphere
%packing. In Euclidian space, this problem is well known to
%crystallization/sedimentation phenomena modelled by hard spheres in a
%gravitational field \autocite{Levin_Crystallizationhardspheres_2000,%Pusey_Hardspherescrystallization_2009}.

%The most widely used interaction potential in chemical and physical sciences is
%the so-called Lennard-Jones (LJ) potential
%\autocite{Jones_DeterminationMolecularFields_1924,Lennard-Jones_Cohesion_1931},
%which is attractive in the long range and repulsive in the short range. In
%reduced units the LJ$(a,b)$ potential takes the form, 
%%
%
%%
%Such a potential maximizes the number of contacts between spheres, and for the
%famous LJ(6,12) case leads to one and only one minimum for the \ac{GNC}
%\autocite{Trombach_stickyhardsphereLennardJonestypeclusters_2018} -- the ideal icosahedron %(shown in figure~\ref{fig:GN}) as
%in the case of the Tammes problem. The icosahedral motif originally proposed by
%Mackay \autocite{Mackay-1962} plays a very important role in cluster physics and
%chemistry
%\autocite{Hoare_Physicalclustermechanics_1975,Klots90,Uppenbrink-1991,vandewaal93,%Wales_Whatcancalculations_1996,Wales_ChangesMorphologyCapping_1996,%Wales_Structuredynamicsspherical_2006}.
%
%For large exponents $(a,b)$ the LJ potential approaches the sticky hard-sphere
%(SHS) limit originally introduced by Baxter \autocite{baxter68,%Gazzillo_AnalyticsolutionsBaxter_2004},
%
%SHS models have been used intensively in many areas, such as crystallization,
%flocculation, colloidal suspensions, micelles, protein solutions, or in the
%exact enumeration of rigid SHS clusters
%\autocite{Stell_Stickyspheresrelated_1991,Jamnik_Spatialcorrelationssolvation_1996,%Santos_Radialdistributionfunctions_1998,Gazzillo_AnalyticsolutionsBaxter_2004,%Hoy_MinimalEnergyPackings_2010,Arkus_Minimalenergyclusters_2009,Arkus-2010,%Arkus_DerivingFiniteSphere_2011,Hoy_Structurefinitesphere_2012,%Hayes_ScienceStickySpheres_2012,Holmes-Cerfon_geometricalapproachcomputing_2013,%Holmes-Cerfon_EnumeratingRigidSphere_2016,Holmes-Cerfon_StickySphereClusters_2017,%Kallus_Freeenergysingular_2017}.
%The SHS model has the advantage that an adjacency matrix $A$ can be introduced
%with entries $A_{ij}=1$ if spheres $i$ and $j$ touch, and 0 otherwise. The
%number of contacts between spheres then simply becomes $N_c=\sum_{i<j}^N
%A_{ij}$. It also opens the way for graph-theoretical treatments of cluster
%structures as we shall see.

%A complete set of non-isomorphic rigid SHS clusters has recently been
%determined for cluster size $N \leq 14$ via exact enumeration studies employing
%geometric rejection
%\autocite{Hoy_Structuredynamicsmodel_2015,Holmes-Cerfon_EnumeratingRigidSphere_2016}
%and rigidity rules \autocite{Holmes-Cerfon_StickySphereClusters_2017}. These
%include the rigid \acp{GNC} as a subset. Recently, Trombach et al. determined
%that the number of non-isomorphic rigid \acp{GNC} $n_{\mathrm GNC}$ is
%surprisingly large ($n_{\mathrm GNC}=737$)
%\autocite{Trombach_stickyhardsphereLennardJonestypeclusters_2018}. An
%explanation for this has not been given yet. Moreover, the condition
%$\lim_{a,b\rightarrow \infty} V_{a,b}^\mathrm{LJ}(r) \rightarrow
%V_\mathrm{SHS}(r)$ implies that at certain $a,b$ values symmetry broken
%solutions away from the ideal icosahedral structure appear. Where exactly this
%happens, and when the icosahedral structure does not survive anymore, is not
%known. In order to close this gap, we decided to analyze the rigid \acp{GNC} and
%corresponding symmetry breaking effects in more detail. This is much in the
%spirit of Wales who already pointed out that the global characteristics of the
%energy landscape can be quite sensitive to the nature of the interatomic
%potential applied \autocite{Wales_MicroscopicBasisGlobal_2001}.

\section{Computational Details}

Coordinates for \ac{GNC} structures have been obtained by searching for
adjacency matrices of the results for $N=13$ from
\citeauthor{Holmes-Cerfon_EnumeratingRigidSphere_2016}
\autocite{Holmes-Cerfon_EnumeratingRigidSphere_2016} with one row or column
containing twelve ``1'' entries. Sub-graph isomorphism was verified using the
\textit{VF2} algorithm \autocite{Cordella_SubGraphIsomorphism_2004} as
implemented in the \textit{boost graph library}
\autocite{Siek_BoostGraphLibrary_2002} using the program package
\textsc{Spheres}. Structural optimisations with \ac{LJ} potentials have been
carried out using the multidimensional function minimiser from the \Cpp library
\textit{dlib} \autocite{King_DlibmlMachineLearning_2009} and an energy
convergence criterion of $10^{-15}$. Results from the optimisation procedure
were analysed based on the \ac{EDM}, which is unique for non-isomorphic
structures apart from permutation, translation, rotation and inversion. For this
the distances were sorted lexicographically.


\section{The Gregory-Newton Problem for Soft Potentials}
\label{sec:thegregorynewtonproblemforsoftpotentials}

The question of the Newton number in three dimensions has been resolved almost
70 years ago\autocite{Schutte_ProblemdreizehnKugeln_1952}. The proof is valid
for hard-sphere short-range potentials, but little is known about the behaviour
of such clusters under long-range potentials such as the Kratzer
potential\autocite{Kratzer_ultrarotenRotationsspektrenHalogenwasserstoffe_1920}.
For unequally sized spheres, some simple results are known; for example, 13 hard
spheres of radius $r_s$ can touch a central sphere of unit radius only if
$r_s\leq 0.9165$.\autocite{phillips12} For clusters bound by the aforementioned
long-range potentials the situation is far more complicated as it requires to
minimise energy rather than distances between neighbouring particles.
Nonetheless, it is important to expand our knowledge on these kind of systems as
they are crucial to understand real systems such as coordination compounds,
which have recently been shown to possess coordination numbers as high as
17\autocite{Kaltsoyannis-2017} or even 20.\autocite{Suresh-2016} 

The optimisation procedures explained in chapter~\ref{sec:theprogramspheres}
were used to minimise the energy of a starting structure consisting of 13
spheres surrounding a center sphere with a fixed distance of one. Generating
such a starting structure where all surrounding spheres are evenly spaced is
impossible since there exists no triangulation of a sphere with 13 vertices,
where every vertex has degree five or
six\autocite{Schwerdtfeger_topologyfullerenes_2015}. To generate an approximate
distribution the Fibonacci sphere
algorithm\autocite{Gonzalez_MeasurementAreasSphere_2010,Keinert_SphericalFibonacciMapping_2015}
was used and the generated structure of size $N=14$ was the starting point for
optimisations with \ac{LJ} potentials with small exponents. The difference
between the largest and smallest \ac{COS} distance was used as a measure for
whether the 13th sphere enters the first coordination shell. A value of zero
would be expected for this to be true.

The results for all positive integer combinations of $m\leq11$ and $n\leq12$
with $m<n$ are depicted in figure~\ref{fig:gregorynewton-N14}. 
%
\begin{figure}[htb]
    \centering
    \includegraphics[width=.8\textwidth]{gregory-newton/N14.pdf}
    \caption{Relation of \acs{LJ} exponents $m$ and $n$ to the difference of
    largest and smallest \acs{COS} distances.  A value of zero would imply that
    all surrounding spheres are touching the center sphere.}
    \label{fig:gregorynewton-N14}
\end{figure}
%
Even for the combination of smallest exponents (1,2) it is clear that the
\ac{COS} distances vary from sphere to sphere. For this potential the largest
\ac{COS} distance is $r_\text{max}=0.882$, while the shortest one is
$r_\text{min}=0.804$. While the longest distance only shows up once, the
shortest distance appears twice. All other ten distances fall in the range
between $r = 0.845$ and $r = 0.861$. The $r_\text{max} /r_\text{min}$ ratio is
$1.097$ and much shorter compared to $r_\text{max} /r_\text{min}= \sqrt{2}$ for
the closed packed lattice, or the shortest distance possible for the \ac{SHS}
system which is $r_{14}^\text{GN} = 1.347$ (see
section~\ref{sec:addinga14thsphere}). Hence, the 13th sphere almost touches the
center sphere.

Note that all \ac{COS} distances for the $N=14$ (1,2)-\ac{LJ} cluster are
significantly shorter than $r=1$, due to the $N(N-1)/2$ attractive two-body
interactions and the softness of the potential.  For infinite (e.g.
body-centred cubic or close-packed) lattices of particles interacting via
$V^\mathrm{LJ}_{mn}(r)$ with $n> m >3$, one can prove
\autocite{Schwerdtfeger_ExtensionLennardJonespotential_2006} that the nearest
neighbour distance is
%
\begin{equation}
    r_\mathrm{NN}(m,n)=\left( L_n L_m^{-1}\right)^\frac{1}{n-m}. 
    \label{eqn:lattice}
\end{equation}%
%
Here $L_n$ is the Lennard-Jones-Ingham lattice coefficient for a specific
lattice determined from 3D lattice sums.  Since $L_n<L_m$ for $n>m$, we see
that $r_\mathrm{NN}<1$, and $\lim\limits_{m,n\rightarrow
\infty}r_\mathrm{NN}(m,n)=1$.  The shortest distances found in (6,12)-\ac{LJ}
clusters $r_\text{min}(N)$ are: $r_\text{min}(8)=0.986767$,
$r_\text{min}(9)=0.964404$, $r_\text{min}(10)=0.964382$,
$r_\text{min}(11)=0.956345$, $r_\text{min}(12)=0.947842$, and
$r_\text{min}(13)=0.952179$.  Surprisingly, $r_\text{min}(12)$ is smaller than
$r_\mathrm{NN}(6,12)$ for typical crystalline lattices; $r_\mathrm{NN}(6,12)$
values are $0.95066$, $0.95186$ and $0.97123$ for simple cubic, body-centred
cubic and close-packed lattices, respectively.  This result shows that stable
clusters do not necessarily have longer bonds compared to the solid state.

\section{Rigid Gregory-Newton Clusters and Corresponding Graphs}

The recent results by Holmes-Cerfon contain a putatively complete set of rigid
SHS clusters of size $N=13$ and $N=14$
\autocite{Holmes-Cerfon_EnumeratingRigidSphere_2016}. The rigid \acp{GNC} can
easily be identified as a subset of the set of all non-isomorphic rigid SHS
clusters, i.e. $\{ S_\mathrm{GN}\}\subset \{ S_\mathrm{SHS}\}$; these have
adjacency matrices $A$ with exactly one column and row containing twelve "1"
entries due to 12 spheres kissing the central sphere. A surprisingly large
number of 737 non-isomorphic $N = 13$ \acp{GNC} out of 98,540 rigid SHS clusters
can be found \autocite{Trombach_stickyhardsphereLennardJonestypeclusters_2018}.
There are four different possible contact numbers $N_c$ with $\{724,10,1,2\}$
rigid \ac{GNC} corresponding to $N_c=\{33,34,35,36\}$
\autocite{Robinson_Arrangement24points_1961}; therefore, none of those clusters
are hypostatic.

For further analysis and without loss of generality the central sphere was
removed and the remaining non-isomorphic shell of spheres\footnote{Note that
rigidity requires the presence of the central sphere.} was analysed, also called
a contact graph according to
\citeauthor{Schuette_AufwelcherKugel_1951}\autocite{Schuette_AufwelcherKugel_1951}
This has the advantage that these shells are related to planar connected graphs.
In the following the corresponding connected planar graph of such a shell of
spheres with the central sphere missing will be referred to as a \textit{GN
graph}. The question arises if all 737 non-isomorphic GN graphs are sub-graphs
of the icosahedral graph, as shown in figure~\ref{fig:GN}. This would make sense
as it is impossible to increase the degree of any vertex beyond five in
the \ac{GN} graph. Note that the icosahedral cluster is completely unjammed and
its space of (infinitesimal) deformations has dimension
24.\autocite{Kusner_ConfigurationSpacesEqual_2018}

Employing the \textit{VF2} algorithm
\autocite{Cordella_SubGraphIsomorphism_2004} as implemented in the \textit{boost
graph library} \autocite{Siek_BoostGraphLibrary_2002} all 737 non-isomorphic
\ac{GN} graphs $G_\mathrm{GN}(N,E')$ (vertex count $|N|=12$, edge count
$|E'|<30$) are found to be (edge-induced) sub-graphs of the icosahedral graph
$G_\mathrm{ico}(N,E)$ ($|N|=12$, $|E|=30$), which implies that their vertices
can all be mapped to vertices of the icosahedral graph with certain edges
deleted such that the sub-graph remains connected ($N_\mathrm{GN} =
N_\mathrm{ico}$ and $E_\mathrm{GN}\subset E_\mathrm{ico}$). An extensive list of
all sub-graphs is included in the appendix (Tables~\ref{tab:icosubgraphs} and
\ref{tab:verticesandfaces}). Note, not all GN graphs are 3-connected and
therefore are not strictly polyhedral according to Steinitz's theorem
\autocite{Steinitz_PolyederundRaumeinteilungen_1916}. These are the graphs which
have vertices of degree 2, i.e. $|N_2|>0$, and there are 304 of them,
table~\ref{tab:verticesandfaces}. As the many non-isomorphic graphs listed in
the appendix are obtained from a certain combination of edge deletions under the
constraint of maintaining rigidity, it is not surprising at all that the number
of non-isomorphic \ac{GN} graphs is so large.

\begin{figure}[htb]
    \centering
    \subfloat[hcp, $|E|=24, \omega =1$.\label{subfig:hcpgraph}]{\includegraphics[width=0.43\columnwidth]{gregory-newton/hcp-new.pdf}\hspace{0.03\textwidth}\includegraphics[width=.37\columnwidth]{gregory-newton/hcp.png}}\\
    \subfloat[fcc, $|E|=24, \omega =2$.\label{subfig:fccgraph}]{\includegraphics[width=0.43\columnwidth]{gregory-newton/fcc-new.pdf}\hspace{0.03\textwidth}\includegraphics[width=.37\columnwidth]{gregory-newton/fcc.png}}
    \caption{\acs{GN} \acs{hcp} (triangular orthobicupola) and \acs{fcc}
    (cuboctahedron) graphs (central sphere removed) as sub-graphs of the
    icosahedral graph and corresponding rigid \acsp{GNC}.  Red lines indicate
    the edges that were removed to create the \acs{GN} graph. The ordinal
    numbers $\omega$ refer to Table~\ref{tab:icosubgraphs} in the appendix.}
    \label{fig:GNshellgraphs}
\end{figure}

The results show, that at least six and up to a maximum of nine edges have to be
removed from the icosahedral graph to create a \ac{GN} graph. Removing six edges
from the icosahedral graph results in 24 edges, or $N_c=36$ if the central
sphere is included. For $N=13$ this is exactly equal to $3N-3$ which is the
maximum contact number observed for this cluster size
\autocite{Hoy_Structuredynamicsmodel_2015,Holmes-Cerfon_EnumeratingRigidSphere_2016}.
Consequently, removing nine edges gives $N_c=33=3N-6$, meaning that rigid
\acp{GNC} cannot be hypostatic (i.e. $N_c < 3N-6$). Interestingly, there are
only two graphs with maximum edge count of $|E|=24$, which are exactly the
fragments of the \acf{fcc} and \acf{hcp} bulk structures, respectively. These
are the result from removing 6 edges in such a way, that exactly one edge is
removed from every vertex in the icosahedral graph (thus the degree of every
vertex is 4), see figure~\ref{fig:GNshellgraphs}. Removing edges in this way
implies that the resulting two graphs consist of triangles and rectangles only.
The difference between the \ac{fcc} and \ac{hcp} clusters is in the way their
square faces are connected; in the \ac{fcc} case the square faces only connect
via vertices (cuboctahedron), while in the \ac{hcp} case the square faces come
in pairs sharing one edge (triangular orthobicupola or Johnson solid
$J_{27}$)\autocite{Kusner_ConfigurationSpacesEqual_2018}.

The construction of hcp and fcc structures by a continuous deformation of an
icosahedron has been described in detail by Kusner et al.
\autocite{Kusner_ConfigurationSpacesEqual_2018} and goes back to Conway and
Sloane in 1988 \autocite{Conway_SpherePackingsLattices_1999}. \ac{hcp} and
\ac{fcc} can both be obtained from a rearrangement of the spheres in an
icosahedron by forming a (zig-zag) cycle (closed path) through six vertices, and
arranging those spheres on the path such that they are in-plane with the central
sphere, which becomes part of the hexagonal plane as in the bulk \ac{fcc} and
\ac{hcp} packing (figure~\ref{fig:ico-fcc-trans}).
%
\begin{figure}[htb]
    \centering
    \includegraphics[width=.8\columnwidth]{gregory-newton/plane.png}
    \caption{Illustration of one zig-zag path (light blue spheres) that needs to
    be deformed such that it aligns with the triangular plane (shown in grey) of
    the \acs{fcc} crystal.}
    \label{fig:ico-fcc-trans}
\end{figure}
%
Additionally, the plane has to be rotated by $\pi/6$ to create the \ac{fcc}
structure. The \ac{hcp} structure can be constructed by also rotating either the
top or the bottom plane by the same amount in either direction parallel to the
hexagonal plane. Kusner noted that a smooth deformation from the icosahedral
configuration to hcp requires 9 moving spheres
\autocite{Kusner_ConfigurationSpacesEqual_2018}. This interesting transition
path may be the key for the icosahedral to closed-packed rearrangements in
larger clusters, which has previously been described in terms of catastrophe
theory as a cusp catastrophe \autocite{Wales_MicroscopicBasisGlobal_2001}. 

Even though the rearrangement from the icosahedral to either the \ac{fcc} or
\ac{hcp} cluster structure can easily be realised for the \ac{GNC}, there should
be clusters where the icosahedral motif is still clearly visible, i.e. only
small rearrangements of the spheres are necessary to break icosahedral symmetry
and form a rigid cluster. These are, for example, the ones with maximum count of
triangles, i.e. according to table~\ref{tab:verticesandfaces} the GN graphs with
$|F_3|=10$ with edge counts of $|E|=$ 22 or 21. Two of these are shown with
their corresponding graphs in figure~\ref{fig:GNicographs}.
%
\begin{figure}[htb]
    \centering
    \subfloat[icosahedral motif, $|E|=22, \omega =4$.\label{subfig:ico4graph}]{\includegraphics[width=0.47\columnwidth]{gregory-newton/4-22-new.pdf}\hspace{0.03\textwidth}\includegraphics[width=.3\columnwidth]{gregory-newton/GNico4-22}}\\
    \subfloat[icosahedral motif, $|E|=22, \omega =7$.\label{subfig:ico7graph}]{\includegraphics[width=0.47\columnwidth]{gregory-newton/7-22-new.pdf}\hspace{0.03\textwidth}\includegraphics[width=.3\columnwidth]{gregory-newton/GNico7-22.png}}
    \caption{Representative \acs{GN} graphs (central sphere removed) with
    $|F_3|=10$ as sub-graphs of the icosahedral graph and corresponding rigid
    \acsp{GNC}. The icosahedral motif in the 3D embedding is clearly visible.
    Red lines indicate the edges that were removed to create the \acs{GN} graph.
    The ordinal numbers $\omega$ refer to Table~\ref{tab:icosubgraphs} in the
    appendix.}
    \label{fig:GNicographs}
\end{figure}

Figure~\ref{fig:GNJohnsongraph} shows the graph with the next highest edge count
after the \ac{fcc} and \ac{hcp} packings. 
%
\begin{figure}[htb]
    \centering
    \subfloat[Distorted elongated pentagonal bipyramid (Johnson solid $J_{16}$), $|E|=23, \omega =3$.\label{subfig:johnsongraph}]{\includegraphics[width=0.47\columnwidth]{gregory-newton/johnson-new.pdf}\hspace{0.03\textwidth}\includegraphics[width=.3\columnwidth]{gregory-newton/johnson.png}}
    \caption{\acs{GN} graph (central sphere removed) as sub-graphs of the
    icosahedral graph and corresponding \acs{GN} Johnson-like solid (with edges
    removed). Red lines indicate the edges that were removed from the
    icosahedral graph to create the \acs{GN} graph. The ordinal number $\omega$
    refers to Table~\ref{tab:icosubgraphs} in the appendix.}
    \label{fig:GNJohnsongraph}
\end{figure}
%
The motif of a distorted elongated pentagonal bipyramid (Johnson solid $J_{16}$)
is clearly visible. Note that the Johnson solid can be obtained by deleting five
edges and rotating the two opposite pentagonal pyramids by $2\pi /5$. One of the
resulting square faces has to be stretched to obey the SHS conditions, which is
achieved by removing two additional edges. In the graph this implies that a
hexagonal face is formed. Note that this \ac{GNC} is also the cluster with the
largest distance $r_\mathrm{max}^\mathrm{RE}= 1.47823719$ that corresponds to a
removed edge (RE) in the \ac{GN} graph. Capping this cluster with one more
sphere over the distorted square face with $r_\mathrm{max}^\mathrm{RE}$ leads to
the structure with the shortest distance a sphere in the second coordination
shell can have to the central sphere
($r^\mathrm{COS}=1.347150628$)\footnote{This will be investigated in detail in
section~\ref{sec:addinga14thsphere}.} out of all 895,478
\ac{GN} clusters with $N=14$
\autocite{Trombach_stickyhardsphereLennardJonestypeclusters_2018}.

If more edges are removed from the icosahedral graph larger $n$-gonal faces
appear, with the largest face being a 12-gon.
%

\section{Symmetry-Broken Lennard-Jones Gregory-Newton Clusters}

All 737 non-isomorphic rigid \acp{GNC} optimise to the ideal icosahedral
symmetry if a (6,12)-\ac{LJ} potential is applied
\autocite{Trombach_stickyhardsphereLennardJonestypeclusters_2018} (however, for
larger sized icosahedral structures many more minima
appear).\autocite{Doye-1995,Wales_TheoreticalPredictionsStructure_1996,Doye_effectrangepotential_1996,Doye_Structuralconsequencesrange_1997}
As mentioned in the introduction, for equally sized hard spheres a cluster with
icosahedral symmetry leaves gaps between the spheres of the outer shell, i.e.
they do not touch, and it is therefore not considered rigid under SHS
conditions. Hence, at certain $(a,b)$ combinations a phase transition must occur
in the $(a,b)$-\ac{LJ} energy landscape where non-icosahedral local minima
appear. In order to determine those $(a,b)$ combinations, all 3D cluster
geometries were optimised with varying exponents $(6\leq a \leq 34$ and $7\leq
b\leq 35)$ with $(b>a)$ and the number of resulting minimum structures was
analysed. The results are shown schematically in figure~\ref{fig:ico-2d}.
%
\begin{figure}[htb]\centering
    \includegraphics[width=.8\columnwidth]{gregory-newton/ico-2d.pdf}
    \caption{Number of unique structures resulting from an optimisation with a
    $(a,b)$-\acs{LJ} potential. The lowest contour line shows the point where more
    than one structure results from the optimisation and the distance between
    contour lines is 1.}
    \label{fig:ico-2d}
\end{figure}

Another interesting limiting case of the \ac{LJ} potential with exponents $a\to
0$ and $b\to \infty$, resembling a constant attractive potential with an
infinite wall, should be mentioned. In such a potential the kissing spheres can
move freely in the available space without change of energy. Indeed, in the
region of low $a$ and high $b$ values and increasing number of unique structures
is found. For example, values of $a=0.6$ and $b=120.0$ result in two distinct
structures that are both derived from the icosahedral motif.

Figure~\ref{fig:no-ico} contains additional information showing three major
phase transitions in the topology of the energy landscape going from low to high
$(a,b)$ exponents. 
%
\begin{figure}[htb]\centering
    \includegraphics[width=.8\columnwidth]{gregory-newton/no-ico.pdf}
    \caption{Different types of energy landscapes arising from combinations of
    the $(a,b)$-\acs{LJ}  exponents. (1) One single (icosahedral) minimum, (2) more
    than one minimum with the icosahedron as the global minimum, (3) more than
    one minimum with the icosahedron becoming a local (and not global) minimum,
    (4) the icosahedral motif disappears completely. The unshaded small area
    in the bottom right corner corresponds to $a>b$, which is excluded. The
    resolution for $a$ is 1.0 and for $b$ 0.25.}
    \label{fig:no-ico}
\end{figure}
%
In the blue shaded area (1), the Mackay icosahedron is the sole minimum in the
potential energy landscape. The first transition occurs when this symmetry can
be broken, and other local minima are supported by the $(a,b)$-\ac{LJ} potential
besides the icosahedron. This is indicated in figure~\ref{fig:no-ico} by the
smallest, orange region (2), which still contains the perfect icosahedron as the
global minimum. At slightly higher exponents, other structures become
energetically more favourable and replace the icosahedron as the global minimum,
region (3). However, the icosahedron remains as a local minimum in the potential
energy surface. The last transition occurs when the \ac{LJ} potential becomes
\ac{SHS}-like, and the icosahedral cluster completely disappears from the
potential energy surface, region (4). The three transition lines are generally
smooth.

Figure~\ref{fig:compareLJ} shows representative \ac{LJ} potentials for
combinations of the $(a,b)$ exponents (with low and high $a$ values) on the
phase transition lines drawn in figure~\ref{fig:no-ico}. 
%
\begin{figure}[htb]\centering
    \includegraphics[width=.8\columnwidth]{gregory-newton/compareLJ.pdf}
    \caption{Comparison of different shapes of \acs{LJ} potentials at the phase
    transition lines shown in figure~\ref{fig:no-ico} with the traditional
    (6,12)-\acs{LJ} potential (black solid line). Dashed lines refer to
    potentials with low $a$ values (left side of figure~\ref{fig:no-ico}), while
    solid lines refer to potentials with high $a$ values (right side of
    figure~\ref{fig:no-ico}).}
    \label{fig:compareLJ}
\end{figure}
%
At these phase transition lines, the corresponding \ac{LJ} potentials show
narrow and steep repulsive potentials compared to the $(6,12)$-\ac{LJ}
potential, which all look very similar in the short range ($r<1$). However,
they differ substantially in the long range ($r>1$).

The $(a,b)$ parameters can be related to the so-called \ac{LJ} hard-sphere radius
$\sigma$ (given by the intersection with the abscissa) through
equation~\eqref{eqn:nmpot},
%
\begin{equation}
    \sigma=\left(\frac{b}{a}\right)^{\frac{1}{a-b}}.
\end{equation}
%
and only the $(a,\sigma)$ combinations shown in figure~\ref{fig:hardsphere}
along the phase transition lines have to be considered.
%
\begin{figure}[htb]\centering
    \includegraphics[width=.8\columnwidth]{gregory-newton/sigma.pdf}
    \caption{Hard-sphere radii $\sigma$ in reduced units for the $(a,b)$-\acs{LJ}
    potentials on the transition lines shown in figure~\ref{fig:no-ico}.}
    \label{fig:hardsphere}
\end{figure}

The variation of $\sigma$ along the phase transitions lines for
(2)$\rightarrow$(3) and (3)$\rightarrow$(4) are rather small. However, all three
transitions clearly show different ranges for $\sigma$ and thus can be
characterized by the \ac{LJ} hard-sphere radius. These are also much larger
compared to the (6,12)-\ac{LJ} hard-sphere radius of $\sigma=0.891$, and close
to the ideal hard sphere radius of 1 in the \ac{SHS} model. This demonstrates
that the shape of the \ac{LJ} potential in the repulsive region has a
significant influence on the position of the transition lines, and therefore on
the topology of the energy landscape. In contrast, these transitions seem to be
far less affected by the shape of the potential in the attractive region. Only
for the transition (1)$\rightarrow$(2) a larger variation in $\sigma$ is
observed.

Finally, the results show that long-range interactions stabilise the icosahedral
cluster. Therefore, the assumption that second-nearest-neighbour interactions
may be important seems to come naturally. However, first-nearest neighbour
interactions are sufficient for stabilising this structure, i.e. if the \ac{GN}
clusters are optimised with a truncated (6,12)-\ac{LJ} potential that ignores
second-nearest-neighbour interactions by setting the range of interactions to
distances below 1.5, it is observed, that the icosahedron is recovered.

\section{Adding a 14th Sphere}
\label{sec:addinga14thsphere}

Finally, the set of \ac{SHS} clusters with $N=14$ from exact enumeration
results\autocite{Holmes-Cerfon_EnumeratingRigidSphere_2016} has been
investigated with respect to the existence of \acp{GNC}. The extra sphere in
these clusters is required to enter the second coordination shell because of
the Newton number. In this case, an even larger number of clusters exists
($14,529$), which is $\approx 0.016|\mathcal{M}_\mathrm{SHS}(14)|$. All of these
structures optimise to just one of two possible (6,12)-LJ minima of \ac{GN}
type. The first is the Mackay icosahedron capped at one of its triangular faces,
and the second is an elongated pentagonal bipyramid (belonging to the class of
Johnson solids) with the 14th sphere capping a square face.

Most of these $N=14$ clusters are minimally rigid ($N_c=3N-6=36$), while only a
few are hyperstatic ($N_c > 3N-6$) and none are hypostatic ($N_c < 3N-6$). There
are $\{14369,144,8,6,2\}$ such clusters with $N_c=\{36,37,38,39,40\}$ and
$N=14$.  The clusters with $N_c=40$ are \ac{hcp} and \ac{fcc} core-shell
structures capped at a square face; these arrangements maximise $N_c$. Most of
the clusters with $N_c=\{38,39\}$ are deformed versions of the elongated
pentagonal bipyramid mentioned above, indicating that this arrangement is a
favoured route to these intermediate-energy structures.  However, $N_c=39$ also
contains \ac{hcp} and \ac{fcc} structures capped at a triangular face.  The
first example of a cluster derived from a perfect icosahedral symmetry shows up
at lower value $N_c=37$.  Representative examples for clusters with high
contact numbers are depicted in figure~\ref{fig:N14}.  
%
\begin{figure}
    \centering
    \subfloat[$r_{14}^\text{GN}=1.34715$, $N_c=39$\label{subfig:short-greg-newton}]{\includegraphics[width=0.4\textwidth]{gregory-newton/short.png}}
    \subfloat[$r_{14}^\text{GN}=1.37515$, $N_c=36$\label{subfig:2ndshort-greg-newton}]{\includegraphics[width=0.4\textwidth]{gregory-newton/2ndshort.png}}\\
    \subfloat[$r_{14}^\text{GN}=\sqrt{2}$, $N_c=40$\label{subfig:sqrt2-greg-newton}]{\includegraphics[width=0.4\textwidth]{gregory-newton/sqrt2.png}}
    \subfloat[$r_{14}^\text{GN}=\sqrt{\frac{8}{3}}$, $N_c=39$\label{subfig:sqrt83-greg-newton}]{\includegraphics[width=0.4\textwidth]{gregory-newton/sqrt83.png}}\\
    \caption{Graphical representations of \acs{SHS} packings with $N=14$, where
    a center sphere is maximally contacting. The orange sphere in each cluster
    is the 14th outer sphere, not able to touch the center sphere (in black).
    \protect\subref{subfig:short-greg-newton} distorted elongated pentagonal
    bipyramid (Johnson solid); \protect\subref{subfig:2ndshort-greg-newton}
    distorted icosahedron; \protect\subref{subfig:sqrt2-greg-newton} hcp capped
    on a square; \protect\subref{subfig:sqrt83-greg-newton} hcp capped on a
    triangle.}
    \label{fig:N14}
\end{figure}

Surprisingly, the $N = 14$ cluster with the closest \ac{COS} distance
$r_\text{min}^\text{COS}$ was not known. Here, this gap is closed by determining
the \ac{COS} distance for all \ac{SHS} Gregory-Newton type clusters with $N=14$.
One single cluster with $r_\text{min}^\text{COS}=1.3471506281091$ is found.  Its
structure (figure~\ref{subfig:short-greg-newton}) is similar to the elongated
pentagonal bipyramid with one of the square faces stretched to form a regular
rectangle. The 14th sphere caps this deformed face, becoming the vertex of a
deformed octahedron and allowing the outer sphere to get closer to the central
sphere.  The next-smallest-$r^\text{COS}$ cluster ($r^\text{COS} = 1.37515$) is
shown in figure~\ref{subfig:2ndshort-greg-newton}. It does not belong to the
category of the clusters derived from the elongated pentagonal bipyramid, but
instead can be described as being icosahedral-like.  The short distance is
achieved by attaching the 14th sphere to 3 spheres that do not form a face of
the cluster (because they are separated by a distance larger than $1$.)

The distribution of $r^\text{COS}$ values for the full set of $N=14$ \ac{GN}
clusters is shown in figure~\ref{fig:greg-newton}.
%
\begin{figure}
    \centering
    \includegraphics[width=0.8\textwidth]{gregory-newton/greg-newton.pdf}
    \caption{Frequency of distances from the cluster center to the most distant
    sphere for all Gregory-Newton-like clusters with $N=14$ contained in the
    structures from
    \citeauthor{Holmes-Cerfon_EnumeratingRigidSphere_2016}\autocite{Holmes-Cerfon_EnumeratingRigidSphere_2016}.
    The width of the bars is $0.01$~$r_e$.}
    \label{fig:greg-newton}
\end{figure}
%
Motifs with larger $r^\text{COS}$ are far more prevalent. For example, the peak
at $r^\text{COS} = 1.41$ corresponds to structures where the 14th sphere is
touching 4 other spheres that are part of a tetragonal pyramid, therefore
forming a regular octahedron with a tip-to-tip distance of $\sqrt{2}$
(figure~\ref{subfig:sqrt2-greg-newton}).  The maximum $r^\text{COS}$ value
($1.63$) corresponds to capping triangular faces, so that the most distant
sphere is part of a regular trigonal bipyramid with a height of $\sqrt{8/3}$
(figure~\ref{subfig:sqrt83-greg-newton}). The structures corresponding to the
bars at $1.60,1.58$ and $1.55$ are derived from the regular trigonal bipyramid
and result from breaking its axial bonds. In these structures, the more bonds
are broken, or the further the axial spheres are separated, the shorter the
\ac{COS} distance becomes.

\section{Conclusion}

Rigid \acp{GNC} have been analysed by graph theoretical means. All 737
non-isomorphic \ac{GN} graphs are sub-graphs of the icosahedral graph obtained by
deleting a minimum of 6 and a maximum of 9 edges. There are only two structures
with maximum edge count of 24 corresponding to the sphere packings of the
\ac{fcc} and \ac{hcp} structures, which can be obtained from the icosahedral
structure by a smooth rearrangement moving the six spheres along a closed
zig-zag path into the (hexagonal) plane. The common (6,12)-\ac{LJ} potential has
only one minimum structure corresponding to the ideal icosahedron where the 12
outer spheres do not touch each other. Symmetry breaking requires a very
repulsive short-range LJ potential. The $(a,b)$-line in the $(a,b)$-\ac{LJ}
potential where the icosahedron completely disappears has also been determined.
While the results shown here depend on the functional form chosen (the
Lennard-Jones potential), similar results are expected for other well known
potentials such as the Morse potential.

It was also shown that for softer potentials, it is still unfavourable for a
13th outer sphere to touch the center sphere. The Gregory-Newton argument still
holds true for even the softest (1,2)-\ac{LJ} potential.

The sphere kissing problem in higher dimensions is a well known
problem\autocite{Conway_SpherePackingsLattices_1999} (in two dimensions there is
only 1 non-isomorphic \ac{GNC}). How many non-isomorphic rigid \acp{GNC} there
are in greater than three dimensions is currently unknown. Moreover, the rigid
kissing sphere problem can be extended to other (convex or not) topologies
instead of a central sphere, e.g. kissing spheres on an ellipsoid. Those would
all be very interesting questions to investigate.