\section{Acknowledgments}

The first author was supported by a Graduate Fellowship from Rensselaer Polytechnic Institute during the completion of this work.
\section{Alternative Models}

\subsection{Steiner Tree Model}

In addition to studying OCNs, we took a preliminary look at comparing OCNs with Steiner trees.

Steiner trees are minimum-weight trees that connect a given subset of nodes in a graph.  The specified nodes are called terminals.  Steiner trees may include non-terminal nodes; these are called Steiner points (or Steiner nodes).

In \cite{Chung1989Steiner}, Steiner trees built on a grid are considered.  The underlying graph is all points in a Euclidean plane, completely connected, and the terminals are the points of a grid.  We began comparing Steiner trees with OCNs by looking at the relative optimality of OCNs minimizing $H_{0.5}(s)$ and normal Steiner trees on grids of the same sizes.

In order to evaluate the energy of a Steiner tree, the tree needs to be directed; this can be done by choosing the outlet to be an arbitrary corner node and directing all the edges towards the outlet.

For the purposes of this study area was defined as the number of terminal nodes in the subtree rooted at a given node.  This way, the terminal nodes have exactly one unit of area associated with each of them, just as the grid nodes do in the OCN.

The OCNs here have realistic link lengths, where diagonal links are length $\sqrt{2}$, to make them comparable to the Steiner trees.

OCNs $s$ of size 2x2, 3x3, 4x4, and 5x5 were constructed and $H_{\gamma}(s)$ for each were compared with $H_{\gamma}(t)$ for Steiner trees $t$ of the same size.  The Steiner trees used were found in \cite{Chung1989Steiner}.  Figure \ref{fig:steinerpics} depicts the pairs of trees that were compared.  On the left of each pair is the Steiner tree, with open circles for the Steiner nodes, and on the right of each is the corresponding OCN.

\begin{figure}[H]
	\begin{subfigure}[b]{0.17\textwidth}
		\centering
		\includegraphics[scale = .5]{images/Steiner/2x2pics.png}
		\caption{2x2}
		\label{fig:2x2pics}
	\end{subfigure}
	\begin{subfigure}[b]{0.21\textwidth}
		\centering
		\includegraphics[scale = .5]{images/Steiner/3x3pics.png}
		\caption{3x3}
		\label{fig:3x3pics}
	\end{subfigure}
	\begin{subfigure}[b]{0.27\textwidth}
		\centering
		\includegraphics[scale = .5]{images/Steiner/4x4pics.png}
		\caption{4x4}
		\label{fig:4x4pics}
	\end{subfigure}
	\begin{subfigure}[b]{0.33\textwidth}
		\centering
		\includegraphics[scale = .5]{images/Steiner/5x5pics.png}
		\caption{5x5}
		\label{fig:5x5pics}
	\end{subfigure}
	\caption{Steiner trees and OCNs for different nxn grids.}
	\label{fig:steinerpics}
\end{figure}

Since Steiner trees minimize total link length, it was clear that they would always have lowest energy for $\gamma = 0$.  For larger values of gamma, however, where area becomes a bigger factor than link length, we hypothesized that the OCNs would have a lower energy value than the Steiner trees.  Finding whether this was true, and where the crossover point would be, was the goal of the comparisons.

\subsection{Analytic Results}

\begin{figure}[H]
	\begin{subfigure}[b]{0.49\textwidth}
		\centering
		\includegraphics[scale = .75]{images/Steiner/2x2graph.pdf}
		\caption{2x2}
		\label{fig:2x2graph}
	\end{subfigure}\quad
	\begin{subfigure}[b]{0.49\textwidth}
		\centering
		\includegraphics[scale = .75]{images/Steiner/3x3graph.pdf}
		\caption{3x3}
		\label{fig:3x3graph}
	\end{subfigure}\\
	\begin{subfigure}[b]{0.49\textwidth}
		\centering
		\includegraphics[scale = .75]{images/Steiner/4x4graph.pdf}
		\caption{4x4}
		\label{fig:4x4graph}
	\end{subfigure}\quad
	\begin{subfigure}[b]{0.49\textwidth}
		\centering
		\includegraphics[scale = .75]{images/Steiner/5x5graph.pdf}
		\caption{5x5}
		\label{fig:5x5graph}
	\end{subfigure}
	\caption{Plots of $H_{\gamma}(s)$ (blue) and $H_{\gamma}(t)$ (red), where $s$ is an OCN and $t$ a Steiner tree of the given size.  In each plot, the horizontal axis is $\gamma$ and the vertical axis is energy.}
	\label{fig:steinergraphs}
\end{figure}

Figure \ref{fig:steinergraphs} shows plots of the energy for corresponding Steiner trees and OCNs, $0 \le \gamma \le 1$.  As expected, the Steiner trees have lower energy than the OCNs for smaller $\gamma$, while for $\gamma$ past a certain crossover point the OCNs are better.  Table 3
lists the crossover values of $\gamma$ for each size.

\begin{table}[H]
	\centering
	\begin{tabular}{| l l |}
		\hline
		size & crossover $\gamma$\\ \hline
		2x2 & 0.55433153188825812933\\
		3x3 & 0.63036154802337965422\\
		4x4 & 0.01907039679465522977\\
		5x5 & 0.32264708810531986705\\
		\hline
	\end{tabular}
	\caption{Crossover values of $\gamma$}
	\label{tab:crossovers}
\end{table}
\medskip

\subsection{Interpretation}

From the sizes studied, it is not possible to detect a trend in crossover points.

The next step would be to combine the two models, constructing Steiner-type trees that minimize $H_{\gamma}(s)$ for varying $\gamma$ and comparing them with real river networks and with normal OCNs.
\section{Analytical Proofs} \label{sec:proofs}

\subsection{Energy bounds}

In \cite{Colaiori1997Analytical}, Colaiori et al. proved lower bounds for the energy of OCNs on an orthogonal grid (one with no diagonal links).  Here, we look at lower bounds for the energy of OCNs on an eight-neighbor grid.  The main result for this section (Theorem 3) is that for an $n \times n$ 8-way grid, $\frac{3}{2}n^2 - \frac{7}{2}n + 1$ is a lower bound on $H_{0.5}(s)$.

\begin{lem} \label{lem3} Let (P) be the optimization problem

	\[ \begin{array}{rl}
		\min        & \displaystyle\sum_{i = 1}^n A_i^{\gamma} \\
		\mbox{s.t.} & \displaystyle\sum_{i = 1}^n A_i = n + m \\
					& A_i \ge 1 \:\: \forall i \\
					& A_i \in \mathbb{Z} \:\: \forall i
		\end{array}
	\]

	where $n$ is the number of entries in the vector $A$, $m$ is some natural number, and $0 \le \gamma \le 1$. \\

	$A_1 = A_2 = ... = A_{n-1} = 1, A_n = m+1$ is an optimal solution to $(P)$.
\end{lem}

\begin{proof}
Suppose A is an optimal solution for $(P)$ and $\exists j$ s. t. $1 < A_j < m + 1$.

$\exists k \:\: s. t. \:\: 1 < A_k < m + 1, j \ne k$.

We can construct $A^*$ such that $A^*_i = A_i \:\: \forall i \ne j, k$; $A^*_j = 1$; $A^*_k = A_j + A_k - 1$.

Then $\sum_{i = 1}^n {A^*_i}^{\gamma} = (\sum_{i = 1}^n A_i^{\gamma}) - A_j^{\gamma} - A_k^{\gamma} + 1 + (A_j + A_k - 1)^{\gamma}$.

