\section{Introduction
Networks are pervasive in the sciences.  For example, they are used
in ecology to represent food webs and in engineering and computer
science to design high quality internet router connections.  Depending
on the application, one particular graph property may be more
important than another. Oftentimes, a desired property is to have a
connected graph or to optimize a particular metric while constrained
to connected graphs \cite{Kincaid11}.

One of these measures, the Randi\'{c} index\xspace of a graph, developed by Milan
Randi\'{c}, was originally used in chemistry \cite{randic75}. The Randi\'{c} index\xspace
of a graph can be thought of as an assortativity measure. A network is
described as disassortative if high-degree nodes are predominantly
attached to low-degree nodes \cite{li2005towards}. Minimizing the
Randi\'{c} index\xspace, in many instances, will produce a graph with disassortativity
\cite{li2005towards}.  Why is this optimization problem of interest?
Li et al. \cite{li2005towards}, in the design of an internet router
network, found that networks that maximized throughput also had small
values for the Randi\'{c} index\xspace.  In addition, the Randi\'{c} index\xspace has been shown
\cite{Kincaid09} to correlate with synchronization, an important
property in many network applications.  We explore this correlation
for directed graphs in section \ref{sec:netsync}. Optimizing the Randi\'{c} index\xspace is also useful in analyzing imaging data of cortical thicknesses which we discuss in section \ref{sec:ct}. Our focus is to investigate
algorithms that minimize the Randi\'{c} index\xspace of a graph over all connected
realizations while keeping the degrees of the nodes fixed.
However, our complexity results and algorithms also extend
to the case of maximizing the Randi\'{c} index\xspace as well.

\subsection{Notation and Definitions}
\label{sec:notation}

We assume the reader to have a knowledge of graph theory (see, e.g.,
\cite{West}). We let $\mathbb{R}$ and $\mathbb{Z}$ denote the set of real numbers and integers, respectively. We consider an undirected graph, $G =
(N,E)$, which consists of nodes, $N$, and edges, $E$. We assume that our graph is simple, i.e., there are no self-loops and no multi-edges.
The degree of a node is defined as $d_i(G) := |\{j:(i,j) \in E\}|$.
We denote the node-node adjacency matrix by $A(G)$. The \textit{degree sequence} is the list of the degrees of all the nodes in a graph, which
we represent as  $d(G)=(d_1(G),d_2(G),...,d_n(G))$.
A sequence of non-negative integers is considered
\textit{graphic} if it is the degree sequence of a graph.
Degree sequences can correspond to more than one adjacency matrix or graph. We call these graphs different
realizations of the degree sequence. When the particular graph is clear
from context, we omit $G$ in the previous definitions.

Let nodes $u,v \in N$. We say that $u$ and $v$ are \textit{connected} if there exists a path from $u$ to $v$.
A graph is connected if for all $u\in N$ there exists a path to every other node.

\vskip .1 true in
\begin{definition}
The Randi\'{c} index\xspace of a graph $G=(N,E)$ is defined as $$R_\alpha(G)=\sum_{(v_i,v_j)\in E} (d_i \cdot d_j)^{\alpha},$$ where $\alpha \in \mathbb{R} - \{ 0 \}$.
\end{definition}
\vskip .1 true in
{\noindent}
\!\!We consider the case when $\alpha = 1$ which has received considerable attention in a variety of contexts~\cite{gutman1973graph,li2005towards,beichl2008measuring}. For convenience, we define $R(G) = R_1(G)$.
A natural optimization problem is: \\

{\noindent}\textbf{Minimum Randi\'{c} Index Problem.} Given a graphic degree sequence
what is a graph realization with the minimum Randi\'{c} index\xspace with $\alpha=1$? \\

We define the {\it connected minimum Randi\'{c} index\xspace problem} as the minimum Randi\'{c} index\xspace
problem with the additional constraint of minimizing over all connected realizations.

\vskip .1 true in
\begin{definition}
\label{def:pbmatch}
For a graph $G = (N,E)$ and a positive integer vector $b=(b_1,\cdots,b_n) \in
\mathbb{Z}^n$, a {\it perfect $b$-matching} is a subset of edges $M
\subseteq E$ such that for node $v_i \in N$, the degree of $v_i$ in the
graph $(N,M)$ is $b_i$.
\end{definition}
\vskip .1 true in

\noindent
An associated optimization problem is: \\

{\noindent}\textbf{Minimum Weight Perfect b-Matching Problem.}  Given a positive integer
vector $b$, a graph $G=(N,E)$ and a set of edge weights $w:E \rightarrow \mathbb{R}$, find a perfect
$b$-matching with minimum weight.
\vskip .1 true in
In section 3, we will see that the minimum Randi\'{c} index\xspace problem is equivalent to the
minimum weight perfect $b$-matching problem on a complete graph $G$ with an appropriate choice of weights.  We will also show that, using an arbitrary graph $G$, by constraining the minimum weight perfect $b$-matching problem to connected perfect $b$-matchings, the problem becomes NP-Hard.

\subsection{Network Measures of Assortativity}
\label{sec:randic}

The Randi\'{c} index\xspace of a graph was originally defined in chemistry.
In 1975, the chemist Milan
Randi\'{c} \cite{randic75} proposed the index $R_{\alpha}(N)$
for the cases $\alpha = -1 \hbox{ and } \alpha = -1/2$ under the name {\it branching index}. He explained the utility of $R$ in measuring the extent of branching of
the carbon-atom skeleton of saturated hydrocarbons.
His approach ``reveals some inherent relationships between [structures] which can be traced to connectivity" \cite{randic75}. It is sometimes
referred to as the connectivity index by scholars in chemistry~\cite{hansen2003variable}.
Bollob\'{a}s and Erd\"{o}s~\cite{Erdos} generalized this index by allowing $\alpha$ to take
on any non-zero real number.  A survey of results for the Randi\'{c} index can be found in \cite{LSW08}.

The Randi\'{c} index\xspace with $\alpha=1$ is important enough to have been ``discovered'' several times in the literature.
For example, the second Zagreb index, defined by Gutman~\cite{gutman1973graph} and also used in chemistry, is actually $R(G)$.  In 2005, Li et al. \cite{li2005towards} investigated what they called the s-metric of a graph which, seemingly unknown to them, is actually $R(G)$. They used $R(G)$ to differentiate between graph realizations of a given
degree sequence following a power law distribution \cite{li2005towards} in the design of internet router networks.
For a fixed degree distribution they plotted $R(G)$ versus throughput for hundreds of graph realizations.  They noted that $R(G)$ measures the ``hub-like core" of a graph and is maximized when high-degree nodes are connected to other high degree nodes (assortative).  Conversely, the minimum values of $R(G)$ were predominantly associated with networks that maximized throughput and were dissassortative. In 2008, Beichl and Cloteaux investigated how well random networks generated with a chosen $R(G)$ can model the structure of real networks such as the
Internet. The graphs produced by optimizing $R(G)$ resulted in better models than the ones that used simple uniform sampling \cite{beichl2008measuring}.

\subsection{Random Graph Classes}
\label{sec:graphs}

Our computational experiments require random graphs.
We make use of three types of graphs: Erd\H{o}s-R\'{e}nyi, geometric and scale-free. The structure
of these graphs depends on the parameters chosen. \\

{\noindent}\textbf{Erd\H{o}s-R\'{e}nyi  Graphs.}
A number of nodes $n$ and a probability of connection $p$ are chosen. A uniform random number on the interval $[0,1]$ is generated for each possible edge. If the
number generated for an edge is less than $p$ then the edge is added. \\

{\noindent}\textbf{Geometric Graphs.}
A number of nodes $n$ and a radius $r$ is chosen.  Each node $v_i$ is placed uniformly at random in the unit square, giving coordinates $(x_{i},y_{i})$.  We connect nodes $v_i$ and $v_j$ if $(x_i - x_j)^2 + (y_i -y_j)^2 \leq r^2$ \cite{taylor2009contest}.  \\

{\noindent}\textbf{Scale-Free Graphs.}
A preferential attachment algorithm is used to create graphs whose degree sequences
follow a power-law distribution. Following the convention in the literature we will refer
to these graphs as ``scale-free''.  A number of nodes $n$ is chosen. New nodes are added and
connected to existing nodes, based on a probability proportional to the current degree of the nodes,
until you reach $n$ nodes, making it more likely that a new node will be connected to a higher degree node \cite{taylor2009contest}.
The algorithm allows a minimum node degree to be specified. \\

\section {Formulation and Complexity}

In this section, we formulate
the minimum Randi\'{c} index\xspace problem as a \textit{minimum weight perfect
b-matching problem}, which is solvable in polynomial time
\cite{schrijver2003combinatorial}. Note that this problem does not
enforce connectivity. We then show that approximating the minimum Randi\'{c} index\xspace
problem with connectivity is NP-Hard.

\subsection{The $b$-Matching Problem }

Consider a graph $G = (N,E)$, a positive integer vector $b=(b_1,\cdots,b_n) \in
\mathbb{Z}^n$ and $M \subseteq E$, a perfect b-matching.  For a given $b$-matching,
$M$, the graph induced by $M$ is $(N,M)$. We denote the set of perfect $b$-matchings of
a graph $G$ by $\mathcal{P}_b(G)$.  For edge weights $w:E \rightarrow \mathbb{R}$, the
{\it minimum weight perfect $b$-matching problem} is finding the perfect
$b$-matching with minimum weight, i.e., to calculate
\begin{equation}
  \label{eq:perfectb}
  M^\ast(G) := \arg\min\left\{\sum_{e \in M} w(e) : M \in \mathcal{P}_b(G)\right\}.
\end{equation}
For example, let G be the undirected, weighted graph below
\begin{center}
\begin{tikzpicture}[node distance=20mm]
        \tikzstyle{every node}=%
        [%
          fill=blue!50!black!20,%
          draw=blue!50!black,%
          minimum size=8mm,%
          circle,%
          thick,%
          inner sep=0pt%
        ]

        \node (1)              {$v_1$};
        \node (2) [right of=1] {$v_2$};
        \node (3) [below of=1] {$v_3$};
        \node (4) [below of=2] {$v_4$};

        \path [thick] [every node/.style={font=\sffamily\small}]
        		      (1) edge node[above] {3} (2)
                         (1) edge node [xshift= -0.7 cm, yshift = 0.4 cm] {2} (4)
                         (3) edge node [xshift= 0.7 cm, yshift = 0.4 cm] {7} (2)
                         (2) edge node[right] {4} (4)
                         (3) edge node [below]{1}(4);
      \end{tikzpicture} \\
\end{center}
and let $b = (2,1,1,2)$ for nodes $v_1$, $v_2$, $v_3$ and $v_4$ respectively.
We select $b_i$ edges that will connect to the $i$th node and that
will produce the minimum weight. Therefore, the matching induces the graph $G'$ below.
\begin{center}
\begin{tikzpicture}[node distance=20mm]
        \tikzstyle{every node}=%
        [%
          fill=blue!50!black!20,%
          draw=blue!50!black,%
          minimum size=8mm,%
          circle,%
          thick,%
          inner sep=0pt%
        ]

        \node (1)              {$v_1$};
        \node (2) [right of=1] {$v_2$};
        \node (3) [below of=1] {$v_3$};
        \node (4) [below of=2] {$v_4$};

        \path [thick][every node/.style={font=\sffamily\small}]
        		      (1) edge node[above] {3} (2)
                         (1) edge node [xshift= -0.35 cm] {2} (4)
                         (3) edge node [below]{1}(4);

      \end{tikzpicture} \\
\end{center}
Note that for this example the solution $G'$ is the only perfect $b$-matching for $G$.  Let an instance of the minimum Randi\'{c} index\xspace problem be given with a positive integer vector $b
\in \mathbb{Z}^n$ and graph $G=(N,E)$.  Then the graphs induced
by the matchings are feasible subgraph realizations for the minimum
Randi\'{c} index\xspace problem. Also, the edge sets
of the feasible subgraph realizations of the minimum Randi\'{c} index\xspace problem are
perfect $b$-matchings. Thus, the set of feasible $b$-matchings on $G$ is identical to
the set of feasible subgraph realizations on $G$ to the minimum Randi\'{c} index\xspace
problem on $G$.
Thus, to formulate an instance of minimum Randi\'{c} index\xspace problem as a minimum weight perfect
$b$-matching problem, set
\begin{equation}
  \label{eq:w}
  w_{ij} = b_i \cdot b_j.
\end{equation}
Therefore, we can create an instance of a minimum weight perfect
$b$-matching to solve the minimum Randi\'{c} index\xspace problem. Since the $b$-matching
problem can be solved in polynomial time, finding the minimum Randi\'{c} index\xspace of
a graph can also be done in polynomial time.  Note that this method
does not enforce connectivity.

We show that even approximating the connected minimum Randi\'{c} index\xspace problem is
NP-Hard. We first define approximation algorithms (see \cite{Vazirani01} for
further details about approximation algorithms). Let $S \subset \mathbb{R}^n$
and $f:S \rightarrow \mathbb{R}$ be a given feasibility set and objective
function, respectively. Define an $\alpha$-approximation algorithm for
the minimization problem $v^* = \min_{x \in S} f(x)$ as a polynomial
time algorithm that finds a solution $y \in S$ with $f(y) \leq \alpha
v^*$. We say that we can {\it approximate} a minimization problem if
there exists an $\alpha$ such that an $\alpha$-approximation algorithm
exists. Note that $\alpha \geq 1$ is implicit with $\alpha=1$ only if
an exact algorithm exists.

\begin{theorem}
  Approximating the minimum Randi\'{c} index\xspace subject to a connectivity constraint is NP-Hard.
\end{theorem}

{\it Proof. } Recall that a Hamiltonian cycle on $G$ is a tour (set of adjacent
edges or, equivalently, nodes in $G$) that visits each node exactly
once, except for the start node which is equal to the last node on the
tour. We claim the existance of a Hamiltonian cycle on a given graph
is equivalent to the feasibility of the minimum Randi\'{c} index\xspace with connectivity
on a related instance. Recall that an instance of Hamiltonian cycle consists of a graph, so
let such an instance be given with $G = (N,E)$. Now define the vector
$b \in \mathbb{R}^{|N|}$ by setting $b_i = 2$ for $i \in \{1,\ldots,|N|\}$ and
consider the resulting connected minimum Randi\'{c} index\xspace instance using the graph
$G$ and vector $b$.

We first show that if the minimum Randi\'{c} index\xspace instance $(G,b)$ is feasible,
then there is a Hamiltonian cycle on $G$. Suppose there is a feasible
solution $H=(N,F)$, which means $F \subseteq E$, each node $u \in N$
has degree 2, and $H$ is connected. As each node has even degree and
$H$ is connected, there is an Eulerian cycle, $T$, on $H$. We claim
that $T$ is a Hamiltonian cycle on $G$, which means each node is
visited exactly once by $T$ except the start node which is visited
exactly twice. Choose $u \in N$ and note that two edges are adjacent to $u$. Then, because $T$ traverses every arc, the node $u$ is visited. Denote the start node of $T$ by $s \in N$ and consider traversing
$T$ beginning at $s$. If the traversal visits a node $u \in N
\setminus \{s\}$ more then once then an edge was traversed into $u$, a
second distinct edge was traversed out of $u$, and a third distinct
edge was traversed into $u$, a contradiction as there are exactly two
distinct edges adjacent to $u$ in $H$. The same argument applies if
$s$ is visited more then once before the traversal is complete. So
each node is visited exactly once by $T$ except the start node, which
is returned to when the traversal is complete, i.e., $T$ is a
Hamiltonian cycle. Thus, if the instance $(G,b)$ is feasible, then $G$
possess a Hamiltonian cycle.

We now show that if there is a Hamiltonian cycle on $G$, then the
minimum Randi\'{c} index\xspace problem on $(G,b)$ is feasible. Consider a Hamiltonian
cycle, $C$ on $G$ and the subgraph induced by $C$. Such a subgraph is
connected as each node is visited. Also, each node has degree 2 as
each node $u \in N$ has one arc used to enter $u$ and exactly one
distinct arc $u$ to exit. Thus, if $G$ possess a Hamiltonian cycle,
then $(G,b)$ must be feasible to the given connected minimum Randi\'{c} index\xspace
instance.

Now suppose there were an $\alpha$-approximation algorithm to the
connected minimum Randi\'{c} index\xspace problem for some $\alpha \geq 1$. If the
algorithm returns a solution to the instance $(G,b)$ then $G$ posesses
a Hamiltonian cycle. If it does not, then the instance $(G,b)$ was not
feasibile and $G$ does not posess a Hamiltonian cycle. Note that the
argument does not rely on what value $\alpha$ is.  $\ddag$

The following corollaries are immediate.
\begin{corollary}
  Approximating $b$-matching when the graph induced by the matching must
  be connected is NP-Hard.
\end{corollary}

Because the proof showed that finding a feasible solution to the minimum Randi\'{c} index\xspace index is NP-Hard, we immediately can state the following corollary.
\begin{corollary}
  Approximating the maximum Randi\'{c} index\xspace subject to a connectivity constraint is NP-Hard.
\end{corollary}

Note that we have not shown what the complexity is when the input
graph $G$ is the complete graph. We leave this as an open problem.

\subsection{Example Transformation} \label{subset:re} Given the degree
sequence $d=(3,2,2,2,2,1)$, what is a graph realization with the minimum Randi\'{c} index\xspace?
We let nodes $v_1,v_2,v_3,v_4,v_5, v_6 \in N$ with $b=(3,2,2,2,2,1)$. Now we can form the complete graph $G$,
with weights corresponding to $b_i\cdot b_j$ for every node $v_{i},v_{j} \in
N$.
\begin{center}
\begin{tikzpicture}[node distance=30mm]
        \tikzstyle{nnode}=%
        [%
          fill=blue!50!black!20,%
          draw=blue!50!black,%
          minimum size=8mm,%
          circle,%
          thick,%
          inner sep=0pt%
        ]

        \node[nnode] (a)              {$3$};
        \node (la) [xshift = -.7 cm]{\textcolor{red}{$v_1$}};
        \node[nnode]  (b) [right of=a] {$2$};
        \node (lb) [right of= a, xshift = .7 cm]{\textcolor{red}{$v_2$}};
        \node [nnode] (c) [below of=a, xshift= -1.5 cm] {$2$};
        \node (lc) [below of=a, xshift = -2.2 cm]{\textcolor{red}{$v_3$}};
         \node (lg) [below of=a, xshift = -3.2 cm]{$G:$};
        \node [nnode] (d) [below of=b, xshift=1.5 cm] {$2$};
        \node (ld) [below of=b, xshift = 2.2 cm]{\textcolor{red}{$v_4$}};
        \node (le) [below of=c, xshift = .7 cm]{\textcolor{red}{$v_5$}};
        \node[nnode]  (e) [below of=c, xshift = 1.5 cm] {$2$};
        \node[nnode]  (f) [below of=d, xshift = -1.5 cm] {$1$};
        \node (lf) [below of=d, xshift = -.7 cm]{\textcolor{red}{$v_6$}};

         \path [thick][every node/.style={font=\sffamily\small}]

        		      (a) edge node[above] {6} (b)
		      (b) edge node[above, xshift=.25 cm] {4} (d)
		      (d) edge node[below, xshift=.25 cm] {2} (f)
		      (e) edge node[below] {2} (f)
		      (e) edge node[below, xshift=-.25 cm] {4} (c)
		      (c) edge node[above, xshift=-.25 cm] {6} (a)

		      (a) edge node[above, xshift=-1.25cm, yshift=.75cm] {6} (d)
		      (a) edge node[xshift=-.8cm, yshift=2cm] {3} (f)
		      (a) edge node[xshift=.2cm, yshift=1.8cm ]{6} (e)

		      (b) edge node[xshift=-.2cm, yshift=1.75cm] {2} (f)
		      (b) edge node[xshift=.8cm, yshift=2.1cm] {4} (e)
		      (b) edge node[xshift=1.25cm, yshift=1.1cm] {4} (c)

		      (c) edge node[xshift=-2cm, yshift=.25cm] {4} (d)
		      (c) edge node[xshift=-1.25cm, yshift=1.1cm] {2} (f)

		      (d) edge node[xshift=1.25cm, yshift=1.1cm] {4} (e);

\end{tikzpicture} \\
\vskip .5 true in
\end{center}
  $$v_1~~ v_2~~ v_3 ~~v_4 ~~v_5 ~~v_6 $$
$$  \begin{bmatrix}
    0 & 6 & 6 & 6 & 6 & 3 \\
    6 & 0 & 4 & 4 & 4 & 2 \\
    6 & 4 & 0 & 4 & 4 & 2 \\
    6 & 4 & 4 & 0 & 4 & 2 \\
    6 & 4 & 4 & 4 & 0 & 2 \\
    3 & 2 & 2 & 2 & 2 & 0 \\
  \end{bmatrix}$$
Now we solve the minimum weight perfect $b$-matching for $G$ and obtain $G'$:
\begin{center}
\begin{tikzpicture}[node distance=30mm]
        \tikzstyle{nnode}=%
        [%
          fill=green!50!black!20,%
          draw=green!50!black,%
          minimum size=8mm,%
          circle,%
          thick,%
          inner sep=0pt%
        ]

        \node[nnode] (a)              {$3$};
        \node (la) [xshift = -.7 cm]{\textcolor{red}{$v_1$}};
        \node[nnode]  (b) [right of=a] {$2$};
        \node (lb) [right of= a, xshift = .7 cm]{\textcolor{red}{$v_2$}};
        \node [nnode] (c) [below of=a, xshift= -1.5 cm] {$2$};
        \node (lc) [below of=a, xshift = -2.2 cm]{\textcolor{red}{$v_3$}};
           \node (lg) [below of=a, xshift = -3.2 cm]{$G':$};
        \node [nnode] (d) [below of=b, xshift=1.5 cm] {$2$};
        \node (ld) [below of=b, xshift = 2.2 cm]{\textcolor{red}{$v_4$}};
        \node (le) [below of=c, xshift = .7 cm]{\textcolor{red}{$v_5$}};
        \node[nnode]  (e) [below of=c, xshift = 1.5 cm] {$2$};
        \node[nnode]  (f) [below of=d, xshift = -1.5 cm] {$1$};
        \node (lf) [below of=d, xshift = -.7 cm]{\textcolor{red}{$v_6$}};

         \path [thick][every node/.style={font=\sffamily\small}]

        		      (a) edge node[above] {6} (b)
		      (c) edge node[above, xshift=-.25 cm] {6} (a)

		      (a) edge node[xshift=-.7cm, yshift=2cm] {3} (f)

		      (b) edge node[xshift=.8cm, yshift=2.1cm] {4} (e)

		      (c) edge node[xshift=-2cm, yshift=.25cm] {4} (d)

		      (d) edge node[xshift=1.25cm, yshift=1.1cm] {4} (e);

\end{tikzpicture} \\
\vskip .5 true in
\end{center}
$$  \begin{bmatrix}
    0 & 1 & 1 & 0 & 0 & 1 \\
    1 & 0 & 0 & 0 & 1 & 0 \\
    1 & 0 & 0 & 1& 0 & 0 \\
    0 & 0 & 1 & 0 & 1 & 0 \\
    0 & 1& 0 & 1 & 0 & 0 \\
    1 & 0 & 0 & 0 & 0 & 0 \\
  \end{bmatrix}$$
$$R(G') = 6+6+4+4+4+3 = 27 $$
$G'$ is a solution for the minimum weight perfect $b$-matching.
The sum of the weights is the minimum Randi\'{c} index\xspace and the unweighted adjacency
matrix is the corresponding graph realization. Note that there are other
solutions to the matching that will produce the minimum Randi\'{c} index\xspace and a different
realization. That is, the solution is not unique.

\section {Algorithms}

Our primary goal is to devise an algorithm to solve the minimum Randi\'{c} index\xspace
problem.  An algorithm that is useful when creating graphs with a
specified degree sequence is the Havel-Hakimi algorithm.  The
Havel-Hakimi algorithm can be used to check if a degree sequence is
graphic and to find a realization of that sequence.

\subsection{Havel-Hakimi Algorithm}
\label{sec:havel}
The Havel-Hakimi algorithm is useful when we have a non-negative integer sequence and we want to know if it is graphic, and if so, what is a realization of the sequence.  Its algorithm is given below:\\

\def\norm#1{\left\|#1\right\|}
\newcommand{\fprod}[1]{\ensuremath{\left\langle #1 \right\rangle}}
\begin{center}
   \fbox{
     \begin{minipage}{.95\linewidth}
       \begin{codebox}
         {\Procname{\underline{\proc{Havel-Hakimi}} \cite{havel1955remark,hakimi1962realizability}}}
         \zi {\bf Inputs: } $d$, a non-negative integer sequence\\
         \zi {\bf Outputs: } $A$, adjacency matrix of a graph realization of $d$ (if graphic) \\
         \zi Initialize a $|d| \times |d|$ adjacency matrix $A$ so that $A_{ij}=A_{ji} = 0$.
         \zi \While $d$ is not the 0 sequence
         \zi \Do
         \zi Pick a random index $i$ of $d$.
         \zi Subtract $1$ from the $d_i$ nodes with largest degrees (not including $d_i$) and set $d_i$ to 0.

         \zi For the $d_i$ nodes of largest degrees in the previous step (call them $v_j$), set
         \zi \ \ \ $A_{ij} = A_{ji} = 1$.
	  \zi \If $|d| == 1$
	  \zi \ \ \ \ \ $d$ is not graphic.
	  \zi \ \ \ \ \ \Return null.
         \zi \End
         \zi \Return $A$
         \label{alg:HH}
       \end{codebox}
     \end{minipage}
   }
   \vspace*{+.1in}
\end{center}

\subsection{Two-Switches and the Metagraph}
\label{sec:twoswitch}
One way to generate a collection of realizations for a degree sequence is to move
from one realization to another by doing a two-switch, with an example as follows:
\begin{center}
\begin{tikzpicture}[node distance=15mm]
        \tikzstyle{every node}=%
        [%
          fill=green!50!black!20,%
          draw=green!50!black,%
          minimum size=8mm,%
          circle,%
          thick,%
          inner sep=0pt%
        ]

        \node (A)              {$a$};
        \node (C) [right of=A] {$c$};
        \node (B) [below of=A] {$b$};
        \node (D) [below of=C] {$d$};

        \path [thick](A) edge (B)
                           (C) edge (D);

      \end{tikzpicture}
~~~~~~~~~~~~~
\begin{tikzpicture}[node distance=15mm]
        \tikzstyle{every node}=%
        [%
          fill=green!50!black!20,%
          draw=green!50!black,%
          minimum size=8mm,%
          circle,%
          thick,%
          inner sep=0pt%
        ]

        \node (A)              {$a$};
        \node (C) [right of=A] {$c$};
        \node (B) [below of=A] {$b$};
        \node (D) [below of=C] {$d$};

        \path [thick](A) edge (B)
                           (C) edge (D);
        \path [dashed](A) edge (D)
                         (C) edge (B);

      \end{tikzpicture}
~~~~~~~~~~~~~
\begin{tikzpicture}[node distance=15mm]
        \tikzstyle{every node}=%
        [%
          fill=green!50!black!20,%
          draw=green!50!black,%
          minimum size=8mm,%
          circle,%
          thick,%
          inner sep=0pt%
        ]

        \node (A)              {$a$};
        \node (C) [right of=A] {$c$};
        \node (B) [below of=A] {$b$};
        \node (D) [below of=C] {$d$};

        \path [thick](A) edge (D)
                         (C) edge (B);

      \end{tikzpicture}
\end{center}
When doing a two-switch, we examine two edges, $(a,b),(c,d)\in E$.  If $(a,d)\notin E$ and $(b,c)\notin E$ then we can remove edges $(a,b)$ and $(c,d)$ and create
edges $(a,d)$ and $(b,c)$. This is not a unique move, since we could also use $(a,c)$ and $(b,d)$ if $(a,c)\notin E$ and $(b,d)
\notin E$. Two-switching is an easy way to obtain a different graph with the same degree sequence after a graph is created using
the Havel-Hakimi algorithm.  The two-switch pseudocode is given below:\\
\begin{center}
   \fbox{
     \begin{minipage}{.85\linewidth}
       \begin{codebox}
         {\Procname{\underline{\proc{Two-Switch}}}}
         \zi {\bf Inputs: }$G$, a graph with degree sequence $d$.
         \zi {\bf Outputs: }$G^{\prime}$, a new graph with degree sequence $d$.
         \zi Let $G^{\prime} := G$.
         \zi Pick a random edge $(a,b)\in E(G^{\prime})$.
         \zi Find a node $c\in N(G^{\prime})$ that is not connected to $b$.
         \zi Find a node $d\in N(G^{\prime})$ that is connected to $c$, but not to $a$.
         \zi Remove edges $(a,b),(c,d)$ from $E(G^\prime)$ and add $(a,d),(c,b)$ to $E(G^\prime)$.
         \zi \Return $G^{\prime}$
         \label{alg:tw}
       \end{codebox}
     \end{minipage}
   }
   \vspace*{+.1in}
\end{center}
We can construct a {\it metagraph} of a degree sequence, where
the metagraph is an undirected graph with each node representing a graph realization of a degree sequence and each edge representing a two-switch.  The following theorem shows that the metagraph is always a connected graph.
\begin{theorem} (Ryser's Theorem \cite{ryser1957combinatorial}.)
Given graphs $G$ and $G'$ such that $d(G)=d(G')$, there exists a sequence of two-switches going from $G$ to $G'$.
\end{theorem}

\subsection{Heuristic for Disconnected Realizations}
\label{sec:heuristic}

Since finding a connected graph realization with minimum Randi\'{c} index\xspace is NP-Hard,
we present a heuristic using two-switches to connect disconnected realizations. The heuristic sequentially performs a two-switch
between pairs of connected components until all the components are connected:
\begin{center}
   \vspace*{+.1in}
   \fbox{
     \begin{minipage}{.85\linewidth}
       \begin{codebox}
         {\Procname{\underline{\proc{Two-switch Heuristic}}}}
         \zi {\bf Inputs: }$G$, a disconnected graph with degree sequence $d$
         \zi {\bf Outputs: }$G^{\prime}$, a connected graph with degree sequence $d$
         \zi
         \zi Let $G^{\prime} := G$.
         \zi \While the number of connected components in $G^{\prime}$ is $\geq 2$
         \zi \Do a two switch with two components to connect them using two randomly
         \zi chosen edges from each component
         \zi \End
         \zi \Return $G^{\prime}$
       \end{codebox}
     \end{minipage}
   }
   \vspace*{+.1in}
\end{center}

Note that the method to connect the disconnected realizations may not produce graphs
with the best structure since there is only 1 edge connecting one component to another.
Also note that we do not need to check whether the randomly chosen edges are adjacent or
not since they are in separate connected components.

\section{Solving the Minimum Randi\'{c} index\xspace Problem}

In this section, we focus on the case where the input graph $G$ is the
complete graph.
To solve the minimum Randi\'{c} index\xspace problem we used code that solves a minimum weight perfect $b$-matching problem.
The code used is for generalized matching problems and was written by Vlad Schogolev, Bert Huang, and Stuart Andrews.
Their code uses the GOBLIN graph library (http://goblin2.sourceforge.net/). Huang's paper on loopy belief propagation for bipartite maximum weight
$b$-matching uses this code \cite{huang2007loopy}.
The code solves a \textit{maximum} weight perfect $b$-matching problem given the weight matrix, $b$ vector and the exact solution algorithm choice.
After transforming our minimum Randi\'{c} index\xspace problem instance into a minimum weight perfect
$b$-matching instance, the $b$ vector will correspond to the degree constraints and the weight matrix to the possible degree products (see Section
\ref{subset:re}). However, since the code solves the $maximum$ matching we transform our problem so that solving for the maximum yields the solution for the minimum. \\

Given a weight matrix $H$ we transform these weights into a matrix $H_2$ such that the maximum matching using $H_2$ will yield the same
solution as the minimum matching using $H$. To do this we take a matrix $M$ with 1s in all positions except for the diagonal which has 0s. We
then multiply every entry by one more than the maximum entry of $H$. $H$ is then subtracted from $M$ yielding $H_2$.
The following pseudocode implements this algorithm and will solve the minimum Randi\'{c} index\xspace problem for a given degree sequence:
\begin{center}
   \vspace*{+.1in}
   \fbox{
     \begin{minipage}{.85\linewidth}
       \begin{codebox}
         {\Procname{\underline{\proc{Algorithm to solve minimum Randi\'{c} index\xspace with $b$-matching }}}}
         \zi {\bf Inputs: }$A$, an adjacency matrix with degree sequence $d$.
         \zi {\bf Outputs: }$A^{\prime}$, the new adjacency matrix with degree sequence $d$
         \zi  ~~~~~~~~~~~~~~~and minimized Randic index $r$.
         \zi
         \zi Create weight matrix $H$ of degree products.
         \zi Let $b := d$.
         \zi Let $H_2 := (1+\max H)\cdot M - H$, where $M$ is the adjacency matrix of all 1s except
         \zi \ \ \ along the diagonal.
         \zi Call b-match solver with $b$ and $H_2$ to get adjacency matrix $A^{\prime}$ of optimal solution
         \zi Calculate $r=R(A^{\prime})$
         \zi \Return $A^{\prime}$ and $r$
       \end{codebox}
     \end{minipage}
   }
   \vspace*{+.1in}
\end{center}
This algorithm returns the minimum Randi\'{c} index\xspace of a graph and a realization.
We know that the $b$-matching code runs in polynomial time (\cite{schrijver2003combinatorial}).
It is easy to see that the transformation steps are done in
polynomial time as well. We used three types of randomly generated graphs to test the algorithm performance:
Erd\H{o}s-R\'{e}nyi, geometric and scale-free.  We limited our computational experiments to degree sequences for which connected realizations
were known to exist. The Randi\'{c} index\xspace before and after the optimization was recorded.   After the optimization we check
if the graph realization with the minimum Randi\'{c} index\xspace is connected.
We generated graphs with 25, 50 and 100 nodes. In addition, 100 of each graph type and size were generated.
Tables \ref{tab:25}, \ref{tab:50} and \ref{tab:100}  present  results from the runs. Note that the number of graphs connected
after the run plus the number of graphs disconnected plus the number of graphs with no connected realizations is 100 for each graph type.
\begin{table}[h!]
  \centering
  \begin{tabular}{c|c|c|c}
{\bf    Graph type} & {\bf connected} & {\bf disconnected} & {\bf no connected realizations}\\
    \hline
    Erd\H{o}s-R\'{e}nyi  & 67 & 1 & 32\\
    Geometric & 61 & 2 & 37\\
    Scale-Free & 93 & 7 & 0\\
    \hline
  \end{tabular}
  \caption{25 node graphs}
  \label{tab:25}
\end{table}
\begin{table}[h!]
  \centering
  \begin{tabular}{c|c|c|c}
{\bf    Graph type}&{\bf connected} & {\bf disconnected} & {\bf no connected realizations}\\
    \hline
    Erd\H{o}s-R\'{e}nyi  & 50 & 5 & 45\\
    Geometric & 57 & 3 & 40\\
    Scale-Free & 85 & 15 & 0\\
    \hline
  \end{tabular}
  \caption{50 node graphs}
  \label{tab:50}
\end{table}
\begin{table}[h!]
  \centering
  \begin{tabular}{c|c|c|c}
{\bf    Graph type}&{\bf connected} &{\bf disconnected} &{\bf no connected realizations}\\
    \hline
    Erd\H{o}s-R\'{e}nyi  & 16 & 2 & 82\\
    Geometric & 30 & 6 & 64\\
    Scale-Free & 91 & 8 & 1\\
    \hline
  \end{tabular}
  \caption{100 node graphs}
  \label{tab:100}
\end{table}
The MATLAB functions used to generate the geometric and scale-free graphs are from CONTEST:
A Controllable Test Matrix Toolbox for MATLAB \cite{taylor2009contest}.
In addition, the necessary and sufficient conditions for a non-negative integer sequence $\{a_i\}$ to be
realizable as the degrees of the nodes of a connected graph are that $a_i \neq 0$ for all $i$
and the sum of the integers $a_i$ is even and not less than $2(n-1).$ This condition was used to
discard graphs with a degree sequence that had no connected realizations \cite{chen1997}.\\

In general from our runs, the realization of the minimum Randi\'{c} index\xspace was connected. There are minimum Randi\'{c} index\xspace
graph realizations that are disconnected and we do not know if there are other realizations with this
Randi\'{c} index\xspace that are connected since the $b$-matching solver only produces one solution. But there were often
a large proportion of graphs that had no connected realization at all. This largely depends on parameters
chosen for the randomly generated graphs.   If the random graph produced has most nodes with large degrees
then it is unlikely that any graph realization would be disconnected. We were interested in generating graphs
that have both connected and disconnected realizations and investigating whether the realization generated with
the minimum Randi\'{c} index\xspace was connected or not.\\

The following parameters were chosen after extensive experimentation so that the number of instances with no
connected realizations was small.
For the Erd\H{o}s-R\'{e}nyi  graphs we used an average degree per node of 4.25. The corresponding $p$ values used were calculated
using $p = \frac{4.25}{n}$ where $n$ is the number of nodes in the graph. Thus $p=.17$ for $n=25$, $p=.085$ for
$n=50$, and $p=.043$ for $n=100$. For the geometric graphs we used an average degree per node of 6. The radii
were calculated using $r=\sqrt{\frac{6}{\pi n}}$. Our corresponding radii were $r=.276$ for $n=25$, $r=.195$
for $n=50$, and $r=.138$ for $n=100$. We used scale-free graphs with a minimum node degree of 2.\\

The left box plots for each of 25, 50 and 100 node graphs in Figures \ref{fig:er3}, \ref{fig:geo3} and \ref{fig:pref3} show the percent difference between
the graph's original Randi\'{c} index\xspace and the minimum Randi\'{c} index\xspace. The percent difference is calculated from $\frac{original - minimum} {minimum} \times 100.$
The right box plots for each of 25, 50 and 100 node graphs in Figures \ref{fig:er3}, \ref{fig:geo3}, and \ref{fig:pref3} show the percent difference between the minimum Randi\'{c} index\xspace and the Randi\'{c} index\xspace after the heuristic algorithm in Section~\ref{sec:heuristic} was applied.
This percent
difference is calculated with $\frac{after~heuristic - minimum} {minimum} \times 100.$ The number of graphs that used the heuristic
depended on the number of optimal graph realizations that were disconnected. Note that this is a different number for each graph type
and size. See Tables \ref{tab:25}, \ref{tab:50} and \ref{tab:100} for those numbers.

 \begin{figure}[!ht]
   \centering
   \resizebox{.6\textwidth}{!}{\includegraphics{er_comp.png}}
      \caption{Comparing percent differences for Erd\H{o}s-R\'{e}nyi  graphs.}
   \label{fig:er3}
\end{figure}

\begin{figure}[!h]
   \centering
   \resizebox{.6\textwidth}{!}{\includegraphics{geo_comp.png}}
      \caption{Comparing percent differences for geometric graphs.}
   \label{fig:geo3}
\end{figure}

\begin{figure}[!htb]
   \centering
   \resizebox{.6\textwidth}{!}{\includegraphics{pref_comp.png}}
      \caption{Comparing percent differences for scale-free graphs.}
   \label{fig:pref3}
\end{figure}

\section{Applications}
\label{sec:apps}

We will now show results for two applications where the effects of network connectivity measures have shown or are hypothesized to play an important role:  neuronal synchronization and dementia.

\subsection{Neuronal synchronization}
\label{sec:netsync}

The first example is exploring the effect of assortativity on the synchronous firing of neurons in the preB\"{o}tzinger complex.  This collection of neurons is responsible for the control of respiration in mammals \cite{Feldman:2006p907}.  As neuronal networks are inherently directed, we must first define a directed Randi\'{c} index\xspace \cite{ZamoraLopez:2008p6831} and extend the $b$-matching algorithm to the directed case.  The extension is similar to \cite{liuNature2011} in which the maximum matching problem is extended to directed graphs for a network controllability problem.

\subsubsection{Extension to directed networks}

Consider the directed graph $\vec{G} = (N,E)$ with vertex set $N = \{v_{1},\ldots,v_{n}\}$ and edge (or arc) set $E = \{(v_{i},v_{j}) \, | \, v_{i} \rightarrow v_{j}\}$.  The degree sequence for $\vec{G}$ is a non-negative integer-pair sequence $d = \{(d^{+}_{i},d^{-}_{i}) \, | \, i=1,\ldots,n\}$, where we denote the out-degree and in-degree sequences by $d^{+}$ and $d^{-}$, respectively.  There are four different Randi\'{c} index\xspace-type measures \cite{ZamoraLopez:2008p6831} given by
\begin{equation}
\label{eq:dirandic}
R^{pq}(G) = \sum_{(v_i,v_j)\in E}d^{p}_{i}d^{q}_{j},
\end{equation}
where $p,q\in\{-,+\}$.  This can be seen as a natural extension of $R(G)$ to the directed case.  In this context, we will now define an extension of perfect $b$-matching as follows.
\begin{definition}[Perfect $b$-Matching for Directed Graphs]
\label{def:dir_pbmatch}
For a directed graph $\vec{G} = (N,E)$ and positive integer-pair sequence $b = (b^{+},b^{-}) = \{(b^{+}_i,b^{-}_{i}) \in \mathbb{Z}^n \times Z^n \, | \, i=1\,\ldots,n\}$, a {\it perfect $b$-matching} is a subset of edges $M
\subseteq E$ such that for node $v_i \in N$, the out and in-degree of $v_i$ in the subgraph $(N,M)$ is $b^{+}_i$ and $b^{-}_{i}$, respectively.
\end{definition}
\noindent In an analogous way as the undirected case, we can define a minimum weight perfect b-matching problem as follows.
\vskip 0.1in
{\noindent}\textbf{Minimum Weight Perfect b-Matching Problem.}  Given a positive integer-pair sequence $b$, a directed graph $\vec{G}=(N,E)$ and a set of edge weights $w:E \rightarrow \mathbb{R}$, find a perfect $b$-matching with minimum weight.
\vskip 0.1in
In order to use the existing algorithm for undirected graphs, we consider the equivalent bipartite form of $\vec{G}$ given by $\vec{G}^{*} = (N^{*},E^{*})$, where $N^{*} = N^{+}\cup N^{-}$, $N^{+} =  \{v_{1}^{+},\ldots,v_{n}^{+}\}$, $N^{-} =  \{v_{1}^{-},\ldots,v_{n}^{-}\}$, and $e^{*} = (v_{i}^{+},v_{j}^{-}) \in E^{*}$ if and only if $e = (v_{i},v_{j}) \in E$ (see Figure~\ref{fig:digraph_example} for an example).
\begin{figure}
\centering
\includegraphics{digraph_example}
\caption{{\bf Left}: Directed graph $\vec{G}=(N,E)$ with node set $N=\{v_{1},v_{2},v_{3},v_{4}\}$ and degree sequence $\{(2,1),(2,1),(1,2),(1,2)\}$.  {\bf Right}:  Bipartite representation $\vec{G}^{*} = (N^{*},E^{*})$, where $N^{*} = N^{+}\cup N^{-}$, $N^{+} =  \{v_{1}^{+},\ldots,v_{n}^{+}\}$, $N^{-} =  \{v_{1}^{-},\ldots,v_{n}^{-}\}$, and $e^{*} = (v_{i}^{+},v_{j}^{-}) \in E^{*}$ if and only if $e = (v_{i},v_{j}) \in E$.}
\label{fig:digraph_example}
\end{figure}
We can define a perfect $b$-matching on the corresponding bipartite graph $\vec{G}^{*} = (N^{*},E^{*})$ as a subset of edges $M^{*}\subseteq E^{*}$ such that for node $v_i \in N^{+}$ or $N^{-}$, the degree of node $v_i$ in $(N^{*},M^{*})$ is $b^{+}_{i}$ or $b^{-}_{i}$, respectively.  This modified definition is a special case of the undirected version in Definition \ref{def:pbmatch} with node set $N = (v_{1}^{+},\ldots,v_{n}^{+},v_{1}^{-},\ldots,v_{n}^{-})$ and the positive vector $b = (b_{1}^{+},\ldots,b_{n}^{+},b_{1}^{-},\ldots,b_{n}^{-})$.  Thus, to find a minimum weight perfect $b$-matching for a directed graph $\vec{G}$, we do the following (see Figure~\ref{fig:bigraph}):
\begin{figure}
\centering
\includegraphics{digraph}
\caption{Bipartite setup of the perfect $b$-matching problem for the directed graph given in Figure~\ref{fig:digraph_example}, with node set $N = (v_{1}^{+},\ldots,v_{4}^{+},v_{1}^{-},\ldots,v_{4}^{-})$ and the positive vector $b = (b_{1}^{+},\ldots,b_{4}^{+},b_{1}^{-},\ldots,b_{4}^{-}) = (2,2,1,1,1,1,2,2)$.  For minimum weight perfect $b$-matching, we would set $w_{ij} = d_{i}^{p}d_{j}^{q}$, where $p,q\in\{+,-\}$.}
\label{fig:bigraph}
\end{figure}
\begin{itemize}
\item Let $B$ be the complete bipartite graph $K_{n,n}$ minus the edges $\{(v_{i}^{+},v_{i}^{-}) \, | \, i=1,\ldots,n\}$.
\item Let $b^{+} = d^{+}$ and $b^{-} = d^{-}$.
\item For $p,q\in\{+,-\}$, let edge weights $w_{ij} = d_{i}^{p}d_{j}^{q}$.
\end{itemize}

\subsubsection{Neuronal networks}

Using NeuronetExperimenter \cite{NeuronetExperimenter:ra0SJb3l}, we simulated 150 rhythmogenic neurons in the preB\"{o}tzinger complex using the Rubin-Hayes neuron model \cite{Rubin:2009ie} (\url{http://senselab.med.yale.edu/modeldb/ShowModel.asp?model=125649}).  It is unknown what degree distribution and network connectivity arise in neuronal networks.  However, it seems reasonable to expect that neurons closer to each other are more likely to be connected.  Hence, we have chosen to model these networks with 3D geometric (directed) graphs.  Raster plots of simulation results for two 3D geometric networks are displayed in Figures~\ref{fig:neuron5} and \ref{fig:neuron20} (center columns), along with simulation results from realizations with the minimum and maximum directed Randi\'{c} indices in \eqref{eq:dirandic}.  These results display an overall tendency for the realization with minimum and maximum Randi\'{c} index\xspace, regardless of type, to lead to a faster and slower breathing rhythm, respectively.  This is not always the case, however, as can be seen with $R^{+\,-}$ in Figure~\ref{fig:neuron20}.  There is also a tendency for more and less synchronous firing from realizations with respective minimum and maximum Randi\'{c} index\xspace.  These results warrant a more thorough quantitative analysis, which is beyond the scope of this work.

\begin{figure}[!h]
\centering
\includegraphics[width=5in]{prebotz5}
\caption{Raster plots for simulations of 150 preB\"{o}tzinger complex neurons connected in a 3D geometric network (middle column, all rows the same), as well as network realizations with the minimum (left column) and maximum (right column) Randi\'{c} index\xspace $R^{pq}$, where $pq = +\,+$, $+\,-$, $-\,+$, and $-\,-$ in rows 1, 2, 3 and 4, respectively.}
\label{fig:neuron5}
\end{figure}

\begin{figure}[!h]
\centering
\includegraphics[width=5in]{prebotz20}
\caption{Setup same as in Figure~\ref{fig:neuron5} for a different 3D geometric network.}
\label{fig:neuron20}
\end{figure}

\vfill\eject

\subsection{Normalizing the Randi\'{c} index\xspace index in cortical thickness networks}
\label{sec:ct}
Recently, researchers have begun using graph measures of connectivity to investigate the difference between structural magnetic resonance images taken from healthy individuals and individuals diagnosed with dementia (for example, \cite{Bai12,Yao10,HeChen09,HeChen08}). Networks in these studies are formed by calculating correlations between the cortical thicknesses of different brain regions, a technique based on the correlation between cortical thicknesss loss and dementia \cite{Fjell06,Gross12}. This method of network creation has been speculated to provide more insight on the functional relationships between brain regions \cite{Gong12}. Details of how the networks are formed can be found in \cite{HeChen08,SoldanMRJP14} and we present one such method below. After networks for each subject population (e.g., normal, subjects with dementia) are formed, network measures (e.g., Randi\'{c} index\xspace) are calculated and evaluated for significance.

A particular challenge to using comparisons between the network measures is complicated by the varying number of edges in each of the networks. In particular, some kind of normalization is required for some of the measures used.  In order to normalize, we propose the following scheme:
\begin{enumerate}
	\item For each network, $N$, determine the underlying degree distribution, $b$.
    \item Calculate $U_b$, the maximum Randi\'{c} index\xspace index for $b$, respectively.
    \item Use $U_b$ to normalize the Randi\'{c} index\xspace for $N$.
\end{enumerate}
We use data from the Alzheimer's Disease Neuroimaging Initiative (ADNI) and further analyzed by FreeSurfer, a technique developed by Fischl and Dale \cite{FischlD00} to measure cortical thicknesses. When applied to the ADNI data, Fischl and Dale discretized the cortical layer into 68 different regions. Soldan, et al.~\cite{SoldanMRJP14} used the following steps to generate networks:
\begin{enumerate}
\item Subset the population into categories based on whether they were diagnosed as normal (NORMAL), diagnosed with mild cognitive impairment for three or more years without disease progression (MCI), diagnosed with mild cognitive impairment and then progressed to Alzheimer's Disease within three years (MCIAD), or diagnosed with Alzheimer's Disease.
\item Within each population, use regression to control for subject age, gender, education level, and interaction effects between age and gender.
\item Use either partial or Pearson's correlations to calculate coefficients and $p$-values between each of the 68 regions.

\label{li:corr}
\item Use False Discover Rate \cite{Benjamini95} calculations to determine significant correlation coefficients with an error rate of 5\%.
\item Use one of the following schemes to determine edge weights: \label{li:weights}
	\begin{enumerate}
    \item Use no weights: simply include edges or not.
    \item Use the absolute value of the correlation coefficients as the edge weight. \label{li:abs}
    \item Use the product of the normalized cortical thicknesses (i.e., the cortical thicknesses divided by the maximimum)
    between the two regions connected by the significant edge. \label{li:prod}
    \item Use both \ref{li:abs} and \ref{li:prod}.
    \end{enumerate}
\end{enumerate}
Given the four diagnostic categories, the two possibilities in step \ref{li:corr}, and the four possibilities in step \ref{li:weights}, there are a total of 8 different networks. After calculating the associated degree distributions, we calculated the maximium Randi\'{c} index\xspace. As an example, we display the results for one network in Table \ref{tab:optimizedMRI}. The results for all eight networks can be seen in Appendix \ref{app:maxct}

\begin{table}
  \centering
  \begin{tabular}{l|c|c|c}
    \hline	Group	&	Original	&	Maximum ($U_b$)	& 	Normalized \\	\hline
AD	&	172910.5	&	182773.6	&	$94.6\%$	\\
MCI-AD	&	182098.8	&	192005.4	&	$94.84\%$	\\
MCI	&	146719	&	155715.5	&	$94.22\%$	\\
Normal	&	231957.8	&	244387.6	&	$94.91\%$	\\	\hline
  \end{tabular}
  \caption{Original, maximum, and normalized Randi\'{c} index\xspace. CT and absolute correlation weighting used with Pierson's correlations.}
  \label{tab:optimizedMRI}
\end{table}

By normalizing the Randi\'{c} index\xspace, comparisons between the groups can more accurately determine whether the assortativity was due to the actual network topology versus other features, such as the total number of edges. Currently, studies instead arbitrarily delete edges in networks, which effectively ignore different brain regions in order to compare non-normalized graph measures \cite{Yao10,HeChen09,HeChen08}. Normalizing by dividing out by the optimized metric allows for comparisons without ignoring network features. To actually use the optimization in such a statistical study, we note that significance testing (e.g., permutation testing to determine whether the comparisons are valid) would be required.

\section{Conclusions and Future Work}
We have shown that the minimum Randi\'{c} index\xspace problem can be solved in polynomial time.
With use of available $b$-matching code we have developed an algorithm that produces
a graph realization with the minimum Randi\'{c} index\xspace for a given degree sequence.  Not all optimal
solutions are connected.  A two-switch heuristic was developed to connect disconnected
optimal solutions.  Although the graph structure of these new connected graphs is fragile,
the Randi\'{c} index\xspace changed relatively little. \\

From our experiments, when generating Erd\H{o}s-R\'{e}nyi, geometric and scale-free graphs, the realization with
the minimum Randi\'{c} index\xspace is generally connected.  Undoubtably, this result is influenced by the parameters chosen
for the randomly generated graphs, but many graphs that will have a disconnected minimum Randi\'{c} index\xspace realization
have no connected realizations at all. \\

There are a number of future topics to explore.  We want to develop a
better way to connect graphs using the two-switch heuristic so that
the structure of the graph is less fragile.  Further experiments with
the parameters of randomly generated graphs are needed to understand
the conditions under which the number of graphs that are disconnected
or have no connected realizations changes.  We also are interested in determining the complexity
of the connected Randi\'{c} index\xspace problem when the input graph is not the complete
graph.

\bibliographystyle{plain}


\section{Background}
The low-temperature behaviour of the heavy-fermion metal YbRh$_2$Si$_2$\, was
studied by means of thermal and electrical transport across its
field-induced quantum critical point (QCP). The thermal
conductivity of metals is typically measured down to temperatures
in the Kelvin range.~\cite{2Smith08,2Wakeham11} Because of the low
energy scales, measurements in the heavy-fermion metals have been
extended down to temperatures as low as 0.06\,K
(Ref.~\onlinecite{2Tanatar07}) or even 0.04\,K
(Ref.~\onlinecite{2Tomokuni11}). In our case, special efforts have
been made to measure the thermal conductivity at $B = 0$ and
0.02\,T down to 0.025\,K due to the very low N\'eel temperature
$T_{\rm N}$ = 0.07\,K (at $B = 0$). Furthermore, we have performed
substantially more temperature scans for $B=0$ in order to reduce
the statistical error of the data. At $B = 0.06$\,T, which is
close to the critical field $B_{\rm c}$, we were able to measure
the thermal conductivity down to 0.04\,K. At higher fields, the
measurements of the thermal conductivity were performed down to
0.06\,K.
\section{Contacts}\label{sec:contact}
Previous thermal transport studies pointed out the strong
influence of thermal contact resistances and demonstrated the
necessity of a careful contact
preparation.\cite{2Tanatar07,2Seyfarth05} Subsequent to polishing
the YbRh$_2$Si$_2$\, sample (sample 1) for optimization of their geometry for
our measurements, the sample surface was cleaned in an ultrasonic
bath and rinsed in ethanol. The contact pads were prepared by
evaporating a gold film (thickness of $\approx 50$\,nm) and
applying a lift-off technique which allowed for optimum pad width
and separation. Gold evaporation was conducted by using electron
beam evaporators in ultra-high vacuum which turned out to provide
superior contacts compared to thermal evaporation. Gold wires
(50\,$\mu$m) were then attached to the contact pads by silver
paint to make use of the pad area. The silver paint covered an
area of width of about 130\,$\mu$m, which is very small compared
to the distance between the contacts (2.9\,mm). This contributes
to a systematic standard error of about 4.5\% in the estimation of
the geometry factor. It is the same for all experiments. We did
not observe any indication of diminished thermal conductivity due
to contact resistances down to 0.025\,K. The same contacts were
used for thermal and electrical conductivity measurements.
\section{Systematic errors}\label{sec:error}
Although most measurements were performed on a sample cut
specifically for thermal transport, sample 1 ($4.2 \times 0.5
\times 0.1$ mm$^{3}$), additional electrical resistivity
measurements were done on a second sample (sample 2) with average
size $1.7 \times 0.41 \times 0.06$ mm$^{3}$ from the same piece of
crystal. This was necessary to systematically study the dependency
of the low-$T$ electrical resistivity on the current, {\it cf.}
Section IV. Because of the inhomogeneous thickness of the samples,
the geometry factor could not be determined with very high
precision.
\begin{figure}[b]
\centering
\includegraphics[width=0.5\textwidth]{figS1.pdf}
\caption{Temperature dependency of the electrical resistivities of
two YbRh$_2$Si$_2$\,\,single crystals from the same batch.
The data of sample 2 have been multiplied by a
factor $1.25 \pm 0.03$ and corrected by a difference
in residual resistivity of 0.22\,$\mu\Omega$cm.}
\label{figS1}
\end{figure}
However, the measured resistivities could perfectly be rescaled by
a factor $1.25 \pm 0.03$ and corrected by a difference in residual
resistivity of 0.22\,$\mu\Omega$cm. The former number allows to
estimate the uncertainty in the determination of the sample
dimensions, {\it i.e.}, the systematic standard error for the
geometry factor. It is about 2.5\%, and it is the same for all
measurements. The thermal and electrical transport coefficients
were measured on the same sample (sample 1) with the same contact
configuration, and the Lorenz ratio $L(T)/L_0 =\rho(T)/w(T)$ is
affected by an additional systematic error of about 4.5\% due to
the finite width of the contacts, see Section~\ref{sec:contact}.
Thus, a total systematic error of 7\% has to be considered, {\it
i.e.}, it will shift systematically all curves of Fig.~3 of the
main text. This explains why in the region of the field-induced
Fermi-liquid in the phase diagram, {\it e.g.}, at $B \ge 0.6$\,T
and below $T$ = 0.15\,K, the difference $w(T) - \rho(T)$ is not
exactly zero, and the Lorenz ratio does not reach exactly one.
\section{Electrical resistivity}
The electrical resistivity $\rho(T)$ was determined by a
four-point ac-technique in a $^3$He-$^4$He dilution refrigerator
down to $T\approx$ 0.02\,K. Figure~S\ref{figS1} shows the
resistivity data of samples 1 and 2 at $B = 0$ and 0.1\,T,
respectively. The data for the second sample, which has a lower
residual resistivity, have been scaled as described above. The
data agree nicely, and we conclude that both samples show the same
overall behaviour. The slightly larger resistivity in the
\begin{figure}[t]
\centering
\includegraphics[width=0.5\textwidth]{figS2.pdf}
\caption{Temperature dependency of the electrical resistivity of
sample 2. The kink at about 0.07\,K indicates the N\'eel
temperature $T_{\rm N}$\,. Below $T_{\rm N}$\, a $T^2$ dependency of the resistivity
is observed as expected~\cite{2Gegenwart02} ({\it cf.} inset where
the individual curves are shifted vertically).}
\label{figS2}
\end{figure}
zero-field curve below 0.05\,K for sample 1 is due to heating
effects. Therefore, the second sample was measured at different
currents to investigate the influence of the current on the
curvature of $\rho(T)$. The results below 0.1\,K are displayed in
the inset of Fig.~S\ref{figS2} plotted as a function of $T^{2}$, a
temperature dependence expected below $T_{\rm N}$\, = 0.07\,K for $B = 0$
(Ref.~\onlinecite{2Gegenwart02}). The heating effect can be
neglected at a current of $5\,\mu$A. We used the resistivity data
with this current to evaluate the Lorenz ratio. As will be shown
in a forthcoming paper,~\cite{2Lausberg} the heating effect is
particularly strong in the vicinity of $B_c$. In addition, it is
increasing with increasing residual resistivity $\rho_0$. For a
single crystal of very high perfection with $\rho_0 \approx 0.5
~\mu\Omega$\,cm, no heating effect could be observed down to $T
\approx 0.02$\,K using a current of 50\,$\mu$A
(Ref.~\onlinecite{2Westerkamp08}). In this case, $\Delta \rho \sim
T$ was observed below $T= 0.1$\,K down to the lowest accessible
temperature of 0.02\,K.

In all samples the magnetoresistivity was measured and compared to
that of the samples investigated in
Refs.~\onlinecite{2Friedemann10,2Friedemann11}, to assure consistency
of the properties observed in our samples with those of the best
single crystals grown so far. As an example, the
\begin{figure}[t]
\centering
\includegraphics[width=0.5\textwidth]{figS3_new.pdf}
\caption{Magnetoresistivity of sample 2 measured with a current of
$5\: \mu$A. The upper inset magnifies the low field, low
temperature region of the main panel to emphasize that the
magnetoresistivity is positive inside the AF phase, but exhibits a
rapid drop across a crossover field $B^{*}$. This is further
illustrated in the lower inset for our lowest temperature $T =$
0.03 K. The rapid crossover across $B^*$ was systematically
analysed \cite{2Friedemann10,2Friedemann11}; in the zero-temperature
limit, it corresponds to a sharp jump because the crossover width
extrapolates to zero. The low-temperature data are also consistent
with the existence of a peak in the residual resistivity as a
secondary feature superimposed on the sharp jump, although its
interpretation is not clearcut as it may also be associated with
the N\'eel transition. }
\label{figS3}
\end{figure}
magnetoresistivity of sample 2 is shown in Fig.~S\ref{figS3}: At
$T = 0.03$\,K the change in resistivity
associated with the crossover is about 10\% (see insets of the
same figure), as also found in the samples studied in
Ref.~\onlinecite{2Friedemann10} ({\it cf.} Fig.~S5 of the
Supporting Information of Ref.~\onlinecite{2Friedemann10}). The
dominant feature is a rapid crossover across $B^*$, corresponding
to a jump in the zero-temperature limit
\cite{2Friedemann10,2Friedemann11}. The isothermal resistivity vs.
the magnetic field at low temperatures are also consistent with
the existence of a peak near $B^*$ as a secondary feature
superimposed on the sharp jump, although it could also be
associated with the (classical) N\'eel transition.

\section{Thermal conductivity}
\subsection{Experimental details}
The thermal conductivity $w(T)$ was measured using a steady-state
two-thermometer, one-heater technique in a $^3$He-$^4$He dilution
refrigerator down to $T\approx$ 0.025~K. The heat loss along the
wires and suspensions for thermometers and heater is estimated to
be a factor of 1000 smaller than the heat flow through the sample.
$\kappa$ was calculated from the electrical power $P$ released by
a resistive heater, the temperature difference $\Delta T$ between
two contacts on the sample and the geometry factor $A$ as $\kappa
= A \cdot P/ \Delta T$. $\Delta T$ was measured by two RuO$_{2}$
resistance thermometers, that were calibrated against the primary
(``cold finger``) sensor during each temperature run, using the
resistances measured at zero heating power. A stability of better
than 0.1\,\% of both sample temperature read-outs was achieved in
the entire measurement range. For each stabilized bath temperature
a set of four different heat currents was applied which results in
temperature gradients $\Delta T/T$ = 1\% - 7\% along the sample.
The proportionality between the applied heater power and the
achieved temperature gradients at a constant bath temperature
proves that the system is in the regime of linear response. The
resulting raw data of $\kappa$ were then averaged. Our measurement
procedure implies the following uncertainties: i) The uncertainty
in $A$ represents a systematic error, that shifts all $\kappa (T)$
curves by a constant factor, as described in
Section~\ref{sec:error}; ii) the uncertainty in $P$ is negligible
because current and voltage at the heater can be measured with
high accuracy; iii) the uncertainty in $\Delta T$ results from the
measurement of the thermometer resistances and from errors in the
calibration. Concerning the last point, two or more temperature
runs have been performed at many fields, each with its own
calibration. For a given field, the calculated $\kappa(T)$ values
fall on top of each other within the scattering of the data.
Therefore, the error in the calibration is negligible compared to
the data scattering. In fact, the largest uncertainty arises from
the measurement of the thermometer resistances due to the limited
excitation current, especially at low $T$. It leads to scattering
of the raw data and is significantly reduced in our final data by
averaging over several points. As error bars we took the standard
deviation of the raw data as displayed in Figs.~2 and 3 of the
main text as well as Figs.~S\ref{figS4} and S\ref{figS5}. The
systematic error due to the sample and contact geometry is not
included, because it results only in a shift of all curves as
explained in Section III. We attribute the observation of $w_0$
being about 1\% larger than $\rho_0$ in the Fermi liquid regime at
high fields (Figs.~2{\sffamily\bfseries d} and
S4{\sffamily\bfseries f}) to this systematic error.

Finally, the noise level in the data of the thermal resistivity
increases substantially with increasing magnetic field, as can be
directly seen in Fig.~S\ref{figS4}. It is presumably due to the
vibrations of the set-up wires in
\begin{figure}[t]
\centering
\includegraphics[width=0.47\textwidth]{supp_resistivity.pdf}
\caption{Thermal and electrical resistivity $w = L_0 T/\kappa_{\rm
el}\approx L_0 T/\kappa$ and $\rho$ plotted as a function of
temperature, at $T \le 0.5$\,K for $B = 0.02$\,T
{\sffamily\bfseries a}, 0.08\,T {\sffamily\bfseries b}, 0.1\,T
{\sffamily\bfseries c}, 0.2\,T {\sffamily\bfseries d}, 0.4\,T
{\sffamily\bfseries e} and 0.6\,T {\sffamily\bfseries f} $(B \perp
c)$. The arrows in {\sffamily\bfseries b} through
{\sffamily\bfseries f} indicate the crossover to Fermi-liquid
behaviour ({\it cf.} Fig.~1{\sffamily\bfseries a} of the main
text). A significant drop is seen in the residual ($T\rightarrow
0$) thermal and electrical resistivities when going from $B =
0.02$\,T {\sffamily\bfseries a} to 0.08\,T {\sffamily\bfseries b}
and 0.1\,T {\sffamily\bfseries c}, which parallels the drop in the
isothermal magnetoresistivity across $B_{\rm c}$ ({\it cf.}
Fig.~S\ref{figS3}) as previously observed and attributed to an
abrupt increase in the charge carrier
concentration.~\cite{2Friedemann10,2Friedemann11} Also, the slight
rise of the residual thermal and electrical resistivities upon
increasing field ({\it cf.} {\sffamily\bfseries d} through
{\sffamily\bfseries f}) confirms the trend observed in
Ref.~\onlinecite{2Friedemann11} and ascribed there to the
magnetoresistivity in the paramagnetic phase of YbRh$_2$Si$_2$.}
\label{figS4}
\end{figure}
magnetic fields. This is the main reason why we could not perform
reliable measurements below 0.06\,K at fields larger than 0.06\,T.
\subsection{Phonon contribution}
The phonon contribution, $\kappa_{\rm ph}(T)$\,, to the measured thermal
conductivity $\kappa(T)$ can be separated from the electronic
part, $\kappa_{\rm el}(T)$\,, by avoiding the low-temperature range where a
significant inelastic scattering of the charge carriers has been
evidenced ({\it cf.} Fig.~1{\sffamily\bfseries b} of the main text
and Ref.~\onlinecite{2Tomokuni11}).

After having subtracted $\kappa_{\rm el}(T)$\, = $\kappa_{WF}(T) = L_0 T/\rho(T)$
from $\kappa(T)$ for 6\,K $< T < 12$\,K, {\it i.e.}, assuming the
validity of the Wiedemann-Franz (WF) law to hold in this
temperature range for the electronic heat transport, the phonon
contribution $\kappa_{\rm ph}(T)$\, is found to follow a $T^{\epsilon}$ dependence
with $\epsilon = 2 \pm 0.2$. As shown in Fig.~1{\sffamily\bfseries
b}, main text, an extrapolation of this power law to lower
temperatures indicates a negligible $\kappa_{\rm ph}(T)$\, below 1\,K, {\it i.e.},
within the temperature range of interest in the present work. The
uncertainty in the exponent has only little effect on the
estimated $\kappa_{\rm ph}(T)$\, below 5\,K. The insignificance of $\kappa_{\rm ph}(T)$\, below 1\,K
is further corroborated by the almost constant Lorenz ratio $L /
L_0 \approx 0.86$ within
\begin{figure}[t]
\centering
\includegraphics[width=0.47\textwidth]{supp_analysis.pdf}
\caption{Difference $(w - \rho)$ vs $T$ at $B = 0.02$\,T, 0.08\,T
and 0.1\,T {\sffamily\bfseries a} as well as 0.2\,T, 0.4\,T and
0.6\,T {\sffamily\bfseries b}. {\sffamily\bfseries c},
{\sffamily\bfseries d} Lorenz ratio $L/L_0 = \rho/w$ vs $T$ for
the same fields as in {\sffamily\bfseries a}, {\sffamily\bfseries
b}. The data for $B = 0.02$\,T are similar to those for $B = 0$
(\textit{cf.} Figs.~3{\sffamily\bfseries a} and
{\sffamily\bfseries c} of the main text). They provide convincing
evidence for an extra heat channel which adds to the one of the
electronic quasiparticles and ist most likely due to AF magnons
(see text). For $B = 0.6$\,T a Fermi-liquid phase forms below $T
\approx 0.15$\,K. This can also be anticipated for $B = 0.4$\,T
and 0.2\,T and even for $B = 0.1$\,T and 0.08\,T at
correspondingly lower crossover temperatures.} \label{figS5}
\end{figure}
the range $0.5$\,K $< T  < 1$\,K, presumably associated with a
measured thermal conductivity of purely electronic origin.
Assuming this Lorenz ratio to extend to above 2\,K, $\kappa_{\rm el}(T)$\, = $L\,
T/\rho(T)$ can easily be estimated. Subtracting this from the
measured $\kappa(T)$, the phonon part is obtained once more. Below
about 5\,K, $\kappa_{\rm ph}(T)$\, obtained by the latter procedure is found to be
in very good agreement with that from the former one.

We, therefore, conclude that  below about 10\,K, $\kappa_{\rm ph}(T)$\, of YbRh$_2$Si$_2$\,
follows a $T^{2}$ dependence, as commonly expected for a metal in
this temperature range due to dominating phonon scattering from
the conduction electrons.~\cite{2Ziman60} A nearly
$T^{2}$-dependence of $\kappa_{\rm ph}(T)$\, has indeed been observed for several
heavy-fermion compounds, e.g., CeB$_{6}$,~\cite{2Peysson85}
CeNiSn~\cite{2Kitagawa02} and CeAuAl$_{3}$~\cite{2Aoki00} in a
similar temperature range. In these systems the dominant phonon
wavelength, $\lambda_{\rm ph}$, is assumed to be shorter than the
mean free path of charge carriers, $l_{\rm el}$,
Ref.~\onlinecite{2Zimmerman59}. This assumption, thus, appears to
hold also for the YbRh$_2$Si$_2$\, single crystal studied here. On the other
hand, for CeCu$_{2}$Si$_{2}$,~\cite{2Franz78}
CeAl$_{3}$~\cite{2Ott84} and CeCu$_{6}$~\cite{2Peysson86} $\kappa_{\rm ph}(T)$\,
depends almost linearly on $T$, which hints at the opposite
relation between $\lambda_{\rm ph}$ and $l_{\rm el}$, {\it i.e.},
$\lambda_{\rm ph} > l_{\rm el}$, Ref.~\onlinecite{2Zimmerman59}.
\subsection{Electronic contribution}
The thermal conductivity obtained by subtracting $\kappa_{\rm ph}(T)$\, from the
measured $\kappa(T)$ is $\kappa_{\rm el}(T)$. In Fig.~S\ref{figS6}
we show that the temperature coefficient, $\kappa_{\rm el}/T$,
increases logarithmically upon cooling from $T = 2$\,K to about
\begin{figure}[t]
\centering
\includegraphics[width=8.4cm]{supp_kappa-vs-logT.pdf}
\caption{Temperature coefficient of the $B = 0$ electronic thermal
conductivity, $\kappa_{\rm el}/T$, of YbRh$_2$Si$_2$\, plotted as a function
of $T$ on a logarithmic scale to emphasize that it follows a
logarithmic behaviour on cooling from 2\,K to about 0.3\,K.}
\label{figS6}
\end{figure}
0.3\,K. A corresponding logarithmic divergence was observed also
for the Sommerfeld coefficient of the electronic specific heat,
$\gamma = C_{\rm el}/T$,~\cite{2Trovarelli00,2Custers03} and the
thermopower coefficient $S/T$.~\cite{2Hartmann10}

In the following, we are interested in the corresponding thermal
resistivity $w(T) = L_0 T/\kappa_{\rm el} (T)$, which is displayed
together with $\rho(T)$, in Fig.~2 of the main text for four
different fields. The results of $w(T)$ and $\rho(T)$ for other
magnetic fields are shown in Fig.~S\ref{figS4}. The data taken at
$B = 0.08$\,T through $0.6$\,T are complementary to those at $B =
0.3$\,T and 1 T presented in Figs.~2{\sffamily\bfseries c} and
2{\sffamily\bfseries d} of the main text. The data at $B =
0.02$\,T are likewise complementary to those at zero field
presented in Fig.~2{\sffamily\bfseries a} of the main text. Like
in the zero-field case, they show a low-temperature downturn in
the AF phase which contains the magnon contribution to the heat
conduction. Figure S\ref{figS5} displays the difference $w(T) -
\rho(T)$ {\sffamily\bfseries a},{\sffamily\bfseries b} and the
ratio $\rho(T)/w(T) = L/L_0$ {\sffamily\bfseries
c},{\sffamily\bfseries d} for the data shown in Fig.~S\ref{figS4}.
The existence of a magnon contribution to the thermal conductivity
at $B = 0.02$\,T is indicated by $(w - \rho)$ being negative and
$L/L_0$ being larger than one below $T \approx 0.03$\,K, very
similar to the results at $B = 0$ (\textit{cf.}
Figs.~3{\sffamily\bfseries a} and 3{\sffamily\bfseries c} of the
main text). Within the experimental uncertainties, the data at $B
= 0.6$\,T indicate Fermi-liquid behaviour and the validity of the
WF law below $T \approx 0.15$\,K. Fermi-liquid behaviour below a
crossover temperature which continuously decreases with decreasing
magnetic field is inferred for both $B = 0.4$\,T and 0.2\,T
{\sffamily\bfseries b},{\sffamily\bfseries d}. This implies
$L/$L$_{0} \rightarrow 1$ as $T \rightarrow 0$, with which also
the results obtained for $B = 0.1$\,T and even 0.08\,T are
compatible {\sffamily\bfseries a},{\sffamily\bfseries c}.

The data for $w$ and $\rho$ measured as a function of temperature
for fixed magnetic fields, shown in Fig.~2 of the main text and
Fig.~S\ref{figS4}, have been used to determine the isothermal
Lorenz ratio as a function of field displayed in
Fig.~3{\sffamily\bfseries e} of the main text. The isothermal
field dependence exhibits a minimum near $B^*$. This striking
behaviour is also seen in preliminary results obtained from direct
isothermal measurements of $w$ and $\rho$ as a function of field
in a sample with lower residual resistivity~\cite{2isothermal_new}.

\subsection{Magnon contribution}
For the measurements performed at $B = 0$
(Fig.~2{\sffamily\bfseries a}, main text) as well as 0.02\,T
(Fig.~S\ref{figS4}{\sffamily\bfseries a}) a downturn in the
thermal resistivity is found below $T \approx 0.1$\,K and $T
\approx 0.08$\,K, respectively, while at $B = 0.06$\,T a similar
but smaller feature occurs below $T \approx 0.07$\,K
(Fig.~2{\sffamily\bfseries b}, main text). One may ascribe this
drop in $w(T)$ to the freezing out of inelastic scatterings
provided by spin fluctuations in the electronic heat transport.
However, while ferromagnetic spin fluctuations remain unchanged,
antiferromagnetic ones grow at sufficiently low temperatures and
magnetic fields, as inferred from NMR Knight shift and
spin-lattice relaxation rate measurements, respectively.
\cite{2Ishida.02} Therefore, the only natural interpretation of
this drop involves heat carriers which add to the electronic ones.
The existence of a corresponding excess contribution to the
thermal conductivity, $\Delta \kappa (T)$, is proven by our
observation that for $B=0$ and $B=0.02$ T the thermal resistivity
becomes smaller than the electrical resistivity, {\it cf.}
Figs.~2{\sffamily\bfseries a} and S4{\sffamily\bfseries a}.
$\Delta \kappa (T)$ adds to $\kappa_{\rm el} (T)$ which we assume
to be given by $L_0 T/\rho$ at temperatures well below $T_{\rm N}$\,, {\it
i.e.}, we assume the Wiedemann-Franz law to describe the
electronic transport in the Fermi-liquid phase well below $T_{\rm
N}$. The observed $\Delta\kappa(T)$ falls into the range $2 -
5\cdot10^{-3}$\,W/Km between 0.025 and 0.03\,K. Because this value
of $\Delta\kappa$ is close to the experimental uncertainty, its
temperature dependence is hard to be experimentally determined.
Nevertheless, its existence has been confirmed by repeated
measurements. Lattice vibrations are unapt to account for this
extra thermal conductivity, as the largest expected $\kappa_{\rm ph}$\,, which is
limited by the sample dimension ($\approx 100\, \mu$m) and,
employing the lattice specific heat,~\cite{2Custers03} is estimated
to be less than $1\cdot 10^{-5}$\,W/Km in this range.

Another potential thermal heat channel is that of AF magnons. The
signature of these spin-wave excitations was clearly observed
\begin{figure}[t]
\centering
\includegraphics[width=0.52\textwidth]{figS7.pdf}
\caption{Specific heat of YbRh$_2$Si$_2$\,, shown as $\Delta C /T$ versus
$T^{2}$. $\Delta C(T) = C(T) - C_{\rm ph}(T) - C_{\rm Q} (T)$,
where $C_{\rm ph}$/$C_{\rm Q}$ denotes the phonon/nuclear
quadrupole contribution. The red line indicates a $T^{3}$
contribution to $\Delta C$ below 0.05\,K.~\cite{2Custers03}}
\label{figS7}
\end{figure}
\cite{2Custers03} in the zero-field specific heat examined down to
$T = 0.018$\,K for an YbRh$_2$Si$_2$\, single crystal being of similar quality
as samples 1 and 2. As seen in Fig.~S\ref{figS7}, below $T =
0.05$\,K the specific heat can be described by $\Delta C = C -
C_{\rm ph} - C_{\rm Q} = \gamma T + \beta T^{3}$, with $C_{\rm
ph}$ and $C_{\rm Q}$ being, respectively, the phonon and the
nuclear quadrupole contributions, $\gamma = 1.64$\,J/(K$^{2}$mol)
and $\beta = 132.2$\,J/(K$^{4}$mol). The huge electronic
contribution $C_{\rm el} = \gamma T$ denotes a heavy Landau FL
phase.~\cite{2Custers03} The Debye-like term $C_m \sim T^{3}$ is
characteristic of the contribution of long-wavelength AF acoustic
magnons. In the framework of the Debye theory and using the
measured $\beta$ value, we can estimate the ``magnetic Debye
temperature", $\Theta_m$, and the group velocity $v_{m}$ of the
magnons to be 4.2\,K and 36\,m/s, respectively. This Debye
temperature corresponds to the AF exchange interaction suitably
averaged over the three main spatial directions. The extracted
group velocity is substantially smaller than that for typical
Ce-based heavy-fermion AF metals, {\it e.g.}, about one-tenth of
the corresponding value for CeAl$_{2}$,~\cite{2Bredl78} reflecting
both the averaging over the spatial directions as well as the weak
N\'eel order of YbRh$_2$Si$_2$\,. The classical kinetic relation, $\kappa_{m}
= (1/3) C_{m} v_{m} l_{m}$ ($l_{m}$: magnon mean free path),
allows to estimate the magnon contribution to the thermal
conductivity, $\kappa_{m}$. In order to obtain a $\kappa_{m}$ of
$2 - 5\cdot10^{-3}$\,W/Km in the range $0.025 - 0.03$\,K as
observed in Fig.~2{\sffamily\bfseries a} of the main text and
Fig.~S\ref{figS4}{\sffamily\bfseries a}, the corresponding magnon
mean free path $l_{m}$ has to be in the range of $2 - 13$\,$\mu$m.
Since in the low-temperature limit $l_{m}$ is expected to be equal
in size to that of the AF domains, the latter are estimated to be
of the order of a few $\mu$m below about 0.03\,K which, indeed, is
a reasonable order of magnitude.~\cite{2Note}

The less pronounced downturn in $w(T)$ observed below $T =
0.07$\,K at $B = 0.06$\,T (Fig.~2{\sffamily\bfseries b}, main
text) is ascribed to overdamped AF magnons, which were shown, via
inelastic neutron scattering
experiments,\cite{2Wiltshire83,2Wiltshire85} to exist for
antiferromagnetically ordered materials substantially above the
N\'eel temperature. Heat transport by short-lived magnon
excitations has been reported, e.g., for the parent compound of
the 214 high-$T_{c}$ cuprates, La$_{2}$CuO$_{4}$, an $S = 1/2$, 2D
antiferromagnet with a N\'eel temperature $T_{\rm N} \approx
310$\,K.\cite{2Hess03} In particular, $S = 1$ chain systems, like
Y$_{2}$BaNiO$_{5}$ (Ref.~\onlinecite{2Kordonis05}) and
Ni$($C$_{2}$H$_{8}$N$_{2})_{2}$NO$_{2}($ClO$_{4})$ (NENP)
(Ref.~\onlinecite{2Sologubenko08}) have served as model systems in
this context.

We would like to note that the magnon contribution cannot be
avoided by performing thermal Hall measurements. In contrast to
the case of phonons, \cite{2Wakeham11} magnons may generate a
transverse thermal gradient. \cite{2Onose10}

\section{Extrapolation of the Lorenz ratio to zero temperature in
quantum critical systems}

We now discuss the various isofield and isothermal scans that are
used for the extrapolation of the Lorenz ratio in the vicinity o a
QCP. We will discuss the special case that is pertinent to YbRh$_2$Si$_2$,
namely a field-induced QCP separating an AF ordered phase and a
paramagnetic Fermi-liquid phase.

Fig.~S\ref{figS8}{\sffamily\bfseries a} illustrates different
temperature scans in different parts of the phase diagram. Scan A
starts from the quantum critical regime, but runs into the ordered
phase; by construction, it passes through the phase boundary. The
zero-temperature limit of $w(T)$, $\rho(T)$ and $L(T)/{L_0}$ so
extrapolated does not reflect the quantum critical behavior, but
instead only captures that of the ordered phase. Scan B also
starts from the quantum critical regime, but runs into the
low-temperature paramagnetic Fermi liquid phase; it passes through
a crossover temperature. Again, the zero-temperature limit so
extrapolated does not reflect the quantum critical behavior, but
instead only captures that of the paramagnetic Fermi liquid phase.
Scan C starts from the quantum critical regime, and goes all the
way to the QCP as the temperature is lowered. The zero-temperature
limit so extrapolated captures the properties of the QCP. In
practice, however, the same purpose can be achieved by carrying
out scans C$_1$ or C$_2$ (provided there is enough dynamical range
in temperature). These start in the quantum critical regime but
stop before running into the N{\' e}el temperature or the
crossover temperature to the paramagnetic Fermi liquid region. The
zero-temperature limit so extrapolated also manifests the behavior
of the QCP.

The different temperature scans will also be manifested in the
isothermal properties of $L/{L_0}$ as a function of the control
parameter, as illustrated in Fig.~S\ref{figS8}{\sffamily\bfseries
b}. Scan a starts from the paramagnetic Fermi-liquid regime, and
runs into the quantum critical regime. As it passes through the
crossover scale between the two regimes, we expect to see a
crossover in $L/{L_0}$ in a way that compliments what is seen in
the temperature scan B. Scan B starts from the paramagnetic
Fermi-liquid regime, passing through the quantum critical regime,
and runs into the ordered phase. Because it passes through both
the crossover temperature and the N\'eel temperature lines, this
isothermal scan will
\begin{figure}[t]
\centering
\includegraphics[width=0.45\textwidth]{figS8_final}
\caption{ Isofield {\sffamily\bfseries a} and isothermal
{\sffamily\bfseries b} scans for the extrapolation of the Lorenz
ratio in the vicinity of a QCP, as described in detail in the
text. } \label{figS8}
\end{figure}
reflect the features of both temperature scan B and temperature
scan A.

In this work, our primary focus is on temperature scans C$_1$
(Figs.\ 2{\sffamily\bfseries a}, 3{\sffamily\bfseries a},
3{\sffamily\bfseries c}), scans close to C (Fig.
2{\sffamily\bfseries b}), and scans C$_2$ ({\it e.g.} Figs.
S4{\sffamily\bfseries b} and {\sffamily\bfseries c}) through B
(Figs. 2{\sffamily\bfseries c}, 2{\sffamily\bfseries d},
3{\sffamily\bfseries b}, 3{\sffamily\bfseries d}). In addition,
isothermal scans a (Fig. 3{\sffamily\bfseries e}) are conducted.
Because of the masking effect of magnons on the electronic heat
transport, we do not attempt to analyze any scan b. Temperature
scans A cannot be performed in YbRh$_2$Si$_2$ because of the small
value of the critical field $B_c$.

\section{Theoretical aspects}
The validity of the WF law in the vicinity of a heavy-fermion
spin-density-wave QCP was briefly discussed in the main text.
Across such a QCP, the heavy quasiparticles remain intact in the
main part of the Fermi surface. In these ``cold'' regions, the
quasiparticles do not experience scattering by the AF spin
fluctuations, and their spectral weight---the quasiparticle
residue $Z$---stays finite as the control parameter moves through
the QCP (Fig.~S9). $L/{L_0}$ must be equal to one in this case.
Even the contributions from the ``hot'' regions, which experience
scatterings by the bosonic collective fluctuations, cannot violate
the WF law, as was illustrated by the case of ZrZn$_2$
(Ref.~\onlinecite{2Smith08}) where ferromagnetic fluctuations
influence the entire Fermi surface.

This is in contrast to the Kondo-destroying local quantum critical
description which is illustrated in Fig.~4 of the main text. In
the paramagnetic Fermi-liquid state, the conduction electron
\begin{figure}[t]
\centering
\includegraphics[width=0.45\textwidth]{figS9_final.pdf}
\caption{The evolution of the quasiparticle weights across a
heavy-fermion spin-density-wave quantum critical points. $Z_L$ of
the large Fermi surface stays nonzero as the system is tuned
through the QCP.} \label{figS9}
\end{figure}
self-energy, $\Sigma ({\boldsymbol k},{\omega},T)$ contains a pole
at the Kondo-resonance energy. It converts the $f$-moments into a
part of the Landau quasiparticles, thereby creating a large Fermi
surface. The quasiparticle residue at the large Fermi momenta,
$\pmb{k}^{\rm L}_{\rm F}$, is nonzero, $Z_L(\pmb{k}^{\rm L}_{\rm
F},\omega=0,T=0) \ne 0$. In the AF Fermi-liquid state, the Kondo
resonance is destroyed. The Fermi surface is given by that of the
conduction electrons only in the presence of a staggered magnetic
field. Away from the ``hot spots", there are Landau quasiparticles
associated with a small Fermi surface. The quasiparticle residue
at such generic small Fermi momenta, $\pmb{k}^{\rm S}_{\rm F}$, is
nonzero, $Z_S(\pmb{k}^{\rm S}_{\rm F},\omega=0,T=0) \ne 0$.

In these Fermi-liquid regimes, both the electrical and electronic
heat currents are predominantly carried by Landau quasiparticles.
At nonzero temperatures, the electronic heat carriers experience
inelastic scatterings. Because the non-Umklapp processes are
considerably more efficient in relaxing electronic heat current
than electrical current, even in anisotropic systems at low
temperatures, the electron-electron scatterings lead to an
electronic Lorenz ratio $L_{\rm el}/{L_0} < 1$ at nonzero
temperatures.~\cite{2Herring67} The effect is similar to the usual
case of electron-phonon scattering.~\cite{2Ziman60} The Umklapp
scatterings are expected to contribute equally efficiently to the
thermal and electrical resistivities. In the zero-temperature
limit, well-defined quasiparticles remain, but the inelastic
scatterers are frozen out; only the elastic scattering processes
remain, and $L_{\rm el}/{L_0}$ is equal to $1$. This is
illustrated in Fig.~S\ref{figS10}, where $L_{\rm el}/{L_0}$
reaches $1$ as the system moves away from the QCP (at $B_c$) into
the Fermi-liquid regimes on both sides.

In the quantum critical regime, low-energy electronic excitations
occur at both the small and large Fermi surfaces (Fig.~4, main
text). The single-electron self-energy at both $\pmb{k}^{\rm
S}_{\rm F}$ and $\pmb{k}^{\rm L}_{\rm F}$ vanishes at the
$\omega=0$ and $T=0$ limit, but has the scaling form as a function
of $\omega$ and $T$ as given in Eq.~1 (main text). This reflects
the fluctuations of the Fermi surfaces, which
\begin{figure}[b]
\centering
\includegraphics[width=0.45\textwidth]{figS10_final.pdf}
\caption{Expected isothermal behaviour of the Lorenz ratio across
a Kondo-destroying QCP. $B^*(T) = B(T^*)$ is the Kondo breakdown
scale at a nonzero but small $T$. At zero temperature an abrupt
dip occurs at $B^* = B_c$} \label{figS10}
\end{figure}
characterize the quantum fluctuations in the {\it entire} quantum
critical regime. We will consider the electrical and heat currents
carried by the electronic excitations at both $\pmb{k}^{\rm
S}_{\rm F}$ and $\pmb{k}^{\rm L}_{\rm F}$, all of which are
non-Fermi liquid in nature as specified by Eq.~1. These electronic
current carriers are subject to inelastic scatterings that are
associated with the quantum criticality. While such scattering
processes are many-body in nature, they can still be divided into
Umklapp and non-Umklapp processes. The non-Umklapp processes will
contribute considerably more to the thermal resistivity than to
the electrical resistivity. As a result, $L_{\rm el}/{L_0}$ will
be less than~1.

The end result is a dip of $L_{\rm el}/{L_0}$ near $B_c$.
Fig.~S\ref{figS10} illustrates the corresponding isothermal
behaviour of $L_{\rm el}/{L_0}$ for a given nonzero but low $T$,
where the dip is expected to be centered at $B^*(T) = B(T^*)$,
{\it cf.} Fig~1{\sffamily\bfseries a} of the main text. As $T
\rightarrow 0$, the dip becomes abrupt at $B_c$.