text:
\section{Introduction}

We consider a network of $n$ sensors, each one placed at a fixed position in
two-dimensional space and capable of communicating with another sensor if and
only if the Euclidean distance between the two is at most $R$, for some constant
radius $R>0$. If $\delta_{ij}$ denotes this distance for sensors $i$ and $j$,
then a graph representation of the network can be obtained by letting each
sensor be a node and creating an edge between any two distinct nodes $i$ and $j$
such that $\delta_{ij}\le R$. Such a representation is, aside from a scale
factor, equivalent to a unit disk graph \cite{disk}.

Often $n$ is a very large integer and the network is essentially unstructured,
in the sense that the sensors' positions, although fixed, are generally unknown.
In domains for which this holds, generalizing the graph representation in such a
way that each node's position is given by random variables becomes a crucial
step, since it opens the way to the investigation of relevant distributions
related to all networks that result from the same deployment process. Such a
generalization, which can be done for any number of dimensions, is known as a
random geometric graph \cite{penrose}. Similarly to the random graphs of
Erd\H{o}s and R\'{e}nyi \cite{erdos} and related structures \cite{nsw01}, many
important properties of random geometric graphs are known, including some
related to connectivity and the appearance of the giant component
\cite{appel_connect,appel_maxdegree,appel_mindegree} and others more closely
related to applications \cite{percolation,gupta,mcdiarmid}.

One curious aspect of random geometric graphs is that, if nodes are positioned
uniformly at random, the expected Euclidean distance between any two nodes is a
constant in the limit of very large $n$, depending only on the number of
dimensions (two, in our case) \cite{box}. In this case, distance-dependent
analyses must necessarily couple the Euclidean distance with some other type of
distance between nodes. The natural candidate is the standard graph-theoretic
distance between two nodes, given by the number of edges on a shortest path
between them \cite{bondy}. For nodes $i$ and $j$, this distance is henceforth
denoted by $d_{ij}$ and referred to simply as the distance between $i$ and $j$.

Given $i$ and $j$, the Euclidean distance $\delta_{ij}$ and the distance
$d_{ij}$ between the two nodes are not independent of each other, but rather
interrelate in a complex way. Our goal in this paper is to explore the
relationship between the two when all sensors are positioned uniformly at random
in a given two-dimensional region. Specifically, for $i$ and $j$ two distinct
nodes chosen at random, we study the probability that $d_{ij}=d$ for some
integer $d>0$, given that $\delta_{ij}=\delta$ for some real number
$\delta\ge 0$. Similarly, we also study the probability density associated with
$\delta_{ij}=\delta$ when $d_{ij}=d$. Our study is analytical whenever feasible,
but is also computational throughout. Depending on the value of $d$, we are in a
few cases capable of providing exact closed-form expressions, but in general
what we give are approximations, either derived mathematically or inferred from
simulation data exclusively.

We remark, before proceeding, that we perceive the study of distance-related
distributions for random geometric graphs as having great applicability in the
field of sensor networks, particularly in domains in which it is important for
each sensor to have a good estimate of its location. In fact, of all possible
applications that we normally envisage for sensor networks \cite{instrumenting},
network localization is crucial in all cases that require the sensed data to be
tagged with reliable indications of where the data come from; it has also been
shown to be important even for routing purposes \cite{routing}. So, although we
do not dwell on the issue of network localization anywhere else in the paper, we
now digress momentarily to clarify what we think the impact of distance-related
distributions may be.

The problem of network localization has been tackled from a variety of
perspectives, including rigidity-theoretic studies \cite{rigid,theory},
approaches that are primarily algorithmic, either centralized
\cite{convex,local_connect} or distributed
\cite{range_free,dv_based,global_coord,robust_noisy}, and others that generalize
on our assumptions by taking advantage of sensor mobility
\cite{local_mobile,robust_noisy} or uneven radii \cite{mds}. In general one
assumes the existence of some anchor sensors (regularly placed \cite{gps-less}
or otherwise), for which positions are known precisely, and then the problem
becomes reduced to the problem of providing, for each of the other sensors, the
Euclidean distances that separate it from three of the anchors (its tripolar
coordinates with respect to those anchors, from which the sensor's position can
be easily calculated \cite{wolfram}).

Finding a sensor's Euclidean distance to an anchor is not simple, though.
Sometimes signal propagation is used for direct or indirect measurement
\cite{radar,spoton,cricket,robust_range,omnet,aps}, but there are approaches
that rely on no such techniques \cite{amorph,gps-less,dv_based}. The latter
include one of the most successful distributed approaches \cite{dv_based}, which
nonetheless suffers from increasing lack of accuracy as sparsity or irregularity
in sensor positioning become more pronounced. The algorithm of \cite{dv_based}
assumes, for each anchor $i$, that each edge on any shortest path to $i$ is
equivalent to a fixed Euclidean distance, which is estimated by $i$ in
communication with the other anchors and by simple proportionality can be used
by any node to infer its Euclidean distance to $i$. We believe that knowledge of
distance-related distributions has an important role to play in replacing this
assumption and perhaps dispelling the algorithm's difficulties in the less
favorable circumstances alluded to above.

We proceed in the following manner. In Section~2 we give some notation and
establish the overall approach to be followed when pursuing the analytical
characterization of distance-related distributions. Then in Sections~3 through
5 we present the mathematical analysis of the $d=1$ through $d=3$ cases. We
continue in Section~6 with computational results related to $d\ge 1$ and close
in Section~7 with some discussion and concluding remarks.

\section{Overall approach}

Let $i$ and $j$ be two distinct, randomly chosen nodes. For $d>0$ an integer and
$\delta\ge 0$ a real number, we use $P_\delta(d)$ to denote the probability,
conditioned on $\delta_{ij}=\delta$, that $d_{ij}=d$. Likewise, we use
$p_d(\delta)$ to denote the probability density, conditioned on $d_{ij}=d$,
associated with $\delta_{ij}=\delta$. These two quantities relate to each other
in the standard way of combining integer and continuous random variables
\cite{trivedi}.

If we assume that $P_\delta(d)$ is known for all applicable values of $d$ and
$\delta$, then it follows from Bayes' theorem that
\begin{equation}
p_d(\delta)
=\frac{P_\delta(d)p(\delta)}{P(d)},
\end{equation}
where $p(\delta)$ is the unconditional probability density associated with the
occurrence of an Euclidean distance of $\delta$ separating two randomly chosen
nodes and $P(d)$ is the unconditional probability that the distance between them
is $d$. Clearly, $P(d)=\int_{r=0}^{dR}P_r(d)p(r)\mathrm{d}r$, since $P_r(d)=0$
for $r>dR$. Moreover, $p(r)$ is proportional to the circumference of a
radius-$r$ circle, $2\pi r$, which yields
\begin{equation}
p_d(\delta)=\frac{P_\delta(d)\delta}{\int_{r=0}^{dR}P_r(d)r\mathrm{d}r}.
\label{eq:pfromP}
\end{equation}

In view of Equation (\ref{eq:pfromP}), our approach henceforth is to concentrate
on calculating $P_\delta(d)$ for all appropriate values of $d$ and $\delta$,
and then to use the equation to obtain $p_d(\delta)$. In order to calculate
$P_\delta(d)$, we fix two nodes $a$ and $b$ such that $\delta_{ab}=\delta$ and
proceed by analyzing how the two radius-$R$ circles (the one centered at $a$ and
the one at $b$) relate to each other. While doing so, we assume that the
two-dimensional region containing the graph has unit area, so that the area of
any of its sub-regions automatically gives the probability that it contains a
randomly chosen node. We assume further that all border effects can be safely
ignored (but see Section~6 for the computational setup that justifies this).

\section{The distance-1 and distance-2 cases}

The case of $d=1$ is straightforward, since $d_{ab}=d$ if and only if
$\delta\le R$. Consequently,
\begin{equation}
P_\delta(1)=\cases{
1,&if $\delta\le R$;\cr
0,&otherwise
}
\end{equation}
and, by Equation (\ref{eq:pfromP}),
\begin{equation}
p_1(\delta)=\cases{
2\delta/R^2,&if $\delta\le R$;\cr
0,&otherwise.
}
\end{equation}

For $d=2$, we have $d_{ab}=d$ if and only if $\delta >R$ and at least one node
$k$ exists, with $k\notin\{a,b\}$, such that $\delta_{ak}\le R$ and
$\delta_{bk}\le R$. The probability that this holds for a randomly chose $k$ is
given by the intersection area of the radius-$R$ circles centered at $a$ and
$b$, here denoted by $\rho_\delta$. From \cite{wolfram}, we have
\begin{equation}
\rho_\delta=\cases{
2R^2\cos^{-1}\left(\delta/2R\right)
-\delta\sqrt{R^2-\delta^2/4},&if $\delta\le 2R$;\cr
0,&otherwise.
}
\label{eq:intersection}
\end{equation}
Because any node that is not $a$ or $b$ may, independently, belong to such
intersection, we have
\begin{equation}
P_\delta(2)=\cases{
1-(1-\rho_\delta)^{n-2},&if $\delta>R$;\cr
0,&otherwise.
}
\end{equation}
As for $p_2(\delta)$, it is as given by Equation (\ref{eq:pfromP}), equaling $0$
if $\delta\le R$ or $\delta>2R$ (we remark that a closed-form expression is
obtainable also in this case, but it is too cumbersome and is for this reason
omitted).

\section{The distance-3 case: exact basis}

The $d=3$ case is substantially more complex than its predecessors in Section~3.
We begin by noting that $d_{ab}=d$ if and only if the following three conditions
hold:
\begin{itemize}
\item[C1.] $\delta>R$.
\item[C2.] No node $i$ exists such that both $\delta_{ai}\le R$ and
$\delta_{bi}\le R$.
\item[C3.] At least one node $k\notin\{a,b\}$ exists, and for this $k$ at least
one node $\ell\notin\{a,b,k\}$, such that $\delta_{ak}\le R$,
$\delta_{k\ell}\le R$, $\delta_{b\ell}\le R$, $\delta_{a\ell}>R$, and finally
$\delta_{bk}>R$.
\end{itemize}
For each fixed $k$ and $\ell$ in Condition~C3, these three conditions result
from the requirement that nodes $a$, $k$, $\ell$, and $b$, in this order,
constitute a shortest path from $a$ to $b$.

If we fix some node $k\notin\{a,b\}$ for which $\delta_{ak}\le R$ and
$\delta_{bk}>R$, the probability that Condition~C3 is satisfied by $k$ and a
randomly chosen $\ell$ is a function of intersection areas of circles that
varies from case to case, depending on the value of $\delta$. There are two
cases to be considered, as illustrated in Figure~\ref{fig:2or3circles}. In the
first case, illustrated in part (a) of the figure, $R<\delta\le 2R$ and node
$\ell$ is to be found in the intersection of the radius-$R$ circles centered at
$b$ and $k$, provided it is not also in the radius-$R$ circle centered at $a$.
The intersection area of interest results from computing the intersection area
of two circles (those centered at $b$ and $k$) and subtracting from it the
intersection area of three circles (those centered at $a$, $b$, and $k$). The
former of these intersection areas is given as in
Equation~(\ref{eq:intersection}), with $\delta_{bk}$ substituting for $\delta$;
as for the latter, closed-form expressions also exist, as given in
\cite{3circles}. The second case, shown in part (b) of
Figure~\ref{fig:2or3circles}, is that of $2R<\delta\le 3R$, and then the
intersection area of interest is the one of the circles centered at $b$ and $k$.
Regardless of which case it is, we use $\sigma_\delta^k$ to denote the resulting
area. Thus, the probability that at least one $\ell$ exists for fixed $k$ is
$1-(1-\sigma_\delta^k)^{n-3}$.

\begin{figure}[t]
\centering
\begin{tabular}{c}
\includegraphics[scale=0.65]{fig1a.eps}\\
{\small (a)}\\
\includegraphics[scale=0.65]{fig1b.eps}\\
{\small (b)}
\end{tabular}
\caption{Regions (shown in shades) whose areas yield the value of
$\sigma_\delta^k$ for $R<\delta\le 2R$ (a) and $2R<\delta\le 3R$ (b).}
\label{fig:2or3circles}
\end{figure}

Now let $P'_\delta(3)$ be the probability that a randomly chosen $k$ satisfies
Condition~C3. Let also $K_\delta$ be the region inside which such a node can be
found with nonzero probability. If $x_k$ and $y_k$ are the Cartesian coordinates
of node $k$, then each possible location of $k$ inside $K_\delta$ contributes to
$P'_\delta(3)$ the infinitesimal probability
$[1-(1-\sigma_\delta^k)^{n-3}]\mathrm{d}x_k\mathrm{d}y_k$. It follows that
\begin{equation}
P'_\delta(3)=\int_{k\in K_\delta}
[1-(1-\sigma_\delta^k)^{n-3}]\mathrm{d}x_k\mathrm{d}y_k.
\end{equation}

There are three possibilities for the region $K_\delta$, shown in parts (a)
through (c) of Figure~\ref{fig:Kdelta} as shaded regions, respectively for
$R<\delta\le R\sqrt{3}$, $R\sqrt{3}<\delta\le 2R$, and $2R<\delta\le 3R$. The
shaded region in part (a) is delimited by four radius-$R$ circles, the ones
centered at nodes $a$ (above and below) and $b$ (on the right) and the ones
centered at points $D$ and $E$ (on the left). As $\delta$ gets increased beyond
$R\sqrt{3}$---and, at the threshold, point $D$ becomes collinear with point $B$
and node $b$---we move into part (b) of the figure, where the shaded region is
now delimited on the left either by the radius-$R$ circles centered at $D$ and
$E$ or by the radius-$2R$ circle centered at $b$, depending on the point of
common tangent between each of the radius-$R$ circles and the radius-$2R$
circle. The next threshold leads $\delta$ beyond $2R$, and in part (c) of the
figure the shaded region is delimited on the left by the radius-$2R$ circle
centered at $b$, on the right by the radius-$R$ circle centered at $a$.