Since $f(x)= x^\gamma$ when $0 \le \gamma \le 1$ is concave, and $A_j, A_k > 1$, $(A_j + A_k - 1)^{\gamma} \le A_j^{\gamma} + A_k^{\gamma} - 1$.  So $\sum_{i = 1}^n {A^*_i}^{\gamma} \le \sum_{i = 1}^n A_i^{\gamma}$.  Thus $A^*$ must be an optimal solution for $(P)$.  By an iteration of this argument, there must be an optimal solution $A^+$ with all entries except for one equal to $1$.  Without loss of generality, the entries can be rearranged so that $A^+_n$ is the largest (only non-unit) entry.
\end{proof} \medskip

\begin{thm} For an $n \times n$ eight-neighbor grid $G$,  the optimal $H_1(s)$ is $\frac{4n^3 - 3n^2 - n}{6}$. \end{thm}
\begin{proof}

\begin{figure}[h]
	\centering
	\includegraphics[scale = .5]{images/5x5stripes.png}
	\caption{The five ``stripes'' of a $5 \times 5$ 8-way grid with the outlet in the bottom left corner.}
	\label{fig:stripes}
\end{figure}

Consider the subsets of vertices $C_k$ where vertex $i \in C_k$ iff the distance from $i$ to the outlet equals $k-1$ (Figure \ref{fig:stripes}).  In an $n \times n$ 8-way grid there are n such stripes.  Assuming the outlet is located in a corner, $|C_k| = 2k - 1$.
Clearly,

\begin{equation} \label{eq:doublesum}
H_\gamma(s) = \sum_{i \in links(s)} A_i^\gamma = \sum_{k=2}^n \sum_{i \in C_k} A_i^\gamma
\end{equation}

Recall that the area of node $i$ is the number of nodes upstream of $i$, plus $1$ for $i$ itself.  So for a given $k$, $\sum_{i \in C_k} A_i^\gamma$ will have to include the number of nodes in $C_k$ and the total number of nodes upstream of all $i \in C_k$, which must include all the nodes $j \in \bigcup_{l > k} C_l$, since their paths to the outlet must pass through stripe $C_k$ at some point.  $\bigcup_{l > k} C_l$ is simply the set of all the nodes in the grid minus those within the $k$x$k$ section enclosed by $C_k$, which is $n^2 - k^2$ nodes.  Thus the total amount of area being passed through $C_k$ is at least
\begin{equation} \label{qu:stripearea}
|C_k| + n^2 - k^2
\end{equation}

Then we have

\[H_1(s) \: =\: \sum_{k=2}^n \sum_{i \in C_k} A_i \:\ge\: \sum_{k=2}^n ( 2k - 1 +  n^2 - k^2 ) \]

The rightmost sum simplifies to $\frac{1}{6}(4n^3 - 3n^2 - n)$; therefore this is a lower bound on $H_1(s)$.

Consider the network $t$ where every link from node $i$ to node $j$ is such that if $i \in C_k$ then $j \in C_{k-1}$.  Here, quantity \ref{qu:stripearea} is exactly the area passing through stripe $C_k$, so $H_1(t) = \frac{1}{6}(4n^3 - 3n^2 - n)$.  Hence the lower bound is achievable.

\end{proof} \medskip

\begin{thm} \label{thm:main} For an $n \times n$ eight-neighbor grid, $\frac{3}{2}n^2 - \frac{7}{2}n + 1$ is a lower bound on $H_{0.5}(s)$. \end{thm}
\begin{proof}

Following the same reasoning as above, we have that

\begin{equation}
H_\gamma(s) = \sum_{k=2}^n \sum_{i \in C_k} A_i^\gamma
\end{equation}

and the total amount of area being passed through $C_k$ is at least $|C_k| + n^2 - k^2$.

By Lemma \ref{lem3}, the optimal way to distribute this area over the vertices in $C_k$ is to send all the upstream area, $n^2-k^2$, through one vertex in $C_k$, leaving each of the other vertices in $C_k$ with area 1.  (Note that this is not usually feasible, but it is a lower bound.) \\*
Thus for a given $k$,
\[ \sum_{i \in C_k} A_i^\gamma \ge |C_k| - 1 + (n^2 - k^2 + 1)^\gamma \]
Since $f(x) = x^\gamma$ is concave,
\[\sum_{k=2}^n (n^2 - k^2 + 1)^\gamma \:\ge\: \sum_{k=2}^n [n^{2\gamma} - (k^2 - 1)^\gamma] \:\ge\: n^{2\gamma + 1} - n^{2\gamma} - \sum_{k=2}^n (k^2)^\gamma  \]
So
\[ \sum_{k=2}^n \sum_{i \in C_k} A_i^\gamma
 \:\ge\: \sum_{k=2}^n [2k-2 + (n^2 - k^2 + 1)^\gamma]
  \:\ge\: n^2 - 3n + n^{2\gamma +1} - n^{2\gamma} - \sum_{k=2}^n k^{2\gamma}  \]
For $\gamma = 0.5$, this evaluates to $2n^2 - 4n - (\frac{n^2-n}{2} - 1) \:=\: \frac{3}{2}n^2 - \frac{7}{2}n + 1$.
\end{proof}

\medskip

\subsection{Volume Scaling}

In \cite{Maritan2002}, Maritan et al. make an analytical prediction for the scaling of volume with area in rivers.  Their analysis depends on the relationship $ \langle L_x \rangle \propto A_x^h$, where $ \langle L_x \rangle $ is the mean distance from vertices upstream of $x$ to $x$, and $h$ is Hack's exponent, from Hack's law (a power-law relationship between basin length and area) \cite{Hack57}.  This is combined with the equation $V_x = A_x \langle L_x \rangle $ to arrive at the relationship $V_x \propto A_x^{1+h}$.  This relationship is verified by their data on real river networks as well as OCNs, where (for both of which) $h$ is typically close to $0.57$ \cite{Maritan2002}.

If $h$ were to equal exactly $0.5$, this would imply isometric scaling of length and area in river basins, as $A_x$ would scale directly with $l_x^2$, and shape would be preserved.  This would line up with the prediction in \cite{size_and_form} that $V \propto A_x^{{D+1}/D}$ in the most efficient networks for $D = 2$.  However, since river basins tend to elongate with growth, getting proportionally narrower as they get larger, $h$ is usually slightly greater than $0.5$, and the lower bound predicted in \cite{size_and_form} is rarely reached.

Shifting focus to three dimensions, we can look for evidence of a corresponding Hack's law for three dimensional OCNs (that is, evidence that area scales with length to some constant power).  If there is a relatively constant scaling exponent $h$, the analysis in \cite{Maritan2002} can be extended to three dimensions, since the steps taken there do not otherwise depend on dimension.  Thus, we might expect to see $V_x \propto A_x^{1+h}$ for some $h$; according to the analysis in \cite{size_and_form},  $V_x$ should scale at least as $A_x^{4/3}$, so we might expect $1/3$ to be a lower bound for $h$, with $h$ slightly greater than $1/3$ if three dimensional OCNs elongate in a similar way to two dimensional ones.  These predictions are borne out in the simulation results that follow.
\section{Definitions and Methods}

Consider an area of land, and embed a grid network in it such that every node in the network is associated with a unit area of land.  Construct a spanning tree on the grid and direct it so that all nodes have a path in the tree to a single root (call it the outlet).  Define a partial order on the nodes such that $x \le y $ iff $x$ is on the path from $y$ to the root.  If $x < y$, we will say $y$ is upstream of $x$.  The terms area, volume, upstream length, link length, and energy are defined here in this context.

Note that each node in the spanning tree, except for the root, has exactly one link directed out of it.  Thus we can name each link according to its corresponding node; that is, when we speak of link $x$ we are referring to the link out of node $x$.  All the properties defined here for nodes will be extended to their corresponding links.  For example, the area of link $x$ is defined to be the area of node $x$.

The area of node $x$, $A_x$, is the total number of nodes whose paths to the outlet include $x$; that is, the number of nodes $y$ such that $x \le y$.  This is equivalent to the sum of the areas of each node directly linked to $x$ plus one for $x$ itself:

\begin{equation} \label{eq:area}
	A_x \:=\: \sum_{y} A_y \:+\: 1
\end{equation}
where y ranges over the nodes directly linked to x \cite{Rinaldo1992Minimum}.

The volume of node $x$, $C_x$, is the sum of all the areas of all nodes upstream of $x$ \cite{size_and_form}:
\begin{equation} \label{eq:volume}
	C_x \:=\: \sum_{y | x < y} A_y
\end{equation}

Equation \ref{eq:volume} can be transformed into a recursive definition for volume:

\begin{equation} \label{eq:recvol}
	 C_x \:=\: \sum_{y} (C_y + A_y)
\end{equation}
where y ranges over the nodes directly linked to x.

Upstream length is another quantity of interest; it is defined as the number of links in the path from a node $x$ to a source node, where at each step upstream the path taken goes to the node with the largest area \cite{Cieplak1998Models}.  Statistically, it is equivalent to define it as the distance (counted in number of links) to the farthest away upstream source \cite{Rinaldo2006Trees}.

Link length is a weight assigned to each link, generally representing the geometric length of the link in the embedding.  The length of link $i$ is denoted by $l_i$.

Energy is a function of the link lengths and areas of all the links.  The energy of a given tree configuration $s$ is given by:
\begin{equation} \label{eq:energy}
	H_{\gamma}(s) = \sum_{i \in links(s)} {A_i^{\gamma} l_i}
\end{equation}
with $0 \le \gamma \le 1$.  A tree network that achieves the minimum $H_{\gamma}(s)$ over all possible tree configurations on a given grid is an optimal channel network \cite{Rinaldo1992Minimum}.

Figure \ref{fig:ocnex} shows a grid and a possible tree configuration $s_1$.  $s_1$ minimizes $H_{1/2}(s)$ for all configurations on this grid.  Area, volume and upstream length values are also depicted for the nodes in $s_1$.

\begin{figure}[H]
	\centering
	\begin{subfigure}[c]{0.2\textwidth}
		\centering
		\includegraphics[scale = 0.5]{images/4x4orthgrid.png}
		\caption{Underlying grid \\ { } }
		\label{fig:ocnexgrid}
	\end{subfigure} \hfill
	\begin{subfigure}[c]{0.2\textwidth}
		\centering
		\includegraphics[scale = 0.5]{images/4x4ortharea.png}
		\caption{Areas of nodes \\ in $s_1$}
		\label{fig:ocnexarea}
	\end{subfigure} \hfill
	\begin{subfigure}[c]{0.2\textwidth}
		\centering
		\includegraphics[scale = 0.5]{images/4x4orthvolume.png}
		\caption{Volumes of nodes \\ in $s_1$}
		\label{fig:ocnexvolume}
	\end{subfigure} \hfill
	\begin{subfigure}[c]{0.2\textwidth}
		\centering
		\includegraphics[scale = 0.5]{images/4x4orthupl.png}
		\caption{Upstream lengths \\ of links in $s_1$}
		\label{fig:ocnexupl}
	\end{subfigure} \hfill
	\caption{The outlets are in the top left corners.  Tree $s_1$ minimizes $H_{1/2}(s)$ for the grid in \ref{fig:ocnexgrid}. $H_{1/2}(s_1) = 24.3024964149$.}
	\label{fig:ocnex}
\end{figure}

Rinaldo et al. \cite{Rinaldo1992Minimum} derived a gamma value of 1/2 from estimations of physical properties of rivers, and hypothesized that natural rivers would minimize $H_{1/2}(s)$.  A remarkable result of \cite{Rinaldo1992Minimum} is that while the global optima of OCNs do not match well statistically with real river networks, the local optima are an excellent match.  The idea that nature tends towards a local optimum makes sense, so this suggests that OCNs are a good model for what actually happens in a river \cite{RodriguezIturbe1997Optimal}.

When $\gamma = 1$ and all links are of unit length, the energy is equivalent to the volume.  When $\gamma = 0$, the energy is equivalent to the total length of all links, so minimizing the energy is equivalent to finding a minimum spanning tree (which, when all links are of unit length, would be \emph{any} spanning tree).

In most studies of OCNs, the network of potential links is a 2-dimensional lattice with additional links between diagonally-adjacent pairs of nodes (Fig. \ref{fig:2dgrid}).  For the sake of simplicity, all the links, diagonal and orthogonal, are here considered to have unit length, because it has been shown (\cite{RodriguezIturbe1997Optimal}) that the properties with which this paper is concerned are independent of whether the diagonal links are given realistic length or not.

\begin{figure}[H]
	\begin{subfigure}[b]{0.5\textwidth}
		\centering
		\includegraphics[scale = .5]{images/2dgrid.png}
		\caption{7x7 network of potential links}
		\label{fig:2dgrid}
	\end{subfigure}
	\begin{subfigure}[b]{0.5\textwidth}
		\centering
		\includegraphics[scale = .5]{images/3dgrid.png}
		\caption{2x2x2 3d network of potential links}
		\label{fig:3dgrid}
	\end{subfigure}
	\caption{}
\end{figure}

In this study we also analyzed a three dimensional version of optimal channel networks, built on a 3-dimensional grid (Fig. \ref{fig:3dgrid}).

To make an OCN, we first generated a random spanning tree on a given grid.  We used Prim's algorithm for minimum spanning trees for this.  Essentially, the algorithm builds the tree from a root node, maintaining a set of frontier edges from the grid that are incident on exactly one node in the tree.  Our algorithm randomly chooses one edge at a time from this set, adds it to the tree, and updates the set.  Rinaldo et al. also began with random spanning trees, though they did not specify their method for generating them \cite{Rinaldo1992Minimum}.

The random tree was then optimized using a version of Lin's algorithm for TSP, as in \cite{Rinaldo1992Minimum}.  That algorithm is as follows: \medskip

\begin{algorithm}[H]
\nonumber
\SetAlgoLined
\KwIn{Undirected graph $G$; random spanning tree $T$, directed so that there exists a path in $T$ from each node to a single root}
\KwOut{Locally optimal spanning tree on $G$}

\While{Convergence has not been reached}{
Step 1: Choose a random node $i$ in $T$

Step 2: Choose a random neighbor of $i$ in $G$ distinct from the node $i$

\Indp \Indp links to in $T$

\Indm \Indm Step 3: Construct $T^*$ from $T$ by redirecting $i$ to link to this new node

\lIf{$T^*$ contains a loop}{
	discard $T^*$
}

\lElseIf{$H_{\gamma}(T^*) < H_{\gamma}(T)$}{
	discard $T$ and let $T^*$ be the new $T$
}

\lElse{discard $T^*$}

}
\label{alg:mine}
\end{algorithm}

The condition for convergence was that the number of improvements made be only 1\% of the total number of iterations (or 2\% for larger networks).  See Figure \ref{fig:steps} for an illustration of the process.

\begin{figure}[H]
	\begin{subfigure}[b]{0.5\textwidth}
		\centering
		\includegraphics[scale = .6]{images/3x3step0.png}
		\caption{Initial tree configuration $T$.  The outlet is in the top left corner. $H_{\gamma}(T) = 9.65028$.}

	\end{subfigure} \quad
	\begin{subfigure}[b]{0.5\textwidth}
		\centering
		\includegraphics[scale = .6]{images/3x3step1.png}
		\caption{Step 1: Node $i$ has been chosen; the unused links to its neighbours in $G$ are highlighted.}

	\end{subfigure} \\
	\begin{subfigure}[b]{0.5\textwidth}
		\centering
		\includegraphics[scale = .6]{images/3x3step3.png}
		\caption{Step 3: A new link for $i$ in the tree has been chosen, creating $T^*$.  $H_{\gamma}(T^*) = 9.8637 > H_{\gamma}(T)$, so this tree will be discarded.}

	\end{subfigure} \quad
	\begin{subfigure}[b]{0.5\textwidth}
		\centering
		\includegraphics[scale = .6]{images/3x3final.png}
		\caption{the next iteration, Step 3:  A different new link for $i$ in the tree has been chosen, creating $T^*$.  $H_{\gamma}(T^*) = 9.44949 < H_{\gamma}(T)$, so this tree will become the new $T$.}

	\end{subfigure}
	\caption{Optimization algorithm}
	\label{fig:steps}
