text:
\subsubsection{Introduction to valence-bond solids}\index{valence-bond solid}

In section~\ref{sec:TOS} we have seen that the classically ordered
N\'eel state is a candidate to describe the groundstate of the
antiferromagnetic Heisenberg model Eq.~\eqref{eq:heisenberg}
with $J>0$ in the thermodynamic limit on a bipartite lattice. The energy
expectation value of this state on a single bond is $e_{\text{N\'eel}} = -J/4$.

The state which minimizes the energy of a single bond is, however, a singlet
state $\ket{S=0}$ formed by the two spins on the bond with energy
$e_{\text{VB}} = -3J/4$, called a valence bond (VB) or dimer. A valence bond
covering of an $N$-site lattice can then be described by a tensor product of
$N/2$ VBs, where each site belongs to exactly one VB\footnote{The set of all
possible valence bond coverings with arbitrary length spans the full
$S_{\text{tot}}=0$ sector of the models Hilbert space and is overcomplete
{\cite{Liang1988,Lhuillier2011}}.}.
Another possible candidate for the thermodynamic groundstate of
Eq.~\eqref{eq:heisenberg} is then a superposition of all possible VB coverings
with only nearest neighbour VBs. Such states do not break the $SU(2)$
spin-rotational symmetry as $S_{\text{tot}}=0$ and are in general not
eigenstates of the Hamiltonian: Acting with the operator $\mathbf{S}_i \cdot
\mathbf{S}_j$ between sites $i$ and $j$ belonging to two different VBs changes
the VB configuration.

This manifold of VB coverings is highly degenerate. As the VB coverings
are in general not eigenstates of the Hamiltonian, they encounter quantum
fluctuations. The energy corrections due to these fluctuations are typically
not equivalent for different coverings, although the bare energies are
identical.
The VB coverings with the largest energy gain are selected by the fluctuations
as the true groundstate configurations. If this \textit{order-by-disorder}
mechanism \cite{Shender1982,Henley1989} selects regular patterns \index{order by disorder}
of VB coverings, the discrete lattice symmetries can be spontaneously broken in the
thermodynamic limit, and a \textit{valence bond solid} (VBS)\index{valence bond
solid} can be
formed. Fig.~\ref{fig:cvbs} and Fig.~\ref{fig:svbs} show two different VBS
states on the square lattice. VBSs show no long-range spin order, but long-range
dimer-correlations $\langle(\mathbf{S}_a \cdot \mathbf{S}_{a'}) (\mathbf{S}_b
\cdot \mathbf{S}_{b'})\rangle$ where $a,a'$ and $b,b'$ label sites on
individual dimers. In section~\ref{sec:sym_analysis_VB} we will see how
different VBS states can be identified and distinguished by the quantum numbers
of the quasi-degenerate groundstate manifold on finite-size systems.

The groundstate of the Heisenberg model Eq.~\eqref{eq:heisenberg} on the square
lattice is not a VBS but a N\'eel state, which has a lower variational energy
already on the classical level. Nevertheless, several models featuring VBS
groundstates are known in 1- and 2-D
\cite{Fouet2001,Lauchli2002,Lauchli2005,Mambrini2006,Gelle2008}. Interestingly,
in \cite{Sandvik2007} a model was proposed, which shows a direct continuous
quantum phase transition between a N\'eel state and a VBS. This transition
exhibits very exotic, non-classical behaviour and is called \textit{deconfined
quantum critical point} \cite{Senthil2004}.

\subsubsection{Identification of VBSs from finite-size spectra}
\label{sec:sym_analysis_VB}

\paragraph{Columnar valence-bond solid}
\index{valence-bond solid!columnar}
A columnar VBS (cVBS) on a square lattice is shown in Fig.~\ref{fig:cvbs}. Four
equivalent states can be found, indicating that there will be a four-fold
quasi-degenerate groundstate manifold. A cVBS breaks the translational
and point-group symmetries of an isotropic SU(2)-invariant Hamiltonian on the
lattice spontaneously but not the continuous spin symmetry group since it
is a singlet and thus invariant under spin rotations.

\begin{figure}[h]
 \centering
 \includegraphics[width=\textwidth]{cvbs}
 \caption{The four columnar VBS coverings of a square lattice. Valence bonds
(spin singlets) are indicated by blue ellipses.}
 \label{fig:cvbs}
\end{figure}

In the following we use Eq.~\eqref{eq:stabilizerformula} to compute the
symmetry sectors of the groundstate manifold. The discrete symmetry group
we consider is
\begin{equation}
 \mathcal{G} = \mathcal{D} = \mathcal{T} \times \text{PG}
\end{equation}
where $\mathcal{T} = \mathbb{Z}_2 \times \mathbb{Z}_2 = \{1, t_x, t_y,
t_x t_y\}$ are the non-trivial lattice translations with translation vectors
\begin{equation}
 \mathbf{t}_1 = (0,0), \quad \mathbf{t}_x = (1,0), \quad
\mathbf{t}_y= (0,1), \quad \mathbf{t}_{xy} = (1,1)
\end{equation}
and $\text{PG} = \text{C}_4$ denotes the point-group of four-fold lattice
rotations\footnote{The dihedral group $\text{D}_4$ is also a symmetry group of
the model. For the sake of simplicity we decided to only consider the subgroup
$\text{C}_4$ in this section.}. To compute the groundstate symmetry sectors we
do not need to consider the full symmetry group $\mathcal{G}$ but only the
stabilizer $\text{Stab}(\ket{\Psi_{cVBS}})$, leaving one of the states in
Fig.~\ref{fig:cvbs} unchanged. Without loss of generality we choose the first
covering as prototype
$\ket{\Psi_{cVBS}}$. The stabilizer is given by
\begin{equation}
 \text{Stab}(\ket{\Psi_{cVBS}}) = \{1\times1\} \cup \{1\times C_2\} \cup
\{t_y \times 1\} \cup \{t_y \times C_2\}
\end{equation}
where $C_2$ denotes the rotation about an angle $\pi$ around the center of a
plaquette.

The irreducible representations (irreps) of the group of lattice translations
$\mathcal{T}$ can be labelled by the allowed momenta $\mathbf{k}$
\begin{equation}
 \mathbf{k} \in \mbox{Irreps}(\mathcal{T}) = \{(0,0), (\pi,0), (0,\pi),
(\pi,\pi)\},
\end{equation}
and the corresponding characters for an element $t \in \mathcal{T}$ are
\begin{equation}
 \chi_{\mathbf{k}}(t) = \mathrm{e}^{i \mathbf{k} \cdot \mathbf{t}}.
\end{equation}
The irreps (called A, B and E, see~\cite{Mulliken1955}) and
characters for the point-group
$\text{C}_4$ are tabulated in Tab.~\ref{tab:char_table_c4}.

\begin{table}[htb]
 \begin{center}
 \begin{tabular}{c|c|c|c|c}
  $\text{C}_4$ & $1$ & $C_4$ & $C_2$ & $(C_4)^3$ \\ \hline
  A & +1 & +1 & +1 & +1 \\
  B & +1 & -1 & +1 & -1 \\
  E$_a$ & +1 & +i & -1 & -i \\
  E$_b$ & +1 & -i & -1 & +i
 \end{tabular}
 \caption{Character table for pointgroup $\mathbf{C}_4$.}
 \label{tab:char_table_c4}
 \end{center}
\end{table}

Using the character-stabilizer formula Eq.~\eqref{eq:stabilizerformula}
we can now reduce the representation induced by the state $\ket{\Psi_{cVBS}}$
to irreducible representations to get the
quantum numbers of the quasi-degenerate groundstate manifold. Let us explicitely
consider $\mathbf{k}=(0,0)$ as an example:
\begin{align}
 n_{(0,0)A} &= n_{(0,0)B} = \frac{1}{|\text{Stab}(\ket{\Psi_{cVBS}})|} \sum_{d \in
 \text{Stab}(\ket{\Psi_{cVBS}})}
 \chi_{A}(d) \chi_{\mathbf{k}=(0,0)}(d) \\
 &= \frac{1}{4} \left[ 1\cdot  e^{i \mathbf{k} \cdot (0,0)} +
 1\cdot  e^{i \mathbf{k} \cdot (0,0)} +
 1\cdot  e^{i \mathbf{k} \cdot (0,1)} +
 1\cdot  e^{i \mathbf{k} \cdot (0,1)} \right] = 1 \\

 n_{(0,0)E_a} &= n_{(0,0)E_b} = \frac{1}{|\text{Stab}(\ket{\Psi_{cVBS}})|} \sum_{d \in
 \text{Stab}(\ket{\Psi_{cVBS}})}
 \chi_{B}(d) \chi_{\mathbf{k}=(0,0)}(d) \\
 &= \frac{1}{4} \left[ 1\cdot  e^{i \mathbf{k} \cdot (0,0)} +
 (-1)\cdot  e^{i \mathbf{k} \cdot (0,0)} +
 1\cdot  e^{i \mathbf{k} \cdot (0,1)} +
 (-1)\cdot  e^{i \mathbf{k} \cdot (0,1)} \right] = 0
\end{align}

Eventually, the cVBS covering will be described by a four-fold quasi-degenerate
groundstate manifold with the following quantum numbers\footnote{The little
group $\rho_{\mathbf{k}}$ for the momenta $\mathbf{k}=(0,\pi)$ and
$\mathbf{k}=(\pi,0)$ is only the subgroup C$_2$ of 2-fold roations of the
symmetry group C$_4$ considered in this example. The irrep called A therefore
denotes the trivial irreducible representation of C$_2$ for these momenta. When
computing the multiplicities for these momenta one should also note, that two
prototype states have to be considered to span the full groundstate manifold
under the symmetry elements of C$_2$.}.
\begin{equation}
 \chi(\ket{\Psi_{cVBS}}) = (0,0) \mbox{A} \oplus (0,0) \mbox{B} \oplus (\pi,0)
\mbox{A}
\oplus (0,\pi) \mbox{A}.
\end{equation}
VBS states are a superposition of spin singlets on the lattice, therefore the
spin quantum number for all levels in the groundstate manifold must be trivial,
$S_{\text{tot}}=0$.

\paragraph{Staggered valence-bond solid}
\index{valence-bond solid!staggered}
The columnar VBS is not the only regular dimer covering of the square lattice.
Another possible regular covering is the staggered VBS (sVBS), where again four
equivalent configurations span the groundstate manifold. One of these
configurations is shown in Fig.~\ref{fig:svbs}.

\begin{figure}[htb]
 \centering
 \includegraphics[width=.25\textwidth]{svbs}
 \caption{One of the four identical staggered VBS coverings on the square
lattice.}
 \label{fig:svbs}
\end{figure}

Similarly, also the sVBS spontaneously breaks the translational and point-group
symmetries of an isotropic Hamiltonian, but not the spin-rotational symmetry.
Following the same steps as before we can compute the quantum
numbers of the four quasi-degenerate groundstates for the sVBS. The stabilizer
turns out to be different to the case of the cVBS and thus also the
decomposition into irreps yields a different result:
\begin{equation}
 \chi(\ket{\Psi_{sVBS}}) = (0,0) \mbox{A} \oplus (0,0) \mbox{B} \oplus
 (\pi,\pi) \mbox{E}_a \oplus  (\pi,\pi) \mbox{E}_b.
\end{equation}

Tab.~\ref{tab:irreps_VBS} shows a comparison of the irreducible representations
in the groundstate manifold of the cVBS and sVBS states.

\begin{table}[ht]
\centering
\begin{tabular}{c|cc}
Irreps & cVBS & sVBS\\ \hline
$(0,0)$ A & 1 & 1\\
$(0,0)$ B & 1 & 1\\
$(\pi,0)$ A & 1 & 0\\
$(0,\pi)$ A & 1 & 0\\
$(\pi,\pi) \text{E}_a$  & 0 & 1\\
$(\pi,\pi) \text{E}_b$  & 0 & 1\\
\end{tabular}
\caption{Multiplicities of the irreducible representations in the four-fold
degenerate groundstate manifolds of the columnar and staggered VBS on a square
lattice.}
\label{tab:irreps_VBS}
\end{table}

By a careful analysis of the quasi-degenerate states and their quantum numbers
on finite systems it is thus possible to identify and distinguish different
VBS phases which spontaneously break the space group
symmetries in the thermodynamic limit.

\section{Introduction} \index{tower of states}
\input{introduction}

\section{Tower of states} \label{sec:TOS} \index{tower of states!scaling}
We start our discussion on spontaneous symmetry breaking by
investigating the Heisenberg model on the square
\index{Heisenberg model}
lattice. Its Hamiltonian is given by
\begin{equation}
  \label{eq:heisenberg}
  H = J\sum\limits_{\left< i, j \right>} \mathbf{S}_i \cdot \mathbf{S}_j
\end{equation}
and is invariant under global SU($2$) spin rotations, i.e. a rotation
of every spin on each site with the same rotational SU($2$) matrix.
Therefore the total spin
\begin{equation}
  \label{eq:totalspin}
  \mathbf{S}_{\text{tot}}^2 = \left(\sum\limits_{i} \mathbf{S}_i\right)^2=
S_{\text{tot}}(S_{\text{tot}}+1)
\end{equation}
is a conserved quantity of this model and every state in the
spectrum of this Hamiltonian can be labeled via its total spin quantum
number $S_{\text{tot}}$. The Heisenberg Hamiltonian on the square
lattice has the property of being \textit{bipartite}: The lattice
can be divided into two sublattices $A$ and $B$ such that every term in
Eq.~\eqref{eq:heisenberg} connects one site from sublattice $A$ to
sublattice $B$.
It was found out early \cite{Anderson1952} that the groundstate
of this model bears resemblance with the classical \textit{N\'{e}el
  state} \index{N\'{e}el antiferromagnet}
\begin{equation}
  \label{eq:neelclass}
  \ket{\text{N\'{e}el class.}} = \ket{\uparrow \downarrow \uparrow \downarrow
  \cdots}
\end{equation}
where the spin-ups live on the $A$ sublattice and the spin-downs live
on the $B$ sublattice. The total spin $S_{\text{tot}}$
is not a good quantum number for this state. From elementary spin
algebra we know that
it is rather a superposition of several states with different total
spin quantum numbers. For example the 2-site state
\begin{equation}
  \label{eq:singletplustriplet}
  \ket{\uparrow \downarrow } = \frac{\ket{\uparrow \downarrow} -
    \ket{\downarrow\uparrow}}{2} + \frac{\ket{\uparrow \downarrow} + \ket{\downarrow\uparrow}}{2}
     = \ket{S_{\text{tot}}=0, m=0} + \ket{S_{\text{tot}}=1, m=0}
\end{equation}
is the superposition of a singlet ($S_{tot}=0$) and a triplet ($S_{tot}=1$).
Therefore if such a state were to be a
groundstate of Eq.~\eqref{eq:heisenberg} several states with different
total spin would have to be degenerate. It turns out that on finite bipartite
lattices this is not the case: The total groundstate of the
Heisenberg model on bipartite lattices can be proven to be a singlet state with
$S_{\text{tot}}=0$.
This result is known as \textit{Marshall's Theorem}
\cite{Marshall1955,Lieb1962,Auerbach1994}.
So how can the N\'{e}el state resemble the singlet groundstate?
To understand this we drastically simplify the Heisenberg model
and investigate a toy model whose spectrum can be fully understood
analytically.

\subsection{Toy model: the Lieb-Mattis model} \index{Lieb-Mattis model}
By introducing the Fourier transformed spin operators
\begin{equation}
  \label{eq:spinoperatorfourier}
  \mathbf{S}_\mathbf{k} = \frac{1}{\sqrt{N}}\sum\limits_{j=0}^N
e^{i\mathbf{k} \cdot \mathbf{x}_j}\mathbf{S}_j
\end{equation}
we can rewrite the original Heisenberg Hamiltonian in terms of these
operators as
\begin{equation}
  \label{eq:heisenbergfourier}
  H = J\sum\limits_{\mathbf{k} \in \text{BZ}} \omega_\mathbf{k}
\mathbf{S_k}\cdot\mathbf{S}_{-\mathbf{k}}
\end{equation}
where $\omega_\mathbf{k} = \cos(k_x) + \cos(k_y)$ and the sum over $\mathbf{k}$
runs over the momenta within the first Brillouin zone~(BZ). Let $\mathbf{k}_0 =
(\pi,\pi)$ be the ordering wavevector which is the dual to the translations
that leave the square N\'{e}el state invariant. We now want to look at the
truncated Hamiltonian
\begin{equation}
  \label{eq:liebmattis}
  H_{\text{LM}} = 2J\left(\mathbf{S}_{(0,0)}^2 -
\mathbf{S}_{\mathbf{k}_0}\cdot\mathbf{S}_{-\mathbf{k}_0}\right)
\end{equation}
where we omit all Fourier components in
Eq.~\eqref{eq:heisenbergfourier} except $\mathbf{k}=(0,0)$ and
$\mathbf{k}_0 = (\pi,\pi)$.
This model is called the Lieb-Mattis model \cite{Lieb1962} and has a
simple analytical solution. To see this we notice that
Eq.~\eqref{eq:liebmattis} can be written as
\begin{equation}
  \label{eq:liebmattis2}
  H_{LM} = \frac{4J}{N}\sum\limits_{i\in A, j\in B} \mathbf{S}_i \cdot
\mathbf{S}_j
\end{equation}
in real space, where $A$ and $B$ denote the two bipartite sublattices of the
square lattice and each spin is only coupled with spins
in the other sublattice. The interaction strength is equal
regardless of the distance between the two spins. Thus this
model is not likely to be experimentally relevant. Yet it will serve
as an illustrative example how breaking the spin-rotational symmetry
manifests itself in the spectrum of a finite size system.
We can rewrite Eq. \eqref{eq:liebmattis2} as
\begin{align}
  \label{eq:liebmattis3}
  H_{LM} &=  \frac{4J}{N}\left(\sum\limits_{i,j \in A\cup B} \mathbf{S}_i \cdot
\mathbf{S}_j
    - \sum\limits_{i,j \in A} \mathbf{S}_i \cdot \mathbf{S}_j -
    \sum\limits_{i,j \in B} \mathbf{S}_i \cdot \mathbf{S}_j \right) \\
  &=  \frac{4J}{N}( \mathbf{S}_{\text{tot}}^2 - \mathbf{S}_A^2 - \mathbf{S}_B^2)
\end{align}
This shows that the Lieb-Mattis model can be considered as the
coupling of two large spins $S_A$ and $S_B$ to a total spin $S_{\text{tot}}$.

We find that the operators $\mathbf{S}_{\text{tot}}^2$, $S_{\text{tot}}^z$,
$\mathbf{S}_A^2$ and $\mathbf{S}_B^2$
commute with this Hamiltonian and therefore the sublattice spins
$S_A$ and $S_B$ as well as the total spin $S_{\text{tot}}$ and its z-component
$m_{\text{tot}}$ are good quantum numbers for this model.
For a lattice with $N$ sites ($N$ even) the sublattice spins can be chosen in
the range $S_{A,B}\in \{0,1,\dots,N/4\}$ and by coupling them
\begin{align}
 S_{\text{tot}} &\in \{|S_A-S_B|,|S_A-S_B|+1,\dots,S_A+S_B\} \\
 m_{\text{tot}} &\in \{-S_{\text{tot}}, -S_{\text{tot}}+1, \dots,
S_{\text{tot}}\}
\end{align}
can be chosen\footnote{This set of states spans the full Hilbertspace of the
model.}. A state $\ket{S_{\text{tot}}, m, S_A, S_B}$ is thus an eigenstate of
the systems with energy
\begin{equation}
 E(S_{\text{tot}}, m, S_A, S_B) = \frac{4J}{N}
 \left[S_{\text{tot}}(S_{\text{tot}}+1) - S_A(S_A+1) - S_B(S_B+1) \right]
 \label{eq:energy_lm}
\end{equation}
independent of $m$, so each state is at least ($2 S_{\text{tot}}+1$)-fold
degenerate.

\paragraph{Tower of states}
We first want to consider only the lowest energy states for each
$S_{\text{tot}}$ sector. These states build the famous \textit{tower of states}
and collapse in the thermodynamic limit to a highly degenerate groundstate
manifold, as we will see in the following.

For a given total spin $S_{\text{tot}}$ the lowest energy states are built by
maximizing the last two terms in Eq.~\eqref{eq:energy_lm} with $S_A=S_B=N/4$ and
\begin{equation}
 E_0(S_{\text{tot}}) = E(S_{\text{tot}},m, N/4,N/4) = \frac{4J}{N}
 S_{\text{tot}}(S_{\text{tot}}+1) - J (\frac{N}{4}+1)
 \label{eq:energy_lm1}
\end{equation}
The groundstate of a finite system will thus be the singlet state with
$S_{\text{tot}}=0$
\footnote{The groundstate of the Heisenberg model
Eq.~\eqref{eq:heisenberg} on a bipartite sublattice with equal sized sublattices
is also proven to be a singlet state $S_{tot}=0$ by Marshall's Theorem
\cite{Auerbach1994, Marshall1955, Lieb1962}.}
. On a finite system the
groundstate is, therefore, totally symmetric under global spin rotations and
does not break the $SU(2)$-symmetry. In the thermodynamic limit $N\rightarrow
\infty$, however, the energy of all these states scales to zero and all of them
constitute to the groundstate manifold.

The classical N\'eel state with fully polarized spins on each
sublattice can be built out of these states by a linear combination of all the
$S_{\text{tot}}$ levels with $m_{tot}=0$ \cite{Lhuillier2005}. All other N\'eel
states pointing in a different direction in spin-space can be equivalently
built out of this groundstate manifold by considering linear combinations with
other $m_{tot}$ quantum numbers. In the thermodynamic limit, any infinitesimal
small field will force the N\'eel state to choose a direction and the
groundstate spontaneously breaks $SU(2)$-symmetry.

The states which constitute the groundstate manifold in the thermodynamic limit
can be readily identified on finite-size systems as well, where their energy
and spin quantum number are given by Eq.~\eqref{eq:energy_lm1}. These states
are called the \textit{tower of states} (TOS) or also \textit{Anderson tower},
\textit{thin spectrum} and \textit{quasi-degenerate joint states}
\cite{Anderson1952, Kaplan1990,Hasenfratz1993,Azaria1993}.

\paragraph{Excitations}
The lowest excitations above the tower of states can be built by lowering the
spin of one sublattice $S_A$ or $S_B$ by one, see Eq.~\eqref{eq:energy_lm}. Let
us set $S_A=N/4$ and $S_B=N/4-1$ which implies that $S_{\text{tot}} \in
\{1,2,\dots,N/2-1\}$. We can directly compute the energy $E_1(S_{\text{tot}})$
of these excited states for each allowed $S_{\text{tot}}$. The energy gap to
the tower of states
\begin{equation}
 E_{\text{exc}}(S_{\text{tot}}) = E_1(S_{\text{tot}}) - E_0(S_{\text{tot}}) = J
\end{equation}
is constant\footnote{This is an artifact of the infinite-range interaction in
the Lieb-Mattis model. In the original Heisenberg model these modes become
gapless magnon excitations.}.
Hence, the lowest excitations of the Lieb-Mattis model
are static spinflips. The next lowest excitations are spinflips on
both sublattices, $S_A=S_B=N/4-1$ with excitation
energy $E_{\text{exc}_2} = 2J$ and $S_{\text{tot}} \in \{0,1,\dots,N/2-2\}$. We
observe that only the energy gap of the TOS levels vanishes
in the thermodynamic limit, so the TOS indeed solely contributes to the
groundstate manifold.

\paragraph{Quantum Fluctuations}
When we introduced the Lieb-Mattis model Eq.~\eqref{eq:liebmattis} from the
Heisenberg model Eq.~\eqref{eq:heisenbergfourier} we neglected all Fourier
components except of $\mathbf{k}=(0,0)$ and $\mathbf{k} = \mathbf{k}_0$. This
was a quite crude approximation and it is not guaranteed that all results for
the Lieb-Mattis model will survive for the short-range Heisenberg model. To get
some first results regarding this question, we can introduce small quantum
fluctuations on top of the N\'eel groundstate of the Lieb-Mattis model and
perform a perturbative spin-wave analysis in first order\footnote{A more
detailed discussion can be found in \cite{Lhuillier2005}.}. This approach does
not affect the scaling of the tower of states levels, but it has an important
effect on the excitations. They are not static particles anymore, but become
spinwaves (magnons) with a dispersion, which is linear around the ordering-wave
vector $\mathbf{k}=\mathbf{k}_0$ and $\mathbf{k}= (0,0)$. On a finite-size
lattice the momentum space is discrete with a distance proportional to $1/L$
between momentum space points, where $L$ is the linear size of the system. The
energy of the lowest excitation above the TOS, the single magnon gap, therefore
scales as $E_{\text{exc}} \propto J/L$ to zero\footnote{In the thermodynamic
    limit the single magnon mode is gapless and has linear dispersion around
    $\mathbf{k}=\mathbf{k}_0$ and $\mathbf{k}= (0,0)$. It corresponds to the
    well-known Goldstone mode which is generated when a continuous symmetry is
spontaneously broken.}.
As the scaling is, however, slower for $d>1$-dimensional systems than the TOS
scaling, these levels do not influence the groundstate manifold in the
thermodynamic limit. Furthermore, the excitation of two magnons results in a
two-particle continuum above the magnon mode.

The properties of the TOS and its excitations are summarized in
Fig.~\ref{fig:TOS_Neel}. The left figure shows the general properties of the
finite-size energy spectrum which can be expected when a continuous symmetry
group is spontaneously broken in the thermodynamic limit. The right figure
depicts the energy spectrum for the Heisenberg model on a square lattice with
$N=32$ sites, obtained with Exact Diagonalization. One can clearly identify the
TOS, the magnon dispersion and the many-particle continuum. The existence of a
N\'eel TOS was not only confirmed numerically for the Heisenberg model on the
square lattice, but also with analytical techniques beyond the simplification
to the Lieb-Mattis model \cite{Anderson1952,Hasenfratz1993,Azaria1993}. The
different symbols in Fig.~\ref{fig:TOS_Neel} represent different quantum numbers
related to the space-group symmetries on the lattice. In the next section we
will see that the structure of these quantum numbers depends on the exact shape
of the symmetry-broken state and we will learn how to compute them.

\begin{figure}[ht]
 \centering
   \includegraphics[width=.4\textwidth]{TOS_Scheme}\hspace{10mm}$\mbox{}$
  \includegraphics[width=.4\textwidth]{TOS_Neel_HB_square_N.32_final}
  \caption{(Left): Schematic finite-size energy spectrum of an antiferromagnet
   breaking SU($2$) spin-rotational symmetry. The TOS levels are the
lowest energy levels for each total spin $S$ and scale with
$1/N$ to the groundstate energy. The low energy magnon
excitations are seperated from the TOS and a continuum of higher
energy states and scale with $1/L$. (Right): Energy spectrum for the
Heisenberg model on a square lattice from ED. The TOS levels are connected by a dashed
line. The single magnon dispersion (green boxes) with $S_{\text{tot}} \in
\{1,2,\dots\}$ are well separated from the TOS and the higher multi-particle
continuum. The different symbols represent quantum numbers related to
space-group symmetries and agree with the expectations for a N\'eel state
(See section~\ref{sec:sym_analysis}).}
 \label{fig:TOS_Neel}
\end{figure}

\section{Symmetry analysis} \label{sec:sym_analysis}

\input{symmetry_analysis}

\section{Examples}\label{sec:examples}
\subsection{Discrete symmetry breaking}

\input{discrete_symm_breaking}

\subsection{Continuous symmetry breaking}\label{sec:cont_symm_breaking}
In this section we give several examples of systems breaking continuous
SO($3$) symmetry. We discuss the introductory example of the square
lattice Heisenberg antiferromagnet, calculate the irreps in the TOS and
compare this to actual energy spectra from Exact Diagonalization on a
finite lattice in section~\ref{sec:squareantif}.
In section~\ref{sec:magtri} we discuss three magnetic orders on the
triangular lattice and an extended Heisenberg model where all of these are
stabilized. We present results from Exact Diagonalization and compare the
representations in these spectra to the predictions from TOS analysis.
Finally, we introduce quadrupolar order and show that also this kind
of symmetry breaking can be analyzed using the TOS technique in
section~\ref{sec:quadru}.

\subsubsection{Heisenberg antiferromagnet on the square lattice}
\label{sec:squareantif}
\index{Heisenberg model}\index{antiferromagnetism}
\input{heisenberg_square_su2_breaking}

\subsubsection{Magnetic order on the triangular lattice}\label{sec:magtri}
\input{magnetic_order_triangular}

\subsubsection{Quadrupolar order}\label{sec:quadru}
\index{quadrupolar order}
\input{quadrupolar}

\section{Outlook}
\input{outlook}

\clearpage
\bibliographystyle{correl}

\subsection{Representation theory for space groups}
\index{representation theory}
For finite discrete groups such as the space group
of a finite lattice the full set of irreducible representations
(irreps) can be worked out.
Let us first discuss some basic groups.
Let's consider a $n\times n$ square lattice with periodic
boundary conditions and a translationally invariant Hamiltonian
like the Heisenberg model on it. In the following we will
set the lattice spacing to $a=1$. The discrete symmetry
group we consider is $\mathcal{T} = \mathbb{Z}_n \times \mathbb{Z}_n$
corresponding to the group of translations on this lattice.
This is an Abelian group of order $n^2$. Its representations
can be labeled by the momentum vectors $\mathbf{k} = (\frac{2\pi
  i}{n}, \frac{2\pi j}{n})$, $i,j \in \{ 0, \cdots, n-1\} $
which just correspond to the reciprocal Bloch vectors defined on
this lattice. Put differently, the vectors $\mathbf{k}$ are the reciprocal lattice points
of the lattice spanned by the simulation torus of our $n\times n$ square lattice.
The character $\chi_\mathbf{k}$ of the $\mathbf{k}$-representation is given by
\begin{equation}
  \label{eq:transchars}
  \chi_\mathbf{k}(\mathbf{t}) = \text{e}^{i\mathbf{k}\cdot \mathbf{t}}
\end{equation}
where $\mathbf{t} \in \mathcal{T}$ is the vector of translation. This is just
the usual Bloch factor for translationally invariant systems.

Let us now consider a (symmorphic) space group of the form
$\mathcal{D} = \mathcal{T}\times \text{PG}$ as
the discrete symmetry group of the lattice
where $\text{PG}$ is the pointgroup of the lattice. For a model on
a $n\times n $ square lattice this could for example be the dihedral
group of order 8, D$_4$, consisting of fourfold rotations
together with reflections. The representation theory and
the character tables of these point groups are well-established.
An example for such character tables can be found in
Tabs.~\ref{tab:char_table_c4} and \ref{tab:d6chartable} for the
cyclic group $C_4$ and the dihedral group $D_6$~\footnote{We follow the
labeling scheme for point group representations according to
Mulliken~\cite{Mulliken1955}.}.
Since $\mathcal{D}$ is now a product of the translation and
the point group we could
think that the irreducible representations of $\mathcal{D}$
are simply given by the product representations $(\mathbf{k} \otimes \rho)$
where $\mathbf{k}$ labels a momentum representation and $\rho$ an irrep
of $\text{PG}$. But here is a small yet important caveat. We have to be careful
since $\mathcal{D}$ is only a semidirect product of groups as
translations and pointgroup symmetries do not necessarily commute.
This alters the representation theory for this product of groups
and the irreps of $\mathcal{D}$ are not just simply the products
of irreps of $\mathcal{T}$ and $\text{PG}$. Instead the full set
of irreps for this group is given by $(\mathbf{k} \otimes \rho_\mathbf{k})$ where
$\rho_\mathbf{k}$ is an irrep of the so called \textit{little group}
$L_\mathbf{k}$ of $\mathbf{k}$ defined as
\index{representation theory:little group}
\begin{equation}
  \label{eq:deflittlegroup}
  L_\mathbf{k} = \left\{ g \in \text{PG}; g(\mathbf{k}) = \mathbf{k} \right\}
\end{equation}
which is just the stabilizer of $\mathbf{k}$ in $\text{PG}$. For example, all
pointgroup elements leave $\mathbf{k}=(0,0)$ invariant, thus the little group
of $\mathbf{k}=(0,0)$ is the full pointgroup PG. In general, this does not hold
for other momenta and only a subgroup of  $\text{PG}$ will be
the little group of $\mathbf{k}$. In Fig.~\ref{fig:3sublatt_triangular}
we show the $\mathbf{k}$-points of a $6\times6 $ triangular
lattice together with its little groups as an example. The $K$ point
in the Brillouin zone has a D$_3$ little group, the $M$ point a D$_2$
little group. Having discussed the represenation theory for
(symmorphic) space groups we state that the characters of these
representations are simply given by
\begin{equation}
  \label{eq:spacechars}
  \chi_{(\mathbf{k},\rho_\mathbf{k})}(\mathbf{t},p) = \text{e}^{i \mathbf{k}\cdot \mathbf{t}} \chi_{\rho_\mathbf{k}}(p)
\end{equation}
where $\mathbf{t} \in \mathcal{T}$, $p \in \text{PG}$ and $\chi_{\rho_\mathbf{k}}$ denotes
the character of the representation $\rho_\mathbf{k}$ of the little group
$L_\mathbf{k}$.

\subsection[Predicting irreducible representations in spontaneous
  symmetry breaking]{Predicting irreducible representations in spontaneous
  symmetry breaking \subsectionmark{Predicting irreducible representations}
}
\subsectionmark{Predicting irreducible representations}
\index{tower of states!quantum numbers}
\index{representation theory:irreducible representations}
Spontaneous symmetry breaking at $T=0$ occurs when
the groundstate $\ket{\psi_{\textnormal{GS}}}$ of $H$ in the
thermodynamic limit is not invariant under the full symmetry group
$\mathcal{G}$ of $H$. We will call a specific groundstate
$\ket{\psi_{\textnormal{GS}}}$ a \textit{prototypical state} and
the \textit{groundstate manifold} is defined by
\begin{equation}
  \label{eq:defdegenerategs}
  V_{\textnormal{GS}} = \spn\left\{ \ket{\psi_{\textnormal{GS}}^i }\right\}
\end{equation}
where $\ket{\psi_{\textnormal{GS}}^i }$ is the set of degenerate
groundstates in the thermodynamic limit.
This groundstate manifold space can be finite or infinite dimensional
depending on the situation. For breaking a discrete finite symmetry,
such as in the example given in section \ref{sec:sym_analysis_VB},
this groundstate manifold will be finite dimensional, for breaking continuous SO($3$)
spin rotational symmetry\footnote{The actual symmetry group of
  Heisenberg antiferromagnets is usually SU($2$). For simplicity
  we only consider the subgroup SO($3$) in these notes which yields
  the same predictions for the case of sublattices with even number of
  sites (corresponding to integer total sublattice spin).}
as in section~\ref{sec:cont_symm_breaking}
it is infinite dimensional in the thermodynamic
limit. For every symmetry $g \in \mathcal{G}$
we denote by $O_g$ the symmetry operator acting on the Hilbert space.
The groundstate manifold becomes degenerate in the thermodynamic
limit and we want to calculate the quantum numbers of the groundstates
in this manifold. Another way of saying this is that we want to compute
the irreducible representations of $\mathcal{G}$ to which the groundstates belong
to.
For this we look at the action $\Gamma$ of the symmetry group $\mathcal{G}$ on
$V_{\textnormal{GS}}$ defined by
\begin{align}
  \label{eq:representationongroundstates}
  \Gamma: &\mathcal{G} \rightarrow \aut(V_{\textnormal{GS}})\\
  & g \mapsto \left(  \braket{\psi_{\textnormal{GS}}^i| O_g |\psi_{\textnormal{GS}}^j}\right)_{i,j}
\end{align}
This is a representation of $\mathcal{G}$ on $V_{\textnormal{GS}}$, so
every group element $g \in \mathcal{G}$ is mapped to an invertible
matrix on $V_{\textnormal{GS}}$. In general this representation is
reducible and can be decomposed into a direct sum of irreducible representations
\begin{equation}
  \label{eq:ssbirrepdecomp}
  \Gamma = \bigoplus_\rho n_\rho \rho
\end{equation}
These irreducible representations $\rho$ are the quantum
numbers of the eigenstates in the groundstate manifold and $n_\rho$
are its respective multiplicities (or degeneracies).
Therefore these irreps constitute the TOS for spontaneous symmetry
breaking \cite{Lhuillier2005}. To compute the
multiplicities we can use a central result from representation theory,
the \textit{character formula}
\begin{equation}
  \label{eq:characterformula}
  n_\rho = \frac{1}{|\mathcal{G}|}\sum\limits_{g\in \mathcal{G}}\overline{\chi_\rho(g)}\tr(\Gamma(g))
\end{equation}
where $\chi_\rho(g)$ is the character of the representation $\rho$ and
$\tr(\Gamma(g))$ denotes the trace over the representation matrix
$\Gamma(g)$ as defined in Eq.~\eqref{eq:representationongroundstates}.
Often we have the case that
\begin{equation}
  \label{eq:vanishingoverlapthermolim}
   \braket{\psi_{\textnormal{GS}} | O_g | \psi_{\textnormal{GS}}'} =
   \begin{cases}
     1 \text{ if } O_g\ket{\psi_{\textnormal{GS}}'} = \ket{\psi_{\textnormal{GS}}} \\
     0 \text{ otherwise}
   \end{cases}
\end{equation}
With this we can simplify Eq.~\eqref{eq:characterformula} to
what we call the \textit{character-stabilizer formula}
\begin{equation}
  \label{eq:stabilizerformula}
  n_{\rho} = \frac{1}{|\text{Stab}(\ket{\psi_{\textnormal{GS}}})|}\sum
  \limits_{g \in \text{Stab}(\ket{\psi_{\textnormal{GS}}})}\chi_\rho(g)
\end{equation}
where
\begin{equation}
  \label{eq:stabilizerdefinition}
  \text{Stab}(\ket{\psi_{\textnormal{GS}}}) \equiv \{ g \in \mathcal{G} :
  \; O_g
  \ket{\psi_{\textnormal{GS}}} = \ket{\psi_{\textnormal{GS}}}\}
\end{equation}
is the stabilizer of a prototypical state
$\ket{\psi_{\textnormal{GS}}}$
~\footnote{In some cases, the orbit of the prototypical state
{$G.\ket{\psi_{GS}} = \{g \in G: \; O_g \ket{\psi_{GS}}\}$} does
not span the full set of degenerate groundstates
{$\ket{\psi^i_{GS}}$}.
In this case, we have to find a set of prototypical states with
different orbits, such that the union of these orbits spans the
full groundstate manifold. Then, Eq.~\unexpanded{\eqref{eq:stabilizerformula}}
has to be applied to each prototypical state, individually, and
the final multiplicity is the sum of the individual results.
}.
We see that for applying the character-stabilizer formula in
Eq.~\eqref{eq:stabilizerformula} only two ingredients are needed:
\begin{itemize}
\item the stabilizer
  $\text{Stab}(\ket{\psi_{\textnormal{GS}}})$ of a prototypical state
  $\ket{\psi_{\textnormal{GS}}}$ in the groundstate manifold
\item the characters of the irreducible representations of the
  symmetry group $\mathcal{G}$
\end{itemize}
We want to remark that in the case
of $\mathcal{G} =  \mathcal{D} \times \mathcal{C}$
where $\mathcal{D}$ is a discrete symmetry group such as the spacegroup of a
lattice and $\mathcal{C}$ is a continuous symmetry group such as
SO($3$) rotations for Heisenberg spins the
Eqs.~\eqref{eq:characterformula} and \eqref{eq:stabilizerformula}
include integrals over Lie groups additionally to the sum over
the elements of the discrete symmetry group $\mathcal{D}$.
Furthermore also the characters for Lie groups like SO($3$)
are well-known. For an element $R \in $ SO($3$) the irreducible
representations are labeled by the spin $S$ and its characters
are given by
\begin{equation}
       \chi_S(R) = \frac{\sin\left[(S+\frac{1}{2})\phi\right]}{\sin\frac{\phi}{2}}
\end{equation}
where $\phi \in [0, 2\pi]$ is the angle of rotation of the spin
rotation $R$. We work out several examples for this case in
section~\ref{sec:cont_symm_breaking}
and compare the results to actual numerical data from Exact Diagonalization.


text:
\section{Introduction}

Field-effect transistors (FETs) with thin films of organic $\pi$-conjugated materials as active semiconductor have developed towards a serious alternative for low-cost and flexible electronics on large areas.\cite{sekitani2,sekitani} The fact that charge transport in those organic films occurs essentially in the first molecular monolayer in close proximity to the dielectric interface \cite{tanase,mottaghi} makes a self-assembled monolayer field-effect transistor (SAMFET) a perfect device for studying charge-carrier transport properties in organic materials.\cite{brondijk} In SAMFETs, the active $\pi$-system is linked to a flexible insulating alkyl-chain that is covalently bound to a suitable surface. Thus, the $\pi$-systems form an almost 2D confined transport channel on molecular scale in thickness \cite{smits,schmaltz} and the n-alkyl linkers compose a hybrid gate dielectric together with the insulating anchor oxide ($AlO_x$).\cite{Klauk} Hole and electron transport have been demonstrated in SAMFETs making the concept a veritable approach in molecular scale electronics.\cite{smits,schmaltz,novak2,ringk}

The experimental advances have triggered theoretical studies of charge transport
in organic thin films within the context of SAMFETs.\cite{brondijk,jaeger} The theoretical description of
charge transport in these systems involves the principle problem
of accounting for structural flexibility and inhomogeneity.
Thus, the characterization of the electronic structure must be combined with conformational
sampling, which in view of the complexity of the systems, represents a significant challenge.
Furthermore, the mechanism of electron transport is {\em a priori} not
obvious. While dephasing due to both structural fluctuations and inelastic processes
may favor hopping-type transport at higher temperatures, band-type coherent transport
can also contribute, in particular in systems with stronger intermolecular coupling and
at lower temperatures.

In previous work, some of the authors investigated charge transport in
$C_{60}$-based SAMs using a combination of molecular-dynamics (MD) simulation, semiempirical electronic
structure calculations and Monte Carlo (MC) transport studies.\cite{jaeger}
In agreement with experimental results,\cite{stubhan,novak} the study revealed the crucial influence of the
morphology of the SAM on the transport properties.
In the present paper, we extend the previous study using a fully quantum mechanical description of transport.
Specifically, we employ Landauer transport theory,\cite{Landauer57} which provides a
coherent transport treatment based on elastic scattering mechanisms and is thus complementary
to the MC-based investigations in Ref.\ \onlinecite{jaeger}.
Moreover, we use the theory of local currents \cite{pecchia} to investigate pathways of
electron transport in the SAM and study the influence of structural fluctuations on
the transport properties.

The article is organized as follows: In Sec.\
\ref{sec:methods_setup}, we introduce  the SAMFET
devices containing the layers of organic molecules,
 which are investigated throughout the paper. The
theoretical methodology  is
outlined in Sec.\ \ref{sec:methods_md} - \ref{sec:methods_loccurr},
including MD simulations, electronic
structure calculations, the Landauer transport approach and the theory of local currents.
In Sec.\ \ref{sec:results}, the results of the simulations obtained
for various SAMs, which differ in their structural properties, are presented. In particular, we
discuss the electronic structure and the transport properties
and their relation to the morphology of the SAMs.
In addition, local pathways for
charge transport in the SAMs are investigated and the influence of structural fluctuations on the
transport properties of the SAMs is analyzed.

\begin{figure*}
	\centering
		\includegraphics[width=.95\textwidth]{Figure1.eps}
		\caption{Schematic and chemical composition
 of SAMFET devices. (a) Chemical structure of the SAM
forming phosphonic acids $C_{10}$-PA ({\bf 1})
and $C_{60}$-$C_{18}$-PA ({\bf 2}).
(b) Schematic SAMFET setup with $Au$ source and
drain electrodes, and $AlO_x$ gate dielectric between SAM and an $Al$ gate.}
\label{fig:device}
\end{figure*}

\section{Systems and Methods}  \label{sec:methods}

The theoretical methodology employed to simulate charge transport in
SAMs involves three steps: (i) a characterization of the structure of
the SAMs with MD simulations, (ii) a determination of the electronic structure
using semiempirical molecular orbital (MO) theory, and (iii) transport calculations
within Landauer theory. In the following, the different methods are described.
We begin with a brief introduction of the systems investigated and their use in SAMFET devices.

\subsection{Setup} \label{sec:methods_setup}

A scheme of the self-assembled monolayer field-effect
 transistor (SAMFET) device and its chemical
composition is shown in Fig.\ \ref{fig:device}. The SAM is formed by fullerene-functionalized
octadecyl-phosphonic acids (PAs) (in the following denoted by
$C_{60}$-$C_{18}$-PA)
 and $C_{10}$-PA in various stoichiometric ratios.
The SAM is separated from the aluminum gate electrode
 at the bottom by an $AlO_x$ layer, which together
with the alkyl chains of the SAM builds the dielectric
 of the device.\cite{tbauer} Lithographically
patterned gold (30 nm) is placed on top of the SAM, serving as
 source and drain electrodes.\cite{jaeger} The
semiconducting $C_{60}$ head groups of the functionalized
 PA in the SAM form the active transistor channel in the device.
In this paper, we focus on the modeling of charge transport within
the SAM. Therefore, the influence of a gate potential and
the $AlO_x$ layer is not taken into account. The coupling to the
gold electrodes is described implicitely using self energies (see below).

\subsection{Molecular-Dynamics Simulations and Semiempirical MO Theory} \label{sec:methods_md}

The underlying conformational sampling of the self-assembled monolayer systems investigated in this publication is based on classical atomistic
molecular-dynamics simulations published and described in detail
previously.\cite{jaeger} Here, the SAM-forming molecules were initially aligned
on a clean (0001) aluminum-oxide surface that was equilibrated prior to
depositing phosphonic acids using an interatomic potential model
parameterized by Sun et al..\cite{sunj} The parameters for the phosphonates are
based on the general Amber force field (GAFF) \cite{wangj} and the MD simulations
were performed with the program DL-POLY.\cite{todorovsmith} Details of the simulation
are given in Ref.\ \onlinecite{jaeger}.

Following the MD simulations, an ensemble of snapshots was extracted and
the $AlO_x$ substrate was removed for further processing. The snapshot
geometries were then treated by semiempirical MO
calculations using the restricted Hartree-Fock formalism and the AM1
Hamiltonian.\cite{dewar} All semiempirical MO calculations were performed using
the parallel EMPIRE program. \cite{empire}

\subsection{Model and Transport Theory} \label{sec:methods_transport}

To study charge transport through the SAM, we use an effective
single-particle model with parameters determined by the
semiempirical MO calculations discussed above. Correspondingly, for a given
snapshot of the structure of the SAM, the
Hamiltonian of the SAM reads
\begin{equation}
H_S= \sum_j |\phi_j\rangle\epsilon_j \langle \phi_j|.
\end{equation}
Thereby, $\epsilon_j$ denotes the energy of an electron
in the $j$th MO  $|\phi_j\rangle$.

Within the Landauer approach, transport through the SAM is described by
the transmission function  $t(E)$, which for an electron with energy
$E$ is given by
\cite{Datta}
\begin{equation}
 t(E)={\rm tr}\left\lbrace\Gamma_RG_S\Gamma_LG_S^{\dagger}\right\rbrace.
\label{eqn:1}
\end{equation}
Hereby, $G_S$ denotes the retarded Green's function of the SAM,
\begin{equation}
 G_S(E)=\left[ E^+ - H_S-\Sigma_L(E)-\Sigma_L(E)\right]^{-1},
\end{equation}
which involves the self energies $\Sigma_{\alpha}(E)$
that describe the coupling of the SAM to the left and
right electrodes ($\alpha=L,R$).
 The self energies are related to the width-function
$\Gamma_{\alpha}(E)$, used in Eq.\ (\ref{eqn:1}), via
\begin{equation}
 \Sigma_{\alpha}(E)= \Delta_{\alpha}(E)-\frac{i}{2}\Gamma_{\alpha}(E),
\end{equation}
where $\Delta_{\alpha}$ denotes the level-shift function.

The calculation of the self energy requires explicit modeling of the gold
electrodes and the gold-SAM interface.\cite{Benesch08,cuevasscheer} In this paper,  we
use a simpler strategy. We assume that the gold electrodes
couple preferentially to selected
carbon atoms of the $C_{60}$ head groups at the
left and right boundaries of the SAM. This choice is suggested by the
geometrical structure of the SAM. Furthermore, we invoke the
wide-band limit, where $\Gamma_{\alpha}$ can be approximated
by a constant value and
$\Delta_{\alpha}$ vanishes. The wide-band limit is a good
approximation for gold electrodes.
Specifically, the matrix elements of the self energies $\Sigma_{\alpha}(E)$
in a local basis (represented by atomic orbitals $|\chi_{\nu}\rangle$)
are given by
\begin{equation}
 (\Sigma_{\alpha}(E))_{\nu\nu}=-\frac{i}{2}(\Gamma_{\alpha})_{\nu\nu}
\end{equation}
 with $(\Gamma_{\alpha})_{\nu\nu}= 1 eV$
for orbitals $\nu$ corresponding to the outermost hexagon of Carbon atoms
of the $C_{60}$ head groups at the left and right boundary of the SAM
and $(\Gamma_{\alpha})_{\nu\nu'}= 0$ otherwise.
The value $(\Gamma_{\alpha})_{\nu\nu}= 1$ eV respresents a reasonable choice
for molecule-gold contacts. It is emphasized
that the results
do not depend significantly on the exact value of this parameter, because
the bottleneck for charge transport in the systems considered is the SAM
itself and not the SAM-gold contact.

Based on the transmission function $t(E)$, the electrical current through
the SAM is given by the expression
\begin{equation}
 I(V) = \frac{2e}{h}\int {\rm d}E \, t(E)(f_L(E)-f_R(E)),
\label{eqn:curr}
\end{equation}
where $f_{L/R}(E)$ denotes the Fermi function for the electrons in the
left/right lead. It is given by
\begin{equation}
f_{\alpha}(E)= \frac{1}{1+e^{(E-\mu_{\alpha})/k_BT}}
\end{equation}
with the Boltzmann constant $k_B$, the electrode temperature $T$ and their chemical potentials $\mu_{\alpha}$ ($\alpha=L,R$). For a symmetric voltage drop around the Fermi energy $E_F$, the chemical potentials are given by
\begin{equation}
 \mu_{\alpha}=E_F\pm \frac{eV}{2},
\end{equation}
where $V$ denotes the bias voltage. In the evaluation of the current, we assume that the transmission function
$t(E)$ depends only weakly on the bias voltage and use its value at
zero bias voltage.

It is noted that Landauer transport theory, as used here,
can only describe  coherent transport based on elastic
scattering mechanisms. Dephasing and inelastic processes
due to electron-phonon or
electron-electron interaction are neglected. The study of these
processes requires an extension of the model
and more advances transport methods, which will be the topic of future work.
The influence of structural fluctuation will be considered in Sec.\
\ref{sec:results_fluc}.

\subsection{Analysis of Transport Pathways Using Local Currents}  \label{sec:methods_loccurr}

Eqs.\ (\ref{eqn:1}) and (\ref{eqn:curr}) describe charge transport in terms
of the overall electrical current and transmission function of the complete SAM.
In order to have a more detailed description and, in particular,
to analyze pathways
of charge transport in the SAM, we use  the method
 of local currents.\cite{Todorov,solomon}
As discussed in detail in Refs.\ \onlinecite{Todorov} and \onlinecite{solomon},
within this technique
the overall current $I$ through a surface perpendicular to the
transport direction
is represented in terms of local currents via
\begin{equation}
 I=\sum_{m\in M_L}\sum_{n\in M_R} I_{mn},\label{eqn:Iloc}
\end{equation}
where $M_{L/R}$ denotes atomic sites left and right of the chosen
surface, respectively.
As can be shown,\cite{solomon,pecchia} the local contributions
to the current from site $n$ to $m$ are given by
\begin{equation}
 I_{mn}=\frac{2e}{h}\int \frac{dE}{2\pi}K_{mn}(E),
\label{eqn:localcurrs}
\end{equation}
where

\begin{equation}
 K_{mn}(E)=\sum_{\nu\in m}\sum_{\mu\in n}
(V_{\nu\mu}G_{\mu\nu}^<(E)-V_{\mu\nu}G_{\nu\mu}^<(E)),
\label{eqn:2}
\end{equation}
and
\begin{equation}
 G^{<}=(if_LG_S\Gamma_L G_S^{\dagger}
+if_RG_S\Gamma_R G_S^{\dagger}).
\end{equation}
Thereby, we have used a representation in the local atomic orbital basis
$|\chi_{\nu}\rangle$, where the Hamiltonian reads
\begin{equation}
 H_S= \sum_{\mu} |\chi_{\mu}\rangle\tilde\epsilon_{\mu} \langle \chi_{\mu}|+\sum_{\mu\ne\nu} |\chi_{\nu}\rangle V_{\nu\mu} \langle \chi_{\mu}|.
\end{equation}
Here, $\tilde\epsilon_{\mu}$ denotes the energy of an electron in orbital
$\mu$ and $V_{\nu\mu}$ the coupling constant between orbitals $\nu$ and $\mu$.
For temperature $T=0$, the expression
$\sum_{m\in M_L}\sum_{n\in M_R} K_{mn}$ can be
identified with the transmission function, Eq.\ (\ref{eqn:1}),
and thus $K_{mn}$ defines local contributions to the transmission
between pairs of atoms.

\begin{figure*}
	\centering
		\includegraphics[width=.95\textwidth]{systems.eps}
		\caption{Top- and sideview of SAM consisting
of 25 $C_{60}$-$C_{18}$-PA (a) mixed with 75 $C_{10}$-PA
after a MD simulation time of 50 ns (b) without the $C_{10}$-PA,
 but with identical geometry as shown in (a),
(c) Pure  25 $C_{60}$-$C_{18}$-PA SAM morphology as
derived from a different simulation.}
\label{fig:systems}
\end{figure*}

In the applications considered below, we use a further coarse graining, where
the indices $m$ and $n$ in Eqs.\ (\ref{eqn:Iloc} - \ref{eqn:2}) refer not
to single atoms but to the $C_{60}$-$C_{18}$-
or $C_{10}$-PA moieties in the SAM. Consequently, the
sums in Eq.\ (\ref{eqn:2}) also extend
over the orbitals of all atoms of each $C_{60}$-$C_{18}$-
or $C_{10}$-PA moiety. This allows a more straightforward analysis
of the pathways of charge transport through the SAM.

\section{Results and Discussion}\label{sec:results}

We have investigated three different SAM configurations
as shown in Fig.\ \ref{fig:systems}, which differ both in the mixing ratio of
$C_{60}$-$C_{18}$ and $C_{10}$-PA groups and in their morphology.
The basic unit representing
the SAM is a cell consisting of 25 $C_{60}$-$C_{18}$ PAs.
Depending on the concentration of $C_{10}$-PA, the MD simulations
reveal pronounced differences in the SAM morphologies.
Fig. \ref{fig:systems} (a) and (c) show snapshots of a 25:75
mixed SAM (containing, in addition, 75 $C_{10}$-PAs) and a pure SAM (without
$C_{10}$-PAs), respectively, after a MD simulation time of 50 ns.
The pure SAM exhibits a considerably reduced thickness and an
enhanced disorder of the $C_{60}$ head groups. Moreover, the $C_{60}$
head groups are much closer to the surface in the pure SAM and,
therefore, favor leakage currents into the substrate. Monte Carlo (MC)
simulations and experiments have proven this morphology
to be less conductive and have suggested a mixing
ratio in the range of 25:75 -- 50:50  to optimize the performance of
the SAMFET devices. \cite{rumpel}

In addition, we have considered the SAM depicted in
 Fig. \ref{fig:systems} (b), which was obtained from the
mixed SAM in panel (a) by removing the insulating $C_{10}$ alkyl chains
after the MD simulation. Although the geometrical structure of this
system is artificial, it can give insight into
the role the $C_{10}$-PA plays in the transport through
the layer. Note that the realistic morphology of the
pure  SAM as derived from  MD simulations is the one shown in (c).

\subsection{Transport Properties and Molecular Orbital Analysis} \label{sec:results_mo}

We first consider the 25:75 SAM depicted in Fig.\ \ref{fig:systems} (a).
The energy spectrum of the MOs of this system in the range of [-10,-2] eV
is shown in Fig \ref{fig:spec}.
\begin{figure}
 \includegraphics[width=1\linewidth]{spectrum.eps}
\caption{ Electronic structure of the 25:75 mixed SAM.
(a) Energy levels shown as black lines. The blue dashed line denotes the Fermi
level of gold used in the transport calculations ($E_F=~-5.1$ eV).
(b) highest occupied
MO (c) lowest unoccupied MO. The MOs are plotted at isovalue $c=0.001$.}
\label{fig:spec}
\end{figure}
\begin{figure*}
\includegraphics[width=.9\textwidth]{figure3a.eps}
\caption{Electronic structure and transport properties of the 25:75 mixed
SAM Top: MO analysis, Center: Energy levels (black lines)
and transmission function
(red lines), Bottom: Current-voltage characteristic}
\label{fig:moanalysis}
\end{figure*}
The spectrum exhibits a dense distribution of
energy levels and a HOMO-LUMO energy
gap of 3.26 eV. In addition, the HOMO and the LUMO of the
SAM are depicted. Due to the pronounced structural disorder,
all MOs are localized. The analysis
reveals that the occupied orbitals close to the HOMO
are mainly localized on the alkyl chains, while the unoccupied
orbitals close to the LUMO are localized on the fullerenes.
As fullerenes are electron conductors, only unoccupied energy states are
relevant for transport processes. Therefore, in the
simulations it is sufficient to take into account only
the unoccupied part of the spectrum.

 Fig.\ \ref{fig:moanalysis} shows the transmission function
 of the 25:75 SAM in an energy range close to the Fermi level
(starting at $E_{\rm LUMO}$). In addition, the energy levels are
depicted as black lines. The transmission
 function exhibits a series of narrow peaks, which can be associated to
 individual energy levels of the molecular layer.
The small widths of the peaks are
 due to the strong localization of the
corresponding molecular orbitals, which results in a weak coupling
to the electrodes. The
current-voltage characteristic of the SAM at $T=300$ K,
depicted in Fig. \ref{fig:moanalysis} (bottom), shows a
non-Ohmic behavior. Thermally broadened steps in the current
appear at voltages,
 where the chemical potential of the left electrode,
$\mu_L=E_F+eV/2$, is in resonance with
energies $E_i$ of the transmission peaks, i.e. at
$V_i=2|E_F-E_i|$.

As in the transmission function, the step structures
in the current can be associated to contributions from
individual MOs. For example, the rise of the current at $V\approx 1.65$ V
is due to contributions of LUMO+1 and LUMO+2. Although their
transmission values are
comparably large, their contributions to the current
are small because of the small widths of their transmission
 peaks. The step at $V\approx 1.82$ V
corresponds to LUMO+5, which
exhibits a small transmission coefficient of $t=0.05$
but a much larger peak width.
At larger voltages, we
find the contributions of several higher unoccupied
MOs. The overall magnitude of the current is on the order of $nA$,
which is in agreement with experimental data.\cite{jaeger}

\subsection{Influence of Morphology and Mixing Ratio} \label{sec:results_mixed}

\begin{figure}
\includegraphics[width=.55\textwidth]{figure3b.eps}
\caption{Electronic structure and transport properties of the pure SAM with
  nonrelaxed geometry (see Fig.\ \ref{fig:systems} (b)) Top:
MO analysis, Center: Energy levels (black lines) and transmission function (red lines), Bottom: Corresponding current}
\label{fig:moanalysis_wo}
\end{figure}

In order to explore the influence of the morphology of the SAMs
and the mixing ratio of the SAM-forming molecules on the
transport properties, we consider the
pure $C_{60}$-$C_{18}$-PA SAMs with the two different
geometries shown in Fig.\ \ref{fig:systems} (b) and (c).
We first focus on the pure SAM with the unrelaxed geometry,
which was obtained by removing the $C_{10}$-PAs groups from
the 25:75 mixed SAM without further MD simulation.
The transmission function and the current-voltage characteristic
depicted in  Fig.\ \ref{fig:moanalysis_wo}
reveal the pronounced influence of the
$C_{10}$-PAs on the conductance properties of the SAMs.
Compared to the mixed SAM, we find that the LUMO
of the pure systems is shifted by almost 1 eV to
higher energies.  This and also a slight shift of the
HOMO to lower values leads to an increase of the bandgap compared to
the mixed SAM.
Due to the  larger bandgap, the onset of the current is shifted to
 higher voltages. However, the current magnitude is
significantly larger than in the mixed SAM.
The pronounced increase of the current at higher
voltages ($V \approx 3.8$ V) can be related to transport through
 LUMO+30. As illustrated in  Fig.\ \ref{fig:moanalysis_wo} (top panel),
LUMO+30 is a delocalized MO, which strongly
couples to the electrodes resulting in a high current
value of about 150 nA. Such delocalized MOs are
 not found in the mixed SAM. Thus, the presence of
the $C_{10}$-PAs in the mixed SAM
appears to decrease the coupling
between the fullerene head groups in the layer, which, in combination
with the disorder of the layer, results in the more localized orbitals
in the mixed SAM. Thus, in principle, the pure SAM with the geometry
considered in Fig.\ \ref{fig:moanalysis_wo} has (at least for
higher voltages) better conductance properties than the mixed SAM.
It should be noted, though, that this is only the case for the
artificial geometry considered in Fig.\ \ref{fig:moanalysis_wo}.
If the system is relaxed by further MD simulation, the morphology
changes drastically to the geometry depicted in Fig.\ \ref{fig:systems} (c),
which will be considered next.

\begin{figure}
\includegraphics[width=.55\textwidth]{figure3c.eps}
\caption{Analysis of pure SAM (equilibrated) Top:
MO analysis, Center: Energy levels (black lines) and transmission function (red lines), Bottom: Corresponding current}
\label{fig:moanalysis_flat}
\end{figure}

The electronic structure and conductance properties of
the pure SAM with the morphology depicted in
Fig.\ \ref{fig:systems} (c) are shown in
Fig.\ \ref{fig:moanalysis_flat}.
The results reveal that the almost flat and strongly disordered morphology
modifies the electronic structure profoundly and results
in a very small current, which confirms conclusions from
previous simulations and experiments.\cite{jaeger}
Evidently, due to the disorder, the
$C_{60}$ head groups in the layer cannot form
a conductive channel. The contributing MOs are
strongly localized, leading to narrow peak structures
 in the transmission function.
Also, we find a shift of the LUMO to higher energies.
The current onset is located at high voltages and
has values about one order of magnitude smaller
than for the mixed SAM. These findings show
that the poor conduction properties of pure
$C_{60}$-$C_{18}$-PA SAMs can clearly be related
 to the decreased order in the molecular
structure. Additionally, when used in
SAMFETs, higher leakage currents can appear
 due to the smaller distance of the fullerenes
to the substrate, as was indeed found in the transport measurements
\cite{burkhardt,jedaa}.

\subsection{Pathways for Electron Transport} \label{sec:results_paths}

So far, we have analyzed the conductance properties of the SAMs
in terms of the the overall transmission function and currents.
To obtain a more detailed understanding of the transport properties
and mechanisms,
we analyze the local pathways for electron transport through the SAM.
Previous investigations using Monte
Carlo path searches based on a local electron affinity potential
in the mixed SAM have suggested
 that there are specific pathways for efficient
electron transport across the fullerenes. \cite{jaeger}
In the following we will extend this analysis based on the
quantum mechanical transport theory employed here using the
theory of local currents and transmissions outlined in Sec.\
\ref{sec:methods_loccurr}.

\begin{figure}[h!]
	\centering
\includegraphics[width=0.75\linewidth]{localtrans.eps}
\caption{Local transmissions between pairs of
$C_{60}$-$C_{18}$-PAs and their relation to the
MO structure. The colors of the arrows are ordered from blue
to orange referring to an increasing transmission amplitude from
0 to 0.8. Black arrows indicate the $C_{60}$ groups contacted by the
electrodes. (a,c,d) Local transmissions at
$E_{\text{LUMO+2}}$ for different contact scenarios,
(b) LUMO+2, (e) Local transmissions at $E_{\text{LUMO+25}}$,
(f) LUMO+25}
\label{fig:loctrans}
\end{figure}

We first consider the local transmission, defined
by Eq.\ (\ref{eqn:2}). Test calculations show that
contributions from the insulating $C_{10}$-PA groups, e.g.
from $C_{10}$-PA to $C_{60}$-$C_{18}$-PA  or
 between pairs of $C_{10}$-PAs, are very small.
Therefore, we restrict the discussions to local
contributions only between $C_{60}$-containing molecules.
The local transmissions were calculated at the energies $E$,
 where the total transmission (cf.\ Fig.\ \ref{fig:moanalysis}) exhibits its maxima.
Each of these maxima can be related to a
single MO in the spectrum. The results for energies $E_{\rm LUMO+2}$ (with
a total transmission of $t= 0.79$) and $E_{\rm LUMO+25}$ ($t= 0.69$)
are shown in Fig. \ref{fig:loctrans}. In the plots,
local transmissions between all $C_{60}$-$C_{18}$-PA
molecules are depicted as arrows. The color code ranging from blue
to orange indicates the local transmission value.
An arrow is only drawn when this value exceeds 0.01.
In this manner, a pathway for the electron through
the monolayer can be visualized. The shape of the pathway can clearly
be related to the shape of the corresponding MO as
shown, e.g.\ in Fig.\ \ref{fig:loctrans} (a,b) for
LUMO+2.
The non-vanishing local transmission components are
mainly located on molecules which exhibit significant
MO density. This is also illustrated for the transmission path
at energy $E_{\rm LUMO+25}$ in Fig.\ \ref{fig:loctrans} (e,f).

By integrating the local transmissions over energy, the
local current is obtained, see Eq.\ (\ref{eqn:localcurrs}).
This corresponds to a weighted sum of all local transmission paths
contributing at a certain voltage. Tuning the bias voltage,
the number of energy levels within the bias
window and thus the number of possible transport channels
is increased. By plotting the local current components, we
can visualize all channels which are open at a specific
 voltage.
Fig.\ \ref{fig:mcpaths} shows the local currents at a bias
voltage of $V=2.4$ V. The local currents are represented
as arrows between pairs of $C_{60}$-$C_{18}$ PAs, where
the color of the arrows indicates the strength of the current.
An arrow is only drawn if the current
is larger than 0.01 nA. At the voltage considered, several transmission
channels contribute. The dominating electron path corresponds to
the contribution of LUMO+25 (cf.\ Fig.\ \ref{fig:loctrans} (e,f)).
In the upper region of the
SAM, smaller contributions of the lower unoccupied orbitals
(especially LUMO+2,+5,+14, cf.\ Fig.\ \ref{fig:loctrans} (a)) are seen.
There are also crossing components of the local currents
 between the different regions.
The pathways obtained based on the local current calculations
are compared in Fig.\ \ref{fig:mcpaths} to the results of
 MC paths searches (yellow colored density) performed using
the local electron affinity energy map
\cite{clark,jaeger} of the system.
In the MC simulation, paths were initialized at the
left boundary of the SAM such that they can start at all five
outermost $C_{60}$ groups in agreement with the contact scenario used in the
local current calculation. In the MC simulation a temperature of
 $T= 2500$ K is used. Computationally this temperature can be attributed to the acceptance
rate of the MC steps and physically it is related to a higher thermal
energy of the electrons (285 meV) experimentally supplied by the
operation voltage of the devices.
We find that the channels of high local currents agree well
with the MC results. There are some minor deviations.
For example, the MC simulation misses the currents that cross
between different regions. On the other hand, the MC simulation
may also find paths localized in between the molecules, which
are by definition not present in the local current that
describes transitions between molecules.
The overall rather good agreement of the results obtained by two
quite different approaches supports the robustness of the results.

\begin{figure*}
	\centering
\includegraphics[width=0.65\linewidth]{current-1.eps}
\caption{Local currents at $V=2.4$ V (scale from blue to
orange in nA) compared to MC transfer paths
through the mixed SAM. The MC paths are shown as
underlying yellow paths
and were calculated at $T= 2500$ K.}
\label{fig:mcpaths}
\end{figure*}

\subsection{Influence of Electrode-SAM Contact Geometry}

We have, furthermore, investigated the influence of
the electrode-SAM contact geometry on the local
transmission paths and the overall current.
This was achieved by varying the self energy (cf.\ Sec.\ \ref{sec:methods_transport}),
which allows the  C$_{60}$ head groups to which the electrodes couple to be selected.
The influence of the contact scenario on transport
is illustrated in Fig.\ \ref{fig:loctrans}
for the local transmission at energy $E_{\rm LUMO+2}$.
The results show that the starting and end
points of the pathways depend strongly on the
C$_{60}$ head group the electrodes couple to. For example,
in Fig.\ \ref{fig:loctrans} (c) and (d), where a smaller number of C$_{60}$ where contacted, the
transmission values are significantly reduced.  Moreover
we find that the transmission is strongly decreased
if the electrodes couple only to molecules far away
from those with MO density.

\begin{figure}[h!]
	\centering
\includegraphics[width=0.95\linewidth]{Plt_curr_contacts_300K.eps}
\caption{ Comparison of currents for the different contact scenarios depicted in Fig. \ref{fig:loctrans} (a, c, d)}
\label{fig:currscont}
\end{figure}

The total current obtained for the three different contact
scenarios depicted in
Fig. \ref{fig:loctrans} (a, c, d) is shown in Fig.\ \ref{fig:currscont}.
At larger voltages, where many transmission channels contribute,
the overall current scales linearly with the
number of molecules contacted.
However, at smaller voltages this is not the case.
For example, the current obtained for contact to only a single
fullerene (cf.\ Fig.\ \ref{fig:loctrans} (d)) is
larger than that for the scenario of Fig. \ref{fig:loctrans}
(c) with three fullerenes in contact.
This shows that for smaller voltage,
the number of contributing channels and thus the current
magnitude depends sensitively on the specific electrode-molecule
contact scenario.

\subsection{Influence of Structural Fluctuations} \label{sec:results_fluc}

All results discussed so far were obtained
for a SAM with fixed geometry, i.e.\ a single snapshot taken at a fixed
time of the MD simulation. However, structural fluctuations may
play an important role for transport processes
in the SAM. For instance, the formation of conductive
pathways depends critically on the molecular arrangement,
e.g.\ on the distance between the fullerene head groups.

To study the importance of structural fluctuations,
we have analyzed the electronic
structure and the  transport
properties for several snapshots taken
at different times along a MD trajectory.
Fig.\ \ref{fig:snaps} shows the energy spectrum, the transmission function and
the current-voltage characteristic for three representative examples,
 taken at $t=16$ ns, $t=25$ ns and $t=44.25$ ns.
The comparison of the  spectra for the different geometries
reveals significant fluctuations of the energy levels.
The energy of the LUMO, for example, fluctuates by about
$0.1$ eV. As a result, the transmission functions change significantly, which,
in turn, results in fluctuations of the current by about one order of magnitude.
These examples show that structural fluctuations can influence the
electronic structure and the transport properties in the C$_{60}$-SAMs significantly.

A more detailed analysis of the influence of
fluctuations on the transport
properties requires, in principle, a theoretical treatment beyond
the methodology employed here, using, e.g.\, time-dependent
density matrix or nonequilibrium Green's function (NEGF)
methods.\cite{Stefanucci2013} This will be the
the subject of future work.

\begin{figure
	\centering
\includegraphics[width=0.7\linewidth]{fluctuations.eps}
		\caption{Energy level spectrum, transmission function
(a) and current-voltage characteristic (b) for different geometries of the
SAM taken at three different
times of the MD simulation.}
\label{fig:snaps}
\end{figure}

\section{Conclusions}

We have investigated charge transport
in $C_{60}$-based SAMs, which are used in SAMFET devises.
The theoretical study was based on a combination of MD simulations to characterize
the structure of the SAMs, semiempirical MO electronic structure calculation and
Landauer transport theory.

The analysis shows that the transport characteristics can be related to
the structural and electronic properties of the SAMs.
Due to the structural disorder of the SAMs, the
molecular orbitals are predominantly localized. As a result, the transmission characteristics
are dominated by narrow peaks, which can
be assigned to the specific molecular orbitals, and give rise to relatively low levels of electrical current,
in qualitative agreement with experiment.
A comparison of SAMs  with different mixing
ratios of $C_{60}$  and alkene groups revealed, furthermore, the strong dependence of the
transport properties on the layer composition. The
conductivity depends crucially  on the formation
of conductive transport channels across the $C_{60}$
head groups, which requires a certain degree of order. Extending previous work,\cite{novak,jaeger}
our studies show  that mixing
$C_{60}$-$C_{18}$-PAs with $C_{10}$-PA increases the
molecular order and facilitates transport. In SAMs consisting of pure
 $C_{60}$-$C_{18}$-PAs, on the other hand, the disordered molecular structure
 hinders the formation of conduction channels and results in a poor performance.

We have also analyzed pathways of charge transport in the SAMs employing the theory of local currents.
The results demonstrate that the electrical current is dominated by few selected pathways,
which facilitate transport through the SAM. The specific contact scenario to the electrodes may influence the
selection of pathways, and thus affect the level of the electrical current.

Another important aspect of transport in organic SAMs is the
effect of structural fluctuation on the transport properties.
Our results show that although the structural fluctuations within the SAM are relatively small,
they have a significant effect on the transport properties and result in fluctuations of
the current by about an order of magnitude. This is
especially due to fluctuations in the alignment of the fullerenes.
To further examine the role of
fluctuations on the transport properties,
a time-dependent treatment of charge transport
\cite{Woiczikowski,croy} is required, using, e.g.,
time-dependent nonequilibrium Green's function (NEGF) theory.\cite{popescu}
This may also open the door to a more complete theoretical treatment, including
inelastic processes due to electron-electron and electron-phonon interaction.

\section*{Acknowledgments}
This work was supported by the
the Deutsche Forschungsgemeinschaft (DFG) through the Cluster
of Excellence 'Engineering of Advanced Materials' (EAM) and
SFB 953. Generous allocation of
computing time at the computing centers in Erlangen (RRZE),
Munich (LRZ), and J\"ulich (JSC) is gratefully acknowledged.

redpajama_set_name: RedPajamaArXiv


text:
\section{Introduction}
\label{sec:intro}
In music production, a song is distributed to the market through a final process called \textit{mastering}, performed by audio engineers \cite{zolzer2011dafx}. During the process, engineers need to consider the musician's desire, sound quality, sound trend, consistency of the album, and distribution format (CD-ROM, half-inch reel tape, PCM 1630 U-Matic tape, etc.). Based on the consideration, they deal with the tone, balance, volume, and space across the song. Examples of possible actions taken during the process is as follows \cite{shelvock2012audio}:

\begin{itemize}
\setlength{\parskip}{-3pt}
    \item Editing minor flaws and applying noise reduction to eliminate clicks, dropouts, hum, and hiss.
    \item Equalizing audio across tracks for optimized frequency distribution.
    \item Adjusting stereo width.
    \item Dynamic range compression or expansion.
    \item Peak limit.
\end{itemize}

To perform the process, engineering knowledge and know-how are crucial, along with an understanding of the overall music production \cite{wyner2013audio}. Hence, quality mastering is difficult to access at the level of general engineers and musicians. In addition, \textit{remastering} is a process of newly mastering past distributed songs to meet current sound trends and distribution formats \cite{wyner2013audio}. For instance, \textit{Queen}'s original vinyl record,  released in 1980, has been periodically remastered to match the sound trends and distribution formats of the time.

Since mastering is necessary for music production but the most intricate part, we seek to lower this entry barrier by using neural networks. Furthermore, we expected that remastering would also be possible if our proposed method can successfully master. To this end, we propose a system that alters the mastering style of a song to the desired reference track as illustrated in Fig.\ref{fig:task objective}.

Examples of the generated audio samples from our model are available at \url{https://dg22302.github.io/MusicRemasteringSystem/}.

\begin{figure}
  \centering
  \includegraphics[width=\linewidth]{Images & Tables/task_outline.png}
  \setlength{\abovecaptionskip}{-5pt}
  \setlength{\belowcaptionskip}{-15pt}
  \caption{Outline of the task objective.}
  \label{fig:task objective}
\end{figure}

\begin{figure*}[t]
  \centering
  \includegraphics[width=\linewidth]{Images & Tables/architecture.png}
  \setlength{\abovecaptionskip}{-5pt}
  \setlength{\belowcaptionskip}{-5pt}
  \caption{Overview of the proposed method. The Mastering Cloner aims to convert the mastering effects of the input track $A1$ to that of the reference track $B2$ by conditioning the encoded feature of the reference track $2'$ extracted from the pre-trained Music Effects Encoder.
  For the training procedure, tracks $A$ and $B$ applied with the same mastering manipulation are used as a ground truth $A2$ and reference track $B2$, where $A1$ is a differently manipulated sample used as the input of the Mastering Cloner.
  The discriminator is applied with the expectation of generating realistic sounds and similar mastering effects to the projected track.}
  \label{fig:model_architecture}
\vspace{-10pt}
\end{figure*}

\section{Related Works}
\label{sec:related}
Intelligent audio mastering is in its infancy due to insufficient training data for the system - lacking of resources, and semantic representation. Figuring out a sound trend requires big data collections and sensibly detectable analysis  \cite{birtchnell2018automating}. Also, aesthetical sound perception is hard to quantify and analyze by computing models \cite{birtchnell2018listening}.
Despite these issues, studies on the intelligent mastering system have been conducted to alleviate the barriers faced by users.
\cite{najduchowski2018automatic} proposed an automatic music mastering system to automatically enhance unprocessed audio signals with the specific parameters obtained from the reference audio. \cite{ramirez2021differentiable} proposed a differentiable approach of using black-box audio plug-ins by estimating its parameters, which allowed users to start with the generalized setup during the mastering procedure.
However, the application of most prior works focuses on applying estimated parameters to unprocessed audio signals.

An end-to-end approach of converting black-box music effects with already processed audio tracks comes in need upon the loss of original dry source tracks or unavailability of replicating the setup of the mastering chain at the time, which may occur especially with old recordings.
\cite{koo2021reverb} introduced a system that interchanges the musical reverberant effects of two differently processed vocal tracks yet required massive storage of already-processed data to train.
We propose an end-to-end remastering system that thoroughly converts the originally mastered effects to the desired style. The model is trained in a self-supervised manner, making it convenient to collect training data where only mastered tracks with its identity are needed.

\section{Methodology}
\label{sec:methodology}
The proposed system shown in Fig.\ref{fig:model_architecture} is trained in a self-supervised manner, where track $A$ and $B$ are two different segments but from the same song. The training procedure is carried out under the assumption that any part of a single song has the equivalent mastering style. Then, by manipulating the mastering style of these tracks, the entire system is trained to transform the input track into the manipulated style of the target. We explain details of each sub modules in this section.

\subsection{Mastering Effects Manipulator}
\label{ssec:mastering manipulator}
Despite the nature of the subjective preference in music mastering and remastering tasks, shared norms and conventions exist \cite{de2017ten}.
Normally, the mastering engineer applies audio effects with analog audio equipment and Virtual Analog or digital plugins on the mixed track. The applied effect sequence is called the \textit{mastering chain} \cite{shelvock2012audio, wyner2013audio}.
The general process of a mastering chain runs as follows. First, the audio engineer makes a tonal adjustment to the entire track. Then, the engineer makes an overall spatial adjustment and harmonious grouping called Gluing \cite{shelvock2012audio}.  Lastly, the final volume of the music is set.
Following this convention, we designed the Mastering Effects Manipulator as the most basic mastering chain by implementing parametric \textit{equalizer} (EQ) for tonal control, multiband \textit{stereo imager} for spatial adjustment, and \textit{maximizer} (or \textit{compressor}) for voluminous control.
For the gain controller, EQ, and maximizer, we follow the implementation to pymixconsole\footnote{\url{https://github.com/csteinmetz1/pymixconsole}} - the collection of audio effects modules presented at \cite{steinmetz2021automatic}, and implement the multiband stereo imager that manipulates stereo width in each band using Linkwitz–Riley crossover filter \cite{linkwitz1976active}.

\subsection{Music Effects Encoder}
\label{ssec:music effects encoder}
Models trained with contrastive learning in a self-supervised manner such as SimCLR \cite{chen2020simple} have recently shown that features extracted from these models contain powerful representations.
The work \cite{spijkervet2021contrastive} applied this methodology to the music domain, where they assigned different sections of the same song as positive pairs and sections from other songs as negative samples as the contrastive objective.
This learning objective makes it extremely convenient for not only training the model but also collecting the dataset since only the song identity of given audio is required.
Accordingly, we trained our Music Effects Encoder with a contrastive objective and expected it would imply the song's overall timbre and mood, including its mastering style.

The Music Effects Encoder consists of multiple 1-dimensional convolutional blocks, where each block includes two convolutional layers with a residual \cite{ he2016deep} connection in between. Each convolutional layer is followed by a batch normalization and a rectified linear unit (ReLU) activation function. The output of the last convolutional layer is time-wise encoded through global average pooling, resulting in the dimensionality of 2048, and is used as the music effects feature $g_{enc}(.)$.
Following the training procedure in \cite{chen2020simple, spijkervet2021contrastive}, $g_{enc}(.)$ is mapped to 512-dimensional features using a linear layer for the contrastive objective, where the loss function used is normalized temperature-scaled cross-entropy loss.

\subsection{Mastering Cloner}
\label{ssec:mastering cloner}
The Mastering Cloner $\psi$ aims to transform the input track $A1$ into $A2'$ by conditioning the encoded reference track $g_{enc}(B2)$.
The network architecture follows that of the Wave-U-Net \cite{stoller2018wave} with some modifications for a finer quality.
First, we adopt the down and upsampling method suggested in \cite{karras2021alias}. The main idea is to prevent aliasing caused by the nonlinear activation functions simply by oversampling. Thus, for each down and upsampling, a leaky ReLU function is applied after the input signal is upsampled (at least by a factor of 2), then is downsampled with a low-pass filter to remove irrelevant frequencies caused during oversampling.
Second, the first layer of the U-Net is computed with a striding factor of 1 without performing skip-connection. This is not only to preserve the high-frequency structure of the input but also to allow opportunities to convert the mastering style with more variation.
Third, conditioning $g_{enc}(B2)$ is performed at the end of each decoder's block with Feature-wise Linear Modulation (FiLM) \cite{perez2018film} operation.
Finally, after the final decoding block, an additional convolutional layer is applied with no activation function, where the output values are trimmed to the range of [-1, +1].

The objective function of the Mastering Cloner is a combination of root mean squared (RMS) loss $\mathcal{L}_{RMS}$ with a weight term $\gamma$ and multi-scale spectral loss $\mathcal{L}_{MSS}$ \cite{engel2020ddsp}.
The purpose of $\mathcal{L}_{RMS}$ is to penalize the volume factor with non-linearity which is defined as
\begin{equation}
\begin{split}
    \gamma &= \rho \cdot \min(\rho^{-1}, |RMS(A2)-RMS(A2')|), \\
    \mathcal{L}_{RMS} &= \mathbb{E}_{A2,A1,B2}[\gamma^{1.5} \cdot \|RMS(A2)-RMS(A2')\|_2^2],  \\
\end{split}
\end{equation}
where $\rho$ is a hyperparameter for $\gamma$.
For $\mathcal{L}_{MSS}$, we utilize the original implementation of \cite{engel2020ddsp} with FFT sizes of (4096, 2048, 1024, 512) and apply it to both stereo-channeled and mid/side-channeled outputs.
Computing mid/side-channeled outputs with multi-scale spectral objective has been proposed as the \textit{stereo loss function} in \cite{steinmetz2021automatic}, which is to separately calculate summation and subtraction of left and right channels with $\mathcal{L}_{MSS}$ so that the phase-related information can also be addressed.
The total objective function for the Mastering Cloner is as follow:
\begin{equation}
\begin{split}
    \mathcal{L}_{\psi} &= \mathcal{L}_{RMS} + \mathcal{L}_{MSS}(A2_{left}, A2'_{left}) + \mathcal{L}_{MSS}(A2_{right}, A2'_{right}) \\
    &+ \mathcal{L}_{MSS}(A2_{mid}, A2'_{mid})
    + \mathcal{L}_{MSS}(A2_{side}, A2'_{side}) .
\end{split}
\end{equation}

\subsection{Projection Discriminator}
\label{ssec:projection discriminator}
Our discriminator $D$ follows the projection discriminator introduced in \cite{miyato2018cgans}.
Similar to the method proposed at \cite{choi2020inference}, we project embeddings extracted from the encoder, the Music Effects Encoder, trained with a self-supervised objective. This encourages the Mastering Cloner to produce more realistic results while conveying the mastering effects of the reference track.

The formation of each encoding block is equivalent to that of the Music Effects Encoder, except it uses 2-dimensional convolutional layers without a residual connection.
The final output of the convolutional blocks is feature-wise averaged pooled, resulting in a 1024-dimensional feature. To match this dimension, the extracted feature from the Music Effects Encoder $g_{enc}(B2)$ is encoded with a 2-layer perceptron and a ReLU activation function in between.

We use hinge version of the standard adversarial loss \cite{lim2017geometric} as adversarial loss:
\begin{equation}
\begin{split}
    \mathcal{L}_{D} &= \mathbb{E}_{A2,B2}[\max(0, 1-D(A2,g_{enc}(B2)))] \\
    &+ \mathbb{E}_{A1,B2}[\max(0, 1+D(A2',g_{enc}(B2)))], \\
    \mathcal{L}_{G} &= -\mathbb{E}_{A1,B2}[D(A2',g_{enc}(B2))].
\end{split}
\end{equation}

\section{Experiment}
\label{sec:typestyle}

\subsection{Dataset}
\label{ssec:subhead}
The Music Effects Encoder and the Mastering Cloner are separately trained with two different music datasets.
Considering the efficiency of contrastive learning, we fit the Music Effects Encoder with a large-scale dataset: the MTG-Jamendo dataset \cite{bogdanov2019mtg}, which includes over 55,000 stereo-channeled audio tracks with a sampling rate of 44.1 kHz.

On the other hand, we evaluated the Mastering Cloner with private collection of already-mastered songs with specifications of stereo-channeled WAV files and a sampling rate of 44.1 kHz.
It includes genres of Pop, Rock, and Hip-hop, released after the 2000s. The number of songs used for the training and validation set is 760 and 28, respectively.
All results shown in this section were computed with the validation set where each sample was manipulated with a fixed random seed for equal comparison between each model.

\begin{figure}
  \centering
  \includegraphics[width=0.85\linewidth]{Images & Tables/volume.png}
  \setlength{\abovecaptionskip}{5pt}
  \setlength{\belowcaptionskip}{-15pt}
  \caption{Voluminal analysis with the Analyzer of Multimeter.}
  \label{fig:volume}
\end{figure}

\begin{figure*}[t]
  \centering
  \includegraphics[width=\linewidth]{Images & Tables/image_low.png}
  \setlength{\abovecaptionskip}{-10pt}
  \setlength{\belowcaptionskip}{-10pt}
  \caption{Four frequency groups - low, low-mid, mid and high - spatial analysis with the Stereo Imager of Ozone 9.}
  \label{fig:stereo}
\end{figure*}

\subsection{Experimental Setups}
\label{ssec:experimental setups}
We first pre-trained Music Effects Encoder following the contrastive training procedure explained at subsection \ref{ssec:music effects encoder}, where two different segments from the same song are used as positive pairs for each batch, and the other segments from other batches are used as negative samples.
The module is trained with a stereo-channeled waveform of 44.1 kHz to capture high-fidelity audio effects and spatial information of the given track.
This is in comparison to the music representation from \cite{spijkervet2021contrastive}, where the input signal of mono-channeled and a sampling rate of 16 kHz is encoded to the dimensionality of 512, which may restrain its performance upon audio effects.
Furthermore, since the network is capable of encoding variable-length inputs into features of fixed size, we randomly segmentized input audio with a duration of $[5, 10]$ seconds.
We set the batch size to 380 using multiple GPUs and trained it for 100 epochs with the Adam optimizer \cite{kingma2014adam} of a learning rate of $2\cdot10^{-4}$. The temperature parameter for the contrastive loss function is set to 0.5.

We freeze Music Effects Encoder for feature extraction of the reference track while training the Mastering Cloner. Although the input length is set to 2.97 seconds during the training procedure, the network can process variable-length inputs due to its fully convolutional nature.
The input for the discriminator is a stereo-channeled log-magnitude spectrogram computed with a Hamming window of 2048 samples with 75\% overlap and a 2048-point fast Fourier transform.
We set the weighting factor $\rho$ used for $\mathcal{L}_{RMS}$ to 100. The batch size during the training procedure is 64, and Adam optimizer with a learning rate of $2\cdot10^{-4}$ is used. We started to train the discriminator after 100 epochs of training the Mastering Cloner.

Detailed configurations of each module, including its pre-trained models and source code, are publicly available\footnote{\url{https://github.com/jhtonyKoo/e2e_music_remastering_system}}.

\subsection{Qualitative Evaluation}
\label{ssec:qual eval}
We qualitatively analyze our results through plug-ins of Logic Pro X (DAW) and Ozone 9 to analyze whether the generated track reflects the sound effects of the reference track.
Specifically, we compare the signals of the network input, output, and ground truth (GT) in terms of three sonic aspects - volume (loudness), space (width), and tone.

For voluminal analysis, we use the volume level meter of the \textit{Analyzer of Multimeter}, with two-volume metrics - RMS and Loudness Unit Full Scale (LUFS).
As shown in Fig.\ref{fig:volume}, both RMS and LUFS of the generated track are almost the same as GT.
The LUFS changes from -6.3 Decibel Full Scale (dBFS) to -8.9 dBFS, where GT's value is -9.0 dBFS.
As we observe the volume change in stereo-channel, the output signal far widens the gap between the loudness of the left and right channel, which mimics its tendency to GT.

For spatial and tonal analysis, we use frequency-band stereo imager - \textit{the Stereo Imager} (Ozone 9).
Here, we divide the input into four frequency groups - low, low-mid, mid, and high. Through analysis, stereo images of the model output in each groups generally resembles that of the GT's than the input. Especially for the frequency group in low, low-mid and mid groups - from 20Hz to 5.1kHz, in Fig.\ref{fig:stereo}, output image spreads more to the left and right spaces than the input's, closer to GT.
For tonal analysis, on the contrary, the output spectral shape above 1kHz vary from GT's, where it drastically reduces.

\begin{table}[t]
\footnotesize
\caption{Model description}
\label{table:model_description}
\centering
\begin{tabular}{|c|c|c|}
\hline
\multirow{2}{*}{\textbf{Model}} & \multicolumn{2}{c|}{\textbf{Method}}                                                                                \\ \cline{2-3}
                                & \begin{tabular}[c]{@{}c@{}}Pre-trained\\ Music Effects Encoder\end{tabular} & \begin{tabular}[c]{@{}c@{}}Projection\\ Discriminator\end{tabular} \\ \hline
model 1                         & $\times$                                                                    & $\times$                                                           \\ \hline
model 2                         & $\circ$                                                                     & $\times$                                                           \\ \hline
model 3                         & $\circ$                                                                     & $\circ$                                                            \\ \hline
\end{tabular}
\end{table}

\begin{table}[t]
\footnotesize
\setlength{\belowcaptionskip}{-5pt}
\setlength{\abovecaptionskip}{+2pt}
\caption{Quantitative results}
\label{table:quant results}
\centering
\begin{tabular}{c|c|c|c|c}
\hline
\textbf{Methods} & \textbf{$\Delta$RMS ($\downarrow$)}   & \textbf{$\Delta$RMS-side ($\downarrow$)}    & \textbf{fw-SNR ($\uparrow$)} & \textbf{STOI ($\uparrow$)}  \\ \hline
input            & 0.0319          & 0.0531          & 19.63           & 86.70\%          \\ \hline
model 1          & 0.0259          & 0.0560          & 14.42           &  78.79\%              \\
model 2          & 0.0226 & 0.0397 & 15.02  & 78.13\% \\
model 3          & \textbf{0.0212}          & \textbf{0.0396}          & \textbf{15.08}           & \textbf{83.26\%}         \\ \hline
\end{tabular}
\vspace{-10pt}
\end{table}

\subsection{Quantitative Evaluation}
\label{ssec:quan eval}
Since there are no other works capable of this particular task, we solely compare three different models of specifications denoted in Table \ref{table:model_description} for the quantitative evaluations.
We measured performance of the volume, stereo width, tone and timbre, and perceptuality through RMS difference of stereo channels, RMS difference of side channels (RMS-side), frequency-weighted segmental Source-to-Noise Ratio (fw-SNR) \cite{ma2009objective}, and short-time objective intelligibility (STOI) \cite{taal2010short}, respectively.
For fw-SNR and STOI, metrics frequently used in speech enhancement tasks, we computed the difference between target mastered track and model output to observe not only the content but also the mastering style.

The results in Table \ref{table:quant results} show the use of pre-trained Music Effects Encoder significantly improved all other metrics except for STOI, whereas training the system with an adversarial term notably enhanced its perceptuality.
From this observation, we can infer that the end-to-end training procedure primarily focuses on preserving the quality of the input track, as conditioning the representations from the pre-trained encoder helps the Mastering Cloner with reproducing the reference track's mastering effects.
This also indicates that pre-training the encoder with the contrastive objective helped capture audio effects-related information.
All of our proposed models could successfully emulate the volume and spatial styles closer to the target than the input, but is less similar in tone and timbre.
Regarding perceptuality, there is a minor difference in STOI metric between model 3 and input, which can be interpreted as being comparable in quality to samples manipulated with the DSP approach.

\begin{figure}[t]
  \centering
  \includegraphics[width=0.65\linewidth]{Images & Tables/listening_test_violin_plot.png}
  \setlength{\abovecaptionskip}{-0pt}
  \setlength{\belowcaptionskip}{-10pt}
  \caption{Listening test violins plots.}
  \label{fig:violinplot}
\vspace{-5pt}
\end{figure}

\subsection{Listening Test}
\label{ssec:listening test}
To account for the perceptual difference of the mastering style, we conducted a listening test using the multiple stimuli with hidden reference and anchor (MUSHRA) protocol using the Web Audio Evaluation Tool \cite{schoeffler2018webmushra}.
17 musicians and sound engineers who are familiar with the concept of mastering participated in the test and were presented with 15 questions.
Each question included a target mastered sample as a reference sample and a hidden reference as a high-anchor (HA), the input of the network as a mid-anchor (MA), an MA sample randomly manipulated with the Mastering Effects Manipulator as a low-anchor (LA), and output of the proposed model (model 3).
Participants were instructed to rate each samples according to similarity of the mastering style to the reference sample on a scale from 0 to 1.

The results of the listening test is plotted as violin plot in Figure \ref{fig:violinplot}.
The median scores in figure order are (0.32, 0.67, 0.71, and 1.00) with a standard deviation of (0.23, 0.23, 0.22, 0.11), respectively.
Even though the scores between our generated samples and MA seems trivial, we found a significant difference with $p<0.001$, with conducted multiple post-hoc paired t-tests with Bonferroni correction \cite{holm1979simple} for each anchor with our generated samples. From this analysis, we believe the model could change the mastering style closer to the target without distorting the input track's content.

\section{Conclusion}
\label{sec:conclusion}
We proposed an end-to-end system where the input audio track is remastered in the mastering style of the reference track. The model was trained without any additional labels except for the identity of the song and in an adversarial manner for a realistic mastering style and a finer quality. We analyzed the results through both qualitative and quantitative approaches and show that the system can successfully generate a sample of results reflecting the mastering style of the reference track.

A couple of ideas can be performed in the future: 1. By adding Gaussian noise to the system's input, the capacity of generating different mastering styles may extend. 2. Including a controllable module upon producing a mastering style that reflects user's preference.

\section{Acknowledgement}
\label{sec:acknowledgement}
This work was supported by Institute of Information \& communications Technology Planning \& Evaluation (IITP) grant funded by the Korea government(MSIT) [NO.2021-0-01343, Artificial Intelligence Graduate School Program (Seoul National University)].

\vfill\pagebreak

\bibliographystyle{IEEEbib}