\begin{figure}[p]
\centering
\begin{tabular}{c}
\includegraphics[scale=0.65]{fig2a.eps}\\
{\small (a)}\\
\includegraphics[scale=0.65]{fig2b.eps}\\
{\small (b)}\\
\includegraphics[scale=0.65]{fig2c.eps}\\
{\small (c)}
\end{tabular}
\caption{Regions (shown in shades) where node $k$ can be found with nonzero
probability for $R<\delta\le R\sqrt{3}$ (a), $R\sqrt{3}<\delta\le 2R$ (b), and
$2R<\delta\le 3R$ (c).}
\label{fig:Kdelta}
\end{figure}

Figure~\ref{fig:Kdelta} is also useful in helping us obtain a more operational
version of the expression for $P'_\delta(3)$, to be used in Section~6. First we
establish a Cartesian coordinate system by placing its origin at node $a$ and
making the positive abscissa axis go through node $b$. In this system, the
shaded regions in all of parts (a) through (c) of the figure are symmetrical
with respect to the abscissa axis. If for each value of $x_k$ we let
$y_k^-(x_k)$ and $y_k^+(x_k)$ be, respectively, the minimum and maximum $y_k$
values in the upper half of the shaded region for the value of $\delta$ at hand,
then
\begin{equation}
P'_\delta(3)=2\int_{x_k=x_k^-}^{x_k^+}\int_{y_k=y_k^-(x_k)}^{y_k^+(x_k)}
[1-(1-\sigma_\delta^k)^{n-3}]\mathrm{d}x_k\mathrm{d}y_k,
\end{equation}
where $x_k^-$ and $x_k^+$ bound the possible values of $x_k$ for the given
$\delta$.

All pertinent values of $x_k^-$ and $x_k^+$, as well as of $y_k^-(x_k)$ and
$y_k^+(x_k)$, are given in Table~\ref{tab:limits}, where $\delta^-$ and
$\delta^+$ indicate, respectively, the lower and upper limit for $\delta$ in
each of the three possible cases. This table's entries make reference to the
abscissae of points $A$, $B$, $C$, and $D$ (respectively $x_A$, $x_B$, $x_C$,
and $x_D$) and to the ordinate of point $D$ ($y_D$). These are given in
Table~\ref{tab:nodes}.

\begin{table}[p]
\centering
{\renewcommand{\arraystretch}{1.5}
\caption{Cartesian coordinates delimiting the upper halves of shaded regions in
Figure~\ref{fig:Kdelta}.}
\label{tab:limits}
{\small\input{tab1}}
}
\end{table}

\begin{table}[p]
\centering
{\renewcommand{\arraystretch}{1.5}
\caption{Cartesian coordinates used in Table~\ref{tab:limits}.}
\label{tab:nodes}
{\small\input{tab2}}
}
\end{table}

\section{The distance-3 case: approximate extension}

Obtaining $P_\delta(3)$ from $P'_\delta(3)$ requires that we fulfill the
remaining requirements set by Conditions~C2 and C3 in Section~4. These are that
no node exists in the intersection of the radius-$R$ circles centered at $a$ and
$b$ and that at least one node $k$ exists with the properties given in
Condition~C3. While the probability of the former requirement is simply
$(1-\rho_\delta)^{n-2}$, expressing the probability of the latter demands that
we make a careful approximation to compensate for the lack of independence of
certain events with respect to one another.