\end{figure}

\begin{figure}[H]
	\begin{subfigure}[b]{0.5\textwidth}
		\centering
		\includegraphics[scale = .4]{images/60x60before.png}
		\caption{}
		\label{fig:ocnbefore}
	\end{subfigure}
	\begin{subfigure}[b]{0.5\textwidth}
		\centering
		\includegraphics[scale = .4]{images/60x60after.png}
		\caption{}
		\label{fig:ocnafter}
	\end{subfigure}
	\caption{A random tree spanning a 60x60 grid (a) and the resulting locally optimal OCN (b).  The outlets are in the top left corners.}
\end{figure}

\begin{figure}[H]
	\centering
	\includegraphics[scale = .5]{images/3dOCN.png}
	\caption{An OCN on a 3x3x3 grid; the red links are in the OCN.}
	\label{fig:3docn}
\end{figure}

Figure \ref{fig:ocnbefore} shows a random tree on a 60x60 grid and Figure \ref{fig:ocnafter} shows the result of using the above algorithm to optimize the tree.  In order to make the structure of the network clearer, links are drawn with a thickness proportional to the log of their area.  Figure \ref{fig:3docn} depicts a three-dimensional OCN on a 3x3x3 grid.

A range of grid sizes were analyzed.  Two-dimensional $n$x$n$ grids with $n =$ $10$, $20$, $30$, $40$, $50$, $60$, $70$, and $80$, and three-dimensional $n$x$n$x$n$ grids with $n = 8$, $10$, $12$, $14$, $16$, $18$, and $20$ were used.  For the smaller grids, data was averaged from at least $4$ realizations of each size, while for the larger grids at least $3$ realizations were used for each.  The quantities measured were distributions of area and length in whole basins, and average length, volume, and energy per area of subbasins and whole basins.

\section{Conclusion}

We have looked at three dimensional OCNs as they compare with metabolite distribution networks in organisms.
In the organism context, volume is interpreted as volume of blood (which is directly proportional to the mass or volume of an organism) and area is proportional to the number of capillaries \cite{Savage2008Sizing}.  In this context, the grid used to create OCNs no longer represents fixed space (as it does in the two-dimensional/river context) and the amount of physical area directly associated with each node may not be the same over different body masses.  This would only have a linear effect on the relationship between area and mass, however, so this has no effect on the scaling exponent.

Though the three-dimensional networks do not actually look like the metabolic networks in organisms, we have shown that they have similar area/volume scaling behavior.  This adds credence to the idea that the observed metabolic scaling is a result of characteristics of the distribution network.

Another way to study the possible similarities between OCNs and metabolic networks would be to derive an energy functional for three dimensional OCNs from metabolic networks themselves, rather than using the one derived from rivers.  This could also include consideration of heterogeneity.  Heterogeneity in the river context is addressed in \cite{RodriguezIturbe1997Optimal} and \cite{Colaiori1997Analytical}; possibly this work could be extended to three dimensions.
\section{Experimental Results}

\subsection{Energy Scaling}

\begin{figure}[H]
	\centering
	\includegraphics[scale = .6]{images/AvsE2D.png}
	\caption{Log-log plot of area vs. energy for 2-dimensional OCNs, showing both analytical lower bound (lower, blue) and observed values (upper, red).}
	\label{fig:energies}
\end{figure}

Fig. \ref{fig:energies} shows how empirically observed minimum possible energy scales with area in two-dimensional OCNs.  The estimated line of best fit is shown, as well as the analytical lower bound derived in Section \ref{sec:proofs}.  The observed values of energy follow a power law with a small exponent, while the lower bound is essentially linear with respect to area.

\subsection{Length Scaling}

Scaling exponents $h$ where $l \propto A^h$ in two and three dimensions is shown in Table \ref{tab:length}.  The observed values of the exponents for two dimensions are within the bounds found in other studies \cite{Hack57, Maritan2002}.
For three dimensions, the exponent is also fairly consistent, though less so than in two dimensions.  It also deviates slightly more from the isometric value of 1/3 than the two dimensional version does from 1/2.  The fact that it is higher than 1/3 fulfils the hypothesis that three-dimensional basins elongate like two-dimensional ones.

\begin{table}[H]
	\centering
	\begin{tabular}{| l c |}
		\hline
		 basins analyzed &	$h$ \\ \hline
		2d whole basins &	0.6066895125 \\
		2d all subbasins &	0.5808520176 \\
		3d whole basins &	0.3645961454 \\
		3d all subbasins &	0.4230582779 \\
		\hline
	\end{tabular}
	\caption{Length scaling exponents}
	\label{tab:length}
\end{table}

\subsection{Volume Scaling}

Fig. \ref{fig:2dvol} shows the scaling of volume with area in two-dimensional OCNs, for whole networks and for all subbasins.    Note that the exponents are very similar, indicating that the scaling behavior is the same within basins as across different sizes of basins.  The exponents are both close to 1.57, as in \cite{Maritan2002}.

\begin{figure}[H]
	\begin{subfigure}[b]{0.5\textwidth}
		\centering
		\includegraphics[scale = .45]{images/AvsC2D/wholes.png}
		\caption{}
		\label{fig:wholes2}
	\end{subfigure}
	\begin{subfigure}[b]{0.5\textwidth}
		\centering
		\includegraphics[scale = .262]{images/AvsC2D/subsall.png}
		\caption{}
		\label{fig:2}
	\end{subfigure}
	\caption{Volume plotted with area for whole nxn basins (\ref{fig:wholes2}) and all subbasins (\ref{fig:2}) of 2-dimensional OCNs.}
	\label{fig:2dvol}
\end{figure}

\begin{figure}[H]
	\begin{subfigure}[b]{0.5\textwidth}
		\centering
		\includegraphics[scale = .45]{images/AvsC3D/wholes.png}
		\caption{Whole nxn basins}
		\label{fig:wholes}
	\end{subfigure}
	\begin{subfigure}[b]{0.5\textwidth}
		\centering
		\includegraphics[scale = .45]{images/AvsC3D/subsall.png}
		\caption{all subbasins}
		\label{fig:3dsubs}
	\end{subfigure}
	\begin{subfigure}[b]{0.5\textwidth}
		\centering
		\includegraphics[scale = .45]{images/AvsC3D/subs10.png}
		\caption{subbasins of 10x10x10 OCNs}
		\label{fig:3dsubs10}
	\end{subfigure} \quad
	\begin{subfigure}[b]{0.5\textwidth}
		\centering
		\includegraphics[scale = .45]{images/AvsC3D/subsallsmall.png}
		\caption{all subbasins with Area $\le$ 1000}
		\label{fig:3dsubssmall}
	\end{subfigure}
	\caption{Volume plotted with area for 3-dimensional OCNs.}
	\label{fig:3dvol}
\end{figure}

\begin{table}[H]
	\centering
	\begin{tabular}{| l c |}
		\hline
		basins analyzed & $\alpha$ \\ \hline
		whole networks &	1.386 \\
		all subbasins &	1.424 \\
		all subbasins of 20x20x20 network &	1.437 \\
		subbasins with area $>$ 1000 &	1.3776 \\
		subbasins of 20x20x20 network with area $>$ 1000 &	1.342 \\
		all subbasins with area $\le$ 1000 &	1.44456 \\
		subbasins of 20x20x20 network with area $\le$ 100 &	1.457 \\
		subbasins of 10x10x10 network &	1.4506 \\
		\hline
	\end{tabular}
	\caption{Volume scaling exponents for three-dimensional OCNs}
	\label{tab:3dvol}