For node $i\notin\{a,b\}$, let $\epsilon_i$ stand for the event that
Condition~C3 does not hold for $k=i$. Let also $Q_\delta(\epsilon_i)$ be the
probability of $\epsilon_i$ and $Q_\delta$ the joint probability of all $n-2$
events.  Clearly, $Q_\delta(\epsilon_i)=1-P'_\delta(3)$ for any $i$ and, for
$\delta>R$, $P_\delta(3)=(1-Q_\delta)(1-\rho_\delta)^{n-2}$. Therefore, if all
the $n-2$ events were independent of one another, we would have
\begin{equation}
Q_\delta=\prod_{i\notin\{a,b\}}Q_\delta(\epsilon_i)=[1-P'_\delta(3)]^{n-2}
\end{equation}
and, consequently,
\begin{equation}
P_\delta(3)=\cases{
\{1-[1-P'_\delta(3)]^{n-2}\}(1-\rho_\delta)^{n-2},&if $\delta>R$;\cr
0,&otherwise.
}
\end{equation}

However, once we know of a certain node $i$ that Condition~C3 does not hold for
it, immediately we reassess as less likely that the condition holds for nodes in
the Euclidean vicinity of $i$. The $n-2$ events introduced above are then not
unconditionally independent of one another, although we do expect whatever
degree of dependence there is to wane progressively as we move away from node
$i$.

We build on this intuition by postulating the existence of an integer $n'<n-2$
such that the independence of the $n'$ events not only holds but is also
sufficient to determine $P_\delta(3)$ as indicated above, provided the
corresponding $n'$ nodes are picked uniformly at random. But since this is
precisely the way in which, by assumption, sensors are positioned, it suffices
that any $n'$ nodes be selected, yielding
\begin{equation}
P_\delta(3)=\cases{
\{1-[1-P'_\delta(3)]^{n'}\}(1-\rho_\delta)^{n-2},&if $\delta>R$;\cr
0,&otherwise.
}
\end{equation}
Similarly to the previous cases, $p_3(\delta)$ is given by
Equation~(\ref{eq:pfromP}) and equals $0$ if $\delta\le R$ or $\delta>3R$.

It remains, of course, for the value of $n'$ to be discovered if our postulate
is to be validated. We have done this empirically, by means of computer
simulations, as discussed in Section~6.

\section{Computational results}

In this section we present simulation results and, for $d=1,2,3$, contrast them
with the analytic predictions of Sections~3 through 5. The latter are obtained
by numerical integration when a closed-form expression is not available (the
case of $d=3$ also requires simulations for finding $n'$; see below). For $d>3$,
we demonstrate that good approximations by Gaussians can be obtained.

We use $n=1\,000$ and a circular region of unit area, therefore of radius
$\sqrt{1/\pi}\approx 0.564$, for the placement of nodes. Node $a$ is always
placed at the circle's center, which has Cartesian coordinates $(0,0)$, and all
results refer to distances to $a$. Our choice for the value of $R$ depends on
the expected number of neighbors (or connectivity) of a node, which we denote by
$z$ and use as the main parameter. Since $z=\pi R^2n$ for large $n$, choosing
the value of $z$ immediately yields the value of $R$ to be used. We use $z=3\pi$
and $z=5\pi$, which yield, respectively, $R\approx 0.055$ and $R\approx 0.071$.
We note that both values of $z$ are significantly above the phase transition
that gives rise to the giant component, which happens at $z\approx 4.52$
\cite{rgg}.  In all our experiments, then, graphs are connected with high
probability.

For each value of $z$, each simulation result we present is an average over
$10^6$ independent trials. Each trial uses a matrix of accumulators having $n-1$
rows (one for each of the possible distance values) and $1\,000\sqrt{1/\pi}$
columns (one for each of the $0.001$-wide bins into which Euclidean distances
are compartmentalized). A trial consists of: placing $n-1$ nodes uniformly at
random in the circle; computing the Euclidean distance between each node and
node $a$; computing the distances between each node and node $a$ (this is done
with Dijkstra's algorithm \cite{cormen}); updating the accumulator that
corresponds to each node, given its two distances. At the end of each trial, its
contributions to $P_\delta(d)$ and $p_d(\delta)$ are computed, with
$d=1,2,\ldots,n-1$ and $\delta$ ranging through the middle points of all bins.
If $M$ is the matrix of accumulators, then these contributions are given,
respectively, by $M(d,\delta)/\sum_{d'}M(d',\delta)$ and
$M(d,\delta)/0.001\sum_{\delta'}M(d,\delta')$.

The case of $d=3$ requires two additional simulation procedures, one for
determining simulation data for $P'_\delta(3)$, the other to determine $n'$ for
use in obtaining analytic predictions for $P_\delta(3)$. The former of these
fixes node $b$ at coordinates $(\delta,0)$ and performs $10^7$ independent
trials. At each trial, two nodes are generated uniformly at random in the
circle. At the end of all trials, the desired probability is computed as the
fraction of trials that resulted in nodes $k$ and $\ell$ as in Section~4.

The simulation for the determination of $n'$ is conducted for $\delta=2R$ only,
whence $\rho_\delta=0$. This is the value of $\delta$ for which the results from
the simulation above for $P_\delta(3)$ and the analytic prediction for
$1-[1-P'_\delta(3)]^{n-2}$ differ the most (data not shown). Moreover, as we
will see shortly, the value of $n'$ we find using this value of $\delta$ is good
for all other values as well. The simulation is aimed at finding the value of
$Q_\delta$ and proceeds in $10^9$ independent trials. Each trial fixes node $b$
at $(\delta,0)$ and places the remaining $n-2$ nodes in the circle uniformly at
random. The fraction of trials resulting in no node qualifying as the node $k$
of Section~4 is the value of $Q_\delta$. We set $n'$ to be the $m<n-2$ that
minimizes $\vert Q_\delta-[1-P'_\delta(3)]^m\vert$, where $P'_\delta(3)$ refers
to the analytic prediction. Our results are $n'=779$ for $z=3\pi$, $n'=780$ for
$z=5\pi$.

Results for $d=1$ are shown in Figure~\ref{fig:d=1}, for $d=2$ in
Figure~\ref{fig:d=2}, for $d=3$ in Figures~\ref{fig:P'} and \ref{fig:d=3}, and
for $d>3$ in Figure~\ref{fig:d>3}. In all figures, both $P_\delta(d)$ and
$p_d(\delta)$ are plotted against $\delta$, since it seems better to visualize
what happens as one gets progressively farther from node $a$ in Euclidean terms.
For this reason, the plots for $P_\delta(d)$ do not constitute a probability
distribution for any fixed value of $d$.

\begin{figure}[p]
\centering
\begin{tabular}{c}
\includegraphics[scale=0.65]{fig3a.eps}\\
{\small (a)}\\
\includegraphics[scale=0.65]{fig3b.eps}\\
{\small (b)}
\end{tabular}
\caption{$P_\delta(1)$ (a) and $p_1(\delta)$ (b). Solid lines give the analytic
predictions.}
\label{fig:d=1}
\end{figure}

\begin{figure}[p]
\centering
\begin{tabular}{c}
\includegraphics[scale=0.65]{fig4a.eps}\\
{\small (a)}\\
\includegraphics[scale=0.65]{fig4b.eps}\\
{\small (b)}
\end{tabular}
\caption{$P_\delta(2)$ (a) and $p_2(\delta)$ (b). Solid lines give the analytic
predictions.}
\label{fig:d=2}
\end{figure}

\begin{figure}[t]
\centering
\begin{tabular}{c}
\includegraphics[scale=0.65]{fig5.eps}
\end{tabular}
\caption{$P'_\delta(3)$. Solid lines give the analytic predictions.}
\label{fig:P'}
\end{figure}

\begin{figure}[p]
\centering
\begin{tabular}{c}
\includegraphics[scale=0.65]{fig6a.eps}\\
{\small (a)}\\
\includegraphics[scale=0.65]{fig6b.eps}\\
{\small (b)}
\end{tabular}
\caption{$P_\delta(3)$ (a) and $p_3(\delta)$ (b). Solid lines give the analytic
predictions.}
\label{fig:d=3}
\end{figure}

\begin{figure}[p]
\centering
\begin{tabular}{c}
\includegraphics[scale=0.65]{fig7a}\\
{\small (a)}\\
\includegraphics[scale=0.65]{fig7b}\\
{\small (b)}
\end{tabular}
\caption{$P_\delta(d)$ and $p_d(\delta)$ for $d>3$, with $z=3\pi$ (a, b) and
$z=5\pi$ (c, d). Solid lines give the Gaussians that best fit some of the
$p_d(\delta)$ data, each of mean $\mu$ and standard deviation $\sigma$ as
indicated.}
\label{fig:d>3}
\end{figure}

\addtocounter{figure}{-1}
\begin{figure}[p]
\centering
\begin{tabular}{c}
\includegraphics[scale=0.65]{fig7c}\\
{\small (c)}\\
\includegraphics[scale=0.65]{fig7d}\\
{\small (d)}
\end{tabular}
\caption{Continued.}
\end{figure}

\section{Discussion and conclusion}

The results summarized in Figures~\ref{fig:d=1} through \ref{fig:d=3} reveal
excellent agreement between the analytic predictions we derived in Sections~3
through 5 and our simulation data. This holds not only for the simple cases of
$d=1$ and $d=2$, but also for the considerably more complex cases of
$P'_\delta(3)$ and $P_\delta(3)$. The latter, in particular, depends on the
empirically determined $n'$. In this respect, it is clear from
Figure~\ref{fig:d=3} that, even though $n'$ could have been calculated for a
greater assortment of $\delta$ values, doing it exclusively for $\delta=2R$
seems to have been sufficient.

Figure~\ref{fig:d>3} contemplates some of the $d>3$ cases, for which we derived
no analytic predictions. The values of $d$ that the figure covers in parts (a,
b) and (c, d), respectively for $z=3\pi$ and $z=5\pi$, are $4,\ldots,11$. Of
these, $d=11$ for $z=5\pi$ in part (d) typifies what happens for larger values
of $d$ as well (omitted for clarity), viz.\ probability densities sharply
concentrated at the border of the radius-$\sqrt{1/\pi}$ circle centered at node
$a$. Note that the same also occurs for $z=3\pi$, but owing to the smaller $R$
it only happens for larger values of $d$ (omitted from part (b), again for
clarity).

For $4\le d\le 11$ with $z=3\pi$, and $4\le d\le 9$ with $z=5\pi$,
Figures~\ref{fig:d>3}(b) and (d) also display Gaussian approximations of
$p_d(\delta)$. Parts (a) and (c) of the figure, in turn, contain the
corresponding simulation data only, and we remark that the absence of some
approximation computed from the Gaussians of part (b) or (d) is not a matter of
difficulty of principle. In fact, the counterpart of Equation~(\ref{eq:pfromP}),
obtained also from Bayes' theorem and such that
\begin{equation}
P_\delta(d)
=\frac{p_d(\delta)P(d)}{p(\delta)}
=\frac{p_d(\delta)P(d)}{\sum_{s=1}^{n-1}p_s(\delta)P(s)},
\end{equation}
can in principle be used with either those Gaussians or the concentrated
densities in place of $p_s(\delta)$ as appropriate for each $s$. What prevents
this, however, is that we lack a characterization of $P(s)$ that is not based on
simulation data only.

Still in regard to Figure~\ref{fig:d>3}, one possible interpretation of the good
fit by Gaussians of the simulation data for $p_d(\delta)$ comes from resorting
to the central limit theorem in its classical form \cite{trivedi}. In order to
do this, we view $\delta$ as valuing the random variable representing the
average Euclidean distance to node $a$ of all nodes that are $d$ edges apart
from $a$. The emergence of $p_d(\delta)$ as a Gaussian for $d>3$ (provided $d$
is small enough that the circle's border is not influential) may then indicate
that, for each value of $d$, the Euclidean distances of those nodes to node $a$
are independent, identically distributed random variables. While we know that
this does not hold for the smaller values of $d$ as a consequence of the
uniformly random positioning of the nodes in the circle (smaller Euclidean
distances to $a$ are less likely to occur for the same value of $d$), it would
appear that it begins to hold as $d$ is increased.

To summarize, we have considered a network of sensors placed uniformly at random
in a two-dimensional region and, for its representation as a random geometric
graph, have studied two distance-related distributions. One of them is the
probability distribution of distances between two randomly chosen nodes,
conditioned on the Euclidean distance between them. The other is the probability
density function associated with the Euclidean distance between two randomly
chosen nodes, given the distance between them. We have provided analytical
characterizations whenever possible, in the simplest cases as closed-form
expressions, and have also validated these predictions through simulations.

While further work related to additional analytical characterizations is worth
undertaking, as is the investigation of the three-dimensional case, we find that
the most promising tracks for future investigation are those that relate to
applications. In Section~1 we illustrated this possibility in the context of
sensor localization, for which it seems that understanding the distance-related
distributions we have studied has the potential to help in the discovery of
better distributed algorithms. Whether there will be success on this front
remains to be seen, as well as whether other applications will be found with the
potential to benefit from the results we have presented.

\subsection*{Acknowledgments}

The authors acknowledge partial support from CNPq, CAPES, and a FAPERJ BBP
grant.

redpajama_set_name: RedPajamaArXiv


text:
\section{INTRODUCTION} \label{sec:Intro}
If low-energy Supersymmetry (SUSY)~\cite{susy} is realized in nature,
the next generation of colliders, the Large Hadron Collider
(LHC) and the International Linear Collider
(ILC) are likely to copiously produce SUSY particles and
will allow for precise measurements of their properties. Once SUSY is
established experimentally, it is the main task to explore the unknown
mechanism of SUSY breaking (SSB). The most general way to do so is to
reconstruct the parameters of the general Minimal Supersymmetric
Standard Model (MSSM)~\cite{mssm} parameter space, which is
significantly more ambitious than a test of specific models of SUSY
breaking, like e.~g.~minimal supergravity (mSUGRA)~\cite{msugra} can
be tested against the data in a relative straight-forward manner (see
e.~g.~\cite{msugrafits}) due to their small number of parameters.

The general MSSM Lagrangian contains more than 100 new parameters
which are only related to each other through the unknown SSB
mechanism. Many of them (like complex phases) are already limited to
some extent by measurements, e.g. by the absence of neutron and
electron electric dipole moments~\cite{edm}. Furthermore universality
of first and second generation appears to be a reasonable
approximation while large Yukawa couplings in the third generation lead to
significant differences. Applying these constraints, the number of
free parameters is reduced to 19, including the Standard Model top
quark mass as a parameter to account for parametric uncertainties. It
is of high interest how well this 19-dimensional parameter space can
be restricted with future measurements from LHC and ILC.

The experimental collaborations ATLAS~\cite{atlas} and CMS~\cite{CMS}
at the LHC have performed detailed simulations of the possibilities to
extract mass information from their future data. At a future electron
positron collider for collision energies up to 1~TeV, the
International Linear Collider (ILC), the kinematically accessible part
of the superpartner spectrum can be studied in great detail due to
favorable background conditions and the well-known initial
state~\cite{ref:LHCILC}. The attempt of an evaluation of the MSSM
Lagrangian parameters and their associated errors is only useful if
the experimental errors of the future measurements are known and under
control. The by far best-studied SUSY scenario to serve as a basis for
such an evaluation is an mSUGRA-inspired benchmark scenario, the SPS1a
scenario~\cite{ref:SPS}.  For this scenario with a relatively light
superpartner spectrum, a wealth of experimental simulations exists and
has recently been compiled in the framework of the international
LHC/ILC study group~\cite{ref:LHCILC}. In this paper, we perform a
global fit of the 19 MSSM parameters to those expected measurements
augmented by possible measurements of topological cross-sections
(i.~e.~cross sections times branching fractions) that will be possible
at the ILC with polarized beams. Since a detailed experimental
simulation for some measurements is lacking, we estimate their
uncertainties conservatively from the predicted cross-sections and
transferring experimental efficiency from well-studied cases.

As theoretical basis for the global fit presented in this paper we use
the loop-level calculations as implemented into the program
SPheno~\cite{ref:SPheno}.  In SPheno, the masses and decay branching
fractions of the superpartners are calculated as well as production cross
sections in $\text{e}^+\text{e}^-$ collisions.

In this paper we present the result of a global fit of the MSSM
Lagrangian parameters at the electro-weak scale for a slightly
modified SPS1a benchmark scenario using the program
Fittino~\cite{ref:Fittino,ref:FittinoProgram}. Previous evaluations of
the errors of those parameters~\cite{msugrafits} did not attempt to
develop a strategy to extract those parameters from data without {\it
a priori} knowledge. Within Fittino, special attention is given
precisely to this task, i.~e.~to find the parameter set which is most
consistent with the data in a $\chi^2$ minimization before a careful
evaluation of errors and correlations is performed. A similar
code~\cite{sfitter} exploits a different strategy for this task.

This paper is organized as follows. In Section~\ref{sec:Fit} we
shortly describe the approach used in Fittino. In
Subsection~\ref{sec:Inputs}, the input observables are explained. The
results of the fit and the error evaluation method are summarized in
Subsection~\ref{sec:Results} and conclusions are drawn in
Section~\ref{sec:conclusions}.

\section{SPS1a' FIT} \label{sec:Fit}
From the numerous fitting options provided by Fittino we have chosen
the following fit procedure to extract the low-energy Lagrangian
parameters. First, in order to be independent of human bias, start
values for the parameters are calculated using tree-level relations
between parameters and a few
observables~\cite{ref:FittinoProgram}. For fits with many parameters
these values are not good enough to allow a fitting tool like
MINUIT~\cite{ref:Minuit} to find the global minimum due to the amount
of loop-level induced cross-dependencies between the individual
sectors of the MSSM. Therefore, in a second step, the parameters are
refined, using a simulated annealing
approach~\cite{ref:FittinoProgram,SimAnn,SimAnn2}. As a result, the
parameter values are close to the global minimum so that a global fit
using MINUIT can find the exact minimum in a third step. To determine
the parameter uncertainties and correlations many individual fits with input
values randomly smeared within their uncertainty range are carried
out. The parameter uncertainties and the correlation
matrix are derived from the spread of the fitted parameter
values.

\subsection{Input Observables from LHC and ILC} \label{sec:Inputs}
A number of anticipated LHC and ILC measurements serve as input
observables to the fit.  For the ILC, running at center-of-mass
energies of 400 GeV, 500 GeV and 1 TeV is considered with 80~\%
electron and 60~\% positron polarization.  The predicted values of the
observables are calculated using the following prescriptions:
\begin{itemize}
\item {\bf Masses:} \newline The experimental uncertainties of the
  mass measurements are taken from \cite{ref:LHCILC}.
\item {\bf Cross-sections:} \newline Only $\text{e}^+ \text{e}^-$
  cross-sections are included in the fit. However, the measurement of
  absolute cross-sections is impossible for many channels, in which
  only a fraction of the final states can be reconstructed. Therefore,
  absolute cross-section measurements are only used for the
  Higgs-strahlung production of the light Higgs boson, which is
  studied in detail in~\cite{ref:HiggsBr}.
\item {\bf Cross-sections times branching fractions:} \newline Since
  no comprehensive study of the precision of cross-section times
  branching fraction measurements is available, the uncertainty is
  assumed to be the error of a counting experiment with the following
  assumptions:
  \begin{itemize}
  \item The selection efficiency amounts to 50~\%.
  \item 80~\% polarization of the electron beam and 60~\%
    polarization of the positron can be achieved.
  \item 500~fb$^{-1}$ per center-of-mass energy and polarization
    is collected.
  \item The relative precision is not allowed to be better than 1~\%
    and the absolute accuracy is at most 0.1~fb to account for
    systematic uncertainties.
  \end{itemize}
  All production processes and decays of SUSY particles and Higgs
  bosons are used which have a cross-section times branching ratio
  value of more than 1~fb in one of the $\text{e}^+ \text{e}^-$
  polarization states LL, RR, LR, RL at $\sqrt{s}=400$ and $500$~GeV
  and LR or RL at $\sqrt{s}=1000$~GeV.
\item {\bf Branching fractions:} \newline The three largest branching
  fractions of the lightest Higgs boson are included. The
  uncertainties on the Higgs branching fractions are taken
  from~\cite{ref:HiggsBr}.
\item {\bf Standard Model parameters:} \newline The present
  uncertainties for $m_{\text{W}}$ and $m_{\text{Z}}$ are used as
  conservative estimates. The experimental uncertainty for the top
  mass $m_{\text{t}}$ is assumed to be 50 MeV~\cite{ref:HiggsBr}.
\end{itemize}

\begin{table}[t]
\caption{The Fittino fit result for the SPS1a' motivated scenario.
  The left column shows the values predicted by SPheno version 2.2.2
  for this scenario, the second column exhibits the central values of
  the fit with unsmeared input observables, the third column displays
  the parameter uncertainties for the fit with experimental
  uncertainties only.  The last column shows the parameter
  uncertainties for the fit with experimental and theoretical
  uncertainties.\vspace{2mm}}
\label{tab:UnsmearedFitResultSPS1aPrime}
\begin{center}
\begin{tabular}{l e{7} e{7} e{7} e{7}}
\hline
\multicolumn{1}{c}{Parameter} & \multicolumn{1}{c}{``True'' value} &
\multicolumn{1}{c}{Fit value} & \multicolumn{1}{r}{Uncertainty} &
\multicolumn{1}{r}{Uncertainty} \\
 & & & \multicolumn{1}{c}{(exp.)} & \multicolumn{1}{c}{(exp.+theor.)}\\
 \hline
$\tan\beta$          &  10.00                &  10.00                 & 0.11 & 0.15 \\
$\mu$                &  400.4  \;\text{GeV}  &  400.4   \;\text{GeV}  & 1.2\;\text{GeV} & 1.3\;\text{GeV} \\
$X_{\tau}$           &  -4449.  \;\text{GeV}  & -4449.    \;\text{GeV}  & 20.\;\text{GeV} & 30. \;\text{GeV} \\
$M_{\tilde{e}_R}$    &  115.60  \;\text{GeV} &  115.60 \;\text{GeV}   &  0.27 \;\text{GeV} & 0.50   \;\text{GeV} \\
$M_{\tilde{\tau}_R}$ &  109.89  \;\text{GeV} &  109.89   \;\text{GeV} &  0.41  \;\text{GeV} & 0.60    \;\text{GeV} \\
$M_{\tilde{e}_L}$    &  181.30  \;\text{GeV} &  181.30   \;\text{GeV} &  0.10  \;\text{GeV} & 0.12    \;\text{GeV} \\
$M_{\tilde{\tau}_L}$ &  179.54  \;\text{GeV} &  179.54   \;\text{GeV} &  0.14  \;\text{GeV} & 0.19    \;\text{GeV} \\
$X_{\text{t}}$     &  -565.7  \;\text{GeV} & -565.7    \;\text{GeV} & 3.1   \;\text{GeV} & 15.4   \;\text{GeV} \\
$X_{\text{b}}$  &  -4935. \;\text{GeV} & -4935.   \;\text{GeV} & 1284. \;\text{GeV} & 1825.  \;\text{GeV} \\
$M_{\tilde{q}_R}$    &  503.   \;\text{GeV} &   503.   \;\text{GeV} & 24.  \;\text{GeV} &  27.  \;\text{GeV} \\
$M_{\tilde{b}_R}$    &  497.   \;\text{GeV} &  497.    \;\text{GeV} & 8.  \;\text{GeV} &  15.  \;\text{GeV} \\
$M_{\tilde{t}_R}$    &  380.9   \;\text{GeV} &  380.9    \;\text{GeV} & 2.5  \;\text{GeV} &  3.9    \;\text{GeV} \\
$M_{\tilde{q}_L}$    &  523.   \;\text{GeV} &  523.    \;\text{GeV} &  10.  \;\text{GeV} & 15.    \;\text{GeV} \\
$M_{\tilde{t}_L}$    &  467.7   \;\text{GeV} &  467.7    \;\text{GeV} & 3.1 \;\text{GeV} & 5.1    \;\text{GeV} \\
$M_1$                &  103.27 \;\text{GeV} &  103.27  \;\text{GeV} & 0.06  \;\text{GeV} & 0.14   \;\text{GeV} \\
$M_2$                &  193.45 \;\text{GeV} &  193.45  \;\text{GeV} & 0.10  \;\text{GeV} & 0.15  \;\text{GeV} \\
$M_3$                &  569.  \;\text{GeV} &  569.    \;\text{GeV} &  7.  \;\text{GeV} &   7.   \;\text{GeV} \\
$m_{\text{A}_{\text{run}}}$ & 312.0 \;\text{GeV} & 311.9 \;\text{GeV} & 4.6  \;\text{GeV} & 6.9   \;\text{GeV} \\
$m_{\text{t}}$     &  178.00  \;\text{GeV} &  178.00   \;\text{GeV} &  0.050  \;\text{GeV} & 0.108    \;\text{GeV} \\
\hline
\multicolumn{5}{c}{$\chi^2$ for unsmeared observables: $5.3\times 10^{-5}$}\\
\hline
\end{tabular}
\end{center}
\end{table}

In order to check the influence of theoretical uncertainties on the
fit results, the fit has been performed twice, once with experimental
uncertainties only and a second time with experimental and theoretical
uncertainties. Theoretical uncertainties for the mass predictions are
scale uncertainties taken from~\cite{ref:SPA}. The experimental and
theoretical contributions are added in quadrature. The assumed
experimental and theoretical uncertainties for the masses are given in
\cite{ref:LHCILC}. For the cross-section and cross-section times
branching fraction measurements, the smallest allowed relative
precision is raised to 2 \% for the fit including theoretical
uncertainties (as opposed to 1 \% for the fit without theoretical
uncertainties). The full list of observables used in the fit and their
uncertainties can be obtained from~\cite{ref:Fittino}.

\subsection{Fit Results} \label{sec:Results}
The input observables described in Section~\ref{sec:Inputs} are used
to determine the SUSY Lagrangian parameters in a global fit under the
assumptions mentioned in Section~\ref{sec:Intro}. In total 18 SUSY
parameters remain to be fitted. In addition to those, the top mass
$m_{\text{t}}$ is fitted, since it strongly influences parts of the
MSSM observables.  Thus 19 Lagrangian parameters are simultaneously
determined in this fit.

As shown in Table~\ref{tab:UnsmearedFitResultSPS1aPrime} all
parameters are perfectly reconstructed at their input values. Due to
the fact that the input observables are unsmeared in this fit, the
final $\chi^2$ is close to zero at $\chi^2 = 5.3 \times 10^{-5}$.

Most parameters are reconstructed to a precision better than or around
1~\%. For the U(1) gaugino mass parameter $M_1$, an accurary below the
per mil level is achieved for the fit including experimental
uncertainties only.  However, its precision is worse by a factor of 2
once theoretical uncertainties are included. The uncertainties on
other parameters increase by up to a factor of 5 if theoretical
uncertainties are taken into account, such as $X_{\text{t}}$. In order
to fully benefit from the precision data provided in particular by the
ILC, work must be invested to reduce uncertainties of theoretical
predications.

After a successful convergence of the simulated annealing algorithm to
the input parameter values, the fit uncertainties are evaluated by
carrying out many individual fits with input values randomly smeared
within their uncertainty range using a Gaussian probability
density. The complete covariance matrix and the correlation matrix are
derived from the spread of the fitted parameter values.  For the fit
without theoretical uncertainties 1002 such fits were made. The
corresponding number for the case including theory uncertainties
amounts to 993. For each of those fits, the parameter set is
minimized. For a large and complex parameter space with large
correlations among the parameters this method has turned out to be
more robust than a MINOS error analysis.  An example of the outcome of
this procedure for the fit (with experimental uncertainties only) is
shown in Figure~\ref{pulldistributions} for the parameters $\tan\beta$
and $M_1$ for the 1002 independent fits.

\begin{figure*}[t]
\centering{
\includegraphics[width=0.45\textwidth]{slac-pub-11295-1.eps}
\includegraphics[width=0.45\textwidth]{slac-pub-11295-2.eps}
}
\caption{Examples for the parameter distributions from the toy
fits. No theoretical uncertainties are included. The agreement of the
parameter distributions with the Gaussian hypothesis is very good.}
\label{pulldistributions}
\end{figure*}

\begin{figure*}[t]
\centering
\includegraphics[width=0.5\textwidth]{slac-pub-11295-3.eps}
\caption{The $\chi^2$ distribution obtained from 1002 toy fits with
smeared input distributions. The agreement with the expectation of a
$\chi^2$ distribution with a mean of $126\pm 0.64$ is very good.}
\label{chisqdistribution}
\end{figure*}

Figure~\ref{chisqdistribution} shows the distribution of the $\chi^2$ values
for these 1002 independent fits for the case without theoretical
uncertainties.  The mean $\chi^2$ obtained from a fit of the $\chi^2$
distribution to the observed distribution is 124.6, in agreement with
the expectation of 126~$\pm$~0.64. This shows that the fits converge
well to the true minimum of the $\chi^2$ surface for the smeared
observables, implying that the uncertainties extracted from the toy
fit value distributions are correct.

\subsection{Relative Importance of Observables}
In order to determine automatically the most important input
observables influencing the precision with which a certain parameter,
the contribution of each observable to $\Delta \chi^2 =
\chi^2_{1\sigma} - \chi^2_{\text{min}}$ is calculated for the
parameter of interest.  $\chi^2_{\text{min}}$ is the $\chi^2$ value at
the minimum of the fit and $\chi^2_{1\sigma}$ is the corresponding
value if the parameter is varied within $\pm 1 \sigma$.  A large
$\Delta \chi^2$ indicates a strong correlation of the parameter with
other parameters. A large contribution $\Delta \chi^2_i$ of an
observable $i$ to the total $\Delta \chi^2$ with respect to the
contributions by the other observables means that the parameter is
mainly determined by the measurement of observable $i$.

\begin{table}[t]
\caption{Individual $\Delta\chi^2$ contributions in the Fittino
    SPS1a' motivated fit. In the first column the parameter and its fitted value
    is shown.  The second column shows the total $\Delta\chi^2$ for
    this parameter, if the parameter is varied by $1\,\sigma$. The
    third and fourth column show the parameters and their contribution
    to the $\Delta\chi^2$ of the parameter in per cent.\vspace{2mm}}\label{tab:IndividualChi2}
\begin{center}
\begin{tabular}{l e{5} l e{7}}
\hline
    Parameter & \multicolumn{1}{c}{Total $\Delta\chi^2$} & Observable & \multicolumn{1}{c}{Contribution to the}  \\
    Value     &                                          &            & \multicolumn{1}{c}{$\Delta\chi^2$ in \%} \\
\hline
\hline
    $\tan\beta$              &  5.0 & $\sigma(\text{e}^-_L\text{e}^+_R
                                       \rightarrow\mathrm{H}^{\pm}\mathrm{H}^{\mp}
                                       \rightarrow\mathrm{t}\bar{\mathrm{b}}
                                       \bar{\mathrm{t}}\mathrm{b})~1$~TeV                 & 31.1 \\
    $10.00 \pm 0.11$         &      & $\sigma(\text{e}^-_L\text{e}^+_R
                                       \rightarrow\mathrm{H}\mathrm{A}
                                       \rightarrow\mathrm{b}\bar{\mathrm{b}}
                                       \mathrm{b}\bar{\mathrm{b}}~1)$~TeV               & 9.61 \\
                             &      & $m_{\mathrm{h}}$         & 8.12 \\
    \hline

    $\mu$                    & 15.2 & $\sigma(\mathrm{e}^-_L\mathrm{e}^+_R
                                       \rightarrow\tilde{\chi}^+_1\tilde{\chi}^-_1
                                       \rightarrow\bar{\nu}_{\tau}\tilde{\chi}^0_1{\tau}^{+}\,
                                       \nu_{\tau}\tilde{\chi}^0_1{\tau}^{-})~400$~GeV   & 14.5 \\
    $400.39 \pm 1.18$~GeV    &      & $\sigma(\mathrm{e}^-_L\mathrm{e}^+_R
                                       \rightarrow\tilde{\chi}^+_1\tilde{\chi}^-_1
                                       \rightarrow\bar{\nu}_{\tau}\tilde{\chi}^0_1{\tau}^{+}\,
                                       \nu_{\tau}\tilde{\chi}^0_1{\tau}^{-})~500$~GeV             & 7.49 \\
                             &      & $\sigma(\mathrm{e}^-_R\mathrm{e}^+_R
                                       \rightarrow\tilde{\chi}^+_1\tilde{\chi}^-_1
                                       \rightarrow\bar{\nu}_{\tau}\tilde{\chi}^0_1{\tau}^{-}\,
                                       \nu_{\tau}\tilde{\chi}^0_1{\tau}^{+})~500$~GeV    &  6.71 \\
    \hline
    $M_{\tilde{\mathrm{q}}_R}$ & 36.6 & $m_{\tilde{\mathrm{e}}_R}$ & 53.6 \\
    $501.6 \pm 23.6$~GeV     &      & $m_{\tilde{\mu}_R}$      &  3.34 \\
                             &      & $\sigma(\mathrm{e}^-_L\mathrm{e}^+_R
                                       \rightarrow\tilde{\chi}^+_1\tilde{\chi}^-_1
                                       \rightarrow\bar{\nu}_{\tau}\tilde{\chi}^0_1{\tau}^{-}\,
                                       \nu_{\tau}\tilde{\chi}^0_1{\tau}^{+})~1$~TeV &  3.25 \\
    \hline
\end{tabular}
\end{center}
\end{table}

Generally the precision of most parameters is determined by ILC data,
as shown as an example for $\tan\beta$ and $\mu$ in
Table~\ref{tab:IndividualChi2}. As a surprise it can be seen that even
for the uncertainty of the squark mass parameters the most
constraining input is provided by ILC slepton mass measurements. It
shows that the small contribution of the squark mass parameter
$M_{\tilde{\text{q}}_{R}}$ to the slepton mass $m_{\tilde{e}_R}$ has a
larger effect compared to the good precision of the measurement of
$m_{\tilde{e}_R}$, than the large contribution of the squark mass
parameter to the much less precisely measured squark mass.

\subsection{LHC and ILC Synergy}
In order to test the dependence of the fit results on the synergy of
LHC and ILC, a fit has been performed with all ILC-influenced
observables removed. While the LHC observables may allow to determine
a further reduced set of parameters (e.g. mSUGRA, if a simplified
model described by fewer parameters still fits the data), the complete
set of parameters studied here is not constrained by the LHC mass
observables which currently enter the fit. For most parameters, the
resulting parameter uncertainties exceed the LHC+ILC uncertainties by
orders of magnitude. Only for 3 of the parameters studied here, namely
$M_{\tilde{\mathrm{q}}_L}$, $M_{\tilde{\mathrm{q}}_R}$ and $M_3$, the
uncertainty is in the same order of magnitude. This exemplifies the
synergy effects of the LHC and the ILC.

\section{CONCLUSIONS} \label{sec:conclusions}
The Lagrangian parameters of the MSSM assuming universality for the
first and second generation and real parameters but without
assumptions on the SUSY breaking mechanism are correctly reconstructed
without usage of {\it a priori} information. This has been achieved using
precision measurements at the LHC and ILC as input to a global fit
exploiting the techniques implemented in the program Fittino.

Most of the Lagrangian parameters can be determined to a precision
around the percent level. For some parameters an accuracy of better
than 1~per mil is achievable if theoretical uncertainties are
neglected. A fit including them revealed that theory uncertainties can
significantly deteriorate the precision of the Lagrangian parameter
determination. More work is needed to improve the accuracy of
theoretical predictions in order to fully benefit from the
experimental precision. The only parameter which cannot be strongly
constrained by the observables used in the presented fit is
$X_{\text{b}}$ resulting in a precision of only 30 \% for this
parameter.

\begin{acknowledgments}
The authors are grateful to Gudrid Moortgat-Pick, Werner Porod, Sven
Heinemeyer and the whole SPA working group for very fruitful
discussions.
\end{acknowledgments}


text:
\section{Introduction}
\label{subgraph}
Over the last decades, the availability of network data on typically
very large scales has created the impetus for the development of new theories and methods for
modeling and describing the properties of large networks. The introduction of exponential random graphs has
aided in this pursuit, as they are able to capture a wide variety of
common network tendencies by representing a complex global
structure through a set of tractable local features. From the point of view of extremal combinatorics and statistical mechanics, investigations have been focused on the variational principle of the limiting normalization constant, concentration of the limiting probability distribution, phase transitions, and asymptotic structures. See for example Chatterjee and Varadhan \cite{CV2011}, Chatterjee and Diaconis \cite{CD2013}, Radin and Yin \cite{RY}, Lubetzky and Zhao \cite{LZ1, LZ2}, Radin and Sadun \cite{RS1, RS2}, Radin et al. \cite{RRS}, Kenyon et al. \cite{KRRS1}, Yin \cite{Yin1}, Kenyon and Yin \cite{KY}, Aristoff and Zhu \cite{AZ2}, and Chatterjee and Dembo \cite{CD2}. The main techniques used in these papers are variants of statistical physics, but the elegant theory of graph limits as developed by Lov\'{a}sz and coauthors (V.T. S\'{o}s, B. Szegedy, C. Borgs, J. Chayes, K. Vesztergombi, ...) \cite{BCLSV1, BCLSV2, BCLSV3, Lov, LS}, also plays an important role in the interdisciplinary inquiry. Building on earlier work of Aldous \cite{Aldous1} and Hoover \cite{Hoover}, the graph limit theory connects sequences of graphs to a unified graphon space equipped with a cut metric. Though the theory itself is tailored to dense graphs, parallel theories for sparse graphs are likewise emerging. See Benjamini and Schramm \cite{BS}, Aldous and Steele \cite{AS}, Aldous and Lyons \cite{AL}, and Lyons \cite{Lyons} where the notion of local weak convergence is discussed and the recent works of Borgs et al. \cite{BCCZ1, BCCZ2} that are making progress towards enriching the existing $L^\infty$ theory of dense graph limits by developing a limiting object for sparse graph sequences based on $L^p$ graphons.

In this paper, we study the ``standard'' family of exponential graph models in the asymptotic regime as the exponential parameters approach infinity. As the model name indicates, we are associating exponential weights to graphical ensembles. For each $n$ let ${\mathcal G}_n$ be the ensemble of simple graphs on $n$ vertices and let ${\mathcal G} = \cup_n \mathcal{G}_n$ be the collection on all simple graphs. The exponential weights are defined in terms of  subgraph densities. For any $H \in {\mathcal G}$ the homomorphism density of $H $ in a graph $G \in {\mathcal G}$ is defined as the probability that a random map on the vertex set of $H$ into the vertex set of $G$, $V(H) \to V(G)$ is edge preserving.
We write the homomorphism density as
\begin{equation}
\label{t} t(H, G)=\frac{|\text{hom}(H,G)|}{|V(G)|^{|V(H)|}}.
\end{equation}
Now, for all $n$ we define a probability distribution on ${\mathcal G}_n$ in terms of the homomorphism density. Let $H_1,..,H_d$ be a given selection of simple graphs, where $H_1 $ is an edge: $H_1 = K_2$.  Let $T: {\mathcal G}_n \to [0,1]^d$, where the components of $T$ are homomorphism densities $\{T(G_n)\}_i=t(H_i, G_n)$. Given $\beta \in {\mathbb R}^{d}$, define the functional on $G_n \in \mathcal{G}_n$:
\begin{equation}\label{ham}
 T^\beta(G_n) := \beta \cdot T(G_n) =  \sum_{i = 1}^{d} \beta_i t(H_i, G_n) ,
\end{equation}
and weight $G_n \in {\mathcal G}_n$ by $e^{n^2 T^\beta(G_n)}$. The normalization constant for the ensemble ${\mathcal G}_n$ is then given by the partition function,
\begin{equation}\label{ptfn}
  Z^\beta_n = \sum_{G_n \in \mathcal{G}_n} \exp \left( n^2 T^{\beta} (G_n)\right).
\end{equation}
The terminology for the partition function is borrowed from thermodynamics. In this context, $|\beta|$ is the inverse temperature.
Renormalizing $T^\beta$ in (\ref{ham}) and fixing $\beta/|\beta|$, one obtains a Hamiltonian $H(G_n)  = - \frac{n^2}{|\beta|} T^{\beta}(G_n)$.   As $n\to \infty$ the major contribution to the partition function  concentrates around the thermal states (labeled by $f^\beta$),  in a manner analogous to standard thermodynamic models (see Section \ref{ierg}). These thermal states are well understood in the large temperature regime $|\beta|\rightarrow 0$. For any selection of subgraphs $H_1,..,H_d$ and $|\beta|$ sufficiently small, the associated thermal state lies in the replica symmetric phase \cite{CD2013}, i.e., $f^\beta \equiv \rho$. On the other hand, the ground states, which are defined as the limit of the thermal states $f^\beta$ as $|\beta | \to \infty$, are not so simple. In some cases the ground state is known to be in the replica symmetric phase \cite{CD2013}, while in other cases the ground state concentrates around a simple graph in ${\mathcal G}$ \cite{YRF}.

Our motivation for this paper comes from the edge-triangle model, obtained by setting $d = 2$ and $H_2 $ a triangle: $H_2 = K_3$. It was shown in \cite{YRF} that there are countably many critical directions of $\beta$, along which the ground state of the model is chosen from finitely many simple graphs with some unknown distribution, and our goal is to develop a mechanism that determines which of these simple graphs is the proper ground state.

\subsection{Graphon topology}
\label{bkgd}
The thermal states $f^\beta$ belong to the space of graph functions ``graphons" which may be understood as generalizations of graphs. The set of graphs $\mathcal{G}$ may be embedded into the space of graphons which consist of symmetric measurable functions from $[0,1]^2$ into $[0,1]$,
\begin{equation}
\mathcal{W} = \{f:[0,1]^2 \to [0,1] \text{ and } f \text{ symmetric}\}.
\end{equation}
For any $n $ and a graph $G_n \in \mathcal{G}_n$, the  graphon representation is the function
\begin{equation}
\label{graphon}
f^{G_n}(x,y)=\left\{%
\begin{array}{ll}
    1, & \hbox{if $(\lceil nx \rceil, \lceil ny \rceil)$ is an edge in $G_n$,} \\
    0, & \hbox{otherwise,} \\
\end{array}%
\right. \quad (x,y) \in [0,1]^2
\end{equation}
where the interval $[0,1] $ may be intuitively thought of as a `continuum' of graph vertices. The distance between graphons is given in terms of the ``cut distance'', defined for $f,h \in {\mathcal W}$, as
\begin{equation}
\label{cut}
d_{\Box}(f, h)=\sup_{S, T \subseteq [0,1]}\left|\int_{S\times
T}\left(f(x, y)-h(x, y)\right)dx\,dy\right|.
\end{equation}
However, a nontrivial difficulty arrises from the arbitrary labeling of vertices as they are embedded in $[0,1]$. Thus we introduce the equivalence $f \sim f_\sigma$ where $f_\sigma(x,y)  =   f(\sigma x, \sigma y)$, where $\sigma: [0,1] \to [0,1]$ is any measure preserving bijection. We write the quotient space of  graphons under the equivalence $\sim$ as $\widetilde {\mathcal W}$, and the equivalance class under $\sim$ of $f\in{\mathcal W}$ as $\widetilde f$. Incorporating the equivalence relation $\sim$ yields a distance
\begin{equation}
\label{cut2}
\delta_{\Box}(\tilde{f},
\tilde{h})=\inf_{\sigma_1, \sigma_2}d_{\Box}(f_{\sigma_1},
h_{\sigma_2}),
\end{equation}
 where the infimum ranges over all measure
preserving bijections $\sigma_1$ and $\sigma_2$, making
$(\widetilde{\mathcal W}, \delta_{\Box})$ a compact metric space (see Section 9.3 of Lov\'{a}sz \cite{Lov}).
With some abuse of notation we also refer to $\delta_{\Box}$ as
the ``cut distance''.

All graphons arise as the limit of some sequence of graphs. Given a graphon $f \in {\mathcal W}$ one may construct a graph $G_n = G(n, f)$ by selecting iid points $x_1,..,x_n$ uniformly from $[0,1]$ which represent the vertices of $G_n$, and then connect vertices $i,j$ with probability $f(x_i, x_j)$. In this context the expected subgraph density is given by the graphon homomorphism density
\begin{equation}
\label{tt} t(H, f)=\int_{[0,1]^k}\prod_{\{i,j\}\in E(H)}f(x_i, x_j)dx_1\cdots dx_k,
\end{equation}
which generalizes (\ref{t}) and is continuous in the metric $\delta_{\Box}$.
  Indeed, for any graph $H$, the subgraph homomorphism density of random graphs selected by the above construction converges almost surely to the graphon homomorphism density,
\begin{equation}
\lim_{n\rightarrow \infty} t(H, G(n,f))  = t(H,f).
\end{equation}

\subsection{Exponential random graphs}\label{ierg}
 As discussed above, we will define  measures on ${\mathcal G}_n$ in terms of subgraph densities. For $G_n \in {\mathcal G}_n$ define the probability
\begin{equation}
\label{pmf2}
{\mathbb P}_n^\beta(G_n)=\exp\left(n^2(T^\beta(G_n)) - \psi_n^\beta\right),
\end{equation}
 where we have introduced the normalization   constant (free energy density),
\begin{equation}
\label{psi2} \psi_n^\beta=\frac{1}{n^2}\log Z_n^\beta.
\end{equation}
By replacing the subgraph homomorphism density (\ref{t}) with graphon homomorphism density (\ref{tt}) in (\ref{ham}), $T^\beta$ extends naturally to $(\widetilde{\mathcal W},\delta_{\Box})$.
Since $T^\beta$ is continuous and bounded on the compact set $\widetilde {\mathcal W}$, there is a nonempty compact subset ${\mathcal K}^\beta$ of $\widetilde {\mathcal W}$ so that $T^\beta$ is maximized on ${\mathcal K}^\beta$. Take $I:[0,1] \to {\mathbb R} $ as
\begin{equation}
\label{I}
I(u) =  \frac12 u\log  u + \frac12(1 - u) \log  (1 - u ),
\end{equation}
and then extend the domain of $I$ to $\widetilde{\mathcal W}$ by
\begin{equation}
I(\tilde{f}) = \int_{[0,1]^2} I (f(x,y)) dxdy,
\end{equation}
where $f$ is any representative element of $\tilde{f}$.
It follows from Lemma 2.1 in
   Chaterjee and Varadhan \cite{CV2011} that $\psi^\beta(\widetilde f):= T^\beta(\widetilde f) - I(\widetilde f)$
   is well defined on $\widetilde {\mathcal W}$ and upper semi-continuous under the cut metric $\delta_\Box$. Let ${\mathcal F}^\beta$ be the subset of $\widetilde {\mathcal W}$ where $\psi^\beta$ is maximized. Then like ${\mathcal K}^\beta$, ${\mathcal F}^\beta$ is a nonempty compact subset of $\widetilde {\mathcal W}$.

For the purpose of this paper, two theorems from Chatterjee and Diaconis \cite{CD2013} (both based on the large deviation result in \cite{CV2011}) merit some
special attention. Together they connect the occurrence of a phase
transition in the exponential random graph model with the solution of a certain
maximization problem. The first theorem (Theorem 3.1 in \cite{CD2013}) states that for any $\beta$, the limiting normalization constant $\lim_{n\to \infty} \psi_n^\beta$ of the exponential random graph always exists and is equal to $\psi^\beta$. The second theorem (Theorem 3.2 in \cite{CD2013}) states that in the large $n$ limit, the quotient image
$\tilde{f}^{G_n}$ of a random graph $G_n$ drawn from (\ref{pmf2})
must lie close to ${\mathcal F}^\beta$ with high probability,
\begin{equation}
\delta_{\Box}(\tilde{f}^{G_n}, {\mathcal F}^\beta)\to 0\hbox{ in
probability as } n\to \infty.
\end{equation}
Since the limiting normalization constant $\psi^\beta$ is the generating function for the limiting
expectations of other random variables on the graph space such as expectations and
correlations of homomorphism densities, a phase transition occurs when $\psi^\beta$ is non-analytic
or when ${\mathcal F}^\beta$ is not a singleton set. Although it is difficult to evaluate $\psi^\beta$ and determine the maximizing set ${\mathcal F}^\beta$ for most $\beta$, we will derive an efficient method to approximate ${\mathcal F}^\beta$ and estimate $\psi^\beta$ for $\beta$ sufficiently far from the origin.

\section{The approximation scheme}
\label{approximation}
Take a finite simple graph $H$ with vertex set $V(H)$ and edge set $E(H)$. Consider a graphon $f \in {\mathcal W}$. For
each $(r, s) \in E(H)$ and each pair of points $x_r, x_s \in [0, 1]$, define
\begin{equation}
\Delta_{H, r, s}f(x_r, x_s):=\int_{[0, 1]^{\left|V(H)\texttt{\char92}\{r, s\}\right|}}\prod_{\substack{(r', s')\in E(H)\\ (r', s')\neq (r, s)}}f(x_{r'}, x_{s'})\prod_{\substack{v\in V(H)\\ v\neq r, s}}dx_v.
\end{equation}
The above definition appears rather complicated, but in essence may be identified with the homomorphism density (\ref{tt}), if we remove edge $(r, s)$ from $H$ and do not integrate out the associated vertices $x_r$ and $x_s$. For $(x, y) \in [0, 1]$, define
\begin{equation}
\label{Delta}
\Delta_H f(x, y):=\sum_{(r, s)\in E(H)} \Delta_{H, r, s} f(x, y),
\end{equation}
which corresponds to the total homomorphism density generated after all possible ways of removing one edge from $H$. We give some examples to illustrate this idea. When $H$ is an edge, $\Delta_H f(x, y) \equiv 1$. When $H$ is a triangle, by symmetry, $\Delta_H f(x, y)=3\int_{[0, 1]}f(x, z)f(y, z)dz$.

\subsection{Motivation}
\label{motivation}
We will consider as a motivating example the edge-triangle model, which is a $2$-parameter exponential random graph model obtained from (\ref{pmf2}) by setting  $H_1$ to  an edge and $H_2$ to  a triangle. Take $\beta=ru$ for $u=(u_1, u_2)\in {\mathbb S}$. Suppose that $u_1>0$ and $u_2<0$, i.e., the $2$-dimensional parameter vector $u$ is pointing towards the $4$th quadrant. Then ${\mathcal K}^u$ (which agrees with ${\mathcal K}^\beta$) consists of graphons $f$ that minimize $(u_1/u_2) e+t$, where $e=t(H_1, f)$ denotes the edge density and $t=t(H_2, f)$ denotes the triangle density of $f$, respectively. This implies that $f$, the maximizers of $T^u$ (and hence of $T^\beta$) must lie on the Razborov curve, which is the lower boundary of the feasible region $R$ of edge-triangle homomorphism densities. As $u_1 e+u_2 t$ is a linear function, $f$ must minimize over the convex hull $P$ of $R$. Since $R$ and $P$ only intersect at the points corresponding to Tur\'{a}n graphons, $f$ must be a Tur\'{a}n graphon. This important fact about the structure of ${\mathcal K}^\beta$ was further used to derive the maximizing graphons of ${\mathcal F}^\beta$ in \cite{YRF}. Using the boundedness of $I$, we found that ${\mathcal F}^\beta$ consists of graphons that can be made arbitrarily close to Tur\'{a}n graphons when the magnitude of $\beta$ is sufficiently large, and exactly which Tur\'{a}n graphon is favored by ${\mathcal F}^\beta$ depends on the direction of $\beta$. However, as nice as these results are, there is ambiguity concerning the optimal Tur\'{a}n graphon along the critical directions of $\beta$, which correspond to normal lines of the convex hull $P$. The subtlety might be due to the fact that Tur\'{a}n graphons, though close to our optimal graphons in cut distance, are not best approximations to the maximizing set ${\mathcal F}^\beta$. Generalizing from Tur\'{a}n graphons \cite{HN}, we say that a graphon $f$ is random-free if $f=\chi_A$ for some symmetric measurable subset $A$ of $[0,1]^2$. When $\beta$ is finite, Chatterjee and Diaconis \cite{CD2013} showed that the maximizing graphons in ${\mathcal F}^\beta$ are almost nowhere random-free -- that is the set $\{x\in[0,1]^2: f(x) \in \{0,1\}\}$ has zero measure. We are thus interested in finding coarse-grained graphons that are not random-free but close to ${\mathcal F}^\beta$ in cut distance, as they are sufficient to distinguish between candidates for the ground state. Resorting to perturbation analysis in the $|\beta| \to \infty$ regime, we will propose a method that keeps track of only the most significant characteristics of the maximizing graphons and  demonstrate its effectiveness.

\subsection{Assumptions}
\label{assumptions}
Let $R$ denote the range of $T$ in ${\mathbb R}^d$ and suppose that the boundary $\partial R$ is piecewise analytic. For $w \in \partial R$, let $F_R(w)$ denote the set of all feasible directions of $R$ at $w$, so that $v \in F_R(w)$ implies that there is $\epsilon > 0$ with $w + \epsilon v \in R$. The (internal) tangent cone $C_R(w)$ of $R$ at $w$ is then given by the closure of $F_R(w)$. See Figure \ref{et} for an illustration of this concept in the edge-triangle model. Let $\beta=r u$ for $u \in {\mathbb S}^{d-1}$, where $r$ is sufficiently large. Suppose that $T^u$ (and hence $T^\beta$ as $|\beta| \to \infty$) is maximized at a set of  random-free graphons $f =\chi_A$, i.e., ${\mathcal K}^u$ (which agrees with ${\mathcal K}^\beta$) consists of random-free graphons only. Further suppose that $f \in {\mathcal K}^u$ has the property that
      \begin{equation} \label{derivative1}
         \Delta_{H_i} f =  a_{i}\chi_{A}+b_{i}\chi_{B},
          \end{equation}
    where $\Delta_{H_i}f$ is defined as in (\ref{Delta}) and $B:=[0, 1]^2 \texttt{\char92} A$. We remark that these nice properties that we assumed are enjoyed beyond the edge-triangle model. Denote by $K_k$ a complete graph on $k$ vertices. Since the vertices of the convex hull $P$ of $R$ for $K_2$ and $K_n$ (any $n$), not just $K_2$ and $K_3$ as in the edge-triangle case, are given by Tur\'{a}n graphons, our argument will run through without much modification in these cases. Utilizing geometry of their respective convex hull \cite{EN}, similar analysis may be extended to more general models. We point out in particular that if $f$ is a Tur\'{a}n graphon, then (\ref{derivative1}) holds for any $H_i$ due to symmetry.

\begin{figure}
\centering
\includegraphics[clip=true, height=3in]{et.pdf}
\caption{Region of attainable edge ($e$) and triangle ($t$) densities for graphons. The upper
boundary is the curve $t=e^{3/2}$ and the lower boundary is a piecewise algebraic curve with
infinitely many concave pieces; see \cite{Razborov}. The tangent cones at complete bipartite graphon (Tur\'{a}n graphon with $2$ classes) and complete tripartite graphon (Tur\'{a}n graphon with $3$ classes) are displayed.} \label{et}
\end{figure}

\subsection{Perturbation of random-free graphons}
\label{perturb}
Consider a maximizing graphon $h \in {\mathcal F}^\beta$. When the magnitude of $\beta$ is sufficiently large, since $T^\beta$ diverges while $I$ stays bounded, ${\mathcal F}^\beta$ and ${\mathcal K}^\beta$ can be made arbitrarily close in the cut metric. From our assumptions in Section \ref{assumptions}, there thus exists a random-free graphon $f=\chi_A$ that is close to $h$ under the cut distance. We will further construct a non-random-free graphon $X_h$ that is close to $h$, simple enough, and yet still retains important information of $h$. Most importantly, we will show that for any finite simple graph $H_i$, $t(H_i, X_h)$ approximates $t(H_i, h)$ at least as well as $t(H_i, f)$ asymptotically. The following proposition is specific for random-free graphons and will be useful for our investigation.

\begin{proposition}
\label{prop}
Let $||\cdot||_1$ denote the $L^1$-norm. For a random-free graphon $f$, $d_\Box(f, h)=O(||f-h||_1)$.
\end{proposition}

\begin{proof}
By (\ref{cut}), it is clear that $d_\Box(f, h) \leq ||f-h||_1$. For the other direction, suppose $f=\chi_A$ for some symmetric measurable subset $A \subset [0, 1]^2$. If $A$ is a rectangle, then $||f-h||_1 \leq 3d_\Box(f, h)$. The general conclusion follows once we recognize that any subset $A$ of $[0, 1]^2$ may be approximated within $\epsilon>0$ by a finite union of disjoint rectangles.
\end{proof}

Let us first expand $t(H_i, h)$ around $f$. Denote a perturbation of $f$ by $g=h-f$. Corresponding to the regime $|\beta| \to \infty$, as explained earlier, $f$ may be chosen so that $|g| \in {\mathcal W}$. We have
\begin{equation}
\label{exph}
t(H_i, h)=t(H_i, f)+\int_{[0, 1]^2}g(x, y)\Delta_{H_i} f(x, y) dxdy+\text{remainder terms}.
\end{equation}
Under our assumption (\ref{derivative1}), we compute the first variation:
\begin{eqnarray}
\label{middle}
\int_{[0, 1]^2}g(x, y)\Delta_{H_i} f(x, y) dxdy&=&\int_{[0, 1]^2} \left(g(x, y)\chi_A+g(x, y)\chi_B\right)\left(a_i\chi_A+b_i\chi_B\right)dxdy\\\nonumber
&=&a_i\int_A g(x, y)dxdy+b_i\int_B g(x, y)dxdy.
\end{eqnarray}
Define
\begin{equation}
\label{ab}
a=\frac{1}{|A|} \int_A h(x, y) dxdy,
               \hspace{.2in}  b=\frac{1}{|B|} \int_B h(x, y) dxdy.
\end{equation}
Using that $f=\chi_A$ and $g=h-f$, (\ref{middle}) reduces to
\begin{equation}
\int_{[0, 1]^2}g(x, y)\Delta_{H_i} f(x, y) dxdy=a_i(a-1)|A|+b_ib|B|.
\end{equation}
The remainder terms (if they exist) give the homomorphism density after all possible ways of removing at least two edges from $H_i$, and so are bounded above by either $\int_{[0, 1]^3} |g(x, y)||g(x, z)|dxdydz$ or $\int_{[0, 1]^4} |g(x, y)||g(u, v)|dxdydudv$. Estimating the latter is easy.
\begin{equation}
\int_{[0, 1]^4} |g(x, y)||g(u, v)|dxdydudv=\left(\int_{[0, 1]^2}|g(x, y)|dxdy\right)^2=||g||_1^2,
\end{equation}
which is of negligible asymptotic order when compared to the first variation (\ref{middle}). The former corresponds to the $2$-star density of a graphon $|g|$ with edge density $||g||_1$. By \cite{AK}, for small $||g||_1$, the $2$-star density is bounded below by $||g||_1^2$ and above by $||g||_1/2+O(||g||_1^2)$, and the upper bound is achieved when $|g|$ is an anticlique of the form
\begin{equation}
|g(x,y)|=\left\{%
\begin{array}{ll}
    1, & \hbox{if $x>c$ and $y>c$,} \\
    0, & \hbox{otherwise,} \\
\end{array}%
\right.
\end{equation}
where $c=1-\sqrt{1-||g||_1}$. The lower bound for this remainder term is of desirable asymptotic order, but the upper bound is of the same order as the first variation. Recall that the maximizing graphon $h$ for a finite $\beta$ is uniformly bounded away from $0$ and $1$ \cite{CD2013}. The graphon $|g|=|h-f|$ is thus likely quite different from an anticlique. This implies that the $2$-star density of $|h-f|$ is of higher order than $||h-f||_1$ and does not achieve the upper bound. The phenomenon was confirmed for example by simulations for the edge-triangle model \cite{KRRS2}.

Define the averaged perturbation by $X_h = a\chi_A + b \chi_B $, where $a$ and $b$ are given in (\ref{ab}). $X_h$ may be viewed as a flattened out version of $h$. Since
\begin{equation}
\label{agree}
||h-f||_1=\int_A (1-h(x, y))dxdy+\int_B h(x, y)dxdy=(1-a)|A|+b|B|=||X_h-f||_1,
\end{equation}
$X_h$ is close to $f$ (and hence $h$) under the $L^1$ distance, and by Proposition \ref{prop}, also under the cut distance. Denote by $g'=X_h-f$. We perform the same expansion for $t(H_i, X_h)$ around $f$ as in the last paragraph:
\begin{equation}
\label{expXh}
t(H_i, X_h)=t(H_i, f)+\int_{[0, 1]^2}g'(x, y)\Delta_{H_i} f(x, y) dxdy+\text{remainder terms}.
\end{equation}
Following similar reasoning as in (\ref{middle}) and using the definition of $X_h$,
\begin{eqnarray}
\int_{[0, 1]^2}g'(x, y)\Delta_{H_i} f(x, y) dxdy&=&a_i \int_A g'(x, y)dxdy+b_i \int_B g'(x, y) dxdy\\\nonumber&=&a_i(a-1)|A|+b_ib|B|.
\end{eqnarray}
This says that $t(H_i, h)$ and $t(H_i, X_h)$ agree except for the remainder terms. As for $t(H_i, h)$, the remainder terms for $t(H_i, X_h)$ (if they exist) are bounded above by either $\int_{[0, 1]^3} |g'(x, y)||g'(x, z)| dxdydz$ or $\int_{[0, 1]^4} |g'(x, y)||g'(u, v)| dxdydudv$. Since $||g'||_1=||g||_1=(1-a)|A|+b|B|$ by (\ref{agree}), the latter is of asymptotic order $||g||_1^2$; while the former gives
\begin{equation}
\int_{[0, 1]^3} |g'(x, y)||g'(x, z)| dxdydz\leq \left(\max\{1-a, b\}\right)^2,
\end{equation}
and so is also of asymptotic order $||g||_1^2$. We conclude that $t(H_i, X_h)$ gives at least as good an asymptotic approximation for $t(H_i, h)$ as $t(H_i, f)$, and a better one when the error term related to the $2$-star density of $|h-f|$ may be dropped.

\subsection{Maximizing graphons}
\label{maximizing}
Let $\beta=r u$ for $u \in {\mathbb R}^d$ and $\|u\|=1$. Under our assumptions, ${\mathcal K}^u$ (which agrees with ${\mathcal K}^\beta$) is a set of random-free graphons. Take $h\in {\mathcal F}^{\beta}$. Then for $\epsilon > 0$ and $r$ large (corresponding to $|\beta| \to \infty$), $d_\Box(f, h)<\epsilon$ for some $f \in {\mathcal K}^u$, i.e., $f=\chi_A$ for some symmetric measurable subset $A$ of $[0, 1]^2$. By (\ref{exph}) and following analysis,
\begin{equation}
\psi^{\beta}(h)=rT^{u}(f)+r\left(\sum_{i=1}^d a_iu_i\right)(a-1)|A|+r\left(\sum_{i=1}^d b_iu_i\right)b|B|-I(h)+\text{remainder terms}.
\end{equation}
Similarly, by (\ref{expXh}) and following analysis,
\begin{equation}
\psi^{\beta}(X_h)=rT^{u}(f)+r\left(\sum_{i=1}^d a_iu_i\right)(a-1)|A|+r\left(\sum_{i=1}^d b_iu_i\right)b|B|-I(X_h)+\text{remainder terms}.
\end{equation}
In both equations above, the entropy $I$ is bounded in contrast with the energy contribution $T^\beta$. Except that $h$ is close to $f$ (and hence $X_h$) in cut distance, we do not have enough information regarding the structure of $h$. Rather than maximizing $\psi^\beta(h)$ over all possible graphons $h$ directly, we will maximize $\psi^\beta(X_h)$ over $2$-parameter families $0 \leq a, b \leq 1$. From the heuristics in Section \ref{perturb}, this is an effective method to approximate the optimal graphon. More than that, in certain situations (for example the edge-triangle model to be discussed in detail in Section \ref{edgetriangle}), we will see that keeping $2$ parameters is not only effective but also sufficient. Notice that
\begin{equation}
\label{entropyapprx}
I(X_h)=I( a\chi_A + b \chi_B   )=  I(a  )  |A| + I(b) |B|.
\end{equation}
We rewrite $X_h$:
\begin{equation}
\label{flat}
\psi^{\beta}(X_h)=rT^{u}(f)+\left(r\left(\sum_{i=1}^d a_iu_i\right)(a-1)-I(a)\right)|A|+\left(r\left(\sum_{i=1}^d b_iu_i\right)b-I(b)\right)|B|+\text{remainder terms}.
\end{equation}
Maximizing each first variation term, we have
\begin{equation}
\label{optimal}
a = \frac{1}{1 + e^{ - 2r \sum_{i=1}^d a_iu_i} }
               , \hspace{.2in}    b = \frac{1}{1 + e^{ - 2r \sum_{i=1}^d b_iu_i} }.
\end{equation}
Under this choice of $a$ and $b$ and provided we can ignore the remainder terms, $\psi^{\beta}(X_h)$ is strictly bigger than $rT^u(f)$, which is the random-free graphon approximation corresponding to $a=1$ and $b=0$. Let us verify that $\psi^\beta(X_h)$ is indeed strictly bigger than $rT^u(f)$ by rigorously managing the error term in (\ref{flat}). As shown earlier in Section \ref{perturb}, the remainder terms are of higher order:
\begin{multline} \label{ruZ3}
      \psi^{\beta}(X_h)=rT^{u}(f)+\left(r\left(\sum_{i=1}^d a_iu_i\right)(a-1)+O(r(a-1)^2)-I(a)\right)|A|\\+\left(r\left(\sum_{i=1}^d b_iu_i\right)b+O(rb^2)-I(b)\right)|B|.
\end{multline}
For any $\epsilon > 0$, there exists large enough $r$ so that $||X_h-f||_1$ is sufficiently small, making $0<1-a<\epsilon$ and $0<b<\epsilon$. Applying these $\epsilon$ bounds in (\ref{ruZ3}) gives
      \begin{equation} \label{ruZ5}
      \psi^{\beta}(X_h) \geq  rT^{u}(f)+\left(r\left(\sum_{i=1}^d a_iu_i+\epsilon\right)(a-1)-I(a)\right)|A|+\left(r\left(\sum_{i=1}^d b_iu_i-\epsilon\right)b-I(b)\right)|B|.
     \end{equation}
Maximizing each first variation term as previously, this yields $\psi^{\beta}(X_h)>rT^u(f)$.  We conclude that, as expected, the addition of one more parameter improves the random-free graphon estimation.

\section{The edge-triangle model}
\label{edgetriangle}
Denote by $K_k$ a complete graph on $k$ vertices. The edge-triangle model is a $2$-parameter exponential random graph model obtained by taking $H_1$ an edge ($K_2$) and $H_2$ a triangle ($K_3$) in (\ref{pmf2}). Consider the set $R=\{(t(K_2, f), t(K_3, f)), f\in {\mathcal W}\}$ of all realizable values of the edge ($e$) and triangle ($t$) homomorphism densities as the graphon $f$ varies over the entire graphon space ${\mathcal W}$. See Figure \ref{et}. The upper boundary curve of $R$ is given by the equation $t=e^{3/2}$, and can be derived using the Kruskal-Katona theorem (see Section 16.3 of \cite{Lov}). The lower boundary curve is trickier. The trivial lower bound of $t=g_1(e):=0$, corresponding to the horizontal segment, is attainable at any $0 \leq e \leq 1/2$ by graphons describing the possibly asymptotic edge density of subgraphs of complete bipartite graphs (Tur\'{a}n graphon with $2$ classes). For $e \geq 1/2$, the optimal bound was obtained by Razborov \cite{Razborov}, who
established, using the flag algebra calculus, that for $(k-1)/k
\leq e \leq k/(k+1)$ with $k \geq 2$,
\begin{equation}
\label{Ra} t \geq g_k(e):=\frac{(k-1) \left(k-2\sqrt{k(k-e(k+1))}\right)
\left(k+ \sqrt{k(k-e(k+1))} \right)^2}{k^2 (k+1)^2}.
\end{equation}
All the curve segments $g_k(e)$ describing the nontrivial part of the lower boundary of $R$ are strictly concave. For $k=1, 2, \ldots$, we set $v_k:=(e_k, t_k)=(t(K_2, f^{K_{k+1}}), t(K_3, f^{K_{k+1}}))$, where explicitly,
\begin{equation}
f^{K_{k+1}}(x, y)=\left\{%
\begin{array}{ll}
    1, & \hbox{if $\lceil (k+1)x\rceil \neq \lceil (k+1)y \rceil$,} \\
    0, & \hbox{otherwise,} \\
\end{array}%
\right. \quad (x,y) \in [0,1]^2
\end{equation}
is the Tur\'{a}n graphon with $k+1$ classes. Thus
\begin{equation}
e_k=\frac{k}{k+1}, \hspace{.2in} t_k=\frac{k(k-1)}{(k+1)^2}.
\end{equation}
For $k=1, 2, \ldots$, let $L_k$ be the line segment joining vertices $v_k$ and $v_{k+1}$ of neighboring Tur\'{a}n graphons. These infinitely many line segments form the convex hull $P$ of $R$, and the length of $L_k$ decreases monotonically to zero as $k$ gets large.

The normal vectors to $L_k$,
\begin{equation}
\label{ok}
o_k=\left(1, -\frac{(k+1)(k+2)}{k(3k+5)}\right)
\end{equation}
are the critical directions of the edge-triangle model. Let $\beta=ro_k$ and take $r\rightarrow \infty$. While the vectors $o_k$ (\ref{ok}) are not normalized as in our derivation (see Section \ref{approximation}), this can be easily adjusted by adapting $r$. Then Tur\'{a}n graphons with $k+1$ and $k+2$ classes both belong to ${\mathcal K}^{o_k}$ (which agrees with ${\mathcal K}^\beta$), and we concluded in \cite{YRF} that a typical graph sampled from the model may behave like either a Tur\'{a}n graphon with $k+1$ classes or a Tur\'{a}n graphon with $k+2$ classes, with no clear preference. Though already quite informative, as explained earlier in Section \ref{motivation}, this result remains somewhat unsatisfactory because it does not indicate whether both such graphons are actually realizable in the limit and in what manner. Notice that underneath our investigation, there is an ordered double asymptotic framework, in the sense that the network size $n$ goes to infinity first followed by the divergence of the parameters $\beta$. In hope of resolving this rather subtle ambiguity within the edge-triangle model, the ``other'' order was also examined in \cite{YRF}, where we first let the magnitude of $\beta$ increase to infinity so as to isolate a simpler sub-model and then study its limiting properties as $n$ grows. Both ordered asymptotics~imply a nearly identical convergence in probability in the cut metric along the noncritical directions. Under the ``original'' as well as the ``reversed'' order, there exist (possibly different) subsequences of the form $\{n_i, \beta_{1,i}, \beta_{2,i}\}$, with $n_i \rightarrow \infty$, $\beta_{1,i}\rightarrow \infty$ and $\beta_{2,i}\rightarrow -\infty$ for $i=1, 2, \ldots$, where the edge-triangle model converges to some Tur\'{a}n graphon specified by the direction of the parameters $(\beta_1, \beta_2)$. Additionally, under the ``reversed'' asymptotics, a detailed categorization of the limiting behavior of the edge-triangle model was obtained: When $\beta$ diverges along the critical direction $o_k$, a typical sampled graph more likely resembles a Tur\'{a}n graphon with $k+2$ classes than with $k+1$ classes. Now that we are equipped with refined perturbation analysis, we would like to sharpen our results under the ``original'' asymptotics and inquire whether the same type of discontinuity in natural parametrization exists.

Let us make a further remark before carrying out the detailed calculations. In the physics literature, people are often interested in cases where the parameter $\beta$ depends on $n$ (some averages need to be satisfied for every $n$). In these models, in place of the normalization constant (free energy density), the relative entropy plays a central role. Analogous (but more complicated) maximization problems and concentration of measure results have been established, which lead to classifications of ensemble equivalence between the microcanonical ensemble and the canonical ensemble. The perturbative methods explored in the current paper are expected to apply in these general parameter situations. In some cases, the perturbation would still be around random-free graphons, and our argument will run through without much adaptation \cite{HMRS}. In some other cases however, the perturbation would be around graphons admitting more intricate structures, and serious future work is needed\cite{GHR,PN}.

\subsection{Perturbation analysis}
Let $\beta=ro_k$. We will compare $\psi^{\beta}(X)$ and $\psi^{\beta}(Y)$, where $X$ is the flattened out graphon close to $f^{K_{k+1}}$ (Tur\'{a}n graphon with $k+1$ classes) and $Y$ is the flattened out graphon close to $f^{K_{k+2}}$ (Tur\'{a}n graphon with $k+2$ classes). Both $X$ and $Y$ are constructed with the optimal perturbation values (\ref{optimal}). We set $u=o_k$ in our calculations below. For $i=1, 2$, we compute $a_i^X$, $b_i^X$, $a_i^Y$, and $b_i^Y$ for $f^{K_{k+1}}$ and $f^{K_{k+2}}$ in (\ref{derivative1}), where $H_1=K_2$ and $H_2=K_3$. As pointed out earlier in Section \ref{approximation}, $a_1^X=b_1^X=a_1^Y=b_1^Y=1$. For notational convenience, from now on we denote $f^{K_{k+1}}$ by $T_{k}$ and $\chi_{[0, 1]^2}-f^{K_{k+1}}$ by $D_{k}$. Then $T_k$ and $D_k$ are indicator functions associated with sets of measure $k/(k+1)$ and $1/(k+1)$, respectively. Reconfirming our assumption in (\ref{derivative1}),
\begin{equation}
\Delta_{K_3} T_k=3\int_{[0, 1]}T_k(x, z)T_k(y, z)dz=\frac{3(k-1)}{k+1}T_k+\frac{3k}{k+1}D_k,
\end{equation}
which gives
\begin{equation}
\label{values}
a_2^X=\frac{3(k-1)}{k+1}, \hspace{.2in} b_2^X=\frac{3k}{k+1},
\end{equation}
\begin{equation*}
a_2^Y=\frac{3k}{k+2}, \hspace{.2in} b_2^Y=\frac{3(k+1)}{k+2}.
\end{equation*}
This yields
\begin{equation}
\label{abs}
\sum_{i=1}^2 a_i^X u_i=\frac{2(k+3)}{k(3k+5)}, \hspace{.2in} \sum_{i=1}^2 b_i^X u_i=-\frac{1}{3k+5},
\end{equation}
\begin{equation*}
\sum_{i=1}^2 a_i^Y u_i=\frac{2}{3k+5}, \hspace{.2in} \sum_{i=1}^2 b_i^Y u_i=-\frac{k+3}{k(3k+5)}.
\end{equation*}

\begin{figure}
\centering
\includegraphics[height=6in, angle=-90]{cone.pdf}
\caption{The (internal) tangent cone $C_R(v_k)$, consisting of graphons of the type $X_h=aT_{k}+bD_k$. Taking $a=1, b=0$ gives the Tur\'{a}n graphon with $k+1$ classes, whose edge-triangle densities are $v_k=(e_k, t_k)$. Here $g_k'(e_k-0)$ and $g_k'(e_k+0)$ indicate the left and right hand side derivatives of $g_k$ (\ref{Ra}) at $e_k$, respectively.}
\label{cone}
\end{figure}

We check the remainder terms for $\psi^\beta(X)$ after first order perturbation. Similar analysis will work for $\psi^\beta(Y)$ and we skip the details. Consider the (internal) tangent cone $C_R(v_k)$, consisting of graphons of the type $X_h=aT_k+bD_k$. See Figure \ref{cone}. Then
\begin{equation}
\label{edge}
t(K_2, X_h)=\int_{[0, 1]^2}\left(aT_k(x, y)+bD_k(x, y)\right)dxdy=\frac{ak+b}{k+1}.
\end{equation}
For $(x, y) \in [0, 1]^2$, we have
\begin{multline}
\int_{[0, 1]}X_h(x, z)X_h(y, z)dz=\int_{[0, 1]}\left((aT_k(x, z)+bD_k(x, z))(aT_k(y, z)+bD_k(y, z))\right)dz\\
=\left(a^2\frac{k-1}{k+1}+2ab\frac{1}{k+1}\right)T_k(x, y)+\left(a^2\frac{k}{k+1}+b^2\frac{1}{k+1}\right)D_k(x, y),
\end{multline}
which implies that
\begin{equation}
\label{triangle}
t(K_3, X_h)=\int_{[0, 1]^3}X_h(x, y)X_h(x, z)X_h(y, z)dxdydz
\end{equation}
\begin{eqnarray*}
&=&\int_{[0, 1]}\left(aT_k(x, y)+bD_k(x, y)\right)\left(\left(a^2\frac{k-1}{k+1}+2ab\frac{1}{k+1}\right)T_k(x, y)+\left(a^2\frac{k}{k+1}+b^2\frac{1}{k+1}\right)D_k(x, y)\right)dxdy\\\nonumber
&=&a^3\frac{k(k-1)}{(k+1)^2}+3a^2b\frac{k}{(k+1)^2}+b^3\frac{1}{(k+1)^2}. \nonumber
\end{eqnarray*}
This yields
\begin{multline}
\label{exact}
\psi^\beta(X_h)=rt(K_2, X_h)-r\frac{(k+1)(k+2)}{k(3k+5)}t(K_3, X_h)-I(a)\frac{k}{k+1}-I(b)\frac{1}{k+1}\\
=r\frac{ak+b}{k+1}-r\frac{(k+1)(k+2)}{k(3k+5)}\left(a^3\frac{k(k-1)}{(k+1)^2}+3a^2b\frac{k}{(k+1)^2}+b^3\frac{1}{(k+1)^2}\right)-I(a)\frac{k}{k+1}-I(b)\frac{1}{k+1}.
\end{multline}
Using (\ref{optimal}) and (\ref{abs}) for $a$ and $b$, (\ref{exact}) gives $\psi^\beta(X)$. After first order perturbation as in (\ref{flat}), the remainder terms are bounded by
\begin{equation}
\text{constant}\cdot r\cdot (\max\{1-a, b\})^2 \approx \text{constant} \cdot r \cdot \max\big \{\exp\left(-\frac{8r(k+3)}{k(3k+5)}\right), \exp\left(-\frac{4r}{3k+5}\right)\big \}.
\end{equation}
The following lemma is useful for our asymptotic derivation.

\begin{lemma}
\label{approx}
Let $I$ be defined as in (\ref{I}). As $x \rightarrow \infty$,
\begin{equation}
\label{simple1}
-\frac{xe^{-x}}{2(1+e^{-x})}-I(\frac{1}{1+e^{-x}}) \approx \frac{e^{-x}}{2},
\end{equation}
\end{lemma}
\noindent and as $x \rightarrow -\infty$,
\begin{equation}
\label{simple2}
\frac{x}{2(1+e^{-x})}-I(\frac{1}{1+e^{-x}}) \approx \frac{e^{x}}{2}.
\end{equation}
\begin{proof}
We recognize that after simplification, the left hand side of (\ref{simple1}) becomes $\log (1+e^{-x})/2$, and the left hand side of (\ref{simple2}) becomes $x/2+\log (1+e^{-x})/2$. The rest is immediate.
\end{proof}

Applying Lemma \ref{approx} to the first order perturbation terms, we see that conforming to our heuristic analysis in Sections \ref{perturb} and \ref{maximizing}, the remainder terms are indeed of negligible order. For $r$ large enough, we have
\begin{equation}
\label{first}
\psi^{\beta}(X) \approx rT^{o_k}(f)+\frac{k}{2(k+1)}\exp\left(-\frac{4r(k+3)}{k(3k+5)}\right)+\frac{1}{2(k+1)}\exp\left(-\frac{2r}{3k+5}\right),
\end{equation}
\begin{equation*}
\psi^{\beta}(Y) \approx rT^{o_k}(f)+\frac{k+1}{2(k+2)}\exp\left(-\frac{4r}{3k+5}\right)+\frac{1}{2(k+2)}\exp\left(-\frac{2r(k+3)}{k(3k+5)}\right).
\end{equation*}
Since $\sum_{i=1}^2 b_i^X u_i=-1/(3k+5)$ has the smallest absolute value among all equations in (\ref{abs}), $\psi^\beta(X)>\psi^{\beta}(Y)$ both in terms of first order perturbation (\ref{first}) and the exact value.

Denote by $\psi^\beta_\text{opt}(X_h)$ and $\psi^\beta_\text{opt}(Y_h)$ the exact optimal value of $\psi^\beta$ within the (internal) tangent cone $C_R(v_k)$ and $C_R(v_{k+1})$, respectively. Let us compare $\psi^\beta_\text{opt}(X_h)$ with $\psi^\beta_\text{opt}(Y_h)$. From our heuristic argument, for any $\epsilon>0$ and $r$ large enough, (\ref{ruZ5}) is satisfied:
\begin{equation}
\psi^{\beta}(X_h) \geq  rT^{o_k}(f)+\left(r\left(\sum_{i=1}^2 a_i^X u_i+\epsilon\right)(a-1)-I(a)\right)\frac{k}{k+1}+\left(r\left(\sum_{i=1}^2 b_i^X u_i-\epsilon\right)b-I(b)\right)\frac{1}{k+1}.
\end{equation}
Following similar reasoning, we also find a bound in the other direction:
\begin{equation}
\psi^{\beta}(X_h) \leq  rT^{o_k}(f)+\left(r\left(\sum_{i=1}^2 a_i^X u_i-\epsilon\right)(a-1)-I(a)\right)\frac{k}{k+1}+\left(r\left(\sum_{i=1}^2 b_i^X u_i+\epsilon\right)b-I(b)\right)\frac{1}{k+1}.
\end{equation}
Maximizing each first variation term and applying Lemma \ref{approx} as previously, this says that $\psi^\beta(X_h)$ is asymptotically bounded below by
\begin{equation}
rT^{o_k}(f)+\frac{k}{2(k+1)}\exp\left(-2r\left(\frac{2(k+3)}{k(3k+5)}+\epsilon\right)\right)+\frac{1}{2(k+1)}\exp\left(-2r\left(\frac{1}{3k+5}+\epsilon\right)\right),
\end{equation}
and above by
\begin{equation}
rT^{o_k}(f)+\frac{k}{2(k+1)}\exp\left(-2r\left(\frac{2(k+3)}{k(3k+5)}-\epsilon\right)\right)+\frac{1}{2(k+1)}\exp\left(-2r\left(\frac{1}{3k+5}-\epsilon\right)\right),
\end{equation}
Similar analysis works for $\psi^\beta(Y_h)$ and we skip the details. Since $\epsilon>0$ can be taken arbitrarily small and $\sum_{i=1}^2 b_i^X u_i=-1/(3k+5)$ has the smallest absolute value among all equations in (\ref{abs}), $\psi^\beta_{\text{opt}}(X_h)>\psi^{\beta}_{\text{opt}}(Y_h)$. As discussed earlier in Section \ref{subgraph}, by Theorems 3.1 and 3.2 in Chatterjee and Diaconis \cite{CD2013}, we conclude that when $\beta$ diverges along the critical direction $o_k$, a typical sampled graph more likely resembles a Tur\'{a}n graphon with $k+1$ classes than with $k+2$ classes. See Table \ref{table1} and Figure \ref{figure1}.

\begin{theorem}
\label{main}
Consider the edge-triangle exponential random graph model, obtained by setting in (\ref{pmf2}) $H_1$ an edge and $H_2$ a triangle. For $k\geq 1$, let $\beta=ro_k$, where $o_k$ is the critical direction (\ref{ok}) and $r$ is sufficiently large. Then in the large $n$ limit, a typical graph drawn from the model behaves like a Tur\'{a}n graphon with $k+1$ classes,
\begin{equation}
\lim_{r\rightarrow \infty} \sup_{\tilde{f} \in {\mathcal F}^\beta} \delta_{\Box}(\tilde{f}, \tilde{f}^{K_{k+1}})=0.
\end{equation}
\end{theorem}

\begin{table}
\begin{center}
\begin{tabular}{c||cccccc}
 & $a$ & $b$ & $\psi^{\beta}$ & $a_{\text{opt}}$ & $b_{\text{opt}}$ & $\psi^\beta_{\text{opt}}$ \\
\hline \hline \\
bipartite & 0.999999998 & $0.0759$ & 5.0197 & \text{almost} 1 & 0.069 & 5.019 \\
tripartite & $0.9933$ & $0.0000454$ & 5.0022 & 0.9943 & 0.000064 & 5.0021 \\ \\
\end{tabular}
\end{center}
\caption{Asymptotic comparison for perturbation around neighboring Tur\'{a}n graphons, where $\beta=(10, -7.5)$ diverges along the critical direction $o_1=(1, -0.75)$. $a$, $b$ and $\psi^{\beta}$ are calculated from (\ref{optimal}) and (\ref{first}) with first order perturbation; $a_{\text{opt}}$, $b_{\text{opt}}$ and $\psi^\beta_{\text{opt}}$ are based on numerical optimization for (\ref{exact}). The bipartite feature is associated with a bigger limiting normalization constant and is favored, matching the asymptotic predictions of Theorem \ref{main}.} \label{table1}
\end{table}

\subsection{Geometric interpretation}

We proceed further and examine the effect of infinitesimal perturbation on the associated edge and triangle densities. From (\ref{edge}) and (\ref{triangle}),
\begin{equation}
\label{tangent}
\left.\frac{\partial}{\partial b} \left(t(K_2, X_h), t(K_3, X_h)\right)\right|_{a=1, b=0} = \left( \frac{1}{k+1},  \frac{3k}{(k+1)^2}\right),
\end{equation}
\begin{equation*}
             \left.\frac{\partial}{\partial a} \left( t(K_2, X_h), t(K_3, X_h)\right)\right|_{a=1, b=0} = \left( \frac{k}{k+1},  \frac{3k(k-1)}{(k+1)^2}\right).
\end{equation*}
Let us present another perspective on this calculation incorporating assumption (\ref{derivative1}), which may be employed to derive infinitesimal variations for more complicated homomorphism densities. The idea appeared in our heuristic analysis before and we make it explicit here. For any $H_i$ so that (\ref{derivative1}) is satisfied,
\begin{multline}
\left.\frac{\partial}{\partial b} t(H_i, X_h)\right|_{a=1, b=0}=\int_{[0, 1]^2} D_k(x, y)\Delta_{H_i}T_k(x, y)dxdy\\
=\int_{[0, 1]^2} D_k(x, y)\left(a_iT_k(x, y)+b_iD_k(x, y)\right)dxdy=b_i\frac{1}{k+1}.
\end{multline}
\begin{multline}
\left.\frac{\partial}{\partial a} t(H_i, X_h)\right|_{a=1, b=0}=\int_{[0, 1]^2} T_k(x, y)\Delta_{H_i}T_k(x, y)dxdy\\
=\int_{[0, 1]^2} T_k(x, y)\left(a_iT_k(x, y)+b_iD_k(x, y)\right)dxdy=a_i\frac{k}{k+1}.
\end{multline}
Using $a_1=b_1=1$, $a_2=3(k-1)/(k+1)$ and $b_2=3k/(k+1)$ (\ref{values}), we recover the partial derivatives calculated above. In particular, we recognize that $-(a_1, a_2)$ points along the left tangent line and $(b_1, b_2)$ points along the right tangent line, while the critical direction $u=o_k$ is the normal vector to the line segment $L_k$ that connects neighboring vertices $v_k$ and $v_{k+1}$. This offers a geometric justification of the ($\pm$) signs in (\ref{abs}).

\begin{figure}
\centering
\includegraphics[width=4in]{simu.pdf}
\caption{A simulated realization of the exponential random graph model on $60$ nodes with edges and triangles as sufficient statistics, where $\beta=(10, -7.5)$ diverges along the critical direction $o_1=(1, -0.75)$. $n$ is of much bigger magnitude than $\beta$ in reflection of the double asymptotic order. The simulated graph displays bipartite feature with edge density $0.557$, matching the asymptotic predictions of Theorem \ref{main}.}
\label{figure1}
\end{figure}

Recall that the lower boundary of attainable edge-triangle densities is a piecewise algebraic curve with infinitely many concave
pieces $g_k$ (\ref{Ra}), and the connection point of $g_k$ and $g_{k+1}$ is $v_k=(e_k, t_k)$. See Figure \ref{et}. We compute
\begin{equation}
g_k'(e)=\frac{3(k-1)}{k(k+1)}\left(k+\sqrt{k(k-e(k+1))}\right),
\end{equation}
which implies that the left and right hand side derivatives at $e_k=k/(k+1)$ are respectively given by
\begin{equation}
g_k'(e_k-0)=\frac{3(k-1)}{k+1}, \hspace{.2in} g_{k+1}'(e_k+0)=\frac{3k}{k+1}.
\end{equation}
The partial derivative vectors in (\ref{tangent}) thus delineate the boundary of the (internal) tangent cone $C_R(v_k)$. In other words, $C_R(v_k)$ spans all possible infinitesimal variations at the Tur\'{a}n graphon with $k+1$ classes. See Figure \ref{cone}. Since the optimizing graphon associated with large enough $\beta$ must lie within the tangent cone of some Tur\'{a}n graphon, it may be represented by a linear combination of Erd\H{o}s-R\'{e}nyi and Tur\'{a}n graphons. Even though the graphon representation may not be unique, optimizing over all possible combinations provides insight into the structure of the maximizing set. Keeping track of the Erd\H{o}s-R\'{e}nyi and Tur\'{a}n characteristics in the edge-triangle model is thus not only an effective but also sufficient method to estimate the normalization constant, and gives evidence of discontinuity of the natural parametrization along the critical directions $o_k$ in the limit as $n$ and then $r$ tend to infinity. This demonstrates the occurrence of discontinuous phase transitions in the edge-triangle model.

\section*{Acknowledgements}
The authors are very grateful to the anonymous referee for the invaluable
suggestions that greatly improved the quality of this paper. Mei Yin thanks Sukhada Fadnavis for helpful conversations. Rajinder Mavi was supported by a postdoctoral fellowship from the Michigan State University Institute for Theoretical and Mathematical Physics. Mei Yin's research was partially supported by NSF grant DMS-1308333.