\end{table}

The three-dimensional results for volume scaling are displayed in Fig. \ref{fig:3dvol} and summarized in Table \ref{tab:3dvol}.    The differences in exponents between the different data sets show higher scaling behavior in the sets of smaller basins; the difference between the exponent for whole basins (\ref{fig:wholes}) and all subbasins (\ref{fig:3dsubs}) is due to the greater concentration of smaller basins in the set of all subbasins, skewing the slope.  Note that the exponents in \ref{fig:3dsubs10} and \ref{fig:3dsubssmall} are similar, showing that scaling behavior in the subbasins is the same for the same distributions of area independent of the size of the encompassing whole basin.  In Section \ref{sec:proofs} we predicted that $\alpha$ would equal $1 + h$, using the analysis in \cite{Maritan2002}, and this relationship is clear in Tables \ref{tab:length} and \ref{tab:3dvol}.

The exponent $\alpha$ ranges from $\approx$1.34 to $\approx$1.46, with a higher exponent for collections of smaller basins.  These results are not inconsistent with the data for metabolic allometry in biological organisms, where a power-law fit to data from a wide range of sizes approximates that metabolic rate scales with mass to the power of 0.70 ($1/\alpha$) \cite{Savage2008Sizing}.  Not only is this value within the range found in the OCNs, but in organisms as well it is found that sets of smaller organisms yield a higher exponent than sets of larger ones \cite{Kolokotrones2010Curvature, Savage2004Predominance}.

\section{Introduction}

Fractal Networks have been widely used to study well known network behaviors.  Some of the properties commonly used in the study of fractal networks are scaling factors (or laws); scaling laws can describe how lengths
are related to area and areas are related to volumes. In this paper, we further study these allometric scaling laws.  We use simulation and elementary mathematical bounding techniques to study such laws.

A common pattern in nature, known as metabolic allometry, is a 3/4-power scaling law between metabolic rate ($M$) and body mass ($B$) in organisms, across and within species: $B \propto M^{3/4}$.  An analogous relationship between water mass ($C$) and contributing area ($A$) in river networks has also been observed: $A \propto C^{2/3}$.

The relationship in biological organisms was first observed by Max Kleiber in 1932 and has garnered a lot of attention in the last decade and a half from biologists, mathematicians, and physicists alike \cite{Kleiber1932Body}.  The areas this relationship finds application in include models of climate change effects in rivers \cite{Renner2012Evaluation} and oceans \cite{ocean_balance}, analyses of river stability \cite{Molnar1998Analysis}, and in studies of information networks \cite{Moses2008Scaling}.  Many theories concerning the origins and nature of this phenomenon have been proposed, and there has been much debate about how ubiquitous the scaling law really is.

In 1997, West, Brown and Enquist (WBE hereafter) introduced a theory purporting to explain metabolic scaling \cite{wbe}. Their theory is based on the idea of a fractal space-filling hierarchical network of resource distribution, with certain assumptions based on simplified biology.  The geometry of the network (e.g. branching ratios, vessel radius relationships) are important to the development of their result. The authors claimed the model is applicable even to organisms without a physical blood vessel network since they can be treated as having a virtual distribution network \cite{West1999Fourth}.

Subsequent studies questioned the WBE results, by showing that data does not appear to strictly follow a 3/4-power scaling law after all; the theoretical approach taken by WBE has also been criticized \cite{Agutter2011Analytic, White2011Manipulative}.  Some authors claim that an exponent of 2/3 fits the data better \cite{Dodds2001Reexamination}, while others assert that no universal scaling law exists at all \cite{Kolokotrones2010Curvature}.
Savage et al. pointed out that the WBE model predicts a curvilinear relationship between $log B$ and $log M$, and only predicts a constant power-law relationship in the infinite limit of body size; however, even with finite size adjustments and variations in the structural properties of the model, they found the results were still inconsistent with data \cite{Savage2008Sizing}.
Some other modifications to the WBE model do yield results somewhat more consistent with data, such as taking into account body temperature \cite{Brown2004Toward} or fluid velocity \cite{Banavar2010General}.
Kolokotrones et al. use regression on large set of organism data to show a quadratic, rather than linear, relationship between $log B$ and $log M$ \cite{Kolokotrones2010Curvature}. Also, their results reported a lot of variation.  They build off the modifications to WBE in \cite{Savage2008Sizing} and find that a model network with a proportional transition between area-increasing and area-preserving branching yields a fit to the data almost as good as their empirically-derived quadratic model.

Other theories and models have been put forth, some based on network properties like the WBE model, and others on other biological properties.  Dodds \cite{Dodds2010Optimal} proposed a model based on virtual networks of metabolite transportation rather than physical networks.

Banavar et al. proposed a general model for any efficient transportation network, predicting $B \propto M^{D/(D+1)}$ (at most) where $D$ is the dimension of the network (i.e. $D = 3$ for organism metabolism, and $D = 2$ for river networks) \cite{size_and_form}.  Thus, they predict that if the network were as efficient as possible, organism metabolism should scale as $B \propto M^{3/4}$, and deviations of the exponent below $3/4$ might be explained by some inefficiency in the metabolic distribution network.  Later the authors added to this model a supply and demand mechanism that accounts for the deviations \cite{Banavar2002Supplydemand}.

Isaac et al. propose a less strict understanding of the ``3/4-power law'': though the 3/4 exponent might not be a consistent, universal constant, it is still an observable trend; if there isn't a single unifying principle explaining why it appears, there still may be multiple reasons for it worth studying \cite{Isaac2010Why}.

In the context of river networks, the analogous properties to metabolic rate and body mass are contributing area (A) and water mass (C).  It is well-observed that the properties obey a power-law relationship with an exponent of slightly more than 2/3: $A \propto C^{2/3}$.  This is consistent with the general resource transportation model of Banavar et al. \cite{size_and_form} for systems in two dimensions, and also with the virtual network model of Dodds \cite{Dodds2010Optimal}.

Rinaldo et al. devised a model of river networks called Optimal Channel Networks (OCNs) \cite{Rinaldo1992Minimum}.  This model uses a grid network to represent the area of a river basin, and constructs a spanning network on the grid that minimizes a functional representing the energy required to drain the area of land through river channels.  The energy functional is derived both from theoretical physical characteristics and observed characteristics of river basins \cite{RodriguezIturbe1997Optimal}.  Statistical properties of Optimal Channel Networks (including metabolic allometry) correspond excellently with observed properties of real river networks \cite{Maritan2002, RodriguezIturbe2011Metabolic}.

This paper explores a three-dimensional version of OCNs, using simulations to compare characteristics of 3d OCNs with 2d OCNs and with organic metabolic networks and study the scaling behavior of area, length, volume, and energy.  We do not attempt to explain 3/4-power scaling in nature, but rather explore some interesting and potentially useful connections, somewhat in the spirit of Isaac et al.'s proposed understanding of the scaling law \cite{Isaac2010Why}.

We further study an alternative model using Steiner trees for OCNs and we compare the scaling behaviors between these two models on small network sizes.

Section 2 describes definitions and methodologies used in our paper.
In Section 3, we derive lower bound results for various scaling laws.
In section 4, we describe our simulation results and show Energy, Length and Volume Scaling with respect to area
and compare it with the lower bound results described in Section 3. In Section 5, we describe an alternative model based on Steiner Trees and derive analytical results for small network sizes. Section 6 gives conclusions and suggestions for future work.


\section{Introduction}

The controllable generation and manipulation of single to few-particle
excitations in mesoscopic conductors has attracted much attention
recently
\cite{hermelin_electrons_2011,mcneil_-demand_2011
FeveMaheBerroirKontosGlattliEtienneSCIENCE316-07
mahe_current_2010,bocquillon_electron_2012
parmentier_current_2012
sherkunov_optimal_2012
gabelli_shaping_2012
VanevicNazarovBelzigPRL99-07,*VanevicNazarovBelzigPRB78-08}
The interest in the subject stems from the emerging field of
electron quantum optics, which has the goal to
develop coherent electronic devices suitable for
quantum communication and quantum computation.
This is based on the analogy with photon quantum optics where pairs
of photons in entangled polarization state are used for
transmission and processing of quantum information. In a
solid-state system, a similar functionality can be achieved using
pairs of quasiparticles with entangled spin or orbital degrees of
freedom
\cite{beenakker_proposal_2003,samuelsson_two-particle_2004}
An important building block of electronic devices aimed at quantum
computation is a tunable electron source capable of creating
quasiparticles in a coherent way. An on-demand coherent electron
source has been realized
recently using a localized electronic
level in a quantum dot, weakly coupled to a conductor, which is
populated and emptied by a modulation of its energy
via a periodic gate
voltage.\cite{FeveMaheBerroirKontosGlattliEtienneSCIENCE316-07}
This results in a sequence of quantized single-electron current
pulses whose properties can be inferred from the measurements
of the current noise spectrum
\cite{parmentier_current_2012,KeelingShytovLevitovPRL101-08}

An altogether different route which does not require electron
confinement has been suggested theoretically by Keeling, Klich,
and Levitov.\cite{KeelingKlichLevitovPRL97-06}
The authors have proposed to make use of time-dependent voltage
pulses to create excitations from a degenerate
Fermi sea in a mesoscopic conductor. It has been found that
Lorentzian-shaped pulses $V(t)$ of a quantized area
$\int eV(t) dt/\hbar = 2\pi N$ ($N$ is integer, $e$ is the electron
charge, and $\hbar$ is the reduced Planck constant,
hereafter $\hbar=1$) create exactly $N$
electrons above the Fermi level with no additional excitations.
The many-body quantum state created by these pulses is a product
state of $N$ particles added to an unperturbed Fermi sea,
independent on the relative position of the pulses, their duration, or
overlap. Conversely, the charge transfer statistics in this case
is binomial, which supports the picture of independent charge
quanta created. The single-particle character of excitations can
be probed by noise measurements: it manifests itself as a
reduction of the current noise power which assumes its minimal
value set by the dc
voltage offset.\cite{IvanovLeeLevitovPRB56-97}
In general, however, a time-dependent field does create additional
excitations in a Fermi sea. These excitations give rise to excess
photon-assisted noise, exceeding the minimal dc
noise level.\cite{LesovikLevitovPRL72-94,PedersenButtikerPRB58-98}
The photon-assisted noise was observed experimentally in
diffusive phase-coherent metallic conductors
\cite{SchoelkopfKozhevnikovProberRooksPRL80-98}
normal-metal--superconductor junctions
\cite{KozhevnikovProberPRL84-99}
and quantum point contacts\cite{ReydelletRocheGlattliEtienneJinPRL90-03}
with harmonic ac-voltage applied.
Noise spectroscopy of a quantum tunnel junction
with a more complex biharmonic voltage drive
has been carried out recently.\cite{gabelli_shaping_2012}

The physical picture behind photon-assisted noise has been
revealed in Ref
~\onlinecite{VanevicNazarovBelzigPRL99-07,*VanevicNazarovBelzigPRB78-08}
by an analysis of the full counting statistics: The elementary
excitations in the system at low temperature are the
{\em electron-hole pairs} created by the ac voltage component,
in addition to electrons injected by the dc voltage offset. This
picture is valid beyond noise measurements and pertains to the
full statistics of the transferred charge.
Since electrons and holes from a pair are created with the same
probability, the pairs give rise only to the excess noise and
higher even-order current correlators, whereas they give no
contribution to the average current and odd-order correlators.
The number of created electron-hole pairs depends on
the shape and the amplitude of the ac voltage applied.
This opens a route toward dynamic control of elementary
excitations in quantum conductors which, if proved feasible,
could be used in the electron-hole sources to produce quasiparticles
with entangled spin or orbital degrees of freedom
\cite{sherkunov_optimal_2012
RychkovPolianskiButtikerPRB72-05
PolianskiSamuelssonBuettikerPRB72-05,samuelsson_dynamic_2005
beenakker_optimal_2005,LebedevLesovikBlatterPRB72-05
splettstoesser_two-particle_2009,moskalets_spectroscopy_2011
grenier_single-electron_2011}
Alternatively, such control can also be used to minimize
excitations present in the system and approach the limit of an
ideal single-electron source using realistic voltage pulses.

In this paper, we investigate the feasibility of a dynamic control
of elementary electron-hole pair excitations.
We analyze experimental data on the photon-assisted noise in
quantum conductors subject to time dependent drive
\cite{ReydelletRocheGlattliEtienneJinPRL90-03,gabelli_shaping_2012}
identify elementary excitations generated in the experiments,
and show how the measured excess noise is composed of the
contributions of electron-hole pairs created. For a quantum
contact with harmonic voltage drive studied in
Ref.~\onlinecite{ReydelletRocheGlattliEtienneJinPRL90-03},
we find that a single electron-hole pair is created per period
with a certain probability, in addition to the electrons that are
injected by dc voltage offset. We find how the probability of pair
creation depends on the amplitude of the ac voltage drive, which
is in agreement with the observed excess noise.

In the recent experiment by Gabelli and
Reulet,\cite{gabelli_shaping_2012}
the authors have studied photon-assisted noise in a quantum tunnel
junction subject to a biharmonic time-dependent voltage, where the
dc offset, the ac amplitudes, and the relative phase are tunable.
The authors have observed the
reduction of the noise when an in-phase second harmonic is present
in the drive. We find that the statistics of the transferred
charge is the simplest when the dc offset is an integer multiple
of the driving frequency. In that case only a few (one or two)
electron-hole pairs are created per period with probabilities that
depend on the shape and the amplitude of the biharmonic voltage
component. In addition,  we relate the reduction of the noise in
the presence of an in-phase second harmonic to the suppressed
probability of the electron-hole pair creation as the drive
voltage approaches the shape of optimal Lorentzian pulses. The
agreement with the observed excess noise corroborates that the
dynamic control of elementary excitations in quantum conductors
has been achieved experimentally.

The system we study is a coherent mesoscopic conductor
characterized by transmission eigenvalues $\{ T_n \}$
where $n$ labels spin-degenerate transport channels. Conductance
of the conductor is $G=(e^2/\pi)\sum_n T_n$ and the
Fano factor $F=(\sum_n T_n R_n)/\sum_n T_n$ where $R_n=1-T_n$ are
reflection probabilities. The conductor is subject to a
periodic voltage drive $V(t) = \bar V + V_\ac(t)$ where
$\bar V$ is dc voltage offset and $V_\ac(t)$ is the ac
voltage component with zero average and period $\gt = 2\pi/\go$.
\begin{figure}[t]
\includegraphics[scale=0.9]{pkV1a.eps}
\caption{\label{fig:pkV1a.eps}
    Probabilities $p_k$ of electron-hole pair creations for
    harmonic drive $V(t) = \bar V + V_1 \cos(\go t)$ as a function
    of the ac amplitude $V_1$ and for different dc voltage
    offsets $e\bar V/\go = N$:
    (a) $N=0$, (b) $N=1$, and (c) $N=2$.
    DC offset $\bar V$ sets the number of electrons injected
    towards the contact while $V_1$ controls the electron-hole pair
    generation.

    Symbols denote experimental data for the excess noise
    $S_\ac$ in units $2GF\go$ taken from
    Ref.~\onlinecite{ReydelletRocheGlattliEtienneJinPRL90-03}
    for a quantum contact with $F=1/2$ and
    ({\scriptsize $\bigcirc$}) $\go/2\pi = 17.32\, \rm GHz$,
    ($\square$) $\go/2\pi = 8.73\, \rm GHz$.
}
\end{figure}
At low temperature $T\ll \go$, the cumulant generating function
$\cS(\gx)$ of the charge transfer statistics can be cast in a
form that manifestly reveals what are the elementary processes
in the system:\cite{VanevicNazarovBelzigPRL99-07}
$\cS(\gx) = \cS_\dc(\gx) + \cS_\ac(\gx)$, where
$\cS_\dc = (t_0 |e\bar V|/\pi)
\sum_n \ln[ 1 + T_n(e^{-i\gk\gx}-1) ]$
and
\begin{align}
\cS_\ac
=&
\frac{2 t_0}{\gt} \sum_{n,k}
\left( (1-\bar v)
\ln[ 1 + p_k^{(N)} T_n R_n (e^{i\gx} + e^{-i\gx}-2)]
\right.  \notag
\\
&+
\left.
\bar v \, \ln[ 1 + p_k^{(N+1)} T_n R_n (e^{i\gx} + e^{-i\gx}-2)]
\right).
\end{align}
Here, $t_0$ is the measurement time which is much larger than the
characteristic time scale on which current fluctuations are
correlated,
$\gk = 1$ ($\gk = -1$) for $e\bar V>0$ ($e\bar V<0$)
is related to direction of the charge transfer,
and $N = \lfloor e\bar V/\go \rfloor$
($\bar v = e\bar V/\go - N$) is the integer (fractional) part
of $e\bar V/\go$. Coefficients
$p_k^{(N)}$ ($0\le p_k^{(N)} \le 1$) in $\cS_\ac$ are the
probabilities of electron-hole pair creations
which depend on the details of the time-dependent voltage applied,
see Fig.~\ref{fig:pkV1a.eps}.
They are given by $p_k^{(N)} = \sin^2(\ga_k/2)$,
where $e^{\pm i \ga_k}$ are the {\em pairs} of complex conjugate
eigenvalues of the matrix
$M_{nm} = \sgn(n+0^+) \sum_{k=-\infty}^\infty
a_{n+k} a^*_{m+k} \sgn(0^+-k-N)$
(cf. Ref.~\onlinecite{VanevicNazarovBelzigPRL99-07,VanevicNazarovBelzigPRB78-08}
for details). The Fourier coefficients $a_n=\gt^{-1}\int_0^\gt
dt\, e^{-i\phi(t)}e^{in\go t}$ characterize the drive, where
$\phi(t)=\int_0^t dt'\, eV_\ac(t')$ is the phase acquired.

\begin{figure}[t]
\includegraphics[scale=0.8]{SIvsVb.eps}
\caption{\label{fig:SIvsVb}
    (a) Current noise $S$ as a function of dc voltage $\bar V$ for
    biharmonic time-dependent drive $V(t) = \bar V + V_1 \cos(\go
    t) + V_2 \cos(2\go t + \gvf)$ with $eV_1/\go = 5.4$, $eV_2/\go
    = 2.7$, and $\gvf = 0$ (thick solid line), $\pi/2$
    (dash-dotted line), and $\pi$ (thin solid line). The excess
    noise $S_\ac$ is shown as a function of $\bar V$ for the same
    voltage drive with (b) $\gvf = 0$ and (c) $\gvf = \pi/2$.
    Symbols in (b) and (c) denote experimental data taken from
    Ref.~\onlinecite{gabelli_shaping_2012} for a tunnel
    junction ($F=1$) at low temperature ($T=0.14\go$).
}
\end{figure}

\begin{figure}[t]
\includegraphics[scale=0.9]{pkPhi0_Pid2a.eps}
\caption{\label{fig:pkPhi0_Pid2a}
    Decomposition of the excess noise $S_\ac$ shown in
    Fig.~\ref{fig:SIvsVb} into elementary electron-hole pair
    excitations [see Eq.~\eqref{eq:Sac-decomp-pk}].
    Probabilities $p_k^{(N)}$ of the pair creations are denoted by
    stacked bars at dc voltages $e\bar V/\go = N$ with $N$
    integer: $k=1$ (blue), $k=2$ (green), and $k=3$ (red bars).
    The elementary excitations comprise a few (one or two)
    electron-hole pairs generated with probabilities
    $p_{1,2}^{(N)}$  per voltage cycle.
}
\end{figure}

The physical interpretation of $\cS(\gx)$ is as follows.
The dc part $\cS_\dc(\gx)$ describes unidirectional single
electron transfers due to a finite voltage offset $\bar V$.
The attempt frequency is $|e\bar V|/\pi$ and the transfer
probability $T_n$ per transport channel.
The ac part of the cumulant generating function $\cS_\ac(\gx)$
describes the events of electron-hole pair creations.
The electron-hole pairs (labelled by $k$) are created by
the ac component of the drive with probability $p_k$ per voltage cycle.
The charge transfer occurs with probability $p_kT_nR_n$
when one particle -- e.g., electron -- is transmitted and the
hole is reflected, or vice versa. Since the electron and hole
from a pair are transmitted with equal probabilities,
the electron-hole pairs give no contribution to the average
current but they do contribute to the current noise power.
We note that the total number of attempts to create an
electron-hole pair is $t_0/\tau$ (per spin-split transport
channel), which gives one attempt per voltage cycle. These
attempts are shared
between processes $p_k^{(N)}$ and $p_k^{(N+1)}$ that correspond
to the nearest integer values of $e\bar V/\go$. The simplest
statistics -- determined by one type of the processes $p_k^{(N)}$
only -- is obtained for the voltage offset $e\bar V = N\go$
which is an integer multiple of the driving frequency.
For a nonzero dc component, the created electron-hole pairs
coexist with $N$ electrons injected towards the contact per period.

From the cumulant generating function we obtain
the current noise power $S=(e^2/t_0)\partial^2_{i\gx}\cS|_{\gx=0}
=S_\dc+S_\ac$ with $S_\dc = GF|e\bar V|$
and
$ S_\ac = 2GF\go \sum_k
\left( (1-\bar v)p_k^{(N)} + \bar v p_k^{(N+1)} \right)$.
Therefore, the excess noise of the electron-hole pairs
is a piecewise linear function of a dc voltage offset.
At $e\bar V/ \go = N$ ($N$ is integer), $S_\ac$
is given by the sum of the probabilities of the pair creations
\begin{equation}\label{eq:Sac-decomp-pk}
S_\ac|_N = 2GF\go \sum_k p_k^{(N)}
\end{equation}
and linearly interpolates in-between for $e\bar V/\go$ noninteger.
The experimentally observed excess nois
\cite{ReydelletRocheGlattliEtienneJinPRL90-03,gabelli_shaping_2012}
and its decomposition into contributions of the electron-hole
pairs is shown in Figs.~\ref{fig:pkV1a.eps}~--~\ref{fig:pkPhi0_Pid2a}.

In Ref.~\onlinecite{ReydelletRocheGlattliEtienneJinPRL90-03}, the
authors measured photon-assisted noise in a quantum
contact ($F=1/2$) subject to harmonic time-dependent drive
$V(t) =\bar V + V_1 \cos(\go t)$. For this drive
$a_n = J_n(eV_1/\go)$ where $J_n(x)$ are Bessel functions of the
first kind. The predicted probabilities $p_k^{(N)}$ of the
electron-hole pair creations are shown in Fig.~\ref{fig:pkV1a.eps}
as a function of the ac amplitude $V_1$ and for different values
$e\bar V/\go =N$ of the dc component. As the drive amplitude $V_1$
is increased, the probability of the pair creation also increases
and more pairs are created per period. The symbols in
Fig.~\ref{fig:pkV1a.eps} denote experimental data for the
excess noise $S_\ac/2GF\go$ at dc offsets
$e\bar V/\go = N$ ($N=0,1,2$), measured for
driving frequencies $\go/2\pi = 17.32\,$GHz and $8.73\,$GHz.
The experimental results are in agreement with
Eq.~\eqref{eq:Sac-decomp-pk}. We find that for the amplitudes
used in the experiment only a single electron-hole pair is
created per period with probability $p_1$. The interval
$eV_1/\go$ where a single pair is excited is broad and is not
restricted to the weak driving limit $eV_1/\go\ll 1$, see
Fig.~\ref{fig:pkV1a.eps}~(a).
For a nonzero dc offset, the pair creation is also
accompanied by the single-electron transfers.
When the ac amplitude is smaller than the dc offset,
the creation of electron-hole pairs is strongly suppressed, see
Fig.~\ref{fig:pkV1a.eps}~(b),(c).
In the situation in which the probability of a pair creation is
small, the noise is dominated by the dc shot noise
$S_\dc \gg S_\ac$. In that case, generated electron-hole pairs can
be probed more directly in a beam splitter geometry where current
{\it cross correlation} occurs due to processes in which the pairs
are split and the two particles enter different outgoing terminals
\cite{VanevicNazarovBelzigPRB78-08,RychkovPolianskiButtikerPRB72-05}

In a recent experiment, Gabelli and
Reulet\cite{gabelli_shaping_2012} studied the noise in a
quantum tunnel junction ($F=1$) with a biharmonic voltage applied,
$V(t) = \bar V + V_1 \cos(\go t) + V_2 \cos(2\go t + \gvf)$,
where the dc offset $\bar V$, the ac amplitudes $V_1$ and $V_2$,
and the phase shift $\gvf$ were tunable.
The coefficients $a_n$ in this case are given by
$a_n=\sum_{m=-\infty}^\infty
J_{n-2m}(eV_1/\go)J_m(eV_2/2\go)e^{-im\gvf}$.
The observed current noise power and the decomposition of the
excess noise into contributions of the electron-hole pairs are
shown in Figs.~\ref{fig:SIvsVb} and~\ref{fig:pkPhi0_Pid2a}.
\begin{figure}[t]
\includegraphics[scale=0.9]{pkV2Phi0_Pid2.eps}
\caption{\label{fig:pkV2Phi0_Pid2}
    Probabilities $p_k$ of the electron-hole pair creations as a
    function of the amplitude $V_2$ of the second harmonic for the
    biharmonic voltage drive with $e\bar V/\go = 3$, $eV_1/\go =
    5.4$, and the phase difference (a) $\gvf=0$ and (b) $\gvf =
    \pi/2$. Dashed lines indicate the excess noise $S_\ac$ in
    units $2GF\go$. The presence of an in-phase second harmonic
    decreases the probability of the electron-hole pair creation
    leading to suppression of $S_\ac$ observed
    experimentally.\cite{gabelli_shaping_2012}
}
\end{figure}
The current noise power as a function of $\bar V$ and for different
phase shifts $\gvf$ is shown in Fig.~\ref{fig:SIvsVb} (a),
where we use the experimental values of the parameters
$eV_1=5.4\go$ and $eV_2=2.7\go$. At low temperature $T\ll \go$,
the current noise power reads\cite{LesovikLevitovPRL72-94}
$S=GF\sum_{n=-\infty}^\infty |e\bar V + n\go|\, |a_n|^2$,
which is a piecewise linear function of the dc voltage offset
$e\bar V$ with kinks at integer multiples of $\go$.
This can be seen in the differential noise
$S'=\partial S/\partial(e\bar V)$ which consists of a series
of steps at $e\bar V/\go = N$ ($N$ is integer)
with the step size
$S'|_{N+0}-S'|_{N-0} = 2GF|a_{-N}|^2$ (cf. Fig.~3 in
Ref.~\onlinecite{gabelli_shaping_2012}).
The excess noise is $S_\ac|_N = 2GF\go \sum_{n=1}^\infty n
|a_{-N\mp n}|^2$ where the upper (lower) sign is taken for
$N\ge 0$ ($N<0$). This analysis quantitatively accounts
for the photon-assisted noise observed experimentally,
see Fig.~\ref{fig:SIvsVb} (b), (c).
However, the information on elementary excitations present in the
system is implicit and encoded in the Fourier components $a_n$
characterizing the drive.

The decomposition of $S_\ac$ into contributions of electron-hole
pair excitations given by Eq.~\eqref{eq:Sac-decomp-pk} is shown in
Fig.~\ref{fig:pkPhi0_Pid2a}. The probabilities $p_k^{(N)}$ of the
pair creations are depicted as stacked bars. We observe that the
excess noise consists of a small number (one or two) of the
electron-hole pairs excited per period with certain probabilities.
The pairs are accompanied with $e\bar V/\go = N$ electrons
injected towards the contact per period. Therefore, for
small offset voltages the number of electron-hole pairs is
significant as compared to the number of electrons injected per
period. As the dc offset is increased, the probability of pair
creation decreases and the excess electrons start to dominate.
The probability of pair creation can further be tuned by
changing the shape of the ac voltage component.

In Ref.~\onlinecite{gabelli_shaping_2012}, the authors have used
a setup where the dc offset and the amplitude of the first harmonic
are kept fixed and the amplitude of the second harmonic $V_2$ is
varied to minimize the noise. In Fig.~\ref{fig:pkV2Phi0_Pid2} we
show the probabilities of the electron-hole pair generation as a
function of $V_2$ for a dc offset $e\bar V/\go = 3$
and the amplitude of the first harmonic $e V_1/\go = 5.4$.
For a simple harmonic drive ($V_2=0$) there is only one
electron-hole pair generated per period with probability $p_1$.
When the second harmonic is in phase with the first one
($\gvf=0$), the probability $p_1$ of the pair generation {\em
decreases} as the amplitude $V_2$ is increased,
Fig.~\ref{fig:pkV2Phi0_Pid2}~(a). As $V_2$ is increased further,
the second electron-hole pair is generated with the
increasing probability $p_2$ per period. For the amplitude
$eV_2/\go \approx 2.6$, the total
probability $p_1+p_2$ of the electron-hole pair creation
exhibits a minimum. This leads to the minimal excess noise in
Eq.~\eqref{eq:Sac-decomp-pk} which has been observed in
Ref.~\onlinecite{gabelli_shaping_2012}.
On the other hand, when the phase difference between the first and
the second harmonic is $\gvf=\pi/2$, the shape of the drive
deviates from the optimal one as $V_2$ is increased. The
probability of the electron-hole pair generation grows
monotonically with $V_2$ [see Fig.~\ref{fig:pkV2Phi0_Pid2}~(b)]
and the excess noise increases.\cite{gabelli_shaping_2012}

In conclusion, we have studied the available experimental data on
quantum noise in mesoscopic conductors with applied time-dependent
voltages and provided an interpretation
in terms of
independent electrons and the electron-hole pairs created by the
drive. This interpretation is valid for a generic quantum
junction and to all orders in the statistics of the transferred
charge. We have shown how the excess photon-assisted noise is
composed of the contributions of the electron-hole pairs, whose
number and probability of creation can be controlled by changing
the shape and amplitude of the applied voltage. The agreement
of the predicted pair creation probabilities and the observed
photon-assisted noise corroborates that the dynamic control
of elementary excitations
has been achieved experimentally. This can be utilized in
electron-hole sources that emit quasiparticles with entangled spin
or orbital degrees of freedom for use in mesoscopic electronics
and electron quantum optics.

We are grateful to Julien Gabelli, Yuli V. Nazarov, and Bertrand
Reulet for fruitful discussions. The research was supported by the
Serbian Ministry of Science, Project No.~171027, and the DFG
through SFB 767 and SPP1285.